From 9b2ad7238f4e3feed4becfac52afba507c7e28f7 Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Thu, 14 Dec 2023 15:19:07 -0700 Subject: [PATCH 1/6] interface changes to prepare PBL for terrain --- .../Diffusion/ComputeTurbulentViscosity.cpp | 6 +- Source/Diffusion/DiffusionSrcForState_N.cpp | 2 +- Source/Diffusion/DiffusionSrcForState_T.cpp | 2 +- Source/Diffusion/EddyViscosity.H | 1 + Source/Diffusion/Make.package | 6 +- .../{ComputeQKESourceTerm.H => PBLModels.H} | 97 +++++++++++++------ Source/Diffusion/PBLModels.cpp | 47 +++------ Source/TimeIntegration/ERF_advance_dycore.cpp | 1 + Source/Utils/TerrainMetrics.H | 20 ++++ 9 files changed, 112 insertions(+), 70 deletions(-) rename Source/Diffusion/{ComputeQKESourceTerm.H => PBLModels.H} (72%) diff --git a/Source/Diffusion/ComputeTurbulentViscosity.cpp b/Source/Diffusion/ComputeTurbulentViscosity.cpp index 1a65ac245..e74fa9f8f 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 da68d7ac9..b1d4065ab 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..5d044fb44 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; 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 72% rename from Source/Diffusion/ComputeQKESourceTerm.H rename to Source/Diffusion/PBLModels.H index 1bfa41ac0..956ca505d 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_N(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,55 @@ 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) +{ + // 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_N(i, j, k, + uvel, vvel, cell_data, izmin, izmax, dz_inv, + 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 a9aec0c5b..bf406655c 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,14 @@ 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); + + // DEBUG DEBUG DEBUG + if (use_terrain) {amrex::Abort("PBL model not yet supported with Terrain");} + // MYNN Level 2.5 PBL Model if (turbChoice.pbl_type == PBLType::MYNN25) { @@ -115,38 +122,12 @@ ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel, // Compute some partial derivatives that we will need (second order) // U and V derivatives are interpolated to account for staggered grid 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_N(i, j, k, + uvel, vvel, cell_data, izmin, izmax, dz_inv, + 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 3ed147cc6..68327b12d 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..f52230a76 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 */ From a3f95afecfbed3148cd6cd0a3275145a5c3d5d6f Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Thu, 14 Dec 2023 17:57:31 -0700 Subject: [PATCH 2/6] terrain metrics for PBL model --- Source/Diffusion/DiffusionSrcForState_T.cpp | 4 +- Source/Diffusion/PBLModels.H | 47 +++++++++--------- Source/Diffusion/PBLModels.cpp | 55 ++++++++++++++++----- 3 files changed, 69 insertions(+), 37 deletions(-) diff --git a/Source/Diffusion/DiffusionSrcForState_T.cpp b/Source/Diffusion/DiffusionSrcForState_T.cpp index 5d044fb44..f427dd38b 100644 --- a/Source/Diffusion/DiffusionSrcForState_T.cpp +++ b/Source/Diffusion/DiffusionSrcForState_T.cpp @@ -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/PBLModels.H b/Source/Diffusion/PBLModels.H index 956ca505d..44762af27 100644 --- a/Source/Diffusion/PBLModels.H +++ b/Source/Diffusion/PBLModels.H @@ -11,22 +11,22 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE void -ComputeVerticalDerivativesPBL_N(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) +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) { 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) @@ -93,7 +93,8 @@ ComputeQKESourceTerms (int i, int j, int k, bool u_ext_dir_on_zlo, bool u_ext_dir_on_zhi, bool v_ext_dir_on_zlo, - bool v_ext_dir_on_zhi) + bool v_ext_dir_on_zhi, + const amrex::Real met_h_zeta=1.0) { // Compute some relevant derivatives amrex::Real dthetadz, dudz, dvdz; @@ -103,12 +104,12 @@ ComputeQKESourceTerms (int i, int j, int k, int izmin = domain.smallEnd(2); int izmax = domain.bigEnd(2); - ComputeVerticalDerivativesPBL_N(i, j, k, - uvel, vvel, cell_data, izmin, izmax, dz_inv, - 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); + 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 bf406655c..c78d42de7 100644 --- a/Source/Diffusion/PBLModels.cpp +++ b/Source/Diffusion/PBLModels.cpp @@ -86,22 +86,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 - { + // 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"); - 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); + 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); @@ -117,17 +140,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; - ComputeVerticalDerivativesPBL_N(i, j, k, - uvel, vvel, cell_data, izmin, izmax, dz_inv, - 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); + 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); From 84a8b598c392535088be2a853b5a28b2eab497c0 Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Thu, 14 Dec 2023 18:02:19 -0700 Subject: [PATCH 3/6] remove barrier --- Source/Diffusion/PBLModels.cpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/Source/Diffusion/PBLModels.cpp b/Source/Diffusion/PBLModels.cpp index c78d42de7..4d99744b6 100644 --- a/Source/Diffusion/PBLModels.cpp +++ b/Source/Diffusion/PBLModels.cpp @@ -30,9 +30,6 @@ ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel, { const bool use_terrain = (z_phys_nd != nullptr); - // DEBUG DEBUG DEBUG - if (use_terrain) {amrex::Abort("PBL model not yet supported with Terrain");} - // MYNN Level 2.5 PBL Model if (turbChoice.pbl_type == PBLType::MYNN25) { From 2d916704d6cc96836c5e0287aad7880638e75b0d Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Thu, 28 Dec 2023 11:42:46 -0700 Subject: [PATCH 4/6] some indent width formatting --- Source/Utils/TerrainMetrics.H | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/Source/Utils/TerrainMetrics.H b/Source/Utils/TerrainMetrics.H index f52230a76..4f8e79ab3 100644 --- a/Source/Utils/TerrainMetrics.H +++ b/Source/Utils/TerrainMetrics.H @@ -344,17 +344,17 @@ 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)); + 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)); + // 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); + return (z_cc - z0_cc); } /** From 5577831c7c3fa0274ebd65043271018d2c0ebe0c Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Sat, 30 Dec 2023 13:55:55 -0800 Subject: [PATCH 5/6] Cleanup of zvel vertical bc's based on Mahesh's bug finds (#1352) --- .../BoundaryConditions_zvel.cpp | 116 +++++++----------- 1 file changed, 46 insertions(+), 70 deletions(-) diff --git a/Source/BoundaryConditions/BoundaryConditions_zvel.cpp b/Source/BoundaryConditions/BoundaryConditions_zvel.cpp index 9ee5f8288..c53a8e860 100644 --- a/Source/BoundaryConditions/BoundaryConditions_zvel.cpp +++ b/Source/BoundaryConditions/BoundaryConditions_zvel.cpp @@ -202,87 +202,63 @@ void ERFPhysBCFunct::impose_vertical_zvel_bcs (const Array4& dest_arr, 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[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[0].lo(2) == ERFBCType::ext_dir && bc_ptr_v[0].lo(2) == ERFBCType::ext_dir) || + (bc_ptr_u[0].lo(2) != ERFBCType::ext_dir && bc_ptr_v[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[0].hi(2) == ERFBCType::ext_dir || + bc_ptr_w[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[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(); } From b649ffbab4ca329b9074c0fa0813816c9f5e2cb7 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Sat, 30 Dec 2023 14:11:23 -0800 Subject: [PATCH 6/6] fix for GPU (#1353) --- .../BoundaryConditions_zvel.cpp | 33 +++++++------------ 1 file changed, 11 insertions(+), 22 deletions(-) diff --git a/Source/BoundaryConditions/BoundaryConditions_zvel.cpp b/Source/BoundaryConditions/BoundaryConditions_zvel.cpp index c53a8e860..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,21 +181,10 @@ 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); @@ -213,7 +202,7 @@ void ERFPhysBCFunct::impose_vertical_zvel_bcs (const Array4& dest_arr, // ******************************************************* // At the bottom boundary we always assert no normal flow - AMREX_ALWAYS_ASSERT(bc_ptr_w[0].lo(2) == ERFBCType::ext_dir); + AMREX_ALWAYS_ASSERT(bc_ptr_w_h[0].lo(2) == ERFBCType::ext_dir); // Moving terrain if (l_use_terrain && l_moving_terrain) @@ -227,8 +216,8 @@ void ERFPhysBCFunct::impose_vertical_zvel_bcs (const Array4& dest_arr, // Static terrain } else if (l_use_terrain) { - AMREX_ALWAYS_ASSERT( (bc_ptr_u[0].lo(2) == ERFBCType::ext_dir && bc_ptr_v[0].lo(2) == ERFBCType::ext_dir) || - (bc_ptr_u[0].lo(2) != ERFBCType::ext_dir && bc_ptr_v[0].lo(2) != ERFBCType::ext_dir) ); + 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); @@ -245,12 +234,12 @@ void ERFPhysBCFunct::impose_vertical_zvel_bcs (const Array4& dest_arr, // Top boundary // ******************************************************* - AMREX_ALWAYS_ASSERT(bc_ptr_w[0].hi(2) == ERFBCType::ext_dir || - bc_ptr_w[0].hi(2) == ERFBCType::foextrap); + AMREX_ALWAYS_ASSERT(bc_ptr_w_h[0].hi(2) == ERFBCType::ext_dir || + bc_ptr_w_h[0].hi(2) == ERFBCType::foextrap); // 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[0].hi(2) == ERFBCType::ext_dir) { + 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) {