ERF_ImplicitDiff_S.cpp File Reference

#include "ERF_Diffusion.H"
#include "ERF_EddyViscosity.H"
#include "ERF_SolveTridiag.H"
#include "ERF_GetRhoAlpha.H"
#include "ERF_GetRhoAlphaForFaces.H"
#include "ERF_WallGradientCoeffs.H"
#include "ERF_NativeScalarDiffusion.H"
#include "ERF_SetupVertDiff.H"

Include dependency graph for ERF_ImplicitDiff_S.cpp:

Macros
#define	INSTANTIATE_IMPLICIT_DIFF_FOR_MOM_LU(STAGDIR)

Functions
void	ImplicitDiffForStateLU_S (const Box &bx, const Box &domain, const int level, const int n, const double dt_d, const GpuArray< Real, AMREX_SPACEDIM *2 > &bc_neumann_vals, const Array4< Real > &cell_data, const Gpu::DeviceVector< Real > &stretched_dz_d, const Array4< const Real > &scalar_zflux, const Array4< const Real > &mu_turb, const SolverChoice &solverChoice, const BCRec *bc_ptr, const bool use_SurfLayer, const Real implicit_fac, const bool use_mrf_countergradient)
template<int stagdir>
void	ImplicitDiffForMomLU_S (const Box &bx, const Box &domain, const int level, const double dt_d, const Array4< const int > &col_kext, const Array4< const Real > &cell_data, const Array4< Real > &face_data, const Array4< const Real > &tau, const Array4< const Real > &tau_corr, const Gpu::DeviceVector< Real > &stretched_dz_d, const Array4< const Real > &mu_turb, const SolverChoice &solverChoice, const BCRec *bc_ptr, const bool use_SurfLayer, const Real implicit_fac, const bool use_ysu_mom_countergradient)
void	ImplicitVerticalFluxForReflux_S (const Box &bx, const Box &domain, const int level, const int n, const Real implicit_fac, const Array4< const Real > &cell_data, const Gpu::DeviceVector< Real > &stretched_dz_d, const Array4< const Real > &mu_turb, const SolverChoice &solverChoice, const Array4< const Real > &mf_mx, const Array4< const Real > &mf_my, const bool use_mrf_countergradient, const Array4< Real > &zflux)

Macro Definition Documentation

INSTANTIATE_IMPLICIT_DIFF_FOR_MOM_LU

#define INSTANTIATE_IMPLICIT_DIFF_FOR_MOM_LU

(

STAGDIR

)

Value:

    template void ImplicitDiffForMomLU_S<STAGDIR> ( \
        const Box&, \
        const Box&, \
        const int, \
        const double, \
        const Array4<const int >&, \
        const Array4<const Real>&, \
        const Array4<      Real>&, \
        const Array4<const Real>&, \
        const Array4<const Real>&, \
        const Gpu::DeviceVector<Real>&, \
        const Array4<const Real>&, \
        const SolverChoice&, \
        const BCRec*, \
        const bool, \
        const Real, \
        const bool);

Function Documentation

ImplicitDiffForMomLU_S()

template<int stagdir>

void ImplicitDiffForMomLU_S	(	const Box &	bx,
		const Box &	domain,
		const int	level,
		const double	dt_d,
		const Array4< const int > &	col_kext,
		const Array4< const Real > &	cell_data,
		const Array4< Real > &	face_data,
		const Array4< const Real > &	tau,
		const Array4< const Real > &	tau_corr,
		const Gpu::DeviceVector< Real > &	stretched_dz_d,
		const Array4< const Real > &	mu_turb,
		const SolverChoice &	solverChoice,
		const BCRec *	bc_ptr,
		const bool	use_SurfLayer,
		const Real	implicit_fac,
		const bool	use_ysu_mom_countergradient
	)

Function for computing the implicit contribution to the vertical diffusion of momentum, with a vertically stretched grid over flat terrain.

This function (explicitly instantiated below) handles staggering in x or y through the template parameter, stagdir.

Parameters

[in]	bx	cell-centered box to loop over
[in]	level	AMR level
[in]	dt_d	time step
[in]	cell_data	conserved cell-centered rho
[in,out]	face_data	conserved momentum
[in]	tau	stress contribution to momentum
[in]	tau_corr	stress contribution to momentum that will be corrected by the implicit solve
[in]	stretched_dz_d	array over z of dz[k]
[in]	mu_turb	turbulent viscosity
[in]	solverChoice	container of parameters
[in]	bc_ptr	container with boundary conditions
[in]	use_SurfLayer	whether we have turned on subgrid diffusion
[in]	implicit_fac	if 1 then fully implicit; if 0 then fully explicit
[in]	use_ysu_mom_countergradient	whether to include YSU momentum countergradient correction

{
    BL_PROFILE_VAR("ImplicitDiffForMom_S()",ImplicitDiffForMom_S);
 
    Real dt = static_cast<Real>(dt_d);
 
    // setup quantities for getRhoAlphaAtFaces()
    DiffChoice dc = solverChoice.diffChoice;
    TurbChoice tc = solverChoice.turbChoice[level];
    bool l_consA  = (dc.molec_diff_type == MolecDiffType::ConstantAlpha);
    bool l_turb   = tc.use_kturb;
    // The correction strains for u/v contain a factor of 1/2, so mu_eff does not carry
    // the factor of 2 of the stresses
    Real mu_eff = (l_consA) ? dc.dynamic_viscosity / dc.rho0_trans
                            : dc.dynamic_viscosity;
 
    // offsets used to average to faces
    constexpr int ioff = (stagdir == 0) ? 1 : 0;
    constexpr int joff = (stagdir == 1) ? 1 : 0;
 
    // Box bounds
    int ilo = bx.smallEnd(0);
    int ihi = bx.bigEnd(0);
    int jlo = bx.smallEnd(1);
    int jhi = bx.bigEnd(1);
    int box_klo = bx.smallEnd(2);
    int box_khi = bx.bigEnd(2);
    amrex::ignore_unused(ilo, ihi, jlo, jhi);
 
    // The vertical extent of the domain.  A column whose top or bottom is not the domain
    // boundary ends at a coarse/fine boundary instead, and must not be given a physical BC row.
    const int dom_klo = domain.smallEnd(2);
    const int dom_khi = domain.bigEnd(2);
 
    // Temporary FABs for tridiagonal solve (allocated on column)
    //   A[k] * x[k-1] + B[k] * x[k] + C[k+1] = RHS[k]
    amrex::FArrayBox RHS_fab, soln_fab, coeffG_fab;
           RHS_fab.resize(bx,1, amrex::The_Async_Arena());
          soln_fab.resize(bx,1, amrex::The_Async_Arena());
        coeffG_fab.resize(bx,1, amrex::The_Async_Arena());
    auto const& RHS_a        =        RHS_fab.array();
    auto const& soln_a       =       soln_fab.array();
    auto const& coeffG_a     =     coeffG_fab.array();
 
    auto dz_ptr = stretched_dz_d.data();
 
    int bc_comp = BCVars::xvel_bc + stagdir;
    bool ext_dir_on_zlo  = (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir ||
                            bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim);
    bool ext_dir_on_zhi  = (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir ||
                            bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim);
    bool foextrap_on_zlo = (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap);
    bool foextrap_on_zhi = (bc_ptr[bc_comp].hi(2) == ERFBCType::foextrap);
    amrex::ignore_unused(foextrap_on_zlo,foextrap_on_zhi);
 
    AMREX_ASSERT_WITH_MESSAGE(foextrap_on_zlo || ext_dir_on_zlo || use_SurfLayer,
                              "Unexpected lower BC for momentum used with implicit vertical diffusion");
    AMREX_ASSERT_WITH_MESSAGE(foextrap_on_zhi || ext_dir_on_zhi,
                              "Unexpected upper BC for momentum used with implicit vertical diffusion");
 
    Real Fact = implicit_fac * dt;
 
#ifdef AMREX_USE_GPU
    ParallelFor(makeSlab(bx,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int)
    {
#else
    for (int j(jlo); j<=jhi; ++j) {
      for (int i(ilo); i<=ihi; ++i) {
#endif
 
          // The face column at (i,j) is shared by the two cell columns on either side of it,
          // and so by the box on either side of a grid seam.  Solve only over the vertical
          // range that both of those cell columns cover: the boxes on the two sides then
          // build the same tridiagonal system and get the same answer, so the duplicated
          // faces along the seam stay consistent.  Faces outside that range sit on a lateral
          // coarse/fine boundary and are set from the coarse level, so it is consistent to
          // leave them to the explicit update.
          const int klo = amrex::max(box_klo, col_kext(i-ioff,j-joff,0,0),
                                              col_kext(i     ,j     ,0,0));
          const int khi = amrex::min(box_khi, col_kext(i-ioff,j-joff,0,1),
                                              col_kext(i     ,j     ,0,1));
 
          // Fewer than two faces in common -- nothing to invert, leave the column explicit
#ifdef AMREX_USE_GPU
          if (khi <= klo) { return; }
#else
          if (khi <= klo) { continue; }
#endif
          // Notes:
          //
          // - In DiffusionSrcForMom (e.g., for x-mom)
          //
          //     Real diffContrib = ...
          //                      + (tau13(i,j,k+1) - tau13(i,j,k)) / dzinv
          //     rho_u_rhs(i,j,k) -= diffContrib;  // note the negative sign
          //
          // - We need to scale the explicit _part_ of `tau13` (for x-mom) by (1 - implicit_fac)
          //   The part that needs to be scaled is stored in `tau_corr`.
          //   E.g., tau13 = myhalf * (du/dz + dw/dx)
          //         tau13_corr = myhalf * du/dz
          //
          // - The momentum (`face_data`) was set to `S_old + S_rhs * dt`
          //   prior to including "ERF_Implicit.H". Recall that S_rhs includes
          //   sources from advection and other forcings, not just diffusion.
          //
          // - To correct momentum, we need to subtract `implicit_fac * diffContrib_corr`
          //   from S_rhs to recover `(1 - implicit_fac) * diffContrib_corr`,
          //   where `diffContrib_corr = -d(tau_corr)/dz`. The negative sign
          //   comes from our convention for the RHS diffusion source.
          //
          //   Subtracting a negative gives the += below; multiply by dt to
          //   get the intermediate momentum on the RHS of the tridiagonal
          //   system.
          //
          // - With a surface_layer BC, tau13/23 holds the vertical flux -d_z(k*u_i)
          //   directly. We must use tau at klo (not tau_corr) with SL BCs.
          //
          // - Finally, the terms ~ RHS += (tau_corr_hi - tau_corr_lo) / dz (below)
          //   essentially undo the explicit diffusion update that will be
          //   handled here implicitly.
 
          // Bottom boundary coefficients and RHS for L decomp
          //===================================================
          Real rhoface, rhoAlpha_lo, rhoAlpha_hi;
          Real dz_inv, dz_inv_lo, dz_inv_hi;
          Real a_tmp, b_tmp, c_tmp, inv_b2_tmp;
          {
              rhoface = myhalf * (cell_data(i,j,klo,Rho_comp) + cell_data(i-ioff,j-joff,klo,Rho_comp));
              getRhoAlphaForFaces(i, j, klo, ioff, joff, rhoAlpha_lo, rhoAlpha_hi,
                                  cell_data, mu_turb, mu_eff,
                                  l_consA, l_turb);
 
              dz_inv    = one / dz_ptr[klo];
              dz_inv_lo = dz_inv;
              dz_inv_hi = two / (dz_ptr[klo] + dz_ptr[klo+1]);
 
              a_tmp = zero;
              c_tmp = -Fact * rhoAlpha_hi * dz_inv_hi * dz_inv;
 
              RHS_a(i,j,klo) = face_data(i,j,klo); // NOTE: this is momenta; solution is velocity
 
              // BCs: Dirichlet (u_i = val), slip wall (w = 0), or surface layer (w = 0) --
              //      but only where the bottom of this column really is the domain boundary
              if (klo != dom_klo) {
                  // Coarse/fine (or grid) boundary: close the system with the known face value
                  // below instead of imposing a physical BC at an interior height
                  dz_inv_lo = two / (dz_ptr[klo] + dz_ptr[klo-1]);
                  a_tmp = -Fact * rhoAlpha_lo * dz_inv_lo * dz_inv;
                  RHS_a(i,j,klo) += Fact * (tau_corr(i,j,klo+1) - tau_corr(i,j,klo)) * dz_inv;
                  const Real rho_below = myhalf * ( cell_data(i     ,j     ,klo-1,Rho_comp)
                                                  + cell_data(i-ioff,j-joff,klo-1,Rho_comp) );
                  RHS_a(i,j,klo) -= a_tmp * (face_data(i,j,klo-1) / rho_below);
              } else if (ext_dir_on_zlo) {
                  RHS_a(i,j,klo) += Fact * (tau_corr(i,j,klo+1) - tau_corr(i,j,klo)) * dz_inv;
                  // Wall stress with the one-sided gradient of ComputeStrain_S,
                  // (c_wall u_wall + c_0 u_0 + c_1 u_1) / dz[klo]: a_tmp multiplies the
                  // wall value, and c_1 adds to the coupling with the cell above
                  const WallGradCoeffs wc = WallGradientCoeffs(dz_ptr[klo], dz_ptr[klo+1]);
                  const Real idz0 = one / dz_ptr[klo];
                  a_tmp  = Fact * rhoAlpha_lo * idz0 * dz_inv * wc.c_wall;
                  c_tmp += Fact * rhoAlpha_lo * idz0 * dz_inv * wc.c_1;
                  const Real rho_wall = myhalf * ( cell_data(i     ,j     ,klo-1,Rho_comp)
                                                 + cell_data(i-ioff,j-joff,klo-1,Rho_comp) );
                  const Real wall_velocity = face_data(i,j,klo-1) / rho_wall;
                  RHS_a(i,j,klo) -= a_tmp * wall_velocity;
              } else if (use_SurfLayer) {
                  // NOTE: tau = -mu*d_z(u_i) w/ SL
                  RHS_a(i,j,klo) += Fact * (tau_corr(i,j,klo+1) - tau(i,j,klo)) * dz_inv;
                  RHS_a(i,j,klo) += Fact * dz_inv * tau(i,j,klo);
              } else {
                  // NOTE: FOEXTRAP has zero lower flux (nothing to add to RHS)
                  RHS_a(i,j,klo) += Fact * (tau_corr(i,j,klo+1) - tau_corr(i,j,klo)) * dz_inv;
              }
 
              // Add YSU momentum countergradient correction at bottom boundary.
              // NOTE: As for the scalars, the lower face at the domain bottom carries no
              //       countergradient flux -- the surface stress is supplied by the
              //       surface layer model or the wall BC above.  At a coarse/fine boundary
              //       the lower face is an interior face and carries it, as in the
              //       interior rows (erf-model/ERF#4822).
              // NOTE: The sign matches the scalar path: tau_i3 = -rho*K*(du_i/dz - gamma_i),
              //       so the countergradient piece of the flux is +rho*K*gamma_i and its
              //       divergence enters the RHS with a minus sign.
              if (use_ysu_mom_countergradient) {
                  const int hgam_comp = (stagdir == 0) ? EddyDiff::HGAMU_v : EddyDiff::HGAMV_v;
                  // Average HGAM* to the staggered face
                  const Real gam_klo = myhalf * (mu_turb(i,j,klo  ,hgam_comp) + mu_turb(i-ioff,j-joff,klo  ,hgam_comp));
                  const Real gam_kp1 = myhalf * (mu_turb(i,j,klo+1,hgam_comp) + mu_turb(i-ioff,j-joff,klo+1,hgam_comp));
                  const Real gam_hi  = myhalf * (gam_klo + gam_kp1);
                  RHS_a(i,j,klo) -= Fact * dz_inv * rhoAlpha_hi * gam_hi;
                  if (klo != dom_klo) {
                      const Real gam_km1 = myhalf * (mu_turb(i,j,klo-1,hgam_comp) + mu_turb(i-ioff,j-joff,klo-1,hgam_comp));
                      const Real gam_lo  = myhalf * (gam_klo + gam_km1);
                      RHS_a(i,j,klo) += Fact * dz_inv * rhoAlpha_lo * gam_lo;
                  }
              }
 
              b_tmp      = rhoface - a_tmp - c_tmp;
              inv_b2_tmp = one;
 
              RHS_a(i,j,klo)    /= b_tmp;         // NOTE: this is now "rho"
              coeffG_a(i,j,klo)  = c_tmp / b_tmp; // NOTE: this is now "gamma"
          }
 
          // Build the coefficients and RHS for L decomp
          //===================================================
          for (int k(klo+1); k < khi; k++) {
              rhoface = myhalf * (cell_data(i,j,k,Rho_comp) + cell_data(i-ioff,j-joff,k,Rho_comp));
              getRhoAlphaForFaces(i, j, k, ioff, joff, rhoAlpha_lo, rhoAlpha_hi,
                                  cell_data, mu_turb, mu_eff,
                                  l_consA, l_turb);
 
              dz_inv    = one / dz_ptr[k];
              dz_inv_lo = two / (dz_ptr[k] + dz_ptr[k-1]);
              dz_inv_hi = two / (dz_ptr[k] + dz_ptr[k+1]);
 
              a_tmp      = -Fact * rhoAlpha_lo * dz_inv_lo * dz_inv;
              c_tmp      = -Fact * rhoAlpha_hi * dz_inv_hi * dz_inv;
              b_tmp      = rhoface - a_tmp - c_tmp;
              inv_b2_tmp = one / (b_tmp - a_tmp * coeffG_a(i,j,k-1));
 
              RHS_a(i,j,k)    = face_data(i,j,k); // NOTE: this is momenta; solution is velocity
              RHS_a(i,j,k)   += Fact * (tau_corr(i,j,k+1) - tau_corr(i,j,k)) * dz_inv;
 
              // Add YSU momentum countergradient correction
              if (use_ysu_mom_countergradient) {
                  const int hgam_comp = (stagdir == 0) ? EddyDiff::HGAMU_v : EddyDiff::HGAMV_v;
                  // AverageHGAM* to the staggered face
                  const Real gam_k   = myhalf * (mu_turb(i,j,k  ,hgam_comp) + mu_turb(i-ioff,j-joff,k  ,hgam_comp));
                  const Real gam_km1 = myhalf * (mu_turb(i,j,k-1,hgam_comp) + mu_turb(i-ioff,j-joff,k-1,hgam_comp));
                  const Real gam_kp1 = myhalf * (mu_turb(i,j,k+1,hgam_comp) + mu_turb(i-ioff,j-joff,k+1,hgam_comp));
                  const Real gam_hi  = myhalf * (gam_k + gam_kp1); // at k+1/2
                  const Real gam_lo  = myhalf * (gam_k + gam_km1); // at k-1/2
                  RHS_a(i,j,k) -= Fact * dz_inv * (rhoAlpha_hi * gam_hi - rhoAlpha_lo * gam_lo);
              }
 
              RHS_a(i,j,k)    = (RHS_a(i,j,k) - a_tmp * RHS_a(i,j,k-1)) * inv_b2_tmp; // NOTE: This is now "rho"
              coeffG_a(i,j,k) = c_tmp * inv_b2_tmp; // NOTE: this is now "gamma"
          } // k
 
          // Top boundary coefficients and RHS for L decomp
          //===================================================
          {
              rhoface = myhalf * (cell_data(i,j,khi,Rho_comp) + cell_data(i-ioff,j-joff,khi,Rho_comp));
              getRhoAlphaForFaces(i, j, khi, ioff, joff, rhoAlpha_lo, rhoAlpha_hi,
                                  cell_data, mu_turb, mu_eff,
                                  l_consA, l_turb);
 
              dz_inv    = one / dz_ptr[khi];
              dz_inv_lo = two / (dz_ptr[khi] + dz_ptr[khi-1]);
              dz_inv_hi = dz_inv;
 
              a_tmp = -Fact * rhoAlpha_lo * dz_inv_lo * dz_inv;
              c_tmp = zero;
 
              RHS_a(i,j,khi)  = face_data(i,j,khi); // NOTE: this is momenta; solution is velocity
              RHS_a(i,j,khi) += Fact * (tau_corr(i,j,khi+1) - tau_corr(i,j,khi)) * dz_inv;
 
              // BCs: Dirichlet (u_i = val), slip wall (w = 0) --
              //      but only where the top of this column really is the domain boundary
              if (khi != dom_khi) {
                  // Coarse/fine (or grid) boundary: close the system with the known face value
                  // above instead of imposing a physical BC at an interior height
                  dz_inv_hi = two / (dz_ptr[khi] + dz_ptr[khi+1]);
                  c_tmp = -Fact * rhoAlpha_hi * dz_inv_hi * dz_inv;
                  const Real rho_above = myhalf * ( cell_data(i     ,j     ,khi+1,Rho_comp)
                                                  + cell_data(i-ioff,j-joff,khi+1,Rho_comp) );
                  RHS_a(i,j,khi) -= c_tmp * (face_data(i,j,khi+1) / rho_above);
              } else if (ext_dir_on_zhi) {
                  // Wall stress with the one-sided gradient of ComputeStrain_S: c_tmp
                  // multiplies the wall value, and c_1 adds to the coupling with the
                  // cell below
                  const WallGradCoeffs wc = WallGradientCoeffs(dz_ptr[khi], dz_ptr[khi-1]);
                  const Real idz0 = one / dz_ptr[khi];
                  c_tmp  = Fact * rhoAlpha_hi * idz0 * dz_inv * wc.c_wall;
                  a_tmp += Fact * rhoAlpha_hi * idz0 * dz_inv * wc.c_1;
                  const Real rho_wall = myhalf * ( cell_data(i     ,j     ,khi+1,Rho_comp)
                                                 + cell_data(i-ioff,j-joff,khi+1,Rho_comp) );
                  const Real wall_velocity = face_data(i,j,khi+1) / rho_wall;
                  RHS_a(i,j,khi) -= c_tmp * wall_velocity;
              }
 
              // Add YSU momentum countergradient correction at top boundary.
              // NOTE: The lower face at khi is the interior face shared with row khi-1, which
              //       has already taken this flux out of the cell below, so it always carries
              //       it; without it here the flux left the column (erf-model/ERF#4822).  The
              //       upper face carries it only at a coarse/fine boundary -- the domain top
              //       is a lid for gamma, as for the scalars.
              if (use_ysu_mom_countergradient) {
                  const int hgam_comp = (stagdir == 0) ? EddyDiff::HGAMU_v : EddyDiff::HGAMV_v;
                  // Average HGAM* to the staggered face
                  const Real gam_khi = myhalf * (mu_turb(i,j,khi  ,hgam_comp) + mu_turb(i-ioff,j-joff,khi  ,hgam_comp));
                  const Real gam_km1 = myhalf * (mu_turb(i,j,khi-1,hgam_comp) + mu_turb(i-ioff,j-joff,khi-1,hgam_comp));
                  const Real gam_lo  = myhalf * (gam_khi + gam_km1);
                  RHS_a(i,j,khi) += Fact * dz_inv * rhoAlpha_lo * gam_lo;
                  if (khi != dom_khi) {
                      const Real gam_kp1 = myhalf * (mu_turb(i,j,khi+1,hgam_comp) + mu_turb(i-ioff,j-joff,khi+1,hgam_comp));
                      const Real gam_hi  = myhalf * (gam_khi + gam_kp1);
                      RHS_a(i,j,khi) -= Fact * dz_inv * rhoAlpha_hi * gam_hi;
                  }
              }
 
              b_tmp      = rhoface - a_tmp - c_tmp;
              inv_b2_tmp = one / (b_tmp - a_tmp * coeffG_a(i,j,khi-1));
 
              // First solve
              soln_a(i,j,khi) = (RHS_a(i,j,khi) - a_tmp * RHS_a(i,j,khi-1)) * inv_b2_tmp;
          }
 
          // Back sweep the U decomp solution
          //===================================================
          for (int k(khi-1); k>=klo; --k) {
              soln_a(i,j,k) = RHS_a(i,j,k) - coeffG_a(i,j,k) * soln_a(i,j,k+1);
          }
 
          // Convert back to momenta
          //===================================================
          for (int k(klo); k<=khi; ++k) {
              rhoface = myhalf * (cell_data(i,j,k,Rho_comp) + cell_data(i-ioff,j-joff,k,Rho_comp));
              face_data(i,j,k) = rhoface * soln_a(i,j,k);
          }
 
#ifdef AMREX_USE_GPU
    });
#else
      } // i
    } // j
#endif
}

Here is the call graph for this function:

ImplicitDiffForStateLU_S()

void ImplicitDiffForStateLU_S	(	const Box &	bx,
		const Box &	domain,
		const int	level,
		const int	n,
		const double	dt_d,
		const GpuArray< Real, AMREX_SPACEDIM *2 > &	bc_neumann_vals,
		const Array4< Real > &	cell_data,
		const Gpu::DeviceVector< Real > &	stretched_dz_d,
		const Array4< const Real > &	scalar_zflux,
		const Array4< const Real > &	mu_turb,
		const SolverChoice &	solverChoice,
		const BCRec *	bc_ptr,
		const bool	use_SurfLayer,
		const Real	implicit_fac,
		const bool	use_mrf_countergradient
	)

Function for computing the implicit contribution to the vertical diffusion of theta, with a vertically stretched grid over flat terrain.

Parameters

[in]	bx	cell-centered box to loop over
[in]	domain	box of the whole domain
[in]	level	AMR level
[in]	n	conserved component index
[in]	dt_d	time step
[in]	bc_neumann_vals	values of derivatives if bc_type == Neumann
[in,out]	cell_data	conserved cell-centered rho, rho theta
[in]	stretched_dz_d	array over z of dz[k]
[in]	scalar_zflux	scalar vertical flux in z-dir
[in]	mu_turb	turbulent viscosity
[in]	solverChoice	container of parameters
[in]	bc_ptr	container with boundary conditions
[in]	use_SurfLayer	whether we have turned on subgrid diffusion
[in]	implicit_fac	if 1 then fully implicit; if 0 then fully explicit
[in]	use_mrf_countergradient	whether to include MRF countergradient correction

{
    BL_PROFILE_VAR("ImplicitDiffForState_S()",ImplicitDiffForState_S);
 
    Real dt = static_cast<Real>(dt_d);
 
    // setup quantities for getRhoAlpha()
#include "ERF_SetupVertDiff.H"
    const int qty_index  = n;
    const int prim_index = qty_index - 1;
    const int prim_scal_index = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? PrimScalar_comp : prim_index;
 
    // Box bounds
    int ilo = bx.smallEnd(0);
    int ihi = bx.bigEnd(0);
    int jlo = bx.smallEnd(1);
    int jhi = bx.bigEnd(1);
    int klo = bx.smallEnd(2);
    int khi = bx.bigEnd(2);
    amrex::ignore_unused(ilo, ihi, jlo, jhi);
 
    // Temporary FABs for tridiagonal solve (allocated on column)
    //   A[k] * x[k-1] + B[k] * x[k] + C[k+1] = RHS[k]
 
    // With LU decomposition, M * x = r is written as L * U * x = r with U * x = rho
    // We then first have L * rho = r and U * x = rho
    amrex::FArrayBox RHS_fab, soln_fab, coeffG_fab;
           RHS_fab.resize(bx,1, amrex::The_Async_Arena());
          soln_fab.resize(bx,1, amrex::The_Async_Arena());
        coeffG_fab.resize(bx,1, amrex::The_Async_Arena());
    auto const& RHS_a        =        RHS_fab.array();
    auto const& soln_a       =       soln_fab.array();
    auto const& coeffG_a     =     coeffG_fab.array();
 
    auto dz_ptr = stretched_dz_d.data();
 
    // BC types are stored by BC variable, not by state component
    const int bc_comp = ResolveNativeScalarDiffusionPolicy(qty_index, solverChoice.diffChoice).bc_comp;
    bool foextrap_on_zlo = (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap);
    bool foextrap_on_zhi = (bc_ptr[bc_comp].hi(2) == ERFBCType::foextrap);
    bool neumann_on_zlo  = (bc_ptr[bc_comp].lo(2) == ERFBCType::neumann);
    bool neumann_on_zhi  = (bc_ptr[bc_comp].hi(2) == ERFBCType::neumann);
    amrex::ignore_unused(foextrap_on_zlo, foextrap_on_zhi);
 
    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(foextrap_on_zlo || neumann_on_zlo || use_SurfLayer,
                                     "Unexpected lower BC for scalars used with implicit vertical diffusion");
    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(foextrap_on_zhi || neumann_on_zhi,
                                     "Unexpected upper BC for scalars used with implicit vertical diffusion");
 
    // k-eqn RANS with a Dirichlet wall value of k: the first cell is held at
    // the value set by the surface layer, so its row reduces to x(klo) = phi(klo)
    // and the row above sees it as a Dirichlet neighbour.
    // NOTE: init_bcs requires zlo.type = surface_layer whenever dirichlet_k is set,
    //       so this matches the l_dirichlet_k pin in erf_slow_rhs_post, which also
    //       requires a surface layer. Do not relax one without the other.
    const bool pin_klo = (qty_index == RhoKE_comp) && (klo == domain.smallEnd(2)) &&
                         (solverChoice.turbChoice[level].rans_type == RANSType::kEqn) &&
                         solverChoice.turbChoice[level].dirichlet_k;
 
    // A column whose bottom or top is not the domain boundary ends at a coarse/fine boundary
    // instead.  There the system must be closed with the value on the other side of that
    // boundary rather than given a physical BC row at an interior height, so that the answer
    // does not depend on where the grid happens to end in the vertical.
    const bool at_zlo = (klo == domain.smallEnd(2));
    const bool at_zhi = (khi == domain.bigEnd(2));
 
    Real Fact = implicit_fac * dt;
 
#ifdef AMREX_USE_GPU
    ParallelFor(makeSlab(bx,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int)
    {
#else
    for (int j(jlo); j<=jhi; ++j) {
        for (int i(ilo); i<=ihi; ++i) {
#endif
            // Bottom boundary coefficients and RHS for L decomp
            //===================================================
            Real rhoAlpha_lo, rhoAlpha_hi;
            Real dz_inv, dz_inv_lo, dz_inv_hi;
            Real a_tmp, b_tmp, c_tmp, inv_b2_tmp;
            {
                getRhoAlpha(i, j, klo, rhoAlpha_lo, rhoAlpha_hi,
                            cell_data, mu_turb, d_alpha_eff, d_eddy_diff_idz,
                            prim_scal_index, l_consA, l_turb);
 
                dz_inv    = one / dz_ptr[klo];
                dz_inv_lo = dz_inv;
                dz_inv_hi = two / (dz_ptr[klo] + dz_ptr[klo+1]);
 
                a_tmp      = zero;
                if (!at_zlo) {
                    dz_inv_lo = two / (dz_ptr[klo] + dz_ptr[klo-1]);
                    a_tmp     = -Fact * rhoAlpha_lo * dz_inv_lo * dz_inv;
                }
                c_tmp      = (pin_klo) ? zero : -Fact * rhoAlpha_hi * dz_inv_hi * dz_inv;
                b_tmp      = cell_data(i,j,klo,Rho_comp) - a_tmp - c_tmp;
                inv_b2_tmp = one;
 
                RHS_a(i,j,klo) = cell_data(i,j,klo,n); // NOTE: this is rho*phi; solution is phi
                if (!at_zlo) {
                    // Coarse/fine boundary: close the system with the known value below
                    RHS_a(i,j,klo) -= a_tmp * (cell_data(i,j,klo-1,n) / cell_data(i,j,klo-1,Rho_comp));
                } else if (pin_klo) {
                    // Dirichlet row: no flux terms
                } else if (use_SurfLayer && scalar_zflux) {
                    RHS_a(i,j,klo) +=  Fact * dz_inv * scalar_zflux(i,j,klo); // NOTE: scalar_zflux = -K*d_z(\phi)
                } else if (neumann_on_zlo) {
                    RHS_a(i,j,klo) += -Fact * dz_inv * rhoAlpha_lo * bc_neumann_vals[2]; // NOTE: N_val = d_z(\phi)
                }
 
                // Add countergradient correction to RHS at bottom boundary.
                // NOTE: The lower face at klo carries no countergradient flux -- the
                //       total surface flux is supplied by the surface layer model or the
                //       Neumann BC, while gamma represents nonlocal transport interior to
                //       the PBL.  Only the upper face contributes here.
                if (use_mrf_countergradient && (n == RhoTheta_comp || n == RhoQ1_comp)) {
                    const int gam_comp = (n == RhoTheta_comp) ? EddyDiff::HGAMT_v : EddyDiff::HGAMQ_v;
                    const Real gam_hi = myhalf * (mu_turb(i, j, klo, gam_comp) + mu_turb(i, j, klo+1, gam_comp));
                    // rhoAlpha*gam is already a flux, so its divergence over cell klo
                    // is scaled by the *cell* spacing, not the face spacing
                    RHS_a(i,j,klo) -= Fact * dz_inv * rhoAlpha_hi * gam_hi;
                    if (!at_zlo) {
                        const Real gam_lo = myhalf * (mu_turb(i, j, klo, gam_comp) + mu_turb(i, j, klo-1, gam_comp));
                        RHS_a(i,j,klo) += Fact * dz_inv * rhoAlpha_lo * gam_lo;
                    }
                }
 
                RHS_a(i,j,klo)    /= b_tmp;         // NOTE: this is now "rho"
                coeffG_a(i,j,klo)  = c_tmp / b_tmp; // NOTE: this is now "gamma"
            }
 
            // Build the coefficients and RHS for L decomp
            //===================================================
            for (int k(klo+1); k < khi; k++) {
                getRhoAlpha(i, j, k, rhoAlpha_lo, rhoAlpha_hi,
                            cell_data, mu_turb, d_alpha_eff, d_eddy_diff_idz,
                            prim_scal_index, l_consA, l_turb);
 
                dz_inv    = one / dz_ptr[k];
                dz_inv_lo = two / (dz_ptr[k] + dz_ptr[k-1]);
                dz_inv_hi = two / (dz_ptr[k] + dz_ptr[k+1]);
 
                a_tmp      = -Fact * rhoAlpha_lo * dz_inv_lo * dz_inv;
                c_tmp      = -Fact * rhoAlpha_hi * dz_inv_hi * dz_inv;
                b_tmp      = cell_data(i,j,k,Rho_comp) - a_tmp - c_tmp;
                inv_b2_tmp = one / (b_tmp - a_tmp * coeffG_a(i,j,k-1));
 
                RHS_a(i,j,k)    = cell_data(i,j,k,n); // NOTE: this is rho*phi; solution is phi
 
                // Add countergradient correction to RHS in interior
                if (use_mrf_countergradient && (n == RhoTheta_comp || n == RhoQ1_comp)) {
                    const int gam_comp = (n == RhoTheta_comp) ? EddyDiff::HGAMT_v : EddyDiff::HGAMQ_v;
                    const Real gam_k   = mu_turb(i, j, k,   gam_comp);
                    const Real gam_km1 = mu_turb(i, j, k-1, gam_comp);
                    const Real gam_kp1 = mu_turb(i, j, k+1, gam_comp);
                    const Real gam_hi  = myhalf * (gam_k + gam_kp1); // at k+½
                    const Real gam_lo  = myhalf * (gam_k + gam_km1); // at k-½
                    // Countergradient flux divergence (implicit contribution to RHS):
                    //   -Fact * dz_inv * [ρα_{k+½}·γ_{k+½} - ρα_{k-½}·γ_{k-½}]
                    // ρα·γ is already a flux, so differencing it across the cell is
                    // scaled by the *cell* spacing dz_inv.  The face spacings
                    // dz_inv_hi/dz_inv_lo belong to the gradient in the diffusion
                    // coefficients above, not here; using them would leave a spurious
                    // tendency for a constant countergradient flux on a stretched grid.
                    RHS_a(i,j,k) -= Fact * dz_inv * (rhoAlpha_hi * gam_hi - rhoAlpha_lo * gam_lo);
                }
 
                RHS_a(i,j,k)    = (RHS_a(i,j,k) - a_tmp * RHS_a(i,j,k-1)) * inv_b2_tmp; // NOTE: This is now "rho"
                coeffG_a(i,j,k) = c_tmp * inv_b2_tmp; // NOTE: this is now "gamma"
            } // k
 
            // Top boundary coefficients and RHS for L decomp
            //===================================================
            {
                getRhoAlpha(i, j, khi, rhoAlpha_lo, rhoAlpha_hi,
                            cell_data, mu_turb, d_alpha_eff, d_eddy_diff_idz,
                            prim_scal_index, l_consA, l_turb);
 
                dz_inv    = one / dz_ptr[khi];
                dz_inv_lo = two / (dz_ptr[khi] + dz_ptr[khi-1]);
                dz_inv_hi = dz_inv;
 
                a_tmp      = -Fact * rhoAlpha_lo * dz_inv_lo * dz_inv;
                c_tmp      = zero;
                if (!at_zhi) {
                    dz_inv_hi = two / (dz_ptr[khi] + dz_ptr[khi+1]);
                    c_tmp     = -Fact * rhoAlpha_hi * dz_inv_hi * dz_inv;
                }
                b_tmp      = cell_data(i,j,khi,Rho_comp) - a_tmp - c_tmp;
                inv_b2_tmp = one / (b_tmp - a_tmp * coeffG_a(i,j,khi-1));
 
                RHS_a(i,j,khi) = cell_data(i,j,khi,n); // NOTE: this is rho*phi; solution is phi
                if (!at_zhi) {
                    // Coarse/fine boundary: close the system with the known value above
                    RHS_a(i,j,khi) -= c_tmp * (cell_data(i,j,khi+1,n) / cell_data(i,j,khi+1,Rho_comp));
                } else if (neumann_on_zhi) {
                    RHS_a(i,j,khi) -= -Fact * dz_inv * rhoAlpha_hi * bc_neumann_vals[5]; // NOTE: N_val = d_z(\phi)
                }
 
                // Add countergradient correction to RHS at top boundary.
                // NOTE: The lower face at khi always carries the countergradient flux, which
                //       row khi-1 has already taken out of the cell below; without it here
                //       that flux would leave the column.  The upper face carries it only at
                //       a coarse/fine boundary -- the domain top is a lid for gamma.
                if (use_mrf_countergradient && (n == RhoTheta_comp || n == RhoQ1_comp)) {
                    const int gam_comp = (n == RhoTheta_comp) ? EddyDiff::HGAMT_v : EddyDiff::HGAMQ_v;
                    const Real gam_lo = myhalf * (mu_turb(i, j, khi, gam_comp) + mu_turb(i, j, khi-1, gam_comp));
                    RHS_a(i,j,khi) += Fact * dz_inv * rhoAlpha_lo * gam_lo;
                    if (!at_zhi) {
                        const Real gam_hi = myhalf * (mu_turb(i, j, khi, gam_comp) + mu_turb(i, j, khi+1, gam_comp));
                        RHS_a(i,j,khi) -= Fact * dz_inv * rhoAlpha_hi * gam_hi;
                    }
                }
 
                // First solve
                soln_a(i,j,khi) = (RHS_a(i,j,khi) - a_tmp * RHS_a(i,j,khi-1)) * inv_b2_tmp;
            }
 
            // Back sweep the U decomp solution
            //===================================================
            for (int k(khi-1); k>=klo; --k) {
                soln_a(i,j,k) = RHS_a(i,j,k) - coeffG_a(i,j,k) * soln_a(i,j,k+1);
            }
 
            // Convert back to rho*theta
            //===================================================
            for (int k(klo); k<=khi; ++k) {
                cell_data(i,j,k,n) = cell_data(i,j,k,Rho_comp) * soln_a(i,j,k);
            }
 
#ifdef AMREX_USE_GPU
    });
#else
      } // i
    } // j
#endif
}

Here is the call graph for this function:

ImplicitVerticalFluxForReflux_S()

void ImplicitVerticalFluxForReflux_S	(	const Box &	bx,
		const Box &	domain,
		const int	level,
		const int	n,
		const Real	implicit_fac,
		const Array4< const Real > &	cell_data,
		const Gpu::DeviceVector< Real > &	stretched_dz_d,
		const Array4< const Real > &	mu_turb,
		const SolverChoice &	solverChoice,
		const Array4< const Real > &	mf_mx,
		const Array4< const Real > &	mf_my,
		const bool	use_mrf_countergradient,
		const Array4< Real > &	zflux
	)

The implicit part of the vertical diffusive flux of component n on the z faces of bx, as the matching ImplicitDiffForStateLU_S solve applied it, for refluxing: in the convention of the advective fluxes the flux registers hold (see AddScalarDiffusionRefluxFluxes_S), i.e. the vertical transfer divided by m_x m_y.

cell_data must hold the solution of that solve, including the ghost values a column that ends at a coarse/fine boundary was closed with. Faces on the domain boundary are set to zero: they are never refluxed.

Parameters

[in]	bx	cell-centered box
[in]	domain	box of the whole domain
[in]	level	AMR level
[in]	n	conserved component index
[in]	implicit_fac	the vert_implicit_fac the solve used
[in]	cell_data	conserved state after the solve
[in]	use_mrf_countergradient	whether the solve added the MRF countergradient flux
[out]	zflux	face flux of component n, in component 0

{
#include "ERF_SetupVertDiff.H"
    const int qty_index  = n;
    const int prim_index = qty_index - 1;
    const int prim_scal_index = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? PrimScalar_comp : prim_index;
    const int khi = bx.bigEnd(2);
    auto dz_ptr = stretched_dz_d.data();
 
    ParallelFor(surroundingNodes(bx,2), [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
    {
        if (k == dom_lo.z || k == dom_hi.z+1) { zflux(i,j,k) = zero; return; }
        // The coefficient on this face, from the cell above it (or below, at the top of bx)
        Real rhoAlpha_lo, rhoAlpha_hi;
        const int kc = (k <= khi) ? k : k-1;
        getRhoAlpha(i, j, kc, rhoAlpha_lo, rhoAlpha_hi, cell_data, mu_turb, d_alpha_eff, d_eddy_diff_idz,
                    prim_scal_index, l_consA, l_turb);
        const Real rhoAlpha = (k <= khi) ? rhoAlpha_lo : rhoAlpha_hi;
        const Real dz_inv_face = two / (dz_ptr[k] + dz_ptr[k-1]);
        const Real phi_hi = cell_data(i,j,k  ,n) / cell_data(i,j,k  ,Rho_comp);
        const Real phi_lo = cell_data(i,j,k-1,n) / cell_data(i,j,k-1,Rho_comp);
        Real G = -implicit_fac * rhoAlpha * dz_inv_face * (phi_hi - phi_lo);
        // The countergradient flux the solve adds through the face, fac rho*alpha gamma
        if (use_mrf_countergradient) {
            const int gam_comp = (n == RhoTheta_comp) ? EddyDiff::HGAMT_v : EddyDiff::HGAMQ_v;
            G += implicit_fac * rhoAlpha * myhalf * (mu_turb(i,j,k,gam_comp) + mu_turb(i,j,k-1,gam_comp));
        }
        zflux(i,j,k) = G / (mf_mx(i,j,0) * mf_my(i,j,0));
    });
}

Here is the call graph for this function: