ERF_ScalarDiffusion.H File Reference

#include <type_traits>
#include "AMReX_Box.H"
#include "AMReX_FArrayBox.H"
#include "AMReX_Gpu.H"
#include "ERF_DiffStruct.H"
#include "ERF_NumericalConstants.H"
#include "ERF_TerrainMetrics.H"
#include "ERF_TerrainImplicitMetric.H"

Include dependency graph for ERF_ScalarDiffusion.H:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Classes
struct	ScalarDiffusionCoefficients
struct	ScalarDiffusionFieldViews
struct	ScalarDiffusionCoefficientMode
struct	ScalarDiffusionPrescribedFluxes
struct	ScalarDiffusionRawFluxDiagnostic
struct	ScalarDiffusionSurfaceFaces
struct	ScalarDiffusionFluxPolicy
struct	StretchedDiffusionBoundaryWeights

Functions
template<bool MultiplyMolecularByDensity, bool AddTurbulence>
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	ScalarDiffusionFaceCoefficient (const amrex::Array4< const amrex::Real > &rho, const int rho_comp, const amrex::Array4< const amrex::Real > &mu_turb, const ScalarDiffusionCoefficients &coefficients, const int i, const int j, const int k, const int di, const int dj, const int dk, const int eddy_comp) noexcept
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE StretchedDiffusionBoundaryWeights	StretchedDiffusionDirichletWeights (const amrex::Real dz_near, const amrex::Real dz_next) noexcept
template<int Direction>
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	ScalarDiffusionFlux_N (const amrex::Array4< const amrex::Real > &scalar, const int scalar_comp, const int i, const int j, const int k, const amrex::Real coeff, const amrex::Real inv_spacing, const amrex::Real map_numerator, const amrex::Real map_denominator, const bool ext_dir_lo, const bool ext_dir_hi) noexcept
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	StretchedScalarGradient (const amrex::Array4< const amrex::Real > &scalar, const int scalar_comp, const int i, const int j, const int k, const amrex::Real *dz, const int klo, const int khi, const bool ext_dir_lo, const bool ext_dir_hi) noexcept
template<int Direction>
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	ScalarDiffusionFlux_S_Horizontal (const amrex::Array4< const amrex::Real > &scalar, const int scalar_comp, const int i, const int j, const int k, const amrex::Real coeff, const amrex::Real inv_spacing, const amrex::Real face_map) noexcept
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	ScalarDiffusionFlux_Tx (const amrex::Real coeff, const amrex::Real face_map, const amrex::Real grad_xi, const amrex::Real h_xi, const amrex::Real grad_z) noexcept
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	ScalarDiffusionFlux_Ty (const amrex::Real coeff, const amrex::Real face_map, const amrex::Real grad_eta, const amrex::Real h_eta, const amrex::Real grad_z) noexcept
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	TerrainScalarDiffusionFluxAtIFace (const amrex::Array4< const amrex::Real > &scalar, const int scalar_comp, const amrex::Real coeff, const int i, const int j, const int k, const amrex::Array4< const amrex::Real > &z_nd, const amrex::Array4< const amrex::Real > &z_cc, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &mf_ux) noexcept
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	TerrainScalarDiffusionFluxAtJFace (const amrex::Array4< const amrex::Real > &scalar, const int scalar_comp, const amrex::Real coeff, const int i, const int j, const int k, const amrex::Array4< const amrex::Real > &z_nd, const amrex::Array4< const amrex::Real > &z_cc, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &mf_vy) noexcept
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	TerrainScalarDiffusionFluxAtKFace (const amrex::Array4< const amrex::Real > &scalar, const int scalar_comp, const amrex::Real coeff, const int i, const int j, const int kface, const amrex::Array4< const amrex::Real > &z_nd, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const bool ext_dir_lo, const bool ext_dir_hi) noexcept
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	TerrainDiffusionGz (const amrex::Real raw_Fz, const amrex::Real mx, const amrex::Real h_xi, const amrex::Real Fbar_x, const amrex::Real my, const amrex::Real h_eta, const amrex::Real Fbar_y) noexcept
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	TerrainDiffusionMappedTx (const amrex::Real raw_Fx, const amrex::Real h_zeta, const amrex::Real my) noexcept
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	TerrainDiffusionMappedTy (const amrex::Real raw_Fy, const amrex::Real h_zeta, const amrex::Real mx) noexcept
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	ScalarDiffusionMappedHorizontalTransfer (const amrex::Real raw_flux, const amrex::Real face_area) noexcept
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	ScalarDiffusionMappedVerticalTransfer (const amrex::Real raw_flux, const amrex::Real mx, const amrex::Real my) noexcept
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	TerrainScalarDiffusionGzAtKFace (const int i, const int j, const int kface, const amrex::Box &domain, const amrex::Array4< const amrex::Real > &raw_xflux, const amrex::Array4< const amrex::Real > &raw_yflux, const amrex::Array4< const amrex::Real > &raw_zflux, const int flux_comp, const amrex::Array4< const amrex::Real > &z_nd, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &mf_mx, const amrex::Array4< const amrex::Real > &mf_my) noexcept
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	TerrainScalarDiffusionMappedTransferAtKFace (const int i, const int j, const int kface, const amrex::Box &domain, const amrex::Array4< const amrex::Real > &raw_xflux, const amrex::Array4< const amrex::Real > &raw_yflux, const amrex::Array4< const amrex::Real > &raw_zflux, const int flux_comp, const amrex::Array4< const amrex::Real > &z_nd, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &mf_mx, const amrex::Array4< const amrex::Real > &mf_my, const bool suppress_transfer=false) noexcept
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	ScalarDiffusionMappedFluxDivergence (const amrex::Real fx_hi, const amrex::Real fx_lo, const amrex::Real fy_hi, const amrex::Real fy_lo, const amrex::Real fz_hi, const amrex::Real fz_lo, const amrex::Real dx_inv, const amrex::Real dy_inv, const amrex::Real dzeta_inv, const amrex::Real mx, const amrex::Real my, const amrex::Real detJ_cell) noexcept
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	TerrainDiffusionDivergence_T (const amrex::Real mapped_Fx_hi, const amrex::Real mapped_Fx_lo, const amrex::Real mapped_Fy_hi, const amrex::Real mapped_Fy_lo, const amrex::Real Gz_hi, const amrex::Real Gz_lo, const amrex::Real dx_inv, const amrex::Real dy_inv, const amrex::Real dz_inv, const amrex::Real mfsq, const amrex::Real detJ_cell) noexcept
template<bool MultiplyMolecularByDensity, bool AddTurbulence>
void	BuildScalarDiffusionFluxes_N_Impl (const amrex::Box &bx, const amrex::Box &domain, const ScalarDiffusionFieldViews &field, const ScalarDiffusionFluxPolicy &policy, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &u, const amrex::Array4< const amrex::Real > &v, const amrex::Array4< const amrex::Real > &mf_ux, const amrex::Array4< const amrex::Real > &mf_uy, const amrex::Array4< const amrex::Real > &mf_vx, const amrex::Array4< const amrex::Real > &mf_vy, const amrex::BCRec *bc_ptr, const int bc_comp)
void	BuildScalarDiffusionFluxes_N (const amrex::Box &bx, const amrex::Box &domain, const ScalarDiffusionFieldViews &field, const ScalarDiffusionFluxPolicy &policy, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &u, const amrex::Array4< const amrex::Real > &v, const amrex::Array4< const amrex::Real > &mf_ux, const amrex::Array4< const amrex::Real > &mf_uy, const amrex::Array4< const amrex::Real > &mf_vx, const amrex::Array4< const amrex::Real > &mf_vy, const amrex::BCRec *bc_ptr, const int bc_comp)
template<bool MultiplyMolecularByDensity, bool AddTurbulence>
void	BuildScalarDiffusionFluxes_S_Impl (const amrex::Box &bx, const amrex::Box &domain, const ScalarDiffusionFieldViews &field, const ScalarDiffusionFluxPolicy &policy, const amrex::Real dx_inv, const amrex::Real dy_inv, const amrex::Real *dz, const int klo, const int khi, const amrex::Array4< const amrex::Real > &mf_ux, const amrex::Array4< const amrex::Real > &mf_uy, const amrex::Array4< const amrex::Real > &mf_vx, const amrex::Array4< const amrex::Real > &mf_vy, const amrex::BCRec *bc_ptr, const int bc_comp)
void	BuildScalarDiffusionFluxes_S (const amrex::Box &bx, const amrex::Box &domain, const ScalarDiffusionFieldViews &field, const ScalarDiffusionFluxPolicy &policy, const amrex::Real dx_inv, const amrex::Real dy_inv, const amrex::Real *dz, const int klo, const int khi, const amrex::Array4< const amrex::Real > &mf_ux, const amrex::Array4< const amrex::Real > &mf_uy, const amrex::Array4< const amrex::Real > &mf_vx, const amrex::Array4< const amrex::Real > &mf_vy, const amrex::BCRec *bc_ptr, const int bc_comp)
template<bool MultiplyMolecularByDensity, bool AddTurbulence>
void	BuildScalarDiffusionFluxes_T_Impl (const amrex::Box &bx, const amrex::Box &domain, const ScalarDiffusionFieldViews &field, const ScalarDiffusionFluxPolicy &policy, const bool rotate, const amrex::Array4< const amrex::Real > &z_nd, const amrex::Array4< const amrex::Real > &z_cc, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &mf_ux, const amrex::Array4< const amrex::Real > &mf_vy, const amrex::BCRec *bc_ptr, const int bc_comp)
void	BuildScalarDiffusionFluxes_T (const amrex::Box &bx, const amrex::Box &domain, const ScalarDiffusionFieldViews &field, const ScalarDiffusionFluxPolicy &policy, const bool rotate, const amrex::Array4< const amrex::Real > &z_nd, const amrex::Array4< const amrex::Real > &z_cc, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &mf_ux, const amrex::Array4< const amrex::Real > &mf_vy, const amrex::BCRec *bc_ptr, const int bc_comp)
void	BuildScalarDiffusionMappedTransfers_N (const amrex::Box &bx, const amrex::Array4< const amrex::Real > &raw_x, const int raw_x_comp, const amrex::Array4< const amrex::Real > &raw_y, const int raw_y_comp, const amrex::Array4< const amrex::Real > &raw_z, const int raw_z_comp, const amrex::Array4< amrex::Real > &mapped_x, const int mapped_x_comp, const amrex::Array4< amrex::Real > &mapped_y, const int mapped_y_comp, const amrex::Array4< amrex::Real > &mapped_z, const int mapped_z_comp, const amrex::Array4< const amrex::Real > &mf_mx, const amrex::Array4< const amrex::Real > &mf_my)
void	BuildScalarDiffusionMappedTransfers_S (const amrex::Box &bx, const amrex::Array4< const amrex::Real > &raw_x, const int raw_x_comp, const amrex::Array4< const amrex::Real > &raw_y, const int raw_y_comp, const amrex::Array4< const amrex::Real > &raw_z, const int raw_z_comp, const amrex::Array4< amrex::Real > &mapped_x, const int mapped_x_comp, const amrex::Array4< amrex::Real > &mapped_y, const int mapped_y_comp, const amrex::Array4< amrex::Real > &mapped_z, const int mapped_z_comp, const amrex::Array4< const amrex::Real > &ax, const amrex::Array4< const amrex::Real > &ay, const amrex::Array4< const amrex::Real > &mf_mx, const amrex::Array4< const amrex::Real > &mf_my)
void	BuildScalarDiffusionMappedTransfers_T (const amrex::Box &bx, const amrex::Box &domain, const amrex::Array4< const amrex::Real > &raw_x, const amrex::Array4< const amrex::Real > &raw_y, const amrex::Array4< const amrex::Real > &raw_z, const int raw_flux_comp, const amrex::Array4< amrex::Real > &mapped_x, const int mapped_x_comp, const amrex::Array4< amrex::Real > &mapped_y, const int mapped_y_comp, const amrex::Array4< amrex::Real > &mapped_z, const int mapped_z_comp, const amrex::Array4< const amrex::Real > &z_nd, const amrex::Array4< const amrex::Real > &ax, const amrex::Array4< const amrex::Real > &ay, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &mf_mx, const amrex::Array4< const amrex::Real > &mf_uy, const amrex::Array4< const amrex::Real > &mf_my, const amrex::Array4< const amrex::Real > &mf_vx, const bool suppress_lower_surface_z_transfer)
void	ApplyScalarMappedFluxDivergence (const amrex::Box &bx, const amrex::Array4< const amrex::Real > &mapped_x, const int mapped_x_comp, const amrex::Array4< const amrex::Real > &mapped_y, const int mapped_y_comp, const amrex::Array4< const amrex::Real > &mapped_z, const int mapped_z_comp, const amrex::Array4< amrex::Real > &rhs, const int rhs_comp, const amrex::Array4< const amrex::Real > &detJ, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &mf_mx, const amrex::Array4< const amrex::Real > &mf_my)
void	ApplyScalarDiffusionMappedDivergence_T (const amrex::Box &bx, const amrex::Box &domain, const ScalarDiffusionFieldViews &field, const amrex::Array4< const amrex::Real > &z_nd, const amrex::Array4< const amrex::Real > &ax, const amrex::Array4< const amrex::Real > &ay, const amrex::Array4< const amrex::Real > &detJ, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &mf_mx, const amrex::Array4< const amrex::Real > &mf_uy, const amrex::Array4< const amrex::Real > &mf_my, const amrex::Array4< const amrex::Real > &mf_vx, const bool suppress_lower_surface_z_flux)
void	ApplyScalarDiffusionFluxDivergence_N (const amrex::Box &bx, const amrex::Array4< const amrex::Real > &xflux, const amrex::Array4< const amrex::Real > &yflux, const amrex::Array4< const amrex::Real > &raw_Fz, const int flux_comp, const amrex::Array4< amrex::Real > &rhs, const int rhs_comp, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &mf_mx, const amrex::Array4< const amrex::Real > &mf_my)
void	ApplyScalarDiffusionFluxDivergence_S (const amrex::Box &bx, const amrex::Array4< const amrex::Real > &xflux, const amrex::Array4< const amrex::Real > &yflux, const amrex::Array4< const amrex::Real > &raw_Fz, const int flux_comp, const amrex::Array4< amrex::Real > &rhs, const int rhs_comp, const amrex::Array4< const amrex::Real > &mf_mx, const amrex::Array4< const amrex::Real > &mf_my, const amrex::Real dx_inv, const amrex::Real dy_inv, const amrex::Real *dz)
template<bool MultiplyMolecularByDensity, bool AddTurbulence>
void	SubtractTerrainMetricImplicitPart_T_Impl (const amrex::Box &zbx, const amrex::Box &domain, const ScalarDiffusionFieldViews &field, const ScalarDiffusionCoefficients &coefficients, const amrex::Array4< const amrex::Real > &z_nd, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &mf_mx, const amrex::Array4< const amrex::Real > &mf_my, const amrex::Real implicit_fac)
void	SubtractTerrainMetricImplicitPart_T (const amrex::Box &bx, const amrex::Box &domain, const ScalarDiffusionFieldViews &field, const ScalarDiffusionFluxPolicy &policy, const amrex::Array4< const amrex::Real > &z_nd, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &mf_mx, const amrex::Array4< const amrex::Real > &mf_my, const amrex::Real implicit_fac)
void	ScaleScalarDiffusionVerticalFlux (const amrex::Box &zbx, const amrex::Array4< amrex::Real > &raw_Fz, const int flux_comp, const amrex::Real explicit_fac)

Function Documentation

ApplyScalarDiffusionFluxDivergence_N()

void ApplyScalarDiffusionFluxDivergence_N ( const amrex::Box &  bx, const amrex::Array4< const amrex::Real > &  xflux, const amrex::Array4< const amrex::Real > &  yflux, const amrex::Array4< const amrex::Real > &  raw_Fz, const int  flux_comp, const amrex::Array4< amrex::Real > &  rhs, const int  rhs_comp, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &  cellSizeInv, const amrex::Array4< const amrex::Real > &  mf_mx, const amrex::Array4< const amrex::Real > &  mf_my  )

inline

Apply the conservative N-grid divergence to a separately selected RHS component. Face scratch is component-selected independently in each view.

{
  const amrex::Real dx_inv = cellSizeInv[0];
  const amrex::Real dy_inv = cellSizeInv[1];
  const amrex::Real dz_inv = cellSizeInv[2];
  amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
    const amrex::Real mfsq = mf_mx(i, j, 0) * mf_my(i, j, 0);
    rhs(i, j, k, rhs_comp) -=
      (xflux(i + 1, j, k, flux_comp) - xflux(i, j, k, flux_comp)) * dx_inv *
        mfsq +
      (yflux(i, j + 1, k, flux_comp) - yflux(i, j, k, flux_comp)) * dy_inv *
        mfsq +
      (raw_Fz(i, j, k + 1, flux_comp) - raw_Fz(i, j, k, flux_comp)) * dz_inv;
  });
}

Referenced by DiffusionSrcForState_N().

Here is the call graph for this function:

Here is the caller graph for this function:

ApplyScalarDiffusionFluxDivergence_S()

void ApplyScalarDiffusionFluxDivergence_S ( const amrex::Box &  bx, const amrex::Array4< const amrex::Real > &  xflux, const amrex::Array4< const amrex::Real > &  yflux, const amrex::Array4< const amrex::Real > &  raw_Fz, const int  flux_comp, const amrex::Array4< amrex::Real > &  rhs, const int  rhs_comp, const amrex::Array4< const amrex::Real > &  mf_mx, const amrex::Array4< const amrex::Real > &  mf_my, const amrex::Real  dx_inv, const amrex::Real  dy_inv, const amrex::Real *  dz  )

inline

S-grid conservative divergence; vertical cell width is the receiver's dz[k].

{
  amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
    const amrex::Real mfsq = mf_mx(i, j, 0) * mf_my(i, j, 0);
    const amrex::Real dzk_inv = amrex::Real(1.0) / dz[k];
    const amrex::Real stateContrib =
      (xflux(i + 1, j, k, flux_comp) - xflux(i, j, k, flux_comp)) * dx_inv *
        mfsq +
      (yflux(i, j + 1, k, flux_comp) - yflux(i, j, k, flux_comp)) * dy_inv *
        mfsq +
      (raw_Fz(i, j, k + 1, flux_comp) - raw_Fz(i, j, k, flux_comp)) * dzk_inv;
    rhs(i, j, k, rhs_comp) -= stateContrib;
  });
}

Referenced by DiffusionSrcForState_S().

Here is the call graph for this function:

Here is the caller graph for this function:

ApplyScalarDiffusionMappedDivergence_T()

void ApplyScalarDiffusionMappedDivergence_T ( const amrex::Box &  bx, const amrex::Box &  domain, const ScalarDiffusionFieldViews &  field, const amrex::Array4< const amrex::Real > &  z_nd, const amrex::Array4< const amrex::Real > &  ax, const amrex::Array4< const amrex::Real > &  ay, const amrex::Array4< const amrex::Real > &  detJ, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &  cellSizeInv, const amrex::Array4< const amrex::Real > &  mf_mx, const amrex::Array4< const amrex::Real > &  mf_uy, const amrex::Array4< const amrex::Real > &  mf_my, const amrex::Array4< const amrex::Real > &  mf_vx, const bool  suppress_lower_surface_z_flux  )

inline

Apply native terrain diffusion in one kernel, mapping each face inline.

{
  using namespace amrex;
  const auto dom_lo = lbound(domain);
  const Real dx_inv = cellSizeInv[0];
  const Real dy_inv = cellSizeInv[1];
  const Real dz_inv = cellSizeInv[2];
  const Array4<const Real> xflux(field.xflux);
  const Array4<const Real> yflux(field.yflux);
  const Array4<const Real> zflux(field.zflux);
  const auto rhs = field.rhs;
  const int flux_comp = field.flux_comp;
  const int rhs_comp = field.rhs_comp;
 
  ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
    const Real mapped_Fx_lo = TerrainDiffusionMappedTx(
      xflux(i, j, k, flux_comp), ax(i, j, k), mf_uy(i, j, 0));
    const Real mapped_Fx_hi = TerrainDiffusionMappedTx(
      xflux(i + 1, j, k, flux_comp), ax(i + 1, j, k), mf_uy(i + 1, j, 0));
    const Real mapped_Fy_lo = TerrainDiffusionMappedTy(
      yflux(i, j, k, flux_comp), ay(i, j, k), mf_vx(i, j, 0));
    const Real mapped_Fy_hi = TerrainDiffusionMappedTy(
      yflux(i, j + 1, k, flux_comp), ay(i, j + 1, k), mf_vx(i, j + 1, 0));
    const Real Gz_lo = (suppress_lower_surface_z_flux && k == dom_lo.z)
                         ? Real(0.0)
                         : TerrainScalarDiffusionGzAtKFace(
                             i, j, k, domain, xflux, yflux, zflux, flux_comp,
                             z_nd, cellSizeInv, mf_mx, mf_my);
    const Real Gz_hi = TerrainScalarDiffusionGzAtKFace(
      i, j, k + 1, domain, xflux, yflux, zflux, flux_comp, z_nd, cellSizeInv,
      mf_mx, mf_my);
    const Real mfsq = mf_mx(i, j, 0) * mf_my(i, j, 0);
    const Real stateContrib = TerrainDiffusionDivergence_T(
      mapped_Fx_hi, mapped_Fx_lo, mapped_Fy_hi, mapped_Fy_lo, Gz_hi, Gz_lo,
      dx_inv, dy_inv, dz_inv, mfsq, detJ(i, j, k));
    rhs(i, j, k, rhs_comp) -= stateContrib;
  });
}

Referenced by DiffusionSrcForState_T().

Here is the call graph for this function:

Here is the caller graph for this function:

ApplyScalarMappedFluxDivergence()

void ApplyScalarMappedFluxDivergence ( const amrex::Box &  bx, const amrex::Array4< const amrex::Real > &  mapped_x, const int  mapped_x_comp, const amrex::Array4< const amrex::Real > &  mapped_y, const int  mapped_y_comp, const amrex::Array4< const amrex::Real > &  mapped_z, const int  mapped_z_comp, const amrex::Array4< amrex::Real > &  rhs, const int  rhs_comp, const amrex::Array4< const amrex::Real > &  detJ, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &  cellSizeInv, const amrex::Array4< const amrex::Real > &  mf_mx, const amrex::Array4< const amrex::Real > &  mf_my  )

inline

Add the divergence of caller-supplied mapped transfers to rhs. Inputs are only the accepted face arrays and the mapped cell geometry; in particular, this operation cannot reconstruct a raw or terrain diffusion candidate. The supplied face arrays and mapped-geometry convention are the same as ApplyScalarAdvectionFluxDivergence. This diffusion operation differs in RHS semantics: it accumulates -divergence into an existing RHS and leaves cells with nonpositive detJ unchanged, whereas advection assigns its output and zeroes such cells.

{
  const amrex::Real dx_inv = cellSizeInv[0];
  const amrex::Real dy_inv = cellSizeInv[1];
  const amrex::Real dzeta_inv = cellSizeInv[2];
  amrex::ParallelFor(
    bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
      if (detJ(i, j, k) > amrex::Real(0.0)) {
        rhs(i, j, k, rhs_comp) -= ScalarDiffusionMappedFluxDivergence(
          mapped_x(i + 1, j, k, mapped_x_comp),
          mapped_x(i, j, k, mapped_x_comp),
          mapped_y(i, j + 1, k, mapped_y_comp),
          mapped_y(i, j, k, mapped_y_comp),
          mapped_z(i, j, k + 1, mapped_z_comp),
          mapped_z(i, j, k, mapped_z_comp), dx_inv, dy_inv, dzeta_inv,
          mf_mx(i, j, 0), mf_my(i, j, 0), detJ(i, j, k));
      }
    });
}

Here is the call graph for this function:

BuildScalarDiffusionFluxes_N()

void BuildScalarDiffusionFluxes_N ( const amrex::Box &  bx, const amrex::Box &  domain, const ScalarDiffusionFieldViews &  field, const ScalarDiffusionFluxPolicy &  policy, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &  cellSizeInv, const amrex::Array4< const amrex::Real > &  u, const amrex::Array4< const amrex::Real > &  v, const amrex::Array4< const amrex::Real > &  mf_ux, const amrex::Array4< const amrex::Real > &  mf_uy, const amrex::Array4< const amrex::Real > &  mf_vx, const amrex::Array4< const amrex::Real > &  mf_vy, const amrex::BCRec *  bc_ptr, const int  bc_comp  )

inline

Build all native N-grid face fluxes for any intensive scalar field.

{
  const auto mode = policy.coefficient_mode;
  if (mode.molecular_is_kinematic) {
    if (mode.add_turbulence) {
      BuildScalarDiffusionFluxes_N_Impl<true, true>(
        bx, domain, field, policy, cellSizeInv, u, v, mf_ux, mf_uy, mf_vx,
        mf_vy, bc_ptr, bc_comp);
    } else {
      BuildScalarDiffusionFluxes_N_Impl<true, false>(
        bx, domain, field, policy, cellSizeInv, u, v, mf_ux, mf_uy, mf_vx,
        mf_vy, bc_ptr, bc_comp);
    }
  } else if (mode.add_turbulence) {
    BuildScalarDiffusionFluxes_N_Impl<false, true>(
      bx, domain, field, policy, cellSizeInv, u, v, mf_ux, mf_uy, mf_vx,
      mf_vy, bc_ptr, bc_comp);
  } else {
    BuildScalarDiffusionFluxes_N_Impl<false, false>(
      bx, domain, field, policy, cellSizeInv, u, v, mf_ux, mf_uy, mf_vx,
      mf_vy, bc_ptr, bc_comp);
  }
}

Referenced by DiffusionSrcForState_N().

Here is the caller graph for this function:

BuildScalarDiffusionFluxes_N_Impl()

template<bool MultiplyMolecularByDensity, bool AddTurbulence>

void BuildScalarDiffusionFluxes_N_Impl ( const amrex::Box &  bx, const amrex::Box &  domain, const ScalarDiffusionFieldViews &  field, const ScalarDiffusionFluxPolicy &  policy, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &  cellSizeInv, const amrex::Array4< const amrex::Real > &  u, const amrex::Array4< const amrex::Real > &  v, const amrex::Array4< const amrex::Real > &  mf_ux, const amrex::Array4< const amrex::Real > &  mf_uy, const amrex::Array4< const amrex::Real > &  mf_vx, const amrex::Array4< const amrex::Real > &  mf_vy, const amrex::BCRec *  bc_ptr, const int  bc_comp  )

inline

Component-explicit N-grid face construction. Boundary-condition selection, map factors, prescribed surface fluxes and raw-flux diagnostics are explicit inputs; this operation does not inspect native conserved-state identity.

{
  using namespace amrex;
  const auto dom_lo = lbound(domain);
  const auto dom_hi = ubound(domain);
  const Box xbx = surroundingNodes(bx, 0);
  const Box ybx = surroundingNodes(bx, 1);
  const Box zbx = surroundingNodes(bx, 2);
  const Real dx_inv = cellSizeInv[0];
  const Real dy_inv = cellSizeInv[1];
  const Real dz_inv = cellSizeInv[2];
  const auto coefficients = policy.coefficients;
  const auto prescribed = policy.prescribed;
  const auto diagnostic = policy.diagnostic;
  const auto surface = policy.surface;
  const int scalar_comp = field.scalar_comp;
  const int rho_comp = field.rho_comp;
  const int flux_comp = field.flux_comp;
  const auto scalar = field.scalar;
  const auto density = field.density;
  const auto mu_turb = field.mu_turb;
  const auto xflux = field.xflux;
  const auto yflux = field.yflux;
  const auto zflux = field.zflux;
 
  ParallelFor(xbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
    const Real coeff = ScalarDiffusionFaceCoefficient<
      MultiplyMolecularByDensity, AddTurbulence>(
      density, rho_comp, mu_turb, coefficients, i, j, k, 1, 0, 0,
      coefficients.eddy_h_comp);
    const bool ext_lo =
      ((bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir) ||
       (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_prim) ||
       (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_upwind &&
        u(dom_lo.x, j, k) >= Real(0.0))) &&
      i == dom_lo.x;
    const bool ext_hi =
      ((bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir) ||
       (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_prim) ||
       (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_upwind &&
        u(dom_hi.x + 1, j, k) <= Real(0.0))) &&
      i == dom_hi.x + 1;
    const bool on_surface =
      (surface.xlo && i == dom_lo.x) ||
      (surface.xhi && i == dom_hi.x + 1);
 
    if (ext_lo) {
      xflux(i, j, k, flux_comp) = ScalarDiffusionFlux_N<0>(
        scalar, scalar_comp, i, j, k, coeff, dx_inv,
        mf_ux(i, j, 0), mf_uy(i, j, 0), true, false);
    } else if (ext_hi) {
      xflux(i, j, k, flux_comp) = ScalarDiffusionFlux_N<0>(
        scalar, scalar_comp, i, j, k, coeff, dx_inv,
        mf_ux(i, j, 0), mf_uy(i, j, 0), false, true);
    } else if (on_surface && prescribed.use_x_on_side) {
      xflux(i, j, k, flux_comp) = prescribed.x(i, j, k, prescribed.x_comp);
    } else {
      xflux(i, j, k, flux_comp) = ScalarDiffusionFlux_N<0>(
        scalar, scalar_comp, i, j, k, coeff, dx_inv,
        mf_ux(i, j, 0), mf_uy(i, j, 0), false, false);
    }
  });
 
  ParallelFor(ybx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
    const Real coeff = ScalarDiffusionFaceCoefficient<
      MultiplyMolecularByDensity, AddTurbulence>(
      density, rho_comp, mu_turb, coefficients, i, j, k, 0, 1, 0,
      coefficients.eddy_h_comp);
    const bool ext_lo =
      ((bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir) ||
       (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_prim) ||
       (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_upwind &&
        v(i, dom_lo.y, k) >= Real(0.0))) &&
      j == dom_lo.y;
    const bool ext_hi =
      ((bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir) ||
       (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_prim) ||
       (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_upwind &&
        v(i, dom_hi.y + 1, k) <= Real(0.0))) &&
      j == dom_hi.y + 1;
    const bool on_surface =
      (surface.ylo && j == dom_lo.y) ||
      (surface.yhi && j == dom_hi.y + 1);
 
    if (ext_lo) {
      yflux(i, j, k, flux_comp) = ScalarDiffusionFlux_N<1>(
        scalar, scalar_comp, i, j, k, coeff, dy_inv,
        mf_vy(i, j, 0), mf_vx(i, j, 0), true, false);
    } else if (ext_hi) {
      yflux(i, j, k, flux_comp) = ScalarDiffusionFlux_N<1>(
        scalar, scalar_comp, i, j, k, coeff, dy_inv,
        mf_vy(i, j, 0), mf_vx(i, j, 0), false, true);
    } else if (on_surface && prescribed.use_y_on_side) {
      yflux(i, j, k, flux_comp) = prescribed.y(i, j, k, prescribed.y_comp);
    } else {
      yflux(i, j, k, flux_comp) = ScalarDiffusionFlux_N<1>(
        scalar, scalar_comp, i, j, k, coeff, dy_inv,
        mf_vy(i, j, 0), mf_vx(i, j, 0), false, false);
    }
  });
 
  ParallelFor(zbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
    const Real coeff = ScalarDiffusionFaceCoefficient<
      MultiplyMolecularByDensity, AddTurbulence>(
      density, rho_comp, mu_turb, coefficients, i, j, k, 0, 0, 1,
      coefficients.eddy_v_comp);
    const bool ext_lo =
      ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) ||
       (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim)) &&
      k == dom_lo.z;
    const bool ext_hi =
      ((bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir) ||
       (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim)) &&
      k == dom_hi.z + 1;
    const bool on_surface =
      (surface.zlo && k == dom_lo.z) ||
      (surface.zhi && k == dom_hi.z + 1);
 
    if (ext_lo) {
      zflux(i, j, k, flux_comp) = ScalarDiffusionFlux_N<2>(
        scalar, scalar_comp, i, j, k, coeff, dz_inv, Real(1.0), Real(1.0),
        true, false);
    } else if (ext_hi) {
      zflux(i, j, k, flux_comp) = ScalarDiffusionFlux_N<2>(
        scalar, scalar_comp, i, j, k, coeff, dz_inv, Real(1.0), Real(1.0),
        false, true);
    } else if (on_surface && prescribed.use_z_on_side) {
      zflux(i, j, k, flux_comp) = prescribed.z(i, j, k, prescribed.z_comp);
    } else if (on_surface) {
      zflux(i, j, k, flux_comp) = Real(0.0);
    } else {
      zflux(i, j, k, flux_comp) = ScalarDiffusionFlux_N<2>(
        scalar, scalar_comp, i, j, k, coeff, dz_inv, Real(1.0), Real(1.0),
        false, false);
    }
 
    if (diagnostic.enabled && (diagnostic.write_on_surface || !on_surface)) {
      diagnostic.z(i, j, k, diagnostic.z_comp) = zflux(i, j, k, flux_comp);
    }
  });
}

Here is the call graph for this function:

BuildScalarDiffusionFluxes_S()

void BuildScalarDiffusionFluxes_S ( const amrex::Box &  bx, const amrex::Box &  domain, const ScalarDiffusionFieldViews &  field, const ScalarDiffusionFluxPolicy &  policy, const amrex::Real  dx_inv, const amrex::Real  dy_inv, const amrex::Real *  dz, const int  klo, const int  khi, const amrex::Array4< const amrex::Real > &  mf_ux, const amrex::Array4< const amrex::Real > &  mf_uy, const amrex::Array4< const amrex::Real > &  mf_vx, const amrex::Array4< const amrex::Real > &  mf_vy, const amrex::BCRec *  bc_ptr, const int  bc_comp  )

inline

Build and map-adjust all native S-grid face fluxes for an intensive field.

{
  const auto mode = policy.coefficient_mode;
  if (mode.molecular_is_kinematic) {
    if (mode.add_turbulence) {
      BuildScalarDiffusionFluxes_S_Impl<true, true>(
        bx, domain, field, policy, dx_inv, dy_inv, dz, klo, khi, mf_ux,
        mf_uy, mf_vx, mf_vy, bc_ptr, bc_comp);
    } else {
      BuildScalarDiffusionFluxes_S_Impl<true, false>(
        bx, domain, field, policy, dx_inv, dy_inv, dz, klo, khi, mf_ux,
        mf_uy, mf_vx, mf_vy, bc_ptr, bc_comp);
    }
  } else if (mode.add_turbulence) {
    BuildScalarDiffusionFluxes_S_Impl<false, true>(
      bx, domain, field, policy, dx_inv, dy_inv, dz, klo, khi, mf_ux,
      mf_uy, mf_vx, mf_vy, bc_ptr, bc_comp);
  } else {
    BuildScalarDiffusionFluxes_S_Impl<false, false>(
      bx, domain, field, policy, dx_inv, dy_inv, dz, klo, khi, mf_ux,
      mf_uy, mf_vx, mf_vy, bc_ptr, bc_comp);
  }
}

Referenced by DiffusionSrcForState_S().

Here is the caller graph for this function:

BuildScalarDiffusionFluxes_S_Impl()

template<bool MultiplyMolecularByDensity, bool AddTurbulence>

void BuildScalarDiffusionFluxes_S_Impl ( const amrex::Box &  bx, const amrex::Box &  domain, const ScalarDiffusionFieldViews &  field, const ScalarDiffusionFluxPolicy &  policy, const amrex::Real  dx_inv, const amrex::Real  dy_inv, const amrex::Real *  dz, const int  klo, const int  khi, const amrex::Array4< const amrex::Real > &  mf_ux, const amrex::Array4< const amrex::Real > &  mf_uy, const amrex::Array4< const amrex::Real > &  mf_vx, const amrex::Array4< const amrex::Real > &  mf_vy, const amrex::BCRec *  bc_ptr, const int  bc_comp  )

inline

{
  using namespace amrex;
  const auto dom_lo = lbound(domain);
  const auto dom_hi = ubound(domain);
  const Box xbx = surroundingNodes(bx, 0);
  const Box ybx = surroundingNodes(bx, 1);
  const Box zbx = surroundingNodes(bx, 2);
  const auto coefficients = policy.coefficients;
  const auto prescribed = policy.prescribed;
  const auto diagnostic = policy.diagnostic;
  const auto surface = policy.surface;
  const int scalar_comp = field.scalar_comp;
  const int rho_comp = field.rho_comp;
  const int flux_comp = field.flux_comp;
  const auto scalar = field.scalar;
  const auto density = field.density;
  const auto mu_turb = field.mu_turb;
  const auto xflux = field.xflux;
  const auto yflux = field.yflux;
  const auto zflux = field.zflux;
 
  ParallelFor(xbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
    const Real coeff = ScalarDiffusionFaceCoefficient<
      MultiplyMolecularByDensity, AddTurbulence>(
      density, rho_comp, mu_turb, coefficients, i, j, k, 1, 0, 0,
      coefficients.eddy_h_comp);
    const bool on_surface =
      (surface.xlo && i == dom_lo.x) ||
      (surface.xhi && i == dom_hi.x + 1);
    if (on_surface && prescribed.use_x_on_side) {
      xflux(i, j, k, flux_comp) = prescribed.x(i, j, k, prescribed.x_comp);
    } else {
      xflux(i, j, k, flux_comp) = ScalarDiffusionFlux_S_Horizontal<0>(
        scalar, scalar_comp, i, j, k, coeff, dx_inv, mf_ux(i, j, 0));
    }
  });
 
  ParallelFor(ybx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
    const Real coeff = ScalarDiffusionFaceCoefficient<
      MultiplyMolecularByDensity, AddTurbulence>(
      density, rho_comp, mu_turb, coefficients, i, j, k, 0, 1, 0,
      coefficients.eddy_h_comp);
    const bool on_surface =
      (surface.ylo && j == dom_lo.y) ||
      (surface.yhi && j == dom_hi.y + 1);
    if (on_surface && prescribed.use_y_on_side) {
      yflux(i, j, k, flux_comp) = prescribed.y(i, j, k, prescribed.y_comp);
    } else {
      yflux(i, j, k, flux_comp) = ScalarDiffusionFlux_S_Horizontal<1>(
        scalar, scalar_comp, i, j, k, coeff, dy_inv, mf_vy(i, j, 0));
    }
  });
 
  ParallelFor(zbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
    const Real coeff = ScalarDiffusionFaceCoefficient<
      MultiplyMolecularByDensity, AddTurbulence>(
      density, rho_comp, mu_turb, coefficients, i, j, k, 0, 0, 1,
      coefficients.eddy_v_comp);
    const bool ext_lo =
      ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) ||
       (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim)) &&
      k == dom_lo.z;
    const bool ext_hi =
      ((bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir) ||
       (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim)) &&
      k == dom_hi.z + 1;
    const bool on_surface =
      (surface.zlo && k == dom_lo.z) ||
      (surface.zhi && k == dom_hi.z + 1);
    const Real grad = StretchedScalarGradient(
      scalar, scalar_comp, i, j, k, dz, klo, khi, ext_lo, ext_hi);
    if (on_surface && prescribed.use_z_on_side) {
      zflux(i, j, k, flux_comp) = prescribed.z(i, j, k, prescribed.z_comp);
    } else if (on_surface) {
      zflux(i, j, k, flux_comp) = Real(0.0);
    } else {
      zflux(i, j, k, flux_comp) = -coeff * grad;
    }
    if (diagnostic.enabled && (diagnostic.write_on_surface || !on_surface)) {
      diagnostic.z(i, j, k, diagnostic.z_comp) = zflux(i, j, k, flux_comp);
    }
  });
 
  // S-grid raw horizontal fluxes are adjusted in the same separate kernels as
  // the native implementation before scaling or applying divergence.
  ParallelFor(xbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
    xflux(i, j, k, flux_comp) /= mf_uy(i, j, 0);
  });
  ParallelFor(ybx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
    yflux(i, j, k, flux_comp) /= mf_vx(i, j, 0);
  });
}

Here is the call graph for this function:

BuildScalarDiffusionFluxes_T()

void BuildScalarDiffusionFluxes_T ( const amrex::Box &  bx, const amrex::Box &  domain, const ScalarDiffusionFieldViews &  field, const ScalarDiffusionFluxPolicy &  policy, const bool  rotate, const amrex::Array4< const amrex::Real > &  z_nd, const amrex::Array4< const amrex::Real > &  z_cc, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &  cellSizeInv, const amrex::Array4< const amrex::Real > &  mf_ux, const amrex::Array4< const amrex::Real > &  mf_vy, const amrex::BCRec *  bc_ptr, const int  bc_comp  )

inline

Build all native terrain-grid face fluxes, including the grown horizontal face regions needed by the mapped z-face transfer.

{
  const auto mode = policy.coefficient_mode;
  if (mode.molecular_is_kinematic) {
    if (mode.add_turbulence) {
      BuildScalarDiffusionFluxes_T_Impl<true, true>(
        bx, domain, field, policy, rotate, z_nd, z_cc, cellSizeInv, mf_ux,
        mf_vy, bc_ptr, bc_comp);
    } else {
      BuildScalarDiffusionFluxes_T_Impl<true, false>(
        bx, domain, field, policy, rotate, z_nd, z_cc, cellSizeInv, mf_ux,
        mf_vy, bc_ptr, bc_comp);
    }
  } else if (mode.add_turbulence) {
    BuildScalarDiffusionFluxes_T_Impl<false, true>(
      bx, domain, field, policy, rotate, z_nd, z_cc, cellSizeInv, mf_ux,
      mf_vy, bc_ptr, bc_comp);
  } else {
    BuildScalarDiffusionFluxes_T_Impl<false, false>(
      bx, domain, field, policy, rotate, z_nd, z_cc, cellSizeInv, mf_ux,
      mf_vy, bc_ptr, bc_comp);
  }
}

Referenced by DiffusionSrcForState_T().

Here is the caller graph for this function:

BuildScalarDiffusionFluxes_T_Impl()

template<bool MultiplyMolecularByDensity, bool AddTurbulence>

void BuildScalarDiffusionFluxes_T_Impl ( const amrex::Box &  bx, const amrex::Box &  domain, const ScalarDiffusionFieldViews &  field, const ScalarDiffusionFluxPolicy &  policy, const bool  rotate, const amrex::Array4< const amrex::Real > &  z_nd, const amrex::Array4< const amrex::Real > &  z_cc, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &  cellSizeInv, const amrex::Array4< const amrex::Real > &  mf_ux, const amrex::Array4< const amrex::Real > &  mf_vy, const amrex::BCRec *  bc_ptr, const int  bc_comp  )

inline

{
  using namespace amrex;
  const auto dom_lo = lbound(domain);
  const auto dom_hi = ubound(domain);
  Box xbx = surroundingNodes(bx, 0);
  Box ybx = surroundingNodes(bx, 1);
  const Box zbx = surroundingNodes(bx, 2);
  if (xbx.smallEnd(2) != dom_lo.z) xbx.growLo(2, 1);
  if (ybx.smallEnd(2) != dom_lo.z) ybx.growLo(2, 1);
  if (xbx.bigEnd(2) != dom_hi.z) xbx.growHi(2, 1);
  if (ybx.bigEnd(2) != dom_hi.z) ybx.growHi(2, 1);
 
  const auto coefficients = policy.coefficients;
  const auto prescribed = policy.prescribed;
  const auto diagnostic = policy.diagnostic;
  const auto surface = policy.surface;
  const int scalar_comp = field.scalar_comp;
  const int rho_comp = field.rho_comp;
  const int flux_comp = field.flux_comp;
  const auto scalar = field.scalar;
  const auto density = field.density;
  const auto mu_turb = field.mu_turb;
  const auto xflux = field.xflux;
  const auto yflux = field.yflux;
  const auto zflux = field.zflux;
 
  ParallelFor(xbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
    const Real coeff = ScalarDiffusionFaceCoefficient<
      MultiplyMolecularByDensity, AddTurbulence>(
      density, rho_comp, mu_turb, coefficients, i, j, k, 1, 0, 0,
      coefficients.eddy_h_comp);
    const bool on_side =
      (surface.xlo && i == dom_lo.x) ||
      (surface.xhi && i == dom_hi.x + 1);
    const bool on_bottom_z =
      surface.zlo && rotate && k == dom_lo.z;
 
    if (on_side && prescribed.use_x_on_side) {
      xflux(i, j, k, flux_comp) = prescribed.x(i, j, k, prescribed.x_comp);
    } else if (on_bottom_z && prescribed.use_x_on_bottom_z) {
      xflux(i, j, k, flux_comp) = prescribed.x(i, j, 0, prescribed.x_comp);
    } else {
      xflux(i, j, k, flux_comp) = TerrainScalarDiffusionFluxAtIFace(
        scalar, scalar_comp, coeff, i, j, k, z_nd, z_cc, cellSizeInv, mf_ux);
    }
  });
 
  ParallelFor(ybx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
    const Real coeff = ScalarDiffusionFaceCoefficient<
      MultiplyMolecularByDensity, AddTurbulence>(
      density, rho_comp, mu_turb, coefficients, i, j, k, 0, 1, 0,
      coefficients.eddy_h_comp);
    const bool on_side =
      (surface.ylo && j == dom_lo.y) ||
      (surface.yhi && j == dom_hi.y + 1);
    const bool on_bottom_z =
      surface.zlo && rotate && k == dom_lo.z;
 
    if (on_side && prescribed.use_y_on_side) {
      yflux(i, j, k, flux_comp) = prescribed.y(i, j, k, prescribed.y_comp);
    } else if (on_bottom_z && prescribed.use_y_on_bottom_z) {
      yflux(i, j, k, flux_comp) = prescribed.y(i, j, 0, prescribed.y_comp);
    } else {
      yflux(i, j, k, flux_comp) = TerrainScalarDiffusionFluxAtJFace(
        scalar, scalar_comp, coeff, i, j, k, z_nd, z_cc, cellSizeInv, mf_vy);
    }
  });
 
  ParallelFor(zbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
    const Real coeff = ScalarDiffusionFaceCoefficient<
      MultiplyMolecularByDensity, AddTurbulence>(
      density, rho_comp, mu_turb, coefficients, i, j, k, 0, 0, 1,
      coefficients.eddy_v_comp);
    const bool ext_lo =
      ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) ||
       (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim)) &&
      k == dom_lo.z;
    const bool ext_hi =
      ((bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir) ||
       (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim)) &&
      k == dom_hi.z + 1;
    const bool on_surface =
      (surface.zlo && k == dom_lo.z) ||
      (surface.zhi && k == dom_hi.z + 1);
 
    const Real raw_flux = TerrainScalarDiffusionFluxAtKFace(
      scalar, scalar_comp, coeff, i, j, k, z_nd, cellSizeInv, ext_lo, ext_hi);
 
    if (on_surface && prescribed.use_z_on_side) {
      zflux(i, j, k, flux_comp) = prescribed.z(i, j, k, prescribed.z_comp);
    } else if (on_surface) {
      zflux(i, j, k, flux_comp) = Real(0.0);
    } else {
      zflux(i, j, k, flux_comp) = raw_flux;
    }
    if (diagnostic.enabled && (diagnostic.write_on_surface || !on_surface)) {
      diagnostic.z(i, j, k, diagnostic.z_comp) = zflux(i, j, k, flux_comp);
    }
  });
}

Here is the call graph for this function:

BuildScalarDiffusionMappedTransfers_N()

void BuildScalarDiffusionMappedTransfers_N ( const amrex::Box &  bx, const amrex::Array4< const amrex::Real > &  raw_x, const int  raw_x_comp, const amrex::Array4< const amrex::Real > &  raw_y, const int  raw_y_comp, const amrex::Array4< const amrex::Real > &  raw_z, const int  raw_z_comp, const amrex::Array4< amrex::Real > &  mapped_x, const int  mapped_x_comp, const amrex::Array4< amrex::Real > &  mapped_y, const int  mapped_y_comp, const amrex::Array4< amrex::Real > &  mapped_z, const int  mapped_z_comp, const amrex::Array4< const amrex::Real > &  mf_mx, const amrex::Array4< const amrex::Real > &  mf_my  )

inline

Materialize the canonical N-grid face transfers. N-grid x/y builder outputs already carry their map-factor ratios; only the vertical flux is divided by the cell-centered horizontal map-factor product.

{
  const amrex::Box xbx = amrex::surroundingNodes(bx, 0);
  const amrex::Box ybx = amrex::surroundingNodes(bx, 1);
  const amrex::Box zbx = amrex::surroundingNodes(bx, 2);
  amrex::ParallelFor(
    xbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
      mapped_x(i, j, k, mapped_x_comp) =
        ScalarDiffusionMappedHorizontalTransfer(
          raw_x(i, j, k, raw_x_comp), amrex::Real(1.0));
    });
  amrex::ParallelFor(
    ybx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
      mapped_y(i, j, k, mapped_y_comp) =
        ScalarDiffusionMappedHorizontalTransfer(
          raw_y(i, j, k, raw_y_comp), amrex::Real(1.0));
    });
  amrex::ParallelFor(
    zbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
      mapped_z(i, j, k, mapped_z_comp) =
        ScalarDiffusionMappedVerticalTransfer(
          raw_z(i, j, k, raw_z_comp), mf_mx(i, j, 0), mf_my(i, j, 0));
    });
}

Here is the call graph for this function:

BuildScalarDiffusionMappedTransfers_S()

void BuildScalarDiffusionMappedTransfers_S ( const amrex::Box &  bx, const amrex::Array4< const amrex::Real > &  raw_x, const int  raw_x_comp, const amrex::Array4< const amrex::Real > &  raw_y, const int  raw_y_comp, const amrex::Array4< const amrex::Real > &  raw_z, const int  raw_z_comp, const amrex::Array4< amrex::Real > &  mapped_x, const int  mapped_x_comp, const amrex::Array4< amrex::Real > &  mapped_y, const int  mapped_y_comp, const amrex::Array4< amrex::Real > &  mapped_z, const int  mapped_z_comp, const amrex::Array4< const amrex::Real > &  ax, const amrex::Array4< const amrex::Real > &  ay, const amrex::Array4< const amrex::Real > &  mf_mx, const amrex::Array4< const amrex::Real > &  mf_my  )

inline

Materialize the canonical S-grid face transfers. ax/ay restore the side Jacobian measure omitted by the native S-grid divergence; detJ and the computational z spacing are supplied separately to the mapped divergence.

{
  const amrex::Box xbx = amrex::surroundingNodes(bx, 0);
  const amrex::Box ybx = amrex::surroundingNodes(bx, 1);
  const amrex::Box zbx = amrex::surroundingNodes(bx, 2);
  amrex::ParallelFor(
    xbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
      mapped_x(i, j, k, mapped_x_comp) =
        ScalarDiffusionMappedHorizontalTransfer(
          raw_x(i, j, k, raw_x_comp), ax(i, j, k));
    });
  amrex::ParallelFor(
    ybx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
      mapped_y(i, j, k, mapped_y_comp) =
        ScalarDiffusionMappedHorizontalTransfer(
          raw_y(i, j, k, raw_y_comp), ay(i, j, k));
    });
  amrex::ParallelFor(
    zbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
      mapped_z(i, j, k, mapped_z_comp) =
        ScalarDiffusionMappedVerticalTransfer(
          raw_z(i, j, k, raw_z_comp), mf_mx(i, j, 0), mf_my(i, j, 0));
    });
}

Here is the call graph for this function:

BuildScalarDiffusionMappedTransfers_T()

void BuildScalarDiffusionMappedTransfers_T ( const amrex::Box &  bx, const amrex::Box &  domain, const amrex::Array4< const amrex::Real > &  raw_x, const amrex::Array4< const amrex::Real > &  raw_y, const amrex::Array4< const amrex::Real > &  raw_z, const int  raw_flux_comp, const amrex::Array4< amrex::Real > &  mapped_x, const int  mapped_x_comp, const amrex::Array4< amrex::Real > &  mapped_y, const int  mapped_y_comp, const amrex::Array4< amrex::Real > &  mapped_z, const int  mapped_z_comp, const amrex::Array4< const amrex::Real > &  z_nd, const amrex::Array4< const amrex::Real > &  ax, const amrex::Array4< const amrex::Real > &  ay, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &  cellSizeInv, const amrex::Array4< const amrex::Real > &  mf_mx, const amrex::Array4< const amrex::Real > &  mf_uy, const amrex::Array4< const amrex::Real > &  mf_my, const amrex::Array4< const amrex::Real > &  mf_vx, const bool  suppress_lower_surface_z_transfer  )

inline

Materialize the canonical terrain-following transfers using ERF's existing face mappings and authoritative Gz construction. The lower suppression switch zeros the complete lower transfer, including terrain cross terms.

Corresponding directional inputs and outputs may alias in place, including the same component. The z transfer is launched first on AMReX's ordered execution stream so it reads raw x/y fluxes before the x/y mapping launches can overwrite them. Cross-direction aliasing is not supported.

For each interior z face mapped from bx, raw_x and raw_y must be valid at both kface and kface-1. Thus, if bx.smallEnd(2) is above the domain bottom, raw_x/raw_y must include bx.smallEnd(2)-1; if bx.bigEnd(2) is below the domain top, they must include bx.bigEnd(2)+1. No outside-domain layer is required at a physical bottom/top face because the Gz helper uses one-sided extrapolation. Raw candidates produced by BuildScalarDiffusionFluxes_T already provide this z coverage.

{
  const amrex::Box xbx = amrex::surroundingNodes(bx, 0);
  const amrex::Box ybx = amrex::surroundingNodes(bx, 1);
  const amrex::Box zbx = amrex::surroundingNodes(bx, 2);
  const amrex::Array4<const amrex::Real> xflux(raw_x);
  const amrex::Array4<const amrex::Real> yflux(raw_y);
  const amrex::Array4<const amrex::Real> zflux(raw_z);
  amrex::ParallelFor(
    zbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
      const bool suppress = suppress_lower_surface_z_transfer &&
                            k == domain.smallEnd(2);
      mapped_z(i, j, k, mapped_z_comp) =
        TerrainScalarDiffusionMappedTransferAtKFace(
          i, j, k, domain, xflux, yflux, zflux, raw_flux_comp, z_nd,
          cellSizeInv, mf_mx, mf_my, suppress);
    });
  amrex::ParallelFor(xbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
    mapped_x(i, j, k, mapped_x_comp) = TerrainDiffusionMappedTx(
      raw_x(i, j, k, raw_flux_comp), ax(i, j, k), mf_uy(i, j, 0));
  });
  amrex::ParallelFor(ybx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
    mapped_y(i, j, k, mapped_y_comp) = TerrainDiffusionMappedTy(
      raw_y(i, j, k, raw_flux_comp), ay(i, j, k), mf_vx(i, j, 0));
  });
}

Here is the call graph for this function:

ScalarDiffusionFaceCoefficient()

template<bool MultiplyMolecularByDensity, bool AddTurbulence>

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real ScalarDiffusionFaceCoefficient ( const amrex::Array4< const amrex::Real > &  rho, const int  rho_comp, const amrex::Array4< const amrex::Real > &  mu_turb, const ScalarDiffusionCoefficients &  coefficients, const int  i, const int  j, const int  k, const int  di, const int  dj, const int  dk, const int  eddy_comp  )

noexcept

Density-weighted face transport coefficient. The booleans are compile-time choices for ConstantAlpha/Constant/None and turbulent/no-turbulence modes; turbulent coefficients are already density weighted and are never multiplied by rho here.

rho_comp and the eddy component are independent of the scalar component.

{
  amrex::Real coeff = coefficients.molecular_coeff;
  if constexpr (MultiplyMolecularByDensity) {
    coeff *= amrex::Real(0.5) *
             (rho(i, j, k, rho_comp) + rho(i - di, j - dj, k - dk, rho_comp));
  }
  if constexpr (AddTurbulence) {
    coeff += amrex::Real(0.5) * (mu_turb(i, j, k, eddy_comp) +
                                 mu_turb(i - di, j - dj, k - dk, eddy_comp));
  }
  return coeff;
}

Referenced by BuildScalarDiffusionFluxes_N_Impl(), BuildScalarDiffusionFluxes_S_Impl(), BuildScalarDiffusionFluxes_T_Impl(), and SubtractTerrainMetricImplicitPart_T_Impl().

Here is the caller graph for this function:

ScalarDiffusionFlux_N()

template<int Direction>

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real ScalarDiffusionFlux_N ( const amrex::Array4< const amrex::Real > &  scalar, const int  scalar_comp, const int  i, const int  j, const int  k, const amrex::Real  coeff, const amrex::Real  inv_spacing, const amrex::Real  map_numerator, const amrex::Real  map_denominator, const bool  ext_dir_lo, const bool  ext_dir_hi  )

noexcept

Compute the existing N-grid signed face flux, F^D=-K grad(chi). scalar_comp selects the intensive scalar and flux_comp is applied by the caller to its independent face-flux view. Horizontal map-factor numerator and denominator remain separate to preserve the native arithmetic association. Physical Dirichlet flags select the legacy one-sided stencil; ext_dir_upwind selection stays in the adapter.

{
  static_assert(Direction >= 0 && Direction < 3);
  amrex::Real grad;
  if constexpr (Direction == 0) {
    if (ext_dir_lo) {
      grad = -(amrex::Real(8.0) / amrex::Real(3.0)) *
               scalar(i - 1, j, k, scalar_comp) +
             amrex::Real(3.0) * scalar(i, j, k, scalar_comp) -
             (amrex::Real(1.0) / amrex::Real(3.0)) *
               scalar(i + 1, j, k, scalar_comp);
    } else if (ext_dir_hi) {
      grad =
        (amrex::Real(8.0) / amrex::Real(3.0)) * scalar(i, j, k, scalar_comp) -
        amrex::Real(3.0) * scalar(i - 1, j, k, scalar_comp) +
        (amrex::Real(1.0) / amrex::Real(3.0)) *
          scalar(i - 2, j, k, scalar_comp);
    } else {
      grad = scalar(i, j, k, scalar_comp) - scalar(i - 1, j, k, scalar_comp);
    }
  } else if constexpr (Direction == 1) {
    if (ext_dir_lo) {
      grad = -(amrex::Real(8.0) / amrex::Real(3.0)) *
               scalar(i, j - 1, k, scalar_comp) +
             amrex::Real(3.0) * scalar(i, j, k, scalar_comp) -
             (amrex::Real(1.0) / amrex::Real(3.0)) *
               scalar(i, j + 1, k, scalar_comp);
    } else if (ext_dir_hi) {
      grad =
        (amrex::Real(8.0) / amrex::Real(3.0)) * scalar(i, j, k, scalar_comp) -
        amrex::Real(3.0) * scalar(i, j - 1, k, scalar_comp) +
        (amrex::Real(1.0) / amrex::Real(3.0)) *
          scalar(i, j - 2, k, scalar_comp);
    } else {
      grad = scalar(i, j, k, scalar_comp) - scalar(i, j - 1, k, scalar_comp);
    }
  } else {
    if (ext_dir_lo) {
      grad = -(amrex::Real(8.0) / amrex::Real(3.0)) *
               scalar(i, j, k - 1, scalar_comp) +
             amrex::Real(3.0) * scalar(i, j, k, scalar_comp) -
             (amrex::Real(1.0) / amrex::Real(3.0)) *
               scalar(i, j, k + 1, scalar_comp);
    } else if (ext_dir_hi) {
      grad =
        (amrex::Real(8.0) / amrex::Real(3.0)) * scalar(i, j, k, scalar_comp) -
        amrex::Real(3.0) * scalar(i, j, k - 1, scalar_comp) +
        (amrex::Real(1.0) / amrex::Real(3.0)) *
          scalar(i, j, k - 2, scalar_comp);
    } else {
      grad = scalar(i, j, k, scalar_comp) - scalar(i, j, k - 1, scalar_comp);
    }
  }
  if constexpr (Direction == 2) {
    return -coeff * grad * inv_spacing;
  } else {
    return -coeff * grad * inv_spacing * map_numerator / map_denominator;
  }
}

ScalarDiffusionFlux_S_Horizontal()

template<int Direction>

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real ScalarDiffusionFlux_S_Horizontal ( const amrex::Array4< const amrex::Real > &  scalar, const int  scalar_comp, const int  i, const int  j, const int  k, const amrex::Real  coeff, const amrex::Real  inv_spacing, const amrex::Real  face_map  )

noexcept

Existing S-grid horizontal raw flux before its separate map-factor adjustment kernel.

{
  static_assert(Direction == 0 || Direction == 1);
  const amrex::Real delta =
    Direction == 0
      ? scalar(i, j, k, scalar_comp) - scalar(i - 1, j, k, scalar_comp)
      : scalar(i, j, k, scalar_comp) - scalar(i, j - 1, k, scalar_comp);
  const amrex::Real gradient = inv_spacing * delta;
  return -coeff * face_map * gradient;
}

ScalarDiffusionFlux_Tx()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real ScalarDiffusionFlux_Tx ( const amrex::Real  coeff, const amrex::Real  face_map, const amrex::Real  grad_xi, const amrex::Real  h_xi, const amrex::Real  grad_z  )

noexcept

{
  return -coeff * face_map * (grad_xi - h_xi * grad_z);
}

Referenced by TerrainScalarDiffusionFluxAtIFace().

Here is the caller graph for this function:

ScalarDiffusionFlux_Ty()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real ScalarDiffusionFlux_Ty ( const amrex::Real  coeff, const amrex::Real  face_map, const amrex::Real  grad_eta, const amrex::Real  h_eta, const amrex::Real  grad_z  )

noexcept

{
  return -coeff * face_map * (grad_eta - h_eta * grad_z);
}

Referenced by TerrainScalarDiffusionFluxAtJFace().

Here is the caller graph for this function:

ScalarDiffusionMappedFluxDivergence()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real ScalarDiffusionMappedFluxDivergence ( const amrex::Real  fx_hi, const amrex::Real  fx_lo, const amrex::Real  fy_hi, const amrex::Real  fy_lo, const amrex::Real  fz_hi, const amrex::Real  fz_lo, const amrex::Real  dx_inv, const amrex::Real  dy_inv, const amrex::Real  dzeta_inv, const amrex::Real  mx, const amrex::Real  my, const amrex::Real  detJ_cell  )

noexcept

Return the positive geometry-weighted mapped flux divergence. For detJ_cell <= 0, return zero without dividing. The diffusion array wrapper subtracts this value from the existing RHS.

{
  if (detJ_cell <= amrex::Real(0.0))
    return amrex::Real(0.0);
  const amrex::Real map_over_jacobian = (mx * my) / detJ_cell;
  return map_over_jacobian *
         ((fx_hi - fx_lo) * dx_inv + (fy_hi - fy_lo) * dy_inv +
          (fz_hi - fz_lo) * dzeta_inv);
}

Referenced by ApplyScalarMappedFluxDivergence().

Here is the caller graph for this function:

ScalarDiffusionMappedHorizontalTransfer()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real ScalarDiffusionMappedHorizontalTransfer ( const amrex::Real  raw_flux, const amrex::Real  face_area  )

noexcept

Pointwise conversion to ERF's mapped conservative face-transfer convention. The transfer is an instantaneous rate and includes the supplied face-area measure; callers must not apply that measure again in the divergence.

{
  return raw_flux * face_area;
}

Referenced by BuildScalarDiffusionMappedTransfers_N(), and BuildScalarDiffusionMappedTransfers_S().

Here is the caller graph for this function:

ScalarDiffusionMappedVerticalTransfer()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real ScalarDiffusionMappedVerticalTransfer ( const amrex::Real  raw_flux, const amrex::Real  mx, const amrex::Real  my  )

noexcept

Convert a raw vertical diffusion flux to the mapped face-transfer space.

{
  return raw_flux / (mx * my);
}

Referenced by BuildScalarDiffusionMappedTransfers_N(), BuildScalarDiffusionMappedTransfers_S(), and TerrainScalarDiffusionMappedTransferAtKFace().

Here is the caller graph for this function:

ScaleScalarDiffusionVerticalFlux()

void ScaleScalarDiffusionVerticalFlux ( const amrex::Box &  zbx, const amrex::Array4< amrex::Real > &  raw_Fz, const int  flux_comp, const amrex::Real  explicit_fac  )

inline

Scale only raw F_z for the explicit integration fraction.

{
  amrex::ParallelFor(zbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
    raw_Fz(i, j, k, flux_comp) *= explicit_fac;
  });
}

Referenced by DiffusionSrcForState_N(), DiffusionSrcForState_S(), and DiffusionSrcForState_T().

Here is the call graph for this function:

Here is the caller graph for this function:

StretchedDiffusionDirichletWeights()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE StretchedDiffusionBoundaryWeights StretchedDiffusionDirichletWeights ( const amrex::Real  dz_near, const amrex::Real  dz_next  )

noexcept

{
  const amrex::Real f = dz_next / dz_near + amrex::Real(2.0);
  const amrex::Real f2 = f * f;
  const amrex::Real c3 = amrex::Real(2.0) / (f - f2);
  return {-(amrex::Real(1.0) - f2) * c3, -f2 * c3, c3};
}

Referenced by StretchedScalarGradient().

Here is the caller graph for this function:

StretchedScalarGradient()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real StretchedScalarGradient ( const amrex::Array4< const amrex::Real > &  scalar, const int  scalar_comp, const int  i, const int  j, const int  k, const amrex::Real *  dz, const int  klo, const int  khi, const bool  ext_dir_lo, const bool  ext_dir_hi  )

noexcept

Existing S-grid vertical-center gradient, including boundary and tile-edge spacing.

{
  if (ext_dir_lo) {
    const auto w = StretchedDiffusionDirichletWeights(dz[k], dz[k + 1]);
    return (amrex::Real(1.0) / dz[k]) *
           (w.c1 * scalar(i, j, k - 1, scalar_comp) +
            w.c2 * scalar(i, j, k, scalar_comp) +
            w.c3 * scalar(i, j, k + 1, scalar_comp));
  }
  if (ext_dir_hi) {
    const auto w = StretchedDiffusionDirichletWeights(dz[k - 1], dz[k - 2]);
    return (amrex::Real(1.0) / dz[k - 1]) *
           (-(
             w.c1 * scalar(i, j, k, scalar_comp) +
             w.c2 * scalar(i, j, k - 1, scalar_comp) +
             w.c3 * scalar(i, j, k - 2, scalar_comp)));
  }
  const amrex::Real inv = (k == klo) ? amrex::Real(1.0) / dz[k]
                          : (k == khi + 1)
                            ? amrex::Real(1.0) / dz[k - 1]
                            : amrex::Real(2.0) / (dz[k] + dz[k - 1]);
  return inv *
         (scalar(i, j, k, scalar_comp) - scalar(i, j, k - 1, scalar_comp));
}

Referenced by BuildScalarDiffusionFluxes_S_Impl().

Here is the call graph for this function:

Here is the caller graph for this function:

SubtractTerrainMetricImplicitPart_T()

void SubtractTerrainMetricImplicitPart_T ( const amrex::Box &  bx, const amrex::Box &  domain, const ScalarDiffusionFieldViews &  field, const ScalarDiffusionFluxPolicy &  policy, const amrex::Array4< const amrex::Real > &  z_nd, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &  cellSizeInv, const amrex::Array4< const amrex::Real > &  mf_mx, const amrex::Array4< const amrex::Real > &  mf_my, const amrex::Real  implicit_fac  )

inline

Dispatch SubtractTerrainMetricImplicitPart_T_Impl on the coefficient mode.

{
  const amrex::Box zbx = amrex::surroundingNodes(bx, 2);
  const auto mode = policy.coefficient_mode;
  if (mode.molecular_is_kinematic) {
    if (mode.add_turbulence) {
      SubtractTerrainMetricImplicitPart_T_Impl<true, true>(
        zbx, domain, field, policy.coefficients, z_nd, cellSizeInv, mf_mx, mf_my, implicit_fac);
    } else {
      SubtractTerrainMetricImplicitPart_T_Impl<true, false>(
        zbx, domain, field, policy.coefficients, z_nd, cellSizeInv, mf_mx, mf_my, implicit_fac);
    }
  } else if (mode.add_turbulence) {
    SubtractTerrainMetricImplicitPart_T_Impl<false, true>(
      zbx, domain, field, policy.coefficients, z_nd, cellSizeInv, mf_mx, mf_my, implicit_fac);
  } else {
    SubtractTerrainMetricImplicitPart_T_Impl<false, false>(
      zbx, domain, field, policy.coefficients, z_nd, cellSizeInv, mf_mx, mf_my, implicit_fac);
  }
}

Referenced by DiffusionSrcForState_T().

Here is the caller graph for this function:

SubtractTerrainMetricImplicitPart_T_Impl()

template<bool MultiplyMolecularByDensity, bool AddTurbulence>

void SubtractTerrainMetricImplicitPart_T_Impl ( const amrex::Box &  zbx, const amrex::Box &  domain, const ScalarDiffusionFieldViews &  field, const ScalarDiffusionCoefficients &  coefficients, const amrex::Array4< const amrex::Real > &  z_nd, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &  cellSizeInv, const amrex::Array4< const amrex::Real > &  mf_mx, const amrex::Array4< const amrex::Real > &  mf_my, const amrex::Real  implicit_fac  )

inline

erf.implicit_terrain_metric: the implicit vertical solve also takes the compact terrain-metric flux -K_h M_s dphi/dz (ERF_TerrainImplicitMetric.H) on the interior zeta faces, so add implicit_fac * K_h M_s dphi/dz to the explicit raw F_z there. K_h is the horizontal coefficient of the x/y fluxes, averaged to the z face; the bottom and top faces, which carry the boundary condition, are left alone.

{
  using namespace amrex;
  const int dom_klo = domain.smallEnd(2);
  const int dom_khi = domain.bigEnd(2);
  const auto scalar = field.scalar;
  const auto density = field.density;
  const auto mu_turb = field.mu_turb;
  const auto zflux = field.zflux;
  const int scalar_comp = field.scalar_comp;
  const int rho_comp = field.rho_comp;
  const int flux_comp = field.flux_comp;
  ParallelFor(zbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
    if (!TerrainMetricImplicitFace(k, dom_klo, dom_khi)) { return; }
    const Real kh = ScalarDiffusionFaceCoefficient<
      MultiplyMolecularByDensity, AddTurbulence>(
      density, rho_comp, mu_turb, coefficients, i, j, k, 0, 0, 1,
      coefficients.eddy_h_comp);
    const Real h_zeta = Compute_h_zeta_AtKface(i, j, k, cellSizeInv, z_nd);
    const Real grad_z = (cellSizeInv[2] / h_zeta) *
      (scalar(i, j, k, scalar_comp) - scalar(i, j, k - 1, scalar_comp));
    const Real metric = TerrainMetricScalarFactor(
      Compute_h_xi_AtKface(i, j, k, cellSizeInv, z_nd),
      Compute_h_eta_AtKface(i, j, k, cellSizeInv, z_nd),
      mf_mx(i, j, 0), mf_my(i, j, 0));
    zflux(i, j, k, flux_comp) += implicit_fac * kh * metric * grad_z;
  });
}

Here is the call graph for this function:

TerrainDiffusionDivergence_T()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real TerrainDiffusionDivergence_T ( const amrex::Real  mapped_Fx_hi, const amrex::Real  mapped_Fx_lo, const amrex::Real  mapped_Fy_hi, const amrex::Real  mapped_Fy_lo, const amrex::Real  Gz_hi, const amrex::Real  Gz_lo, const amrex::Real  dx_inv, const amrex::Real  dy_inv, const amrex::Real  dz_inv, const amrex::Real  mfsq, const amrex::Real  detJ_cell  )

noexcept

Native terrain divergence in the legacy arithmetic association.

{
  amrex::Real stateContrib =
    (mapped_Fx_hi - mapped_Fx_lo) * dx_inv * mfsq +
    (mapped_Fy_hi - mapped_Fy_lo) * dy_inv * mfsq +
    (Gz_hi - Gz_lo) * dz_inv;
  stateContrib /= detJ_cell;
  return stateContrib;
}

Referenced by ApplyScalarDiffusionMappedDivergence_T().

Here is the caller graph for this function:

TerrainDiffusionGz()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real TerrainDiffusionGz ( const amrex::Real  raw_Fz, const amrex::Real  mx, const amrex::Real  h_xi, const amrex::Real  Fbar_x, const amrex::Real  my, const amrex::Real  h_eta, const amrex::Real  Fbar_y  )

noexcept

Raw terrain-coordinate F_z and its conservative transformed transfer are distinct quantities. G_z is the numerator F_z - m_x h_xi Fbar_x - m_y h_eta Fbar_y; T_z=G_z/(m_x m_y).

{
  return raw_Fz - h_xi * mx * Fbar_x - h_eta * my * Fbar_y;
}

Referenced by TerrainScalarDiffusionGzAtKFace().

Here is the caller graph for this function:

TerrainDiffusionMappedTx()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real TerrainDiffusionMappedTx ( const amrex::Real  raw_Fx, const amrex::Real  h_zeta, const amrex::Real  my  )

noexcept

{
  return raw_Fx * h_zeta / my;
}

Referenced by ApplyScalarDiffusionMappedDivergence_T(), and BuildScalarDiffusionMappedTransfers_T().

Here is the caller graph for this function:

TerrainDiffusionMappedTy()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real TerrainDiffusionMappedTy ( const amrex::Real  raw_Fy, const amrex::Real  h_zeta, const amrex::Real  mx  )

noexcept

{
  return raw_Fy * h_zeta / mx;
}

Referenced by ApplyScalarDiffusionMappedDivergence_T(), and BuildScalarDiffusionMappedTransfers_T().

Here is the caller graph for this function:

TerrainScalarDiffusionFluxAtIFace()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real TerrainScalarDiffusionFluxAtIFace ( const amrex::Array4< const amrex::Real > &  scalar, const int  scalar_comp, const amrex::Real  coeff, const int  i, const int  j, const int  k, const amrex::Array4< const amrex::Real > &  z_nd, const amrex::Array4< const amrex::Real > &  z_cc, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &  cellSizeInv, const amrex::Array4< const amrex::Real > &  mf_ux  )

noexcept

Terrain x-face raw diffusion flux, including ERF's x-face map factor. coeff is the already selected effective face coefficient. Prescribed surface flux replacement is boundary policy and is intentionally outside this geometry primitive.

{
  const amrex::Real h_xi = Compute_h_xi_AtIface(i, j, k, cellSizeInv, z_nd);
  const amrex::Real idz_hi =
    amrex::Real(1.0) / (z_cc(i, j, k + 1) - z_cc(i, j, k - 1));
  const amrex::Real idz_lo =
    amrex::Real(1.0) / (z_cc(i - 1, j, k + 1) - z_cc(i - 1, j, k - 1));
  const amrex::Real grad_z =
    amrex::Real(0.5) * (scalar(i, j, k + 1, scalar_comp) * idz_hi +
                        scalar(i - 1, j, k + 1, scalar_comp) * idz_lo -
                        scalar(i, j, k - 1, scalar_comp) * idz_hi -
                        scalar(i - 1, j, k - 1, scalar_comp) * idz_lo);
  const amrex::Real grad_xi =
    cellSizeInv[0] *
    (scalar(i, j, k, scalar_comp) - scalar(i - 1, j, k, scalar_comp));
  return ScalarDiffusionFlux_Tx(coeff, mf_ux(i, j, 0), grad_xi, h_xi, grad_z);
}

Referenced by BuildScalarDiffusionFluxes_T_Impl().

Here is the call graph for this function:

Here is the caller graph for this function:

TerrainScalarDiffusionFluxAtJFace()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real TerrainScalarDiffusionFluxAtJFace ( const amrex::Array4< const amrex::Real > &  scalar, const int  scalar_comp, const amrex::Real  coeff, const int  i, const int  j, const int  k, const amrex::Array4< const amrex::Real > &  z_nd, const amrex::Array4< const amrex::Real > &  z_cc, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &  cellSizeInv, const amrex::Array4< const amrex::Real > &  mf_vy  )

noexcept

Terrain y-face raw diffusion flux, including ERF's y-face map factor.

{
  const amrex::Real h_eta = Compute_h_eta_AtJface(i, j, k, cellSizeInv, z_nd);
  const amrex::Real idz_hi =
    amrex::Real(1.0) / (z_cc(i, j, k + 1) - z_cc(i, j, k - 1));
  const amrex::Real idz_lo =
    amrex::Real(1.0) / (z_cc(i, j - 1, k + 1) - z_cc(i, j - 1, k - 1));
  const amrex::Real grad_z =
    amrex::Real(0.5) * (scalar(i, j, k + 1, scalar_comp) * idz_hi +
                        scalar(i, j - 1, k + 1, scalar_comp) * idz_lo -
                        scalar(i, j, k - 1, scalar_comp) * idz_hi -
                        scalar(i, j - 1, k - 1, scalar_comp) * idz_lo);
  const amrex::Real grad_eta =
    cellSizeInv[1] *
    (scalar(i, j, k, scalar_comp) - scalar(i, j - 1, k, scalar_comp));
  return ScalarDiffusionFlux_Ty(coeff, mf_vy(i, j, 0), grad_eta, h_eta, grad_z);
}

Referenced by BuildScalarDiffusionFluxes_T_Impl().

Here is the call graph for this function:

Here is the caller graph for this function:

TerrainScalarDiffusionFluxAtKFace()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real TerrainScalarDiffusionFluxAtKFace ( const amrex::Array4< const amrex::Real > &  scalar, const int  scalar_comp, const amrex::Real  coeff, const int  i, const int  j, const int  kface, const amrex::Array4< const amrex::Real > &  z_nd, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &  cellSizeInv, const bool  ext_dir_lo, const bool  ext_dir_hi  )

noexcept

Terrain raw vertical diffusion flux at one k face. Dirichlet boundary stencils are selected explicitly by the caller; surface-layer replacement remains separate boundary policy.

{
  amrex::Real grad_z;
  if (ext_dir_lo) {
    const amrex::Real zm = Compute_Z_AtWFace(i, j, kface + 1, z_nd);
    const amrex::Real dz0 = zm - Compute_Z_AtWFace(i, j, kface, z_nd);
    const amrex::Real dz1 = Compute_Z_AtWFace(i, j, kface + 2, z_nd) - zm;
    const amrex::Real idz0 = amrex::Real(1.0) / dz0;
    const amrex::Real f = dz1 / dz0 + amrex::Real(2.0);
    const amrex::Real f2 = f * f;
    const amrex::Real c3 = amrex::Real(2.0) / (f - f2);
    const amrex::Real c2 = -f2 * c3;
    const amrex::Real c1 = -(amrex::Real(1.0) - f2) * c3;
    grad_z = idz0 * (c1 * scalar(i, j, kface - 1, scalar_comp) +
                     c2 * scalar(i, j, kface, scalar_comp) +
                     c3 * scalar(i, j, kface + 1, scalar_comp));
  } else if (ext_dir_hi) {
    const amrex::Real zm = Compute_Z_AtWFace(i, j, kface - 1, z_nd);
    const amrex::Real dz0 = Compute_Z_AtWFace(i, j, kface, z_nd) - zm;
    const amrex::Real dz1 = zm - Compute_Z_AtWFace(i, j, kface - 2, z_nd);
    const amrex::Real idz0 = amrex::Real(1.0) / dz0;
    const amrex::Real f = dz1 / dz0 + amrex::Real(2.0);
    const amrex::Real f2 = f * f;
    const amrex::Real c3 = amrex::Real(2.0) / (f - f2);
    const amrex::Real c2 = -f2 * c3;
    const amrex::Real c1 = -(amrex::Real(1.0) - f2) * c3;
    grad_z = idz0 * (-(
                      c1 * scalar(i, j, kface, scalar_comp) +
                      c2 * scalar(i, j, kface - 1, scalar_comp) +
                      c3 * scalar(i, j, kface - 2, scalar_comp)));
  } else {
    const amrex::Real h_zeta =
      Compute_h_zeta_AtKface(i, j, kface, cellSizeInv, z_nd);
    grad_z = (cellSizeInv[2] / h_zeta) * (scalar(i, j, kface, scalar_comp) -
                                          scalar(i, j, kface - 1, scalar_comp));
  }
  return -coeff * grad_z;
}

Referenced by BuildScalarDiffusionFluxes_T_Impl().

Here is the call graph for this function:

Here is the caller graph for this function:

TerrainScalarDiffusionGzAtKFace()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real TerrainScalarDiffusionGzAtKFace ( const int  i, const int  j, const int  kface, const amrex::Box &  domain, const amrex::Array4< const amrex::Real > &  raw_xflux, const amrex::Array4< const amrex::Real > &  raw_yflux, const amrex::Array4< const amrex::Real > &  raw_zflux, const int  flux_comp, const amrex::Array4< const amrex::Real > &  z_nd, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &  cellSizeInv, const amrex::Array4< const amrex::Real > &  mf_mx, const amrex::Array4< const amrex::Real > &  mf_my  )

noexcept

Construct the native terrain numerator G_z at one vertical face. This helper owns x/y-face interpolation, physical-boundary extrapolation, and terrain cross terms; native divergence consumes G_z directly.

{
  const int klo = domain.smallEnd(2);
  const int khi = domain.bigEnd(2) + 1;
  amrex::Real xbar, ybar;
 
  if (kface == klo) {
    const amrex::Real xlo = myhalf * (raw_xflux(i, j, klo, flux_comp) +
                                    raw_xflux(i + 1, j, klo, flux_comp));
    const amrex::Real xhi = myhalf * (raw_xflux(i, j, klo + 1, flux_comp) +
                                    raw_xflux(i + 1, j, klo + 1, flux_comp));
    xbar = amrex::Real(1.5) * xlo - myhalf * xhi;
    const amrex::Real ylo = myhalf * (raw_yflux(i, j, klo, flux_comp) +
                                    raw_yflux(i, j + 1, klo, flux_comp));
    const amrex::Real yhi = myhalf * (raw_yflux(i, j, klo + 1, flux_comp) +
                                    raw_yflux(i, j + 1, klo + 1, flux_comp));
    ybar = amrex::Real(1.5) * ylo - myhalf * yhi;
  } else if (kface == khi) {
    const amrex::Real xlo = myhalf * (raw_xflux(i, j, khi - 2, flux_comp) +
                                    raw_xflux(i + 1, j, khi - 2, flux_comp));
    const amrex::Real xhi = myhalf * (raw_xflux(i, j, khi - 1, flux_comp) +
                                    raw_xflux(i + 1, j, khi - 1, flux_comp));
    xbar = amrex::Real(1.5) * xhi - myhalf * xlo;
    const amrex::Real ylo = myhalf * (raw_yflux(i, j, khi - 2, flux_comp) +
                                    raw_yflux(i, j + 1, khi - 2, flux_comp));
    const amrex::Real yhi = myhalf * (raw_yflux(i, j, khi - 1, flux_comp) +
                                    raw_yflux(i, j + 1, khi - 1, flux_comp));
    ybar = amrex::Real(1.5) * yhi - myhalf * ylo;
  } else {
    xbar = fourth * (raw_xflux(i, j, kface, flux_comp) +
                     raw_xflux(i + 1, j, kface, flux_comp) +
                     raw_xflux(i, j, kface - 1, flux_comp) +
                     raw_xflux(i + 1, j, kface - 1, flux_comp));
    ybar = fourth * (raw_yflux(i, j, kface, flux_comp) +
                     raw_yflux(i, j + 1, kface, flux_comp) +
                     raw_yflux(i, j, kface - 1, flux_comp) +
                     raw_yflux(i, j + 1, kface - 1, flux_comp));
  }
 
  const amrex::Real mx = mf_mx(i, j, 0);
  const amrex::Real my = mf_my(i, j, 0);
  const amrex::Real h_xi = Compute_h_xi_AtKface(i, j, kface, cellSizeInv, z_nd);
  const amrex::Real h_eta =
    Compute_h_eta_AtKface(i, j, kface, cellSizeInv, z_nd);
  const amrex::Real Gz = TerrainDiffusionGz(
    raw_zflux(i, j, kface, flux_comp), mx, h_xi, xbar, my, h_eta, ybar);
  return Gz;
}

Referenced by ApplyScalarDiffusionMappedDivergence_T(), and TerrainScalarDiffusionMappedTransferAtKFace().

Here is the call graph for this function:

Here is the caller graph for this function:

TerrainScalarDiffusionMappedTransferAtKFace()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real TerrainScalarDiffusionMappedTransferAtKFace ( const int  i, const int  j, const int  kface, const amrex::Box &  domain, const amrex::Array4< const amrex::Real > &  raw_xflux, const amrex::Array4< const amrex::Real > &  raw_yflux, const amrex::Array4< const amrex::Real > &  raw_zflux, const int  flux_comp, const amrex::Array4< const amrex::Real > &  z_nd, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &  cellSizeInv, const amrex::Array4< const amrex::Real > &  mf_mx, const amrex::Array4< const amrex::Real > &  mf_my, const bool  suppress_transfer = false  )

noexcept

Complete terrain vertical transfer at one k face. Terrain interpolation, boundary extrapolation and cross terms remain owned by the Gz helper above. Set suppress_transfer when the caller's resolved boundary policy suppresses the complete lower terrain-normal transfer.

At an interior k face, raw_xflux and raw_yflux must be valid at kface and kface-1. At the physical bottom or top face, the authoritative Gz helper uses one-sided extrapolation and needs no out-of-domain z layer.

{
  if (suppress_transfer) return amrex::Real(0.0);
  const amrex::Real Gz = TerrainScalarDiffusionGzAtKFace(
    i, j, kface, domain, raw_xflux, raw_yflux, raw_zflux, flux_comp, z_nd,
    cellSizeInv, mf_mx, mf_my);
  return ScalarDiffusionMappedVerticalTransfer(
    Gz, mf_mx(i, j, 0), mf_my(i, j, 0));
}

Referenced by BuildScalarDiffusionMappedTransfers_T().

Here is the call graph for this function:

Here is the caller graph for this function: