Skip to content

Commit

Permalink
First pass at adding land/sea masking to evolution
Browse files Browse the repository at this point in the history 
- Read masks in from the NetCDF grid file rather than initial file
- Masks are set to 1 if ocean, 0 if land
- Add psi-point masks
  • Loading branch information
@hklion
hklion committed Jul 30, 2024
1 parent 6c8e4e3 commit df2cfbe
Show file tree
Hide file tree
Showing 20 changed files with 268 additions and 96 deletions.
31 changes: 25 additions & 6 deletions Source/IO/ReadFromInitNetcdf.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,9 +11,7 @@ read_data_from_netcdf (int /*lev*/,
const std::string& fname,
FArrayBox& NC_temp_fab, FArrayBox& NC_salt_fab,
FArrayBox& NC_xvel_fab, FArrayBox& NC_yvel_fab,
FArrayBox& NC_ubar_fab, FArrayBox& NC_vbar_fab,
FArrayBox& NC_mskr_fab, FArrayBox& NC_msku_fab,
FArrayBox& NC_mskv_fab)
FArrayBox& NC_ubar_fab, FArrayBox& NC_vbar_fab)
{
amrex::Print() << "Loading initial solution data from NetCDF file " << fname << std::endl;

Expand All @@ -27,9 +25,6 @@ read_data_from_netcdf (int /*lev*/,
NC_fabs.push_back(&NC_yvel_fab); NC_names.push_back("v"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 3
NC_fabs.push_back(&NC_ubar_fab), NC_names.push_back("ubar"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); // 4
NC_fabs.push_back(&NC_vbar_fab); NC_names.push_back("vbar"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); // 5
NC_fabs.push_back(&NC_mskr_fab); NC_names.push_back("mask_rho"); NC_dim_types.push_back(NC_Data_Dims_Type::SN_WE); // 6
NC_fabs.push_back(&NC_msku_fab); NC_names.push_back("mask_u"); NC_dim_types.push_back(NC_Data_Dims_Type::SN_WE); // 7
NC_fabs.push_back(&NC_mskv_fab); NC_names.push_back("mask_v"); NC_dim_types.push_back(NC_Data_Dims_Type::SN_WE); // 8

// Read the netcdf file and fill these FABs
BuildFABsFromNetCDFFile<FArrayBox,Real>(domain, fname, NC_names, NC_dim_types, NC_fabs);
Expand Down Expand Up @@ -91,4 +86,28 @@ read_coriolis_from_netcdf (int /*lev*/,
// Read the netcdf file and fill these FABs
BuildFABsFromNetCDFFile<FArrayBox,Real>(domain, fname, NC_names, NC_dim_types, NC_fabs);
}

void
read_masks_from_netcdf (int /*lev*/,
const Box& domain,
const std::string& fname,
FArrayBox& NC_mskr_fab,
FArrayBox& NC_msku_fab,
FArrayBox& NC_mskv_fab,
FArrayBox& NC_mskp_fab)
{
amrex::Print() << "Loading masks from NetCDF file " << fname << std::endl;

Vector<FArrayBox*> NC_fabs;
Vector<std::string> NC_names;
Vector<enum NC_Data_Dims_Type> NC_dim_types;

NC_fabs.push_back(&NC_mskr_fab ) ; NC_names.push_back("mask_rho") ; NC_dim_types.push_back(NC_Data_Dims_Type::SN_WE); // 0
NC_fabs.push_back(&NC_msku_fab ) ; NC_names.push_back("mask_u") ; NC_dim_types.push_back(NC_Data_Dims_Type::SN_WE); // 1
NC_fabs.push_back(&NC_mskv_fab ) ; NC_names.push_back("mask_v") ; NC_dim_types.push_back(NC_Data_Dims_Type::SN_WE); // 2
NC_fabs.push_back(&NC_mskp_fab ) ; NC_names.push_back("mask_psi") ; NC_dim_types.push_back(NC_Data_Dims_Type::SN_WE); // 3

// Read the netcdf file and fill these FABs
BuildFABsFromNetCDFFile<FArrayBox,Real>(domain, fname, NC_names, NC_dim_types, NC_fabs);
}
#endif // ROMSX_USE_NETCDF
6 changes: 4 additions & 2 deletions Source/Initialization/REMORA_init.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -189,11 +189,13 @@ REMORA::set_2darrays (int lev)
});
}

print_state(*vec_mskv[lev].get(), IntVect(10,0,0));
FillPatch(lev, t_new[lev], *vec_ubar[lev], GetVecOfPtrs(vec_ubar), BdyVars::ubar,0,false,false);
FillPatch(lev, t_new[lev], *vec_vbar[lev], GetVecOfPtrs(vec_vbar), BdyVars::vbar,0,false,false);

FillPatch(lev, t_new[lev], *vec_msku[lev], GetVecOfPtrs(vec_ubar), BdyVars::null,0,true,false);
FillPatch(lev, t_new[lev], *vec_mskv[lev], GetVecOfPtrs(vec_vbar), BdyVars::null,0,true,false);
FillPatchNoBC(lev, t_new[lev], *vec_msku[lev], GetVecOfPtrs(vec_msku), BdyVars::null,0,true,false);
FillPatchNoBC(lev, t_new[lev], *vec_mskv[lev], GetVecOfPtrs(vec_mskv), BdyVars::null,0,true,false);
print_state(*vec_mskv[lev].get(), IntVect(10,0,0));
}

void
Expand Down
93 changes: 60 additions & 33 deletions Source/Initialization/REMORA_init_from_netcdf.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,9 +16,12 @@ void
read_data_from_netcdf (int /*lev*/, const Box& domain, const std::string& fname,
FArrayBox& NC_temp_fab, FArrayBox& NC_salt_fab,
FArrayBox& NC_xvel_fab, FArrayBox& NC_yvel_fab,
FArrayBox& NC_ubar_fab, FArrayBox& NC_vbar_fab,
FArrayBox& NC_mskr_fab,
FArrayBox& NC_msku_fab, FArrayBox& NC_mskv_fab);
FArrayBox& NC_ubar_fab, FArrayBox& NC_vbar_fab);

void
read_masks_from_netcdf (int /*lev*/, const Box& domain, const std::string& fname,
FArrayBox& NC_mskr_fab, FArrayBox& NC_msku_fab,
FArrayBox& NC_mskv_fab, FArrayBox& NC_mskp_fab);

Real
read_bdry_from_netcdf (const Box& domain, const std::string& fname,
Expand All @@ -33,17 +36,12 @@ init_state_from_netcdf (int lev,
FArrayBox& temp_fab, FArrayBox& salt_fab,
FArrayBox& x_vel_fab, FArrayBox& y_vel_fab,
FArrayBox& ubar_fab, FArrayBox& vbar_fab,
FArrayBox& mskr_fab,
FArrayBox& msku_fab, FArrayBox& mskv_fab,
const Vector<FArrayBox>& NC_temp_fab,
const Vector<FArrayBox>& NC_salt_fab,
const Vector<FArrayBox>& NC_xvel_fab,
const Vector<FArrayBox>& NC_yvel_fab,
const Vector<FArrayBox>& NC_ubar_fab,
const Vector<FArrayBox>& NC_vbar_fab,
const Vector<FArrayBox>& NC_mskr_fab,
const Vector<FArrayBox>& NC_msku_fab,
const Vector<FArrayBox>& NC_mskv_fab);
const Vector<FArrayBox>& NC_vbar_fab);

void
read_bathymetry_from_netcdf (int lev, const Box& domain, const std::string& fname,
Expand Down Expand Up @@ -81,18 +79,13 @@ REMORA::init_data_from_netcdf (int lev)
Vector<FArrayBox> NC_yvel_fab ; NC_yvel_fab.resize(num_boxes_at_level[lev]);
Vector<FArrayBox> NC_ubar_fab ; NC_ubar_fab.resize(num_boxes_at_level[lev]);
Vector<FArrayBox> NC_vbar_fab ; NC_vbar_fab.resize(num_boxes_at_level[lev]);
Vector<FArrayBox> NC_mskr_fab ; NC_mskr_fab.resize(num_boxes_at_level[lev]);
Vector<FArrayBox> NC_msku_fab ; NC_msku_fab.resize(num_boxes_at_level[lev]);
Vector<FArrayBox> NC_mskv_fab ; NC_mskv_fab.resize(num_boxes_at_level[lev]);

for (int idx = 0; idx < num_boxes_at_level[lev]; idx++)
{
read_data_from_netcdf(lev, boxes_at_level[lev][idx], nc_init_file[lev][idx],
NC_temp_fab[idx], NC_salt_fab[idx],
NC_xvel_fab[idx], NC_yvel_fab[idx],
NC_ubar_fab[idx], NC_vbar_fab[idx],
NC_mskr_fab[idx], NC_msku_fab[idx],
NC_mskv_fab[idx]);
NC_ubar_fab[idx], NC_vbar_fab[idx]);
}


Expand All @@ -113,20 +106,13 @@ REMORA::init_data_from_netcdf (int lev)
FArrayBox &yvel_fab = (*yvel_new[lev])[mfi];
FArrayBox &ubar_fab = (*vec_ubar[lev])[mfi];
FArrayBox &vbar_fab = (*vec_vbar[lev])[mfi];
FArrayBox &mskr_fab = (*vec_mskr[lev])[mfi];
FArrayBox &msku_fab = (*vec_msku[lev])[mfi];
FArrayBox &mskv_fab = (*vec_mskv[lev])[mfi];

init_state_from_netcdf(lev, temp_fab, salt_fab,
xvel_fab, yvel_fab,
ubar_fab, vbar_fab,
mskr_fab,
msku_fab, mskv_fab,
NC_temp_fab, NC_salt_fab,
NC_xvel_fab, NC_yvel_fab,
NC_ubar_fab, NC_vbar_fab,
NC_mskr_fab,
NC_msku_fab, NC_mskv_fab);
NC_ubar_fab, NC_vbar_fab);
} // mf
} // omp
}
Expand Down Expand Up @@ -318,6 +304,56 @@ REMORA::init_coriolis_from_netcdf (int lev)
vec_fcor[lev]->FillBoundary(geom[lev].periodicity());
}

/**
* REMORA function that initializes land/sea masks from netcdf file
*
* @param lev Integer specifying the current level
*/
void
REMORA::init_masks_from_netcdf (int lev)
{
// *** FArrayBox's at this level for holding the INITIAL data
Vector<FArrayBox> NC_mskr_fab ; NC_mskr_fab.resize(num_boxes_at_level[lev]);
Vector<FArrayBox> NC_msku_fab ; NC_msku_fab.resize(num_boxes_at_level[lev]);
Vector<FArrayBox> NC_mskv_fab ; NC_mskv_fab.resize(num_boxes_at_level[lev]);
Vector<FArrayBox> NC_mskp_fab ; NC_mskp_fab.resize(num_boxes_at_level[lev]);

for (int idx = 0; idx < num_boxes_at_level[lev]; idx++)
{
read_masks_from_netcdf(lev,boxes_at_level[lev][idx], nc_grid_file[lev][idx],
NC_mskr_fab[idx],NC_msku_fab[idx],
NC_mskv_fab[idx],NC_mskp_fab[idx]);

#ifdef _OPENMP
#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
#endif
{
// Don't tile this since we are operating on full FABs in this routine
for ( MFIter mfi(*cons_new[lev], false); mfi.isValid(); ++mfi )
{
FArrayBox &mskr_fab = (*vec_mskr[lev])[mfi];
FArrayBox &msku_fab = (*vec_msku[lev])[mfi];
FArrayBox &mskv_fab = (*vec_mskv[lev])[mfi];
FArrayBox &mskp_fab = (*vec_mskp[lev])[mfi];

//
// FArrayBox to FArrayBox copy does "copy on intersection"
// This only works here because we have broadcast the FArrayBox of data from the netcdf file to all ranks
//

mskr_fab.template copy<RunOn::Device>(NC_mskr_fab[idx]);
msku_fab.template copy<RunOn::Device>(NC_msku_fab[idx]);
mskv_fab.template copy<RunOn::Device>(NC_mskv_fab[idx]);
mskp_fab.template copy<RunOn::Device>(NC_mskp_fab[idx]);
} // mf
} // omp
} // idx
vec_mskr[lev]->FillBoundary(geom[lev].periodicity());
vec_msku[lev]->FillBoundary(geom[lev].periodicity());
vec_mskv[lev]->FillBoundary(geom[lev].periodicity());
vec_mskp[lev]->FillBoundary(geom[lev].periodicity());
}

/**
* REMORA function that initializes time series of boundary data from a netcdf file
*
Expand Down Expand Up @@ -369,17 +405,12 @@ init_state_from_netcdf (int /*lev*/,
FArrayBox& temp_fab, FArrayBox& salt_fab,
FArrayBox& x_vel_fab, FArrayBox& y_vel_fab,
FArrayBox& ubar_fab, FArrayBox& vbar_fab,
FArrayBox& mskr_fab,
FArrayBox& msku_fab, FArrayBox& mskv_fab,
const Vector<FArrayBox>& NC_temp_fab,
const Vector<FArrayBox>& NC_salt_fab,
const Vector<FArrayBox>& NC_xvel_fab,
const Vector<FArrayBox>& NC_yvel_fab,
const Vector<FArrayBox>& NC_ubar_fab,
const Vector<FArrayBox>& NC_vbar_fab,
const Vector<FArrayBox>& NC_mskr_fab,
const Vector<FArrayBox>& NC_msku_fab,
const Vector<FArrayBox>& NC_mskv_fab)
const Vector<FArrayBox>& NC_vbar_fab)
{
int nboxes = NC_xvel_fab.size();
for (int idx = 0; idx < nboxes; idx++)
Expand All @@ -394,10 +425,6 @@ init_state_from_netcdf (int /*lev*/,
y_vel_fab.template copy<RunOn::Device>(NC_yvel_fab[idx]);
ubar_fab.template copy<RunOn::Device>(NC_ubar_fab[idx],0,0,1);
vbar_fab.template copy<RunOn::Device>(NC_vbar_fab[idx],0,0,1);
mskr_fab.template copy<RunOn::Device>(NC_mskr_fab[idx],0,0,1);
msku_fab.template copy<RunOn::Device>(NC_msku_fab[idx],0,0,1);
mskv_fab.template copy<RunOn::Device>(NC_mskv_fab[idx],0,0,1);

} // idx
}

Expand Down
3 changes: 3 additions & 0 deletions Source/Initialization/REMORA_make_new_level.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -336,6 +336,7 @@ void REMORA::resize_stuff(int lev)
vec_mskr.resize(lev+1);
vec_msku.resize(lev+1);
vec_mskv.resize(lev+1);
vec_mskp.resize(lev+1);
vec_sstore.resize(lev+1);

vec_pm.resize(lev+1);
Expand Down Expand Up @@ -446,6 +447,7 @@ void REMORA::init_stuff (int lev, const BoxArray& ba, const DistributionMapping&
vec_mskr[lev].reset(new MultiFab(ba2d,dm,1,IntVect(NGROW+1,NGROW+1,0)));
vec_msku[lev].reset(new MultiFab(convert(ba2d,IntVect(1,0,0)),dm,1,IntVect(NGROW,NGROW,0)));
vec_mskv[lev].reset(new MultiFab(convert(ba2d,IntVect(0,1,0)),dm,1,IntVect(NGROW,NGROW,0)));
vec_mskp[lev].reset(new MultiFab(convert(ba2d,IntVect(1,1,0)),dm,1,IntVect(NGROW,NGROW,0)));

vec_pm[lev].reset(new MultiFab(ba2d,dm,1,IntVect(NGROW+1,NGROW+2,0)));
vec_pn[lev].reset(new MultiFab(ba2d,dm,1,IntVect(NGROW+2,NGROW+1,0)));
Expand Down Expand Up @@ -477,6 +479,7 @@ void REMORA::init_stuff (int lev, const BoxArray& ba, const DistributionMapping&
vec_mskr[lev]->setVal(1.0_rt);
vec_msku[lev]->setVal(1.0_rt);
vec_mskv[lev]->setVal(1.0_rt);
vec_mskp[lev]->setVal(1.0_rt);

// Initialize these vars even if we aren't using GLS to
// avoid issues on e.g. checkpoint
Expand Down
26 changes: 26 additions & 0 deletions Source/REMORA.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -294,6 +294,9 @@ public:
/** land/sea mask at y-faces (2D) */
amrex::Vector<std::unique_ptr<amrex::MultiFab>> vec_mskv;

/** land/sea mask at cell corners (2D) */
amrex::Vector<std::unique_ptr<amrex::MultiFab>> vec_mskp;

/** map factor (2D) */
amrex::Vector<std::unique_ptr<amrex::MultiFab>> vec_pm;

Expand Down Expand Up @@ -357,6 +360,10 @@ public:
amrex::MultiFab const* mf_fcor,
amrex::MultiFab const* mf_visc2_p,
amrex::MultiFab const* mf_visc2_r,
amrex::MultiFab const* mf_mskr,
amrex::MultiFab const* mf_msku,
amrex::MultiFab const* mf_mskv,
amrex::MultiFab const* mf_mskp,
const amrex::Real dtfast_lev,
const bool predictor_2d_step,
const bool first_2d_step, int my_iif,
Expand All @@ -383,6 +390,8 @@ public:
amrex::MultiFab const* mf_h,
amrex::MultiFab const* mf_pm,
amrex::MultiFab const* mf_pn,
amrex::MultiFab const* mf_msku,
amrex::MultiFab const* mf_mskv,
const int N,
const amrex::Real dt_lev);

Expand All @@ -405,6 +414,8 @@ public:
const amrex::MultiFab* mf_svstr,
const amrex::MultiFab* mf_bustr,
const amrex::MultiFab* mf_bvstr,
const amrex::MultiFab* mf_msku,
const amrex::MultiFab* mf_mskv,
const int iic, const int nfirst,
const int nnew, int nstp, int nrhs,
int N, const amrex::Real dt_lev);
Expand All @@ -425,6 +436,8 @@ public:
const amrex::Array4<amrex::Real const>& h,
const amrex::Array4<amrex::Real const>& pm,
const amrex::Array4<amrex::Real const>& pn,
const amrex::Array4<amrex::Real const>& msku,
const amrex::Array4<amrex::Real const>& mskv,
int iic, int ntfirst, int nrhs, int N,
const amrex::Real dt_lev);

Expand All @@ -439,6 +452,8 @@ public:
const amrex::Array4<amrex::Real const>& pm,
const amrex::Array4<amrex::Real const>& W,
const amrex::Array4<amrex::Real >& FC,
const amrex::Array4<amrex::Real const>& msku,
const amrex::Array4<amrex::Real const>& mskv,
int nrhs, int nnew, int N, const amrex::Real dt_lev);

/** Calculation of the RHS */
Expand Down Expand Up @@ -482,6 +497,7 @@ public:
const amrex::Array4<amrex::Real const>& z_w,
const amrex::Array4<amrex::Real const>& z_r,
const amrex::Array4<amrex::Real const>& h,
const amrex::Array4<amrex::Real const>& mskr,
const int N);

void prsgrd (const amrex::Box& bx,
Expand All @@ -497,6 +513,8 @@ public:
const amrex::Array4<amrex::Real const>& Hz,
const amrex::Array4<amrex::Real const>& z_r,
const amrex::Array4<amrex::Real const>& z_w,
const amrex::Array4<amrex::Real const>& msku,
const amrex::Array4<amrex::Real const>& mskv,
const int nrhs, const int N);

/** Update velocities or tracers with diffusion/viscosity
Expand Down Expand Up @@ -544,6 +562,7 @@ public:
const amrex::Array4<amrex::Real const>& Dphi2,
const amrex::Array4<amrex::Real >& DC,
const amrex::Array4<amrex::Real >& FC,
const amrex::Array4<amrex::Real const>& msk,
const int nnew);

void vert_mean_3d (const amrex::Box& bx,
Expand All @@ -554,6 +573,7 @@ public:
const amrex::Array4<amrex::Real >& DC,
const amrex::Array4<amrex::Real >& CF,
const amrex::Array4<amrex::Real const>& dxlen,
const amrex::Array4<amrex::Real const>& msk,
const int nnew, const int N);

/** Harmonic viscosity */
Expand All @@ -570,6 +590,7 @@ public:
const amrex::Array4<amrex::Real const>& Hz,
const amrex::Array4<amrex::Real const>& pm,
const amrex::Array4<amrex::Real const>& pn,
const amrex::Array4<amrex::Real const>& mskp,
int nrhs, int nnew,
const amrex::Real dt_lev);

Expand All @@ -581,6 +602,8 @@ public:
const amrex::Array4<amrex::Real const>& Hz,
const amrex::Array4<amrex::Real const>& pm,
const amrex::Array4<amrex::Real const>& pn,
const amrex::Array4<amrex::Real const>& msku,
const amrex::Array4<amrex::Real const>& mskv,
const amrex::Real dt_lev,
const int ncomp);

Expand Down Expand Up @@ -611,11 +634,13 @@ public:
void init_gls_vmix (int lev, SolverChoice solver_choice);
void gls_prestep (int lev, amrex::MultiFab* mf_gls, amrex::MultiFab* mf_tke,
amrex::MultiFab& mf_W,
amrex::MultiFab* mf_msku, amrex::MultiFab* mf_mskv,
const int nstp, const int nnew, const int iic, const int ntfirst,
const int N, const amrex::Real dt_lev);
void gls_corrector (int lev, amrex::MultiFab* mf_gls, amrex::MultiFab* mf_tke,
amrex::MultiFab& mf_W, amrex::MultiFab* mf_Akv, amrex::MultiFab* mf_Akt,
amrex::MultiFab* mf_Akk, amrex::MultiFab* mf_Akp,
amrex::MultiFab* mf_msku, amrex::MultiFab* mf_mskv,
const int nstp, const int nnew,
const int N, const amrex::Real dt_lev);

Expand Down Expand Up @@ -659,6 +684,7 @@ public:

void init_data_from_netcdf (int lev);
void init_bdry_from_netcdf ();
void init_masks_from_netcdf (int lev);
void init_bathymetry_from_netcdf (int lev);
void init_zeta_from_netcdf (int lev);
void init_coriolis_from_netcdf (int lev);
Expand Down
Loading

0 comments on commit df2cfbe

Please sign in to comment.