ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Macros | Functions
ERF_ImplicitDiff_T.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_TerrainImplicitMetric.H"
#include "ERF_NativeScalarDiffusion.H"
#include "ERF_SetupVertDiff.H"
Include dependency graph for ERF_ImplicitDiff_T.cpp:

Macros

#define INSTANTIATE_IMPLICIT_DIFF_FOR_MOM_LU(STAGDIR)
 

Functions

void ImplicitDiffForStateLU_T (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 Array4< const Real > &z_nd, const Array4< const Real > &detJ, const GpuArray< Real, AMREX_SPACEDIM > &cellSizeInv, 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, const Array4< const Real > &mf_mx, const Array4< const Real > &mf_my, const bool implicit_metric)
 
template<int stagdir>
void ImplicitDiffForMomLU_T (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 Array4< const Real > &z_nd, const Array4< const Real > &detJ, const GpuArray< Real, AMREX_SPACEDIM > &cellSizeInv, 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, const Array4< const Real > &mf_x, const Array4< const Real > &mf_y, const bool implicit_metric)
 
void ImplicitVerticalFluxForReflux_T (const Box &bx, const Box &domain, const int level, const int n, const Real implicit_fac, const Array4< const Real > &cell_data, const Array4< const Real > &z_nd, const GpuArray< Real, AMREX_SPACEDIM > &cellSizeInv, const Array4< const Real > &mu_turb, const SolverChoice &solverChoice, const Array4< const Real > &mf_mx, const Array4< const Real > &mf_my, const bool implicit_metric, 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_T<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 Array4<const Real>&, \
const Array4<const Real>&, \
const GpuArray<Real, AMREX_SPACEDIM>&, \
const Array4<const Real>&, \
const SolverChoice&, \
const BCRec*, \
const bool, \
const Real, \
const bool, \
const Array4<const Real>&, \
const Array4<const Real>&, \
const bool);
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
SolverChoice
Definition: ERF_DataStruct.H:715

Function Documentation

◆ ImplicitDiffForMomLU_T()

template<int stagdir>
void ImplicitDiffForMomLU_T ( 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 Array4< const Real > &  z_nd,
const Array4< const Real > &  detJ,
const GpuArray< Real, AMREX_SPACEDIM > &  cellSizeInv,
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,
const Array4< const Real > &  mf_x,
const Array4< const Real > &  mf_y,
const bool  implicit_metric 
)

Function for computing the implicit contribution to the vertical diffusion of momentum, over terrain.

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

Parameters
[in]bxcell-centered box to loop over
[in]levelAMR level
[in]dt_dtime step
[in]cell_dataconserved cell-centered rho
[in,out]face_dataconserved momentum
[in]taustress contribution to momentum
[in]tau_corrstress contribution to momentum that will be corrected by the implicit solve
[in]z_ndnodal array of z
[in]detJJacobian determinant
[in]cellSizeInvinverse cell size array
[in]mu_turbturbulent viscosity
[in]solverChoicecontainer of parameters
[in]bc_ptrcontainer with boundary conditions
[in]use_SurfLayerwhether we have turned on subgrid diffusion
[in]implicit_facif 1 then fully implicit; if 0 then fully explicit
[in]use_ysu_mom_countergradientwhether to include YSU momentum countergradient correction
[in]mf_xx map factor at the faces of this component (u: mf_ux, v: mf_vx)
[in]mf_yy map factor at the faces of this component (u: mf_uy, v: mf_vy)
[in]implicit_metricerf.implicit_terrain_metric: add the compact terrain-metric coefficient K_h*M on the interior faces (u and v only)
499 {
500  BL_PROFILE_VAR("ImplicitDiffForMom_T()",ImplicitDiffForMom_T);
501 
502  Real dt = static_cast<Real>(dt_d);
503 
504  // setup quantities for getRhoAlphaAtFaces()
505  DiffChoice dc = solverChoice.diffChoice;
506  TurbChoice tc = solverChoice.turbChoice[level];
507  bool l_consA = (dc.molec_diff_type == MolecDiffType::ConstantAlpha);
508  bool l_turb = tc.use_kturb;
509  // The correction strains for u/v contain a factor of 1/2, so mu_eff does not carry
510  // the factor of 2 of the stresses
511  Real mu_eff = (l_consA) ? dc.dynamic_viscosity / dc.rho0_trans
512  : dc.dynamic_viscosity;
513 
514  // offsets used to average to faces
515  constexpr int ioff = (stagdir == 0) ? 1 : 0;
516  constexpr int joff = (stagdir == 1) ? 1 : 0;
517 
518  // Box bounds
519  int ilo = bx.smallEnd(0);
520  int ihi = bx.bigEnd(0);
521  int jlo = bx.smallEnd(1);
522  int jhi = bx.bigEnd(1);
523  int box_klo = bx.smallEnd(2);
524  int box_khi = bx.bigEnd(2);
525  amrex::ignore_unused(ilo, ihi, jlo, jhi);
526 
527  // The vertical extent of the domain. A column whose top or bottom is not the domain
528  // boundary ends at a coarse/fine boundary instead, and must not be given a physical BC row.
529  const int dom_klo = domain.smallEnd(2);
530  const int dom_khi = domain.bigEnd(2);
531 
532  // Temporary FABs for tridiagonal solve (allocated on column)
533  // A[k] * x[k-1] + B[k] * x[k] + C[k+1] = RHS[k]
534  amrex::FArrayBox RHS_fab, soln_fab, coeffG_fab;
535  RHS_fab.resize(bx,1, amrex::The_Async_Arena());
536  soln_fab.resize(bx,1, amrex::The_Async_Arena());
537  coeffG_fab.resize(bx,1, amrex::The_Async_Arena());
538  auto const& RHS_a = RHS_fab.array();
539  auto const& soln_a = soln_fab.array();
540  auto const& coeffG_a = coeffG_fab.array();
541 
542  Real dz_inv = cellSizeInv[2];
543 
544  int bc_comp = BCVars::xvel_bc + stagdir;
545  bool ext_dir_on_zlo = (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir ||
546  bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim);
547  bool ext_dir_on_zhi = (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir ||
548  bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim);
549  bool foextrap_on_zlo = (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap);
550  bool foextrap_on_zhi = (bc_ptr[bc_comp].hi(2) == ERFBCType::foextrap);
551  amrex::ignore_unused(foextrap_on_zlo,foextrap_on_zhi);
552 
553  AMREX_ASSERT_WITH_MESSAGE(foextrap_on_zlo || ext_dir_on_zlo || use_SurfLayer,
554  "Unexpected lower BC for momentum used with implicit vertical diffusion");
555  AMREX_ASSERT_WITH_MESSAGE(foextrap_on_zhi || ext_dir_on_zhi,
556  "Unexpected upper BC for momentum used with implicit vertical diffusion");
557 
558  Real Fact = implicit_fac * dt * dz_inv;
559 
560 #ifdef AMREX_USE_GPU
561  ParallelFor(makeSlab(bx,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int)
562  {
563 #else
564  for (int j(jlo); j<=jhi; ++j) {
565  for (int i(ilo); i<=ihi; ++i) {
566 #endif
567 
568  // The face column at (i,j) is shared by the two cell columns on either side of it,
569  // and so by the box on either side of a grid seam. Solve only over the vertical
570  // range that both of those cell columns cover: the boxes on the two sides then
571  // build the same tridiagonal system and get the same answer, so the duplicated
572  // faces along the seam stay consistent. Faces outside that range sit on a lateral
573  // coarse/fine boundary and are set from the coarse level, so it is consistent to
574  // leave them to the explicit update.
575  const int klo = amrex::max(box_klo, col_kext(i-ioff,j-joff,0,0),
576  col_kext(i ,j ,0,0));
577  const int khi = amrex::min(box_khi, col_kext(i-ioff,j-joff,0,1),
578  col_kext(i ,j ,0,1));
579 
580  // Fewer than two faces in common -- nothing to invert, leave the column explicit
581 #ifdef AMREX_USE_GPU
582  if (khi <= klo) { return; }
583 #else
584  if (khi <= klo) { continue; }
585 #endif
586  // Notes:
587  //
588  // - In DiffusionSrcForMom (e.g., for x-mom)
589  //
590  // Real diffContrib = ...
591  // + (tau13(i,j,k+1) - tau13(i,j,k)) / dzinv
592  // rho_u_rhs(i,j,k) -= diffContrib; // note the negative sign
593  //
594  // - We need to scale the explicit _part_ of `tau13` (for x-mom) by (1 - implicit_fac)
595  // The part that needs to be scaled is stored in `tau_corr`.
596  // E.g., tau13 = myhalf * (du/dz + dw/dx)
597  // tau13_corr = myhalf * du/dz
598  //
599  // - The momentum (`face_data`) was set to `S_old + S_rhs * dt`
600  // prior to including "ERF_Implicit.H". Recall that S_rhs includes
601  // sources from advection and other forcings, not just diffusion.
602  //
603  // - To correct momentum, we need to subtract `implicit_fac * diffContrib_corr`
604  // from S_rhs to recover `(1 - implicit_fac) * diffContrib_corr`,
605  // where `diffContrib_corr = -d(tau_corr)/dz`. The negative sign
606  // comes from our convention for the RHS diffusion source.
607  //
608  // Subtracting a negative gives the += below; multiply by dt to
609  // get the intermediate momentum on the RHS of the tridiagonal
610  // system.
611  //
612  // - With a surface_layer BC, tau13/23 holds the vertical flux -d_z(k*u_i)
613  // directly. We must use tau at klo (not tau_corr) with SL BCs.
614  //
615  // - The detJ for divergence was multiplied through.
616  // Therefore, it doesn't show up in the A/B denominator,
617  // but it does modify B and the RHS.
618  //
619  // - Finally, the terms ~ RHS += (tau_corr_hi - tau_corr_lo) / dz (below)
620  // essentially undo the explicit diffusion update that will be
621  // handled here implicitly.
622 
623  // Bottom boundary coefficients and RHS for L decomp
624  //===================================================
625  Real rhoface, rhoAlpha_lo, rhoAlpha_hi;
626  Real detJface, met_h_zeta_lo, met_h_zeta_hi;
627  Real a_tmp, b_tmp, c_tmp, inv_b2_tmp;
628  {
629  detJface = myhalf * (detJ(i,j,klo) + detJ(i-ioff,j-joff,klo));
630  rhoface = myhalf * (cell_data(i,j,klo,Rho_comp) + cell_data(i-ioff,j-joff,klo,Rho_comp));
631  getRhoAlphaForFaces(i, j, klo, ioff, joff, rhoAlpha_lo, rhoAlpha_hi,
632  cell_data, mu_turb, mu_eff,
633  l_consA, l_turb);
634  Real metric_lo = zero, metric_hi = zero;
635  if (implicit_metric) {
636  TerrainMetricMomCoeffs<stagdir>(i, j, klo, dom_klo, dom_khi, metric_lo, metric_hi,
637  cell_data, mu_turb, mu_eff, l_consA, l_turb,
638  z_nd, cellSizeInv, mf_x, mf_y);
639  }
640 
641  met_h_zeta_lo = MomFaceHZeta<stagdir>(i,j,klo,cellSizeInv,z_nd);
642  met_h_zeta_hi = MomFaceHZeta<stagdir>(i,j,klo+1,cellSizeInv,z_nd);
643 
644  a_tmp = zero;
645  c_tmp = -Fact * (rhoAlpha_hi + metric_hi) * dz_inv / met_h_zeta_hi;
646 
647  RHS_a(i,j,klo) = detJface * face_data(i,j,klo); // NOTE: this is momenta; solution is velocity
648 
649  // BCs: Dirichlet (u_i = val), slip wall (w = 0), or surface layer (w = 0) --
650  // but only where the bottom of this column really is the domain boundary
651  if (klo != dom_klo) {
652  // Coarse/fine (or grid) boundary: close the system with the known face value
653  // below instead of imposing a physical BC at an interior height
654  a_tmp = -Fact * (rhoAlpha_lo + metric_lo) * dz_inv / met_h_zeta_lo;
655  RHS_a(i,j,klo) += Fact * (tau_corr(i,j,klo+1) - tau_corr(i,j,klo));
656  const Real rho_below = myhalf * ( cell_data(i ,j ,klo-1,Rho_comp)
657  + cell_data(i-ioff,j-joff,klo-1,Rho_comp) );
658  RHS_a(i,j,klo) -= a_tmp * (face_data(i,j,klo-1) / rho_below);
659  } else if (ext_dir_on_zlo) {
660  RHS_a(i,j,klo) += Fact * (tau_corr(i,j,klo+1) - tau_corr(i,j,klo));
661  // Wall stress with the one-sided gradient of ComputeStrain_T,
662  // (c_wall u_wall + c_0 u_0 + c_1 u_1) / dz0 with the physical heights
663  // of the first two cells: a_tmp multiplies the wall value, and c_1 adds
664  // to the coupling with the cell above
665  const Real dz0 = MomFaceCellHeight<stagdir>(i, j, klo , z_nd);
666  const Real dz1 = MomFaceCellHeight<stagdir>(i, j, klo+1, z_nd);
667  const WallGradCoeffs wc = WallGradientCoeffs(dz0, dz1);
668  const Real idz0 = one / dz0;
669  a_tmp = Fact * rhoAlpha_lo * idz0 * wc.c_wall;
670  c_tmp += Fact * rhoAlpha_lo * idz0 * wc.c_1;
671  const Real rho_wall = myhalf * ( cell_data(i ,j ,klo-1,Rho_comp)
672  + cell_data(i-ioff,j-joff,klo-1,Rho_comp) );
673  const Real wall_velocity = face_data(i,j,klo-1) / rho_wall;
674  RHS_a(i,j,klo) -= a_tmp * wall_velocity;
675  } else if (use_SurfLayer) {
676  // NOTE: tau = -mu*d_z(u_i) w/ SL
677  RHS_a(i,j,klo) += Fact * (tau_corr(i,j,klo+1) - tau(i,j,klo));
678  RHS_a(i,j,klo) += Fact * tau(i,j,klo);
679  } else {
680  // NOTE: FOEXTRAP has zero lower flux (nothing to add to RHS)
681  RHS_a(i,j,klo) += Fact * (tau_corr(i,j,klo+1) - tau_corr(i,j,klo));
682  }
683 
684  // Add YSU momentum countergradient correction at bottom boundary.
685  // NOTE: As for the scalars, the lower face at the domain bottom carries no
686  // countergradient flux -- the surface stress is supplied by the
687  // surface layer model or the wall BC above. At a coarse/fine boundary
688  // the lower face is an interior face and carries it, as in the
689  // interior rows (erf-model/ERF#4822).
690  // NOTE: The sign matches the scalar path: tau_i3 = -rho*K*(du_i/dz - gamma_i),
691  // so the countergradient piece of the flux is +rho*K*gamma_i and its
692  // divergence enters the RHS with a minus sign.
693  if (use_ysu_mom_countergradient) {
694  const int hgam_comp = (stagdir == 0) ? EddyDiff::HGAMU_v : EddyDiff::HGAMV_v;
695  // Average HGAM* to the staggered face
696  const Real gam_klo = myhalf * (mu_turb(i,j,klo ,hgam_comp) + mu_turb(i-ioff,j-joff,klo ,hgam_comp));
697  const Real gam_kp1 = myhalf * (mu_turb(i,j,klo+1,hgam_comp) + mu_turb(i-ioff,j-joff,klo+1,hgam_comp));
698  const Real gam_hi = myhalf * (gam_klo + gam_kp1);
699  RHS_a(i,j,klo) -= Fact * rhoAlpha_hi * gam_hi;
700  if (klo != dom_klo) {
701  const Real gam_km1 = myhalf * (mu_turb(i,j,klo-1,hgam_comp) + mu_turb(i-ioff,j-joff,klo-1,hgam_comp));
702  const Real gam_lo = myhalf * (gam_klo + gam_km1);
703  RHS_a(i,j,klo) += Fact * rhoAlpha_lo * gam_lo;
704  }
705  }
706 
707  b_tmp = detJface * rhoface - a_tmp - c_tmp;
708  inv_b2_tmp = one;
709 
710  RHS_a(i,j,klo) /= b_tmp; // NOTE: this is now "rho"
711  coeffG_a(i,j,klo) = c_tmp / b_tmp; // NOTE: this is now "gamma"
712  }
713 
714  // Build the coefficients and RHS for L decomp
715  //===================================================
716  for (int k(klo+1); k < khi; k++) {
717  detJface = myhalf * (detJ(i,j,k) + detJ(i-ioff,j-joff,k));
718  rhoface = myhalf * (cell_data(i,j,k,Rho_comp) + cell_data(i-ioff,j-joff,k,Rho_comp));
719  getRhoAlphaForFaces(i, j, k, ioff, joff, rhoAlpha_lo, rhoAlpha_hi,
720  cell_data, mu_turb, mu_eff,
721  l_consA, l_turb);
722  Real metric_lo = zero, metric_hi = zero;
723  if (implicit_metric) {
724  TerrainMetricMomCoeffs<stagdir>(i, j, k, dom_klo, dom_khi, metric_lo, metric_hi,
725  cell_data, mu_turb, mu_eff, l_consA, l_turb,
726  z_nd, cellSizeInv, mf_x, mf_y);
727  }
728 
729  met_h_zeta_lo = MomFaceHZeta<stagdir>(i,j,k,cellSizeInv,z_nd);
730  met_h_zeta_hi = MomFaceHZeta<stagdir>(i,j,k+1,cellSizeInv,z_nd);
731 
732  a_tmp = -Fact * (rhoAlpha_lo + metric_lo) * dz_inv / met_h_zeta_lo;
733  c_tmp = -Fact * (rhoAlpha_hi + metric_hi) * dz_inv / met_h_zeta_hi;
734  b_tmp = detJface * rhoface - a_tmp - c_tmp;
735  inv_b2_tmp = one / (b_tmp - a_tmp * coeffG_a(i,j,k-1));
736 
737  RHS_a(i,j,k) = detJface * face_data(i,j,k); // NOTE: this is momenta; solution is velocity
738  RHS_a(i,j,k) += Fact * (tau_corr(i,j,k+1) - tau_corr(i,j,k));
739 
740  // Add YSU momentum countergradient correction
741  if (use_ysu_mom_countergradient) {
742  const int hgam_comp = (stagdir == 0) ? EddyDiff::HGAMU_v : EddyDiff::HGAMV_v;
743  // Average HGAM* to the staggered face
744  const Real gam_k = myhalf * (mu_turb(i,j,k ,hgam_comp) + mu_turb(i-ioff,j-joff,k ,hgam_comp));
745  const Real gam_km1 = myhalf * (mu_turb(i,j,k-1,hgam_comp) + mu_turb(i-ioff,j-joff,k-1,hgam_comp));
746  const Real gam_kp1 = myhalf * (mu_turb(i,j,k+1,hgam_comp) + mu_turb(i-ioff,j-joff,k+1,hgam_comp));
747  const Real gam_hi = myhalf * (gam_k + gam_kp1); // at k+1/2
748  const Real gam_lo = myhalf * (gam_k + gam_km1); // at k-1/2
749  RHS_a(i,j,k) -= Fact * (rhoAlpha_hi * gam_hi - rhoAlpha_lo * gam_lo);
750  }
751 
752  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"
753  coeffG_a(i,j,k) = c_tmp * inv_b2_tmp; // NOTE: this is now "gamma"
754  } // k
755 
756  // Top boundary coefficients and RHS for L decomp
757  //===================================================
758  {
759  detJface = myhalf * (detJ(i,j,khi) + detJ(i-ioff,j-joff,khi));
760  rhoface = myhalf * (cell_data(i,j,khi,Rho_comp) + cell_data(i-ioff,j-joff,khi,Rho_comp));
761  getRhoAlphaForFaces(i, j, khi, ioff, joff, rhoAlpha_lo, rhoAlpha_hi,
762  cell_data, mu_turb, mu_eff,
763  l_consA, l_turb);
764  Real metric_lo = zero, metric_hi = zero;
765  if (implicit_metric) {
766  TerrainMetricMomCoeffs<stagdir>(i, j, khi, dom_klo, dom_khi, metric_lo, metric_hi,
767  cell_data, mu_turb, mu_eff, l_consA, l_turb,
768  z_nd, cellSizeInv, mf_x, mf_y);
769  }
770 
771  met_h_zeta_lo = MomFaceHZeta<stagdir>(i,j,khi,cellSizeInv,z_nd);
772  met_h_zeta_hi = MomFaceHZeta<stagdir>(i,j,khi+1,cellSizeInv,z_nd);
773 
774  a_tmp = -Fact * (rhoAlpha_lo + metric_lo) * dz_inv / met_h_zeta_lo;
775  c_tmp = zero;
776 
777  RHS_a(i,j,khi) = detJface * face_data(i,j,khi); // NOTE: this is momenta; solution is velocity
778  RHS_a(i,j,khi) += Fact * (tau_corr(i,j,khi+1) - tau_corr(i,j,khi));
779 
780  // BCs: Dirichlet (u_i = val), slip wall (w = 0) --
781  // but only where the top of this column really is the domain boundary
782  if (khi != dom_khi) {
783  // Coarse/fine (or grid) boundary: close the system with the known face value
784  // above instead of imposing a physical BC at an interior height
785  c_tmp = -Fact * (rhoAlpha_hi + metric_hi) * dz_inv / met_h_zeta_hi;
786  const Real rho_above = myhalf * ( cell_data(i ,j ,khi+1,Rho_comp)
787  + cell_data(i-ioff,j-joff,khi+1,Rho_comp) );
788  RHS_a(i,j,khi) -= c_tmp * (face_data(i,j,khi+1) / rho_above);
789  } else if (ext_dir_on_zhi) {
790  // Wall stress with the one-sided gradient of ComputeStrain_T: c_tmp
791  // multiplies the wall value, and c_1 adds to the coupling with the
792  // cell below
793  const Real dz0 = MomFaceCellHeight<stagdir>(i, j, khi , z_nd);
794  const Real dz1 = MomFaceCellHeight<stagdir>(i, j, khi-1, z_nd);
795  const WallGradCoeffs wc = WallGradientCoeffs(dz0, dz1);
796  const Real idz0 = one / dz0;
797  c_tmp = Fact * rhoAlpha_hi * idz0 * wc.c_wall;
798  a_tmp += Fact * rhoAlpha_hi * idz0 * wc.c_1;
799  const Real rho_wall = myhalf * ( cell_data(i ,j ,khi+1,Rho_comp)
800  + cell_data(i-ioff,j-joff,khi+1,Rho_comp) );
801  const Real wall_velocity = face_data(i,j,khi+1) / rho_wall;
802  RHS_a(i,j,khi) -= c_tmp * wall_velocity;
803  }
804 
805  // Add YSU momentum countergradient correction at top boundary.
806  // NOTE: The lower face at khi is the interior face shared with row khi-1, which
807  // has already taken this flux out of the cell below, so it always carries
808  // it; without it here the flux left the column (erf-model/ERF#4822). The
809  // upper face carries it only at a coarse/fine boundary -- the domain top
810  // is a lid for gamma, as for the scalars.
811  if (use_ysu_mom_countergradient) {
812  const int hgam_comp = (stagdir == 0) ? EddyDiff::HGAMU_v : EddyDiff::HGAMV_v;
813  // Average HGAM* to the staggered face
814  const Real gam_khi = myhalf * (mu_turb(i,j,khi ,hgam_comp) + mu_turb(i-ioff,j-joff,khi ,hgam_comp));
815  const Real gam_km1 = myhalf * (mu_turb(i,j,khi-1,hgam_comp) + mu_turb(i-ioff,j-joff,khi-1,hgam_comp));
816  const Real gam_lo = myhalf * (gam_khi + gam_km1);
817  RHS_a(i,j,khi) += Fact * rhoAlpha_lo * gam_lo;
818  if (khi != dom_khi) {
819  const Real gam_kp1 = myhalf * (mu_turb(i,j,khi+1,hgam_comp) + mu_turb(i-ioff,j-joff,khi+1,hgam_comp));
820  const Real gam_hi = myhalf * (gam_khi + gam_kp1);
821  RHS_a(i,j,khi) -= Fact * rhoAlpha_hi * gam_hi;
822  }
823  }
824 
825  b_tmp = detJface * rhoface - a_tmp - c_tmp;
826  inv_b2_tmp = one / (b_tmp - a_tmp * coeffG_a(i,j,khi-1));
827 
828  // First solve
829  soln_a(i,j,khi) = (RHS_a(i,j,khi) - a_tmp * RHS_a(i,j,khi-1)) * inv_b2_tmp;
830  }
831 
832  // Back sweep the U decomp solution
833  //===================================================
834  for (int k(khi-1); k>=klo; --k) {
835  soln_a(i,j,k) = RHS_a(i,j,k) - coeffG_a(i,j,k) * soln_a(i,j,k+1);
836  }
837 
838  // Convert back to momenta
839  //===================================================
840  for (int k(klo); k<=khi; ++k) {
841  rhoface = myhalf * (cell_data(i,j,k,Rho_comp) + cell_data(i-ioff,j-joff,k,Rho_comp));
842  face_data(i,j,k) = rhoface * soln_a(i,j,k);
843  }
844 
845 #ifdef AMREX_USE_GPU
846  });
847 #else
848  } // i
849  } // j
850 #endif
851 }
MolecDiffType::ConstantAlpha
@ ConstantAlpha
getRhoAlphaForFaces
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void getRhoAlphaForFaces(int i, int j, int k, int ioff, int joff, amrex::Real &rhoAlpha_lo, amrex::Real &rhoAlpha_hi, const amrex::Array4< const amrex::Real > &cell_data, const amrex::Array4< const amrex::Real > &mu_turb, const amrex::Real mu_eff, bool l_consA, bool l_turb, const int eddy_comp=EddyDiff::Mom_v)
Definition: ERF_GetRhoAlphaForFaces.H:23
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:53
khi
const int khi
Definition: ERF_InitCustomPert_Bubble.H:25
klo
const int klo
Definition: ERF_InitCustomPert_KE.H:9
ParallelFor
ParallelFor(fab_box, [=] AMREX_GPU_DEVICE(int i, int j, int k) { qrcuten_arr(i, j, k)=Real(0);qscuten_arr(i, j, k)=Real(0);qicuten_arr(i, j, k)=Real(0);})
one
constexpr amrex::Real one
Definition: ERF_NumericalConstants.H:37
zero
constexpr amrex::Real zero
Definition: ERF_NumericalConstants.H:36
myhalf
constexpr amrex::Real myhalf
Definition: ERF_NumericalConstants.H:41
AMREX_ASSERT_WITH_MESSAGE
AMREX_ASSERT_WITH_MESSAGE(wbar_cutoff_min > wbar_cutoff_max, "ERROR: wbar_cutoff_min < wbar_cutoff_max")
WallGradientCoeffs
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE WallGradCoeffs WallGradientCoeffs(const amrex::Real dz0, const amrex::Real dz1) noexcept
Definition: ERF_WallGradientCoeffs.H:33
BCVars::xvel_bc
@ xvel_bc
Definition: ERF_IndexDefines.H:126
ERFBCType::foextrap
@ foextrap
Definition: ERF_IndexDefines.H:339
ERFBCType::ext_dir
@ ext_dir
Definition: ERF_IndexDefines.H:340
ERFBCType::ext_dir_prim
@ ext_dir_prim
Definition: ERF_IndexDefines.H:343
EddyDiff::HGAMU_v
@ HGAMU_v
Definition: ERF_IndexDefines.H:302
EddyDiff::HGAMV_v
@ HGAMV_v
Definition: ERF_IndexDefines.H:303
DiffChoice
Definition: ERF_DiffStruct.H:24
DiffChoice::molec_diff_type
MolecDiffType molec_diff_type
Selected molecular transport model.
Definition: ERF_DiffStruct.H:96
DiffChoice::dynamic_viscosity
amrex::Real dynamic_viscosity
Dynamic viscosity for momentum diffusion [kg/(m-s)].
Definition: ERF_DiffStruct.H:108
DiffChoice::rho0_trans
amrex::Real rho0_trans
Reference density used to compute dynamic diffusion coefficients [kg/m3].
Definition: ERF_DiffStruct.H:103
SolverChoice::turbChoice
amrex::Vector< TurbChoice > turbChoice
Turbulence options for each AMR level.
Definition: ERF_DataStruct.H:2620
SolverChoice::diffChoice
DiffChoice diffChoice
Diffusion-related options.
Definition: ERF_DataStruct.H:2617
TurbChoice
Definition: ERF_TurbStruct.H:119
TurbChoice::use_kturb
bool use_kturb
Whether any turbulence model is active.
Definition: ERF_TurbStruct.H:977
WallGradCoeffs
Definition: ERF_WallGradientCoeffs.H:25
WallGradCoeffs::c_wall
amrex::Real c_wall
Definition: ERF_WallGradientCoeffs.H:26
WallGradCoeffs::c_1
amrex::Real c_1
Definition: ERF_WallGradientCoeffs.H:28
Here is the call graph for this function:

◆ ImplicitDiffForStateLU_T()

void ImplicitDiffForStateLU_T ( 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 Array4< const Real > &  z_nd,
const Array4< const Real > &  detJ,
const GpuArray< Real, AMREX_SPACEDIM > &  cellSizeInv,
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,
const Array4< const Real > &  mf_mx,
const Array4< const Real > &  mf_my,
const bool  implicit_metric 
)

Function for computing the implicit contribution to the vertical diffusion of theta, with terrain.

Parameters
[in]bxcell-centered box to loop over
[in]domainbox of the whole domain
[in]levelAMR level
[in]nconserved component index
[in]dt_dtime step
[in]bc_neumann_valsvalues of derivatives if bc_type == Neumann
[in,out]cell_dataconserved cell-centered rho, rho theta
[in]z_ndnodal array of z
[in]detJJacobian determinant
[in]cellSizeInvinverse cell size array
[in]scalar_zfluxscalar vertical flux in z-dir
[in]mu_turbturbulent viscosity
[in]solverChoicecontainer of parameters
[in]bc_ptrcontainer with boundary conditions
[in]use_SurfLayerwhether we have turned on subgrid diffusion
[in]implicit_facif 1 then fully implicit; if 0 then fully explicit
[in]use_mrf_countergradientwhether to include MRF countergradient correction
[in]mf_mxx map factor at cell centers
[in]mf_myy map factor at cell centers
[in]implicit_metricerf.implicit_terrain_metric: add the compact terrain-metric coefficient K_h*M_s on the interior faces
99 {
100  BL_PROFILE_VAR("ImplicitDiffForState_T()",ImplicitDiffForState_T);
101 
102  Real dt = static_cast<Real>(dt_d);
103 
104  // setup quantities for getRhoAlpha()
105 #include "ERF_SetupVertDiff.H"
106  const int qty_index = n;
107  const int prim_index = qty_index - 1;
108  const int prim_scal_index = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? PrimScalar_comp : prim_index;
109 
110  // erf.implicit_terrain_metric: the horizontal eddy diffusivity of this quantity, from the
111  // selector the explicit fluxes use, so that the two sides cannot drift apart
112  const int eddy_h_comp =
113  ResolveNativeScalarDiffusionPolicy(qty_index, solverChoice.diffChoice).coefficients.eddy_h_comp;
114  const int dom_klo = domain.smallEnd(2);
115  const int dom_khi = domain.bigEnd(2);
116 
117  // Box bounds
118  int ilo = bx.smallEnd(0);
119  int ihi = bx.bigEnd(0);
120  int jlo = bx.smallEnd(1);
121  int jhi = bx.bigEnd(1);
122  int klo = bx.smallEnd(2);
123  int khi = bx.bigEnd(2);
124  amrex::ignore_unused(ilo, ihi, jlo, jhi);
125 
126  // Temporary FABs for tridiagonal solve (allocated on column)
127  // A[k] * x[k-1] + B[k] * x[k] + C[k+1] = RHS[k]
128 
129  // With LU decomposition, M * x = r is written as L * U * x = r with U * x = rho
130  // We then first have L * rho = r and U * x = rho
131  amrex::FArrayBox RHS_fab, soln_fab, coeffG_fab;
132  RHS_fab.resize(bx,1, amrex::The_Async_Arena());
133  soln_fab.resize(bx,1, amrex::The_Async_Arena());
134  coeffG_fab.resize(bx,1, amrex::The_Async_Arena());
135  auto const& RHS_a = RHS_fab.array();
136  auto const& soln_a = soln_fab.array();
137  auto const& coeffG_a = coeffG_fab.array();
138 
139  Real dz_inv = cellSizeInv[2];
140 
141  // BC types are stored by BC variable, not by state component
142  const int bc_comp = ResolveNativeScalarDiffusionPolicy(qty_index, solverChoice.diffChoice).bc_comp;
143  bool foextrap_on_zlo = (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap);
144  bool foextrap_on_zhi = (bc_ptr[bc_comp].hi(2) == ERFBCType::foextrap);
145  bool neumann_on_zlo = (bc_ptr[bc_comp].lo(2) == ERFBCType::neumann);
146  bool neumann_on_zhi = (bc_ptr[bc_comp].hi(2) == ERFBCType::neumann);
147  amrex::ignore_unused(foextrap_on_zlo, foextrap_on_zhi);
148 
149  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(foextrap_on_zlo || neumann_on_zlo || use_SurfLayer,
150  "Unexpected lower BC for scalars used with implicit vertical diffusion");
151  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(foextrap_on_zhi || neumann_on_zhi,
152  "Unexpected upper BC for scalars used with implicit vertical diffusion");
153 
154  // k-eqn RANS with a Dirichlet wall value of k: the first cell is held at
155  // the value set by the surface layer, so its row reduces to x(klo) = phi(klo)
156  // and the row above sees it as a Dirichlet neighbour.
157  // NOTE: init_bcs requires zlo.type = surface_layer whenever dirichlet_k is set,
158  // so this matches the l_dirichlet_k pin in erf_slow_rhs_post, which also
159  // requires a surface layer. Do not relax one without the other.
160  const bool pin_klo = (qty_index == RhoKE_comp) && (klo == domain.smallEnd(2)) &&
161  (solverChoice.turbChoice[level].rans_type == RANSType::kEqn) &&
162  solverChoice.turbChoice[level].dirichlet_k;
163 
164  // A column whose bottom or top is not the domain boundary ends at a coarse/fine boundary
165  // instead. There the system must be closed with the value on the other side of that
166  // boundary rather than given a physical BC row at an interior height, so that the answer
167  // does not depend on where the grid happens to end in the vertical.
168  const bool at_zlo = (klo == domain.smallEnd(2));
169  const bool at_zhi = (khi == domain.bigEnd(2));
170 
171  Real Fact = implicit_fac * dt * dz_inv;
172 
173 #ifdef AMREX_USE_GPU
174  ParallelFor(makeSlab(bx,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int)
175  {
176 #else
177  for (int j(jlo); j<=jhi; ++j) {
178  for (int i(ilo); i<=ihi; ++i) {
179 #endif
180  // Notes: The detJ for divergence was multiplied through.
181  // Therefore, it doesn't show up in the A/B denominator,
182  // but it does modify B and the RHS.
183 
184  // Bottom boundary coefficients and RHS for L decomp
185  //===================================================
186  Real rhoAlpha_lo, rhoAlpha_hi;
187  Real met_h_zeta_lo, met_h_zeta_hi;
188  Real a_tmp, b_tmp, c_tmp, inv_b2_tmp;
189  {
190  getRhoAlpha(i, j, klo, rhoAlpha_lo, rhoAlpha_hi,
191  cell_data, mu_turb, d_alpha_eff, d_eddy_diff_idz,
192  prim_scal_index, l_consA, l_turb);
193  Real metric_lo = zero, metric_hi = zero;
194  if (implicit_metric) {
195  const Real alpha_mol = d_alpha_eff[prim_scal_index];
196  TerrainMetricScalarCoeffs(i, j, klo, dom_klo, dom_khi, metric_lo, metric_hi,
197  cell_data, mu_turb, alpha_mol, eddy_h_comp,
198  l_consA, l_turb, z_nd, cellSizeInv, mf_mx, mf_my);
199  }
200 
201  met_h_zeta_hi = Compute_h_zeta_AtKface(i,j,klo+1,cellSizeInv,z_nd);
202 
203  a_tmp = zero;
204  if (!at_zlo) {
205  met_h_zeta_lo = Compute_h_zeta_AtKface(i,j,klo,cellSizeInv,z_nd);
206  a_tmp = -Fact * (rhoAlpha_lo + metric_lo) * dz_inv / met_h_zeta_lo;
207  }
208  c_tmp = (pin_klo) ? zero : -Fact * (rhoAlpha_hi + metric_hi) * dz_inv / met_h_zeta_hi;
209  b_tmp = detJ(i,j,klo) * cell_data(i,j,klo,Rho_comp) - a_tmp - c_tmp;
210  inv_b2_tmp = one;
211 
212  RHS_a(i,j,klo) = detJ(i,j,klo) * cell_data(i,j,klo,n); // NOTE: this is rho*phi; solution is phi
213  if (!at_zlo) {
214  // Coarse/fine boundary: close the system with the known value below
215  RHS_a(i,j,klo) -= a_tmp * (cell_data(i,j,klo-1,n) / cell_data(i,j,klo-1,Rho_comp));
216  } else if (pin_klo) {
217  // Dirichlet row: no flux terms
218  } else if (use_SurfLayer && scalar_zflux) {
219  RHS_a(i,j,klo) += Fact * scalar_zflux(i,j,klo); // NOTE: scalar_zflux = -K*d_z(\phi)
220  } else if (neumann_on_zlo) {
221  RHS_a(i,j,klo) += -Fact * rhoAlpha_lo * bc_neumann_vals[2]; // NOTE: N_val = d_z(\phi)
222  }
223 
224  // Add countergradient correction to RHS at bottom boundary.
225  // NOTE: The lower face at klo carries no countergradient flux -- the
226  // total surface flux is supplied by the surface layer model or the
227  // Neumann BC, while gamma represents nonlocal transport interior to
228  // the PBL. Only the upper face contributes here.
229  if (use_mrf_countergradient && (n == RhoTheta_comp || n == RhoQ1_comp)) {
230  const int gam_comp = (n == RhoTheta_comp) ? EddyDiff::HGAMT_v : EddyDiff::HGAMQ_v;
231  const Real gam_hi = myhalf * (mu_turb(i, j, klo, gam_comp) + mu_turb(i, j, klo+1, gam_comp));
232  RHS_a(i,j,klo) -= Fact * rhoAlpha_hi * gam_hi;
233  if (!at_zlo) {
234  const Real gam_lo = myhalf * (mu_turb(i, j, klo, gam_comp) + mu_turb(i, j, klo-1, gam_comp));
235  RHS_a(i,j,klo) += Fact * rhoAlpha_lo * gam_lo;
236  }
237  }
238 
239  RHS_a(i,j,klo) /= b_tmp; // NOTE: this is now "rho"
240  coeffG_a(i,j,klo) = c_tmp / b_tmp; // NOTE: this is now "gamma"
241  }
242 
243  // Build the coefficients and RHS for L decomp
244  //===================================================
245  for (int k(klo+1); k < khi; k++) {
246  getRhoAlpha(i, j, k, rhoAlpha_lo, rhoAlpha_hi,
247  cell_data, mu_turb, d_alpha_eff, d_eddy_diff_idz,
248  prim_scal_index, l_consA, l_turb);
249  Real metric_lo = zero, metric_hi = zero;
250  if (implicit_metric) {
251  const Real alpha_mol = d_alpha_eff[prim_scal_index];
252  TerrainMetricScalarCoeffs(i, j, k, dom_klo, dom_khi, metric_lo, metric_hi,
253  cell_data, mu_turb, alpha_mol, eddy_h_comp,
254  l_consA, l_turb, z_nd, cellSizeInv, mf_mx, mf_my);
255  }
256 
257  met_h_zeta_lo = Compute_h_zeta_AtKface(i,j,k ,cellSizeInv,z_nd);
258  met_h_zeta_hi = Compute_h_zeta_AtKface(i,j,k+1,cellSizeInv,z_nd);
259 
260  a_tmp = -Fact * (rhoAlpha_lo + metric_lo) * dz_inv / met_h_zeta_lo;
261  c_tmp = -Fact * (rhoAlpha_hi + metric_hi) * dz_inv / met_h_zeta_hi;
262  b_tmp = detJ(i,j,k) * cell_data(i,j,k,Rho_comp) - a_tmp - c_tmp;
263  inv_b2_tmp = one / (b_tmp - a_tmp * coeffG_a(i,j,k-1));
264 
265  RHS_a(i,j,k) = detJ(i,j,k) * cell_data(i,j,k,n); // NOTE: this is rho*phi; solution is phi
266 
267  // Add countergradient correction to RHS in interior
268  if (use_mrf_countergradient && (n == RhoTheta_comp || n == RhoQ1_comp)) {
269  const int gam_comp = (n == RhoTheta_comp) ? EddyDiff::HGAMT_v : EddyDiff::HGAMQ_v;
270  const Real gam_k = mu_turb(i, j, k, gam_comp);
271  const Real gam_km1 = mu_turb(i, j, k-1, gam_comp);
272  const Real gam_kp1 = mu_turb(i, j, k+1, gam_comp);
273  const Real gam_hi = myhalf * (gam_k + gam_kp1); // at k+½
274  const Real gam_lo = myhalf * (gam_k + gam_km1); // at k-½
275  // Countergradient flux divergence (implicit contribution to RHS):
276  // -Fact * [ρα_{k+½}·γ_{k+½} - ρα_{k-½}·γ_{k-½}]
277  // NOTE: no 1/met_h_zeta here. met_h_zeta appears in a_tmp/c_tmp only
278  // because the diffusive flux contains an inner vertical
279  // derivative d(phi)/dz ~ dz_inv/h_zeta * delta_phi. The
280  // countergradient piece contains no derivative, so it takes only
281  // the outer dz_inv already carried in Fact.
282  RHS_a(i,j,k) -= Fact * (rhoAlpha_hi * gam_hi - rhoAlpha_lo * gam_lo);
283  }
284 
285  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"
286  coeffG_a(i,j,k) = c_tmp * inv_b2_tmp; // NOTE: this is now "gamma"
287  } // k
288 
289  // Top boundary coefficients and RHS for L decomp
290  //===================================================
291  {
292  getRhoAlpha(i, j, khi, rhoAlpha_lo, rhoAlpha_hi,
293  cell_data, mu_turb, d_alpha_eff, d_eddy_diff_idz,
294  prim_scal_index, l_consA, l_turb);
295  Real metric_lo = zero, metric_hi = zero;
296  if (implicit_metric) {
297  const Real alpha_mol = d_alpha_eff[prim_scal_index];
298  TerrainMetricScalarCoeffs(i, j, khi, dom_klo, dom_khi, metric_lo, metric_hi,
299  cell_data, mu_turb, alpha_mol, eddy_h_comp,
300  l_consA, l_turb, z_nd, cellSizeInv, mf_mx, mf_my);
301  }
302 
303  // Lower-face metric shared with row khi-1.
304  met_h_zeta_lo = Compute_h_zeta_AtKface(i,j,khi ,cellSizeInv,z_nd);
305 
306  a_tmp = -Fact * (rhoAlpha_lo + metric_lo) * dz_inv / met_h_zeta_lo;
307  c_tmp = zero;
308  if (!at_zhi) {
309  met_h_zeta_hi = Compute_h_zeta_AtKface(i,j,khi+1,cellSizeInv,z_nd);
310  c_tmp = -Fact * (rhoAlpha_hi + metric_hi) * dz_inv / met_h_zeta_hi;
311  }
312  b_tmp = detJ(i,j,khi) * cell_data(i,j,khi,Rho_comp) - a_tmp - c_tmp;
313  inv_b2_tmp = one / (b_tmp - a_tmp * coeffG_a(i,j,khi-1));
314 
315  RHS_a(i,j,khi) = detJ(i,j,khi) * cell_data(i,j,khi,n); // NOTE: this is rho*phi; solution is phi
316  if (!at_zhi) {
317  // Coarse/fine boundary: close the system with the known value above
318  RHS_a(i,j,khi) -= c_tmp * (cell_data(i,j,khi+1,n) / cell_data(i,j,khi+1,Rho_comp));
319  } else if (neumann_on_zhi) {
320  RHS_a(i,j,khi) -= -Fact * rhoAlpha_hi * bc_neumann_vals[5]; // NOTE: N_val = d_z(\phi)
321  }
322 
323  // Add countergradient correction to RHS at top boundary.
324  // NOTE: The lower face at khi always carries the countergradient flux, which
325  // row khi-1 has already taken out of the cell below; without it here
326  // that flux would leave the column. The upper face carries it only at
327  // a coarse/fine boundary -- the domain top is a lid for gamma.
328  if (use_mrf_countergradient && (n == RhoTheta_comp || n == RhoQ1_comp)) {
329  const int gam_comp = (n == RhoTheta_comp) ? EddyDiff::HGAMT_v : EddyDiff::HGAMQ_v;
330  const Real gam_lo = myhalf * (mu_turb(i, j, khi, gam_comp) + mu_turb(i, j, khi-1, gam_comp));
331  RHS_a(i,j,khi) += Fact * rhoAlpha_lo * gam_lo;
332  if (!at_zhi) {
333  const Real gam_hi = myhalf * (mu_turb(i, j, khi, gam_comp) + mu_turb(i, j, khi+1, gam_comp));
334  RHS_a(i,j,khi) -= Fact * rhoAlpha_hi * gam_hi;
335  }
336  }
337 
338  // First solve
339  soln_a(i,j,khi) = (RHS_a(i,j,khi) - a_tmp * RHS_a(i,j,khi-1)) * inv_b2_tmp;
340  }
341 
342  // Back sweep the U decomp solution
343  //===================================================
344  for (int k(khi-1); k>=klo; --k) {
345  soln_a(i,j,k) = RHS_a(i,j,k) - coeffG_a(i,j,k) * soln_a(i,j,k+1);
346  }
347 
348  // Convert back to rho*theta
349  //===================================================
350  for (int k(klo); k<=khi; ++k) {
351  cell_data(i,j,k,n) = cell_data(i,j,k,Rho_comp) * soln_a(i,j,k);
352  }
353 
354 #ifdef AMREX_USE_GPU
355  });
356 #else
357  } // i
358  } // j
359 #endif
360 }
getRhoAlpha
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void getRhoAlpha(int i, int j, int k, amrex::Real &rhoAlpha_lo, amrex::Real &rhoAlpha_hi, const amrex::Array4< const amrex::Real > &cell_data, const amrex::Array4< const amrex::Real > &mu_turb, const amrex::Real *d_alpha_eff, const int *d_eddy_diff_idz, int prim_scal_index, bool l_consA, bool l_turb)
Definition: ERF_GetRhoAlpha.H:20
RhoScalar_comp
#define RhoScalar_comp
Definition: ERF_IndexDefines.H:57
RhoTheta_comp
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:54
NSCALARS
#define NSCALARS
Definition: ERF_IndexDefines.H:26
RhoQ1_comp
#define RhoQ1_comp
Definition: ERF_IndexDefines.H:66
PrimScalar_comp
#define PrimScalar_comp
Definition: ERF_IndexDefines.H:81
RhoKE_comp
#define RhoKE_comp
Definition: ERF_IndexDefines.H:55
AMREX_ALWAYS_ASSERT_WITH_MESSAGE
AMREX_ALWAYS_ASSERT_WITH_MESSAGE(axis >=0 &&axis< AMREX_SPACEDIM, "prob.axis must be 0, 1 or 2")
ResolveNativeScalarDiffusionPolicy
NativeScalarDiffusionPolicy ResolveNativeScalarDiffusionPolicy(const int qty_comp, const DiffChoice &choice)
Definition: ERF_NativeScalarDiffusion.H:26
ERF_SetupVertDiff.H
Compute_h_zeta_AtKface
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_zeta_AtKface(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:420
ERFBCType::neumann
@ neumann
Definition: ERF_IndexDefines.H:345
EddyDiff::HGAMQ_v
@ HGAMQ_v
Definition: ERF_IndexDefines.H:301
EddyDiff::HGAMT_v
@ HGAMT_v
Definition: ERF_IndexDefines.H:300
NativeScalarDiffusionPolicy::bc_comp
int bc_comp
Definition: ERF_NativeScalarDiffusion.H:16
NativeScalarDiffusionPolicy::coefficients
ScalarDiffusionCoefficients coefficients
Definition: ERF_NativeScalarDiffusion.H:12
ScalarDiffusionCoefficients::eddy_h_comp
int eddy_h_comp
Definition: ERF_ScalarDiffusion.H:25
Here is the call graph for this function:

◆ ImplicitVerticalFluxForReflux_T()

void ImplicitVerticalFluxForReflux_T ( const Box &  bx,
const Box &  domain,
const int  level,
const int  n,
const Real  implicit_fac,
const Array4< const Real > &  cell_data,
const Array4< const Real > &  z_nd,
const GpuArray< Real, AMREX_SPACEDIM > &  cellSizeInv,
const Array4< const Real > &  mu_turb,
const SolverChoice &  solverChoice,
const Array4< const Real > &  mf_mx,
const Array4< const Real > &  mf_my,
const bool  implicit_metric,
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_T solve applied it, for refluxing: in the convention of the advective fluxes the flux registers hold (see AddScalarDiffusionRefluxFluxes_T), 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]bxcell-centered box
[in]domainbox of the whole domain
[in]levelAMR level
[in]nconserved component index
[in]implicit_facthe vert_implicit_fac the solve used
[in]cell_dataconserved state after the solve
[in]use_mrf_countergradientwhether the solve added the MRF countergradient flux
[out]zfluxface flux of component n, in component 0
915 {
916 #include "ERF_SetupVertDiff.H"
917  const int qty_index = n;
918  const int prim_index = qty_index - 1;
919  const int prim_scal_index = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? PrimScalar_comp : prim_index;
920  const int eddy_h_comp =
921  ResolveNativeScalarDiffusionPolicy(qty_index, solverChoice.diffChoice).coefficients.eddy_h_comp;
922  const int dom_klo = domain.smallEnd(2);
923  const int dom_khi = domain.bigEnd(2);
924  const int khi = bx.bigEnd(2);
925  const Real dz_inv = cellSizeInv[2];
926 
927  ParallelFor(surroundingNodes(bx,2), [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
928  {
929  if (k == dom_lo.z || k == dom_hi.z+1) { zflux(i,j,k) = zero; return; }
930  // The coefficient on this face, from the cell above it (or below, at the top of bx)
931  Real rhoAlpha_lo, rhoAlpha_hi;
932  const int kc = (k <= khi) ? k : k-1;
933  getRhoAlpha(i, j, kc, rhoAlpha_lo, rhoAlpha_hi, cell_data, mu_turb, d_alpha_eff, d_eddy_diff_idz,
934  prim_scal_index, l_consA, l_turb);
935  Real metric_lo = zero, metric_hi = zero;
936  if (implicit_metric) {
937  const Real alpha_mol = d_alpha_eff[prim_scal_index];
938  TerrainMetricScalarCoeffs(i, j, kc, dom_klo, dom_khi, metric_lo, metric_hi,
939  cell_data, mu_turb, alpha_mol, eddy_h_comp,
940  l_consA, l_turb, z_nd, cellSizeInv, mf_mx, mf_my);
941  }
942  const Real rhoAlpha = (k <= khi) ? rhoAlpha_lo : rhoAlpha_hi;
943  const Real coef = rhoAlpha + ((k <= khi) ? metric_lo : metric_hi);
944  const Real h_zeta = Compute_h_zeta_AtKface(i, j, k, cellSizeInv, z_nd);
945  const Real phi_hi = cell_data(i,j,k ,n) / cell_data(i,j,k ,Rho_comp);
946  const Real phi_lo = cell_data(i,j,k-1,n) / cell_data(i,j,k-1,Rho_comp);
947  Real G = -implicit_fac * coef * dz_inv / h_zeta * (phi_hi - phi_lo);
948  // The countergradient flux the solve adds through the face, fac rho*alpha gamma
949  if (use_mrf_countergradient) {
950  const int gam_comp = (n == RhoTheta_comp) ? EddyDiff::HGAMT_v : EddyDiff::HGAMQ_v;
951  G += implicit_fac * rhoAlpha * myhalf * (mu_turb(i,j,k,gam_comp) + mu_turb(i,j,k-1,gam_comp));
952  }
953  zflux(i,j,k) = G / (mf_mx(i,j,0) * mf_my(i,j,0));
954  });
955 }
Here is the call graph for this function: