From af7f0f8392d96d8363df66e4573572b0b0639c5b Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Mon, 1 Apr 2024 13:34:11 -0700 Subject: [PATCH] clean up open bc stuff (#1546) --- Source/Advection/Advection.H | 9 +- Source/Advection/AdvectionSrcForMom.cpp | 89 ++++++++++++- Source/Advection/AdvectionSrcForState.cpp | 49 ++++++- Source/TimeIntegration/ERF_slow_rhs_inc.cpp | 100 +++++++++----- Source/TimeIntegration/ERF_slow_rhs_post.cpp | 41 +----- Source/TimeIntegration/ERF_slow_rhs_pre.cpp | 133 +++---------------- Source/TimeIntegration/TI_slow_headers.H | 1 - Source/TimeIntegration/TI_slow_rhs_fun.H | 7 +- 8 files changed, 230 insertions(+), 199 deletions(-) diff --git a/Source/Advection/Advection.H b/Source/Advection/Advection.H index aef42edf9..c5e9e4334 100644 --- a/Source/Advection/Advection.H +++ b/Source/Advection/Advection.H @@ -47,12 +47,15 @@ void AdvectionSrcForScalars (const amrex::Box& bx, const AdvType horiz_adv_type, const AdvType vert_adv_type, const amrex::Real horiz_upw_frac, const amrex::Real vert_upw_frac, const bool use_terrain, - const amrex::GpuArray, AMREX_SPACEDIM>& flx_arr); + const amrex::GpuArray, AMREX_SPACEDIM>& flx_arr, + const amrex::Box& domain, + const amrex::BCRec* bc_ptr_h); /** Compute advection tendencies for all components of momentum */ void AdvectionSrcForMom (const amrex::Box& bxx, const amrex::Box& bxy, const amrex::Box& bxz, const amrex::Array4< amrex::Real>& rho_u_rhs, const amrex::Array4< amrex::Real>& rho_v_rhs, const amrex::Array4< amrex::Real>& rho_w_rhs, + const amrex::Array4& rho, const amrex::Array4& u , const amrex::Array4& v, const amrex::Array4& w , const amrex::Array4& rho_u , const amrex::Array4& rho_v, @@ -64,7 +67,9 @@ void AdvectionSrcForMom (const amrex::Box& bxx, const amrex::Box& bxy, const amr const amrex::Array4& mf_v, const AdvType horiz_adv_type, const AdvType vert_adv_type, const amrex::Real horiz_upw_frac, const amrex::Real vert_upw_frac, - const bool use_terrain, const int lo_z_face, const int hi_z_face); + const bool use_terrain, const int lo_z_face, const int hi_z_face, + const amrex::Box& domain, + const amrex::BCRec* bc_ptr_h); /** Compute advection tendencies for all normal components of momentum */ void diff --git a/Source/Advection/AdvectionSrcForMom.cpp b/Source/Advection/AdvectionSrcForMom.cpp index 7d1cd2c5e..3875c40a9 100644 --- a/Source/Advection/AdvectionSrcForMom.cpp +++ b/Source/Advection/AdvectionSrcForMom.cpp @@ -1,3 +1,6 @@ +#include "AMReX_BCRec.H" + +#include #include #include @@ -30,13 +33,13 @@ using namespace amrex; * @param[in] horiz_adv_type sets the spatial order to be used for lateral derivatives * @param[in] vert_adv_type sets the spatial order to be used for vertical derivatives * @param[in] use_terrain if true, use the terrain-aware derivatives (with metric terms) - * @param[in] domhi_z maximum k value in the domain */ void AdvectionSrcForMom (const Box& bxx, const Box& bxy, const Box& bxz, const Array4< Real>& rho_u_rhs, const Array4< Real>& rho_v_rhs, const Array4< Real>& rho_w_rhs, + const Array4& cell_data, const Array4& u, const Array4& v, const Array4& w, @@ -54,7 +57,9 @@ AdvectionSrcForMom (const Box& bxx, const Box& bxy, const Box& bxz, const Real horiz_upw_frac, const Real vert_upw_frac, const bool use_terrain, - const int lo_z_face, const int hi_z_face) + const int lo_z_face, const int hi_z_face, + const Box& domain, + const BCRec* bc_ptr_h) { BL_PROFILE_VAR("AdvectionSrcForMom", AdvectionSrcForMom); @@ -313,5 +318,85 @@ AdvectionSrcForMom (const Box& bxx, const Box& bxy, const Box& bxz, } } } + + // Open bc will be imposed upon all vars (we only access cons here for simplicity) + const bool xlo_open = (bc_ptr_h[BCVars::cons_bc].lo(0) == ERFBCType::open); + const bool xhi_open = (bc_ptr_h[BCVars::cons_bc].hi(0) == ERFBCType::open); + const bool ylo_open = (bc_ptr_h[BCVars::cons_bc].lo(1) == ERFBCType::open); + const bool yhi_open = (bc_ptr_h[BCVars::cons_bc].hi(1) == ERFBCType::open); + + // Only advection operations in bndry normal direction with OPEN BC + const int domhi_z = domain.bigEnd(2); + Box tbx_xlo, tbx_xhi, tbx_ylo, tbx_yhi; + Box tby_xlo, tby_xhi, tby_ylo, tby_yhi; + Box tbz_xlo, tbz_xhi, tbz_ylo, tbz_yhi; + if (xlo_open) { + if (bxx.smallEnd(0) == domain.smallEnd(0)) { tbx_xlo = makeSlab(bxx,0,domain.smallEnd(0)); tbx_xlo.growLo(0,-1); } + if (bxy.smallEnd(0) == domain.smallEnd(0)) { tby_xlo = makeSlab(bxy,0,domain.smallEnd(0)); } + if (bxz.smallEnd(0) == domain.smallEnd(0)) { tbz_xlo = makeSlab(bxz,0,domain.smallEnd(0)); } + } + if (xhi_open) { + if (bxx.bigEnd(0) == domain.bigEnd(0)+1) { tbx_xhi = makeSlab(bxx,0,domain.bigEnd(0)+1); tbx_xhi.growHi(0,-1); } + if (bxy.bigEnd(0) == domain.bigEnd(0)) { tby_xhi = makeSlab(bxy,0,domain.bigEnd(0) ); } + if (bxz.bigEnd(0) == domain.bigEnd(0)) { tbz_xhi = makeSlab(bxz,0,domain.bigEnd(0) ); } + } + if (ylo_open) { + if (bxx.smallEnd(1) == domain.smallEnd(1)) { tbx_ylo = makeSlab(bxx,1,domain.smallEnd(1)); } + if (bxy.smallEnd(1) == domain.smallEnd(1)) { tby_ylo = makeSlab(bxy,1,domain.smallEnd(1)); tby_ylo.growLo(1,-1); } + if (bxz.smallEnd(1) == domain.smallEnd(1)) { tbz_ylo = makeSlab(bxz,1,domain.smallEnd(1)); } + } + if (yhi_open) { + if (bxx.bigEnd(1) == domain.bigEnd(1)) { tbx_yhi = makeSlab(bxx,1,domain.bigEnd(1) ); } + if (bxy.bigEnd(1) == domain.bigEnd(1)+1) { tby_yhi = makeSlab(bxy,1,domain.bigEnd(1)+1); tby_yhi.growHi(1,-1); } + if (bxz.bigEnd(1) == domain.bigEnd(1)) { tbz_yhi = makeSlab(bxz,1,domain.bigEnd(1) ); } + } + + // Special advection operator for open BC (bndry normal/tangent operations) + if (xlo_open) { + bool do_lo = true; + AdvectionSrcForOpenBC_Normal(tbx_xlo, 0, rho_u_rhs, u, cell_data, cellSizeInv, do_lo); + AdvectionSrcForOpenBC_Tangent_Ymom(tby_xlo, 0, rho_v_rhs, v, + rho_u, rho_v, Omega, + z_nd, detJ, cellSizeInv, + use_terrain, do_lo); + AdvectionSrcForOpenBC_Tangent_Zmom(tbz_xlo, 0, rho_w_rhs, w, + rho_u, rho_v, Omega, + z_nd, detJ, cellSizeInv, + use_terrain, domhi_z, do_lo); + } + if (xhi_open) { + AdvectionSrcForOpenBC_Normal(tbx_xhi, 0, rho_u_rhs, u, cell_data, cellSizeInv); + AdvectionSrcForOpenBC_Tangent_Ymom(tby_xhi, 0, rho_v_rhs, v, + rho_u, rho_v, Omega, + z_nd, detJ, cellSizeInv, + use_terrain); + AdvectionSrcForOpenBC_Tangent_Zmom(tbz_xhi, 0, rho_w_rhs, w, + rho_u, rho_v, Omega, + z_nd, detJ, cellSizeInv, + use_terrain, domhi_z); + } + if (ylo_open) { + bool do_lo = true; + AdvectionSrcForOpenBC_Tangent_Xmom(tbx_ylo, 1, rho_u_rhs, u, + rho_u, rho_v, Omega, + z_nd, detJ, cellSizeInv, + use_terrain, do_lo); + AdvectionSrcForOpenBC_Normal(tby_ylo, 1, rho_v_rhs, v, cell_data, cellSizeInv, do_lo); + AdvectionSrcForOpenBC_Tangent_Zmom(tbz_ylo, 1, rho_w_rhs, w, + rho_u, rho_v, Omega, + z_nd, detJ, cellSizeInv, + use_terrain, domhi_z, do_lo); + } + if (yhi_open) { + AdvectionSrcForOpenBC_Tangent_Xmom(tbx_yhi, 1, rho_u_rhs, u, + rho_u, rho_v, Omega, + z_nd, detJ, cellSizeInv, + use_terrain); + AdvectionSrcForOpenBC_Normal(tby_yhi, 1, rho_v_rhs, v, cell_data, cellSizeInv); + AdvectionSrcForOpenBC_Tangent_Zmom(tbz_yhi, 1, rho_w_rhs, w, + rho_u, rho_v, Omega, + z_nd, detJ, cellSizeInv, + use_terrain, domhi_z); + } } diff --git a/Source/Advection/AdvectionSrcForState.cpp b/Source/Advection/AdvectionSrcForState.cpp index 732052e24..263927312 100644 --- a/Source/Advection/AdvectionSrcForState.cpp +++ b/Source/Advection/AdvectionSrcForState.cpp @@ -168,7 +168,9 @@ AdvectionSrcForScalars (const Box& bx, const int icomp, const int ncomp, const Real horiz_upw_frac, const Real vert_upw_frac, const bool use_terrain, - const GpuArray, AMREX_SPACEDIM>& flx_arr) + const GpuArray, AMREX_SPACEDIM>& flx_arr, + const Box& domain, + const BCRec* bc_ptr_h) { BL_PROFILE_VAR("AdvectionSrcForScalars", AdvectionSrcForScalars); auto dxInv = cellSizeInv[0], dyInv = cellSizeInv[1], dzInv = cellSizeInv[2]; @@ -177,6 +179,27 @@ AdvectionSrcForScalars (const Box& bx, const int icomp, const int ncomp, const Box ybx = surroundingNodes(bx,1); const Box zbx = surroundingNodes(bx,2); + // Open bc will be imposed upon all vars (we only access cons here for simplicity) + const bool xlo_open = (bc_ptr_h[BCVars::cons_bc].lo(0) == ERFBCType::open); + const bool xhi_open = (bc_ptr_h[BCVars::cons_bc].hi(0) == ERFBCType::open); + const bool ylo_open = (bc_ptr_h[BCVars::cons_bc].lo(1) == ERFBCType::open); + const bool yhi_open = (bc_ptr_h[BCVars::cons_bc].hi(1) == ERFBCType::open); + + // Only advection operations in bndry normal direction with OPEN BC + Box bx_xlo, bx_xhi, bx_ylo, bx_yhi; + if (xlo_open) { + if ( bx.smallEnd(0) == domain.smallEnd(0)) { bx_xlo = makeSlab( bx,0,domain.smallEnd(0));} + } + if (xhi_open) { + if ( bx.bigEnd(0) == domain.bigEnd(0)) { bx_xhi = makeSlab( bx,0,domain.bigEnd(0) );} + } + if (ylo_open) { + if ( bx.smallEnd(1) == domain.smallEnd(1)) { bx_ylo = makeSlab( bx,1,domain.smallEnd(1));} + } + if (yhi_open) { + if ( bx.bigEnd(1) == domain.bigEnd(1)) { bx_yhi = makeSlab( bx,1,domain.bigEnd(1) );} + } + // Inline with 2nd order for efficiency // NOTE: we don't need to weight avg_xmom, avg_ymom, avg_zmom with terrain metrics // (or with EB area fractions) @@ -279,4 +302,28 @@ AdvectionSrcForScalars (const Box& bx, const int icomp, const int ncomp, ( (flx_arr[1])(i,j+1,k,cons_index) - (flx_arr[1])(i,j ,k,cons_index) ) * dyInv + ( (flx_arr[2])(i,j,k+1,cons_index) - (flx_arr[2])(i,j,k ,cons_index) ) * dzInv ); }); + + // Special advection operator for open BC (bndry tangent operations) + if (xlo_open) { + bool do_lo = true; + AdvectionSrcForOpenBC_Tangent_Cons(bx_xlo, 0, icomp, ncomp, advectionSrc, cell_prim, + avg_xmom, avg_ymom, avg_zmom, + detJ, cellSizeInv, use_terrain, do_lo); + } + if (xhi_open) { + AdvectionSrcForOpenBC_Tangent_Cons(bx_xhi, 0, icomp, ncomp, advectionSrc, cell_prim, + avg_xmom, avg_ymom, avg_zmom, + detJ, cellSizeInv, use_terrain); + } + if (ylo_open) { + bool do_lo = true; + AdvectionSrcForOpenBC_Tangent_Cons(bx_ylo, 1, icomp, ncomp, advectionSrc, cell_prim, + avg_xmom, avg_ymom, avg_zmom, + detJ, cellSizeInv, use_terrain, do_lo); + } + if (yhi_open) { + AdvectionSrcForOpenBC_Tangent_Cons(bx_yhi, 1, icomp, ncomp, advectionSrc, cell_prim, + avg_xmom, avg_ymom, avg_zmom, + detJ, cellSizeInv, use_terrain); + } } diff --git a/Source/TimeIntegration/ERF_slow_rhs_inc.cpp b/Source/TimeIntegration/ERF_slow_rhs_inc.cpp index a0a52d7e1..2bc27505e 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_inc.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_inc.cpp @@ -50,7 +50,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] dptr_rhotheta_src custom temperature source term * @param[in] dptr_u_geos custom geostrophic wind profile * @param[in] dptr_v_geos custom geostrophic wind profile * @param[in] d_rayleigh_ptrs_at_lev Vector of {strength of Rayleigh damping, reference value for xvel/yvel/zvel/theta} used to define Rayleigh damping @@ -86,7 +85,6 @@ void erf_slow_rhs_inc (int level, int nrk, std::unique_ptr& mapfac_m, std::unique_ptr& mapfac_u, std::unique_ptr& mapfac_v, - const Real* dptr_rhotheta_src, const Real* dptr_u_geos, const Real* dptr_v_geos, const Real* dptr_wbar_sub, @@ -138,15 +136,15 @@ void erf_slow_rhs_inc (int level, int nrk, Real* vbar = d_rayleigh_ptrs_at_lev[Rayleigh::vbar]; Real* wbar = d_rayleigh_ptrs_at_lev[Rayleigh::wbar]; - // ************************************************************************* + // ***************************************************************************** // Combine external forcing terms - // ************************************************************************* + // ***************************************************************************** const Array grav{0.0, 0.0, -solverChoice.gravity}; const GpuArray grav_gpu{grav[0], grav[1], grav[2]}; - // ************************************************************************* + // ***************************************************************************** // Pre-computed quantities - // ************************************************************************* + // ***************************************************************************** int nvars = S_data[IntVars::cons].nComp(); const BoxArray& ba = S_data[IntVars::cons].boxArray(); const DistributionMapping& dm = S_data[IntVars::cons].DistributionMap(); @@ -169,7 +167,7 @@ void erf_slow_rhs_inc (int level, int nrk, : 2.0 * dc.dynamicViscosity; #ifdef _OPENMP -#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#pragma omp parallel if (Gpu::notInLaunchRegion()) #endif for ( MFIter mfi(S_data[IntVars::cons],TileNoZ()); mfi.isValid(); ++mfi) { @@ -218,7 +216,8 @@ void erf_slow_rhs_inc (int level, int nrk, 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 + + // 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)); @@ -249,9 +248,9 @@ void erf_slow_rhs_inc (int level, int nrk, Array4 s21 = S21.array(); Array4 s31 = S31.array(); Array4 s32 = S32.array(); Array4 tau21 = Tau21->array(mfi); Array4 tau31 = Tau31->array(mfi); Array4 tau32 = Tau32->array(mfi); - //----------------------------------------- + // ***************************************************************************** // Expansion rate compute terrain - //----------------------------------------- + // ***************************************************************************** { BL_PROFILE("slow_rhs_making_er_T"); // First create Omega using velocity (not momentum) @@ -283,9 +282,9 @@ void erf_slow_rhs_inc (int level, int nrk, }); } // end profile - //----------------------------------------- + // ***************************************************************************** // Strain tensor compute terrain - //----------------------------------------- + // ***************************************************************************** { BL_PROFILE("slow_rhs_making_strain_T"); ComputeStrain_T(bxcc, tbxxy, tbxxz, tbxyz, @@ -308,9 +307,21 @@ void erf_slow_rhs_inc (int level, int nrk, }); } - //----------------------------------------- +#ifdef ERF_EXPLICIT_MOST_STRESS + // We've updated the strains at all locations including the + // surface. This is required to get the correct strain-rate + // magnitude. Now, update the stress everywhere but the surface + // to retain the values set by MOST. + if (use_most) { + // Don't overwrite modeled total stress value at boundary + tbxxz.setSmall(2,1); + tbxyz.setSmall(2,1); + } +#endif + + // ***************************************************************************** // Stress tensor compute terrain - //----------------------------------------- + // ***************************************************************************** { BL_PROFILE("slow_rhs_making_stress_T"); @@ -369,9 +380,9 @@ void erf_slow_rhs_inc (int level, int nrk, } else { - //----------------------------------------- + // ***************************************************************************** // Expansion rate compute no terrain - //----------------------------------------- + // ***************************************************************************** { BL_PROFILE("slow_rhs_making_er_N"); ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { @@ -383,9 +394,9 @@ void erf_slow_rhs_inc (int level, int nrk, } // end profile - //----------------------------------------- + // ***************************************************************************** // Strain tensor compute no terrain - //----------------------------------------- + // ***************************************************************************** { BL_PROFILE("slow_rhs_making_strain_N"); ComputeStrain_N(bxcc, tbxxy, tbxxz, tbxyz, @@ -406,9 +417,9 @@ void erf_slow_rhs_inc (int level, int nrk, }); } - //----------------------------------------- + // ***************************************************************************** // Stress tensor compute no terrain - //----------------------------------------- + // ***************************************************************************** { BL_PROFILE("slow_rhs_making_stress_N"); @@ -466,8 +477,15 @@ void erf_slow_rhs_inc (int level, int nrk, // ************************************************************************* // Define updates and fluxes in the current RK stage // ************************************************************************* + + // Open bc will be imposed upon all vars (we only access cons here for simplicity) + const bool xlo_open = (bc_ptr_h[BCVars::cons_bc].lo(0) == ERFBCType::open); + const bool xhi_open = (bc_ptr_h[BCVars::cons_bc].hi(0) == ERFBCType::open); + const bool ylo_open = (bc_ptr_h[BCVars::cons_bc].lo(1) == ERFBCType::open); + const bool yhi_open = (bc_ptr_h[BCVars::cons_bc].hi(1) == ERFBCType::open); + #ifdef _OPENMP -#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#pragma omp parallel if (Gpu::notInLaunchRegion()) #endif std::array flux; @@ -478,6 +496,12 @@ void erf_slow_rhs_inc (int level, int nrk, Box tby = mfi.nodaltilebox(1); Box tbz = mfi.nodaltilebox(2); + // If we are imposing open bc's then don't add rhs terms at the boundary locations + if ( xlo_open && (tbx.smallEnd(0) == domain.smallEnd(0)) ) {tbx.growLo(0,-1);} + if ( xhi_open && (tbx.bigEnd(0) == domain.bigEnd(0)+1) ) {tbx.growHi(0,-1);} + if ( ylo_open && (tby.smallEnd(1) == domain.smallEnd(1)) ) {tby.growLo(1,-1);} + if ( yhi_open && (tby.bigEnd(1) == domain.bigEnd(1)+1) ) {tby.growHi(1,-1);} + // We don't compute a source term for z-momentum on the bottom or top boundary tbz.growLo(2,-1); tbz.growHi(2,-1); @@ -546,8 +570,10 @@ void erf_slow_rhs_inc (int level, int nrk, } const GpuArray, AMREX_SPACEDIM> flx_arr{{AMREX_D_DECL(flux[0].array(), flux[1].array(), flux[2].array())}}; + + // ***************************************************************************** // Contravariant flux field - //----------------------------------------- + // ***************************************************************************** { BL_PROFILE("slow_rhs_making_omega"); Box gbxo = surroundingNodes(bx,2); gbxo.grow(IntVect(1,1,0)); @@ -558,9 +584,9 @@ void erf_slow_rhs_inc (int level, int nrk, } // end profile - //----------------------------------------- + // ***************************************************************************** // Diffusive terms (pre-computed above) - //----------------------------------------- + // ***************************************************************************** // Expansion Array4 er_arr; if (expr) { @@ -610,7 +636,7 @@ void erf_slow_rhs_inc (int level, int nrk, cell_prim, cell_rhs, detJ_arr, dxInv, mf_m, l_horiz_adv_type, l_vert_adv_type, l_horiz_upw_frac, l_vert_upw_frac, - l_use_terrain, flx_arr); + l_use_terrain, flx_arr, domain, bc_ptr_h); if (l_use_diff) { Array4 diffflux_x = dflux_x->array(mfi); @@ -657,9 +683,9 @@ void erf_slow_rhs_inc (int level, int nrk, }); } - // ********************************************************************* + // ***************************************************************************** // Define updates in the RHS of {x, y, z}-momentum equations - // ********************************************************************* + // ***************************************************************************** int lo_z_face; int hi_z_face; if (level == 0) { @@ -670,12 +696,14 @@ void erf_slow_rhs_inc (int level, int nrk, hi_z_face = mfi.validbox().bigEnd(2)+1; } AdvectionSrcForMom(tbx, tby, tbz, - rho_u_rhs, rho_v_rhs, rho_w_rhs, u, v, w, - rho_u , rho_v , omega_arr, + rho_u_rhs, rho_v_rhs, rho_w_rhs, + cell_data, u, v, w, + rho_u, rho_v, omega_arr, z_nd, detJ_arr, dxInv, mf_m, mf_u, mf_v, l_horiz_adv_type, l_vert_adv_type, l_horiz_upw_frac, l_vert_upw_frac, - l_use_terrain, lo_z_face, hi_z_face); + l_use_terrain, lo_z_face, hi_z_face, + domain, bc_ptr_h); if (l_use_diff) { DiffusionSrcForMom_N(tbx, tby, tbz, @@ -699,9 +727,9 @@ void erf_slow_rhs_inc (int level, int nrk, { BL_PROFILE("slow_rhs_inc_xmom"); - // ****************************************************************** + // ***************************************************************************** // NON-TERRAIN VERSION - // ****************************************************************** + // ***************************************************************************** ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // x-momentum equation @@ -742,9 +770,9 @@ void erf_slow_rhs_inc (int level, int nrk, { BL_PROFILE("slow_rhs_inc_ymom"); - // ****************************************************************** + // ***************************************************************************** // NON-TERRAIN VERSION - // ****************************************************************** + // ***************************************************************************** ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // y-momentum equation @@ -796,9 +824,9 @@ void erf_slow_rhs_inc (int level, int nrk, { BL_PROFILE("slow_rhs_pre_zmom"); - // ****************************************************************** + // ***************************************************************************** // NON-TERRAIN VERSION - // ****************************************************************** + // ***************************************************************************** ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // z-momentum equation diff --git a/Source/TimeIntegration/ERF_slow_rhs_post.cpp b/Source/TimeIntegration/ERF_slow_rhs_post.cpp index bb16be0a3..46406304e 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_post.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_post.cpp @@ -192,21 +192,6 @@ void erf_slow_rhs_post (int level, int finest_level, Box tbx = mfi.tilebox(); - // Only advection operations in bndry normal direction with OPEN BC - Box tbx_xlo, tbx_xhi, tbx_ylo, tbx_yhi; - if (xlo_open) { - if (tbx.smallEnd(0) == domain.smallEnd(0)) { tbx_xlo = makeSlab(tbx,0,domain.smallEnd(0)); } - } - if (xhi_open) { - if (tbx.bigEnd(0) == domain.bigEnd(0)) { tbx_xhi = makeSlab(tbx,0,domain.bigEnd(0)); } - } - if (ylo_open) { - if (tbx.smallEnd(1) == domain.smallEnd(1)) { tbx_ylo = makeSlab(tbx,1,domain.smallEnd(1)); } - } - if (yhi_open) { - if (tbx.bigEnd(1) == domain.bigEnd(1)) { tbx_yhi = makeSlab(tbx,1,domain.bigEnd(1)); } - } - // ************************************************************************* // Define flux arrays for use in advection // ************************************************************************* @@ -359,31 +344,7 @@ void erf_slow_rhs_post (int level, int finest_level, dxInv, mf_m, horiz_adv_type, vert_adv_type, horiz_upw_frac, vert_upw_frac, - l_use_terrain, flx_arr); - - // Special advection operator for open BC (bndry normal operations) - if (xlo_open) { - bool do_lo = true; - AdvectionSrcForOpenBC_Tangent_Cons(tbx_xlo, 0, ivar, num_comp, cell_rhs, cur_prim, - avg_xmom, avg_ymom, avg_zmom, - detJ_arr, dxInv, l_use_terrain, do_lo); - } - if (xhi_open) { - AdvectionSrcForOpenBC_Tangent_Cons(tbx_xhi, 0, ivar, num_comp, cell_rhs, cur_prim, - avg_xmom, avg_ymom, avg_zmom, - detJ_arr, dxInv, l_use_terrain); - } - if (ylo_open) { - bool do_lo = true; - AdvectionSrcForOpenBC_Tangent_Cons(tbx_ylo, 1, ivar, num_comp, cell_rhs, cur_prim, - avg_xmom, avg_ymom, avg_zmom, - detJ_arr, dxInv, l_use_terrain, do_lo); - } - if (yhi_open) { - AdvectionSrcForOpenBC_Tangent_Cons(tbx_yhi, 1, ivar, num_comp, cell_rhs, cur_prim, - avg_xmom, avg_ymom, avg_zmom, - detJ_arr, dxInv, l_use_terrain); - } + l_use_terrain, flx_arr, domain, bc_ptr_h); if (l_use_diff) { diff --git a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp index 5c476e523..c86835ee0 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp @@ -138,12 +138,6 @@ void erf_slow_rhs_pre (int level, int finest_level, const bool use_moisture = (solverChoice.moisture_type != MoistureType::None); const bool use_most = (most != nullptr); - // Open bc will be imposed upon all vars (we only access cons here for simplicity) - const bool xlo_open = (bc_ptr_h[BCVars::cons_bc].lo(0) == ERFBCType::open); - const bool xhi_open = (bc_ptr_h[BCVars::cons_bc].hi(0) == ERFBCType::open); - const bool ylo_open = (bc_ptr_h[BCVars::cons_bc].lo(1) == ERFBCType::open); - const bool yhi_open = (bc_ptr_h[BCVars::cons_bc].hi(1) == ERFBCType::open); - const Box& domain = geom.Domain(); const int domhi_z = domain.bigEnd(2); const int domlo_z = domain.smallEnd(2); @@ -555,6 +549,13 @@ void erf_slow_rhs_pre (int level, int finest_level, // ***************************************************************************** // Define updates and fluxes in the current RK stage // ***************************************************************************** + + // Open bc will be imposed upon all vars (we only access cons here for simplicity) + const bool xlo_open = (bc_ptr_h[BCVars::cons_bc].lo(0) == ERFBCType::open); + const bool xhi_open = (bc_ptr_h[BCVars::cons_bc].hi(0) == ERFBCType::open); + const bool ylo_open = (bc_ptr_h[BCVars::cons_bc].lo(1) == ERFBCType::open); + const bool yhi_open = (bc_ptr_h[BCVars::cons_bc].hi(1) == ERFBCType::open); + #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif @@ -568,40 +569,16 @@ void erf_slow_rhs_pre (int level, int finest_level, Box tby = mfi.nodaltilebox(1); Box tbz = mfi.nodaltilebox(2); + // If we are imposing open bc's then don't add rhs terms at the boundary locations + if ( xlo_open && (tbx.smallEnd(0) == domain.smallEnd(0)) ) {tbx.growLo(0,-1);} + if ( xhi_open && (tbx.bigEnd(0) == domain.bigEnd(0)+1) ) {tbx.growHi(0,-1);} + if ( ylo_open && (tby.smallEnd(1) == domain.smallEnd(1)) ) {tby.growLo(1,-1);} + if ( yhi_open && (tby.bigEnd(1) == domain.bigEnd(1)+1) ) {tby.growHi(1,-1);} + // We don't compute a source term for z-momentum on the bottom or top boundary tbz.growLo(2,-1); tbz.growHi(2,-1); - // Only advection operations in bndry normal direction with OPEN BC - Box bx_xlo, bx_xhi, bx_ylo, bx_yhi; - Box tbx_xlo, tbx_xhi, tbx_ylo, tbx_yhi; - Box tby_xlo, tby_xhi, tby_ylo, tby_yhi; - Box tbz_xlo, tbz_xhi, tbz_ylo, tbz_yhi; - if (xlo_open) { - if ( bx.smallEnd(0) == domain.smallEnd(0)) { bx_xlo = makeSlab( bx,0,domain.smallEnd(0)); } - if (tbx.smallEnd(0) == domain.smallEnd(0)) { tbx_xlo = makeSlab(tbx,0,domain.smallEnd(0)); tbx.growLo(0,-1); } - if (tby.smallEnd(0) == domain.smallEnd(0)) { tby_xlo = makeSlab(tby,0,domain.smallEnd(0)); } - if (tbz.smallEnd(0) == domain.smallEnd(0)) { tbz_xlo = makeSlab(tbz,0,domain.smallEnd(0)); } - } - if (xhi_open) { - if ( bx.bigEnd(0) == domain.bigEnd(0)) { bx_xhi = makeSlab( bx,0,domain.bigEnd(0) ); } - if (tbx.bigEnd(0) == domain.bigEnd(0)+1) { tbx_xhi = makeSlab(tbx,0,domain.bigEnd(0)+1); tbx.growHi(0,-1); } - if (tby.bigEnd(0) == domain.bigEnd(0)) { tby_xhi = makeSlab(tby,0,domain.bigEnd(0) ); } - if (tbz.bigEnd(0) == domain.bigEnd(0)) { tbz_xhi = makeSlab(tbz,0,domain.bigEnd(0) ); } - } - if (ylo_open) { - if ( bx.smallEnd(1) == domain.smallEnd(1)) { bx_ylo = makeSlab( bx,1,domain.smallEnd(1)); } - if (tbx.smallEnd(1) == domain.smallEnd(1)) { tbx_ylo = makeSlab(tbx,1,domain.smallEnd(1)); } - if (tby.smallEnd(1) == domain.smallEnd(1)) { tby_ylo = makeSlab(tby,1,domain.smallEnd(1)); tby.growLo(1,-1); } - if (tbz.smallEnd(1) == domain.smallEnd(1)) { tbz_ylo = makeSlab(tbz,1,domain.smallEnd(1)); } - } - if (yhi_open) { - if ( bx.bigEnd(1) == domain.bigEnd(1)) { bx_yhi = makeSlab( bx,1,domain.bigEnd(1) ); } - if (tbx.bigEnd(1) == domain.bigEnd(1)) { tbx_yhi = makeSlab(tbx,1,domain.bigEnd(1) ); } - if (tby.bigEnd(1) == domain.bigEnd(1)+1) { tby_yhi = makeSlab(tby,1,domain.bigEnd(1)+1); tby.growHi(1,-1); } - if (tbz.bigEnd(1) == domain.bigEnd(1)) { tbz_yhi = makeSlab(tbz,1,domain.bigEnd(1) ); } - } - const Array4 & cell_data = S_data[IntVars::cons].array(mfi); const Array4 & cell_prim = S_prim.array(mfi); const Array4 & cell_rhs = S_rhs[IntVars::cons].array(mfi); @@ -793,33 +770,7 @@ void erf_slow_rhs_pre (int level, int finest_level, dxInv, mf_m, l_horiz_adv_type, l_vert_adv_type, l_horiz_upw_frac, l_vert_upw_frac, - l_use_terrain, flx_arr); - - - // Special advection operator for open BC (bndry tangent operations) - if (xlo_open) { - bool do_lo = true; - AdvectionSrcForOpenBC_Tangent_Cons(bx_xlo, 0, RhoTheta_comp, 1, cell_rhs, cell_prim, - avg_xmom, avg_ymom, avg_zmom, - detJ_arr, dxInv, l_use_terrain, do_lo); - } - if (xhi_open) { - AdvectionSrcForOpenBC_Tangent_Cons(bx_xhi, 0, RhoTheta_comp, 1, cell_rhs, cell_prim, - avg_xmom, avg_ymom, avg_zmom, - detJ_arr, dxInv, l_use_terrain); - } - if (ylo_open) { - bool do_lo = true; - AdvectionSrcForOpenBC_Tangent_Cons(bx_ylo, 1, RhoTheta_comp, 1, cell_rhs, cell_prim, - avg_xmom, avg_ymom, avg_zmom, - detJ_arr, dxInv, l_use_terrain, do_lo); - } - if (yhi_open) { - AdvectionSrcForOpenBC_Tangent_Cons(bx_yhi, 1, RhoTheta_comp, 1, cell_rhs, cell_prim, - avg_xmom, avg_ymom, avg_zmom, - detJ_arr, dxInv, l_use_terrain); - } - + l_use_terrain, flx_arr, domain, bc_ptr_h); if (l_use_diff) { Array4 diffflux_x = dflux_x->array(mfi); @@ -882,60 +833,14 @@ void erf_slow_rhs_pre (int level, int finest_level, hi_z_face = mfi.validbox().bigEnd(2)+1; } AdvectionSrcForMom(tbx, tby, tbz, - rho_u_rhs, rho_v_rhs, rho_w_rhs, u, v, w, - rho_u , rho_v , omega_arr, + rho_u_rhs, rho_v_rhs, rho_w_rhs, + cell_data, u, v, w, + rho_u, rho_v, omega_arr, z_nd, detJ_arr, dxInv, mf_m, mf_u, mf_v, l_horiz_adv_type, l_vert_adv_type, l_horiz_upw_frac, l_vert_upw_frac, - l_use_terrain, lo_z_face, hi_z_face); - - // Special advection operator for open BC (bndry normal/tangent operations) - if (xlo_open) { - bool do_lo = true; - AdvectionSrcForOpenBC_Normal(tbx_xlo, 0, rho_u_rhs, u, cell_data, dxInv, do_lo); - AdvectionSrcForOpenBC_Tangent_Ymom(tby_xlo, 0, rho_v_rhs, v, - rho_u, rho_v, omega_arr, - z_nd, detJ_arr, dxInv, - l_use_terrain, do_lo); - AdvectionSrcForOpenBC_Tangent_Zmom(tbz_xlo, 0, rho_w_rhs, w, - rho_u, rho_v, omega_arr, - z_nd, detJ_arr, dxInv, - l_use_terrain, domhi_z, do_lo); - } - if (xhi_open) { - AdvectionSrcForOpenBC_Normal(tbx_xhi, 0, rho_u_rhs, u, cell_data, dxInv); - AdvectionSrcForOpenBC_Tangent_Ymom(tby_xhi, 0, rho_v_rhs, v, - rho_u, rho_v, omega_arr, - z_nd, detJ_arr, dxInv, - l_use_terrain); - AdvectionSrcForOpenBC_Tangent_Zmom(tbz_xhi, 0, rho_w_rhs, w, - rho_u, rho_v, omega_arr, - z_nd, detJ_arr, dxInv, - l_use_terrain, domhi_z); - } - if (ylo_open) { - bool do_lo = true; - AdvectionSrcForOpenBC_Tangent_Xmom(tbx_ylo, 1, rho_u_rhs, u, - rho_u, rho_v, omega_arr, - z_nd, detJ_arr, dxInv, - l_use_terrain, do_lo); - AdvectionSrcForOpenBC_Normal(tby_ylo, 1, rho_v_rhs, v, cell_data, dxInv, do_lo); - AdvectionSrcForOpenBC_Tangent_Zmom(tbz_ylo, 1, rho_w_rhs, w, - rho_u, rho_v, omega_arr, - z_nd, detJ_arr, dxInv, - l_use_terrain, domhi_z, do_lo); - } - if (yhi_open) { - AdvectionSrcForOpenBC_Tangent_Xmom(tbx_yhi, 1, rho_u_rhs, u, - rho_u, rho_v, omega_arr, - z_nd, detJ_arr, dxInv, - l_use_terrain); - AdvectionSrcForOpenBC_Normal(tby_yhi, 1, rho_v_rhs, v, cell_data, dxInv); - AdvectionSrcForOpenBC_Tangent_Zmom(tbz_yhi, 1, rho_w_rhs, w, - rho_u, rho_v, omega_arr, - z_nd, detJ_arr, dxInv, - l_use_terrain, domhi_z); - } + l_use_terrain, lo_z_face, hi_z_face, + domain, bc_ptr_h); if (l_use_diff) { // Note: tau** were calculated with calls to diff --git a/Source/TimeIntegration/TI_slow_headers.H b/Source/TimeIntegration/TI_slow_headers.H index a7ade7ec0..031871904 100644 --- a/Source/TimeIntegration/TI_slow_headers.H +++ b/Source/TimeIntegration/TI_slow_headers.H @@ -165,7 +165,6 @@ void erf_slow_rhs_inc (int level, int nrk, std::unique_ptr& mapfac_m, std::unique_ptr& mapfac_u, std::unique_ptr& mapfac_v, - const amrex::Real* dptr_rhotheta_src, const amrex::Real* dptr_u_geos, const amrex::Real* dptr_v_geos, const amrex::Real* dptr_wbar_sub, diff --git a/Source/TimeIntegration/TI_slow_rhs_fun.H b/Source/TimeIntegration/TI_slow_rhs_fun.H index a4322d7be..3f10fbced 100644 --- a/Source/TimeIntegration/TI_slow_rhs_fun.H +++ b/Source/TimeIntegration/TI_slow_rhs_fun.H @@ -379,9 +379,10 @@ erf_slow_rhs_inc(level, nrk, slow_dt, S_rhs, S_old, S_data, S_prim, S_scratch, xvel_new, yvel_new, zvel_new, - z_t_rk[level], Omega, source, buoyancy, Tau11, Tau22, Tau33, Tau12, - Tau13, Tau21, Tau23, Tau31, Tau32, SmnSmn, eddyDiffs, - Hfx3, Diss, + z_t_rk[level], Omega, source, buoyancy, + Tau11_lev[level].get(), Tau22_lev[level].get(), Tau33_lev[level].get(), Tau12_lev[level].get(), + Tau13_lev[level].get(), Tau21_lev[level].get(), Tau23_lev[level].get(), Tau31_lev[level].get(), + Tau32_lev[level].get(), SmnSmn, eddyDiffs, Hfx3, Diss, fine_geom, solverChoice, m_most, domain_bcs_type_d, domain_bcs_type, z_phys_nd[level], detJ_cc[level], p0, mapfac_m[level], mapfac_u[level], mapfac_v[level],