From df2cfbe16e8a205832f03dc14f408361bddaa1af Mon Sep 17 00:00:00 2001 From: Hannah Klion Date: Tue, 30 Jul 2024 09:57:55 -0700 Subject: [PATCH 1/2] First pass at adding land/sea masking to evolution - 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 --- Source/IO/ReadFromInitNetcdf.cpp | 31 +++++-- Source/Initialization/REMORA_init.cpp | 6 +- .../REMORA_init_from_netcdf.cpp | 93 ++++++++++++------- .../Initialization/REMORA_make_new_level.cpp | 3 + Source/REMORA.H | 26 ++++++ Source/TimeIntegration/REMORA_advance_2d.cpp | 35 ++++--- .../REMORA_advance_2d_onestep.cpp | 4 + Source/TimeIntegration/REMORA_advance_3d.cpp | 28 ++++-- .../TimeIntegration/REMORA_advance_3d_ml.cpp | 16 ++++ Source/TimeIntegration/REMORA_gls.cpp | 40 ++++---- Source/TimeIntegration/REMORA_prestep.cpp | 8 +- .../REMORA_prestep_t_advection.cpp | 6 +- Source/TimeIntegration/REMORA_prsgrd.cpp | 10 +- Source/TimeIntegration/REMORA_rho_eos.cpp | 6 +- Source/TimeIntegration/REMORA_rhs_t_3d.cpp | 8 +- Source/TimeIntegration/REMORA_setup_step.cpp | 26 ++++-- Source/TimeIntegration/REMORA_t3dmix.cpp | 4 + .../REMORA_update_massflux_3d.cpp | 6 ++ Source/TimeIntegration/REMORA_uv3dmix.cpp | 6 +- .../TimeIntegration/REMORA_vert_mean_3d.cpp | 2 + 20 files changed, 268 insertions(+), 96 deletions(-) diff --git a/Source/IO/ReadFromInitNetcdf.cpp b/Source/IO/ReadFromInitNetcdf.cpp index 15af7a3..34a2113 100644 --- a/Source/IO/ReadFromInitNetcdf.cpp +++ b/Source/IO/ReadFromInitNetcdf.cpp @@ -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; @@ -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(domain, fname, NC_names, NC_dim_types, NC_fabs); @@ -91,4 +86,28 @@ read_coriolis_from_netcdf (int /*lev*/, // Read the netcdf file and fill these FABs BuildFABsFromNetCDFFile(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 NC_fabs; + Vector NC_names; + Vector 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(domain, fname, NC_names, NC_dim_types, NC_fabs); +} #endif // ROMSX_USE_NETCDF diff --git a/Source/Initialization/REMORA_init.cpp b/Source/Initialization/REMORA_init.cpp index c4fb2c8..e13288d 100644 --- a/Source/Initialization/REMORA_init.cpp +++ b/Source/Initialization/REMORA_init.cpp @@ -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 diff --git a/Source/Initialization/REMORA_init_from_netcdf.cpp b/Source/Initialization/REMORA_init_from_netcdf.cpp index 7d07753..5273646 100644 --- a/Source/Initialization/REMORA_init_from_netcdf.cpp +++ b/Source/Initialization/REMORA_init_from_netcdf.cpp @@ -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, @@ -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& NC_temp_fab, const Vector& NC_salt_fab, const Vector& NC_xvel_fab, const Vector& NC_yvel_fab, const Vector& NC_ubar_fab, - const Vector& NC_vbar_fab, - const Vector& NC_mskr_fab, - const Vector& NC_msku_fab, - const Vector& NC_mskv_fab); + const Vector& NC_vbar_fab); void read_bathymetry_from_netcdf (int lev, const Box& domain, const std::string& fname, @@ -81,18 +79,13 @@ REMORA::init_data_from_netcdf (int lev) Vector NC_yvel_fab ; NC_yvel_fab.resize(num_boxes_at_level[lev]); Vector NC_ubar_fab ; NC_ubar_fab.resize(num_boxes_at_level[lev]); Vector NC_vbar_fab ; NC_vbar_fab.resize(num_boxes_at_level[lev]); - Vector NC_mskr_fab ; NC_mskr_fab.resize(num_boxes_at_level[lev]); - Vector NC_msku_fab ; NC_msku_fab.resize(num_boxes_at_level[lev]); - Vector 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]); } @@ -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 } @@ -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 NC_mskr_fab ; NC_mskr_fab.resize(num_boxes_at_level[lev]); + Vector NC_msku_fab ; NC_msku_fab.resize(num_boxes_at_level[lev]); + Vector NC_mskv_fab ; NC_mskv_fab.resize(num_boxes_at_level[lev]); + Vector 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(NC_mskr_fab[idx]); + msku_fab.template copy(NC_msku_fab[idx]); + mskv_fab.template copy(NC_mskv_fab[idx]); + mskp_fab.template copy(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 * @@ -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& NC_temp_fab, const Vector& NC_salt_fab, const Vector& NC_xvel_fab, const Vector& NC_yvel_fab, const Vector& NC_ubar_fab, - const Vector& NC_vbar_fab, - const Vector& NC_mskr_fab, - const Vector& NC_msku_fab, - const Vector& NC_mskv_fab) + const Vector& NC_vbar_fab) { int nboxes = NC_xvel_fab.size(); for (int idx = 0; idx < nboxes; idx++) @@ -394,10 +425,6 @@ init_state_from_netcdf (int /*lev*/, y_vel_fab.template copy(NC_yvel_fab[idx]); ubar_fab.template copy(NC_ubar_fab[idx],0,0,1); vbar_fab.template copy(NC_vbar_fab[idx],0,0,1); - mskr_fab.template copy(NC_mskr_fab[idx],0,0,1); - msku_fab.template copy(NC_msku_fab[idx],0,0,1); - mskv_fab.template copy(NC_mskv_fab[idx],0,0,1); - } // idx } diff --git a/Source/Initialization/REMORA_make_new_level.cpp b/Source/Initialization/REMORA_make_new_level.cpp index 867189f..8adb3ba 100644 --- a/Source/Initialization/REMORA_make_new_level.cpp +++ b/Source/Initialization/REMORA_make_new_level.cpp @@ -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); @@ -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))); @@ -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 diff --git a/Source/REMORA.H b/Source/REMORA.H index 5da2e0d..d625eb3 100644 --- a/Source/REMORA.H +++ b/Source/REMORA.H @@ -294,6 +294,9 @@ public: /** land/sea mask at y-faces (2D) */ amrex::Vector> vec_mskv; + /** land/sea mask at cell corners (2D) */ + amrex::Vector> vec_mskp; + /** map factor (2D) */ amrex::Vector> vec_pm; @@ -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, @@ -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); @@ -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); @@ -425,6 +436,8 @@ public: const amrex::Array4& h, const amrex::Array4& pm, const amrex::Array4& pn, + const amrex::Array4& msku, + const amrex::Array4& mskv, int iic, int ntfirst, int nrhs, int N, const amrex::Real dt_lev); @@ -439,6 +452,8 @@ public: const amrex::Array4& pm, const amrex::Array4& W, const amrex::Array4& FC, + const amrex::Array4& msku, + const amrex::Array4& mskv, int nrhs, int nnew, int N, const amrex::Real dt_lev); /** Calculation of the RHS */ @@ -482,6 +497,7 @@ public: const amrex::Array4& z_w, const amrex::Array4& z_r, const amrex::Array4& h, + const amrex::Array4& mskr, const int N); void prsgrd (const amrex::Box& bx, @@ -497,6 +513,8 @@ public: const amrex::Array4& Hz, const amrex::Array4& z_r, const amrex::Array4& z_w, + const amrex::Array4& msku, + const amrex::Array4& mskv, const int nrhs, const int N); /** Update velocities or tracers with diffusion/viscosity @@ -544,6 +562,7 @@ public: const amrex::Array4& Dphi2, const amrex::Array4& DC, const amrex::Array4& FC, + const amrex::Array4& msk, const int nnew); void vert_mean_3d (const amrex::Box& bx, @@ -554,6 +573,7 @@ public: const amrex::Array4& DC, const amrex::Array4& CF, const amrex::Array4& dxlen, + const amrex::Array4& msk, const int nnew, const int N); /** Harmonic viscosity */ @@ -570,6 +590,7 @@ public: const amrex::Array4& Hz, const amrex::Array4& pm, const amrex::Array4& pn, + const amrex::Array4& mskp, int nrhs, int nnew, const amrex::Real dt_lev); @@ -581,6 +602,8 @@ public: const amrex::Array4& Hz, const amrex::Array4& pm, const amrex::Array4& pn, + const amrex::Array4& msku, + const amrex::Array4& mskv, const amrex::Real dt_lev, const int ncomp); @@ -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); @@ -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); diff --git a/Source/TimeIntegration/REMORA_advance_2d.cpp b/Source/TimeIntegration/REMORA_advance_2d.cpp index 6df41e8..c5e39f7 100644 --- a/Source/TimeIntegration/REMORA_advance_2d.cpp +++ b/Source/TimeIntegration/REMORA_advance_2d.cpp @@ -26,6 +26,10 @@ using namespace amrex; * @param[in ] mf_h * @param[inout] mf_visc2_p * @param[inout] mf_visc2_f + * @param[in ] mf_mskr + * @param[in ] mf_msku + * @param[in ] mf_mskv + * @param[in ] mf_mskp * @param[in ] dtfast_lev * @param[in ] predictor_2d_step * @param[in ] first_2d_step @@ -58,6 +62,10 @@ REMORA::advance_2d (int lev, MultiFab const* mf_fcor, MultiFab const* mf_visc2_p, MultiFab const* mf_visc2_r, + MultiFab const* mf_mskr, + MultiFab const* mf_msku, + MultiFab const* mf_mskv, + MultiFab const* mf_mskp, Real dtfast_lev, bool predictor_2d_step, bool first_2d_step, int my_iif, @@ -190,6 +198,11 @@ REMORA::advance_2d (int lev, Array4 const& pn = mf_pn->const_array(mfi); Array4 const& fcor = mf_fcor->const_array(mfi); + Array4 const& mskr = mf_mskr->const_array(mfi); + Array4 const& msku = mf_msku->const_array(mfi); + Array4 const& mskv = mf_mskv->const_array(mfi); + Array4 const& mskp = mf_mskp->const_array(mfi); + Box bx = mfi.tilebox(); Box gbx = mfi.growntilebox(); Box gbx1 = mfi.growntilebox(IntVect(NGROW-1,NGROW-1,0)); @@ -372,7 +385,7 @@ REMORA::advance_2d (int lev, { rhs_zeta(i,j,0) = (DUon(i,j,0)-DUon(i+1,j,0))+ (DVom(i,j,0)-DVom(i,j+1,0)); - zeta_new(i,j,0) = zeta(i,j,0,kstp)+ pm(i,j,0)*pn(i,j,0)*cff1*rhs_zeta(i,j,0); + zeta_new(i,j,0) = (zeta(i,j,0,kstp)+ pm(i,j,0)*pn(i,j,0)*cff1*rhs_zeta(i,j,0)) * mskr(i,j,0); Dnew(i,j,0) = zeta_new(i,j,0)+h(i,j,0); //Pressure gradient terms: @@ -392,8 +405,8 @@ REMORA::advance_2d (int lev, { rhs_zeta(i,j,0)=(DUon(i,j,0)-DUon(i+1,j,0))+ (DVom(i,j,0)-DVom(i,j+1,0)); - zeta_new(i,j,0)=zeta(i,j,0,kstp)+ - pm(i,j,0)*pn(i,j,0)*cff1*rhs_zeta(i,j,0); + zeta_new(i,j,0)=(zeta(i,j,0,kstp)+ + pm(i,j,0)*pn(i,j,0)*cff1*rhs_zeta(i,j,0)) * mskr(i,j,0); Dnew(i,j,0)=zeta_new(i,j,0)+h(i,j,0); //Pressure gradient terms zwrk(i,j,0)=cff5*zeta(i,j,0,krhs)+ @@ -419,6 +432,7 @@ REMORA::advance_2d (int lev, pm(i,j,0)*pn(i,j,0)*(cff+ cff2*rzeta(i,j,0,kstp)- cff3*rzeta(i,j,0,ptsk)); + zeta_new(i,j,0) *= mskr(i,j,0); Dnew(i,j,0)=zeta_new(i,j,0)+h(i,j,0); //Pressure gradient terms zwrk(i,j,0)=cff5*zeta_new(i,j,0)+cff4*zeta(i,j,0,krhs); @@ -515,14 +529,13 @@ REMORA::advance_2d (int lev, // coriolis(xbxD, ybxD, ubar_const, vbar_const, rhs_ubar, rhs_vbar, Drhs, fomn, krhs, 0); } - //----------------------------------------------------------------------- //Add in horizontal harmonic viscosity. // Consider generalizing or copying uv3dmix, where Drhs is used instead of Hz and u=>ubar v=>vbar, drop dt terms //----------------------------------------------------------------------- uv3dmix(xbxD, ybxD, ubar, vbar, ubar, vbar, rhs_ubar, rhs_vbar, visc2_p, visc2_r, Drhs_const, - pm, pn, krhs, nnew, 0.0_rt); + pm, pn, mskp, krhs, nnew, 0.0_rt); //----------------------------------------------------------------------- // Coupling from 3d to 2d @@ -636,7 +649,7 @@ REMORA::advance_2d (int lev, Real Dnew_avg =1.0_rt/(Dnew(i,j,0)+Dnew(i-1,j,0)); ubar(i,j,0,knew)=(ubar(i,j,0,kstp)* (Dstp(i,j,0)+Dstp(i-1,j,0))+ - cff*cff1*rhs_ubar(i,j,0))*Dnew_avg; + cff*cff1*rhs_ubar(i,j,0))*Dnew_avg * msku(i,j,0); }); ParallelFor(ybxD, [=] AMREX_GPU_DEVICE (int i, int j, int ) @@ -645,7 +658,7 @@ REMORA::advance_2d (int lev, Real Dnew_avg=1.0_rt/(Dnew(i,j,0)+Dnew(i,j-1,0)); vbar(i,j,0,knew)=(vbar(i,j,0,kstp)* (Dstp(i,j,0)+Dstp(i,j-1,0))+ - cff*cff1*rhs_vbar(i,j,0))*Dnew_avg; + cff*cff1*rhs_vbar(i,j,0))*Dnew_avg * mskv(i,j,0); }); } else if (predictor_2d_step) { @@ -658,7 +671,7 @@ REMORA::advance_2d (int lev, Real Dnew_avg=1.0_rt/(Dnew(i,j,0)+Dnew(i-1,j,0)); ubar(i,j,0,knew)=(ubar(i,j,0,kstp)* (Dstp(i,j,0)+Dstp(i-1,j,0))+ - cff*cff1*rhs_ubar(i,j,0))*Dnew_avg; + cff*cff1*rhs_ubar(i,j,0))*Dnew_avg * msku(i,j,0); }); ParallelFor(ybxD, [=] AMREX_GPU_DEVICE (int i, int j, int ) @@ -667,7 +680,7 @@ REMORA::advance_2d (int lev, Real Dnew_avg=1.0_rt/(Dnew(i,j,0)+Dnew(i,j-1,0)); vbar(i,j,0,knew)=(vbar(i,j,0,kstp)* (Dstp(i,j,0)+Dstp(i,j-1,0))+ - cff*cff1*rhs_vbar(i,j,0))*Dnew_avg; + cff*cff1*rhs_vbar(i,j,0))*Dnew_avg * mskv(i,j,0); }); } else if ((!predictor_2d_step)) { @@ -684,7 +697,7 @@ REMORA::advance_2d (int lev, (Dstp(i,j,0)+Dstp(i-1,j,0))+ cff*(cff1*rhs_ubar(i,j,0)+ cff2*rubar(i,j,0,kstp)- - cff3*rubar(i,j,0,ptsk)))*Dnew_avg; + cff3*rubar(i,j,0,ptsk)))*Dnew_avg * msku(i,j,0); }); ParallelFor(ybxD, [=] AMREX_GPU_DEVICE (int i, int j, int ) @@ -695,7 +708,7 @@ REMORA::advance_2d (int lev, (Dstp(i,j,0)+Dstp(i,j-1,0))+ cff*(cff1*rhs_vbar(i,j,0)+ cff2*rvbar(i,j,0,kstp)- - cff3*rvbar(i,j,0,ptsk)))*Dnew_avg; + cff3*rvbar(i,j,0,ptsk)))*Dnew_avg * mskv(i,j,0); }); } diff --git a/Source/TimeIntegration/REMORA_advance_2d_onestep.cpp b/Source/TimeIntegration/REMORA_advance_2d_onestep.cpp index a521544..647802f 100644 --- a/Source/TimeIntegration/REMORA_advance_2d_onestep.cpp +++ b/Source/TimeIntegration/REMORA_advance_2d_onestep.cpp @@ -23,6 +23,8 @@ void REMORA::advance_2d_onestep (int lev, Real /*dt_lev*/, Real dtfast_lev, int vec_pm[lev].get(), vec_pn[lev].get(), vec_fcor[lev].get(), vec_visc2_p[lev].get(), vec_visc2_r[lev].get(), + vec_mskr[lev].get(), vec_msku[lev].get(), vec_mskv[lev].get(), + vec_mskp[lev].get(), dtfast_lev, predictor_2d_step, first_2d_step, my_iif, next_indx1); //Corrector. Skip it on last fast step @@ -44,6 +46,8 @@ void REMORA::advance_2d_onestep (int lev, Real /*dt_lev*/, Real dtfast_lev, int vec_pm[lev].get(), vec_pn[lev].get(), vec_fcor[lev].get(), vec_visc2_p[lev].get(), vec_visc2_r[lev].get(), + vec_mskr[lev].get(), vec_msku[lev].get(), vec_mskv[lev].get(), + vec_mskp[lev].get(), dtfast_lev, predictor_2d_step, first_2d_step, my_iif, next_indx1); } } diff --git a/Source/TimeIntegration/REMORA_advance_3d.cpp b/Source/TimeIntegration/REMORA_advance_3d.cpp index 69b95bf..3c76fe0 100644 --- a/Source/TimeIntegration/REMORA_advance_3d.cpp +++ b/Source/TimeIntegration/REMORA_advance_3d.cpp @@ -26,6 +26,8 @@ REMORA::advance_3d (int lev, MultiFab& mf_cons, MultiFab const* mf_h, MultiFab const* mf_pm, MultiFab const* mf_pn, + MultiFab const* mf_msku, + MultiFab const* mf_mskv, const int N, Real dt_lev) { const int nrhs = 0; @@ -70,6 +72,9 @@ REMORA::advance_3d (int lev, MultiFab& mf_cons, Array4 const& pm = mf_pm->const_array(mfi); Array4 const& pn = mf_pn->const_array(mfi); + Array4 const& msku = mf_msku->const_array(mfi); + Array4 const& mskv = mf_mskv->const_array(mfi); + Box bx = mfi.tilebox(); Box gbx2 = mfi.growntilebox(IntVect(NGROW,NGROW,0)); Box gbx21 = mfi.growntilebox(IntVect(NGROW,NGROW,NGROW-1)); @@ -136,13 +141,13 @@ REMORA::advance_3d (int lev, MultiFab& mf_cons, mf_DC[mfi].template setVal(0.,xbx); fab_CF.template setVal(0.,xbx); - vert_mean_3d(xbx,1,0,u,Hz,DU_avg1,DC,CF,pn,nnew,N); + vert_mean_3d(xbx,1,0,u,Hz,DU_avg1,DC,CF,pn,msku,nnew,N); // Reset to zero on the box on which they'll be used mf_DC[mfi].template setVal(0.,ybx); fab_CF.template setVal(0.,ybx); - vert_mean_3d(ybx,0,1,v,Hz,DV_avg1,DC,CF,pm,nnew,N); + vert_mean_3d(ybx,0,1,v,Hz,DV_avg1,DC,CF,pm,mskv,nnew,N); } // Apply physical boundary conditions to u and v @@ -179,8 +184,10 @@ REMORA::advance_3d (int lev, MultiFab& mf_cons, Array4 const& Huon = mf_Huon->array(mfi); Array4 const& Hvom = mf_Hvom->array(mfi); - Array4 const& pm = mf_pm->const_array(mfi); - Array4 const& pn = mf_pn->const_array(mfi); + Array4 const& pm = mf_pm->const_array(mfi); + Array4 const& pn = mf_pn->const_array(mfi); + Array4 const& msku = mf_msku->const_array(mfi); + Array4 const& mskv = mf_mskv->const_array(mfi); Box gbx2 = mfi.growntilebox(IntVect(NGROW,NGROW,0)); @@ -203,12 +210,12 @@ REMORA::advance_3d (int lev, MultiFab& mf_cons, // Reset to zero on the box on which they'll be used fab_FC.template setVal(0.,gbx2); mf_DC[mfi].template setVal(0.,grow(gbx2,IntVect(0,0,1))); - update_massflux_3d(gbx2,1,0,u,ubar,Huon,Hz,pn,DU_avg1,DU_avg2,DC,FC,nnew); + update_massflux_3d(gbx2,1,0,u,ubar,Huon,Hz,pn,DU_avg1,DU_avg2,DC,FC,msku,nnew); // Reset to zero on the box on which they'll be used fab_FC.template setVal(0.,gbx2); mf_DC[mfi].template setVal(0.,grow(gbx2,IntVect(0,0,1))); - update_massflux_3d(gbx2,0,1,v,vbar,Hvom,Hz,pm,DV_avg1,DV_avg2,DC,FC,nnew); + update_massflux_3d(gbx2,0,1,v,vbar,Hvom,Hz,pm,DV_avg1,DV_avg2,DC,FC,mskv,nnew); #endif } @@ -296,7 +303,9 @@ REMORA::advance_3d (int lev, MultiFab& mf_cons, nnew = 1-nstp; if (solverChoice.vert_mixing_type == VertMixingType::GLS) { gls_corrector(lev, vec_gls[lev].get(), vec_tke[lev].get(), mf_W, vec_Akv[lev].get(), - vec_Akt[lev].get(),vec_Akk[lev].get(), vec_Akp[lev].get(), nstp, nnew, N, dt_lev); + vec_Akt[lev].get(),vec_Akk[lev].get(), vec_Akp[lev].get(), + vec_msku[lev].get(), vec_mskv[lev].get(), + nstp, nnew, N, dt_lev); } nnew = 0; @@ -311,6 +320,8 @@ REMORA::advance_3d (int lev, MultiFab& mf_cons, Array4 const& h = mf_h->const_array(mfi); Array4 const& pm = mf_pm->const_array(mfi); Array4 const& pn = mf_pn->const_array(mfi); + Array4 const& msku = mf_msku->const_array(mfi); + Array4 const& mskv = mf_mskv->const_array(mfi); Box bx = mfi.tilebox(); Box gbx = mfi.growntilebox(); @@ -340,7 +351,7 @@ REMORA::advance_3d (int lev, MultiFab& mf_cons, { Array4 const& sstore = mf_sstore->array(mfi, i_comp); rhs_t_3d(bx, gbx, mf_cons.array(mfi,i_comp), sstore, Huon, Hvom, - Hz, pn, pm, W, FC, nrhs, nnew, N,dt_lev); + Hz, pn, pm, W, FC, msku, mskv, nrhs, nnew, N,dt_lev); } } // mfi @@ -354,6 +365,7 @@ REMORA::advance_3d (int lev, MultiFab& mf_cons, Array4 const& Hzk = mf_Hzk.array(mfi); Array4 const& Hz = mf_Hz->const_array(mfi); + Array4 const& cons = mf_cons.array(mfi); Box bx = mfi.tilebox(); diff --git a/Source/TimeIntegration/REMORA_advance_3d_ml.cpp b/Source/TimeIntegration/REMORA_advance_3d_ml.cpp index 7f0d833..0eb132d 100644 --- a/Source/TimeIntegration/REMORA_advance_3d_ml.cpp +++ b/Source/TimeIntegration/REMORA_advance_3d_ml.cpp @@ -25,6 +25,7 @@ void REMORA::advance_3d_ml (int lev, Real dt_lev) vec_Akv[lev], vec_Akt[lev], vec_Hz[lev], vec_Huon[lev], vec_Hvom[lev], vec_z_w[lev], vec_hOfTheConfusingName[lev].get(), vec_pm[lev].get(), vec_pn[lev].get(), + vec_msku[lev].get(), vec_mskv[lev].get(), N, dt_lev); // Ideally can roll all of these into a single call @@ -36,6 +37,7 @@ void REMORA::advance_3d_ml (int lev, Real dt_lev) FillPatchNoBC(lev, t_old[lev], *vec_vbar[lev], GetVecOfPtrs(vec_vbar), BdyVars::vbar,2,false,false); FillPatch(lev, t_old[lev], *vec_sstore[lev], GetVecOfPtrs(vec_sstore), BdyVars::t); + // Fill in three ways: 1) interpolate from coarse grid if lev > 0; 2) fill from physical boundaries; // 3) fine-fine fill of ghost cells with FillBoundary call // Note that we need the fine-fine and physical bc's in order to correctly move the particles @@ -44,6 +46,20 @@ void REMORA::advance_3d_ml (int lev, Real dt_lev) yvel_new[lev]->FillBoundary(geom[lev].periodicity()); FillPatch(lev, t_old[lev], *zvel_new[lev], zvel_new, BdyVars::null); + // Apply land/sea mask to tracers + for ( MFIter mfi(*cons_new[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi ) + { + Array4 const& cons = cons_new[lev]->array(mfi); + Array4 const& mskr = vec_mskr[lev]->array(mfi); + + Box bx = mfi.tilebox(); + + ParallelFor(bx, NCONS, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) + { + cons(i,j,k,n) *= mskr(i,j,0); + }); + } + #ifdef REMORA_USE_PARTICLES // ************************************************************************************** // Update the particle positions diff --git a/Source/TimeIntegration/REMORA_gls.cpp b/Source/TimeIntegration/REMORA_gls.cpp index a510f34..215ba28 100644 --- a/Source/TimeIntegration/REMORA_gls.cpp +++ b/Source/TimeIntegration/REMORA_gls.cpp @@ -4,7 +4,8 @@ using namespace amrex; void REMORA::gls_prestep (int lev, MultiFab* mf_gls, MultiFab* mf_tke, - MultiFab& mf_W, const int nstp, const int nnew, + MultiFab& mf_W, MultiFab* mf_msku, MultiFab* mf_mskv, + const int nstp, const int nnew, const int iic, const int ntfirst, const int N, const Real dt_lev) { // temps: grad, gradL, XF, FX, FXL, EF, FE, FEL @@ -18,6 +19,8 @@ REMORA::gls_prestep (int lev, MultiFab* mf_gls, MultiFab* mf_tke, Array4 const& Hz = vec_Hz[lev]->array(mfi); Array4 const& pm = vec_pm[lev]->array(mfi); Array4 const& pn = vec_pn[lev]->array(mfi); + Array4 const& msku = vec_msku[lev]->array(mfi); + Array4 const& mskv = vec_mskv[lev]->array(mfi); Box bx = mfi.tilebox(); Box xbx = surroundingNodes(bx,0); @@ -66,11 +69,11 @@ REMORA::gls_prestep (int lev, MultiFab* mf_gls, MultiFab* mf_tke, // need XF/FX/FXL from [xlo to xhi] by [ylo to yhi ] on u points ParallelFor(grow(xbx,IntVect(0,0,-1)), [=] AMREX_GPU_DEVICE (int i, int j, int k) { - Real grad_im1 = tke(i-1,j,k,nstp) - tke(i-2,j,k,nstp); - Real grad_ip1 = tke(i+1,j,k,nstp) - tke(i ,j,k,nstp); + Real grad_im1 = (tke(i-1,j,k,nstp) - tke(i-2,j,k,nstp)) * msku(i-1,j,0); + Real grad_ip1 = (tke(i+1,j,k,nstp) - tke(i ,j,k,nstp)) * msku(i+1,j,0); - Real gradL_im1 = gls(i-1,j,k,nstp) - gls(i-2,j,k,nstp); - Real gradL_ip1 = gls(i+1,j,k,nstp) - gls(i ,j,k,nstp); + Real gradL_im1 = (gls(i-1,j,k,nstp) - gls(i-2,j,k,nstp)) * msku(i-1,j,0); + Real gradL_ip1 = (gls(i+1,j,k,nstp) - gls(i ,j,k,nstp)) * msku(i+1,j,0); // Adjust boundaries // TODO: Make sure indices match with what ROMS does @@ -93,11 +96,11 @@ REMORA::gls_prestep (int lev, MultiFab* mf_gls, MultiFab* mf_tke, // need EF/FE/FEL from [xlo to xhi ] by [ylo to yhi+1] ParallelFor(grow(ybx,IntVect(0,0,-1)), [=] AMREX_GPU_DEVICE (int i, int j, int k) { - Real grad_jm1 = tke(i,j-1,k,nstp) - tke(i,j-2,k,nstp); - Real grad_jp1 = tke(i,j+1,k,nstp) - tke(i,j ,k,nstp); + Real grad_jm1 = (tke(i,j-1,k,nstp) - tke(i,j-2,k,nstp)) * mskv(i,j-1,0); + Real grad_jp1 = (tke(i,j+1,k,nstp) - tke(i,j ,k,nstp)) * mskv(i,j+1,0); - Real gradL_jm1 = gls(i,j-1,k,nstp) - gls(i,j-2,k,nstp); - Real gradL_jp1 = gls(i,j+1,k,nstp) - gls(i,j ,k,nstp); + Real gradL_jm1 = (gls(i,j-1,k,nstp) - gls(i,j-2,k,nstp)) * mskv(i,j-1,0); + Real gradL_jp1 = (gls(i,j+1,k,nstp) - gls(i,j ,k,nstp)) * mskv(i,j+1,0); // Adjust boundaries // TODO: Make sure indices match with what ROMS does @@ -212,6 +215,7 @@ void REMORA::gls_corrector (int lev, MultiFab* mf_gls, MultiFab* mf_tke, MultiFab& mf_W, MultiFab* mf_Akv, MultiFab* mf_Akt, MultiFab* mf_Akk, MultiFab* mf_Akp, + MultiFab* mf_msku, MultiFab* mf_mskv, const int nstp, const int nnew, const int N, const Real dt_lev) { @@ -462,6 +466,8 @@ REMORA::gls_corrector (int lev, MultiFab* mf_gls, MultiFab* mf_tke, Array4 const& svstr = vec_svstr[lev]->array(mfi); Array4 const& bustr = vec_bustr[lev]->array(mfi); Array4 const& bvstr = vec_bvstr[lev]->array(mfi); + Array4 const& msku = vec_msku[lev]->array(mfi); + Array4 const& mskv = vec_mskv[lev]->array(mfi); FArrayBox fab_tmp_buoy(bx_growloxy,1, amrex::The_Async_Arena()); fab_tmp_buoy.template setVal(0.); FArrayBox fab_tmp_shear(bx_growloxy,1, amrex::The_Async_Arena()); fab_tmp_shear.template setVal(0.); @@ -531,10 +537,10 @@ REMORA::gls_corrector (int lev, MultiFab* mf_gls, MultiFab* mf_tke, gradP = gls(i ,j,k,2)-gls(i-1,j,k,2); gradP_ip1 = gradP; } else { - gradK = tke(i ,j,k,2)-tke(i-1,j,k,2); - gradK_ip1 = tke(i+1,j,k,2)-tke(i ,j,k,2); - gradP = gls(i ,j,k,2)-gls(i-1,j,k,2); - gradP_ip1 = gls(i+1,j,k,2)-gls(i ,j,k,2); + gradK = (tke(i ,j,k,2)-tke(i-1,j,k,2)) * msku(i ,j,0); + gradK_ip1 = (tke(i+1,j,k,2)-tke(i ,j,k,2)) * msku(i+1,j,0); + gradP = (gls(i ,j,k,2)-gls(i-1,j,k,2)) * msku(i ,j,0); + gradP_ip1 = (gls(i+1,j,k,2)-gls(i ,j,k,2)) * msku(i+1,j,0); } curvK(i,j,k) = gradK_ip1 - gradK; @@ -553,10 +559,10 @@ REMORA::gls_corrector (int lev, MultiFab* mf_gls, MultiFab* mf_tke, //Time step advective terms ParallelFor(growLo(grow(ybx,IntVect(0,0,-1)),1,1), [=] AMREX_GPU_DEVICE (int i, int j, int k) { - Real gradK = tke(i,j ,k,2)-tke(i,j-1,k,2); - Real gradK_jp1 = tke(i,j+1,k,2)-tke(i,j ,k,2); - Real gradP = gls(i,j ,k,2)-gls(i,j-1,k,2); - Real gradP_jp1 = gls(i,j+1,k,2)-gls(i,j ,k,2); + Real gradK = (tke(i,j ,k,2)-tke(i,j-1,k,2)) * mskv(i,j ,0); + Real gradK_jp1 = (tke(i,j+1,k,2)-tke(i,j ,k,2)) * mskv(i,j+1,0); + Real gradP = (gls(i,j ,k,2)-gls(i,j-1,k,2)) * mskv(i,j ,0); + Real gradP_jp1 = (gls(i,j+1,k,2)-gls(i,j ,k,2)) * mskv(i,j+1,0); if (j == dlo.y-1 && (bcr_y.lo(1) == REMORABCType::ext_dir or ic_bc_type==IC_BC_Type::Real)) { gradK = gradK_jp1; diff --git a/Source/TimeIntegration/REMORA_prestep.cpp b/Source/TimeIntegration/REMORA_prestep.cpp index 9504b30..791fb81 100644 --- a/Source/TimeIntegration/REMORA_prestep.cpp +++ b/Source/TimeIntegration/REMORA_prestep.cpp @@ -21,6 +21,8 @@ using namespace amrex; * @param[in ] mf_svstr * @param[in ] mf_bustr * @param[in ] mf_bvstr + * @param[in ] mf_msku + * @param[in ] mf_mskv * @param[in ] iic * @param[in ] ntfirst * @param[in ] nnew @@ -47,6 +49,8 @@ REMORA::prestep (int lev, const MultiFab* mf_svstr, const MultiFab* mf_bustr, const MultiFab* mf_bvstr, + const MultiFab* mf_msku, + const MultiFab* mf_mskv, const int iic, const int ntfirst, const int nnew, int nstp, int nrhs, int N, const Real dt_lev) @@ -91,6 +95,8 @@ REMORA::prestep (int lev, Array4 const& svstr = mf_svstr->const_array(mfi); Array4 const& bustr = mf_bustr->const_array(mfi); Array4 const& bvstr = mf_bvstr->const_array(mfi); + Array4 const& msku = mf_msku->const_array(mfi); + Array4 const& mskv = mf_mskv->const_array(mfi); Real lambda = 1.0_rt; @@ -152,7 +158,7 @@ REMORA::prestep (int lev, Array4 const& sstore = (vec_sstore[lev])->array(mfi,i_comp); prestep_t_advection(bx, gbx, S_old.array(mfi,i_comp), mf_scalarcache.array(mfi,i_comp), Hz, Huon, Hvom, - W, DC, FC, sstore, z_w, h, pm, pn, iic, ntfirst, + W, DC, FC, sstore, z_w, h, pm, pn, msku, mskv, iic, ntfirst, nrhs, N, dt_lev); } diff --git a/Source/TimeIntegration/REMORA_prestep_t_advection.cpp b/Source/TimeIntegration/REMORA_prestep_t_advection.cpp index 067a0ef..4126b89 100644 --- a/Source/TimeIntegration/REMORA_prestep_t_advection.cpp +++ b/Source/TimeIntegration/REMORA_prestep_t_advection.cpp @@ -21,6 +21,8 @@ REMORA::prestep_t_advection (const Box& tbx, const Box& gbx, const Array4& h, const Array4& pm, const Array4& pn, + const Array4& msku, + const Array4& mskv, int iic, int ntfirst, int nrhs, int N, Real dt_lev) { @@ -152,13 +154,13 @@ REMORA::prestep_t_advection (const Box& tbx, const Box& gbx, ParallelFor(utbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k) { //should be t index 3 - FX(i,j,k)=tempold(i,j,k,nrhs)-tempold(i-1,j,k,nrhs); + FX(i,j,k)=(tempold(i,j,k,nrhs)-tempold(i-1,j,k,nrhs)) * msku(i,j,0); }); ParallelFor(vtbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k) { //should be t index 3 - FE(i,j,k)=tempold(i,j,k,nrhs)-tempold(i,j-1,k,nrhs); + FE(i,j,k)=(tempold(i,j,k,nrhs)-tempold(i,j-1,k,nrhs)) * mskv(i,j,0); }); Real cffa=1.0_rt/6.0_rt; diff --git a/Source/TimeIntegration/REMORA_prsgrd.cpp b/Source/TimeIntegration/REMORA_prsgrd.cpp index 7d941bd..63fb6b2 100644 --- a/Source/TimeIntegration/REMORA_prsgrd.cpp +++ b/Source/TimeIntegration/REMORA_prsgrd.cpp @@ -14,6 +14,8 @@ REMORA::prsgrd (const Box& phi_bx, const Box& phi_gbx, const Array4& Hz, const Array4& z_r, const Array4& z_w, + const Array4& msku, + const Array4& mskv, const int nrhs, const int N) { auto phi_bxD=phi_bx; @@ -101,8 +103,8 @@ REMORA::prsgrd (const Box& phi_bx, const Box& phi_gbx, ParallelFor(phi_ubx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - FC(i,j,k)=rho(i,j,k)-rho(i-1,j,k); - aux(i,j,k)=z_r(i,j,k)-z_r(i-1,j,k); + FC(i,j,k)=(rho(i,j,k)-rho(i-1,j,k)) * msku(i,j,0); + aux(i,j,k)=(z_r(i,j,k)-z_r(i-1,j,k)) * msku(i,j,0); }); //This should be nodal aux and FC need wider boxes above @@ -149,8 +151,8 @@ REMORA::prsgrd (const Box& phi_bx, const Box& phi_gbx, //This should be nodal ParallelFor(phi_vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - FC(i,j,k)= rho(i,j,k)-rho(i,j-1,k); - aux(i,j,k)= z_r(i,j,k)-z_r(i,j-1,k); + FC(i,j,k)= (rho(i,j,k)-rho(i,j-1,k)) * mskv(i,j,0); + aux(i,j,k)= (z_r(i,j,k)-z_r(i,j-1,k)) * mskv(i,j,0); }); //This should be nodal aux and FC need wider boxes above diff --git a/Source/TimeIntegration/REMORA_rho_eos.cpp b/Source/TimeIntegration/REMORA_rho_eos.cpp index 0c197c4..13512a7 100644 --- a/Source/TimeIntegration/REMORA_rho_eos.cpp +++ b/Source/TimeIntegration/REMORA_rho_eos.cpp @@ -14,6 +14,7 @@ using namespace amrex; * @param[in ] z_w * @param[in ] z_r * @param[in ] h + * @param[in ] mskr * @param[in ] N */ @@ -28,6 +29,7 @@ REMORA::rho_eos (const Box& bx, const Array4& z_w, const Array4& z_r, const Array4& h, + const Array4& mskr, const int N) { // @@ -51,9 +53,9 @@ REMORA::rho_eos (const Box& bx, ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - rho(i,j,k) = R0 - R0*Tcoef*(state(i,j,k,Temp_comp)-T0) + rho(i,j,k) = (R0 - R0*Tcoef*(state(i,j,k,Temp_comp)-T0) + R0*Scoef*(state(i,j,k,Salt_comp)-S0) - - 1000.0_rt; + - 1000.0_rt) * mskr(i,j,0); }); // diff --git a/Source/TimeIntegration/REMORA_rhs_t_3d.cpp b/Source/TimeIntegration/REMORA_rhs_t_3d.cpp index 282e095..da321c0 100644 --- a/Source/TimeIntegration/REMORA_rhs_t_3d.cpp +++ b/Source/TimeIntegration/REMORA_rhs_t_3d.cpp @@ -15,6 +15,8 @@ using namespace amrex; * @param[in ] pm * @param[in ] W * @param[inout] FC + * @param[in ] msku + * @param[in ] mskv * @param[in ] nrhs * @param[in ] nnew * @param[in ] N @@ -32,6 +34,8 @@ REMORA::rhs_t_3d (const Box& bx, const Box& gbx, const Array4& pm, const Array4& W , const Array4& FC, + const Array4& msku, + const Array4& mskv, int nrhs, int nnew, int N, Real dt_lev) { //copy the tilebox @@ -75,7 +79,7 @@ REMORA::rhs_t_3d (const Box& bx, const Box& gbx, ParallelFor(utbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k) { //should be t index 3 - FX(i,j,k)=sstore(i,j,k,nrhs)-sstore(i-1,j,k,nrhs); + FX(i,j,k)=(sstore(i,j,k,nrhs)-sstore(i-1,j,k,nrhs)) * msku(i,j,0); }); Real cffa=1.0_rt/6.0_rt; @@ -158,7 +162,7 @@ REMORA::rhs_t_3d (const Box& bx, const Box& gbx, ParallelFor(vtbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k) { //should be t index 3 - FE(i,j,k)=sstore(i,j,k,nrhs)-sstore(i,j-1,k,nrhs); + FE(i,j,k)=(sstore(i,j,k,nrhs)-sstore(i,j-1,k,nrhs)) * mskv(i,j,0); }); cffa=1.0_rt/6.0_rt; diff --git a/Source/TimeIntegration/REMORA_setup_step.cpp b/Source/TimeIntegration/REMORA_setup_step.cpp index 26175c9..696ff72 100644 --- a/Source/TimeIntegration/REMORA_setup_step.cpp +++ b/Source/TimeIntegration/REMORA_setup_step.cpp @@ -74,6 +74,11 @@ REMORA::setup_step (int lev, Real time, Real dt_lev) std::unique_ptr& mf_bustr = vec_bustr[lev]; std::unique_ptr& mf_bvstr = vec_bvstr[lev]; + std::unique_ptr& mf_mskr = vec_mskr[lev]; + std::unique_ptr& mf_msku = vec_msku[lev]; + std::unique_ptr& mf_mskv = vec_mskv[lev]; + std::unique_ptr& mf_mskp = vec_mskp[lev]; + MultiFab mf_rw(ba,dm,1,IntVect(NGROW,NGROW,0)); std::unique_ptr& mf_visc2_p = vec_visc2_p[lev]; @@ -129,6 +134,8 @@ REMORA::setup_step (int lev, Real time, Real dt_lev) Array4 const& pm = mf_pm->const_array(mfi); Array4 const& pn = mf_pn->const_array(mfi); + Array4 const& mskr = mf_mskr->const_array(mfi); + Box bx = mfi.tilebox(); Box gbx1 = mfi.growntilebox(IntVect(NGROW-1,NGROW-1,0)); Box gbx2 = mfi.growntilebox(IntVect(NGROW,NGROW,0)); @@ -171,7 +178,7 @@ REMORA::setup_step (int lev, Real time, Real dt_lev) }); Array4 const& state_old = S_old.const_array(mfi); - rho_eos(gbx2,state_old,rho,rhoA,rhoS,bvf,Hz,z_w,z_r,h,N); + rho_eos(gbx2,state_old,rho,rhoA,rhoS,bvf,Hz,z_w,z_r,h,mskr,N); } if (solverChoice.vert_mixing_type == VertMixingType::analytical) { @@ -184,6 +191,7 @@ REMORA::setup_step (int lev, Real time, Real dt_lev) MultiFab mf_W(convert(ba,IntVect(0,0,1)),dm,1,IntVect(NGROW+1,NGROW+1,0)); mf_W.setVal(0.0_rt); + if (solverChoice.use_prestep) { const int nnew = 0; prestep(lev, U_old, V_old, U_new, V_new, @@ -191,6 +199,7 @@ REMORA::setup_step (int lev, Real time, Real dt_lev) S_old, S_new, mf_W, mf_DC, mf_z_r, mf_z_w, mf_h, mf_pm, mf_pn, mf_sustr.get(), mf_svstr.get(), mf_bustr.get(), mf_bvstr.get(), + mf_msku.get(), mf_mskv.get(), iic, ntfirst, nnew, nstp, nrhs, N, dt_lev); } @@ -229,6 +238,10 @@ REMORA::setup_step (int lev, Real time, Real dt_lev) Array4 const& pn = mf_pn->const_array(mfi); Array4 const& fcor = mf_fcor->const_array(mfi); + Array4 const& msku = mf_msku->const_array(mfi); + Array4 const& mskv = mf_mskv->const_array(mfi); + Array4 const& mskp = mf_mskp->const_array(mfi); + Box bx = mfi.tilebox(); Box tbxp1 = bx; @@ -280,7 +293,7 @@ REMORA::setup_step (int lev, Real time, Real dt_lev) FC(i,j,k)=0.0_rt; }); - prsgrd(tbxp1,gbx1,utbx,vtbx,ru,rv,pn,pm,rho,FC,Hz,z_r,z_w,nrhs,N); + prsgrd(tbxp1,gbx1,utbx,vtbx,ru,rv,pn,pm,rho,FC,Hz,z_r,z_w,msku,mskv,nrhs,N); // Apply mixing to temperature and, if use_salt, salt int ncomp = solverChoice.use_salt ? 2 : 1; @@ -288,10 +301,10 @@ REMORA::setup_step (int lev, Real time, Real dt_lev) Array4 const& s_arr_rhs = S_old.array(mfi); Array4 const& diff2_arr = vec_diff2[lev]->array(mfi); - t3dmix(bx, s_arr, s_arr_rhs, diff2_arr, Hz, pm, pn, dt_lev, ncomp); + t3dmix(bx, s_arr, s_arr_rhs, diff2_arr, Hz, pm, pn, msku, mskv, dt_lev, ncomp); Array4 const& diff2_arr_scalar = vec_diff2[lev]->array(mfi,Scalar_comp); - t3dmix(bx, S_new.array(mfi,Scalar_comp), S_old.array(mfi,Scalar_comp), diff2_arr_scalar, Hz, pm, pn, dt_lev, 1); + t3dmix(bx, S_new.array(mfi,Scalar_comp), S_old.array(mfi,Scalar_comp), diff2_arr_scalar, Hz, pm, pn, msku, mskv, dt_lev, 1); if (solverChoice.use_coriolis) { //----------------------------------------------------------------------- @@ -316,14 +329,15 @@ REMORA::setup_step (int lev, Real time, Real dt_lev) if(solverChoice.use_uv3dmix) { const int nnew = 0; - uv3dmix(xbx, ybx, u, v, uold, vold, rufrc, rvfrc, visc2_p, visc2_r, Hz, pm, pn, nrhs, nnew, dt_lev); + uv3dmix(xbx, ybx, u, v, uold, vold, rufrc, rvfrc, visc2_p, visc2_r, Hz, pm, pn, mskp, nrhs, nnew, dt_lev); } } // MFIter int nnew = (iic +1)% 2; nstp = iic % 2; if (solverChoice.vert_mixing_type == VertMixingType::GLS) { - gls_prestep(lev, mf_gls, mf_tke, mf_W, nstp, nnew, iic, ntfirst, N, dt_lev); + gls_prestep(lev, mf_gls, mf_tke, mf_W, mf_msku.get(), mf_mskv.get(), + nstp, nnew, iic, ntfirst, N, dt_lev); } nstp = 0; diff --git a/Source/TimeIntegration/REMORA_t3dmix.cpp b/Source/TimeIntegration/REMORA_t3dmix.cpp index 9126521..252626c 100644 --- a/Source/TimeIntegration/REMORA_t3dmix.cpp +++ b/Source/TimeIntegration/REMORA_t3dmix.cpp @@ -10,6 +10,8 @@ REMORA::t3dmix (const Box& bx, const Array4& Hz, const Array4& pm, const Array4& pn, + const Array4& msku, + const Array4& mskv, const Real dt_lev, const int ncomp) { //----------------------------------------------------------------------- @@ -31,6 +33,7 @@ REMORA::t3dmix (const Box& bx, const Real cff = 0.25_rt * (diff2(i,j,n) + diff2(i-1,j,n)) * pmon_u; FX(i,j,k,n) = cff * (Hz(i,j,k) + Hz(i-1,j,k)) * (state_rhs(i,j,k,n)-state_rhs(i-1,j,k,n)); + FX(i,j,k,n) *= msku(i,j,0); }); ParallelFor(ybx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) @@ -39,6 +42,7 @@ REMORA::t3dmix (const Box& bx, const Real cff = 0.25_rt*(diff2(i,j,n)+diff2(i,j-1,n)) * pnom_v; FE(i,j,k,n) = cff * (Hz(i,j,k) + Hz(i,j-1,k)) * (state_rhs(i,j,k,n) - state_rhs(i,j-1,k,n)); + FE(i,j,k,n) *= mskv(i,j,0); }); /* diff --git a/Source/TimeIntegration/REMORA_update_massflux_3d.cpp b/Source/TimeIntegration/REMORA_update_massflux_3d.cpp index ab9fe96..0b7dfe8 100644 --- a/Source/TimeIntegration/REMORA_update_massflux_3d.cpp +++ b/Source/TimeIntegration/REMORA_update_massflux_3d.cpp @@ -16,6 +16,7 @@ using namespace amrex; * @param[in ] Dphi_avg2 * @param[inout] DC * @param[inout] FC + * @param[in ] msk * @param[in ] nnew component of velocity */ @@ -31,6 +32,7 @@ REMORA::update_massflux_3d (const Box& bx, const Array4& Dphi_avg2, const Array4& DC, const Array4& FC, + const Array4& msk, const int nnew) { const Box& domain = geom[0].Domain(); @@ -85,14 +87,18 @@ REMORA::update_massflux_3d (const Box& bx, for (int k=0; k<=N; k++) { if (i == dlo.x-joff && (bcr_x.lo(0) == REMORABCType::ext_dir or ic_bc_type==IC_BC_Type::Real)) { phi(i,j,k,nnew) -= CF(i,j,0); + phi(i,j,k,nnew) *= msk(i,j,0); } else if (i == dhi.x+1 && (bcr_x.hi(0) == REMORABCType::ext_dir or ic_bc_type==IC_BC_Type::Real)) { phi(i,j,k,nnew) -= CF(i,j,0); + phi(i,j,k,nnew) *= msk(i,j,0); } if (j == dlo.y-ioff && (bcr_y.lo(1) == REMORABCType::ext_dir or ic_bc_type==IC_BC_Type::Real)) { phi(i,j,k,nnew) -= CF(i,j,0); + phi(i,j,k,nnew) *= msk(i,j,0); } else if (j == dhi.y+1 && (bcr_y.hi(1) == REMORABCType::ext_dir or ic_bc_type==IC_BC_Type::Real)) { phi(i,j,k,nnew) -= CF(i,j,0); + phi(i,j,k,nnew) *= msk(i,j,0); } } diff --git a/Source/TimeIntegration/REMORA_uv3dmix.cpp b/Source/TimeIntegration/REMORA_uv3dmix.cpp index bd2acf4..e2cf119 100644 --- a/Source/TimeIntegration/REMORA_uv3dmix.cpp +++ b/Source/TimeIntegration/REMORA_uv3dmix.cpp @@ -15,6 +15,7 @@ REMORA::uv3dmix (const Box& xbx, const Box& ybx, const Array4& Hz, const Array4& pm, const Array4& pn, + const Array4& mskp, int nrhs, int nnew, const Real dt_lev) { @@ -63,7 +64,8 @@ REMORA::uv3dmix (const Box& xbx, const Box& ybx, (pn(i-1,j-1,0)+pn(i-1,j,0)+pn(i,j-1,0)+pn(i,j,0)); const Real pnom_p = (pn(i-1,j-1,0)+pn(i-1,j,0)+pn(i,j-1,0)+pn(i,j,0)) / (pm(i-1,j-1,0)+pm(i-1,j,0)+pm(i,j-1,0)+pm(i,j,0)); - const Real cff = 0.125_rt * (Hz(i-1,j ,k) + Hz(i,j ,k)+ Hz(i-1,j-1,k) + Hz(i,j-1,k))* + const Real cff = mskp(i,j,0) * 0.125_rt * + (Hz(i-1,j ,k) + Hz(i,j ,k)+ Hz(i-1,j-1,k) + Hz(i,j-1,k))* (pmon_p* ((pn(i ,j-1,0)+pn(i ,j,0))*vold(i ,j,k,nrhs)- (pn(i-1,j-1,0)+pn(i-1,j,0))*vold(i-1,j,k,nrhs))+ @@ -118,7 +120,7 @@ REMORA::uv3dmix (const Box& xbx, const Box& ybx, (pn(i-1,j-1,0)+pn(i-1,j,0)+pn(i,j-1,0)+pn(i,j,0)); const Real pnom_p = (pn(i-1,j-1,0)+pn(i-1,j,0)+pn(i,j-1,0)+pn(i,j,0)) / (pm(i-1,j-1,0)+pm(i-1,j,0)+pm(i,j-1,0)+pm(i,j,0)); - const Real cff = 0.125_rt * (Hz(i-1,j ,k)+Hz(i,j ,k)+ + const Real cff = mskp(i,j,0) * 0.125_rt * (Hz(i-1,j ,k)+Hz(i,j ,k)+ Hz(i-1,j-1,k)+Hz(i,j-1,k))* (pmon_p* ((pn(i ,j-1,0)+pn(i ,j,0))*vold(i ,j,k,nrhs)- diff --git a/Source/TimeIntegration/REMORA_vert_mean_3d.cpp b/Source/TimeIntegration/REMORA_vert_mean_3d.cpp index 1ba9b34..b7431b4 100644 --- a/Source/TimeIntegration/REMORA_vert_mean_3d.cpp +++ b/Source/TimeIntegration/REMORA_vert_mean_3d.cpp @@ -14,6 +14,7 @@ REMORA::vert_mean_3d (const Box& phi_bx, const int ioff, const int joff, const Array4& DC, const Array4& CF, const Array4& dxlen, + const Array4& msk, const int nnew, const int N) { @@ -41,5 +42,6 @@ REMORA::vert_mean_3d (const Box& phi_bx, const int ioff, const int joff, ParallelFor(phi_bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { phi(i,j,k) -= DC(i,j,-1); + phi(i,j,k) *= msk(i,j,0); }); } From 6900b4c0efeb4fd321edebe428916f92fee42c7a Mon Sep 17 00:00:00 2001 From: Hannah Klion Date: Fri, 2 Aug 2024 17:18:35 -0700 Subject: [PATCH 2/2] IdealMiniGrid agrees with ROMS - Recalculate psi mask as in ROMS - Make sure masks are actually read from file - When read from file, one ghost cell of bathymetry is maintained, to match with ROMS - Pass masks around as needed --- .../BoundaryConditions_cons.cpp | 35 +++++++------- .../BoundaryConditions_netcdf.cpp | 19 ++++---- .../BoundaryConditions_xvel.cpp | 18 +++---- .../BoundaryConditions_yvel.cpp | 18 +++---- .../BoundaryConditions_zvel.cpp | 18 +++---- .../BoundaryConditions/REMORA_FillPatch.cpp | 14 ++++-- .../BoundaryConditions/REMORA_PhysBCFunct.H | 15 +++--- .../BoundaryConditions/REMORA_PhysBCFunct.cpp | 14 +++--- Source/IO/NCFile.H | 6 +++ Source/Initialization/REMORA_init.cpp | 22 --------- .../REMORA_init_from_netcdf.cpp | 7 ++- .../Initialization/REMORA_make_new_level.cpp | 47 +++++++++++++++++-- Source/REMORA.H | 8 +++- Source/REMORA.cpp | 8 ++++ Source/TimeIntegration/REMORA_advance_3d.cpp | 10 ++-- Source/TimeIntegration/REMORA_gls.cpp | 3 +- Source/TimeIntegration/REMORA_setup_step.cpp | 2 +- Source/TimeIntegration/REMORA_uv3dmix.cpp | 3 +- 18 files changed, 162 insertions(+), 105 deletions(-) diff --git a/Source/BoundaryConditions/BoundaryConditions_cons.cpp b/Source/BoundaryConditions/BoundaryConditions_cons.cpp index e4e547b..65c3c24 100644 --- a/Source/BoundaryConditions/BoundaryConditions_cons.cpp +++ b/Source/BoundaryConditions/BoundaryConditions_cons.cpp @@ -13,10 +13,13 @@ using namespace amrex; // time is the time at which the data should be filled // bccomp is the index into both domain_bcs_type_bcr and bc_extdir_vals for icomp = 0 -- // so this follows the BCVars enum +// n_not_fill perimeter of cells in x and y where BCs are not applied for conditions other than ext_dir. +// Foextrap is done based on the values at dom_lo-n_not_fill and dom_hi+n_not_fill. Reflecting conditions +// are untested. // void REMORAPhysBCFunct::impose_cons_bcs (const Array4& dest_arr, const Box& bx, const Box& domain, - const GpuArray /*dxInv*/, - int icomp, int ncomp, Real /*time*/, int bccomp) + const GpuArray /*dxInv*/, const Array4& mskr, + int icomp, int ncomp, Real /*time*/, int bccomp, int n_not_fill) { BL_PROFILE_VAR("impose_cons_bcs()",impose_cons_bcs); const auto& dom_lo = amrex::lbound(domain); @@ -64,12 +67,12 @@ void REMORAPhysBCFunct::impose_cons_bcs (const Array4& dest_arr, const Box ParallelFor( bx_xlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { if (bc_ptr[n].lo(0) == REMORABCType::ext_dir) { - dest_arr(i,j,k,icomp+n) = l_bc_extdir_vals_d[n][0]; + dest_arr(i,j,k,icomp+n) = l_bc_extdir_vals_d[n][0] * mskr(i,j,0); } }, bx_xhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { if (bc_ptr[n].hi(0) == REMORABCType::ext_dir) { - dest_arr(i,j,k,icomp+n) = l_bc_extdir_vals_d[n][3]; + dest_arr(i,j,k,icomp+n) = l_bc_extdir_vals_d[n][3] * mskr(i,j,0); } } ); @@ -82,12 +85,12 @@ void REMORAPhysBCFunct::impose_cons_bcs (const Array4& dest_arr, const Box ParallelFor( bx_ylo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { if (bc_ptr[n].lo(1) == REMORABCType::ext_dir) { - dest_arr(i,j,k,icomp+n) = l_bc_extdir_vals_d[n][1]; + dest_arr(i,j,k,icomp+n) = l_bc_extdir_vals_d[n][1] * mskr(i,j,0); } }, bx_yhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { if (bc_ptr[n].hi(1) == REMORABCType::ext_dir) { - dest_arr(i,j,k,icomp+n) = l_bc_extdir_vals_d[n][4]; + dest_arr(i,j,k,icomp+n) = l_bc_extdir_vals_d[n][4] * mskr(i,j,0); } } ); @@ -99,12 +102,12 @@ void REMORAPhysBCFunct::impose_cons_bcs (const Array4& dest_arr, const Box ParallelFor( bx_zlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { if (bc_ptr[n].lo(2) == REMORABCType::ext_dir) { - dest_arr(i,j,k,icomp+n) = l_bc_extdir_vals_d[n][2]; + dest_arr(i,j,k,icomp+n) = l_bc_extdir_vals_d[n][2] * mskr(i,j,0); } }, bx_zhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { if (bc_ptr[n].hi(2) == REMORABCType::ext_dir) { - dest_arr(i,j,k,icomp+n) = l_bc_extdir_vals_d[n][5]; + dest_arr(i,j,k,icomp+n) = l_bc_extdir_vals_d[n][5] * mskr(i,j,0); } } ); @@ -115,16 +118,16 @@ void REMORAPhysBCFunct::impose_cons_bcs (const Array4& dest_arr, const Box if (!is_periodic_in_x) { // Populate ghost cells on lo-x and hi-x domain boundaries - Box bx_xlo(bx); bx_xlo.setBig (0,dom_lo.x-1); + Box bx_xlo(bx); bx_xlo.setBig (0,dom_lo.x-1-n_not_fill); bx_xlo.setSmall(2,std::max(dom_lo.z,bx.smallEnd(2))); bx_xlo.setBig (2,std::min(dom_hi.z,bx.bigEnd(2))); - Box bx_xhi(bx); bx_xhi.setSmall(0,dom_hi.x+1); + Box bx_xhi(bx); bx_xhi.setSmall(0,dom_hi.x+1+n_not_fill); bx_xhi.setSmall(2,std::max(dom_lo.z,bx.smallEnd(2))); bx_xhi.setBig (2,std::min(dom_hi.z,bx.bigEnd(2))); ParallelFor(bx_xlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { int iflip = dom_lo.x - 1 - i; if (bc_ptr[n].lo(0) == REMORABCType::foextrap) { - dest_arr(i,j,k,icomp+n) = dest_arr(dom_lo.x,j,k,icomp+n); + dest_arr(i,j,k,icomp+n) = dest_arr(dom_lo.x-n_not_fill,j,k,icomp+n); } else if (bc_ptr[n].lo(0) == REMORABCType::reflect_even) { dest_arr(i,j,k,icomp+n) = dest_arr(iflip,j,k,icomp+n); } else if (bc_ptr[n].lo(0) == REMORABCType::reflect_odd) { @@ -134,7 +137,7 @@ void REMORAPhysBCFunct::impose_cons_bcs (const Array4& dest_arr, const Box bx_xhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { int iflip = 2*dom_hi.x + 1 - i; if (bc_ptr[n].hi(0) == REMORABCType::foextrap) { - dest_arr(i,j,k,icomp+n) = dest_arr(dom_hi.x,j,k,icomp+n); + dest_arr(i,j,k,icomp+n) = dest_arr(dom_hi.x+n_not_fill,j,k,icomp+n); } else if (bc_ptr[n].hi(0) == REMORABCType::reflect_even) { dest_arr(i,j,k,icomp+n) = dest_arr(iflip,j,k,icomp+n); } else if (bc_ptr[n].hi(0) == REMORABCType::reflect_odd) { @@ -147,17 +150,17 @@ void REMORAPhysBCFunct::impose_cons_bcs (const Array4& dest_arr, const Box if (!is_periodic_in_y) { // Populate ghost cells on lo-y and hi-y domain boundaries - Box bx_ylo(bx); bx_ylo.setBig (1,dom_lo.y-1); + Box bx_ylo(bx); bx_ylo.setBig (1,dom_lo.y-1-n_not_fill); bx_ylo.setSmall(2,std::max(dom_lo.z,bx.smallEnd(2))); bx_ylo.setBig (2,std::min(dom_hi.z,bx.bigEnd(2))); - Box bx_yhi(bx); bx_yhi.setSmall(1,dom_hi.y+1); + Box bx_yhi(bx); bx_yhi.setSmall(1,dom_hi.y+1+n_not_fill); bx_yhi.setSmall(2,std::max(dom_lo.z,bx.smallEnd(2))); bx_yhi.setBig (2,std::min(dom_hi.z,bx.bigEnd(2))); ParallelFor( bx_ylo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { int jflip = dom_lo.y - 1 - j; if (bc_ptr[n].lo(1) == REMORABCType::foextrap) { - dest_arr(i,j,k,icomp+n) = dest_arr(i,dom_lo.y,k,icomp+n); + dest_arr(i,j,k,icomp+n) = dest_arr(i,dom_lo.y-n_not_fill,k,icomp+n); } else if (bc_ptr[n].lo(1) == REMORABCType::reflect_even) { dest_arr(i,j,k,icomp+n) = dest_arr(i,jflip,k,icomp+n); } else if (bc_ptr[n].lo(1) == REMORABCType::reflect_odd) { @@ -167,7 +170,7 @@ void REMORAPhysBCFunct::impose_cons_bcs (const Array4& dest_arr, const Box bx_yhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { int jflip = 2*dom_hi.y + 1 - j; if (bc_ptr[n].hi(1) == REMORABCType::foextrap) { - dest_arr(i,j,k,icomp+n) = dest_arr(i,dom_hi.y,k,icomp+n); + dest_arr(i,j,k,icomp+n) = dest_arr(i,dom_hi.y+n_not_fill,k,icomp+n); } else if (bc_ptr[n].hi(1) == REMORABCType::reflect_even) { dest_arr(i,j,k,icomp+n) = dest_arr(i,jflip,k,icomp+n); } else if (bc_ptr[n].hi(1) == REMORABCType::reflect_odd) { diff --git a/Source/BoundaryConditions/BoundaryConditions_netcdf.cpp b/Source/BoundaryConditions/BoundaryConditions_netcdf.cpp index a9ecfbc..0430efd 100644 --- a/Source/BoundaryConditions/BoundaryConditions_netcdf.cpp +++ b/Source/BoundaryConditions/BoundaryConditions_netcdf.cpp @@ -11,7 +11,7 @@ using namespace amrex; */ void -REMORA::fill_from_bdyfiles (MultiFab& mf_to_fill, const Real time, const int bdy_var_type, const int icomp_to_fill) +REMORA::fill_from_bdyfiles (MultiFab& mf_to_fill, const MultiFab& mf_mask, const Real time, const int bdy_var_type, const int icomp_to_fill) { int lev = 0; @@ -99,36 +99,37 @@ REMORA::fill_from_bdyfiles (MultiFab& mf_to_fill, const Real time, const int bdy Box xhi_yhi = xhi & yhi; const Array4& dest_arr = mf_to_fill.array(mfi); + const Array4& mask_arr = mf_mask.array(mfi); if (!xlo.isEmpty()) { ParallelFor(xlo, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - dest_arr(i,j,k,icomp+icomp_to_fill) = oma * bdatxlo_n (ubound(xlo).x,j,k,0) - + alpha * bdatxlo_np1(ubound(xlo).x,j,k,0); + dest_arr(i,j,k,icomp+icomp_to_fill) = (oma * bdatxlo_n (ubound(xlo).x,j,k,0) + + alpha * bdatxlo_np1(ubound(xlo).x,j,k,0)) * mask_arr(i,j,0); }); } if (!xhi.isEmpty()) { ParallelFor(xhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - dest_arr(i,j,k,icomp+icomp_to_fill) = oma * bdatxhi_n (lbound(xhi).x,j,k,0) - + alpha * bdatxhi_np1(lbound(xhi).x,j,k,0); + dest_arr(i,j,k,icomp+icomp_to_fill) = (oma * bdatxhi_n (lbound(xhi).x,j,k,0) + + alpha * bdatxhi_np1(lbound(xhi).x,j,k,0)) * mask_arr(i,j,0); }); } if (!ylo.isEmpty()) { ParallelFor(ylo, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - dest_arr(i,j,k,icomp+icomp_to_fill) = oma * bdatylo_n (i,ubound(ylo).y,k,0) - + alpha * bdatylo_np1(i,ubound(ylo).y,k,0); + dest_arr(i,j,k,icomp+icomp_to_fill) = (oma * bdatylo_n (i,ubound(ylo).y,k,0) + + alpha * bdatylo_np1(i,ubound(ylo).y,k,0)) * mask_arr(i,j,0); }); } if (!yhi.isEmpty()) { ParallelFor(yhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - dest_arr(i,j,k,icomp+icomp_to_fill) = oma * bdatyhi_n (i,lbound(yhi).y,k,0) - + alpha * bdatyhi_np1(i,lbound(yhi).y,k,0); + dest_arr(i,j,k,icomp+icomp_to_fill) = (oma * bdatyhi_n (i,lbound(yhi).y,k,0) + + alpha * bdatyhi_np1(i,lbound(yhi).y,k,0)) * mask_arr(i,j,0); }); } diff --git a/Source/BoundaryConditions/BoundaryConditions_xvel.cpp b/Source/BoundaryConditions/BoundaryConditions_xvel.cpp index 50f5267..eae456a 100644 --- a/Source/BoundaryConditions/BoundaryConditions_xvel.cpp +++ b/Source/BoundaryConditions/BoundaryConditions_xvel.cpp @@ -10,7 +10,7 @@ using namespace amrex; // so this follows the BCVars enum // void REMORAPhysBCFunct::impose_xvel_bcs (const Array4& dest_arr, const Box& bx, const Box& domain, - const GpuArray /*dxInv*/, + const GpuArray /*dxInv*/, const Array4& msku, Real /*time*/, int bccomp) { BL_PROFILE_VAR("impose_xvel_bcs()",impose_xvel_bcs); @@ -60,7 +60,7 @@ void REMORAPhysBCFunct::impose_xvel_bcs (const Array4& dest_arr, const Box bx_xlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { int iflip = dom_lo.x - i; if (bc_ptr[n].lo(0) == REMORABCType::ext_dir) { - dest_arr(i,j,k) = l_bc_extdir_vals_d[n][0]; + dest_arr(i,j,k) = l_bc_extdir_vals_d[n][0]*msku(i,j,0); } else if (bc_ptr[n].lo(0) == REMORABCType::foextrap) { dest_arr(i,j,k) = dest_arr(dom_lo.x,j,k); } else if (bc_ptr[n].lo(0) == REMORABCType::reflect_even) { @@ -72,14 +72,14 @@ void REMORAPhysBCFunct::impose_xvel_bcs (const Array4& dest_arr, const Box // We only set the values on the domain faces themselves if EXT_DIR bx_xlo_face, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { if (bc_ptr[n].lo(0) == REMORABCType::ext_dir) - dest_arr(i,j,k) = l_bc_extdir_vals_d[n][0]; + dest_arr(i,j,k) = l_bc_extdir_vals_d[n][0]*msku(i,j,0); } ); ParallelFor( bx_xhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { int iflip = 2*(dom_hi.x + 1) - i; if (bc_ptr[n].hi(0) == REMORABCType::ext_dir) { - dest_arr(i,j,k) = l_bc_extdir_vals_d[n][3]; + dest_arr(i,j,k) = l_bc_extdir_vals_d[n][3]*msku(i,j,0); } else if (bc_ptr[n].hi(0) == REMORABCType::foextrap) { dest_arr(i,j,k) = dest_arr(dom_hi.x+1,j,k); } else if (bc_ptr[n].hi(0) == REMORABCType::reflect_even) { @@ -91,7 +91,7 @@ void REMORAPhysBCFunct::impose_xvel_bcs (const Array4& dest_arr, const Box // We only set the values on the domain faces themselves if EXT_DIR bx_xhi_face, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { if (bc_ptr[n].lo(3) == REMORABCType::ext_dir) - dest_arr(i,j,k) = l_bc_extdir_vals_d[n][3]; + dest_arr(i,j,k) = l_bc_extdir_vals_d[n][3]*msku(i,j,0); } ); } // not periodic in x @@ -105,7 +105,7 @@ void REMORAPhysBCFunct::impose_xvel_bcs (const Array4& dest_arr, const Box bx_ylo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { int jflip = dom_lo.y - 1 - j; if (bc_ptr[n].lo(1) == REMORABCType::ext_dir) { - dest_arr(i,j,k) = l_bc_extdir_vals_d[n][1]; + dest_arr(i,j,k) = l_bc_extdir_vals_d[n][1]*msku(i,j,0); } else if (bc_ptr[n].lo(1) == REMORABCType::foextrap) { dest_arr(i,j,k) = dest_arr(i,dom_lo.y,k); } else if (bc_ptr[n].lo(1) == REMORABCType::reflect_even) { @@ -117,7 +117,7 @@ void REMORAPhysBCFunct::impose_xvel_bcs (const Array4& dest_arr, const Box bx_yhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { int jflip = 2*dom_hi.y + 1 - j; if (bc_ptr[n].hi(1) == REMORABCType::ext_dir) { - dest_arr(i,j,k) = l_bc_extdir_vals_d[n][4]; + dest_arr(i,j,k) = l_bc_extdir_vals_d[n][4]*msku(i,j,0); } else if (bc_ptr[n].hi(1) == REMORABCType::foextrap) { dest_arr(i,j,k) = dest_arr(i,dom_hi.y,k); } else if (bc_ptr[n].hi(1) == REMORABCType::reflect_even) { @@ -137,7 +137,7 @@ void REMORAPhysBCFunct::impose_xvel_bcs (const Array4& dest_arr, const Box bx_zlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { int kflip = dom_lo.z - 1 - k; if (bc_ptr[n].lo(2) == REMORABCType::ext_dir) { - dest_arr(i,j,k) = l_bc_extdir_vals_d[n][2]; + dest_arr(i,j,k) = l_bc_extdir_vals_d[n][2]*msku(i,j,0); } else if (bc_ptr[n].lo(2) == REMORABCType::foextrap) { dest_arr(i,j,k) = dest_arr(i,j,dom_lo.z); } else if (bc_ptr[n].lo(2) == REMORABCType::reflect_even) { @@ -149,7 +149,7 @@ void REMORAPhysBCFunct::impose_xvel_bcs (const Array4& dest_arr, const Box bx_zhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { int kflip = 2*dom_hi.z + 1 - k; if (bc_ptr[n].hi(2) == REMORABCType::ext_dir) { - dest_arr(i,j,k) = l_bc_extdir_vals_d[n][5]; + dest_arr(i,j,k) = l_bc_extdir_vals_d[n][5]*msku(i,j,0); } else if (bc_ptr[n].hi(2) == REMORABCType::foextrap) { dest_arr(i,j,k) = dest_arr(i,j,dom_hi.z); } else if (bc_ptr[n].hi(2) == REMORABCType::reflect_even) { diff --git a/Source/BoundaryConditions/BoundaryConditions_yvel.cpp b/Source/BoundaryConditions/BoundaryConditions_yvel.cpp index 4505e6c..ec58555 100644 --- a/Source/BoundaryConditions/BoundaryConditions_yvel.cpp +++ b/Source/BoundaryConditions/BoundaryConditions_yvel.cpp @@ -10,7 +10,7 @@ using namespace amrex; // so this follows the BCVars enum // void REMORAPhysBCFunct::impose_yvel_bcs (const Array4& dest_arr, const Box& bx, const Box& domain, - const GpuArray /*dxInv*/, + const GpuArray /*dxInv*/, const Array4& mskv, Real /*time*/, int bccomp) { BL_PROFILE_VAR("impose_yvel_bcs()",impose_yvel_bcs); @@ -59,7 +59,7 @@ void REMORAPhysBCFunct::impose_yvel_bcs (const Array4& dest_arr, const Box bx_xlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { int iflip = dom_lo.x - 1- i; if (bc_ptr[n].lo(0) == REMORABCType::ext_dir) { - dest_arr(i,j,k) = l_bc_extdir_vals_d[n][0]; + dest_arr(i,j,k) = l_bc_extdir_vals_d[n][0]*mskv(i,j,0); } else if (bc_ptr[n].lo(0) == REMORABCType::foextrap) { dest_arr(i,j,k) = dest_arr(dom_lo.x,j,k); } else if (bc_ptr[n].lo(0) == REMORABCType::reflect_even) { @@ -71,7 +71,7 @@ void REMORAPhysBCFunct::impose_yvel_bcs (const Array4& dest_arr, const Box bx_xhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { int iflip = 2*dom_hi.x + 1 - i; if (bc_ptr[n].hi(0) == REMORABCType::ext_dir) { - dest_arr(i,j,k) = l_bc_extdir_vals_d[n][3]; + dest_arr(i,j,k) = l_bc_extdir_vals_d[n][3]*mskv(i,j,0); } else if (bc_ptr[n].hi(0) == REMORABCType::foextrap) { dest_arr(i,j,k) = dest_arr(dom_hi.x,j,k); } else if (bc_ptr[n].hi(0) == REMORABCType::reflect_even) { @@ -95,7 +95,7 @@ void REMORAPhysBCFunct::impose_yvel_bcs (const Array4& dest_arr, const Box bx_ylo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { int jflip = dom_lo.y-j; if (bc_ptr[n].lo(1) == REMORABCType::ext_dir) { - dest_arr(i,j,k) = l_bc_extdir_vals_d[n][1]; + dest_arr(i,j,k) = l_bc_extdir_vals_d[n][1]*mskv(i,j,0); } else if (bc_ptr[n].lo(1) == REMORABCType::foextrap) { dest_arr(i,j,k) = dest_arr(i,dom_lo.y,k); } else if (bc_ptr[n].lo(1) == REMORABCType::reflect_even) { @@ -107,14 +107,14 @@ void REMORAPhysBCFunct::impose_yvel_bcs (const Array4& dest_arr, const Box // We only set the values on the domain faces themselves if EXT_DIR bx_ylo_face, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { if (bc_ptr[n].lo(1) == REMORABCType::ext_dir) - dest_arr(i,j,k) = l_bc_extdir_vals_d[n][1]; + dest_arr(i,j,k) = l_bc_extdir_vals_d[n][1]*mskv(i,j,0); } ); ParallelFor( bx_yhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { int jflip = 2*(dom_hi.y + 1) - j; if (bc_ptr[n].hi(1) == REMORABCType::ext_dir) { - dest_arr(i,j,k) = l_bc_extdir_vals_d[n][4]; + dest_arr(i,j,k) = l_bc_extdir_vals_d[n][4]*mskv(i,j,0); } else if (bc_ptr[n].hi(1) == REMORABCType::foextrap) { dest_arr(i,j,k) = dest_arr(i,dom_hi.y+1,k); } else if (bc_ptr[n].hi(1) == REMORABCType::reflect_even) { @@ -126,7 +126,7 @@ void REMORAPhysBCFunct::impose_yvel_bcs (const Array4& dest_arr, const Box // We only set the values on the domain faces themselves if EXT_DIR bx_yhi_face, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { if (bc_ptr[n].lo(4) == REMORABCType::ext_dir) - dest_arr(i,j,k) = l_bc_extdir_vals_d[n][4]; + dest_arr(i,j,k) = l_bc_extdir_vals_d[n][4]*mskv(i,j,0); } ); } @@ -139,7 +139,7 @@ void REMORAPhysBCFunct::impose_yvel_bcs (const Array4& dest_arr, const Box bx_zlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { int kflip = dom_lo.z - 1 - k; if (bc_ptr[n].lo(2) == REMORABCType::ext_dir) { - dest_arr(i,j,k) = l_bc_extdir_vals_d[n][2]; + dest_arr(i,j,k) = l_bc_extdir_vals_d[n][2]*mskv(i,j,0); } else if (bc_ptr[n].lo(2) == REMORABCType::foextrap) { dest_arr(i,j,k) = dest_arr(i,j,dom_lo.z); } else if (bc_ptr[n].lo(2) == REMORABCType::reflect_even) { @@ -151,7 +151,7 @@ void REMORAPhysBCFunct::impose_yvel_bcs (const Array4& dest_arr, const Box bx_zhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { int kflip = 2*dom_hi.z + 1 - k; if (bc_ptr[n].hi(2) == REMORABCType::ext_dir) { - dest_arr(i,j,k) = l_bc_extdir_vals_d[n][5]; + dest_arr(i,j,k) = l_bc_extdir_vals_d[n][5]*mskv(i,j,0); } else if (bc_ptr[n].hi(2) == REMORABCType::foextrap) { dest_arr(i,j,k) = dest_arr(i,j,dom_hi.z); } else if (bc_ptr[n].hi(2) == REMORABCType::reflect_even) { diff --git a/Source/BoundaryConditions/BoundaryConditions_zvel.cpp b/Source/BoundaryConditions/BoundaryConditions_zvel.cpp index b293822..09aca4c 100644 --- a/Source/BoundaryConditions/BoundaryConditions_zvel.cpp +++ b/Source/BoundaryConditions/BoundaryConditions_zvel.cpp @@ -11,7 +11,7 @@ using namespace amrex; // so this follows the BCVars enum // void REMORAPhysBCFunct::impose_zvel_bcs (const Array4& dest_arr, const Box& bx, const Box& domain, - const GpuArray /*dxInv*/, + const GpuArray /*dxInv*/,const Array4& mskr, Real /*time*/, int bccomp) { const auto& dom_lo = amrex::lbound(domain); @@ -63,7 +63,7 @@ void REMORAPhysBCFunct::impose_zvel_bcs (const Array4& dest_arr, const Box bx_xlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { int iflip = dom_lo.x - 1 - i; if (bc_ptr[n].lo(0) == REMORABCType::ext_dir) { - dest_arr(i,j,k) = l_bc_extdir_vals_d[n][0]; + dest_arr(i,j,k) = l_bc_extdir_vals_d[n][0]*mskr(i,j,0); } else if (bc_ptr[n].lo(0) == REMORABCType::foextrap) { dest_arr(i,j,k) = dest_arr(dom_lo.x,j,k); } else if (bc_ptr[n].lo(0) == REMORABCType::reflect_even) { @@ -75,7 +75,7 @@ void REMORAPhysBCFunct::impose_zvel_bcs (const Array4& dest_arr, const Box bx_xhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { int iflip = 2*dom_hi.x + 1 - i; if (bc_ptr[n].hi(0) == REMORABCType::ext_dir) { - dest_arr(i,j,k) = l_bc_extdir_vals_d[n][3]; + dest_arr(i,j,k) = l_bc_extdir_vals_d[n][3]*mskr(i,j,0); } else if (bc_ptr[n].hi(0) == REMORABCType::foextrap) { dest_arr(i,j,k) = dest_arr(dom_hi.x,j,k); } else if (bc_ptr[n].hi(0) == REMORABCType::reflect_even) { @@ -95,7 +95,7 @@ void REMORAPhysBCFunct::impose_zvel_bcs (const Array4& dest_arr, const Box ParallelFor(bx_ylo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { int jflip = dom_lo.y - 1 - j; if (bc_ptr[n].lo(1) == REMORABCType::ext_dir) { - dest_arr(i,j,k) = l_bc_extdir_vals_d[n][1]; + dest_arr(i,j,k) = l_bc_extdir_vals_d[n][1]*mskr(i,j,0); } else if (bc_ptr[n].lo(1) == REMORABCType::foextrap) { dest_arr(i,j,k) = dest_arr(i,dom_lo.y,k); } else if (bc_ptr[n].lo(1) == REMORABCType::reflect_even) { @@ -107,7 +107,7 @@ void REMORAPhysBCFunct::impose_zvel_bcs (const Array4& dest_arr, const Box bx_yhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { int jflip = 2*dom_hi.y + 1 - j; if (bc_ptr[n].hi(1) == REMORABCType::ext_dir) { - dest_arr(i,j,k) = l_bc_extdir_vals_d[n][4]; + dest_arr(i,j,k) = l_bc_extdir_vals_d[n][4]*mskr(i,j,0); } else if (bc_ptr[n].hi(1) == REMORABCType::foextrap) { dest_arr(i,j,k) = dest_arr(i,dom_hi.y,k); } else if (bc_ptr[n].hi(1) == REMORABCType::reflect_even) { @@ -130,7 +130,7 @@ void REMORAPhysBCFunct::impose_zvel_bcs (const Array4& dest_arr, const Box int n = 0; int kflip = dom_lo.z - k; if (bc_ptr[n].lo(2) == REMORABCType::ext_dir) { - dest_arr(i,j,k) = l_bc_extdir_vals_d[n][2]; + dest_arr(i,j,k) = l_bc_extdir_vals_d[n][2]*mskr(i,j,0); } else if (bc_ptr[n].lo(2) == REMORABCType::foextrap) { dest_arr(i,j,k) = dest_arr(i,j,dom_lo.z); } else if (bc_ptr[n].lo(2) == REMORABCType::reflect_even) { @@ -142,7 +142,7 @@ void REMORAPhysBCFunct::impose_zvel_bcs (const Array4& dest_arr, const Box ParallelFor(bx_zlo_face, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { if (bc_ptr[n].lo(2) == REMORABCType::ext_dir) { - dest_arr(i,j,k) = l_bc_extdir_vals_d[n][2]; + dest_arr(i,j,k) = l_bc_extdir_vals_d[n][2]*mskr(i,j,0); } } ); @@ -151,7 +151,7 @@ void REMORAPhysBCFunct::impose_zvel_bcs (const Array4& dest_arr, const Box bx_zhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { int kflip = 2*(dom_hi.z + 1) - k; if (bc_ptr[n].hi(5) == REMORABCType::ext_dir) { - dest_arr(i,j,k) = l_bc_extdir_vals_d[n][5]; + dest_arr(i,j,k) = l_bc_extdir_vals_d[n][5]*mskr(i,j,0); } else if (bc_ptr[n].hi(5) == REMORABCType::foextrap) { dest_arr(i,j,k) = dest_arr(i,j,dom_hi.z+1); } else if (bc_ptr[n].hi(5) == REMORABCType::reflect_even) { @@ -163,7 +163,7 @@ void REMORAPhysBCFunct::impose_zvel_bcs (const Array4& dest_arr, const Box }, bx_zhi_face, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { if (bc_ptr[n].hi(2) == REMORABCType::ext_dir) { - dest_arr(i,j,k) = l_bc_extdir_vals_d[n][5]; + dest_arr(i,j,k) = l_bc_extdir_vals_d[n][5]*mskr(i,j,0); } } ); diff --git a/Source/BoundaryConditions/REMORA_FillPatch.cpp b/Source/BoundaryConditions/REMORA_FillPatch.cpp index bb85e14..36400bc 100644 --- a/Source/BoundaryConditions/REMORA_FillPatch.cpp +++ b/Source/BoundaryConditions/REMORA_FillPatch.cpp @@ -21,7 +21,8 @@ REMORA::FillPatch (int lev, Real time, MultiFab& mf_to_fill, Vector c #endif const int icomp, const bool fill_all, - const bool fill_set) + const bool fill_set, + const int n_not_fill) { BL_PROFILE_VAR("REMORA::FillPatch()",REMORA_FillPatch); int bccomp; @@ -30,6 +31,8 @@ REMORA::FillPatch (int lev, Real time, MultiFab& mf_to_fill, Vector c Box mf_box(mf_to_fill.boxArray()[0]); bool is_2d = mf_box.length(2) == 1; + MultiFab* mask; + // // *************************************************************************** // The first thing we do is interpolate the momenta on the "valid" faces of @@ -71,20 +74,24 @@ REMORA::FillPatch (int lev, Real time, MultiFab& mf_to_fill, Vector c { bccomp = 0; mapper = &cell_cons_interp; + mask = vec_mskr[lev].get(); } else if (mf_box.ixType() == IndexType(IntVect(1,0,0))) { bccomp = BCVars::xvel_bc; mapper = &face_linear_interp; + mask = vec_msku[lev].get(); } else if (mf_box.ixType() == IndexType(IntVect(0,1,0))) { bccomp = BCVars::yvel_bc; mapper = &face_linear_interp; + mask = vec_mskv[lev].get(); } else { bccomp = BCVars::zvel_bc; mapper = &face_linear_interp; + mask = vec_mskr[lev].get(); } if (lev == 0) @@ -115,14 +122,14 @@ REMORA::FillPatch (int lev, Real time, MultiFab& mf_to_fill, Vector c // *************************************************************************** // Enforce physical boundary conditions - (*physbcs[lev])(mf_to_fill,icomp,ncomp,mf_to_fill.nGrowVect(),time,bccomp); + (*physbcs[lev])(mf_to_fill,*mask,icomp,ncomp,mf_to_fill.nGrowVect(),time,bccomp,n_not_fill); #ifdef REMORA_USE_NETCDF // Fill the data which is stored in the boundary data read from netcdf files if ( (solverChoice.ic_bc_type == IC_BC_Type::Real) && (lev==0) && (bdy_var_type != BdyVars::null) ) { - fill_from_bdyfiles (mf_to_fill,time,bdy_var_type, icomp); + fill_from_bdyfiles (mf_to_fill,*mask,time,bdy_var_type, icomp); } #endif @@ -194,7 +201,6 @@ REMORA::FillPatchNoBC (int lev, Real time, MultiFab& mf_to_fill, Vector const& /*dest*/, + amrex::Array4 const& /*mask*/, const int /*dcomp*/, const int /*numcomp*/, amrex::GeometryData const& /*geom*/, const amrex::Real /*time*/, const amrex::BCRec* /*bcr*/, const int /*bcomp*/, @@ -57,26 +58,26 @@ public: // bccomp is the index into both domain_bcs_type_bcr and bc_extdir_vals for icomp = 0 -- // so this follows the BCVars enum // - void operator() (amrex::MultiFab& mf, int icomp, int ncomp, amrex::IntVect const& nghost, - amrex::Real time, int bccomp); + void operator() (amrex::MultiFab& mf, const amrex::MultiFab& mask, int icomp, int ncomp, amrex::IntVect const& nghost, + amrex::Real time, int bccomp, int n_not_fill=0); void impose_xvel_bcs (const amrex::Array4& dest_arr, const amrex::Box& bx, const amrex::Box& domain, - const amrex::GpuArray dxInv, + const amrex::GpuArray dxInv, const amrex::Array4& msku, amrex::Real time, int bccomp); void impose_yvel_bcs (const amrex::Array4& dest_arr, const amrex::Box& bx, const amrex::Box& domain, - const amrex::GpuArray dxInv, + const amrex::GpuArray dxInv, const amrex::Array4& mskv, amrex::Real time, int bccomp); void impose_zvel_bcs (const amrex::Array4& dest_arr, const amrex::Box& bx, const amrex::Box& domain, - const amrex::GpuArray dxInv, + const amrex::GpuArray dxInv, const amrex::Array4& mskr, amrex::Real time, int bccomp); void impose_cons_bcs (const amrex::Array4& mf, const amrex::Box& bx, const amrex::Box& domain, - const amrex::GpuArray dxInv, - int icomp, int ncomp, amrex::Real time, int bccomp); + const amrex::GpuArray dxInv, const amrex::Array4& mskr, + int icomp, int ncomp, amrex::Real time, int bccomp, int n_not_fill); // For backward compatibility // void FillBoundary (amrex::MultiFab& mf, int dcomp, int ncomp, amrex::IntVect const& nghost, diff --git a/Source/BoundaryConditions/REMORA_PhysBCFunct.cpp b/Source/BoundaryConditions/REMORA_PhysBCFunct.cpp index 83261a8..a37e0ec 100644 --- a/Source/BoundaryConditions/REMORA_PhysBCFunct.cpp +++ b/Source/BoundaryConditions/REMORA_PhysBCFunct.cpp @@ -16,8 +16,8 @@ using namespace amrex; // bccomp is the index into both domain_bcs_type_bcr and bc_extdir_vals for icomp = 0 -- // so this follows the BCVars enum // -void REMORAPhysBCFunct::operator() (MultiFab& mf, int icomp, int ncomp, IntVect const& nghost, - Real time, int bccomp) +void REMORAPhysBCFunct::operator() (MultiFab& mf, const MultiFab& msk, int icomp, int ncomp, IntVect const& nghost, + Real time, int bccomp,int n_not_fill) { if (m_geom.isAllPeriodic()) return; @@ -44,10 +44,11 @@ void REMORAPhysBCFunct::operator() (MultiFab& mf, int icomp, int ncomp, IntVect for (MFIter mfi(mf); mfi.isValid(); ++mfi) { const Array4& dest_arr = mf.array(mfi); + const Array4& msk_arr = msk.array(mfi); Box bx = mfi.validbox(); bx.grow(nghost); if (!gdomain.contains(bx)) { - impose_cons_bcs(dest_arr,bx,domain,dxInv,icomp,ncomp,time,bccomp); + impose_cons_bcs(dest_arr,bx,domain,dxInv,msk_arr,icomp,ncomp,time,bccomp,n_not_fill); } } // mfi @@ -57,16 +58,17 @@ void REMORAPhysBCFunct::operator() (MultiFab& mf, int icomp, int ncomp, IntVect for (MFIter mfi(mf); mfi.isValid(); ++mfi) { Box bx = mfi.validbox(); bx.grow(nghost); + const Array4& msk_arr = msk.array(mfi); if (!gdomain.contains(bx)) { if(bx.ixType() == IndexType(IntVect(1,0,0))) { const Array4& dest_arr = mf.array(mfi,icomp); - impose_xvel_bcs(dest_arr,bx,domain,dxInv,time,bccomp); + impose_xvel_bcs(dest_arr,bx,domain,dxInv,msk_arr,time,bccomp); } else if (bx.ixType() == IndexType(IntVect(0,1,0))) { const Array4& dest_arr = mf.array(mfi,icomp); - impose_yvel_bcs(dest_arr,bx,domain,dxInv,time,bccomp); + impose_yvel_bcs(dest_arr,bx,domain,dxInv,msk_arr,time,bccomp); } else if (bx.ixType() == IndexType(IntVect(0,0,1))) { const Array4& dest_arr = mf.array(mfi,icomp); - impose_zvel_bcs(dest_arr,bx,domain,dxInv,time,bccomp); + impose_zvel_bcs(dest_arr,bx,domain,dxInv,msk_arr,time,bccomp); } } } // mfi diff --git a/Source/IO/NCFile.H b/Source/IO/NCFile.H index b17d59f..2bc1e94 100644 --- a/Source/IO/NCFile.H +++ b/Source/IO/NCFile.H @@ -219,6 +219,12 @@ fill_fab_from_arrays (int iv, my_box.setBig(amrex::IntVect(nsx-2,nsy-1,nsz-1)); my_box.setType(amrex::IndexType(amrex::IntVect(0,1,0))); } + else if (var_name == "mask_psi") + { + my_box.setSmall(amrex::IntVect(0,0,0)); + my_box.setBig(amrex::IntVect(nsx-1,nsy-1,nsz-1)); + my_box.setType(amrex::IndexType(amrex::IntVect(1,1,0))); + } else // everything cell-centered -- these have one ghost cell { my_box.setSmall(amrex::IntVect(-1,-1,0)); diff --git a/Source/Initialization/REMORA_init.cpp b/Source/Initialization/REMORA_init.cpp index e13288d..07ab397 100644 --- a/Source/Initialization/REMORA_init.cpp +++ b/Source/Initialization/REMORA_init.cpp @@ -118,7 +118,6 @@ REMORA::set_2darrays (int lev) x_r(i,j,0) = prob_lo[0] + (i + 0.5_rt) * dx[0]; y_r(i,j,0) = prob_lo[1] + (j + 0.5_rt) * dx[1]; - // const Real z = prob_lo[2] + (k + 0.5_rt) * dx[2]; }); } @@ -130,9 +129,6 @@ REMORA::set_2darrays (int lev) MultiFab* V_old = yvel_new[lev]; std::unique_ptr& mf_ubar = vec_ubar[lev]; std::unique_ptr& mf_vbar = vec_vbar[lev]; - std::unique_ptr& mf_mskr = vec_mskr[lev]; - std::unique_ptr& mf_msku = vec_msku[lev]; - std::unique_ptr& mf_mskv = vec_mskv[lev]; std::unique_ptr& mf_Hz = vec_Hz[lev]; int nstp = 0; @@ -141,10 +137,6 @@ REMORA::set_2darrays (int lev) Array4 const& ubar = (mf_ubar)->array(mfi); Array4 const& vbar = (mf_vbar)->array(mfi); - Array4 const& mskr = (mf_mskr)->array(mfi); - Array4 const& msku = (mf_msku)->array(mfi); - Array4 const& mskv = (mf_mskv)->array(mfi); - Array4 const& Hz = mf_Hz->const_array(mfi); Array4 const& u = U_old->const_array(mfi); Array4 const& v = V_old->const_array(mfi); @@ -153,11 +145,6 @@ REMORA::set_2darrays (int lev) Box ubx2 = mfi.nodaltilebox(0); ubx2.grow(IntVect(NGROW ,NGROW ,0)); // x-face-centered, grown by 2 Box vbx2 = mfi.nodaltilebox(1); vbx2.grow(IntVect(NGROW ,NGROW ,0)); // y-face-centered, grown by 2 - ParallelFor(makeSlab(bx2,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int ) - { - mskr(i,j,0,0) = 1.0_rt; - }); - ParallelFor(makeSlab(ubx2,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int ) { Real CF = 0.; @@ -169,8 +156,6 @@ REMORA::set_2darrays (int lev) CF += avg_hz*u(i,j,k,nstp); } ubar(i,j,0,0) = CF / sum_of_hz; - - msku(i,j,0,0) = 1.0_rt; }); ParallelFor(makeSlab(vbx2,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int ) @@ -184,18 +169,11 @@ REMORA::set_2darrays (int lev) CF += avg_hz*v(i,j,k,nstp); } vbar(i,j,0,0) = CF / sum_of_hz; - - mskv(i,j,0,0) = 1.0_rt; }); } - 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); - - 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 diff --git a/Source/Initialization/REMORA_init_from_netcdf.cpp b/Source/Initialization/REMORA_init_from_netcdf.cpp index 5273646..630476e 100644 --- a/Source/Initialization/REMORA_init_from_netcdf.cpp +++ b/Source/Initialization/REMORA_init_from_netcdf.cpp @@ -152,7 +152,7 @@ REMORA::init_zeta_from_netcdf (int lev) } // omp } // idx vec_zeta[lev]->FillBoundary(geom[lev].periodicity()); - (*physbcs[lev])(*vec_zeta[lev],0,3,vec_zeta[lev]->nGrowVect(),t_old[lev],BCVars::cons_bc); + (*physbcs[lev])(*vec_zeta[lev],*vec_mskr[lev].get(),0,3,vec_zeta[lev]->nGrowVect(),t_old[lev],BCVars::cons_bc); } /** * REMORA function that initializes bathymetry from a netcdf file @@ -200,7 +200,8 @@ REMORA::init_bathymetry_from_netcdf (int lev) } // idx const double dummy_time = 0.0_rt; - FillPatch(lev,dummy_time,*vec_hOfTheConfusingName[lev],GetVecOfPtrs(vec_hOfTheConfusingName)); + FillPatch(lev,dummy_time,*vec_hOfTheConfusingName[lev],GetVecOfPtrs(vec_hOfTheConfusingName), + BdyVars::null,0,true,true,1); int ng = vec_pm[lev]->nGrow(); @@ -348,6 +349,8 @@ REMORA::init_masks_from_netcdf (int lev) } // mf } // omp } // idx + + update_mskp(lev); vec_mskr[lev]->FillBoundary(geom[lev].periodicity()); vec_msku[lev]->FillBoundary(geom[lev].periodicity()); vec_mskv[lev]->FillBoundary(geom[lev].periodicity()); diff --git a/Source/Initialization/REMORA_make_new_level.cpp b/Source/Initialization/REMORA_make_new_level.cpp index 8adb3ba..fa87b1a 100644 --- a/Source/Initialization/REMORA_make_new_level.cpp +++ b/Source/Initialization/REMORA_make_new_level.cpp @@ -445,9 +445,9 @@ void REMORA::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& vec_zeta[lev].reset(new MultiFab(ba2d,dm,3,IntVect(NGROW+1,NGROW+1,0))); // 2d free surface 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_msku[lev].reset(new MultiFab(convert(ba2d,IntVect(1,0,0)),dm,1,IntVect(NGROW+1,NGROW+1,0))); + vec_mskv[lev].reset(new MultiFab(convert(ba2d,IntVect(0,1,0)),dm,1,IntVect(NGROW+1,NGROW+1,0))); + vec_mskp[lev].reset(new MultiFab(convert(ba2d,IntVect(1,1,0)),dm,1,IntVect(NGROW+1,NGROW+1,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))); @@ -541,3 +541,44 @@ REMORA::set_zeta_to_Ztavg (int lev) } } + +void +REMORA::update_mskp (int lev) +{ + for ( MFIter mfi(*vec_mskr[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi ) + { + Array4 const& mskr = vec_mskr[lev]->const_array(mfi); + Array4< Real> const& mskp = vec_mskp[lev]->array(mfi); + + Box bx = mfi.tilebox(); bx.grow(IntVect(1,1,0)); bx.makeSlab(2,0); + + Real cff1 = 1.0_rt; + Real cff2 = 2.0_rt; + + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int) + { + if ((mskr(i-1,j,0) > 0.5) and (mskr(i,j,0) > 0.5) and (mskr(i-1,j-1,0) > 0.5) and (mskr(i,j-1,0) > 0.5)) { + mskp(i,j,0) = 1.0_rt; + } else if ((mskr(i-1,j,0) < 0.5) and (mskr(i,j,0) > 0.5) and (mskr(i-1,j-1,0) > 0.5) and (mskr(i,j-1,0) > 0.5)) { + mskp(i,j,0) = cff1; + } else if ((mskr(i-1,j,0) > 0.5) and (mskr(i,j,0) < 0.5) and (mskr(i-1,j-1,0) > 0.5) and (mskr(i,j-1,0) > 0.5)) { + mskp(i,j,0) = cff1; + } else if ((mskr(i-1,j,0) > 0.5) and (mskr(i,j,0) > 0.5) and (mskr(i-1,j-1,0) < 0.5) and (mskr(i,j-1,0) > 0.5)) { + mskp(i,j,0) = cff1; + } else if ((mskr(i-1,j,0) > 0.5) and (mskr(i,j,0) > 0.5) and (mskr(i-1,j-1,0) > 0.5) and (mskr(i,j-1,0) < 0.5)) { + mskp(i,j,0) = cff1; + } else if ((mskr(i-1,j,0) > 0.5) and (mskr(i,j,0) < 0.5) and (mskr(i-1,j-1,0) > 0.5) and (mskr(i,j-1,0) < 0.5)) { + mskp(i,j,0) = cff2; + } else if ((mskr(i-1,j,0) < 0.5) and (mskr(i,j,0) > 0.5) and (mskr(i-1,j-1,0) < 0.5) and (mskr(i,j-1,0) > 0.5)) { + mskp(i,j,0) = cff2; + } else if ((mskr(i-1,j,0) > 0.5) and (mskr(i,j,0) > 0.5) and (mskr(i-1,j-1,0) < 0.5) and (mskr(i,j-1,0) < 0.5)) { + mskp(i,j,0) = cff2; + } else if ((mskr(i-1,j,0) < 0.5) and (mskr(i,j,0) < 0.5) and (mskr(i-1,j-1,0) > 0.5) and (mskr(i,j-1,0) > 0.5)) { + mskp(i,j,0) = cff2; + } else { + mskp(i,j,0) = 0.0_rt; + } + + }); + } +} diff --git a/Source/REMORA.H b/Source/REMORA.H index d625eb3..d77b041 100644 --- a/Source/REMORA.H +++ b/Source/REMORA.H @@ -140,6 +140,9 @@ public: /** Set zeta components to be equal to time-averaged Zt_avg1 */ void set_zeta_to_Ztavg (int lev); + /** Set psi-point mask to be consistent with rho-point mask */ + void update_mskp (int lev); + /** compute dt from CFL considerations */ amrex::Real estTimeStep (int lev) const; @@ -640,6 +643,7 @@ public: 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_mskr, amrex::MultiFab* mf_msku, amrex::MultiFab* mf_mskv, const int nstp, const int nnew, const int N, const amrex::Real dt_lev); @@ -665,7 +669,8 @@ public: const int bdy_var_type = BdyVars::null, const int icomp=0, const bool fill_all=true, - const bool fill_set=true); + const bool fill_set=true, + const int n_not_fill=0); void FillPatchNoBC (int lev, amrex::Real time, amrex::MultiFab& mf_to_be_filled, @@ -676,6 +681,7 @@ public: const bool fill_set=true); void fill_from_bdyfiles (amrex::MultiFab& mf_to_fill, + const amrex::MultiFab& mf_mask, const amrex::Real time, const int bdy_var_type, const int icomp_to_fill); diff --git a/Source/REMORA.cpp b/Source/REMORA.cpp index 5209fc0..d69c150 100644 --- a/Source/REMORA.cpp +++ b/Source/REMORA.cpp @@ -623,6 +623,10 @@ REMORA::init_only (int lev, Real time) set_zeta_to_Ztavg(lev); amrex::Print() << "Initial data loaded from netcdf file \n " << std::endl; + amrex::Print() << "Calling init_masks_from_netcdf " << std::endl; + init_masks_from_netcdf(lev); + amrex::Print() << "Masks loaded from netcdf file \n " << std::endl; + amrex::Print() << "Calling init_bdry_from_netcdf " << std::endl; init_bdry_from_netcdf(); amrex::Print() << "Boundary data loaded from netcdf file \n " << std::endl; @@ -649,6 +653,10 @@ REMORA::init_only (int lev, Real time) set_zeta_to_Ztavg(lev); amrex::Print() << "Initial data loaded from netcdf file \n " << std::endl; + amrex::Print() << "Calling init_masks_from_netcdf " << std::endl; + init_masks_from_netcdf(lev); + amrex::Print() << "Masks loaded from netcdf file \n " << std::endl; + amrex::Print() << "Calling init_bdry_from_netcdf " << std::endl; init_bdry_from_netcdf(); amrex::Print() << "Boundary data loaded from netcdf file \n " << std::endl; diff --git a/Source/TimeIntegration/REMORA_advance_3d.cpp b/Source/TimeIntegration/REMORA_advance_3d.cpp index 3c76fe0..8eb5d5a 100644 --- a/Source/TimeIntegration/REMORA_advance_3d.cpp +++ b/Source/TimeIntegration/REMORA_advance_3d.cpp @@ -151,15 +151,15 @@ REMORA::advance_3d (int lev, MultiFab& mf_cons, } // Apply physical boundary conditions to u and v - (*physbcs[lev])(mf_u,0,1,mf_u.nGrowVect(),t_old[lev],BCVars::xvel_bc); - (*physbcs[lev])(mf_v,0,1,mf_v.nGrowVect(),t_old[lev],BCVars::yvel_bc); + (*physbcs[lev])(mf_u,*mf_msku,0,1,mf_u.nGrowVect(),t_old[lev],BCVars::xvel_bc); + (*physbcs[lev])(mf_v,*mf_mskv,0,1,mf_v.nGrowVect(),t_old[lev],BCVars::yvel_bc); #ifdef REMORA_USE_NETCDF // Fill the data which is stored in the boundary data read from netcdf files if ( (solverChoice.ic_bc_type == IC_BC_Type::Real) && (lev==0) ) { - fill_from_bdyfiles (mf_u,t_old[lev],BdyVars::u,0); - fill_from_bdyfiles (mf_v,t_old[lev],BdyVars::v,0); + fill_from_bdyfiles (mf_u,*mf_msku,t_old[lev],BdyVars::u,0); + fill_from_bdyfiles (mf_v,*mf_mskv,t_old[lev],BdyVars::v,0); } #endif @@ -303,7 +303,7 @@ REMORA::advance_3d (int lev, MultiFab& mf_cons, nnew = 1-nstp; if (solverChoice.vert_mixing_type == VertMixingType::GLS) { gls_corrector(lev, vec_gls[lev].get(), vec_tke[lev].get(), mf_W, vec_Akv[lev].get(), - vec_Akt[lev].get(),vec_Akk[lev].get(), vec_Akp[lev].get(), + vec_Akt[lev].get(),vec_Akk[lev].get(), vec_Akp[lev].get(), vec_mskr[lev].get(), vec_msku[lev].get(), vec_mskv[lev].get(), nstp, nnew, N, dt_lev); } diff --git a/Source/TimeIntegration/REMORA_gls.cpp b/Source/TimeIntegration/REMORA_gls.cpp index 215ba28..c06b9c2 100644 --- a/Source/TimeIntegration/REMORA_gls.cpp +++ b/Source/TimeIntegration/REMORA_gls.cpp @@ -215,6 +215,7 @@ void REMORA::gls_corrector (int lev, MultiFab* mf_gls, MultiFab* mf_tke, MultiFab& mf_W, MultiFab* mf_Akv, MultiFab* mf_Akt, MultiFab* mf_Akk, MultiFab* mf_Akp, + MultiFab* mf_mskr, MultiFab* mf_msku, MultiFab* mf_mskv, const int nstp, const int nnew, const int N, const Real dt_lev) @@ -416,7 +417,7 @@ REMORA::gls_corrector (int lev, MultiFab* mf_gls, MultiFab* mf_tke, }); } - (*physbcs[lev])(mf_shear2,0,1,mf_shear2.nGrowVect(),t_new[lev],BCVars::cons_bc); + (*physbcs[lev])(mf_shear2,*mf_mskr,0,1,mf_shear2.nGrowVect(),t_new[lev],BCVars::cons_bc); mf_CF.setVal(0.0_rt); for ( MFIter mfi(*mf_gls, TilingIfNotGPU()); mfi.isValid(); ++mfi ) diff --git a/Source/TimeIntegration/REMORA_setup_step.cpp b/Source/TimeIntegration/REMORA_setup_step.cpp index 696ff72..5e6e6cc 100644 --- a/Source/TimeIntegration/REMORA_setup_step.cpp +++ b/Source/TimeIntegration/REMORA_setup_step.cpp @@ -205,7 +205,7 @@ REMORA::setup_step (int lev, Real time, Real dt_lev) // We use FillBoundary not FillPatch here since mf_W is single-level scratch space mf_W.FillBoundary(geom[lev].periodicity()); - (*physbcs[lev])(mf_W,0,1,mf_W.nGrowVect(),t_new[lev],BCVars::zvel_bc); + (*physbcs[lev])(mf_W,*mf_mskr.get(),0,1,mf_W.nGrowVect(),t_new[lev],BCVars::zvel_bc); for ( MFIter mfi(S_old, TilingIfNotGPU()); mfi.isValid(); ++mfi ) { diff --git a/Source/TimeIntegration/REMORA_uv3dmix.cpp b/Source/TimeIntegration/REMORA_uv3dmix.cpp index e2cf119..2efca16 100644 --- a/Source/TimeIntegration/REMORA_uv3dmix.cpp +++ b/Source/TimeIntegration/REMORA_uv3dmix.cpp @@ -133,13 +133,14 @@ REMORA::uv3dmix (const Box& xbx, const Box& ybx, VFx(i,j,k) = on_p*on_p*visc2_p(i,j,0)*cff; }); - ParallelFor(makeSlab(ybx,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int k) + ParallelFor(makeSlab(ybx,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int ) { for (int k=0; k<=N; k++) { const Real cff=dt_lev*0.25_rt*(pm(i,j,0)+pm(i,j-1,0))*(pn(i,j,0)+pn(i,j-1,0)); const Real cff1=0.5_rt*(pn(i,j-1,0)+pn(i,j,0))*(VFx(i+1,j,k)-VFx(i,j ,k)); const Real cff2=0.5_rt*(pm(i,j-1,0)+pm(i,j,0))*(VFe(i ,j,k)-VFe(i,j-1,k)); const Real cff3=cff*(cff1-cff2); + v(i,j,k,nnew)=v(i,j,k,nnew)+cff3; rvfrc(i,j,0) += cff1-cff2; }