ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Functions
ERF_SolveWithEBMLMG.cpp File Reference
#include "ERF.H"
#include "ERF_EB.H"
#include "ERF_Utils.H"
#include "ERF_SolverUtils.H"
#include <AMReX_MLMG.H>
#include <AMReX_MLEBABecLap.H>
Include dependency graph for ERF_SolveWithEBMLMG.cpp:

Functions

void FillZeroAreaFaceFluxes (MultiFab &phi, Array< MultiFab, AMREX_SPACEDIM > &fluxes, const Geometry &geom, EBFArrayBoxFactory const &ebfact, eb_aux_ const &ebfact_u, eb_aux_ const &ebfact_v, eb_aux_ const &ebfact_w)
 
void solve_with_EB_mlmg (int lev, Vector< MultiFab > &rhs, Vector< MultiFab > &phi, Vector< Array< MultiFab, AMREX_SPACEDIM >> &fluxes, EBFArrayBoxFactory const &ebfact, eb_aux_ const &ebfact_u, eb_aux_ const &ebfact_v, eb_aux_ const &ebfact_w, const Geometry &geom, const Vector< amrex::IntVect > &ref_ratio, Array< std::string, 2 *AMREX_SPACEDIM > domain_bc_type, int mg_verbose, Real reltol, Real abstol)
 

Function Documentation

◆ FillZeroAreaFaceFluxes()

void FillZeroAreaFaceFluxes ( MultiFab &  phi,
Array< MultiFab, AMREX_SPACEDIM > &  fluxes,
const Geometry &  geom,
EBFArrayBoxFactory const &  ebfact,
eb_aux_ const &  ebfact_u,
eb_aux_ const &  ebfact_v,
eb_aux_ const &  ebfact_w 
)

Compute phi gradients where the area of the face is zero

13 {
14  BL_PROFILE("ERF::FillZeroAreaFaceFluxes()");
15 
16  const GpuArray<Real, AMREX_SPACEDIM> dxInv = geom.InvCellSizeArray();
17 
18  for (MFIter mfi(phi,TileNoZ()); mfi.isValid(); ++mfi)
19  {
20  const Box& tbx = mfi.tilebox();
21  const Box& xbx = mfi.nodaltilebox(0);
22  const Box& ybx = mfi.nodaltilebox(1);
23  const Box& zbx = mfi.nodaltilebox(2);
24 
25  // EBCellFlagFab const& cflag_fab = (ebfact.get_const_factory())->getMultiEBCellFlagFab()[mfi];
26  EBCellFlagFab const& cflag_fab = (ebfact.getMultiEBCellFlagFab())[mfi];
27  Array4<const EBCellFlag> cflag = cflag_fab.const_array();
28 
29  if (cflag_fab.getType(tbx) == FabType::singlevalued)
30  {
31  Array4<const Real> apx = ebfact.getAreaFrac()[0]->const_array(mfi);
32  Array4<const Real> apy = ebfact.getAreaFrac()[1]->const_array(mfi);
33  Array4<const Real> apz = ebfact.getAreaFrac()[2]->const_array(mfi);
34 
35  Array4<const EBCellFlag> u_cflag = ebfact_u.getMultiEBCellFlagFab()[mfi].const_array();
36  Array4<const EBCellFlag> v_cflag = ebfact_v.getMultiEBCellFlagFab()[mfi].const_array();
37  Array4<const EBCellFlag> w_cflag = ebfact_w.getMultiEBCellFlagFab()[mfi].const_array();
38 
39  Array4<Real const> const& p_arr = phi.const_array(mfi);
40  Array4<Real> const& fx = fluxes[0].array(mfi);
41  Array4<Real> const& fy = fluxes[1].array(mfi);
42  Array4<Real> const& fz = fluxes[2].array(mfi);
43 
44  ParallelFor(xbx, ybx, zbx,
45  // x-face
46  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
47  {
48  if (apx(i,j,k) == zero) {
49  if (!u_cflag(i,j,k).isCovered()) {
50  if (cflag(i,j,k).isCovered() && !cflag(i-1,j,k).isCovered()) {
51  fx(i,j,k) = dxInv[0] * (p_arr(i-3,j,k) - three*p_arr(i-2,j,k) + two*p_arr(i-1,j,k));
52  } else if (cflag(i-1,j,k).isCovered() && !cflag(i,j,k).isCovered()) {
53  fx(i,j,k) = dxInv[0] * (three*p_arr(i+1,j,k) - p_arr(i+2,j,k) - two*p_arr(i,j,k));
54  }
55  }
56  }
57  },
58  // y-face
59  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
60  {
61  if (apy(i,j,k) == zero) {
62  if (!v_cflag(i,j,k).isCovered()) {
63  if (cflag(i,j,k).isCovered() && !cflag(i,j-1,k).isCovered()) {
64  fy(i,j,k) = dxInv[1] * (p_arr(i,j-3,k) - three*p_arr(i,j-2,k) + two*p_arr(i,j-1,k));
65  } else if (cflag(i,j-1,k).isCovered() && !cflag(i,j,k).isCovered()) {
66  fy(i,j,k) = dxInv[1] * (three*p_arr(i,j+1,k) - p_arr(i,j+2,k) - two*p_arr(i,j,k));
67  }
68  }
69  }
70  },
71  // z-face
72  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
73  {
74  if (apz(i,j,k) == zero) {
75  if (!w_cflag(i,j,k).isCovered()) {
76  if (cflag(i,j,k).isCovered() && !cflag(i,j,k-1).isCovered()) {
77  fz(i,j,k) = dxInv[2] * (p_arr(i,j,k-3) - three*p_arr(i,j,k-2) + two*p_arr(i,j,k-1));
78  } else if (cflag(i,j,k-1).isCovered() && !cflag(i,j,k).isCovered()) {
79  fz(i,j,k) = dxInv[2] * (three*p_arr(i,j,k+1) - p_arr(i,j,k+2) - two*p_arr(i,j,k));
80  }
81  }
82  }
83  });
84  } // single-valued
85  } // mfi
86 }
three
constexpr amrex::Real three
Definition: ERF_Constants.H:9
two
constexpr amrex::Real two
Definition: ERF_Constants.H:8
zero
constexpr amrex::Real zero
Definition: ERF_Constants.H:6
dxInv
amrex::GpuArray< Real, AMREX_SPACEDIM > dxInv
Definition: ERF_InitCustomPertVels_ParticleTests.H:17
ParallelFor
ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { const Real *dx=geomdata.CellSize();const Real x=(i+0.5) *dx[0];const Real y=(j+0.5) *dx[1];const Real Omg=erf_vortex_Gaussian(x, y, xc, yc, R, beta, sigma);const Real deltaT=-(gamma - 1.0)/(2.0 *sigma *sigma) *Omg *Omg;const Real rho_norm=std::pow(1.0+deltaT, inv_gm1);const Real T=(1.0+deltaT) *T_inf;const Real p=std::pow(rho_norm, Gamma)/Gamma *rho_0 *a_inf *a_inf;const Real rho_theta=rho_0 *rho_norm *(T *std::pow(p_0/p, rdOcp));state_pert(i, j, k, RhoTheta_comp)=rho_theta - getRhoThetagivenP(p_hse(i, j, k));const Real r2d_xy=std::sqrt((x-xc) *(x-xc)+(y-yc) *(y-yc));state_pert(i, j, k, RhoScalar_comp)=0.25 *(1.0+std::cos(PI *std::min(r2d_xy, R)/R));})
zbx
const Box zbx
Definition: ERF_SetupDiff.H:9
xbx
const Box xbx
Definition: ERF_SetupDiff.H:7
ybx
const Box ybx
Definition: ERF_SetupDiff.H:8
TileNoZ
AMREX_FORCE_INLINE amrex::IntVect TileNoZ()
Definition: ERF_TileNoZ.H:11

Referenced by solve_with_EB_mlmg().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ solve_with_EB_mlmg()

void solve_with_EB_mlmg ( int  lev,
Vector< MultiFab > &  rhs,
Vector< MultiFab > &  phi,
Vector< Array< MultiFab, AMREX_SPACEDIM >> &  fluxes,
EBFArrayBoxFactory const &  ebfact,
eb_aux_ const &  ebfact_u,
eb_aux_ const &  ebfact_v,
eb_aux_ const &  ebfact_w,
const Geometry &  geom,
const Vector< amrex::IntVect > &  ref_ratio,
Array< std::string, 2 *AMREX_SPACEDIM >  domain_bc_type,
int  mg_verbose,
Real  reltol,
Real  abstol 
)

Solve the Poisson equation using EB_enabled MLMG Note that the level may or may not be level zero

Important: we solve on the whole level even if there are disjoint regions

38 {
39  BL_PROFILE("ERF::solve_with_EB_mlmg()");
40 
41  auto const dom_lo = lbound(geom.Domain());
42  auto const dom_hi = ubound(geom.Domain());
43 
44  LPInfo info;
45  // Allow a hidden direction if the domain is one cell wide in any lateral direction
46  if (dom_lo.x == dom_hi.x) {
47  info.setHiddenDirection(0);
48  } else if (dom_lo.y == dom_hi.y) {
49  info.setHiddenDirection(1);
50  }
51 
52  // Make sure the solver only sees the levels over which we are solving
53  Vector<BoxArray> ba_tmp; ba_tmp.push_back(rhs[0].boxArray());
54  Vector<DistributionMapping> dm_tmp; dm_tmp.push_back(rhs[0].DistributionMap());
55  Vector<Geometry> geom_tmp; geom_tmp.push_back(geom);
56 
57  auto bclo = get_lo_projection_bc(geom,domain_bc_type);
58  auto bchi = get_hi_projection_bc(geom,domain_bc_type);
59 
60  amrex::Print() << "BCLO " << bclo[0] << " " << bclo[1] << " " << bclo[2] << std::endl;
61  amrex::Print() << "BCHI " << bchi[0] << " " << bchi[1] << " " << bchi[2] << std::endl;
62 
63  // ****************************************************************************
64  // Multigrid solve
65  // ****************************************************************************
66 
67  MLEBABecLap mleb (geom_tmp, ba_tmp, dm_tmp, info, {&ebfact});
68 
69  mleb.setMaxOrder(2);
70  mleb.setDomainBC(bclo, bchi);
71  if (lev > 0) {
72  mleb.setCoarseFineBC(nullptr, ref_ratio[lev-1], LinOpBCType::Neumann);
73  }
74  mleb.setLevelBC(0, nullptr);
75 
76  //
77  // This sets A = 0, B = 1 so that
78  // the operator A alpha - b del dot beta grad to b
79  // becomes - del dot beta grad
80  //
81  mleb.setScalars(zero, one);
82 
83  Array<MultiFab,AMREX_SPACEDIM> bcoef;
84  for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
85  bcoef[idim].define(convert(ba_tmp[0],IntVect::TheDimensionVector(idim)),
86  dm_tmp[0], 1, 0, MFInfo(), ebfact);
87  bcoef[idim].setVal(-1.0);
88  }
89  mleb.setBCoeffs(0, amrex::GetArrOfConstPtrs(bcoef));
90 
91  MLMG mlmg(mleb);
92 
93  int max_iter = 100;
94  mlmg.setMaxIter(max_iter);
95  mlmg.setVerbose(mg_verbose);
96  mlmg.setBottomVerbose(0);
97 
98  mlmg.solve(GetVecOfPtrs(phi), GetVecOfConstPtrs(rhs), reltol, abstol);
99 
100  mlmg.getFluxes(GetVecOfArrOfPtrs(fluxes));
101 
102  // Add the flux values (gradient phi) to the face-centered cell with zero area fraction
103  // (mlmg.getFluxes does not fill gradient phi at faces with zero area fraction)
104  FillZeroAreaFaceFluxes(phi[0], fluxes[0], geom, ebfact, ebfact_u, ebfact_v, ebfact_w);
105 
106  // ImposeBCsOnPhi(lev,phi[0], geom[lev].Domain());
107 
108  //
109  // This arises because we solve MINUS del dot beta grad phi = div (rho u)
110  //
111  fluxes[0][0].mult(-one);
112  fluxes[0][1].mult(-one);
113  fluxes[0][2].mult(-one);
114 }
one
constexpr amrex::Real one
Definition: ERF_Constants.H:7
dom_hi
const auto & dom_hi
Definition: ERF_SetupVertDiff.H:2
dom_lo
const auto & dom_lo
Definition: ERF_SetupVertDiff.H:1
FillZeroAreaFaceFluxes
void FillZeroAreaFaceFluxes(MultiFab &phi, Array< MultiFab, AMREX_SPACEDIM > &fluxes, const Geometry &geom, EBFArrayBoxFactory const &ebfact, eb_aux_ const &ebfact_u, eb_aux_ const &ebfact_v, eb_aux_ const &ebfact_w)
Definition: ERF_FillZeroAreaFaceFluxes.cpp:10
get_lo_projection_bc
Array< LinOpBCType, AMREX_SPACEDIM > get_lo_projection_bc(Geometry const &lev_geom, Array< std::string, 2 *AMREX_SPACEDIM > l_domain_bc_type)
Definition: ERF_SolverUtils.H:13
get_hi_projection_bc
Array< LinOpBCType, AMREX_SPACEDIM > get_hi_projection_bc(Geometry const &lev_geom, Array< std::string, 2 *AMREX_SPACEDIM > l_domain_bc_type)
Definition: ERF_SolverUtils.H:34
Here is the call graph for this function: