Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add terrain/stretched grid capability for PBL models #1348

Merged
merged 8 commits into from
Dec 28, 2023
Merged

Add terrain/stretched grid capability for PBL models #1348

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
9b2ad72
interface changes to prepare PBL for terrain
baperry2 Dec 14, 2023
a3f95af
terrain metrics for PBL model
baperry2 Dec 15, 2023
84a8b59
remove barrier
baperry2 Dec 15, 2023
786d305
Merge branch 'development' into pbl-terrain
asalmgren Dec 15, 2023
e2f00d3
Merge branch 'development' into pbl-terrain
baperry2 Dec 15, 2023
bd99c24
Merge branch 'development' into pbl-terrain
baperry2 Dec 28, 2023
2cfccf6
Merge branch 'pbl-terrain' of github.com:baperry2/ERF into pbl-terrain
baperry2 Dec 28, 2023
2d91670
some indent width formatting
baperry2 Dec 28, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 4 additions & 2 deletions Source/Diffusion/ComputeTurbulentViscosity.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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);
}
}
2 changes: 1 addition & 1 deletion Source/Diffusion/DiffusionSrcForState_N.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
#include <Diffusion.H>
#include <EddyViscosity.H>
#include <ComputeQKESourceTerm.H>
#include <PBLModels.H>

using namespace amrex;

Expand Down
6 changes: 4 additions & 2 deletions Source/Diffusion/DiffusionSrcForState_T.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
#include <Diffusion.H>
#include <EddyViscosity.H>
#include <TerrainMetrics.H>
#include <ComputeQKESourceTerm.H>
#include <PBLModels.H>

using namespace amrex;

Expand Down Expand Up @@ -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);
});
}
}
1 change: 1 addition & 0 deletions Source/Diffusion/EddyViscosity.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
6 changes: 3 additions & 3 deletions Source/Diffusion/Make.package
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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
98 changes: 68 additions & 30 deletions Source/Diffusion/ComputeQKESourceTerm.H → Source/Diffusion/PBLModels.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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);
Expand Down
101 changes: 54 additions & 47 deletions Source/Diffusion/PBLModels.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
#include "Diffusion.H"
#include "ERF_Constants.H"
#include "TurbStruct.H"
#include "PBLModels.H"

/**
* Function to compute turbulent viscosity with PBL.
Expand All @@ -24,8 +25,11 @@ 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)
{
const bool use_terrain = (z_phys_nd != nullptr);

// MYNN Level 2.5 PBL Model
if (turbChoice.pbl_type == PBLType::MYNN25) {

Expand Down Expand Up @@ -79,22 +83,45 @@ ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel,
const amrex::Array4<amrex::Real> qvel = qturb.array();
const amrex::Array4<amrex::Real> 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<amrex::Real const> &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);

Expand All @@ -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<amrex::Real const> 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);
Expand Down
1 change: 1 addition & 0 deletions Source/TimeIntegration/ERF_advance_dycore.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
}

Expand Down
20 changes: 20 additions & 0 deletions Source/Utils/TerrainMetrics.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<const amrex::Real>& 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
*/
Expand Down
Loading