Merge branch 'development' into AddMoistScalars

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<Real>& dest_arr,
+void ERFPhysBCFunct::impose_lateral_zvel_bcs (const Array4<Real >& dest_arr,
  const Array4<Real const>& xvel_arr,
  const Array4<Real const>& yvel_arr,
  const Box& bx, const Box& domain,
@@ -181,108 +181,73 @@ void ERFPhysBCFunct::impose_vertical_zvel_bcs (const Array4<Real>& 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<BCRec> bcrs_u_d(ncomp);
- amrex::Gpu::DeviceVector<BCRec> bcrs_v_d(ncomp);
- amrex::Gpu::DeviceVector<BCRec> 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<GpuArray<Real, AMREX_SPACEDIM*2>,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();
 }

Source/Diffusion/ComputeTurbulentViscosity.cpp

@@ -16,7 +16,8 @@ ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel,
  const TurbChoice& turbChoice,
  std::unique_ptr<ABLMost>& most,
  const amrex::BCRec* bc_ptr,
- bool /*vert_only*/);
+ bool /*vert_only*/,
+ const std::unique_ptr<amrex::MultiFab>& 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<amrex::MultiFab>& z_phys_nd,
  const TurbChoice& turbChoice, const Real const_grav,
  std::unique_ptr<ABLMost>& 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);
  }
 }

Source/Diffusion/DiffusionSrcForState_N.cpp

@@ -1,6 +1,6 @@
 #include <Diffusion.H>
 #include <EddyViscosity.H>
-#include <ComputeQKESourceTerm.H>
+#include <PBLModels.H>
 
 using namespace amrex;
 

Source/Diffusion/DiffusionSrcForState_T.cpp

@@ -1,7 +1,7 @@
 #include <Diffusion.H>
 #include <EddyViscosity.H>
 #include <TerrainMetrics.H>
-#include <ComputeQKESourceTerm.H>
+#include <PBLModels.H>
 
 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);
  });
  }
 }

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<amrex::MultiFab>& z_phys_nd,
  const TurbChoice& turbChoice, const amrex::Real const_grav,
  std::unique_ptr<ABLMost>& most,
  const amrex::BCRec* bc_ptr,

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

Source/Diffusion/ComputeQKESourceTerm.H → 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<const amrex::Real>& uvel,
- const amrex::Array4<const amrex::Real>& vvel,
- const amrex::Array4<const amrex::Real>& cell_data,
- const amrex::Array4<const amrex::Real>& cell_prim,
- const amrex::Array4<const amrex::Real>& K_turb,
- const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& 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<const amrex::Real>& uvel,
+  const amrex::Array4<const amrex::Real>& vvel,
+  const amrex::Array4<const amrex::Real>& 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<const amrex::Real>& uvel,
+ const amrex::Array4<const amrex::Real>& vvel,
+ const amrex::Array4<const amrex::Real>& cell_data,
+ const amrex::Array4<const amrex::Real>& cell_prim,
+ const amrex::Array4<const amrex::Real>& K_turb,
+ const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& 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);