diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake index e386237c0..4b0dc4cee 100644 --- a/CMake/BuildERFExe.cmake +++ b/CMake/BuildERFExe.cmake @@ -153,6 +153,7 @@ function(build_erf_lib erf_lib_name) ${SRC_DIR}/IO/console_io.cpp ${SRC_DIR}/SourceTerms/ERF_ApplySpongeZoneBCs.cpp ${SRC_DIR}/SourceTerms/ERF_make_buoyancy.cpp + ${SRC_DIR}/SourceTerms/ERF_add_thin_body_sources.cpp ${SRC_DIR}/SourceTerms/ERF_make_mom_sources.cpp ${SRC_DIR}/SourceTerms/ERF_make_sources.cpp ${SRC_DIR}/SourceTerms/ERF_moist_set_rhs.cpp diff --git a/Source/BoundaryConditions/ERF_FillPatch.cpp b/Source/BoundaryConditions/ERF_FillPatch.cpp index dfd6b6422..d0b786c31 100644 --- a/Source/BoundaryConditions/ERF_FillPatch.cpp +++ b/Source/BoundaryConditions/ERF_FillPatch.cpp @@ -561,7 +561,9 @@ ERF::FillBdyCCVels (Vector& mf_cc_vel) for (MFIter mfi(mf_cc_vel[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi) { - const Box& bx = mfi.growntilebox(1); + // Note that we don't fill corners here -- only the cells that share a face + // with interior cells -- this is all that is needed to calculate vorticity + const Box& bx = mfi.tilebox(); const Array4& vel_arr = mf_cc_vel[lev].array(mfi); if (!Geom(lev).isPeriodic(0)) { diff --git a/Source/Derive.H b/Source/Derive.H index d93077a54..f40bdc432 100644 --- a/Source/Derive.H +++ b/Source/Derive.H @@ -122,5 +122,15 @@ void erf_dervortz ( const int* bcrec, const int level); +void erf_dermagvel ( + const amrex::Box& bx, + amrex::FArrayBox& derfab, + int dcomp, + int ncomp, + const amrex::FArrayBox& datfab, + const amrex::Geometry& geomdata, + amrex::Real time, + const int* bcrec, + const int level); } #endif diff --git a/Source/Derive.cpp b/Source/Derive.cpp index 8a6fdf006..6df9c5b33 100644 --- a/Source/Derive.cpp +++ b/Source/Derive.cpp @@ -278,5 +278,31 @@ erf_dervortz ( }); } +void +erf_dermagvel ( + const amrex::Box& bx, + amrex::FArrayBox& derfab, + int dcomp, + int ncomp, + const amrex::FArrayBox& datfab, + const amrex::Geometry& /*geomdata*/, + amrex::Real /*time*/, + const int* /*bcrec*/, + const int /*level*/) +{ + AMREX_ALWAYS_ASSERT(dcomp == 0); + AMREX_ALWAYS_ASSERT(ncomp == 1); + + auto const dat = datfab.array(); // cell-centered velocity + auto tfab = derfab.array(); // cell-centered magvel + + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + Real u = dat(i,j,k,0); + Real v = dat(i,j,k,1); + Real w = dat(i,j,k,2); + tfab(i,j,k,dcomp) = std::sqrt(u*u + v*v + w*w); + }); +} } // namespace diff --git a/Source/ERF.H b/Source/ERF.H index e699fa4e9..117aa9b80 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -767,7 +767,7 @@ private: // Note that the order of variable names here must match the order in Derive.cpp const amrex::Vector derived_names {"soundspeed", "temp", "theta", "KE", "QKE", "scalar", "vorticity_x","vorticity_y","vorticity_z", - "divU", + "magvel", "divU", "pres_hse", "dens_hse", "pressure", "pert_pres", "pert_dens", "eq_pot_temp", "num_turb", "dpdx", "dpdy", "pres_hse_x", "pres_hse_y", "z_phys", "detJ" , "mapfac", diff --git a/Source/IO/Plotfile.cpp b/Source/IO/Plotfile.cpp index 5a0c5cc42..9e9dc184d 100644 --- a/Source/IO/Plotfile.cpp +++ b/Source/IO/Plotfile.cpp @@ -198,6 +198,7 @@ ERF::WritePlotFile (int which, Vector plot_var_names) if (containerHasElement(plot_var_names, "x_velocity" ) || containerHasElement(plot_var_names, "y_velocity" ) || containerHasElement(plot_var_names, "z_velocity" ) || + containerHasElement(plot_var_names, "magvel" ) || containerHasElement(plot_var_names, "vorticity_x") || containerHasElement(plot_var_names, "vorticity_y") || containerHasElement(plot_var_names, "vorticity_z") ) { @@ -294,6 +295,7 @@ ERF::WritePlotFile (int which, Vector plot_var_names) calculate_derived("vorticity_x", mf_cc_vel[lev] , derived::erf_dervortx); calculate_derived("vorticity_y", mf_cc_vel[lev] , derived::erf_dervorty); calculate_derived("vorticity_z", mf_cc_vel[lev] , derived::erf_dervortz); + calculate_derived("magvel" , mf_cc_vel[lev] , derived::erf_dermagvel); if (containerHasElement(plot_var_names, "divU")) { @@ -1256,10 +1258,10 @@ ERF::WriteMultiLevelPlotfileWithTerrain (const std::string& plotfilename, int nl { BL_PROFILE("WriteMultiLevelPlotfileWithTerrain()"); - BL_ASSERT(nlevels <= mf.size()); - BL_ASSERT(nlevels <= ref_ratio.size()+1); - BL_ASSERT(nlevels <= level_steps.size()); - BL_ASSERT(mf[0]->nComp() == varnames.size()); + AMREX_ALWAYS_ASSERT(nlevels <= mf.size()); + AMREX_ALWAYS_ASSERT(nlevels <= ref_ratio.size()+1); + AMREX_ALWAYS_ASSERT(nlevels <= level_steps.size()); + AMREX_ALWAYS_ASSERT(mf[0]->nComp() == varnames.size()); bool callBarrier(false); PreBuildDirectorHierarchy(plotfilename, levelPrefix, nlevels, callBarrier); @@ -1338,78 +1340,78 @@ ERF::WriteGenericPlotfileHeaderWithTerrain (std::ostream &HeaderFile, const std::string &levelPrefix, const std::string &mfPrefix) const { - BL_ASSERT(nlevels <= bArray.size()); - BL_ASSERT(nlevels <= ref_ratio.size()+1); - BL_ASSERT(nlevels <= level_steps.size()); + AMREX_ALWAYS_ASSERT(nlevels <= bArray.size()); + AMREX_ALWAYS_ASSERT(nlevels <= ref_ratio.size()+1); + AMREX_ALWAYS_ASSERT(nlevels <= level_steps.size()); - HeaderFile.precision(17); + HeaderFile.precision(17); - // ---- this is the generic plot file type name - HeaderFile << versionName << '\n'; + // ---- this is the generic plot file type name + HeaderFile << versionName << '\n'; - HeaderFile << varnames.size() << '\n'; + HeaderFile << varnames.size() << '\n'; - for (int ivar = 0; ivar < varnames.size(); ++ivar) { - HeaderFile << varnames[ivar] << "\n"; - } - HeaderFile << AMREX_SPACEDIM << '\n'; - HeaderFile << time << '\n'; - HeaderFile << finest_level << '\n'; - for (int i = 0; i < AMREX_SPACEDIM; ++i) { - HeaderFile << geom[0].ProbLo(i) << ' '; - } - HeaderFile << '\n'; - for (int i = 0; i < AMREX_SPACEDIM; ++i) { - HeaderFile << geom[0].ProbHi(i) << ' '; - } - HeaderFile << '\n'; - for (int i = 0; i < finest_level; ++i) { - HeaderFile << ref_ratio[i][0] << ' '; - } - HeaderFile << '\n'; - for (int i = 0; i <= finest_level; ++i) { - HeaderFile << geom[i].Domain() << ' '; - } - HeaderFile << '\n'; - for (int i = 0; i <= finest_level; ++i) { - HeaderFile << level_steps[i] << ' '; + for (int ivar = 0; ivar < varnames.size(); ++ivar) { + HeaderFile << varnames[ivar] << "\n"; + } + HeaderFile << AMREX_SPACEDIM << '\n'; + HeaderFile << time << '\n'; + HeaderFile << finest_level << '\n'; + for (int i = 0; i < AMREX_SPACEDIM; ++i) { + HeaderFile << geom[0].ProbLo(i) << ' '; + } + HeaderFile << '\n'; + for (int i = 0; i < AMREX_SPACEDIM; ++i) { + HeaderFile << geom[0].ProbHi(i) << ' '; + } + HeaderFile << '\n'; + for (int i = 0; i < finest_level; ++i) { + HeaderFile << ref_ratio[i][0] << ' '; + } + HeaderFile << '\n'; + for (int i = 0; i <= finest_level; ++i) { + HeaderFile << geom[i].Domain() << ' '; + } + HeaderFile << '\n'; + for (int i = 0; i <= finest_level; ++i) { + HeaderFile << level_steps[i] << ' '; + } + HeaderFile << '\n'; + for (int i = 0; i <= finest_level; ++i) { + for (int k = 0; k < AMREX_SPACEDIM; ++k) { + HeaderFile << geom[i].CellSize()[k] << ' '; } HeaderFile << '\n'; - for (int i = 0; i <= finest_level; ++i) { - for (int k = 0; k < AMREX_SPACEDIM; ++k) { - HeaderFile << geom[i].CellSize()[k] << ' '; - } - HeaderFile << '\n'; - } - HeaderFile << (int) geom[0].Coord() << '\n'; - HeaderFile << "0\n"; + } + HeaderFile << (int) geom[0].Coord() << '\n'; + HeaderFile << "0\n"; - for (int level = 0; level <= finest_level; ++level) { - HeaderFile << level << ' ' << bArray[level].size() << ' ' << time << '\n'; - HeaderFile << level_steps[level] << '\n'; + for (int level = 0; level <= finest_level; ++level) { + HeaderFile << level << ' ' << bArray[level].size() << ' ' << time << '\n'; + HeaderFile << level_steps[level] << '\n'; - const IntVect& domain_lo = geom[level].Domain().smallEnd(); - for (int i = 0; i < bArray[level].size(); ++i) - { - // Need to shift because the RealBox ctor we call takes the - // physical location of index (0,0,0). This does not affect - // the usual cases where the domain index starts with 0. - const Box& b = shift(bArray[level][i], -domain_lo); - RealBox loc = RealBox(b, geom[level].CellSize(), geom[level].ProbLo()); - for (int n = 0; n < AMREX_SPACEDIM; ++n) { - HeaderFile << loc.lo(n) << ' ' << loc.hi(n) << '\n'; - } + const IntVect& domain_lo = geom[level].Domain().smallEnd(); + for (int i = 0; i < bArray[level].size(); ++i) + { + // Need to shift because the RealBox ctor we call takes the + // physical location of index (0,0,0). This does not affect + // the usual cases where the domain index starts with 0. + const Box& b = shift(bArray[level][i], -domain_lo); + RealBox loc = RealBox(b, geom[level].CellSize(), geom[level].ProbLo()); + for (int n = 0; n < AMREX_SPACEDIM; ++n) { + HeaderFile << loc.lo(n) << ' ' << loc.hi(n) << '\n'; } - - HeaderFile << MultiFabHeaderPath(level, levelPrefix, mfPrefix) << '\n'; - } - HeaderFile << "1" << "\n"; - HeaderFile << "3" << "\n"; - HeaderFile << "amrexvec_nu_x" << "\n"; - HeaderFile << "amrexvec_nu_y" << "\n"; - HeaderFile << "amrexvec_nu_z" << "\n"; - std::string mf_nodal_prefix = "Nu_nd"; - for (int level = 0; level <= finest_level; ++level) { - HeaderFile << MultiFabHeaderPath(level, levelPrefix, mf_nodal_prefix) << '\n'; } + + HeaderFile << MultiFabHeaderPath(level, levelPrefix, mfPrefix) << '\n'; + } + HeaderFile << "1" << "\n"; + HeaderFile << "3" << "\n"; + HeaderFile << "amrexvec_nu_x" << "\n"; + HeaderFile << "amrexvec_nu_y" << "\n"; + HeaderFile << "amrexvec_nu_z" << "\n"; + std::string mf_nodal_prefix = "Nu_nd"; + for (int level = 0; level <= finest_level; ++level) { + HeaderFile << MultiFabHeaderPath(level, levelPrefix, mf_nodal_prefix) << '\n'; + } } diff --git a/Source/SourceTerms/ERF_add_thin_body_sources.cpp b/Source/SourceTerms/ERF_add_thin_body_sources.cpp new file mode 100644 index 000000000..eb7b210d1 --- /dev/null +++ b/Source/SourceTerms/ERF_add_thin_body_sources.cpp @@ -0,0 +1,131 @@ +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +using namespace amrex; + +/** + * Function for computing the slow RHS for the evolution equations for the density, potential temperature and momentum. + * + * @param[in] xmom_src source terms for x-momentum + * @param[in] ymom_src source terms for y-momentum + * @param[in] zmom_src source terms for z-momentum + * @param[in] xflux_imask_lev thin-body mask on x-faces + * @param[in] yflux_imask_lev thin-body mask on y-faces + * @param[in] zflux_imask_lev thin-body mask on z-faces + * @param[in] thin_xforce_lev x-component of forces on thin immersed bodies + * @param[in] thin_yforce_lev y-component of forces on thin immersed bodies + * @param[in] thin_zforce_lev z-component of forces on thin immersed bodies + */ + +void add_thin_body_sources ( MultiFab & xmom_src, + MultiFab & ymom_src, + MultiFab & zmom_src, + std::unique_ptr& xflux_imask_lev, + std::unique_ptr& yflux_imask_lev, + std::unique_ptr& zflux_imask_lev, + std::unique_ptr& thin_xforce_lev, + std::unique_ptr& thin_yforce_lev, + std::unique_ptr& thin_zforce_lev) +{ + BL_PROFILE_REGION("erf_add_thin_body_sources()"); + + const bool l_have_thin_xforce = (thin_xforce_lev != nullptr); + const bool l_have_thin_yforce = (thin_yforce_lev != nullptr); + const bool l_have_thin_zforce = (thin_zforce_lev != nullptr); + + // ***************************************************************************** + // If a thin immersed body is present, add forcing terms + // ***************************************************************************** + if (l_have_thin_xforce) { + MultiFab::Copy(*thin_xforce_lev, xmom_src, 0, 0, 1, 0); + thin_xforce_lev->mult(-1., 0, 1, 0); + ApplyInvertedMask(*thin_xforce_lev, *xflux_imask_lev, 0); + MultiFab::Add(xmom_src, *thin_xforce_lev, 0, 0, 1, 0); + } + + if (l_have_thin_yforce) { + MultiFab::Copy(*thin_yforce_lev, ymom_src, 0, 0, 1, 0); + thin_yforce_lev->mult(-1., 0, 1, 0); + ApplyInvertedMask(*thin_yforce_lev, *yflux_imask_lev, 0); + MultiFab::Add(ymom_src, *thin_yforce_lev, 0, 0, 1, 0); + } + + if (l_have_thin_zforce) { + MultiFab::Copy(*thin_zforce_lev, zmom_src, 0, 0, 1, 0); + thin_zforce_lev->mult(-1., 0, 1, 0); + ApplyInvertedMask(*thin_zforce_lev, *zflux_imask_lev, 0); + MultiFab::Add(zmom_src, *thin_zforce_lev, 0, 0, 1, 0); + } + +#if 0 +#ifndef AMREX_USE_GPU + if (l_have_thin_xforce) { + // TODO: Implement particles to better track and output these data + if (nrk==2) { + for ( MFIter mfi(S_data[IntVars::cons],TileNoZ()); mfi.isValid(); ++mfi) + { + const Box& tbx = mfi.nodaltilebox(0); + const Array4 & fx = thin_xforce_lev->const_array(mfi); + const Array4 & mask = xflux_imask_lev->const_array(mfi); + ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { + if (mask(i,j,k)==0) { + amrex::AllPrint() << "thin body fx"< & fy = thin_yforce_lev->const_array(mfi); + const Array4 & mask = yflux_imask_lev->const_array(mfi); + ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k) { + if (mask(i,j,k)==0) { + amrex::AllPrint() << "thin body fy"< & fz = thin_zforce_lev->const_array(mfi); + const Array4 & mask = zflux_imask_lev->const_array(mfi); + ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { + if (mask(i,j,k)==0) { + amrex::AllPrint() << "thin body fz"< d_rayleigh_ptrs_at_lev, const int n_qstate); +void add_thin_body_sources (amrex::MultiFab& xmom_source, + amrex::MultiFab& ymom_source, + amrex::MultiFab& zmom_source, + std::unique_ptr& xflux_imask_lev, + std::unique_ptr& yflux_imask_lev, + std::unique_ptr& zflux_imask_lev, + std::unique_ptr& thin_xforce_lev, + std::unique_ptr& thin_yforce_lev, + std::unique_ptr& thin_zforce_lev); + #if defined(ERF_USE_NETCDF) void moist_set_rhs ( diff --git a/Source/TimeIntegration/ERF_slow_rhs_inc.cpp b/Source/TimeIntegration/ERF_slow_rhs_inc.cpp index 924ed4675..81b111876 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_inc.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_inc.cpp @@ -55,12 +55,6 @@ using namespace amrex; * @param[in] mapfac_m map factor at cell centers * @param[in] mapfac_u map factor at x-faces * @param[in] mapfac_v map factor at y-faces - * @param[in] xflux_imask_lev thin-body mask on x-faces - * @param[in] yflux_imask_lev thin-body mask on y-faces - * @param[in] zflux_imask_lev thin-body mask on z-faces - * @param[in] thin_xforce_lev x-component of forces on thin immersed bodies - * @param[in] thin_yforce_lev y-component of forces on thin immersed bodies - * @param[in] thin_zforce_lev z-component of forces on thin immersed bodies */ void erf_slow_rhs_inc (int level, int nrk, @@ -97,16 +91,13 @@ void erf_slow_rhs_inc (int level, int nrk, const MultiFab* p0, std::unique_ptr& mapfac_m, std::unique_ptr& mapfac_u, - std::unique_ptr& mapfac_v, - std::unique_ptr& xflux_imask_lev, - std::unique_ptr& yflux_imask_lev, - std::unique_ptr& zflux_imask_lev, - std::unique_ptr& thin_xforce_lev, - std::unique_ptr& thin_yforce_lev, - std::unique_ptr& thin_zforce_lev) + std::unique_ptr& mapfac_v) { BL_PROFILE_REGION("erf_slow_rhs_pre_inc()"); + const BCRec* bc_ptr_d = domain_bcs_type_d.data(); + const BCRec* bc_ptr_h = domain_bcs_type_h.data(); + DiffChoice dc = solverChoice.diffChoice; TurbChoice tc = solverChoice.turbChoice[level]; @@ -137,19 +128,12 @@ void erf_slow_rhs_inc (int level, int nrk, const bool use_most = (most != nullptr); - const BCRec* bc_ptr_d = domain_bcs_type_d.data(); - const BCRec* bc_ptr_h = domain_bcs_type_h.data(); - const Box& domain = geom.Domain(); const int domhi_z = domain.bigEnd(2); const int domlo_z = domain.smallEnd(2); const GpuArray dxInv = geom.InvCellSizeArray(); - const bool l_have_thin_xforce = (thin_xforce_lev != nullptr); - const bool l_have_thin_yforce = (thin_yforce_lev != nullptr); - const bool l_have_thin_zforce = (thin_zforce_lev != nullptr); - // ***************************************************************************** // Combine external forcing terms // ***************************************************************************** @@ -712,7 +696,7 @@ void erf_slow_rhs_inc (int level, int nrk, lo_z_face = mfi.validbox().smallEnd(2); hi_z_face = mfi.validbox().bigEnd(2)+1; } - AdvectionSrcForMom(tbx, tby, tbz, + AdvectionSrcForMom(bx, tbx, tby, tbz, rho_u_rhs, rho_v_rhs, rho_w_rhs, cell_data, u, v, w, rho_u, rho_v, omega_arr, @@ -756,6 +740,41 @@ void erf_slow_rhs_inc (int level, int nrk, } }); + // ***************************************************************************** + // Zero out source terms for x- and y- momenta if at walls or inflow + // We need to do this -- even though we call the boundary conditions later -- + // because the slow source is used to update the state in the fast interpolater. + // ***************************************************************************** + { + if ( (bx.smallEnd(0) == domain.smallEnd(0)) && + (bc_ptr_h[BCVars::xvel_bc].lo(0) == ERFBCType::ext_dir) ) { + Box lo_x_dom_face(bx); lo_x_dom_face.setBig(0,bx.smallEnd(0)); + ParallelFor(lo_x_dom_face, [=] AMREX_GPU_DEVICE (int i, int j, int k) { + rho_u_rhs(i,j,k) = 0.; + }); + } + if ( (bx.bigEnd(0) == domain.bigEnd(0)) && + (bc_ptr_h[BCVars::xvel_bc].hi(0) == ERFBCType::ext_dir) ) { + Box hi_x_dom_face(bx); hi_x_dom_face.setSmall(0,bx.bigEnd(0)+1); hi_x_dom_face.setBig(0,bx.bigEnd(0)+1); + ParallelFor(hi_x_dom_face, [=] AMREX_GPU_DEVICE (int i, int j, int k) { + rho_u_rhs(i,j,k) = 0.; + }); + } + if ( (bx.smallEnd(1) == domain.smallEnd(1)) && + (bc_ptr_h[BCVars::yvel_bc].lo(1) == ERFBCType::ext_dir) ) { + Box lo_y_dom_face(bx); lo_y_dom_face.setBig(1,bx.smallEnd(1)); + ParallelFor(lo_y_dom_face, [=] AMREX_GPU_DEVICE (int i, int j, int k) { + rho_v_rhs(i,j,k) = 0.; + }); + } + if ( (bx.bigEnd(1) == domain.bigEnd(1)) && + (bc_ptr_h[BCVars::yvel_bc].hi(1) == ERFBCType::ext_dir) ) { + Box hi_y_dom_face(bx); hi_y_dom_face.setSmall(1,bx.bigEnd(1)+1); hi_y_dom_face.setBig(1,bx.bigEnd(1)+1);; + ParallelFor(hi_y_dom_face, [=] AMREX_GPU_DEVICE (int i, int j, int k) { + rho_v_rhs(i,j,k) = 0.; + }); + } + amrex::Box b2d = tbz; b2d.setSmall(2,0); b2d.setBig(2,0); diff --git a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp index 7891900a0..e1d57d88e 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp @@ -1,7 +1,6 @@ #include #include #include -#include #include #include @@ -60,12 +59,6 @@ using namespace amrex; * @param[in] mapfac_v map factor at y-faces * @param[inout] fr_as_crse YAFluxRegister at level l at level l / l+1 interface * @param[inout] fr_as_fine YAFluxRegister at level l at level l-1 / l interface - * @param[in] xflux_imask_lev thin-body mask on x-faces - * @param[in] yflux_imask_lev thin-body mask on y-faces - * @param[in] zflux_imask_lev thin-body mask on z-faces - * @param[in] thin_xforce_lev x-component of forces on thin immersed bodies - * @param[in] thin_yforce_lev y-component of forces on thin immersed bodies - * @param[in] thin_zforce_lev z-component of forces on thin immersed bodies */ void erf_slow_rhs_pre (int level, int finest_level, @@ -108,13 +101,7 @@ void erf_slow_rhs_pre (int level, int finest_level, EBFArrayBoxFactory const& ebfact, #endif YAFluxRegister* fr_as_crse, - YAFluxRegister* fr_as_fine, - std::unique_ptr& xflux_imask_lev, - std::unique_ptr& yflux_imask_lev, - std::unique_ptr& zflux_imask_lev, - std::unique_ptr& thin_xforce_lev, - std::unique_ptr& thin_yforce_lev, - std::unique_ptr& thin_zforce_lev) + YAFluxRegister* fr_as_fine) { BL_PROFILE_REGION("erf_slow_rhs_pre()"); @@ -164,10 +151,6 @@ void erf_slow_rhs_pre (int level, int finest_level, const GpuArray dxInv = geom.InvCellSizeArray(); const Real* dx = geom.CellSize(); - const bool l_have_thin_xforce = (thin_xforce_lev != nullptr); - const bool l_have_thin_yforce = (thin_yforce_lev != nullptr); - const bool l_have_thin_zforce = (thin_zforce_lev != nullptr); - // ***************************************************************************** // Combine external forcing terms // ***************************************************************************** @@ -1011,76 +994,5 @@ void erf_slow_rhs_pre (int level, int finest_level, } // two-way coupling } // end profile } // mfi - - // Enforce thin immersed body interface condition, save forces - if (l_have_thin_xforce) { - MultiFab::Copy(*thin_xforce_lev, S_rhs[IntVars::xmom], 0, 0, 1, 0); - thin_xforce_lev->mult(-1., 0, 1, 0); - ApplyInvertedMask(*thin_xforce_lev, *xflux_imask_lev, 0); - MultiFab::Add(S_rhs[IntVars::xmom], *thin_xforce_lev, 0, 0, 1, 0); - - // TODO: Implement particles to better track and output these data -#ifndef AMREX_USE_GPU - if (nrk==2) { - for ( MFIter mfi(S_data[IntVars::cons],TileNoZ()); mfi.isValid(); ++mfi) - { - const Box& tbx = mfi.nodaltilebox(0); - const Array4 & fx = thin_xforce_lev->const_array(mfi); - const Array4 & mask = xflux_imask_lev->const_array(mfi); - ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - if (mask(i,j,k)==0) { - amrex::AllPrint() << "thin body fx"<mult(-1., 0, 1, 0); - ApplyInvertedMask(*thin_yforce_lev, *yflux_imask_lev, 0); - MultiFab::Add(S_rhs[IntVars::ymom], *thin_yforce_lev, 0, 0, 1, 0); - - // TODO: Implement particles to better track and output these data -#ifndef AMREX_USE_GPU - if (nrk==2) { - for ( MFIter mfi(S_data[IntVars::cons],TileNoZ()); mfi.isValid(); ++mfi) - { - const Box& tby = mfi.nodaltilebox(1); - const Array4 & fy = thin_yforce_lev->const_array(mfi); - const Array4 & mask = yflux_imask_lev->const_array(mfi); - ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - if (mask(i,j,k)==0) { - amrex::AllPrint() << "thin body fy"<mult(-1., 0, 1, 0); - ApplyInvertedMask(*thin_zforce_lev, *zflux_imask_lev, 0); - MultiFab::Add(S_rhs[IntVars::zmom], *thin_zforce_lev, 0, 0, 1, 0); - - // TODO: Implement particles to better track and output these data -#ifndef AMREX_USE_GPU - if (nrk==2) { - for ( MFIter mfi(S_data[IntVars::cons],TileNoZ()); mfi.isValid(); ++mfi) - { - const Box& tbz = mfi.nodaltilebox(2); - const Array4 & fz = thin_zforce_lev->const_array(mfi); - const Array4 & mask = zflux_imask_lev->const_array(mfi); - ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - if (mask(i,j,k)==0) { - amrex::AllPrint() << "thin body fz"<& S_old, Vector& F_slow, const Real time_for_fp, const Real slow_dt, +#ifdef ERF_USE_POISSON_SOLVE + const int nrk) +#else const int /*nrk*/) +#endif { BL_PROFILE("no_substep_fun"); int n_data = IntVars::NumTypes; diff --git a/Source/TimeIntegration/TI_slow_headers.H b/Source/TimeIntegration/TI_slow_headers.H index 997651b45..7d8621dbf 100644 --- a/Source/TimeIntegration/TI_slow_headers.H +++ b/Source/TimeIntegration/TI_slow_headers.H @@ -70,13 +70,7 @@ void erf_slow_rhs_pre (int level, int finest_level, int nrk, amrex::EBFArrayBoxFactory const& ebfact, #endif amrex::YAFluxRegister* fr_as_crse, - amrex::YAFluxRegister* fr_as_fine, - std::unique_ptr& xflux_imask_lev, - std::unique_ptr& yflux_imask_lev, - std::unique_ptr& zflux_imask_lev, - std::unique_ptr& thin_xforce_lev, - std::unique_ptr& thin_yforce_lev, - std::unique_ptr& thin_zforce_lev); + amrex::YAFluxRegister* fr_as_fine); /** * Function for computing the slow RHS for the evolution equations for the scalars other than density or potential temperature @@ -177,12 +171,6 @@ void erf_slow_rhs_inc (int level, int nrk, const amrex::MultiFab* p0, std::unique_ptr& mapfac_m, std::unique_ptr& mapfac_u, - std::unique_ptr& mapfac_v, - std::unique_ptr& xflux_imask_lev, - std::unique_ptr& yflux_imask_lev, - std::unique_ptr& zflux_imask_lev, - std::unique_ptr& thin_xforce_lev, - std::unique_ptr& thin_yforce_lev, - std::unique_ptr& thin_zforce_lev); + std::unique_ptr& mapfac_v); #endif #endif diff --git a/Source/TimeIntegration/TI_slow_rhs_fun.H b/Source/TimeIntegration/TI_slow_rhs_fun.H index 70182fb73..b07d4c237 100644 --- a/Source/TimeIntegration/TI_slow_rhs_fun.H +++ b/Source/TimeIntegration/TI_slow_rhs_fun.H @@ -37,6 +37,9 @@ } } + Real* dptr_u_geos = solverChoice.custom_geostrophic_profile ? d_u_geos[level].data(): nullptr; + Real* dptr_v_geos = solverChoice.custom_geostrophic_profile ? d_v_geos[level].data(): nullptr; + // Construct the source terms for the cell-centered (conserved) variables make_sources(level, nrk, slow_dt, S_data, S_prim, cc_src, #if defined(ERF_USE_RRTMGP) @@ -113,9 +116,8 @@ xmom_src, ymom_src, zmom_src, r0_new, fine_geom, solverChoice, mapfac_m[level], mapfac_u[level], mapfac_v[level], - dptr_rhotheta_src, dptr_rhoqt_src, - dptr_wbar_sub, d_rayleigh_ptrs_at_lev, - n_qstate); + dptr_u_geos, dptr_v_geos, dptr_wbar_sub, + d_rayleigh_ptrs_at_lev, n_qstate); erf_slow_rhs_pre(level, finest_level, nrk, slow_dt, S_rhs, S_data, S_prim, S_scratch, xvel_new, yvel_new, zvel_new, @@ -129,9 +131,11 @@ #ifdef ERF_USE_EB EBFactory(level), #endif - fr_as_crse, fr_as_fine, - xflux_imask[level], yflux_imask[level], zflux_imask[level], - thin_xforce[level], thin_yforce[level], thin_zforce[level]); + fr_as_crse, fr_as_fine); + + add_thin_body_sources(xmom_src, ymom_src, zmom_src, + xflux_imask[level], yflux_imask[level], zflux_imask[level], + thin_xforce[level], thin_yforce[level], thin_zforce[level]); // We define and evolve (rho theta)_0 in order to re-create p_0 in a way that is consistent // with our update of (rho theta) but does NOT maintain dp_0 / dz = -rho_0 g. This is why @@ -154,7 +158,6 @@ #endif for ( MFIter mfi(*p0,TilingIfNotGPU()); mfi.isValid(); ++mfi) { - const Array4 rt0_arr = rt0.array(mfi); const Array4 rt0_tmp_arr = rt0_new.array(mfi); @@ -213,9 +216,8 @@ xmom_src, ymom_src, zmom_src, r0, fine_geom, solverChoice, mapfac_m[level], mapfac_u[level], mapfac_v[level], - dptr_rhotheta_src, dptr_rhoqt_src, - dptr_wbar_sub, d_rayleigh_ptrs_at_lev, - n_qstate); + dptr_u_geos, dptr_v_geos, dptr_wbar_sub, + d_rayleigh_ptrs_at_lev, n_qstate); erf_slow_rhs_pre(level, finest_level, nrk, slow_dt, S_rhs, S_data, S_prim, S_scratch, xvel_new, yvel_new, zvel_new, @@ -229,9 +231,11 @@ #ifdef ERF_USE_EB EBFactory(level), #endif - fr_as_crse, fr_as_fine, - xflux_imask[level], yflux_imask[level], zflux_imask[level], - thin_xforce[level], thin_yforce[level], thin_zforce[level]); + fr_as_crse, fr_as_fine); + + add_thin_body_sources(xmom_src, ymom_src, zmom_src, + xflux_imask[level], yflux_imask[level], zflux_imask[level], + thin_xforce[level], thin_yforce[level], thin_zforce[level]); } #ifdef ERF_USE_NETCDF @@ -382,6 +386,9 @@ Real slow_dt = new_stage_time - old_step_time; + Real* dptr_u_geos = solverChoice.custom_geostrophic_profile ? d_u_geos[level].data(): nullptr; + Real* dptr_v_geos = solverChoice.custom_geostrophic_profile ? d_v_geos[level].data(): nullptr; + make_sources(level, nrk, slow_dt, S_data, S_prim, cc_src, #if defined(ERF_USE_RRTMGP) qheating_rates[level], @@ -397,9 +404,8 @@ xmom_src, ymom_src, zmom_src, r0, fine_geom, solverChoice, mapfac_m[level], mapfac_u[level], mapfac_v[level], - dptr_rhotheta_src, dptr_rhoqt_src, - dptr_wbar_sub, d_rayleigh_ptrs_at_lev, - n_qstate); + dptr_u_geos, dptr_v_geos, dptr_wbar_sub, + d_rayleigh_ptrs_at_lev, n_qstate); erf_slow_rhs_inc(level, nrk, slow_dt, S_rhs, S_old, S_data, S_prim, S_scratch, @@ -410,10 +416,11 @@ Tau32_lev[level].get(), SmnSmn, eddyDiffs, Hfx3, Diss, fine_geom, solverChoice, m_most, domain_bcs_type_d, domain_bcs_type, z_phys_nd[level], ax[level], ay[level], az[level], detJ_cc[level], p0, - mapfac_m[level], mapfac_u[level], mapfac_v[level], - xflux_imask[level], yflux_imask[level], zflux_imask[level], - thin_xforce[level], thin_yforce[level], thin_zforce[level]); + mapfac_m[level], mapfac_u[level], mapfac_v[level]); + add_thin_body_sources(xmom_src, ymom_src, zmom_src, + xflux_imask[level], yflux_imask[level], zflux_imask[level], + thin_xforce[level], thin_yforce[level], thin_zforce[level]); #ifdef ERF_USE_NETCDF // Populate RHS for relaxation zones if using real bcs @@ -428,16 +435,5 @@ } } #endif - -#if 0 - // HACK -- NO RELAXATION INSIDE FINE GRIDS - // Compute RHS for fine interior ghost - if (level > 0 && cf_width>0) { - 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], - domain_bcs_type, S_rhs, S_data); - } -#endif }; // end slow_rhs_fun_inc #endif