diff --git a/Exec/RegTests/DynamicRefinement/inputs_onelevel b/Exec/RegTests/DynamicRefinement/inputs_threelevel similarity index 70% rename from Exec/RegTests/DynamicRefinement/inputs_onelevel rename to Exec/RegTests/DynamicRefinement/inputs_threelevel index 2d5e66693..bd1681aa5 100644 --- a/Exec/RegTests/DynamicRefinement/inputs_onelevel +++ b/Exec/RegTests/DynamicRefinement/inputs_threelevel @@ -6,15 +6,13 @@ amrex.fpe_trap_invalid = 1 fabarray.mfiter_tile_size = 1024 1024 1024 # PROBLEM SIZE & GEOMETRY -geometry.prob_lo = -12 -12 -1 -geometry.prob_hi = 12 12 1 -amr.n_cell = 96 96 4 - geometry.prob_lo = -6 -6 -1 geometry.prob_hi = 6 6 1 amr.n_cell = 96 96 4 -amr.blocking_factor = 32 32 4 +amr.blocking_factor_x = 16 16 16 +amr.blocking_factor_y = 16 16 16 +amr.blocking_factor_z = 1 1 1 geometry.is_periodic = 1 1 0 @@ -31,16 +29,20 @@ erf.v = 1 # verbosity in ERF.cpp amr.v = 1 # verbosity in Amr.cpp # REFINEMENT / REGRIDDING -amr.max_level = 1 # maximum level number allowed -amr.ref_ratio_vect = 2 2 1 -erf.refinement_indicators = hi_scal1 +amr.max_level = 2 # maximum level number allowed +amr.ref_ratio_vect = 2 2 1 4 4 1 +erf.refinement_indicators = hi_scal1 hi_scal2 erf.hi_scal1.max_level = 1 erf.hi_scal1.field_name = scalar erf.hi_scal1.value_greater = 1. +erf.hi_scal2.max_level = 2 +erf.hi_scal2.field_name = scalar +erf.hi_scal2.value_greater = 2. amr.n_error_buf = 4 erf.regrid_int = 2 -erf.cf_width = 0 # Internal relaxation -erf.cf_set_width = 0 # Internal Dirichlet +erf.cf_width = 4 # Internal relaxation (not ready w/ multiple boxes at a level) +erf.cf_set_width = 1 # Internal Dirichlet (not ready w/ multiple boxes at a level) +amr.iterate_grids = true # Prevents level 1 BA decomposition (use true to decompose) # CHECKPOINT FILES erf.check_file = chk # root name of checkpoint file @@ -48,8 +50,8 @@ erf.check_int = -100 # number of timesteps between checkpoints # PLOTFILES erf.plot_file_1 = plt # prefix of plotfile name -erf.plot_int_1 = 100 # number of timesteps between plotfiles -erf.plot_vars_1 = density x_velocity y_velocity z_velocity pressure theta temp scalar pres_hse dens_hse pert_pres pert_dens +erf.plot_int_1 = 1 # number of timesteps between plotfiles +erf.plot_vars_1 = density rhotheta x_velocity y_velocity z_velocity pressure theta temp scalar pres_hse dens_hse pert_pres pert_dens mapfac # SOLVER CHOICE erf.alpha_T = 0.0 diff --git a/Exec/RegTests/DynamicRefinement/inputs_twolevel b/Exec/RegTests/DynamicRefinement/inputs_twolevel index f75a006f2..d7886f73d 100644 --- a/Exec/RegTests/DynamicRefinement/inputs_twolevel +++ b/Exec/RegTests/DynamicRefinement/inputs_twolevel @@ -6,15 +6,13 @@ amrex.fpe_trap_invalid = 1 fabarray.mfiter_tile_size = 1024 1024 1024 # PROBLEM SIZE & GEOMETRY -geometry.prob_lo = -12 -12 -1 -geometry.prob_hi = 12 12 1 -amr.n_cell = 96 96 4 - geometry.prob_lo = -6 -6 -1 geometry.prob_hi = 6 6 1 amr.n_cell = 96 96 4 -amr.blocking_factor = 32 32 4 +amr.blocking_factor_x = 32 32 +amr.blocking_factor_y = 32 32 +amr.blocking_factor_z = 1 1 geometry.is_periodic = 1 1 0 @@ -31,20 +29,16 @@ erf.v = 1 # verbosity in ERF.cpp amr.v = 1 # verbosity in Amr.cpp # REFINEMENT / REGRIDDING -amr.max_level = 2 # maximum level number allowed -amr.ref_ratio_vect = 2 2 1 4 4 1 -erf.refinement_indicators = hi_scal1 hi_scal2 +amr.max_level = 1 # maximum level number allowed +amr.ref_ratio_vect = 2 2 1 +erf.refinement_indicators = hi_scal1 erf.hi_scal1.max_level = 1 erf.hi_scal1.field_name = scalar erf.hi_scal1.value_greater = 1. -erf.hi_scal2.max_level = 2 -erf.hi_scal2.field_name = scalar -erf.hi_scal2.value_greater = 2. amr.n_error_buf = 4 erf.regrid_int = 2 -erf.cf_width = 0 # Internal relaxation (not ready w/ multiple boxes at a level) -erf.cf_set_width = 0 # Internal Dirichlet (not ready w/ multiple boxes at a level) -amr.iterate_grids = false # Prevents level 1 BA decomposition (use true to decompose) +erf.cf_width = 0 # Internal relaxation +erf.cf_set_width = 0 # Internal Dirichlet # CHECKPOINT FILES erf.check_file = chk # root name of checkpoint file @@ -52,7 +46,7 @@ erf.check_int = -100 # number of timesteps between checkpoints # PLOTFILES erf.plot_file_1 = plt # prefix of plotfile name -erf.plot_int_1 = 100 # number of timesteps between plotfiles +erf.plot_int_1 = 50 # number of timesteps between plotfiles erf.plot_vars_1 = density x_velocity y_velocity z_velocity pressure theta temp scalar pres_hse dens_hse pert_pres pert_dens # SOLVER CHOICE diff --git a/Source/BoundaryConditions/ERF_FillPatch.cpp b/Source/BoundaryConditions/ERF_FillPatch.cpp index 0b86c94eb..9890fe906 100644 --- a/Source/BoundaryConditions/ERF_FillPatch.cpp +++ b/Source/BoundaryConditions/ERF_FillPatch.cpp @@ -17,7 +17,7 @@ PhysBCFunctNoOp null_bc; */ void -ERF::FillPatch (int lev, Real time, const Vector& mfs) +ERF::FillPatch (int lev, Real time, const Vector& mfs, bool fillset) { BL_PROFILE_VAR("ERF::FillPatch()",ERF_FillPatch); int bccomp; @@ -73,11 +73,11 @@ ERF::FillPatch (int lev, Real time, const Vector& mfs) } // var_idx // Coarse-Fine set region - if (lev>0 && coupling_type=="OneWay" && cf_set_width>0) { - FPr_c[lev-1].fill(*mfs[Vars::cons], time, null_bc, domain_bcs_type, true); - FPr_u[lev-1].fill(*mfs[Vars::xvel], time, null_bc, domain_bcs_type, true); - FPr_v[lev-1].fill(*mfs[Vars::yvel], time, null_bc, domain_bcs_type, true); - FPr_w[lev-1].fill(*mfs[Vars::zvel], time, null_bc, domain_bcs_type, true); + if (lev>0 && coupling_type=="OneWay" && cf_set_width>0 && fillset) { + FPr_c[lev-1].FillSet(*mfs[Vars::cons], time, null_bc, domain_bcs_type); + FPr_u[lev-1].FillSet(*mfs[Vars::xvel], time, null_bc, domain_bcs_type); + FPr_v[lev-1].FillSet(*mfs[Vars::yvel], time, null_bc, domain_bcs_type); + FPr_w[lev-1].FillSet(*mfs[Vars::zvel], time, null_bc, domain_bcs_type); } // *************************************************************************** @@ -228,10 +228,14 @@ ERF::FillIntermediatePatch (int lev, Real time, // Coarse-Fine set region if (lev>0 && coupling_type=="OneWay" && cf_set_width>0) { - FPr_c[lev-1].fill(*mfs[Vars::cons], time, null_bc, domain_bcs_type, true); - FPr_u[lev-1].fill(*mfs[Vars::xvel], time, null_bc, domain_bcs_type, true); - FPr_v[lev-1].fill(*mfs[Vars::yvel], time, null_bc, domain_bcs_type, true); - FPr_w[lev-1].fill(*mfs[Vars::zvel], time, null_bc, domain_bcs_type, true); + if (cons_only) { + FPr_c[lev-1].FillSet(*mfs[Vars::cons], time, null_bc, domain_bcs_type); + } else { + FPr_c[lev-1].FillSet(*mfs[Vars::cons], time, null_bc, domain_bcs_type); + FPr_u[lev-1].FillSet(*mfs[Vars::xvel], time, null_bc, domain_bcs_type); + FPr_v[lev-1].FillSet(*mfs[Vars::yvel], time, null_bc, domain_bcs_type); + FPr_w[lev-1].FillSet(*mfs[Vars::zvel], time, null_bc, domain_bcs_type); + } } // *************************************************************************** @@ -251,10 +255,7 @@ ERF::FillIntermediatePatch (int lev, Real time, (*physbcs[lev])(mfs,icomp_cons,ncomp_cons,ngvect_cons,ngvect_vels,init_type,cons_only,BCVars::cons_bc,time); // *************************************************************************** - // - // It is important that we apply the MOST bcs after we have imposed all the others - // so that we have enough information in the ghost cells to calculate the viscosity - // + // MOST boundary conditions if (!(cons_only && ncomp_cons == 1) && m_most && allow_most_bcs) m_most->impose_most_bcs(lev,mfs,eddyDiffs); } diff --git a/Source/BoundaryConditions/ERF_FillPatcher.H b/Source/BoundaryConditions/ERF_FillPatcher.H index b52afe349..8a9f5b10b 100644 --- a/Source/BoundaryConditions/ERF_FillPatcher.H +++ b/Source/BoundaryConditions/ERF_FillPatcher.H @@ -2,37 +2,57 @@ #define ERF_FILLPATCHER_H_ #include +#include +#include class ERFFillPatcher { public: - ERFFillPatcher (amrex::BoxArray const& fba, amrex::DistributionMapping fdm, + ERFFillPatcher (amrex::BoxArray const& fba, amrex::DistributionMapping const& fdm, amrex::Geometry const& fgeom, - amrex::BoxArray cba, amrex::DistributionMapping cdm, + amrex::BoxArray const& cba, amrex::DistributionMapping const& cdm, amrex::Geometry const& cgeom, - int nghost, int nghost_subset, int ncomp, amrex::InterpBase* interp); - - ~ERFFillPatcher ( ) - { - delete m_cf_fine_data; - delete m_cf_fine_subset_data; - delete m_cf_crse_data[0]; - delete m_cf_crse_data[1]; - } + int nghost, int nghost_set, int ncomp, amrex::InterpBase* interp); - void Define (amrex::BoxArray const& fba, amrex::DistributionMapping fdm, + void Define (amrex::BoxArray const& fba, amrex::DistributionMapping const& fdm, amrex::Geometry const& fgeom, - amrex::BoxArray cba, amrex::DistributionMapping cdm, + amrex::BoxArray const& cba, amrex::DistributionMapping const& cdm, amrex::Geometry const& cgeom, - int nghost, int nghost_subset, int ncomp, amrex::InterpBase* interp); + int nghost, int nghost_sset, int ncomp, + amrex::InterpBase* interp); + + void BuildMask (amrex::BoxArray const& fba, int nghost, int nghost_set); - void registerCoarseData (amrex::Vector const& crse_data, + void RegisterCoarseData (amrex::Vector const& crse_data, amrex::Vector const& crse_time); + void InterpFace (amrex::MultiFab& fine, + amrex::MultiFab const& crse, + int mask_val); + + void InterpCell (amrex::MultiFab& fine, + amrex::MultiFab const& crse, + amrex::Vector const& bcr, + int mask_val); + + int GetSetMaskVal () { return m_set_mask; } + + int GetRelaxMaskVal () { return m_relax_mask; } + + amrex::iMultiFab* GetMask () { return m_cf_mask.get(); } + template - void fill (amrex::MultiFab& mf, amrex::Real time, - BC& cbc, amrex::Vector const& bcs, bool fill_subset=false); + void FillSet (amrex::MultiFab& mf, amrex::Real time, + BC& cbc, amrex::Vector const& bcs); + + template + void FillRelax (amrex::MultiFab& mf, amrex::Real time, + BC& cbc, amrex::Vector const& bcs); + + template + void Fill (amrex::MultiFab& mf, amrex::Real time, + BC& cbc, amrex::Vector const& bcs, int mask_val); private: @@ -47,15 +67,17 @@ private: int m_ncomp; amrex::InterpBase* m_interp; amrex::IntVect m_ratio; - amrex::Vector m_cf_crse_data; - amrex::MultiFab* m_cf_fine_data; - amrex::MultiFab* m_cf_fine_subset_data; + std::unique_ptr m_cf_crse_data_old; + std::unique_ptr m_cf_crse_data_new; + std::unique_ptr m_cf_mask; amrex::Vector m_crse_times; amrex::Real m_dt_crse; + int m_set_mask{2}; + int m_relax_mask{1}; }; /* - * Fill data at the coarse-fine boundary + * Fill fine data in the set region * * @param[out] mf MultiFab to be filled * @param[in] time Time at which to fill data @@ -64,8 +86,41 @@ private: */ template void -ERFFillPatcher::fill (amrex::MultiFab& mf, amrex::Real time, - BC& cbc, amrex::Vector const& bcs, bool fill_subset) +ERFFillPatcher::FillSet (amrex::MultiFab& mf, amrex::Real time, + BC& cbc, amrex::Vector const& bcs) +{ + Fill(mf,time,cbc,bcs,m_set_mask); +} + +/* + * Fill fine data in the relax region + * + * @param[out] mf MultiFab to be filled + * @param[in] time Time at which to fill data + * @param[in] cbc Coarse boundary condition + * @param[in] bcs Vector of boundary conditions + */ +template +void +ERFFillPatcher::FillRelax (amrex::MultiFab& mf, amrex::Real time, + BC& cbc, amrex::Vector const& bcs) +{ + Fill(mf,time,cbc,bcs,m_relax_mask); +} + +/* + * Fill fine data in the relax region + * + * @param[out] mf MultiFab to be filled + * @param[in] time Time at which to fill data + * @param[in] cbc Coarse boundary condition + * @param[in] bcs Vector of boundary conditions + * @param[in] mask_val Value to assign mask array + */ +template +void +ERFFillPatcher::Fill (amrex::MultiFab& mf, amrex::Real time, + BC& cbc, amrex::Vector const& bcs, int mask_val) { constexpr amrex::Real eps = std::numeric_limits::epsilon(); AMREX_ALWAYS_ASSERT((time >= m_crse_times[0]-eps) && (time <= m_crse_times[1]+eps)); @@ -75,32 +130,28 @@ ERFFillPatcher::fill (amrex::MultiFab& mf, amrex::Real time, amrex::Real fac_old = 1.0 - fac_new; // Boundary condition operator - cbc(*(m_cf_crse_data[0]), 0, m_ncomp, amrex::IntVect(0), time, 0); + cbc(*(m_cf_crse_data_old), 0, m_ncomp, amrex::IntVect(0), time, 0); // Coarse MF to hold time interpolated data - amrex::MultiFab crse_data_time_interp(m_cf_crse_data[0]->boxArray(), m_cf_crse_data[0]->DistributionMap(), + amrex::MultiFab crse_data_time_interp(m_cf_crse_data_old->boxArray(), + m_cf_crse_data_old->DistributionMap(), m_ncomp, amrex::IntVect{0}); // Time interpolate the coarse data amrex::MultiFab::LinComb(crse_data_time_interp, - fac_old, *(m_cf_crse_data[0]), 0, - fac_new, *(m_cf_crse_data[1]), 0, + fac_old, *(m_cf_crse_data_old), 0, + fac_new, *(m_cf_crse_data_new), 0, 0, m_ncomp, amrex::IntVect{0}); - // Ensure fine domain box is correct index type - amrex::Box fdest_dom = amrex::convert(m_fgeom.Domain(),m_cf_fine_data->boxArray().ixType()); - - // Spatially interpolate the time-interpolated coarse data - amrex::FillPatchInterp(*m_cf_fine_data, 0, crse_data_time_interp, 0, m_ncomp, amrex::IntVect(0), - m_cgeom, m_fgeom, fdest_dom, m_ratio, m_interp, bcs, 0); - - // Fill whole region or subset? - if (fill_subset) { - //amrex::MultiFab::Copy(m_cf_fine_subset_data, m_cf_fine_data, 0, 0, m_ncomp, amrex::IntVect{0}); - m_cf_fine_subset_data->ParallelCopy(*m_cf_fine_data, 0, 0, m_ncomp, amrex::IntVect{0}, amrex::IntVect{0}); - mf.ParallelCopy(*m_cf_fine_subset_data, 0, 0, m_ncomp, amrex::IntVect{0}, amrex::IntVect{0}); + // Call correct spatial interpolation type + amrex::IndexType m_ixt = mf.boxArray().ixType(); + int ixt_sum = m_ixt[0]+m_ixt[1]+m_ixt[2]; + if (ixt_sum == 0) { + InterpCell(mf,crse_data_time_interp,bcs,mask_val); + } else if (ixt_sum == 1) { + InterpFace(mf,crse_data_time_interp,mask_val); } else { - mf.ParallelCopy(*m_cf_fine_data, 0, 0, m_ncomp, amrex::IntVect{0}, amrex::IntVect{0}); + amrex::Abort("ERF_FillPatcher only supports face linear and cell cons linear interp!"); } } #endif diff --git a/Source/BoundaryConditions/ERF_FillPatcher.cpp b/Source/BoundaryConditions/ERF_FillPatcher.cpp index 381e68778..7721c8368 100644 --- a/Source/BoundaryConditions/ERF_FillPatcher.cpp +++ b/Source/BoundaryConditions/ERF_FillPatcher.cpp @@ -16,11 +16,11 @@ using namespace amrex; * @param[in] interp interpolation operator to be used */ -ERFFillPatcher::ERFFillPatcher (BoxArray const& fba, DistributionMapping fdm, +ERFFillPatcher::ERFFillPatcher (BoxArray const& fba, DistributionMapping const& fdm, Geometry const& fgeom, - BoxArray cba, DistributionMapping cdm, + BoxArray const& cba, DistributionMapping const& cdm, Geometry const& cgeom, - int nghost, int nghost_subset, + int nghost, int nghost_set, int ncomp, InterpBase* interp) : m_fba(fba), m_cba(cba), @@ -33,15 +33,10 @@ ERFFillPatcher::ERFFillPatcher (BoxArray const& fba, DistributionMapping fdm, // Vector to hold times for coarse data m_crse_times.resize(2); - m_cf_crse_data.resize(2); - - // Init MF patches - m_cf_fine_data = nullptr; m_cf_fine_subset_data = nullptr; - m_cf_crse_data[0] = nullptr; m_cf_crse_data[1] = nullptr; // Define the coarse and fine MFs Define(fba, fdm, fgeom, cba, cdm, cgeom, - nghost, nghost_subset, ncomp, interp); + nghost, nghost_set, ncomp, interp); } @@ -58,91 +53,107 @@ ERFFillPatcher::ERFFillPatcher (BoxArray const& fba, DistributionMapping fdm, * @param[in] ncomp number of components to be filled * @param[in] interp interpolation operator to be used */ -void ERFFillPatcher::Define (BoxArray const& fba, DistributionMapping fdm, +void ERFFillPatcher::Define (BoxArray const& fba, DistributionMapping const& fdm, Geometry const& fgeom, - BoxArray cba, DistributionMapping cdm, + BoxArray const& cba, DistributionMapping const& cdm, Geometry const& cgeom, - int nghost, int nghost_subset, + int nghost, int nghost_set, int ncomp, InterpBase* interp) { AMREX_ALWAYS_ASSERT(nghost < 0); - AMREX_ALWAYS_ASSERT(nghost_subset <= 0); - AMREX_ALWAYS_ASSERT(nghost <= nghost_subset); + AMREX_ALWAYS_ASSERT(nghost_set <= 0); + AMREX_ALWAYS_ASSERT(nghost <= nghost_set); // Set data memebers m_fba = fba; m_cba = cba; m_fdm = fdm; m_cdm = cdm; - m_fgeom = fgeom; m_cgeom = cgeom; - m_nghost = nghost; m_nghost_subset = nghost_subset; + m_fgeom = fgeom; m_cgeom = cgeom; + m_nghost = nghost; m_nghost_subset = nghost_set; m_ncomp = ncomp; m_interp = interp; // Delete old MFs if they exist - if (m_cf_fine_data) { - delete m_cf_fine_data; delete m_cf_fine_subset_data; - delete m_cf_crse_data[0]; delete m_cf_crse_data[1]; - } + if (m_cf_crse_data_old) m_cf_crse_data_old.reset(); + if (m_cf_crse_data_new) m_cf_crse_data_new.reset(); + if (m_cf_mask) m_cf_mask.reset(); // Index type for the BL/BA IndexType m_ixt = fba.ixType(); - // Box bounding all the fine boxes - Box const fine_valid_box = amrex::grow(fba.minimalBox(), IntVect(nghost,nghost,0)); - Box const fine_valid_box_subset = amrex::grow(fba.minimalBox(), IntVect(nghost_subset,nghost_subset,0)); - // Refinement ratios for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { m_ratio[idim] = m_fgeom.Domain().length(idim) / m_cgeom.Domain().length(idim); } - // Fine box list (interior halo regions exactly) - BoxList bl; - bl.set(m_ixt); - for (int ibox = 0; ibox < fba.size(); ++ibox) { - Box vbx = fba[ibox]; - BoxList const& bndry = amrex::boxDiff(vbx, fine_valid_box); - if (bndry.isNotEmpty()) { - bl.join(bndry); - } - } - - // Fine subset box list (interior halo regions exactly) - BoxList bl_subset; - bl_subset.set(m_ixt); - for (int ibox = 0; ibox < fba.size(); ++ibox) { - Box vbx = fba[ibox]; - BoxList const& bndry = amrex::boxDiff(vbx, fine_valid_box_subset); - if (bndry.isNotEmpty()) { - bl_subset.join(bndry); - } - } - - // Coarse box list (coincides with grown fine box list) + // Coarse box list BoxList cbl; cbl.set(m_ixt); - cbl.reserve(bl.size()); - for (auto const& b : bl) { - cbl.push_back(interp->CoarseBox(b, m_ratio)); + cbl.reserve(fba.size()); + for (int i(0); iCoarseBox(fba[i], m_ratio)); } - // Box arrays for fine and coarse - BoxArray cf_fba(std::move(bl)); + // Box arrays for the coarse data BoxArray cf_cba(std::move(cbl)); - BoxArray cf_fba_s(std::move(bl_subset)); - DistributionMapping cf_dm(cf_fba); - DistributionMapping cf_dm_s(cf_fba_s); + // Two coarse patches to hold the data to be interpolated + m_cf_crse_data_old = std::make_unique (cf_cba, fdm, m_ncomp, 0); + m_cf_crse_data_new = std::make_unique (cf_cba, fdm, m_ncomp, 0); - // Fine patch to hold the time-interpolated state - m_cf_fine_data = new MultiFab (cf_fba, cf_dm, m_ncomp, 0); - - // Fine subset patch to hold the time-interpolated state - m_cf_fine_subset_data = new MultiFab (cf_fba_s, cf_dm, m_ncomp, 0); + // Integer masking array + m_cf_mask = std::make_unique (fba, fdm, 1, 0); - // Two coarse patches to hold the data to be interpolated - m_cf_crse_data[0] = new MultiFab (cf_cba, cf_dm, m_ncomp, 0); - m_cf_crse_data[1] = new MultiFab (cf_cba, cf_dm, m_ncomp, 0); + // Populate mask array + if (nghost_set < 0) { + m_cf_mask->setVal(m_set_mask); + BuildMask(fba,nghost_set,m_set_mask-1); + } else { + m_cf_mask->setVal(m_relax_mask); + } + BuildMask(fba,nghost,m_relax_mask-1); } +void ERFFillPatcher::BuildMask (BoxArray const& fba, + int nghost, + int mask_val) +{ + // Minimal bounding box of fine BA plus a halo cell + Box fba_bnd = amrex::grow(fba.minimalBox(), IntVect(1,1,0)); + + // BoxList and BoxArray to store complement + BoxList com_bl; BoxArray com_ba; + + // Compute the complement + fba.complementIn(com_bl,fba_bnd); + + // com_bl cannot be null since we grew with halo cells + AMREX_ALWAYS_ASSERT(com_bl.size() > 0); + + // Grow the complement boxes and trim with the bounding box + Vector& com_bl_v = com_bl.data(); + for (int i(0); i& mask_arr = m_cf_mask->array(mfi); + + for (auto const& b : com_bl) { + Box com_bx = vbx & b; + amrex::ParallelFor(com_bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + mask_arr(i,j,k) = mask_val; + }); + } + } +} /* * Register the coarse data to be used by the ERFFillPatcher @@ -151,7 +162,7 @@ void ERFFillPatcher::Define (BoxArray const& fba, DistributionMapping fdm, * @param[in] crse_time times at which crse_data is defined */ -void ERFFillPatcher::registerCoarseData (Vector const& crse_data, +void ERFFillPatcher::RegisterCoarseData (Vector const& crse_data, Vector const& crse_time) { AMREX_ALWAYS_ASSERT(crse_data.size() == 2); // old and new @@ -162,14 +173,107 @@ void ERFFillPatcher::registerCoarseData (Vector const& crse_dat // m_cf_crse_data into ghost cells in the z-dir. So we need // to include ghost cells for crse_data when doing the copy IntVect src_ng = crse_data[0]->nGrowVect(); - IntVect dst_ng = m_cf_crse_data[0]->nGrowVect(); - m_cf_crse_data[0]->ParallelCopy(*crse_data[0], 0, 0, m_ncomp, - src_ng, dst_ng, m_cgeom.periodicity()); // old data - m_cf_crse_data[1]->ParallelCopy(*crse_data[1], 0, 0, m_ncomp, - src_ng, dst_ng, m_cgeom.periodicity()); // new data + IntVect dst_ng = m_cf_crse_data_old->nGrowVect(); + m_cf_crse_data_old->ParallelCopy(*(crse_data[0]), 0, 0, m_ncomp, + src_ng, dst_ng, m_cgeom.periodicity()); // old data + m_cf_crse_data_new->ParallelCopy(*(crse_data[1]), 0, 0, m_ncomp, + src_ng, dst_ng, m_cgeom.periodicity()); // new data m_crse_times[0] = crse_time[0]; // time of "old" coarse data m_crse_times[1] = crse_time[1]; // time of "new" coarse data m_dt_crse = crse_time[1] - crse_time[0]; } + + +void ERFFillPatcher::InterpFace (MultiFab& fine, + MultiFab const& crse, + int mask_val) +{ + int ncomp = m_ncomp; + IntVect ratio = m_ratio; + + for (MFIter mfi(fine); mfi.isValid(); ++mfi) { + Box const& fbx = mfi.validbox(); + + Array4 const& fine_arr = fine.array(mfi); + Array4 const& crse_arr = crse.const_array(mfi); + Array4 const& mask_arr = m_cf_mask->const_array(mfi); + + if (fbx.type(0) == IndexType::NODE) { + AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(RunOn::Gpu,fbx,ncomp,i,j,k,n, + { + if (mask_arr(i,j,k) == mask_val) face_linear_interp_x(i,j,k,n,fine_arr,crse_arr,ratio); + }); + } else if (fbx.type(1) == IndexType::NODE) { + AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(RunOn::Gpu,fbx,ncomp,i,j,k,n, + { + if (mask_arr(i,j,k) == mask_val) face_linear_interp_y(i,j,k,n,fine_arr,crse_arr,ratio); + }); + } else { + AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(RunOn::Gpu,fbx,ncomp,i,j,k,n, + { + if (mask_arr(i,j,k) == mask_val) face_linear_interp_z(i,j,k,n,fine_arr,crse_arr,ratio); + }); + } // IndexType::NODE + } // MFiter +} + +void ERFFillPatcher::InterpCell (MultiFab& fine, + MultiFab const& crse, + Vector const& bcr, + int mask_val) +{ + + int ncomp = m_ncomp; + IntVect ratio = m_ratio; + IndexType m_ixt = fine.boxArray().ixType(); + Box const& cdomain = amrex::convert(m_cgeom.Domain(), m_ixt); + + for (MFIter mfi(fine); mfi.isValid(); ++mfi) { + Box const& fbx = mfi.validbox(); + + Array4 const& fine_arr = fine.array(mfi); + Array4 const& crse_arr = crse.const_array(mfi); + Array4 const& mask_arr = m_cf_mask->const_array(mfi); + + bool run_on_gpu = Gpu::inLaunchRegion(); + amrex::ignore_unused(run_on_gpu); + + amrex::ignore_unused(m_fgeom); + + const Box& crse_region = m_interp->CoarseBox(fbx,ratio); + Box cslope_bx(crse_region); + for (int dim = 0; dim < AMREX_SPACEDIM; dim++) { + if (ratio[dim] > 1) { + cslope_bx.grow(dim,-1); + } + } + + FArrayBox ccfab(cslope_bx, ncomp*AMREX_SPACEDIM); + Array4 const& tmp = ccfab.array(); + Array4 const& ctmp = ccfab.const_array(); + +#ifdef AMREX_USE_GPU + AsyncArray async_bcr(bcr.data(), (run_on_gpu) ? ncomp : 0); + BCRec const* bcrp = (run_on_gpu) ? async_bcr.data() : bcr.data(); + + Elixir cceli; + if (run_on_gpu) { cceli = ccfab.elixir(); } +#else + BCRec const* bcrp = bcr.data(); +#endif + + AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(RunOn::Gpu, cslope_bx, ncomp, i, j, k, n, + { + mf_cell_cons_lin_interp_mcslope(i,j,k,n, tmp, crse_arr, 0, ncomp, + cdomain, ratio, bcrp); + }); + + AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(RunOn::Gpu, fbx, ncomp, i, j, k, n, + { + if (mask_arr(i,j,k) == mask_val) mf_cell_cons_lin_interp(i,j,k,n, fine_arr, 0, ctmp, + crse_arr, 0, ncomp, ratio); + }); + } // MFIter +} diff --git a/Source/ERF.H b/Source/ERF.H index 4c73186d1..7eb43bb7f 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -321,8 +321,6 @@ private: void update_terrain_arrays (int lev, amrex::Real time); - void define_grids_to_evolve (int lev, const amrex::BoxArray& ba); // NOLINT - void Construct_ERFFillPatchers (int lev); void Define_ERFFillPatchers (int lev); @@ -349,7 +347,7 @@ private: // Compute a vector of new MultiFabs by copying from valid region and filling ghost cells - // works for single level and 2-level cases (fill fine grid ghost by interpolating from coarse) - void FillPatch (int lev, amrex::Real time, const amrex::Vector& mf); + void FillPatch (int lev, amrex::Real time, const amrex::Vector& mf, bool fillset=true); #if defined(ERF_USE_MOISTURE) // Compute a new MultiFab by copying from valid region and filling ghost cells - @@ -513,9 +511,6 @@ private: amrex::Vector > > > mri_integrator_mem; amrex::Vector> physbcs; - // BoxArray at each level to define where we actually evolve the solution - amrex::Vector grids_to_evolve; - // Store Theta variable for MOST BC amrex::Vector> Theta_prim; diff --git a/Source/ERF.cpp b/Source/ERF.cpp index bcb358dd8..a5ebc4af1 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -133,8 +133,6 @@ ERF::ERF () nsubsteps[lev] = MaxRefRatio(lev-1); } - grids_to_evolve.resize(nlevs_max); - t_new.resize(nlevs_max, 0.0); t_old.resize(nlevs_max, -1.e100); dt.resize(nlevs_max, 1.e100); @@ -480,16 +478,6 @@ ERF::InitData () #endif } - // ************************************************************************* - // At this point, init_only has been called (from ERF::MakeNewLevelFromScratch), - // which sets lev_new and base_state - // ************************************************************************* - - // Define after wrfbdy_width is known - for (int lev = 0; lev <= finest_level; lev++) { - define_grids_to_evolve(lev, grids[lev]); - } - if (input_bndry_planes) { // Create the ReadBndryPlanes object so we can handle reading of boundary plane data amrex::Print() << "Defining r2d for the first time " << std::endl; @@ -521,7 +509,7 @@ ERF::InitData () // If not restarting we need to fill qmoist given qt and qp. if (restart_chkfile.empty()) { micro.Init(vars_new[lev][Vars::cons], qmoist[lev], - grids_to_evolve[lev], Geom(lev), 0.0); // dummy value, not needed just to diagnose + grids[lev], Geom(lev), 0.0); // dummy value, not needed just to diagnose micro.Update(vars_new[lev][Vars::cons], qmoist[lev]); } } @@ -1205,50 +1193,6 @@ ERF::AverageDownTo (int crse_lev) // NOLINT } } -void -ERF::define_grids_to_evolve (int lev, const BoxArray& ba) // NOLINT -{ - - Box domain(geom[lev].Domain()); - if (lev == 0 && ( init_type == "real" || init_type == "metgrid" ) ) - { - int width = wrfbdy_set_width; - Box shrunk_domain(domain); - shrunk_domain.grow(0,-width); - shrunk_domain.grow(1,-width); - BoxList bl; - int N = static_cast(ba.size()); - for (int i = 0; i < N; ++i) bl.push_back(ba[i] & shrunk_domain); - grids_to_evolve[lev].define(std::move(bl)); - } else if (lev == 1 && regrid_int < 0) { - int width = cf_set_width; - Box shrunk_domain(boxes_at_level[lev][0]); - shrunk_domain.grow(0,-width); - shrunk_domain.grow(1,-width); - BoxList bl; - int N = static_cast(ba.size()); - for (int i = 0; i < N; ++i) bl.push_back(ba[i] & shrunk_domain); - grids_to_evolve[lev].define(std::move(bl)); -#if 0 - if (num_boxes_at_level[lev] > 1) { - for (int i = 1; i < num_boxes_at_level[lev]; i++) { - int width = cf_set_width; - Box shrunk_domain(boxes_at_level[lev][i]); - shrunk_domain.grow(0,-width); - shrunk_domain.grow(1,-width); - BoxList bl; - int N = static_cast(ba.size()); - for (int i = 0; i < N; ++i) bl.push_back(ba[i] & shrunk_domain); - grids_to_evolve[lev].define(std::move(bl)); - } - } -#endif - } else { - // Just copy grids... - grids_to_evolve[lev] = ba; - } -} - void ERF::Construct_ERFFillPatchers (int lev) { @@ -1304,10 +1248,10 @@ ERF::Register_ERFFillPatchers (int lev) { auto& lev_new = vars_new[lev-1]; auto& lev_old = vars_old[lev-1]; - FPr_c[lev-1].registerCoarseData({&lev_old[Vars::cons], &lev_new[Vars::cons]}, {t_old[lev-1], t_new[lev-1]}); - FPr_u[lev-1].registerCoarseData({&lev_old[Vars::xvel], &lev_new[Vars::xvel]}, {t_old[lev-1], t_new[lev-1]}); - FPr_v[lev-1].registerCoarseData({&lev_old[Vars::yvel], &lev_new[Vars::yvel]}, {t_old[lev-1], t_new[lev-1]}); - FPr_w[lev-1].registerCoarseData({&lev_old[Vars::zvel], &lev_new[Vars::zvel]}, {t_old[lev-1], t_new[lev-1]}); + FPr_c[lev-1].RegisterCoarseData({&lev_old[Vars::cons], &lev_new[Vars::cons]}, {t_old[lev-1], t_new[lev-1]}); + FPr_u[lev-1].RegisterCoarseData({&lev_old[Vars::xvel], &lev_new[Vars::xvel]}, {t_old[lev-1], t_new[lev-1]}); + FPr_v[lev-1].RegisterCoarseData({&lev_old[Vars::yvel], &lev_new[Vars::yvel]}, {t_old[lev-1], t_new[lev-1]}); + FPr_w[lev-1].RegisterCoarseData({&lev_old[Vars::zvel], &lev_new[Vars::zvel]}, {t_old[lev-1], t_new[lev-1]}); } #ifdef ERF_USE_MULTIBLOCK @@ -1360,8 +1304,6 @@ ERF::ERF (const amrex::RealBox& rb, int max_level_in, nsubsteps[lev] = MaxRefRatio(lev-1); } - grids_to_evolve.resize(nlevs_max); - t_new.resize(nlevs_max, 0.0); t_old.resize(nlevs_max, -1.e100); dt.resize(nlevs_max, 1.e100); diff --git a/Source/ERF_make_new_level.cpp b/Source/ERF_make_new_level.cpp index a7ee11f43..e2b4ef88c 100644 --- a/Source/ERF_make_new_level.cpp +++ b/Source/ERF_make_new_level.cpp @@ -235,8 +235,6 @@ ERF::MakeNewLevelFromCoarse (int lev, Real time, const BoxArray& ba, // Build the data structures for calculating diffusive/turbulent terms update_arrays(lev, ba, dm); - define_grids_to_evolve(lev, grids[lev]); - FillCoarsePatch(lev, time, {&lev_new[Vars::cons],&lev_new[Vars::xvel], &lev_new[Vars::yvel],&lev_new[Vars::zvel]}); @@ -250,8 +248,6 @@ ERF::MakeNewLevelFromCoarse (int lev, Real time, const BoxArray& ba, void ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapping& dm) { - define_grids_to_evolve(lev, ba); - Vector temp_lev_new(Vars::NumTypes); Vector temp_lev_old(Vars::NumTypes); @@ -286,7 +282,7 @@ ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapp // This will fill the temporary MultiFabs with data from vars_new // ******************************************************************************************** FillPatch(lev, time, {&temp_lev_new[Vars::cons],&temp_lev_new[Vars::xvel], - &temp_lev_new[Vars::yvel],&temp_lev_new[Vars::zvel]}); + &temp_lev_new[Vars::yvel],&temp_lev_new[Vars::zvel]}, false); // ******************************************************************************************** // Copy from new into old just in case @@ -482,6 +478,4 @@ ERF::ClearLevel (int lev) // Clears the integrator memory mri_integrator_mem[lev].reset(); physbcs[lev].reset(); - - grids_to_evolve[lev].clear(); } diff --git a/Source/Microphysics/Init.cpp b/Source/Microphysics/Init.cpp index 15cf82404..9c2cac5de 100644 --- a/Source/Microphysics/Init.cpp +++ b/Source/Microphysics/Init.cpp @@ -15,17 +15,17 @@ using namespace amrex; * @param[in] qc_in Cloud variables input * @param[in,out] qv_in Vapor variables input * @param[in] qi_in Ice variables input - * @param[in] grids_to_evolve The boxes on which we will evolve the solution + * @param[in] grids The boxes on which we will evolve the solution * @param[in] geom Geometry associated with these MultiFabs and grids * @param[in] dt_advance Timestep for the advance */ void Microphysics::Init (const MultiFab& cons_in, MultiFab& qmoist, - const BoxArray& grids_to_evolve, + const BoxArray& grids, const Geometry& geom, const Real& dt_advance) { m_geom = geom; - m_gtoe = grids_to_evolve; + m_gtoe = grids; auto dz = m_geom.CellSize(2); auto lowz = m_geom.ProbLo(2); diff --git a/Source/Microphysics/Microphysics.H b/Source/Microphysics/Microphysics.H index 6a0473b98..d6a5cf660 100644 --- a/Source/Microphysics/Microphysics.H +++ b/Source/Microphysics/Microphysics.H @@ -90,7 +90,7 @@ class Microphysics { // init void Init (const amrex::MultiFab& cons_in, amrex::MultiFab& qmoist, - const amrex::BoxArray& grids_to_evolve, + const amrex::BoxArray& grids, const amrex::Geometry& geom, const amrex::Real& dt_advance); diff --git a/Source/Radiation/Radiation.H b/Source/Radiation/Radiation.H index 3dc812719..cfa7a7f6c 100644 --- a/Source/Radiation/Radiation.H +++ b/Source/Radiation/Radiation.H @@ -36,52 +36,52 @@ class Radiation { ~Radiation() = default; // init - void initialize(const amrex::MultiFab& cons_in, - const amrex::MultiFab& qc_in, - const amrex::MultiFab& qv_in, - const amrex::MultiFab& qi_in, - const amrex::BoxArray& grids_to_evolve, - const amrex::Geometry& geom, - const amrex::Real& dt_advance, - const bool& do_sw_rad, - const bool& do_lw_rad, - const bool& do_aero_rad, - const bool& do_snow_opt, - const bool& is_cmip6_volcano); + void initialize (const amrex::MultiFab& cons_in, + const amrex::MultiFab& qc_in, + const amrex::MultiFab& qv_in, + const amrex::MultiFab& qi_in, + const amrex::BoxArray& grids, + const amrex::Geometry& geom, + const amrex::Real& dt_advance, + const bool& do_sw_rad, + const bool& do_lw_rad, + const bool& do_aero_rad, + const bool& do_snow_opt, + const bool& is_cmip6_volcano); // run radiation model - void run(); + void run (); // call back - void on_complete(); - - void radiation_driver_lw(int ncol, int nlev, - real3d& gas_vmr, - real2d& pmid, real2d& pint, real2d& tmid, real2d& tint, - real3d& cld_tau_gpt, real3d& aer_tau_bnd, - FluxesByband& fluxes_clrsky, FluxesByband& fluxes_allsky, - real2d& qrl, real2d& qrlc); - - void radiation_driver_sw(int ncol, - real3d& gas_vmr, real2d& pmid, real2d& pint, real2d& tmid, - real2d& albedo_dir, real2d& albedo_dif, real1d& coszrs, - real3d& cld_tau_gpt, real3d& cld_ssa_gpt, real3d& cld_asm_gpt, - real3d& aer_tau_bnd, real3d& aer_ssa_bnd, real3d& aer_asm_bnd, - FluxesByband& fluxes_clrsky, FluxesByband& fluxes_allsky, - real2d& qrs, real2d& qrsc); - - void set_daynight_indices(real1d& coszrs, - int1d& day_indices, - int1d& night_indices); - - void get_gas_vmr(std::vector& gas_names, - real3d& gas_vmr); - - void calculate_heating_rate(real2d& flux_up, - real2d& flux_dn, - real2d& pint, - real2d& heating_rate); - private: + void on_complete (); + + void radiation_driver_lw (int ncol, int nlev, + real3d& gas_vmr, + real2d& pmid, real2d& pint, real2d& tmid, real2d& tint, + real3d& cld_tau_gpt, real3d& aer_tau_bnd, + FluxesByband& fluxes_clrsky, FluxesByband& fluxes_allsky, + real2d& qrl, real2d& qrlc); + + void radiation_driver_sw (int ncol, + real3d& gas_vmr, real2d& pmid, real2d& pint, real2d& tmid, + real2d& albedo_dir, real2d& albedo_dif, real1d& coszrs, + real3d& cld_tau_gpt, real3d& cld_ssa_gpt, real3d& cld_asm_gpt, + real3d& aer_tau_bnd, real3d& aer_ssa_bnd, real3d& aer_asm_bnd, + FluxesByband& fluxes_clrsky, FluxesByband& fluxes_allsky, + real2d& qrs, real2d& qrsc); + + void set_daynight_indices (real1d& coszrs, + int1d& day_indices, + int1d& night_indices); + + void get_gas_vmr (std::vector& gas_names, + real3d& gas_vmr); + + void calculate_heating_rate (real2d& flux_up, + real2d& flux_dn, + real2d& pint, + real2d& heating_rate); +private: // geometry amrex::Geometry m_geom; diff --git a/Source/Radiation/Radiation.cpp b/Source/Radiation/Radiation.cpp index b4b72e037..d60bdad54 100644 --- a/Source/Radiation/Radiation.cpp +++ b/Source/Radiation/Radiation.cpp @@ -76,7 +76,7 @@ void Radiation::initialize(const amrex::MultiFab& cons_in, const amrex::MultiFab& qc_in, const amrex::MultiFab& qv_in, const amrex::MultiFab& qi_in, - const amrex::BoxArray& grids_to_evolve, + const amrex::BoxArray& grids, const amrex::Geometry& geom, const amrex::Real& dt_advance, const bool& do_sw_rad, @@ -85,7 +85,7 @@ void Radiation::initialize(const amrex::MultiFab& cons_in, const bool& do_snow_opt, const bool& is_cmip6_volcano) { m_geom = geom; - m_gtoe = grids_to_evolve; + m_gtoe = grids; auto dz = m_geom.CellSize(2); auto lowz = m_geom.ProbLo(2); diff --git a/Source/TimeIntegration/ERF_TimeStep.cpp b/Source/TimeIntegration/ERF_TimeStep.cpp index 2c26ab0d2..5d5361c20 100644 --- a/Source/TimeIntegration/ERF_TimeStep.cpp +++ b/Source/TimeIntegration/ERF_TimeStep.cpp @@ -39,6 +39,10 @@ ERF::timeStep (int lev, Real time, int iteration) if (coupling_type=="OneWay" && cf_width>0) { Define_ERFFillPatchers(lev+1); Register_ERFFillPatchers(lev+1); + if (lev < max_level-1) { + Define_ERFFillPatchers(lev+2); + Register_ERFFillPatchers(lev+2); + } } // mark that we have regridded this level already diff --git a/Source/TimeIntegration/ERF_advance_dycore.cpp b/Source/TimeIntegration/ERF_advance_dycore.cpp index 5f2867adf..c700c67ca 100644 --- a/Source/TimeIntegration/ERF_advance_dycore.cpp +++ b/Source/TimeIntegration/ERF_advance_dycore.cpp @@ -313,10 +313,10 @@ void ERF::advance_dycore(int level, // Register coarse data for coarse-fine fill if (level0) { - FPr_c[level].registerCoarseData({&cons_old, &cons_new}, {old_time, old_time + dt_advance}); - FPr_u[level].registerCoarseData({&xvel_old, &xvel_new}, {old_time, old_time + dt_advance}); - FPr_v[level].registerCoarseData({&yvel_old, &yvel_new}, {old_time, old_time + dt_advance}); - FPr_w[level].registerCoarseData({&zvel_old, &zvel_new}, {old_time, old_time + dt_advance}); + FPr_c[level].RegisterCoarseData({&cons_old, &cons_new}, {old_time, old_time + dt_advance}); + FPr_u[level].RegisterCoarseData({&xvel_old, &xvel_new}, {old_time, old_time + dt_advance}); + FPr_v[level].RegisterCoarseData({&yvel_old, &yvel_new}, {old_time, old_time + dt_advance}); + FPr_w[level].RegisterCoarseData({&zvel_old, &zvel_new}, {old_time, old_time + dt_advance}); } if (verbose) Print() << "Done with advance_dycore at level " << level << std::endl; diff --git a/Source/TimeIntegration/ERF_advance_microphysics.cpp b/Source/TimeIntegration/ERF_advance_microphysics.cpp index 37de746b1..090a01935 100644 --- a/Source/TimeIntegration/ERF_advance_microphysics.cpp +++ b/Source/TimeIntegration/ERF_advance_microphysics.cpp @@ -3,12 +3,12 @@ using namespace amrex; #if defined(ERF_USE_MOISTURE) -void ERF::advance_microphysics(int lev, - MultiFab& cons, - const Real& dt_advance) +void ERF::advance_microphysics (int lev, + MultiFab& cons, + const Real& dt_advance) { micro.Init(cons, qmoist[lev], - grids_to_evolve[lev], + grids[lev], Geom(lev), dt_advance); diff --git a/Source/TimeIntegration/ERF_fast_rhs_MT.cpp b/Source/TimeIntegration/ERF_fast_rhs_MT.cpp index bf6b052d0..afb303590 100644 --- a/Source/TimeIntegration/ERF_fast_rhs_MT.cpp +++ b/Source/TimeIntegration/ERF_fast_rhs_MT.cpp @@ -17,7 +17,6 @@ using namespace amrex; * * @param[in] step which fast time step * @param[in] level level of resolution - * @param[in] grids_to_evolve the region in the domain excluding the relaxation and specified zones * @param[in] S_slow_rhs slow RHS computed in erf_slow_rhs_pre * @param[in] S_prev previous solution * @param[in] S_stg_data solution at previous RK stage @@ -47,7 +46,6 @@ using namespace amrex; */ void erf_fast_rhs_MT (int step, int /*level*/, - BoxArray& grids_to_evolve, Vector& S_slow_rhs, // the slow RHS already computed const Vector& S_prev, // if step == 0, this is S_old, else the previous solution Vector& S_stg_data, // at last RK stg: S^n, S^* or S^** @@ -117,10 +115,7 @@ void erf_fast_rhs_MT (int step, int /*level*/, // will require additional changes for ( MFIter mfi(S_stg_data[IntVar::cons],false); mfi.isValid(); ++mfi) { - // Construct intersection of current tilebox and valid region for updating - Box valid_bx = grids_to_evolve[mfi.index()]; - Box bx = mfi.tilebox() & valid_bx; - + Box bx = mfi.tilebox(); Box tbx = surroundingNodes(bx,0); Box tby = surroundingNodes(bx,1); Box tbz = surroundingNodes(bx,2); @@ -129,7 +124,7 @@ void erf_fast_rhs_MT (int step, int /*level*/, const Array4 & stg_xmom = S_stg_data[IntVar::xmom].const_array(mfi); const Array4 & stg_ymom = S_stg_data[IntVar::ymom].const_array(mfi); const Array4 & stg_zmom = S_stg_data[IntVar::zmom].const_array(mfi); - const Array4 & prim = S_stg_prim.const_array(mfi); + const Array4 & prim = S_stg_prim.const_array(mfi); const Array4& slow_rhs_cons = S_slow_rhs[IntVar::cons].const_array(mfi); const Array4& slow_rhs_rho_u = S_slow_rhs[IntVar::xmom].const_array(mfi); @@ -175,9 +170,9 @@ void erf_fast_rhs_MT (int step, int /*level*/, // Note: it is important to grow the tilebox rather than use growntilebox because // we need to fill the ghost cells of the tilebox so we can use them below - Box gbx = mfi.tilebox() & valid_bx; gbx.grow(1); - Box gtbx = mfi.nodaltilebox(0) & surroundingNodes(valid_bx,0); gtbx.grow(1); gtbx.setSmall(2,0); - Box gtby = mfi.nodaltilebox(1) & surroundingNodes(valid_bx,1); gtby.grow(1); gtby.setSmall(2,0); + Box gbx = mfi.tilebox(); gbx.grow(1); + Box gtbx = mfi.nodaltilebox(0); gtbx.grow(1); gtbx.setSmall(2,0); + Box gtby = mfi.nodaltilebox(1); gtby.grow(1); gtby.setSmall(2,0); { BL_PROFILE("fast_rhs_copies_0"); @@ -326,7 +321,7 @@ void erf_fast_rhs_MT (int step, int /*level*/, // This must be done before we set cur_xmom and cur_ymom, since those // in fact point to the same array as prev_xmom and prev_ymom // ********************************************************************* - Box gbxo = mfi.nodaltilebox(2) & surroundingNodes(valid_bx,2); + Box gbxo = mfi.nodaltilebox(2); { BL_PROFILE("fast_MT_making_omega"); Box gbxo_lo = gbxo; gbxo_lo.setBig(2,0); @@ -549,10 +544,6 @@ void erf_fast_rhs_MT (int step, int /*level*/, // ************************************************************************** // Define updates in the RHS of rho and (rho theta) // ************************************************************************** - - // We note that valid_bx is the actual grid, while bx may be a tile within that grid - // const auto& vbx_hi = amrex::ubound(valid_bx); - { BL_PROFILE("fast_rho_final_update"); amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept diff --git a/Source/TimeIntegration/ERF_fast_rhs_N.cpp b/Source/TimeIntegration/ERF_fast_rhs_N.cpp index c7d10a425..98cde76a9 100644 --- a/Source/TimeIntegration/ERF_fast_rhs_N.cpp +++ b/Source/TimeIntegration/ERF_fast_rhs_N.cpp @@ -16,7 +16,6 @@ using namespace amrex; * * @param[in] step which fast time step * @param[in] level level of resolution - * @param[in] grids_to_evolve the region in the domain excluding the relaxation and specified zones * @param[in] S_slow_rhs slow RHS computed in erf_slow_rhs_pre * @param[in] S_prev previous solution * @param[in] S_stage_data solution at previous RK stage @@ -36,7 +35,6 @@ using namespace amrex; */ void erf_fast_rhs_N (int step, int /*level*/, - BoxArray& grids_to_evolve, Vector& S_slow_rhs, // the slow RHS already computed const Vector& S_prev, // if step == 0, this is S_old, else the previous solution Vector& S_stage_data, // S_bar = S^n, S^* or S^** @@ -108,7 +106,7 @@ void erf_fast_rhs_N (int step, int /*level*/, const Array4 & cur_cons = S_data[IntVar::cons].array(mfi); const Array4& prev_cons = S_prev[IntVar::cons].const_array(mfi); const Array4& stage_cons = S_stage_data[IntVar::cons].const_array(mfi); - const Array4& lagged_delta_rt = S_scratch[IntVar::cons].array(mfi); + const Array4& lagged_delta_rt = S_scratch[IntVar::cons].array(mfi); const Array4& old_drho = Delta_rho.array(mfi); const Array4& old_drho_w = Delta_rho_w.array(mfi); @@ -117,18 +115,17 @@ void erf_fast_rhs_N (int step, int /*level*/, const Array4& prev_zmom = S_prev[IntVar::zmom].const_array(mfi); const Array4& stage_zmom = S_stage_data[IntVar::zmom].const_array(mfi); - Box valid_bx = grids_to_evolve[mfi.index()]; - Box gbx = mfi.tilebox() & valid_bx; gbx.grow(1); + Box gbx = mfi.tilebox(); gbx.grow(1); if (step == 0) { amrex::ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - cur_cons(i,j,k,Rho_comp) = prev_cons(i,j,k,Rho_comp); - cur_cons(i,j,k,RhoTheta_comp) = prev_cons(i,j,k,RhoTheta_comp); + cur_cons(i,j,k,Rho_comp) = prev_cons(i,j,k,Rho_comp); + cur_cons(i,j,k,RhoTheta_comp) = prev_cons(i,j,k,RhoTheta_comp); }); } // step = 0 - Box gtbz = mfi.nodaltilebox(2) & surroundingNodes(valid_bx,2); + Box gtbz = mfi.nodaltilebox(2); gtbz.grow(IntVect(1,1,0)); amrex::ParallelFor(gtbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { @@ -136,7 +133,6 @@ void erf_fast_rhs_N (int step, int /*level*/, }); const Array4& theta_extrap = extrap.array(mfi); - amrex::ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { old_drho(i,j,k) = cur_cons(i,j,k,Rho_comp) - stage_cons(i,j,k,Rho_comp); old_drho_theta(i,j,k) = cur_cons(i,j,k,RhoTheta_comp) - stage_cons(i,j,k,RhoTheta_comp); @@ -156,11 +152,10 @@ void erf_fast_rhs_N (int step, int /*level*/, for ( MFIter mfi(S_stage_data[IntVar::cons],TilingIfNotGPU()); mfi.isValid(); ++mfi) { // We define lagged_delta_rt for our next step as the current delta_rt - Box valid_bx = grids_to_evolve[mfi.index()]; - Box gbx = mfi.tilebox() & valid_bx; gbx.grow(1); + Box gbx = mfi.tilebox(); gbx.grow(1); - const Array4& lagged_delta_rt = S_scratch[IntVar::cons].array(mfi); - const Array4& old_drho_theta = Delta_rho_theta.array(mfi); + const Array4& lagged_delta_rt = S_scratch[IntVar::cons].array(mfi); + const Array4& old_drho_theta = Delta_rho_theta.array(mfi); amrex::ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { lagged_delta_rt(i,j,k,RhoTheta_comp) = old_drho_theta(i,j,k); @@ -176,11 +171,8 @@ void erf_fast_rhs_N (int step, int /*level*/, #endif for ( MFIter mfi(S_stage_data[IntVar::cons],TilingIfNotGPU()); mfi.isValid(); ++mfi) { - const Box& valid_bx = grids_to_evolve[mfi.index()]; - - // Construct intersection of current tilebox and valid region for updating - Box tbx = mfi.nodaltilebox(0) & surroundingNodes(valid_bx,0); - Box tby = mfi.nodaltilebox(1) & surroundingNodes(valid_bx,1); + Box tbx = mfi.nodaltilebox(0); + Box tby = mfi.nodaltilebox(1); const Array4 & stage_xmom = S_stage_data[IntVar::xmom].const_array(mfi); const Array4 & stage_ymom = S_stage_data[IntVar::ymom].const_array(mfi); @@ -275,9 +267,7 @@ void erf_fast_rhs_N (int step, int /*level*/, #endif for ( MFIter mfi(S_stage_data[IntVar::cons],TileNoZ()); mfi.isValid(); ++mfi) { - // Construct intersection of current tilebox and valid region for updating - Box bx = mfi.tilebox() & grids_to_evolve[mfi.index()]; - + Box bx = mfi.tilebox(); Box tbz = surroundingNodes(bx,2); const Array4 & stage_xmom = S_stage_data[IntVar::xmom].const_array(mfi); @@ -494,14 +484,14 @@ void erf_fast_rhs_N (int step, int /*level*/, Real zflux_lo = beta_2 * soln_a(i,j,k ) + beta_1 * old_drho_w(i,j,k ); Real zflux_hi = beta_2 * soln_a(i,j,k+1) + beta_1 * old_drho_w(i,j,k+1); - avg_zmom(i,j,k ) += facinv*zflux_lo; + avg_zmom(i,j,k) += facinv*zflux_lo; // Note that in the solve we effectively impose soln_a(i,j,vbx_hi.z+1)=0 // so we don't update avg_zmom at k=vbx_hi.z+1 - temp_rhs_arr(i,j,k,0) += ( zflux_hi - zflux_lo ) * dzi; - temp_rhs_arr(i,j,k,1) += 0.5 * dzi * - ( zflux_hi * (prim(i,j,k) + prim(i,j,k+1)) - zflux_lo * (prim(i,j,k) + prim(i,j,k-1)) ); + temp_rhs_arr(i,j,k,Rho_comp ) += dzi * ( zflux_hi - zflux_lo ); + temp_rhs_arr(i,j,k,RhoTheta_comp) += 0.5 * dzi * ( zflux_hi * (prim(i,j,k) + prim(i,j,k+1)) + - zflux_lo * (prim(i,j,k) + prim(i,j,k-1)) ); }); } // end profile } // mfi @@ -511,23 +501,23 @@ void erf_fast_rhs_N (int step, int /*level*/, #endif for ( MFIter mfi(S_stage_data[IntVar::cons],TilingIfNotGPU()); mfi.isValid(); ++mfi) { - const Box& bx = mfi.tilebox() & grids_to_evolve[mfi.index()]; + const Box& bx = mfi.tilebox(); - const Array4& cur_cons = S_data[IntVar::cons].array(mfi); - auto const& temp_rhs_arr = temp_rhs.const_array(mfi); - auto const& slow_rhs_cons = S_slow_rhs[IntVar::cons].const_array(mfi); + int cons_dycore{2}; + const Array4& cur_cons = S_data[IntVar::cons].array(mfi); + auto const& temp_rhs_arr = temp_rhs.const_array(mfi); + auto const& slow_rhs_cons = S_slow_rhs[IntVar::cons].const_array(mfi); - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + amrex::ParallelFor(bx, cons_dycore, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - cur_cons(i,j,k,0) += dtau * (slow_rhs_cons(i,j,k,0) - temp_rhs_arr(i,j,k,0)); - cur_cons(i,j,k,1) += dtau * (slow_rhs_cons(i,j,k,1) - temp_rhs_arr(i,j,k,1)); + cur_cons(i,j,k,n) += dtau * (slow_rhs_cons(i,j,k,n) - temp_rhs_arr(i,j,k,n)); }); - const Array4& cur_xmom = S_data[IntVar::xmom].array(mfi); - const Array4& cur_ymom = S_data[IntVar::ymom].array(mfi); + const Array4& cur_xmom = S_data[IntVar::xmom].array(mfi); + const Array4& cur_ymom = S_data[IntVar::ymom].array(mfi); - const Array4& temp_cur_xmom_arr = temp_cur_xmom.const_array(mfi); - const Array4& temp_cur_ymom_arr = temp_cur_ymom.const_array(mfi); + const Array4& temp_cur_xmom_arr = temp_cur_xmom.const_array(mfi); + const Array4& temp_cur_ymom_arr = temp_cur_ymom.const_array(mfi); Box tbx = surroundingNodes(bx,0); Box tby = surroundingNodes(bx,1); diff --git a/Source/TimeIntegration/ERF_fast_rhs_T.cpp b/Source/TimeIntegration/ERF_fast_rhs_T.cpp index 89812fefe..51a3134f2 100644 --- a/Source/TimeIntegration/ERF_fast_rhs_T.cpp +++ b/Source/TimeIntegration/ERF_fast_rhs_T.cpp @@ -17,7 +17,6 @@ using namespace amrex; * * @param[in] step which fast time step * @param[in] level level of resolution - * @param[in] grids_to_evolve the region in the domain excluding the relaxation and specified zones * @param[in] S_slow_rhs slow RHS computed in erf_slow_rhs_pre * @param[in] S_prev previous solution * @param[in] S_stage_data solution at previous RK stage @@ -40,7 +39,6 @@ using namespace amrex; */ void erf_fast_rhs_T (int step, int /*level*/, - BoxArray& grids_to_evolve, Vector& S_slow_rhs, // the slow RHS already computed const Vector& S_prev, // if step == 0, this is S_old, else the previous solution Vector& S_stage_data, // S_bar = S^n, S^* or S^** @@ -127,20 +125,19 @@ void erf_fast_rhs_T (int step, int /*level*/, const Array4& stage_ymom = S_stage_data[IntVar::ymom].const_array(mfi); const Array4& stage_zmom = S_stage_data[IntVar::zmom].const_array(mfi); - Box valid_bx = grids_to_evolve[mfi.index()]; - Box gbx = mfi.tilebox() & valid_bx; gbx.grow(1); + Box gbx = mfi.tilebox(); gbx.grow(1); if (step == 0) { amrex::ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - cur_cons(i,j,k,Rho_comp) = prev_cons(i,j,k,Rho_comp); - cur_cons(i,j,k,RhoTheta_comp) = prev_cons(i,j,k,RhoTheta_comp); + cur_cons(i,j,k,Rho_comp) = prev_cons(i,j,k,Rho_comp); + cur_cons(i,j,k,RhoTheta_comp) = prev_cons(i,j,k,RhoTheta_comp); }); } // step = 0 - Box gtbx = mfi.nodaltilebox(0) & surroundingNodes(valid_bx,0); gtbx.grow(IntVect(1,1,0)); - Box gtby = mfi.nodaltilebox(1) & surroundingNodes(valid_bx,1); gtby.grow(IntVect(1,1,0)); - Box gtbz = mfi.nodaltilebox(2) & surroundingNodes(valid_bx,2); gtbz.grow(IntVect(1,1,0)); + Box gtbx = mfi.nodaltilebox(0); gtbx.grow(IntVect(1,1,0)); + Box gtby = mfi.nodaltilebox(1); gtby.grow(IntVect(1,1,0)); + Box gtbz = mfi.nodaltilebox(2); gtbz.grow(IntVect(1,1,0)); amrex::ParallelFor(gtbx, gtby, gtbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { @@ -173,10 +170,9 @@ void erf_fast_rhs_T (int step, int /*level*/, for ( MFIter mfi(S_stage_data[IntVar::cons],TilingIfNotGPU()); mfi.isValid(); ++mfi) { // We define lagged_delta_rt for our next step as the current delta_rt - Box valid_bx = grids_to_evolve[mfi.index()]; - Box gbx = mfi.tilebox() & valid_bx; gbx.grow(1); - const Array4& old_drho_theta = Delta_rho_theta.array(mfi); - const Array4& lagged_delta_rt = S_scratch[IntVar::cons].array(mfi); + Box gbx = mfi.tilebox(); gbx.grow(1); + const Array4& old_drho_theta = Delta_rho_theta.array(mfi); + const Array4& lagged_delta_rt = S_scratch[IntVar::cons].array(mfi); amrex::ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { lagged_delta_rt(i,j,k,RhoTheta_comp) = old_drho_theta(i,j,k); }); @@ -191,10 +187,8 @@ void erf_fast_rhs_T (int step, int /*level*/, #endif for ( MFIter mfi(S_stage_data[IntVar::cons],TilingIfNotGPU()); mfi.isValid(); ++mfi) { - // Construct intersection of current tilebox and valid region for updating - Box valid_bx = grids_to_evolve[mfi.index()]; - Box tbx = mfi.nodaltilebox(0) & surroundingNodes(valid_bx,0); - Box tby = mfi.nodaltilebox(1) & surroundingNodes(valid_bx,1); + Box tbx = mfi.nodaltilebox(0); + Box tby = mfi.nodaltilebox(1); const Array4 & stage_xmom = S_stage_data[IntVar::xmom].const_array(mfi); const Array4 & stage_ymom = S_stage_data[IntVar::ymom].const_array(mfi); @@ -314,10 +308,7 @@ void erf_fast_rhs_T (int step, int /*level*/, #endif for ( MFIter mfi(S_stage_data[IntVar::cons],TileNoZ()); mfi.isValid(); ++mfi) { - // Construct intersection of current tilebox and valid region for updating - Box valid_bx = grids_to_evolve[mfi.index()]; - Box bx = mfi.tilebox() & valid_bx; - + Box bx = mfi.tilebox(); Box tbz = surroundingNodes(bx,2); const Array4 & stage_zmom = S_stage_data[IntVar::zmom].const_array(mfi); @@ -417,7 +408,7 @@ void erf_fast_rhs_T (int step, int /*level*/, } // end profile // ********************************************************************* - Box gbxo = mfi.nodaltilebox(2) & surroundingNodes(valid_bx,2); + Box gbxo = mfi.nodaltilebox(2); { BL_PROFILE("fast_T_making_omega"); Box gbxo_lo = gbxo; gbxo_lo.setBig(2,0); diff --git a/Source/TimeIntegration/ERF_make_buoyancy.cpp b/Source/TimeIntegration/ERF_make_buoyancy.cpp index c348095b7..7ad36127c 100644 --- a/Source/TimeIntegration/ERF_make_buoyancy.cpp +++ b/Source/TimeIntegration/ERF_make_buoyancy.cpp @@ -21,7 +21,6 @@ using namespace amrex; * equation for the z-component of momentum in the slow integrator. There * are three options for how buoyancy is computed (two are the same in the absence of moisture). * - * @param[in] grids_to_evolve the region in the domain excluding the relaxation and specified zones * @param[in] S_data current solution * @param[in] S_prim primitive variables (i.e. conserved variables divided by density) * @param[out] buoyancy the buoyancy term computed here @@ -36,8 +35,7 @@ using namespace amrex; * @param[in] r0 Reference (hydrostatically stratified) density */ -void make_buoyancy (BoxArray& grids_to_evolve, - Vector& S_data, +void make_buoyancy (Vector& S_data, const MultiFab& S_prim, MultiFab& buoyancy, #if defined(ERF_USE_MOISTURE) @@ -72,10 +70,7 @@ void make_buoyancy (BoxArray& grids_to_evolve, if (solverChoice.buoyancy_type == 1) { for ( MFIter mfi(buoyancy,TilingIfNotGPU()); mfi.isValid(); ++mfi) { - const Box& valid_bx = surroundingNodes(grids_to_evolve[mfi.index()],2); - - // Construct intersection of current tilebox and valid region for updating - Box tbz = mfi.tilebox() & valid_bx; + Box tbz = mfi.tilebox(); // We don't compute a source term for z-momentum on the bottom or top domain boundary if (tbz.smallEnd(2) == klo) tbz.growLo(2,-1); @@ -122,10 +117,7 @@ void make_buoyancy (BoxArray& grids_to_evolve, #endif for ( MFIter mfi(buoyancy,TilingIfNotGPU()); mfi.isValid(); ++mfi) { - const Box& valid_bx = surroundingNodes(grids_to_evolve[mfi.index()],2); - - // Construct intersection of current tilebox and valid region for updating - Box tbz = mfi.tilebox() & valid_bx; + Box tbz = mfi.tilebox(); // We don't compute a source term for z-momentum on the bottom or top boundary if (tbz.smallEnd(2) == klo) tbz.growLo(2,-1); @@ -162,10 +154,7 @@ void make_buoyancy (BoxArray& grids_to_evolve, for ( MFIter mfi(buoyancy,TilingIfNotGPU()); mfi.isValid(); ++mfi) { - const Box& valid_bx = surroundingNodes(grids_to_evolve[mfi.index()],2); - - // Construct intersection of current tilebox and valid region for updating - Box tbz = mfi.tilebox() & valid_bx; + Box tbz = mfi.tilebox(); // We don't compute a source term for z-momentum on the bottom or top domain boundary if (tbz.smallEnd(2) == klo) tbz.growLo(2,-1); @@ -234,10 +223,7 @@ void make_buoyancy (BoxArray& grids_to_evolve, #endif for ( MFIter mfi(buoyancy,TilingIfNotGPU()); mfi.isValid(); ++mfi) { - const Box& valid_bx = surroundingNodes(grids_to_evolve[mfi.index()],2); - - // Construct intersection of current tilebox and valid region for updating - Box tbz = mfi.tilebox() & valid_bx; + Box tbz = mfi.tilebox(); // We don't compute a source term for z-momentum on the bottom or top domain boundary if (tbz.smallEnd(2) == klo) tbz.growLo(2,-1); @@ -285,10 +271,7 @@ void make_buoyancy (BoxArray& grids_to_evolve, #endif for ( MFIter mfi(buoyancy,TilingIfNotGPU()); mfi.isValid(); ++mfi) { - const Box& valid_bx = surroundingNodes(grids_to_evolve[mfi.index()],2); - - // Construct intersection of current tilebox and valid region for updating - Box tbz = mfi.tilebox() & valid_bx; + Box tbz = mfi.tilebox(); // We don't compute a source term for z-momentum on the bottom or top domain boundary if (tbz.smallEnd(2) == klo) tbz.growLo(2,-1); @@ -332,10 +315,7 @@ void make_buoyancy (BoxArray& grids_to_evolve, for ( MFIter mfi(buoyancy,TilingIfNotGPU()); mfi.isValid(); ++mfi) { - const Box& valid_bx = surroundingNodes(grids_to_evolve[mfi.index()],2); - - // Construct intersection of current tilebox and valid region for updating - Box tbz = mfi.tilebox() & valid_bx; + Box tbz = mfi.tilebox(); // We don't compute a source term for z-momentum on the bottom or top domain boundary if (tbz.smallEnd(2) == klo) tbz.growLo(2,-1); diff --git a/Source/TimeIntegration/ERF_make_fast_coeffs.cpp b/Source/TimeIntegration/ERF_make_fast_coeffs.cpp index 68b36e011..c728636b6 100644 --- a/Source/TimeIntegration/ERF_make_fast_coeffs.cpp +++ b/Source/TimeIntegration/ERF_make_fast_coeffs.cpp @@ -12,7 +12,6 @@ using namespace amrex; * integrator (the acoustic substepping). * * @param[in] level level of refinement - * @param[in] grids_to_evolve the region in the domain excluding the relaxation and specified zones * @param[out] fast_coeffs the coefficients for the tridiagonal solver computed here * @param[in] S_stage_data solution at the last stage * @param[in] S_stage_prim primitive variables (i.e. conserved variables divided by density) at the last stage @@ -28,7 +27,6 @@ using namespace amrex; */ void make_fast_coeffs (int /*level*/, - BoxArray& grids_to_evolve, MultiFab& fast_coeffs, Vector& S_stage_data, // S_bar = S^n, S^* or S^** const MultiFab& S_stage_prim, @@ -73,11 +71,7 @@ void make_fast_coeffs (int /*level*/, for ( MFIter mfi(S_stage_data[IntVar::cons],TileNoZ()); mfi.isValid(); ++mfi) { - const Box& valid_bx = grids_to_evolve[mfi.index()]; - - // Construct intersection of current tilebox and valid region for updating - Box bx = mfi.tilebox() & valid_bx; - + Box bx = mfi.tilebox(); Box tbz = surroundingNodes(bx,2); const Array4 & stage_cons = S_stage_data[IntVar::cons].const_array(mfi); diff --git a/Source/TimeIntegration/ERF_slow_rhs_inc.cpp b/Source/TimeIntegration/ERF_slow_rhs_inc.cpp index 706be525e..86f658973 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_inc.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_inc.cpp @@ -23,7 +23,6 @@ using namespace amrex; * @param[in] level level of resolution * @param[in] nrk which RK stage * @param[in] dt slow time step - * @param[in] grids_to_evolve the region in the domain excluding the relaxation and specified zones * @param[out] S_rhs RHS computed here * @param[in] S_old old-time solution * @param[in] S_data current solution @@ -71,7 +70,6 @@ using namespace amrex; void erf_slow_rhs_inc (int /*level*/, int nrk, amrex::Real dt, - BoxArray& grids_to_evolve, Vector& S_rhs, Vector& S_old, Vector& S_data, @@ -170,9 +168,8 @@ void erf_slow_rhs_inc (int /*level*/, int nrk, #endif for ( MFIter mfi(S_data[IntVar::cons],TileNoZ()); mfi.isValid(); ++mfi) { - // Construct intersection of current tilebox and valid region for updating - const Box& valid_bx = grids_to_evolve[mfi.index()]; - Box bx = mfi.tilebox() & valid_bx; + const Box& bx = mfi.tilebox(); + const Box& valid_bx = mfi.validbox(); // Velocities const Array4 & u = xvel.array(mfi); @@ -211,11 +208,10 @@ void erf_slow_rhs_inc (int /*level*/, int nrk, //------------------------------------------------------------------------------- // Strain/Stress tile boxes - Box vbx = valid_bx; - Box bxcc = mfi.tilebox() & valid_bx; - Box tbxxy = mfi.tilebox(IntVect(1,1,0)) & vbx.convert(IntVect(1,1,0)); - Box tbxxz = mfi.tilebox(IntVect(1,0,1)) & vbx.convert(IntVect(1,0,1)); - Box tbxyz = mfi.tilebox(IntVect(0,1,1)) & vbx.convert(IntVect(0,1,1)); + Box bxcc = mfi.tilebox(); + Box tbxxy = mfi.tilebox(IntVect(1,1,0)); + Box tbxxz = mfi.tilebox(IntVect(1,0,1)); + Box tbxyz = mfi.tilebox(IntVect(0,1,1)); // We need a halo cells for terrain bxcc.grow(IntVect(1,1,0)); tbxxy.grow(IntVect(1,1,0)); @@ -470,14 +466,11 @@ void erf_slow_rhs_inc (int /*level*/, int nrk, #endif for ( MFIter mfi(S_data[IntVar::cons],TileNoZ()); mfi.isValid(); ++mfi) { - const Box& valid_bx = grids_to_evolve[mfi.index()]; - - // Construct intersection of current tilebox and valid region for updating - Box bx = mfi.tilebox() & valid_bx; - - Box tbx = mfi.nodaltilebox(0) & surroundingNodes(valid_bx,0); - Box tby = mfi.nodaltilebox(1) & surroundingNodes(valid_bx,1); - Box tbz = mfi.nodaltilebox(2) & surroundingNodes(valid_bx,2); + Box bx = mfi.tilebox(); + Box tbx = mfi.nodaltilebox(0); + Box tby = mfi.nodaltilebox(1); + Box tbz = mfi.nodaltilebox(2); + Box valid_bx = mfi.validbox(); // We don't compute a source term for z-momentum on the bottom or top boundary tbz.growLo(2,-1); diff --git a/Source/TimeIntegration/ERF_slow_rhs_post.cpp b/Source/TimeIntegration/ERF_slow_rhs_post.cpp index 98e3dc03a..62f7e2020 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_post.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_post.cpp @@ -397,7 +397,7 @@ void erf_slow_rhs_post (int level, #elif defined(ERF_USE_WARM_NO_PRECIP) icomp = RhoQv_comp; #endif - zero_RHS_in_set_region(icomp, 1, bx_xlo, bx_xhi, bx_ylo, bx_yhi, cell_rhs); + wrfbdy_zero_rhs_in_set_region(icomp, 1, bx_xlo, bx_xhi, bx_ylo, bx_yhi, cell_rhs); } #endif diff --git a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp index 61f7700cf..58a768be7 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp @@ -23,7 +23,6 @@ using namespace amrex; * @param[in] level level of resolution * @param[in] nrk which RK stage * @param[in] dt slow time step - * @param[in] grids_to_evolve the region in the domain excluding the relaxation and specified zones * @param[out] S_rhs RHS computed here * @param[in] S_data current solution * @param[in] S_prim primitive variables (i.e. conserved variables divided by density) @@ -71,7 +70,6 @@ using namespace amrex; void erf_slow_rhs_pre (int level, int nrk, amrex::Real dt, - BoxArray& grids_to_evolve, Vector& S_rhs, Vector& S_data, const MultiFab& S_prim, @@ -170,9 +168,8 @@ void erf_slow_rhs_pre (int level, #endif for ( MFIter mfi(S_data[IntVar::cons],TileNoZ()); mfi.isValid(); ++mfi) { - // Construct intersection of current tilebox and valid region for updating - const Box& valid_bx = grids_to_evolve[mfi.index()]; - Box bx = mfi.tilebox() & valid_bx; + const Box& bx = mfi.tilebox(); + const Box& valid_bx = mfi.validbox(); // Velocities const Array4 & u = xvel.array(mfi); @@ -211,12 +208,11 @@ void erf_slow_rhs_pre (int level, //------------------------------------------------------------------------------- // Strain/Stress tile boxes - Box vbx = valid_bx; - Box bxcc = mfi.tilebox() & valid_bx; - Box tbxxy = mfi.tilebox(IntVect(1,1,0)) & vbx.convert(IntVect(1,1,0)); - Box tbxxz = mfi.tilebox(IntVect(1,0,1)) & vbx.convert(IntVect(1,0,1)); - Box tbxyz = mfi.tilebox(IntVect(0,1,1)) & vbx.convert(IntVect(0,1,1)); - // We need a halo cells for terrain + Box bxcc = mfi.tilebox(); + Box tbxxy = mfi.tilebox(IntVect(1,1,0)); + Box tbxxz = mfi.tilebox(IntVect(1,0,1)); + Box tbxyz = mfi.tilebox(IntVect(0,1,1)); + // We need a halo cell for terrain bxcc.grow(IntVect(1,1,0)); tbxxy.grow(IntVect(1,1,0)); tbxxz.grow(IntVect(1,1,0)); @@ -468,14 +464,11 @@ void erf_slow_rhs_pre (int level, #endif for ( MFIter mfi(S_data[IntVar::cons],TileNoZ()); mfi.isValid(); ++mfi) { - const Box& valid_bx = grids_to_evolve[mfi.index()]; - - // Construct intersection of current tilebox and valid region for updating - Box bx = mfi.tilebox() & valid_bx; - - Box tbx = mfi.nodaltilebox(0) & surroundingNodes(valid_bx,0); - Box tby = mfi.nodaltilebox(1) & surroundingNodes(valid_bx,1); - Box tbz = mfi.nodaltilebox(2) & surroundingNodes(valid_bx,2); + Box bx = mfi.tilebox(); + Box tbx = mfi.nodaltilebox(0); + Box tby = mfi.nodaltilebox(1); + Box tbz = mfi.nodaltilebox(2); + Box valid_bx = mfi.validbox(); // We don't compute a source term for z-momentum on the bottom or top boundary tbz.growLo(2,-1); @@ -532,7 +525,7 @@ void erf_slow_rhs_pre (int level, //----------------------------------------- // Perturbational pressure field //----------------------------------------- - Box gbx = mfi.tilebox() & grids_to_evolve[mfi.index()]; gbx.grow(IntVect(1,1,0)); + Box gbx = mfi.tilebox(); gbx.grow(IntVect(1,1,0)); FArrayBox pprime; pprime.resize(gbx,1); Elixir pp_eli = pprime.elixir(); const Array4 & pp_arr = pprime.array(); diff --git a/Source/TimeIntegration/TI_fast_rhs_fun.H b/Source/TimeIntegration/TI_fast_rhs_fun.H index c0cd4bf3e..9d28b4dde 100644 --- a/Source/TimeIntegration/TI_fast_rhs_fun.H +++ b/Source/TimeIntegration/TI_fast_rhs_fun.H @@ -2,15 +2,15 @@ * Wrapper for calling the routine that creates the fast RHS */ auto fast_rhs_fun = [&](int fast_step, int n_sub, - Vector& S_slow_rhs, - const Vector& S_old, - Vector& S_stage, - Vector& S_data, - Vector& S_scratch, - const Real dtau, - const Real inv_fac, - const Real old_substep_time, - const Real new_substep_time) + Vector& S_slow_rhs, + const Vector& S_old, + Vector& S_stage, + Vector& S_data, + Vector& S_scratch, + const Real dtau, + const Real inv_fac, + const Real old_substep_time, + const Real new_substep_time) { BL_PROFILE("fast_rhs_fun"); if (verbose) amrex::Print() << "Calling fast rhs at level " << level << " with dt = " << dtau << std::endl; @@ -67,13 +67,13 @@ auto fast_rhs_fun = [&](int fast_step, int n_sub, // We have to call this each step since it depends on the substep time now // Note we pass in the *old* detJ here - make_fast_coeffs(level, grids_to_evolve[level], fast_coeffs, S_stage, S_prim, pi_stage, fine_geom, + make_fast_coeffs(level, fast_coeffs, S_stage, S_prim, pi_stage, fine_geom, solverChoice.use_terrain, solverChoice.gravity, solverChoice.c_p, detJ_cc[level], r0, pi0, dtau, beta_s); if (fast_step == 0) { // If this is the first substep we pass in S_old as the previous step's solution - erf_fast_rhs_MT(fast_step, level, grids_to_evolve[level], + erf_fast_rhs_MT(fast_step, level, S_slow_rhs, S_old, S_stage, S_prim, pi_stage, fast_coeffs, S_data, S_scratch, fine_geom, solverChoice.gravity, solverChoice.use_lagged_delta_rt, @@ -84,7 +84,7 @@ auto fast_rhs_fun = [&](int fast_step, int n_sub, mapfac_m[level], mapfac_u[level], mapfac_v[level]); } else { // If this is not the first substep we pass in S_data as the previous step's solution - erf_fast_rhs_MT(fast_step, level, grids_to_evolve[level], + erf_fast_rhs_MT(fast_step, level, S_slow_rhs, S_data, S_stage, S_prim, pi_stage, fast_coeffs, S_data, S_scratch, fine_geom, solverChoice.gravity, solverChoice.use_lagged_delta_rt, @@ -98,19 +98,19 @@ auto fast_rhs_fun = [&](int fast_step, int n_sub, if (fast_step == 0) { // If this is the first substep we make the coefficients since they are based only on stage data - make_fast_coeffs(level, grids_to_evolve[level], fast_coeffs, S_stage, S_prim, pi_stage, fine_geom, + make_fast_coeffs(level, fast_coeffs, S_stage, S_prim, pi_stage, fine_geom, solverChoice.use_terrain, solverChoice.gravity, solverChoice.c_p, detJ_cc[level], r0, pi0, dtau, beta_s); // If this is the first substep we pass in S_old as the previous step's solution - erf_fast_rhs_T(fast_step, level, grids_to_evolve[level], + erf_fast_rhs_T(fast_step, level, S_slow_rhs, S_old, S_stage, S_prim, pi_stage, fast_coeffs, S_data, S_scratch, fine_geom, solverChoice.gravity, Omega, z_phys_nd[level], detJ_cc[level], dtau, beta_s, inv_fac, mapfac_m[level], mapfac_u[level], mapfac_v[level]); } else { // If this is not the first substep we pass in S_data as the previous step's solution - erf_fast_rhs_T(fast_step, level, grids_to_evolve[level], + erf_fast_rhs_T(fast_step, level, S_slow_rhs, S_data, S_stage, S_prim, pi_stage, fast_coeffs, S_data, S_scratch, fine_geom, solverChoice.gravity, Omega, z_phys_nd[level], detJ_cc[level], dtau, beta_s, inv_fac, @@ -120,19 +120,19 @@ auto fast_rhs_fun = [&](int fast_step, int n_sub, if (fast_step == 0) { // If this is the first substep we make the coefficients since they are based only on stage data - make_fast_coeffs(level, grids_to_evolve[level], fast_coeffs, S_stage, S_prim, pi_stage, fine_geom, + make_fast_coeffs(level, fast_coeffs, S_stage, S_prim, pi_stage, fine_geom, solverChoice.use_terrain, solverChoice.gravity, solverChoice.c_p, detJ_cc[level], r0, pi0, dtau, beta_s); // If this is the first substep we pass in S_old as the previous step's solution - erf_fast_rhs_N(fast_step, level, grids_to_evolve[level], + erf_fast_rhs_N(fast_step, level, S_slow_rhs, S_old, S_stage, S_prim, pi_stage, fast_coeffs, S_data, S_scratch, fine_geom, solverChoice.gravity, dtau, beta_s, inv_fac, mapfac_m[level], mapfac_u[level], mapfac_v[level]); } else { // If this is not the first substep we pass in S_data as the previous step's solution - erf_fast_rhs_N(fast_step, level, grids_to_evolve[level], + erf_fast_rhs_N(fast_step, level, S_slow_rhs, S_data, S_stage, S_prim, pi_stage, fast_coeffs, S_data, S_scratch, fine_geom, solverChoice.gravity, dtau, beta_s, inv_fac, @@ -141,13 +141,15 @@ auto fast_rhs_fun = [&](int fast_step, int n_sub, } + // TODO: DELETE ALL OF THIS; FILL BCs PICKS THIS UP W/O G2E + /* #ifdef ERF_USE_NETCDF // Update vars in set zone (relaxation already updated) if (init_type=="real" && level==0 && wrfbdy_set_width>0) { if (fast_step == 0 ) { - update_interior_ghost(dtau, wrfbdy_set_width, boxes_at_level[level], S_slow_rhs, S_old, S_data); + update_interior_ghost(level, dtau, wrfbdy_set_width, boxes_at_level[level], S_slow_rhs, S_old, S_data); } else { - update_interior_ghost(dtau, wrfbdy_set_width, boxes_at_level[level], S_slow_rhs, S_data, S_data); + update_interior_ghost(level, dtau, wrfbdy_set_width, boxes_at_level[level], S_slow_rhs, S_data, S_data); } } #endif @@ -155,11 +157,12 @@ auto fast_rhs_fun = [&](int fast_step, int n_sub, // Update vars in set zone (relaxation already updated) if (level>0 && cf_set_width>0) { if (fast_step == 0 ) { - update_interior_ghost(dtau, cf_set_width, boxes_at_level[level], S_slow_rhs, S_old, S_data); + update_interior_ghost(level, dtau, cf_set_width, boxes_at_level[level], S_slow_rhs, S_old, S_data, &FPr_c[level-1], &FPr_u[level-1], &FPr_v[level-1], &FPr_w[level-1]); } else { - update_interior_ghost(dtau, cf_set_width, boxes_at_level[level], S_slow_rhs, S_data, S_data); + update_interior_ghost(level, dtau, cf_set_width, boxes_at_level[level], S_slow_rhs, S_data, S_data, &FPr_c[level-1], &FPr_u[level-1], &FPr_v[level-1], &FPr_w[level-1]); } } + */ // Even if we update all the conserved variables we don't need // to fillpatch the slow ones every acoustic substep diff --git a/Source/TimeIntegration/TI_headers.H b/Source/TimeIntegration/TI_headers.H index 443c266ac..f1eb5fbaf 100644 --- a/Source/TimeIntegration/TI_headers.H +++ b/Source/TimeIntegration/TI_headers.H @@ -14,7 +14,6 @@ */ void erf_slow_rhs_pre (int level, int nrk, amrex::Real dt, - amrex::BoxArray& grids_to_evolve, amrex::Vector& S_rhs, amrex::Vector& S_data, const amrex::MultiFab& S_prim, @@ -108,7 +107,6 @@ void erf_slow_rhs_post (int level, int nrk, * */ void erf_fast_rhs_N (int step, int level, - amrex::BoxArray& grids_to_evolve, amrex::Vector& S_slow_rhs, const amrex::Vector& S_prev, amrex::Vector& S_stage_data, @@ -130,7 +128,6 @@ void erf_fast_rhs_N (int step, int level, * */ void erf_fast_rhs_T (int step, int level, - amrex::BoxArray& grids_to_evolve, amrex::Vector& S_slow_rhs, const amrex::Vector& S_prev, amrex::Vector& S_stage_data, @@ -155,7 +152,6 @@ void erf_fast_rhs_T (int step, int level, * */ void erf_fast_rhs_MT (int step, int level, - amrex::BoxArray& grids_to_evolve, amrex::Vector& S_slow_rhs, const amrex::Vector& S_prev, amrex::Vector& S_stg_data, @@ -187,7 +183,6 @@ void erf_fast_rhs_MT (int step, int level, * integrator (the acoustic substepping). */ void make_fast_coeffs (int level, - amrex::BoxArray& grids_to_evolve, amrex::MultiFab& fast_coeffs, amrex::Vector& S_stage_data, const amrex::MultiFab& S_stage_prim, @@ -207,19 +202,18 @@ void make_fast_coeffs (int level, * equation for the z-component of momentum in the slow integrator. There * are three options for how buoyancy is computed (two are the same in the absence of moisture). */ -void make_buoyancy (amrex::BoxArray& grids_to_evolve, - amrex::Vector< amrex::MultiFab>& S_data, - const amrex::MultiFab & S_prim, - amrex::MultiFab& buoyancy, +void make_buoyancy (amrex::Vector< amrex::MultiFab>& S_data, + const amrex::MultiFab & S_prim, + amrex::MultiFab& buoyancy, #if defined(ERF_USE_MOISTURE) - const amrex::MultiFab& qmoist, - const amrex::Gpu::DeviceVector qv_d, - const amrex::Gpu::DeviceVector qc_d, - const amrex::Gpu::DeviceVector qi_d, + const amrex::MultiFab& qmoist, + const amrex::Gpu::DeviceVector qv_d, + const amrex::Gpu::DeviceVector qc_d, + const amrex::Gpu::DeviceVector qi_d, #endif - const amrex::Geometry geom, - const SolverChoice& solverChoice, - const amrex::MultiFab* r0); + const amrex::Geometry geom, + const SolverChoice& solverChoice, + const amrex::MultiFab* r0); #endif #ifdef ERF_USE_POISSON_SOLVE @@ -228,68 +222,67 @@ void make_buoyancy (amrex::BoxArray& grids_to_evolve, * */ void erf_slow_rhs_inc (int level, int nrk, - amrex::Real dt, - amrex::BoxArray& grids_to_evolve, - amrex::Vector& S_rhs, - amrex::Vector& S_old, - amrex::Vector& S_data, - const amrex::MultiFab& S_prim, - amrex::Vector& S_scratch, - const amrex::MultiFab& xvel, - const amrex::MultiFab& yvel, - const amrex::MultiFab& zvel, + amrex::Real dt, + amrex::Vector& S_rhs, + amrex::Vector& S_old, + amrex::Vector& S_data, + const amrex::MultiFab& S_prim, + amrex::Vector& S_scratch, + const amrex::MultiFab& xvel, + const amrex::MultiFab& yvel, + const amrex::MultiFab& zvel, #if defined(ERF_USE_MOISTURE) - const amrex::MultiFab& qvapor, + const amrex::MultiFab& qvapor, #endif - std::unique_ptr& z_t, - amrex::MultiFab& Omega, - const amrex::MultiFab& source, - const amrex::MultiFab& buoyancy, - amrex::MultiFab* Tau11, - amrex::MultiFab* Tau22, - amrex::MultiFab* Tau33, - amrex::MultiFab* Tau12, - amrex::MultiFab* Tau13, - amrex::MultiFab* Tau21, - amrex::MultiFab* Tau23, - amrex::MultiFab* Tau31, - amrex::MultiFab* Tau32, - amrex::MultiFab* SmnSmn, - amrex::MultiFab* eddyDiffs, - amrex::MultiFab* Hfx3, - amrex::MultiFab* Diss, - const amrex::Geometry geom, - const SolverChoice& solverChoice, - std::unique_ptr& most, - const amrex::Gpu::DeviceVector& domain_bcs_type_d, - const amrex::Vector& domain_bcs_type, - std::unique_ptr& z_phys_nd, - std::unique_ptr& dJ, - const amrex::MultiFab* p0, - std::unique_ptr& mapfac_m, - std::unique_ptr& mapfac_u, - std::unique_ptr& mapfac_v, - const amrex::Real* dptr_rayleigh_tau, - const amrex::Real* dptr_rayleigh_ubar, - const amrex::Real* dptr_rayleigh_vbar, - const amrex::Real* dptr_rayleigh_wbar, - const amrex::Real* dptr_rayleigh_thetabar); + std::unique_ptr& z_t, + amrex::MultiFab& Omega, + const amrex::MultiFab& source, + const amrex::MultiFab& buoyancy, + amrex::MultiFab* Tau11, + amrex::MultiFab* Tau22, + amrex::MultiFab* Tau33, + amrex::MultiFab* Tau12, + amrex::MultiFab* Tau13, + amrex::MultiFab* Tau21, + amrex::MultiFab* Tau23, + amrex::MultiFab* Tau31, + amrex::MultiFab* Tau32, + amrex::MultiFab* SmnSmn, + amrex::MultiFab* eddyDiffs, + amrex::MultiFab* Hfx3, + amrex::MultiFab* Diss, + const amrex::Geometry geom, + const SolverChoice& solverChoice, + std::unique_ptr& most, + const amrex::Gpu::DeviceVector& domain_bcs_type_d, + const amrex::Vector& domain_bcs_type, + std::unique_ptr& z_phys_nd, + std::unique_ptr& dJ, + const amrex::MultiFab* p0, + std::unique_ptr& mapfac_m, + std::unique_ptr& mapfac_u, + std::unique_ptr& mapfac_v, + const amrex::Real* dptr_rayleigh_tau, + const amrex::Real* dptr_rayleigh_ubar, + const amrex::Real* dptr_rayleigh_vbar, + const amrex::Real* dptr_rayleigh_wbar, + const amrex::Real* dptr_rayleigh_thetabar); #endif -void -ApplySpongeZoneBCs( - const SpongeChoice& spongeChoice, - const amrex::Geometry geom, - const amrex::Box& tbx, - const amrex::Box& tby, - const amrex::Box& tbz, - const amrex::Array4& rho_u_rhs, - const amrex::Array4& rho_v_rhs, - const amrex::Array4& rho_w_rhs, - const amrex::Array4& rho_u, - const amrex::Array4& rho_v, - const amrex::Array4& rho_w, - const amrex::Box& bx, - const amrex::Array4& cell_rhs, - const amrex::Array4& cell_data); +void ApplySpongeZoneBCs (const SpongeChoice& spongeChoice, + const amrex::Geometry geom, + const amrex::Box& tbx, + const amrex::Box& tby, + const amrex::Box& tbz, + const amrex::Array4& rho_u_rhs, + const amrex::Array4& rho_v_rhs, + const amrex::Array4& rho_w_rhs, + const amrex::Array4& rho_u, + const amrex::Array4& rho_v, + const amrex::Array4& rho_w, + const amrex::Box& bx, + const amrex::Array4& cell_rhs, + const amrex::Array4& cell_data); + + diff --git a/Source/TimeIntegration/TI_slow_rhs_fun.H b/Source/TimeIntegration/TI_slow_rhs_fun.H index f5c081f0e..08fc4ce2e 100644 --- a/Source/TimeIntegration/TI_slow_rhs_fun.H +++ b/Source/TimeIntegration/TI_slow_rhs_fun.H @@ -74,13 +74,13 @@ MultiFab* r0_new = &r_hse_new; MultiFab* p0_new = &p_hse_new; - make_buoyancy(grids_to_evolve[level], S_data, S_prim, buoyancy, + make_buoyancy(S_data, S_prim, buoyancy, #if defined(ERF_USE_MOISTURE) qmoist[level], qv_d, qc_d, qi_d, #endif fine_geom, solverChoice, r0_new); - erf_slow_rhs_pre(level, nrk, slow_dt, grids_to_evolve[level], S_rhs, S_data, S_prim, S_scratch, + erf_slow_rhs_pre(level, nrk, slow_dt, S_rhs, S_data, S_prim, S_scratch, xvel_new, yvel_new, zvel_new, #if defined(ERF_USE_MOISTURE) qmoist[level], @@ -170,13 +170,13 @@ } else { // If not moving_terrain - make_buoyancy(grids_to_evolve[level], S_data, S_prim, buoyancy, + make_buoyancy(S_data, S_prim, buoyancy, #if defined(ERF_USE_MOISTURE) qmoist[level], qv_d, qc_d, qi_d, #endif fine_geom, solverChoice, r0); - erf_slow_rhs_pre(level, nrk, slow_dt, grids_to_evolve[level], S_rhs, S_data, S_prim, S_scratch, + erf_slow_rhs_pre(level, nrk, slow_dt, S_rhs, S_data, S_prim, S_scratch, xvel_new, yvel_new, zvel_new, #if defined(ERF_USE_MOISTURE) qmoist[level], @@ -195,7 +195,7 @@ #ifdef ERF_USE_NETCDF // Populate RHS for relaxation zones if (init_type=="real" && level==0) { - wrfbdy_compute_interior_ghost_RHS(bdy_time_interval, start_bdy_time, new_stage_time, slow_dt, + wrfbdy_compute_interior_ghost_rhs(bdy_time_interval, start_bdy_time, new_stage_time, slow_dt, wrfbdy_width-1, wrfbdy_set_width, fine_geom, S_rhs, S_data, bdy_data_xlo, bdy_data_xhi, bdy_data_ylo, bdy_data_yhi); @@ -204,11 +204,10 @@ // Compute RHS for fine interior ghost if (level>0 && coupling_type=="OneWay" && cf_width>0) { - fine_compute_interior_ghost_RHS(new_stage_time, slow_dt, - cf_width-1, cf_set_width, fine_geom, + fine_compute_interior_ghost_rhs(new_stage_time, slow_dt, + cf_width, cf_set_width, fine_geom, &FPr_c[level-1], &FPr_u[level-1], &FPr_v[level-1], &FPr_w[level-1], - boxes_at_level[level], domain_bcs_type, - S_rhs, S_data); + domain_bcs_type, S_rhs, S_data); } // S_rhs[IntVar::cons].FillBoundary(fine_geom.periodicity()); @@ -318,13 +317,13 @@ Real slow_dt = new_stage_time - old_step_time; // If not moving_terrain - make_buoyancy(grids_to_evolve[level], S_data, S_prim, buoyancy, + make_buoyancy(S_data, S_prim, buoyancy, #if defined(ERF_USE_MOISTURE) qmoist[level], qv_d, qc_d, qi_d, #endif fine_geom, solverChoice, r0); - erf_slow_rhs_inc(level, nrk, slow_dt, grids_to_evolve[level], + erf_slow_rhs_inc(level, nrk, slow_dt, S_rhs, S_old, S_data, S_prim, S_scratch, xvel_new, yvel_new, zvel_new, #if defined(ERF_USE_MOISTURE) @@ -343,7 +342,7 @@ #ifdef ERF_USE_NETCDF // Populate RHS for relaxation zones if (init_type=="real" && level==0) { - wrfbdy_compute_interior_ghost_RHS(bdy_time_interval, start_bdy_time, new_stage_time, slow_dt, + wrfbdy_compute_interior_ghost_rhs(bdy_time_interval, start_bdy_time, new_stage_time, slow_dt, wrfbdy_width-1, wrfbdy_set_width, fine_geom, S_rhs, S_data, bdy_data_xlo, bdy_data_xhi, bdy_data_ylo, bdy_data_yhi); @@ -352,8 +351,8 @@ // Compute RHS for fine interior ghost if (level>0 && coupling_type=="OneWay" && cf_width>0) { - fine_compute_interior_ghost_RHS(new_stage_time, slow_dt, - cf_width-1, cf_set_width, + fine_compute_interior_ghost_rhs(new_stage_time, slow_dt, + cf_width, cf_set_width, &FPr_c[level-1], &FPr_u[level-1], &FPr_v[level-1], &FPr_w[level-1], boxes_at_level[level], domain_bcs_type, S_rhs, S_data); diff --git a/Source/Utils/InteriorGhostCells.cpp b/Source/Utils/InteriorGhostCells.cpp index b0ae3ff1e..9ad271b00 100644 --- a/Source/Utils/InteriorGhostCells.cpp +++ b/Source/Utils/InteriorGhostCells.cpp @@ -19,18 +19,17 @@ PhysBCFunctNoOp void_bc; * @param[in] ng_vect number of ghost cells in each direction * @param[in] get_int_ng flag to get ghost cells inside the domain */ - void -compute_interior_ghost_bxs_xy(const Box& bx, - const Box& domain, - const int& width, - const int& set_width, - Box& bx_xlo, - Box& bx_xhi, - Box& bx_ylo, - Box& bx_yhi, - const IntVect& ng_vect, - const bool get_int_ng) +compute_interior_ghost_bxs_xy (const Box& bx, + const Box& domain, + const int& width, + const int& set_width, + Box& bx_xlo, + Box& bx_xhi, + Box& bx_ylo, + Box& bx_yhi, + const IntVect& ng_vect, + const bool get_int_ng) { AMREX_ALWAYS_ASSERT(bx.ixType() == domain.ixType()); @@ -105,19 +104,19 @@ compute_interior_ghost_bxs_xy(const Box& bx, * @param[in] start_bdy_time time of the first boundary data read in */ void -wrfbdy_compute_interior_ghost_RHS(const Real& bdy_time_interval, - const Real& start_bdy_time, - const Real& time, - const Real& delta_t, - const int& width, - const int& set_width, - const Geometry& geom, - Vector& S_rhs, - Vector& S_data, - Vector>& bdy_data_xlo, - Vector>& bdy_data_xhi, - Vector>& bdy_data_ylo, - Vector>& bdy_data_yhi) +wrfbdy_compute_interior_ghost_rhs (const Real& bdy_time_interval, + const Real& start_bdy_time, + const Real& time, + const Real& delta_t, + const int& width, + const int& set_width, + const Geometry& geom, + Vector& S_rhs, + Vector& S_data, + Vector>& bdy_data_xlo, + Vector>& bdy_data_xhi, + Vector>& bdy_data_ylo, + Vector>& bdy_data_yhi) { BL_PROFILE_REGION("wrfbdy_compute_interior_ghost_RHS()"); @@ -389,7 +388,7 @@ wrfbdy_compute_interior_ghost_RHS(const Real& bdy_time_interval, compute_interior_ghost_bxs_xy(tbx, domain, set_width, 0, bx_xlo, bx_xhi, bx_ylo, bx_yhi); - zero_RHS_in_set_region(Rho_comp, 2, bx_xlo, bx_xhi, bx_ylo, bx_yhi, rhs_cons); + wrfbdy_zero_rhs_in_set_region(Rho_comp, 2, bx_xlo, bx_xhi, bx_ylo, bx_yhi, rhs_cons); } { @@ -400,7 +399,7 @@ wrfbdy_compute_interior_ghost_RHS(const Real& bdy_time_interval, compute_interior_ghost_bxs_xy(tbx, domain, set_width, 0, bx_xlo, bx_xhi, bx_ylo, bx_yhi); - zero_RHS_in_set_region(0, 1, bx_xlo, bx_xhi, bx_ylo, bx_yhi, rhs_xmom); + wrfbdy_zero_rhs_in_set_region(0, 1, bx_xlo, bx_xhi, bx_ylo, bx_yhi, rhs_xmom); } { @@ -411,7 +410,7 @@ wrfbdy_compute_interior_ghost_RHS(const Real& bdy_time_interval, compute_interior_ghost_bxs_xy(tbx, domain, set_width, 0, bx_xlo, bx_xhi, bx_ylo, bx_yhi); - zero_RHS_in_set_region(0, 1, bx_xlo, bx_xhi, bx_ylo, bx_yhi, rhs_ymom); + wrfbdy_zero_rhs_in_set_region(0, 1, bx_xlo, bx_xhi, bx_ylo, bx_yhi, rhs_ymom); } { @@ -422,7 +421,7 @@ wrfbdy_compute_interior_ghost_RHS(const Real& bdy_time_interval, compute_interior_ghost_bxs_xy(tbx, domain, set_width, 0, bx_xlo, bx_xhi, bx_ylo, bx_yhi); - zero_RHS_in_set_region(0, 1, bx_xlo, bx_xhi, bx_ylo, bx_yhi, rhs_zmom); + wrfbdy_zero_rhs_in_set_region(0, 1, bx_xlo, bx_xhi, bx_ylo, bx_yhi, rhs_zmom); } } // mfi @@ -480,10 +479,10 @@ wrfbdy_compute_interior_ghost_RHS(const Real& bdy_time_interval, continue; } - compute_Laplacian_relaxation(delta_t, icomp, 1, width, set_width, dom_lo, dom_hi, F1, F2, - tbx_xlo, tbx_xhi, tbx_ylo, tbx_yhi, - arr_xlo, arr_xhi, arr_ylo, arr_yhi, - data_arr, rhs_arr); + wrfbdy_compute_laplacian_relaxation(delta_t, icomp, 1, width, set_width, dom_lo, dom_hi, F1, F2, + tbx_xlo, tbx_xhi, tbx_ylo, tbx_yhi, + arr_xlo, arr_xhi, arr_ylo, arr_yhi, + data_arr, rhs_arr); } // mfi } // ivar } @@ -505,19 +504,18 @@ wrfbdy_compute_interior_ghost_RHS(const Real& bdy_time_interval, * @param[in] S_data current value of the solution */ void -fine_compute_interior_ghost_RHS(const Real& time, - const Real& delta_t, - const int& width, - const int& set_width, - const Geometry& geom, - ERFFillPatcher* FPr_c, - ERFFillPatcher* FPr_u, - ERFFillPatcher* FPr_v, - ERFFillPatcher* FPr_w, - Vector& boxes_at_level, - Vector& domain_bcs_type, - Vector& S_rhs_f, - Vector& S_data_f) +fine_compute_interior_ghost_rhs (const Real& time, + const Real& delta_t, + const int& width, + const int& set_width, + const Geometry& geom, + ERFFillPatcher* FPr_c, + ERFFillPatcher* FPr_u, + ERFFillPatcher* FPr_v, + ERFFillPatcher* FPr_w, + Vector& domain_bcs_type, + Vector& S_rhs_f, + Vector& S_data_f) { BL_PROFILE_REGION("fine_compute_interior_ghost_RHS()"); @@ -552,152 +550,99 @@ fine_compute_interior_ghost_RHS(const Real& time, MultiFab& fmf_p = fmf_p_v[ivar_idx]; amrex::MultiFab::Copy(fmf_p,fmf, 0, 0, num_var, fmf.nGrowVect()); + // Integer mask MF + int set_mask_val; + int relax_mask_val; + iMultiFab* mask; + // Fill fine patch on interior halo region //========================================================== if (ivar_idx == IntVar::cons) { - FPr_c->fill(fmf_p, time, void_bc, domain_bcs_type); + FPr_c->FillRelax(fmf_p, time, void_bc, domain_bcs_type); + mask = FPr_c->GetMask(); + set_mask_val = FPr_c->GetSetMaskVal(); + relax_mask_val = FPr_c->GetRelaxMaskVal(); } else if (ivar_idx == IntVar::xmom) { - FPr_u->fill(fmf_p, time, void_bc, domain_bcs_type); + FPr_u->FillRelax(fmf_p, time, void_bc, domain_bcs_type); + mask = FPr_u->GetMask(); + set_mask_val = FPr_u->GetSetMaskVal(); + relax_mask_val = FPr_u->GetRelaxMaskVal(); #ifdef _OPENMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) #endif for ( MFIter mfi(fmf_p,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { - Box tbx = mfi.tilebox(); + Box vbx = mfi.validbox(); + const Array4& prim_arr = fmf_p.array(mfi); - const Array4& rho_arr = fmf_p_v[0].const_array(mfi); + const Array4& rho_arr = fmf_p_v[0].const_array(mfi); + const Array4& mask_arr = mask->const_array(mfi); - for (int g_ind(0); g_indfill(fmf_p, time, void_bc, domain_bcs_type); + FPr_v->FillRelax(fmf_p, time, void_bc, domain_bcs_type); + mask = FPr_v->GetMask(); + set_mask_val = FPr_v->GetSetMaskVal(); + relax_mask_val = FPr_v->GetRelaxMaskVal(); #ifdef _OPENMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) #endif for ( MFIter mfi(fmf_p,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { - Box tbx = mfi.tilebox(); + Box vbx = mfi.validbox(); + const Array4& prim_arr = fmf_p.array(mfi); - const Array4& rho_arr = fmf_p_v[0].const_array(mfi); + const Array4& rho_arr = fmf_p_v[0].const_array(mfi); + const Array4& mask_arr = mask->const_array(mfi); - for (int g_ind(0); g_indfill(fmf_p, time, void_bc, domain_bcs_type); + FPr_w->FillRelax(fmf_p, time, void_bc, domain_bcs_type); + mask = FPr_w->GetMask(); + set_mask_val = FPr_w->GetSetMaskVal(); + relax_mask_val = FPr_w->GetRelaxMaskVal(); #ifdef _OPENMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) #endif for ( MFIter mfi(fmf_p,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { - Box tbx = mfi.tilebox(); + Box vbx = mfi.validbox(); + const Array4& prim_arr = fmf_p.array(mfi); - const Array4& rho_arr = fmf_p_v[0].const_array(mfi); + const Array4& rho_arr = fmf_p_v[0].const_array(mfi); + const Array4& mask_arr = mask->const_array(mfi); - for (int g_ind(0); g_ind& rhs_arr = rhs.array(mfi); - for (int g_ind(0); g_ind& mask_arr = mask->const_array(mfi); - zero_RHS_in_set_region(0, num_var, bx_xlo, bx_xhi, bx_ylo, bx_yhi, rhs_arr); - } // g_ind + amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + if (mask_arr(i,j,k) == set_mask_val) { + rhs_arr(i,j,k) = 0.0; + } + }); } // mfi // For Laplacian stencil rhs.FillBoundary(geom.periodicity()); + // Compute RHS in relaxation region //========================================================== #ifdef _OPENMP @@ -737,98 +679,118 @@ fine_compute_interior_ghost_RHS(const Real& time, #endif for ( MFIter mfi(fmf_p,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { - Box tbx = mfi.tilebox(); + Box vbx = mfi.validbox(); + const Array4& rhs_arr = rhs.array(mfi); const Array4& fine_arr = fmf_p.const_array(mfi); const Array4& data_arr = fmf.const_array(mfi); - const Array4& rhs_arr = rhs.array(mfi); - - for (int g_ind(0); g_ind& mask_arr = mask->const_array(mfi); + + const auto& vbx_lo = lbound(vbx); + const auto& vbx_hi = ubound(vbx); + + int icomp = 0; + + int Spec_z = set_width; + int Relax_z = width - Spec_z; + amrex::Real num = amrex::Real(Spec_z + Relax_z); + amrex::Real denom = amrex::Real(Relax_z - 1); + amrex::ParallelFor(vbx, num_var, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + { + if (mask_arr(i,j,k) == relax_mask_val) { + + // Indices + Real n_ind; + int ii{width-1}; int jj{width-1}; + bool near_x_lo_wall{false}; bool near_x_hi_wall{false}; + bool near_y_lo_wall{false}; bool near_y_hi_wall{false}; + bool mask_x_found{false}; bool mask_y_found{false}; + + // Near x-wall + if ((i-vbx_lo.x) < width) { + near_x_lo_wall = true; + ii = i-vbx_lo.x; + if (mask_arr(vbx_lo.x,j,k) == 2) mask_x_found = true; + } else if ((vbx_hi.x-i) < width) { + near_x_hi_wall = true; + ii = vbx_hi.x-i; + if (mask_arr(vbx_hi.x,j,k) == 2) mask_x_found = true; + } + + // Near y-wall + if ((j-vbx_lo.y) < width) { + near_y_lo_wall = true; + jj = j-vbx_lo.y; + if (mask_arr(i,vbx_lo.y,k) == 2) mask_y_found = true; + } else if ((vbx_hi.y-j) < width) { + near_y_hi_wall = true; + jj = vbx_hi.y-j; + if (mask_arr(i,vbx_hi.y,k) == 2) mask_y_found = true; + } + + // Found a nearby masked cell (valid n_ind) + if (mask_x_found && mask_y_found) { + n_ind = std::min(ii,jj) + 1.0; + } else if (mask_x_found) { + n_ind = ii + 1.0; + } else if (mask_y_found) { + n_ind = jj + 1.0; + // Pesky corner cell + } else { + if (near_x_lo_wall || near_x_hi_wall) { + Real dj_min{width-1.0}; + int j_lb = std::max(vbx_lo.y,j-width); + int j_ub = std::min(vbx_hi.y,j+width); + int li = (near_x_lo_wall) ? vbx_lo.x : vbx_hi.x; + for (int lj(j_lb); lj<=j_ub; ++lj) { + if (mask_arr(li,lj,k) == 2) { + mask_y_found = true; + dj_min = std::min(dj_min,(Real) std::abs(lj-j)); + } + } + if (mask_y_found) { + Real mag = sqrt( amrex::Real(dj_min*dj_min + ii*ii) ); + n_ind = std::min(mag,width-1.0) + 1.0; + } else { + amrex::Abort("Mask not found near x wall!"); + } + } else if (near_y_lo_wall || near_y_hi_wall) { + Real di_min{width-1.0}; + int i_lb = std::max(vbx_lo.x,i-width); + int i_ub = std::min(vbx_hi.x,i+width); + int lj = (near_y_lo_wall) ? vbx_lo.y : vbx_hi.y; + for (int li(i_lb); li<=i_ub; ++li) { + if (mask_arr(li,lj,k) == 2) { + mask_x_found = true; + di_min = std::min(di_min,(Real) std::abs(li-i)); + } + } + if (mask_x_found) { + Real mag = sqrt( amrex::Real(di_min*di_min + jj*jj) ); + n_ind = std::min(mag,width-1.0) + 1.0; + } else { + amrex::Abort("Mask not found near y wall!"); + } + } else { + amrex::Abort("Relaxation cell must be near a wall!"); + } + } + + amrex::Real Factor = (num - n_ind)/denom; + amrex::Real d = data_arr(i ,j ,k ,n+icomp) + delta_t*rhs_arr(i , j , k ,n+icomp); + amrex::Real d_ip1 = data_arr(i+1,j ,k ,n+icomp) + delta_t*rhs_arr(i+1, j , k ,n+icomp); + amrex::Real d_im1 = data_arr(i-1,j ,k ,n+icomp) + delta_t*rhs_arr(i-1, j , k ,n+icomp); + amrex::Real d_jp1 = data_arr(i ,j+1,k ,n+icomp) + delta_t*rhs_arr(i , j+1, k ,n+icomp); + amrex::Real d_jm1 = data_arr(i ,j-1,k ,n+icomp) + delta_t*rhs_arr(i , j-1, k ,n+icomp); + amrex::Real delta = fine_arr(i ,j ,k,n) - d; + amrex::Real delta_xp = fine_arr(i+1,j ,k,n) - d_ip1; + amrex::Real delta_xm = fine_arr(i-1,j ,k,n) - d_im1; + amrex::Real delta_yp = fine_arr(i ,j+1,k,n) - d_jp1; + amrex::Real delta_ym = fine_arr(i ,j-1,k,n) - d_jm1; + amrex::Real Laplacian = delta_xp + delta_xm + delta_yp + delta_ym - 4.0*delta; + rhs_arr(i,j,k,n) += (F1*delta - F2*Laplacian) * Factor; + } + }); } // mfi } // ivar_idx } -/** - * Update the solution in the relaxation zone - * - * @param[in] delta_t timestep - * @param[in] width number of cells in (relaxation+specified) zone - * @param[in] geom container for geometric information - * @param[in] S_rhs RHS to be added here - * @param[in] S_old previous value of the solution - * @param[out] S_data new value of the solution defined here - */ - -void -update_interior_ghost(const Real& delta_t, - const int& width, - Vector& boxes_at_level, - Vector& S_rhs, - const Vector& S_old, - Vector& S_data) -{ - BL_PROFILE_REGION("update_interior_ghost()"); - - Vector ncomp_map = {2, 1, 1, 1}; - for (int ivar(IntVar::cons); ivar < IntVar::NumVars; ivar++) - { - int ncomp = ncomp_map[ivar]; - IndexType m_ixt = S_data[ivar].boxArray().ixType(); - -#ifdef _OPENMP -#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) -#endif - for ( MFIter mfi(S_data[ivar],amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - Box tbx = mfi.tilebox(); - Array4 rhs_arr = S_rhs[ivar].array(mfi); - Array4 old_arr = S_old[ivar].array(mfi); - Array4 new_arr = S_data[ivar].array(mfi); - - for (int g_ind(0); g_ind& S_rhs, + amrex::Vector& S_data, + amrex::Vector>& bdy_data_xlo, + amrex::Vector>& bdy_data_xhi, + amrex::Vector>& bdy_data_ylo, + amrex::Vector>& bdy_data_yhi); + +/* + * Compute relaxation region RHS at fine-crse interface + */ +void +fine_compute_interior_ghost_rhs (const amrex::Real& time, + const amrex::Real& delta_t, + const int& width, + const int& set_width, + const amrex::Geometry& geom, + ERFFillPatcher* FPr_c, + ERFFillPatcher* FPr_u, + ERFFillPatcher* FPr_v, + ERFFillPatcher* FPr_w, + amrex::Vector& domain_bcs_type, + amrex::Vector& S_rhs_f, + amrex::Vector& S_data_f); /** * Zero RHS in the set region @@ -80,13 +113,13 @@ void compute_interior_ghost_bxs_xy (const amrex::Box& bx, AMREX_GPU_HOST AMREX_FORCE_INLINE void -zero_RHS_in_set_region (const int& icomp, - const int& num_var, - const amrex::Box& bx_xlo, - const amrex::Box& bx_xhi, - const amrex::Box& bx_ylo, - const amrex::Box& bx_yhi, - const amrex::Array4& rhs_arr) +wrfbdy_zero_rhs_in_set_region (const int& icomp, + const int& num_var, + const amrex::Box& bx_xlo, + const amrex::Box& bx_xhi, + const amrex::Box& bx_ylo, + const amrex::Box& bx_yhi, + const amrex::Array4& rhs_arr) { amrex::ParallelFor(bx_xlo, num_var, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { @@ -107,7 +140,6 @@ zero_RHS_in_set_region (const int& icomp, }); } - /** * Compute the Laplacian RHS in the relaxation zone * @@ -134,25 +166,25 @@ zero_RHS_in_set_region (const int& icomp, AMREX_GPU_HOST AMREX_FORCE_INLINE void -compute_Laplacian_relaxation (const amrex::Real& delta_t, - const int& icomp, - const int& num_var, - const int& width, - const int& set_width, - const amrex::Dim3& dom_lo, - const amrex::Dim3& dom_hi, - const amrex::Real& F1, - const amrex::Real& F2, - const amrex::Box& bx_xlo, - const amrex::Box& bx_xhi, - const amrex::Box& bx_ylo, - const amrex::Box& bx_yhi, - const amrex::Array4& arr_xlo, - const amrex::Array4& arr_xhi, - const amrex::Array4& arr_ylo, - const amrex::Array4& arr_yhi, - const amrex::Array4& data_arr, - const amrex::Array4& rhs_arr) +wrfbdy_compute_laplacian_relaxation (const amrex::Real& delta_t, + const int& icomp, + const int& num_var, + const int& width, + const int& set_width, + const amrex::Dim3& dom_lo, + const amrex::Dim3& dom_hi, + const amrex::Real& F1, + const amrex::Real& F2, + const amrex::Box& bx_xlo, + const amrex::Box& bx_xhi, + const amrex::Box& bx_ylo, + const amrex::Box& bx_yhi, + const amrex::Array4& arr_xlo, + const amrex::Array4& arr_xhi, + const amrex::Array4& arr_ylo, + const amrex::Array4& arr_yhi, + const amrex::Array4& data_arr, + const amrex::Array4& rhs_arr) { // RHS computation int Spec_z = set_width; @@ -255,69 +287,4 @@ compute_Laplacian_relaxation (const amrex::Real& delta_t, } }); } - -/* - * Compute relaxation region RHS with wrfbdy - */ -void wrfbdy_compute_interior_ghost_RHS (const amrex::Real& bdy_time_interval, - const amrex::Real& start_bdy_time, - const amrex::Real& time, - const amrex::Real& delta_t, - const int& width, - const int& set_width, - const amrex::Geometry& geom, - amrex::Vector& S_rhs, - amrex::Vector& S_data, - amrex::Vector>& bdy_data_xlo, - amrex::Vector>& bdy_data_xhi, - amrex::Vector>& bdy_data_ylo, - amrex::Vector>& bdy_data_yhi); - -/* - * Compute moist relaxation region RHS with wrfbdy - */ -#if defined(ERF_USE_NETCDF) && (defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP)) -void wrfbdy_compute_moist_interior_ghost_RHS (amrex::Box tbx, - const amrex::Real& bdy_time_interval, - const amrex::Real& start_bdy_time, - const amrex::Real& time, - const amrex::Real& delta_t, - const int& width, - const int& set_width, - const amrex::Geometry& geom, - const amrex::Array4& rhs_arr, - const amrex::Array4& data_arr, - amrex::Vector>& bdy_data_xlo, - amrex::Vector>& bdy_data_xhi, - amrex::Vector>& bdy_data_ylo, - amrex::Vector>& bdy_data_yhi); -#endif - -/* - * Compute relaxation region RHS at fine-crse interface - */ -void -fine_compute_interior_ghost_RHS (const amrex::Real& time, - const amrex::Real& delta_t, - const int& width, - const int& set_width, - const amrex::Geometry& geom, - ERFFillPatcher* FPr_c, - ERFFillPatcher* FPr_u, - ERFFillPatcher* FPr_v, - ERFFillPatcher* FPr_w, - amrex::Vector& boxes_at_level, - amrex::Vector& domain_bcs_type, - amrex::Vector& S_rhs_f, - amrex::Vector& S_data_f); - -/* - * Update the variables in the relaxation/set region - */ -void update_interior_ghost (const amrex::Real& delta_t, - const int& width, - amrex::Vector& boxes_at_level, - amrex::Vector& S_rhs, - const amrex::Vector& S_old, - amrex::Vector& S_data); #endif diff --git a/Source/main.cpp b/Source/main.cpp index 74432a4f7..16cbb229c 100644 --- a/Source/main.cpp +++ b/Source/main.cpp @@ -37,7 +37,9 @@ void add_par () { pp.add("n_proper",2); pp.add("max_grid_size",2048); - pp.add("blocking_factor",1); + + int bf=1; + pp.queryAdd("blocking_factor",bf); pp.add("n_error_buf",0); }