diff --git a/Source/BoundaryConditions/BoundaryConditions_zvel.cpp b/Source/BoundaryConditions/BoundaryConditions_zvel.cpp index 9ee5f8288..b41f5c279 100644 --- a/Source/BoundaryConditions/BoundaryConditions_zvel.cpp +++ b/Source/BoundaryConditions/BoundaryConditions_zvel.cpp @@ -13,7 +13,7 @@ using namespace amrex; * @param[in] z_phys_nd height coordinate at nodes * @param[in] bccomp index into m_domain_bcs_type */ -void ERFPhysBCFunct::impose_lateral_zvel_bcs (const Array4& dest_arr, +void ERFPhysBCFunct::impose_lateral_zvel_bcs (const Array4& dest_arr, const Array4& xvel_arr, const Array4& yvel_arr, const Box& bx, const Box& domain, @@ -181,108 +181,73 @@ void ERFPhysBCFunct::impose_vertical_zvel_bcs (const Array4& dest_arr, amrex::setBC(bx, domain, bccomp_v, 0, 1, m_domain_bcs_type, bcrs_v); amrex::setBC(bx, domain, bccomp_w, 0, 1, m_domain_bcs_type, bcrs_w); - amrex::Gpu::DeviceVector bcrs_u_d(ncomp); - amrex::Gpu::DeviceVector bcrs_v_d(ncomp); - amrex::Gpu::DeviceVector bcrs_w_d(ncomp); -#ifdef AMREX_USE_GPU - Gpu::htod_memcpy_async(bcrs_u_d.data(), bcrs_u.data(), sizeof(BCRec)); - Gpu::htod_memcpy_async(bcrs_v_d.data(), bcrs_v.data(), sizeof(BCRec)); - Gpu::htod_memcpy_async(bcrs_w_d.data(), bcrs_w.data(), sizeof(BCRec)); -#else - std::memcpy(bcrs_u_d.data(), bcrs_u.data(), sizeof(BCRec)); - std::memcpy(bcrs_v_d.data(), bcrs_v.data(), sizeof(BCRec)); - std::memcpy(bcrs_w_d.data(), bcrs_w.data(), sizeof(BCRec)); -#endif - const amrex::BCRec* bc_ptr_u = bcrs_u_d.data(); - const amrex::BCRec* bc_ptr_v = bcrs_v_d.data(); - const amrex::BCRec* bc_ptr_w = bcrs_w_d.data(); + // We use these for the asserts below + const amrex::BCRec* bc_ptr_u_h = bcrs_u.data(); + const amrex::BCRec* bc_ptr_v_h = bcrs_v.data(); + const amrex::BCRec* bc_ptr_w_h = bcrs_w.data(); bool l_use_terrain = (m_z_phys_nd != nullptr); bool l_moving_terrain = (terrain_type == TerrainType::Moving); GpuArray,AMREX_SPACEDIM+NVAR_max> l_bc_extdir_vals_d; - for (int i = 0; i < ncomp; i++) - for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) + for (int i = 0; i < ncomp; i++) { + for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) { l_bc_extdir_vals_d[i][ori] = m_bc_extdir_vals[bccomp_w+i][ori]; + } + } + // ******************************************************* + // Bottom boundary + // ******************************************************* + + // At the bottom boundary we always assert no normal flow + AMREX_ALWAYS_ASSERT(bc_ptr_w_h[0].lo(2) == ERFBCType::ext_dir); + + // Moving terrain + if (l_use_terrain && l_moving_terrain) { - Box bx_zlo(bx); bx_zlo.setBig (2,dom_lo.z-1); - Box bx_zhi(bx); bx_zhi.setSmall(2,dom_hi.z+2); + //************************************************************ + // NOTE: z_t depends on the time interval in which it is + // evaluated so we can't arbitrarily define it at a + // given time, we must specify an interval + //************************************************************ - // Populate face values on z-boundaries themselves only if EXT_DIR - Box bx_zlo_face(bx); bx_zlo_face.setSmall(2,dom_lo.z ); bx_zlo_face.setBig(2,dom_lo.z ); - Box bx_zhi_face(bx); bx_zhi_face.setSmall(2,dom_hi.z+1); bx_zhi_face.setBig(2,dom_hi.z+1); + // Static terrain + } else if (l_use_terrain) { - ParallelFor(bx_zlo, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - int n = 0; - int kflip = dom_lo.z - k; - if (bc_ptr_w[n].lo(2) == ERFBCType::ext_dir) { - dest_arr(i,j,k) = l_bc_extdir_vals_d[n][2]; - } else if (bc_ptr_w[n].lo(2) == ERFBCType::foextrap) { - dest_arr(i,j,k) = dest_arr(i,j,dom_lo.z); - } else if (bc_ptr_w[n].lo(2) == ERFBCType::reflect_even) { - dest_arr(i,j,k) = dest_arr(i,j,kflip); - } else if (bc_ptr_w[n].lo(2) == ERFBCType::reflect_odd) { - dest_arr(i,j,k) = -dest_arr(i,j,kflip); - } + AMREX_ALWAYS_ASSERT( (bc_ptr_u_h[0].lo(2) == ERFBCType::ext_dir && bc_ptr_v_h[0].lo(2) == ERFBCType::ext_dir) || + (bc_ptr_u_h[0].lo(2) != ERFBCType::ext_dir && bc_ptr_v_h[0].lo(2) != ERFBCType::ext_dir) ); + + ParallelFor(makeSlab(bx,2,dom_lo.z), [=] AMREX_GPU_DEVICE (int i, int j, int k) { + dest_arr(i,j,k) = WFromOmega(i,j,k,l_bc_extdir_vals_d[0][2],xvel_arr,yvel_arr,z_phys_nd,dxInv); }); - if (l_use_terrain && l_moving_terrain) - { - //************************************************************ - // NOTE: z_t depends on the time interval in which it is - // evaluated so we can't arbitrarily define it at a - // given time, we must specify an interval - //************************************************************ - } else if (l_use_terrain) { - ParallelFor(bx_zlo_face, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { - if (bc_ptr_w[n].lo(2) == ERFBCType::ext_dir) { - if (bc_ptr_u[n].lo(2) == ERFBCType::ext_dir && - bc_ptr_v[n].lo(2) == ERFBCType::ext_dir) { - dest_arr(i,j,k) = WFromOmega(i,j,k,l_bc_extdir_vals_d[n][2],xvel_arr,yvel_arr,z_phys_nd,dxInv); + // No terrain + } else { + ParallelFor(makeSlab(bx,2,dom_lo.z), [=] AMREX_GPU_DEVICE (int i, int j, int k) { + dest_arr(i,j,k) = l_bc_extdir_vals_d[0][2]; + }); + } - } else if (bc_ptr_u[n].lo(2) != ERFBCType::ext_dir && - bc_ptr_v[n].lo(2) != ERFBCType::ext_dir) { - dest_arr(i,j,k) = WFromOmega(i,j,k,l_bc_extdir_vals_d[n][2],xvel_arr,yvel_arr,z_phys_nd,dxInv); - } else { -#ifndef AMREX_USE_GPU - amrex::Abort("Bad combination of boundary conditions"); -#endif - } - } - }); - } else { - ParallelFor(bx_zlo_face, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { - if (bc_ptr_w[n].lo(2) == ERFBCType::ext_dir) { - dest_arr(i,j,k) = l_bc_extdir_vals_d[n][2]; - } - }); - } + // ******************************************************* + // Top boundary + // ******************************************************* - ParallelFor( - bx_zhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { - int kflip = 2*(dom_hi.z + 1) - k; - if (bc_ptr_w[n].hi(5) == ERFBCType::ext_dir) { - dest_arr(i,j,k) = l_bc_extdir_vals_d[n][5]; - } else if (bc_ptr_w[n].hi(5) == ERFBCType::foextrap) { - dest_arr(i,j,k) = dest_arr(i,j,dom_hi.z+1); - } else if (bc_ptr_w[n].hi(5) == ERFBCType::reflect_even) { - dest_arr(i,j,k) = dest_arr(i,j,kflip); - } else if (bc_ptr_w[n].hi(5) == ERFBCType::reflect_odd) { - dest_arr(i,j,k) = -dest_arr(i,j,kflip); - } + AMREX_ALWAYS_ASSERT(bc_ptr_w_h[0].hi(2) == ERFBCType::ext_dir || + bc_ptr_w_h[0].hi(2) == ERFBCType::foextrap); - }, - bx_zhi_face, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { - if (bc_ptr_w[n].hi(2) == ERFBCType::ext_dir) { - if (l_use_terrain) - dest_arr(i,j,k) = WFromOmega(i,j,k,l_bc_extdir_vals_d[n][5],xvel_arr,yvel_arr,z_phys_nd,dxInv); - else - dest_arr(i,j,k) = l_bc_extdir_vals_d[n][5]; + // NOTE: if we set SlipWall at top, that generates ERFBCType::ext_dir which sets w=0 here + // NOTE: if we set Outflow at top, that generates ERFBCType::foextrap which doesn't touch w here + if (bc_ptr_w_h[0].hi(2) == ERFBCType::ext_dir) { + ParallelFor(makeSlab(bx,2,dom_hi.z+1), [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + if (l_use_terrain) { + dest_arr(i,j,k) = WFromOmega(i,j,k,l_bc_extdir_vals_d[0][5],xvel_arr,yvel_arr,z_phys_nd,dxInv); + } else { + dest_arr(i,j,k) = l_bc_extdir_vals_d[0][5]; } - } - ); + }); } Gpu::streamSynchronize(); } diff --git a/Source/Diffusion/ComputeTurbulentViscosity.cpp b/Source/Diffusion/ComputeTurbulentViscosity.cpp index b919595ec..2a0389277 100644 --- a/Source/Diffusion/ComputeTurbulentViscosity.cpp +++ b/Source/Diffusion/ComputeTurbulentViscosity.cpp @@ -16,7 +16,8 @@ ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel, const TurbChoice& turbChoice, std::unique_ptr& most, const amrex::BCRec* bc_ptr, - bool /*vert_only*/); + bool /*vert_only*/, + const std::unique_ptr& z_phys_nd); /** * Function for computing the turbulent viscosity with LES. @@ -368,6 +369,7 @@ void ComputeTurbulentViscosity (const amrex::MultiFab& xvel , const amrex::Multi amrex::MultiFab& Hfx1, amrex::MultiFab& Hfx2, amrex::MultiFab& Hfx3, amrex::MultiFab& Diss, const amrex::Geometry& geom, const amrex::MultiFab& mapfac_u, const amrex::MultiFab& mapfac_v, + const std::unique_ptr& z_phys_nd, const TurbChoice& turbChoice, const Real const_grav, std::unique_ptr& most, const amrex::BCRec* bc_ptr, @@ -407,6 +409,6 @@ void ComputeTurbulentViscosity (const amrex::MultiFab& xvel , const amrex::Multi if (turbChoice.pbl_type != PBLType::None) { // NOTE: state_new is passed in for Cons_old (due to ptr swap in advance) ComputeTurbulentViscosityPBL(xvel, yvel, cons_in, eddyViscosity, - geom, turbChoice, most, bc_ptr, vert_only); + geom, turbChoice, most, bc_ptr, vert_only, z_phys_nd); } } diff --git a/Source/Diffusion/DiffusionSrcForState_N.cpp b/Source/Diffusion/DiffusionSrcForState_N.cpp index 14fcceddd..cfbd7c078 100644 --- a/Source/Diffusion/DiffusionSrcForState_N.cpp +++ b/Source/Diffusion/DiffusionSrcForState_N.cpp @@ -1,6 +1,6 @@ #include #include -#include +#include using namespace amrex; diff --git a/Source/Diffusion/DiffusionSrcForState_T.cpp b/Source/Diffusion/DiffusionSrcForState_T.cpp index 7bbb113dd..f427dd38b 100644 --- a/Source/Diffusion/DiffusionSrcForState_T.cpp +++ b/Source/Diffusion/DiffusionSrcForState_T.cpp @@ -1,7 +1,7 @@ #include #include #include -#include +#include using namespace amrex; @@ -565,11 +565,13 @@ DiffusionSrcForState_T (const amrex::Box& bx, const amrex::Box& domain, bool v_ext_dir_on_zhi = ( (bc_ptr[BCVars::yvel_bc].lo(5) == ERFBCType::ext_dir) ); amrex::ParallelFor(bx,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + const Real met_h_zeta = Compute_h_zeta_AtCellCenter(i,j,k,dxInv,z_nd); cell_rhs(i, j, k, qty_index) += ComputeQKESourceTerms(i,j,k,u,v,cell_data,cell_prim, mu_turb,dxInv,domain,pbl_B1_l,tm_arr(i,j,0), c_ext_dir_on_zlo, c_ext_dir_on_zhi, u_ext_dir_on_zlo, u_ext_dir_on_zhi, - v_ext_dir_on_zlo, v_ext_dir_on_zhi); + v_ext_dir_on_zlo, v_ext_dir_on_zhi, + met_h_zeta); }); } } diff --git a/Source/Diffusion/EddyViscosity.H b/Source/Diffusion/EddyViscosity.H index ddfca3c00..24332aec4 100644 --- a/Source/Diffusion/EddyViscosity.H +++ b/Source/Diffusion/EddyViscosity.H @@ -16,6 +16,7 @@ ComputeTurbulentViscosity (const amrex::MultiFab& xvel , const amrex::MultiFab& amrex::MultiFab& Hfx1, amrex::MultiFab& Hfx2, amrex::MultiFab& Hfx3, amrex::MultiFab& Diss, const amrex::Geometry& geom, const amrex::MultiFab& mapfac_u, const amrex::MultiFab& mapfac_v, + const std::unique_ptr& z_phys_nd, const TurbChoice& turbChoice, const amrex::Real const_grav, std::unique_ptr& most, const amrex::BCRec* bc_ptr, diff --git a/Source/Diffusion/Make.package b/Source/Diffusion/Make.package index fb2aa4d33..1b665ad5d 100644 --- a/Source/Diffusion/Make.package +++ b/Source/Diffusion/Make.package @@ -3,7 +3,7 @@ CEXE_sources += DiffusionSrcForMom_T.cpp CEXE_sources += DiffusionSrcForState_N.cpp CEXE_sources += DiffusionSrcForState_T.cpp - + CEXE_sources += ComputeStress_N.cpp CEXE_sources += ComputeStress_T.cpp @@ -13,8 +13,8 @@ CEXE_sources += ComputeStrain_T.cpp CEXE_sources += PBLModels.cpp CEXE_sources += NumericalDiffusion.cpp CEXE_sources += ComputeTurbulentViscosity.cpp - + CEXE_headers += Diffusion.H CEXE_headers += EddyViscosity.H CEXE_headers += NumericalDiffusion.H -CEXE_headers += ComputeQKESourceTerm.H +CEXE_headers += PBLModels.H diff --git a/Source/Diffusion/ComputeQKESourceTerm.H b/Source/Diffusion/PBLModels.H similarity index 71% rename from Source/Diffusion/ComputeQKESourceTerm.H rename to Source/Diffusion/PBLModels.H index 3586c1d15..4405a3ab3 100644 --- a/Source/Diffusion/ComputeQKESourceTerm.H +++ b/Source/Diffusion/PBLModels.H @@ -2,44 +2,32 @@ #define _PBLMODELS_H_ /** - * Function for computing the QKE source terms. + * Function for computing vertical derivatives for use in PBL model * * @param[in] u velocity in x-dir * @param[in] v velocity in y-dir * @param[in] cell_data conserved cell center vars - * @param[in] cell_prim primitive cell center vars - * @param[in] K_turb turbulent viscosity - * @param[in] cellSizeInv inverse cell size array - * @param[in] domain box of the whole domain - * @param[in] pbl_B1_l a parameter - * @param[in] theta_mean average theta */ AMREX_GPU_DEVICE AMREX_FORCE_INLINE -amrex::Real -ComputeQKESourceTerms (int i, int j, int k, - const amrex::Array4& uvel, - const amrex::Array4& vvel, - const amrex::Array4& cell_data, - const amrex::Array4& cell_prim, - const amrex::Array4& K_turb, - const amrex::GpuArray& cellSizeInv, - const amrex::Box& domain, - amrex::Real pbl_B1_l, - const amrex::Real theta_mean, - bool c_ext_dir_on_zlo, - bool c_ext_dir_on_zhi, - bool u_ext_dir_on_zlo, - bool u_ext_dir_on_zhi, - bool v_ext_dir_on_zlo, - bool v_ext_dir_on_zhi) +void +ComputeVerticalDerivativesPBL(int i, int j, int k, + const amrex::Array4& uvel, + const amrex::Array4& vvel, + const amrex::Array4& cell_data, + const int izmin, + const int izmax, + const amrex::Real& dz_inv, + const bool c_ext_dir_on_zlo, + const bool c_ext_dir_on_zhi, + const bool u_ext_dir_on_zlo, + const bool u_ext_dir_on_zhi, + const bool v_ext_dir_on_zlo, + const bool v_ext_dir_on_zhi, + amrex::Real& dthetadz, + amrex::Real& dudz, + amrex::Real& dvdz) { - // Compute some relevant derivatives - amrex::Real dz_inv = cellSizeInv[2]; - int izmin = domain.smallEnd(2); - int izmax = domain.bigEnd(2); - amrex::Real dthetadz, dudz, dvdz; - amrex::Real source_term = 0.0; if ( k==izmax && c_ext_dir_on_zhi ) { dthetadz = (1.0/3.0)*(-cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp) - 3.0 * cell_data(i,j,k ,RhoTheta_comp)/cell_data(i,j,k ,Rho_comp) @@ -72,6 +60,56 @@ ComputeQKESourceTerms (int i, int j, int k, } else { dvdz = 0.25*( vvel(i,j,k+1) - vvel(i,j,k-1) + vvel(i,j+1,k+1) - vvel(i,j+1,k-1) )*dz_inv; } +} + +/** + * Function for computing the QKE source terms. + * + * @param[in] u velocity in x-dir + * @param[in] v velocity in y-dir + * @param[in] cell_data conserved cell center vars + * @param[in] cell_prim primitive cell center vars + * @param[in] K_turb turbulent viscosity + * @param[in] cellSizeInv inverse cell size array + * @param[in] domain box of the whole domain + * @param[in] pbl_B1_l a parameter + * @param[in] theta_mean average theta + */ +AMREX_GPU_DEVICE +AMREX_FORCE_INLINE +amrex::Real +ComputeQKESourceTerms (int i, int j, int k, + const amrex::Array4& uvel, + const amrex::Array4& vvel, + const amrex::Array4& cell_data, + const amrex::Array4& cell_prim, + const amrex::Array4& K_turb, + const amrex::GpuArray& cellSizeInv, + const amrex::Box& domain, + amrex::Real pbl_B1_l, + const amrex::Real theta_mean, + bool c_ext_dir_on_zlo, + bool c_ext_dir_on_zhi, + bool u_ext_dir_on_zlo, + bool u_ext_dir_on_zhi, + bool v_ext_dir_on_zlo, + bool v_ext_dir_on_zhi, + const amrex::Real met_h_zeta=1.0) +{ + // Compute some relevant derivatives + amrex::Real dthetadz, dudz, dvdz; + amrex::Real source_term = 0.0; + + amrex::Real dz_inv = cellSizeInv[2]; + int izmin = domain.smallEnd(2); + int izmax = domain.bigEnd(2); + + ComputeVerticalDerivativesPBL(i, j, k, + uvel, vvel, cell_data, izmin, izmax, dz_inv/met_h_zeta, + c_ext_dir_on_zlo, c_ext_dir_on_zhi, + u_ext_dir_on_zlo, u_ext_dir_on_zhi, + v_ext_dir_on_zlo, v_ext_dir_on_zhi, + dthetadz, dudz, dvdz); // Production (We store mu_turb, which is 0.5*K_turb) source_term += 4.0*K_turb(i,j,k,EddyDiff::Mom_v) * (dudz*dudz + dvdz*dvdz); diff --git a/Source/Diffusion/PBLModels.cpp b/Source/Diffusion/PBLModels.cpp index 1c0b21628..dccbb9bb2 100644 --- a/Source/Diffusion/PBLModels.cpp +++ b/Source/Diffusion/PBLModels.cpp @@ -3,6 +3,7 @@ #include "Diffusion.H" #include "ERF_Constants.H" #include "TurbStruct.H" +#include "PBLModels.H" /** * Function to compute turbulent viscosity with PBL. @@ -24,8 +25,11 @@ ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel, const TurbChoice& turbChoice, std::unique_ptr& most, const amrex::BCRec* bc_ptr, - bool /*vert_only*/) + bool /*vert_only*/, + const std::unique_ptr& z_phys_nd) { + const bool use_terrain = (z_phys_nd != nullptr); + // MYNN Level 2.5 PBL Model if (turbChoice.pbl_type == PBLType::MYNN25) { @@ -79,22 +83,45 @@ ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel, const amrex::Array4 qvel = qturb.array(); const amrex::Array4 qvel_old = qturb_old.array(); - amrex::ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), bx, - [=] AMREX_GPU_DEVICE (int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept - { - qvel(i,j,k) = std::sqrt(cell_data(i,j,k,RhoQKE_comp) / cell_data(i,j,k,Rho_comp)); - qvel_old(i,j,k) = std::sqrt(cell_data(i,j,k,RhoQKE_comp) / cell_data(i,j,k,Rho_comp) + eps); - AMREX_ASSERT_WITH_MESSAGE(qvel(i,j,k) > 0.0, "QKE must have a positive value"); - AMREX_ASSERT_WITH_MESSAGE(qvel_old(i,j,k) > 0.0, "Old QKE must have a positive value"); - - const amrex::Real Zval = gdata.ProbLo(2) + (k + 0.5)*gdata.CellSize(2); - if (sbx.contains(i,j,k)) { - amrex::Gpu::deviceReduceSum(&qint(i,j,0,0), Zval*qvel(i,j,k), handler); - amrex::Gpu::deviceReduceSum(&qint(i,j,0,1), qvel(i,j,k), handler); - } - }); + // vertical integrals to compute lengthscale + if (use_terrain) { + const amrex::Array4 &z_nd = z_phys_nd->array(mfi); + const auto invCellSize = geom.InvCellSizeArray(); + amrex::ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), bx, + [=] AMREX_GPU_DEVICE (int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept + { + qvel(i,j,k) = std::sqrt(cell_data(i,j,k,RhoQKE_comp) / cell_data(i,j,k,Rho_comp)); + qvel_old(i,j,k) = std::sqrt(cell_data(i,j,k,RhoQKE_comp) / cell_data(i,j,k,Rho_comp) + eps); + AMREX_ASSERT_WITH_MESSAGE(qvel(i,j,k) > 0.0, "QKE must have a positive value"); + AMREX_ASSERT_WITH_MESSAGE(qvel_old(i,j,k) > 0.0, "Old QKE must have a positive value"); + + if (sbx.contains(i,j,k)) { + const amrex::Real Zval = Compute_Zrel_AtCellCenter(i,j,k,z_nd); + const amrex::Real dz = Compute_h_zeta_AtCellCenter(i,j,k,invCellSize,z_nd); + amrex::Gpu::deviceReduceSum(&qint(i,j,0,0), Zval*qvel(i,j,k)*dz, handler); + amrex::Gpu::deviceReduceSum(&qint(i,j,0,1), qvel(i,j,k)*dz, handler); + } + }); + } else { + amrex::ParallelFor(amrex::Gpu::KernelInfo().setReduction(true), bx, + [=] AMREX_GPU_DEVICE (int i, int j, int k, amrex::Gpu::Handler const& handler) noexcept + { + qvel(i,j,k) = std::sqrt(cell_data(i,j,k,RhoQKE_comp) / cell_data(i,j,k,Rho_comp)); + qvel_old(i,j,k) = std::sqrt(cell_data(i,j,k,RhoQKE_comp) / cell_data(i,j,k,Rho_comp) + eps); + AMREX_ASSERT_WITH_MESSAGE(qvel(i,j,k) > 0.0, "QKE must have a positive value"); + AMREX_ASSERT_WITH_MESSAGE(qvel_old(i,j,k) > 0.0, "Old QKE must have a positive value"); + + // Not multiplying by dz: its constant and would fall out when we divide qint0/qint1 anyway + if (sbx.contains(i,j,k)) { + const amrex::Real Zval = gdata.ProbLo(2) + (k + 0.5)*gdata.CellSize(2); + amrex::Gpu::deviceReduceSum(&qint(i,j,0,0), Zval*qvel(i,j,k), handler); + amrex::Gpu::deviceReduceSum(&qint(i,j,0,1), qvel(i,j,k), handler); + } + }); + } amrex::Real dz_inv = geom.InvCellSize(2); + const auto& dxInv = geom.InvCellSizeArray(); int izmin = geom.Domain().smallEnd(2); int izmax = geom.Domain().bigEnd(2); @@ -110,43 +137,23 @@ ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel, const auto& u_star_arr = u_star_mf->array(mfi); const auto& t_star_arr = t_star_mf->array(mfi); + const amrex::Array4 z_nd; + if (use_terrain) { + const auto& z_nd = z_phys_nd->array(mfi); + } + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { // Compute some partial derivatives that we will need (second order) // U and V derivatives are interpolated to account for staggered grid + const amrex::Real met_h_zeta = use_terrain ? Compute_h_zeta_AtCellCenter(i,j,k,dxInv,z_nd) : 1.0; amrex::Real dthetadz, dudz, dvdz; - if ( k==izmax && c_ext_dir_on_zhi ) { - dthetadz = (1.0/3.0)*(-cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp) - - 3.0 * cell_data(i,j,k ,RhoTheta_comp)/cell_data(i,j,k ,Rho_comp) - + 4.0 * cell_data(i,j,k+1,RhoTheta_comp)/cell_data(i,j,k+1,Rho_comp) )*dz_inv; - } else if ( k==izmin && c_ext_dir_on_zlo ) { - dthetadz = (1.0/3.0)*( cell_data(i,j,k+1,RhoTheta_comp)/cell_data(i,j,k+1,Rho_comp) - + 3.0 * cell_data(i,j,k ,RhoTheta_comp)/cell_data(i,j,k ,Rho_comp) - - 4.0 * cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp) )*dz_inv; - } else { - dthetadz = 0.5*(cell_data(i,j,k+1,RhoTheta_comp)/cell_data(i,j,k+1,Rho_comp) - - cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp))*dz_inv; - } - - if ( k==izmax && u_ext_dir_on_zhi ) { - dudz = (1.0/6.0)*( (-uvel(i ,j,k-1) - 3.0 * uvel(i ,j,k ) + 4.0 * uvel(i ,j,k+1)) - + (-uvel(i+1,j,k-1) - 3.0 * uvel(i+1,j,k ) + 4.0 * uvel(i+1,j,k+1)) )*dz_inv; - } else if ( k==izmin && u_ext_dir_on_zlo ) { - dudz = (1.0/6.0)*( (uvel(i ,j,k+1) + 3.0 * uvel(i ,j,k ) - 4.0 * uvel(i ,j,k-1)) - + (uvel(i+1,j,k+1) + 3.0 * uvel(i+1,j,k ) - 4.0 * uvel(i+1,j,k-1)) )*dz_inv; - } else { - dudz = 0.25*(uvel(i,j,k+1) - uvel(i,j,k-1) + uvel(i+1,j,k+1) - uvel(i+1,j,k-1))*dz_inv; - } - - if ( k==izmax && v_ext_dir_on_zhi ) { - dvdz = (1.0/6.0)*( (-vvel(i,j ,k-1) - 3.0 * vvel(i,j ,k ) + 4.0 * vvel(i,j ,k+1)) - + (-vvel(i,j+1,k-1) - 3.0 * vvel(i,j+1,k ) + 4.0 * vvel(i,j+1,k+1)) )*dz_inv; - } else if ( k==izmin && v_ext_dir_on_zlo ) { - dvdz = (1.0/6.0)*( (vvel(i,j ,k+1) + 3.0 * vvel(i,j ,k ) - 4.0 * vvel(i,j ,k-1)) - + (vvel(i,j+1,k+1) + 3.0 * vvel(i,j+1,k ) - 4.0 * vvel(i,j+1,k-1)) )*dz_inv; - } else { - dvdz = 0.25*(vvel(i,j,k+1) - vvel(i,j,k-1) + vvel(i,j+1,k+1) - vvel(i,j+1,k-1))*dz_inv; - } + ComputeVerticalDerivativesPBL(i, j, k, + uvel, vvel, cell_data, izmin, izmax, dz_inv/met_h_zeta, + c_ext_dir_on_zlo, c_ext_dir_on_zhi, + u_ext_dir_on_zlo, u_ext_dir_on_zhi, + v_ext_dir_on_zlo, v_ext_dir_on_zhi, + dthetadz, dudz, dvdz); // Spatially varying MOST amrex::Real surface_heat_flux = -u_star_arr(i,j,0) * t_star_arr(i,j,0); diff --git a/Source/TimeIntegration/ERF_advance_dycore.cpp b/Source/TimeIntegration/ERF_advance_dycore.cpp index f0ec29c24..80f74eecb 100644 --- a/Source/TimeIntegration/ERF_advance_dycore.cpp +++ b/Source/TimeIntegration/ERF_advance_dycore.cpp @@ -220,6 +220,7 @@ void ERF::advance_dycore(int level, state_old[IntVar::cons], *eddyDiffs, *Hfx1, *Hfx2, *Hfx3, *Diss, // to be updated fine_geom, *mapfac_u[level], *mapfac_v[level], + z_phys_nd[level], tc, solverChoice.gravity, m_most, bc_ptr_d); } diff --git a/Source/Utils/TerrainMetrics.H b/Source/Utils/TerrainMetrics.H index 0d8faa9af..4f8e79ab3 100644 --- a/Source/Utils/TerrainMetrics.H +++ b/Source/Utils/TerrainMetrics.H @@ -337,6 +337,26 @@ Compute_h_eta_AtEdgeCenterI (const int &i, const int &j, const int &k, return met_h_eta; } +// Relative height above terrain surface at cell center from z_nd (nodal absolute height) +AMREX_GPU_DEVICE +AMREX_FORCE_INLINE +amrex::Real +Compute_Zrel_AtCellCenter (const int &i, const int &j, const int &k, + const amrex::Array4& z_nd) +{ + const amrex::Real z_cc = 0.125*( z_nd(i ,j ,k ) + z_nd(i ,j ,k+1) + + + z_nd(i+1,j ,k ) + z_nd(i ,j ,k+1) + + z_nd(i ,j+1,k ) + z_nd(i ,j+1,k+1) + + z_nd(i+1,j+1,k ) + z_nd(i ,j+1,k+1)); + + // Note: we assume the z_nd array spans from the bottom to top of the domain + // i.e. no domain decomposition across processors in vertical direction + const amrex::Real z0_cc = 0.25*( z_nd(i ,j ,0) + z_nd(i ,j+1,0) + + z_nd(i+1,j ,0) + z_nd(i+1,j+1,0)); + + return (z_cc - z0_cc); +} + /** * Define omega given u,v and w */