ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Functions
ERF_MakeGradP.cpp File Reference
#include <AMReX_MultiFab.H>
#include <AMReX_ArrayLim.H>
#include "AMReX_BCRec.H"
#include "ERF.H"
#include "ERF_SrcHeaders.H"
#include "ERF_DataStruct.H"
#include "ERF_Utils.H"
#include "ERF_EB.H"
#include <ERF_EBSlopes.H>
Include dependency graph for ERF_MakeGradP.cpp:

Functions

void make_gradp_pert (int level, const SolverChoice &solverChoice, const Geometry &geom, Vector< MultiFab > &S_data, const MultiFab &p0, const MultiFab &z_phys_nd, const MultiFab &z_phys_cc, Vector< std::unique_ptr< MultiFab >> &mapfac, const eb_ &ebfact, Vector< MultiFab > &gradp)
 
void compute_gradp (const MultiFab &p, const Geometry &geom, const MultiFab &z_phys_nd, const MultiFab &z_phys_cc, Vector< std::unique_ptr< MultiFab >> &mapfac, const eb_ &ebfact, Vector< MultiFab > &gradp, const SolverChoice &solverChoice)
 Compute the full pressure gradient. More...
 
void compute_gradp_xy (const MultiFab &p, const Geometry &geom, const MultiFab &z_phys_cc, Vector< std::unique_ptr< MultiFab >> &mapfac, const eb_ &ebfact, Vector< MultiFab > &gradp, const SolverChoice &solverChoice)
 Compute the horizontal components of the pressure gradient. More...
 
void compute_gradp_z (const MultiFab &p, const Geometry &geom, const MultiFab &z_phys_nd, const eb_ &ebfact, Vector< MultiFab > &gradp, const SolverChoice &solverChoice)
 Compute the vertical component of the pressure gradient. More...
 
void compute_gradp_interpz (const MultiFab &p, const Geometry &geom, const MultiFab &z_phys_nd, const MultiFab &z_phys_cc, Vector< std::unique_ptr< MultiFab >> &mapfac, Vector< MultiFab > &gradp, const SolverChoice &solverChoice)
 Compute the pressure gradient using vertical interpolation. More...
 

Function Documentation

◆ compute_gradp()

void compute_gradp ( const MultiFab &  p,
const Geometry &  geom,
const MultiFab &  z_phys_nd,
const MultiFab &  z_phys_cc,
Vector< std::unique_ptr< MultiFab >> &  mapfac,
const eb_ &  ebfact,
Vector< MultiFab > &  gradp,
const SolverChoice &  solverChoice 
)

Compute the full pressure gradient.

Parameters
[in]pPressure field.
[in]geomGeometry container.
[in]z_phys_ndPhysical height on nodes.
[in]z_phys_ccPhysical height on cell centers.
[in]mapfacMap factors.
[in]ebfactEB factory.
[out]gradpPressure gradient components.
[in]solverChoiceSolver options.
132 {
133  compute_gradp_xy(p,geom,z_phys_cc,mapfac,ebfact,gradp,solverChoice);
134  compute_gradp_z(p,geom,z_phys_nd,ebfact,gradp,solverChoice);
135 }
compute_gradp_z
void compute_gradp_z(const MultiFab &p, const Geometry &geom, const MultiFab &z_phys_nd, const eb_ &ebfact, Vector< MultiFab > &gradp, const SolverChoice &solverChoice)
Compute the vertical component of the pressure gradient.
Definition: ERF_MakeGradP.cpp:406
compute_gradp_xy
void compute_gradp_xy(const MultiFab &p, const Geometry &geom, const MultiFab &z_phys_cc, Vector< std::unique_ptr< MultiFab >> &mapfac, const eb_ &ebfact, Vector< MultiFab > &gradp, const SolverChoice &solverChoice)
Compute the horizontal components of the pressure gradient.
Definition: ERF_MakeGradP.cpp:148
WSM6Ind::p
@ p
Definition: ERF_WSM6.H:280

Referenced by ERF::compute_max_pressure_gradient_diagnostic(), and ERF::FillPlot3DVars().

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

◆ compute_gradp_interpz()

void compute_gradp_interpz ( const MultiFab &  p,
const Geometry &  geom,
const MultiFab &  z_phys_nd,
const MultiFab &  z_phys_cc,
Vector< std::unique_ptr< MultiFab >> &  mapfac,
Vector< MultiFab > &  gradp,
const SolverChoice &  solverChoice 
)

Compute the pressure gradient using vertical interpolation.

Parameters
[in]pPressure field.
[in]geomGeometry container.
[in]z_phys_ndPhysical height on nodes.
[in]z_phys_ccPhysical height on cell centers.
[in]mapfacMap factors.
[out]gradpPressure gradient components.
[in]solverChoiceSolver options.
552 {
553  const bool l_use_terrain_fitted_coords = (solverChoice.mesh_type != MeshType::ConstantDz);
554 
555  const Box domain = geom.Domain();
556  const int domain_klo = domain.smallEnd(2);
557  const int domain_khi = domain.bigEnd(2);
558 
559  const GpuArray<Real, AMREX_SPACEDIM> dxInv = geom.InvCellSizeArray();
560 
561  // *****************************************************************************
562  // Take gradient of relevant quantity (p0, pres, or pert_pres = pres - p0)
563  // *****************************************************************************
564  for ( MFIter mfi(p); mfi.isValid(); ++mfi)
565  {
566  Box tbx = mfi.nodaltilebox(0);
567  Box tby = mfi.nodaltilebox(1);
568  Box tbz = mfi.nodaltilebox(2);
569 
570  // We don't compute gpz on the bottom or top domain boundary
571  if (tbz.smallEnd(2) == domain_klo) {
572  tbz.growLo(2,-1);
573  }
574  if (tbz.bigEnd(2) == domain_khi+1) {
575  tbz.growHi(2,-1);
576  }
577 
578  // Terrain metrics
579  const Array4<const Real>& z_nd_arr = z_phys_nd.const_array(mfi);
580  const Array4<const Real>& z_cc_arr = z_phys_cc.const_array(mfi);
581 
582  const Array4<const Real>& p_arr = p.const_array(mfi);
583 
584  const Array4< Real>& gpx_arr = gradp[GpVars::gpx].array(mfi);
585  const Array4< Real>& gpy_arr = gradp[GpVars::gpy].array(mfi);
586  const Array4< Real>& gpz_arr = gradp[GpVars::gpz].array(mfi);
587 
588  const Array4<const Real>& mf_ux_arr = mapfac[MapFacType::u_x]->const_array(mfi);
589  const Array4<const Real>& mf_vy_arr = mapfac[MapFacType::v_y]->const_array(mfi);
590 
591  ParallelFor(tbx, tby, tbz,
592  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
593  {
594  if (l_use_terrain_fitted_coords) {
595  Real p_lo = p_arr(i-1,j,k);
596  Real p_hi = p_arr(i,j,k);
597  Real dz_int = myhalf * (z_cc_arr(i,j,k) - z_cc_arr(i-1,j,k));
598  if (dz_int > 0) {
599  // Klemp 2011, Eqn. 16: s = 1/2
600  if (k==domain_klo) {
601  p_hi = quad_interp_1d(z_cc_arr(i,j,k) - dz_int,
602  z_cc_arr(i,j,k ), p_arr(i,j,k ),
603  z_cc_arr(i,j,k+1), p_arr(i,j,k+1),
604  z_cc_arr(i,j,k+2), p_arr(i,j,k+2));
605  } else {
606  p_hi -= dz_int * ( ( p_arr(i ,j,k ) - p_arr(i ,j,k-1))
607  / (z_cc_arr(i ,j,k ) - z_cc_arr(i ,j,k-1)) );
608  }
609  if (k==domain_khi) {
610  p_lo = quad_interp_1d(z_cc_arr(i-1,j,k) + dz_int,
611  z_cc_arr(i-1,j,k-2), p_arr(i-1,j,k-2),
612  z_cc_arr(i-1,j,k-1), p_arr(i-1,j,k-1),
613  z_cc_arr(i-1,j,k ), p_arr(i-1,j,k ));
614  } else {
615  p_lo += dz_int * ( ( p_arr(i-1,j,k+1) - p_arr(i-1,j,k ))
616  / (z_cc_arr(i-1,j,k+1) - z_cc_arr(i-1,j,k )) );
617  }
618  } else if (dz_int < 0) {
619  // Klemp 2011, Eqn. 16: s = -1/2
620  if (k==domain_khi) {
621  p_hi = quad_interp_1d(z_cc_arr(i,j,k) - dz_int,
622  z_cc_arr(i,j,k-2), p_arr(i,j,k-2),
623  z_cc_arr(i,j,k-1), p_arr(i,j,k-1),
624  z_cc_arr(i,j,k ), p_arr(i,j,k ));
625  } else {
626  p_hi -= dz_int * ( ( p_arr(i ,j,k+1) - p_arr(i ,j,k ))
627  / (z_cc_arr(i ,j,k+1) - z_cc_arr(i ,j,k )) );
628  }
629  if (k==domain_klo) {
630  p_lo = quad_interp_1d(z_cc_arr(i-1,j,k) + dz_int,
631  z_cc_arr(i-1,j,k ), p_arr(i-1,j,k ),
632  z_cc_arr(i-1,j,k+1), p_arr(i-1,j,k+1),
633  z_cc_arr(i-1,j,k+2), p_arr(i-1,j,k+2));
634  } else {
635  p_lo += dz_int * ( ( p_arr(i-1,j,k ) - p_arr(i-1,j,k-1))
636  / (z_cc_arr(i-1,j,k ) - z_cc_arr(i-1,j,k-1)) );
637  }
638  }
639  gpx_arr(i,j,k) = dxInv[0] * (p_hi - p_lo);
640  } else {
641  gpx_arr(i,j,k) = dxInv[0] * (p_arr(i,j,k) - p_arr(i-1,j,k));
642  }
643 
644  // NOTE that the gradp array now carries the map factor!
645  gpx_arr(i,j,k) *= mf_ux_arr(i,j,0);
646  },
647  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
648  {
649  if (l_use_terrain_fitted_coords) {
650  Real p_lo = p_arr(i,j-1,k);
651  Real p_hi = p_arr(i,j,k);
652  Real dz_int = myhalf * (z_cc_arr(i,j,k) - z_cc_arr(i,j-1,k));
653  if (dz_int > 0) {
654  // Klemp 2011, Eqn. 16: s = 1/2
655  if (k==domain_klo) {
656  p_hi = quad_interp_1d(z_cc_arr(i,j,k) - dz_int,
657  z_cc_arr(i,j,k ), p_arr(i,j,k ),
658  z_cc_arr(i,j,k+1), p_arr(i,j,k+1),
659  z_cc_arr(i,j,k+2), p_arr(i,j,k+2));
660  } else {
661  p_hi -= dz_int * ( ( p_arr(i,j ,k ) - p_arr(i,j ,k-1))
662  / (z_cc_arr(i,j ,k ) - z_cc_arr(i,j ,k-1)) );
663  }
664  if (k==domain_khi) {
665  p_lo = quad_interp_1d(z_cc_arr(i,j-1,k) + dz_int,
666  z_cc_arr(i,j-1,k-2), p_arr(i,j-1,k-2),
667  z_cc_arr(i,j-1,k-1), p_arr(i,j-1,k-1),
668  z_cc_arr(i,j-1,k ), p_arr(i,j-1,k ));
669  } else {
670  p_lo += dz_int * ( ( p_arr(i,j-1,k+1) - p_arr(i,j-1,k ))
671  / (z_cc_arr(i,j-1,k+1) - z_cc_arr(i,j-1,k )) );
672  }
673  } else if (dz_int < 0) {
674  // Klemp 2011, Eqn. 16: s = -1/2
675  if (k==domain_khi) {
676  p_hi = quad_interp_1d(z_cc_arr(i,j,k) - dz_int,
677  z_cc_arr(i,j,k-2), p_arr(i,j,k-2),
678  z_cc_arr(i,j,k-1), p_arr(i,j,k-1),
679  z_cc_arr(i,j,k ), p_arr(i,j,k ));
680  } else {
681  p_hi -= dz_int * ( ( p_arr(i,j ,k+1) - p_arr(i,j ,k ))
682  / (z_cc_arr(i,j ,k+1) - z_cc_arr(i,j ,k )) );
683  }
684  if (k==domain_klo) {
685  p_lo = quad_interp_1d(z_cc_arr(i,j-1,k) + dz_int,
686  z_cc_arr(i,j-1,k ), p_arr(i,j-1,k ),
687  z_cc_arr(i,j-1,k+1), p_arr(i,j-1,k+1),
688  z_cc_arr(i,j-1,k+2), p_arr(i,j-1,k+2));
689  } else {
690  p_lo += dz_int * ( ( p_arr(i,j-1,k ) - p_arr(i,j-1,k-1))
691  / (z_cc_arr(i,j-1,k ) - z_cc_arr(i,j-1,k-1)) );
692  }
693  }
694  gpy_arr(i,j,k) = dxInv[1] * (p_hi - p_lo);
695  } else {
696  gpy_arr(i,j,k) = dxInv[1] * (p_arr(i,j,k) - p_arr(i,j-1,k));
697  }
698 
699  // NOTE that the gradp array now carries the map factor!
700  gpy_arr(i,j,k) *= mf_vy_arr(i,j,0);
701  },
702  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
703  {
704  // Note: identical to gradp_type == 0
705  Real met_h_zeta = (l_use_terrain_fitted_coords) ? Compute_h_zeta_AtKface(i, j, k, dxInv, z_nd_arr) : 1;
706  gpz_arr(i,j,k) = dxInv[2] * ( p_arr(i,j,k)-p_arr(i,j,k-1) ) / met_h_zeta;
707  });
708  } // mfi
709 }
v_y
@ v_y
Definition: ERF_DataStruct.H:30
u_x
@ u_x
Definition: ERF_DataStruct.H:29
dxInv
amrex::GpuArray< Real, AMREX_SPACEDIM > dxInv
Definition: ERF_InitCustomPertVels_ParticleTests.H:17
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);})
myhalf
constexpr amrex::Real myhalf
Definition: ERF_NumericalConstants.H:34
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
Compute_h_zeta_AtKface
AMREX_GPU_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:409
quad_interp_1d
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real quad_interp_1d(amrex::Real z, amrex::Real z0, amrex::Real p0, amrex::Real z1, amrex::Real p1, amrex::Real z2, amrex::Real p2)
Definition: ERF_Utils.H:753
GpVars::gpz
@ gpz
Definition: ERF_IndexDefines.H:226
GpVars::gpy
@ gpy
Definition: ERF_IndexDefines.H:225
GpVars::gpx
@ gpx
Definition: ERF_IndexDefines.H:224
SolverChoice::mesh_type
static MeshType mesh_type
Vertical mesh representation.
Definition: ERF_DataStruct.H:1958

Referenced by make_gradp_pert().

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

◆ compute_gradp_xy()

void compute_gradp_xy ( const MultiFab &  p,
const Geometry &  geom,
const MultiFab &  z_phys_cc,
Vector< std::unique_ptr< MultiFab >> &  mapfac,
const eb_ &  ebfact,
Vector< MultiFab > &  gradp,
const SolverChoice &  solverChoice 
)

Compute the horizontal components of the pressure gradient.

Parameters
[in]pPressure field.
[in]geomGeometry container.
[in]z_phys_ccPhysical height on cell centers.
[in]mapfacMap factors.
[in]ebfactEB factory.
[out]gradpPressure gradient components.
[in]solverChoiceSolver options.
155 {
156  const bool l_use_terrain_fitted_coords = (solverChoice.mesh_type != MeshType::ConstantDz);
157 
158  const Box domain = geom.Domain();
159  const int domain_klo = domain.smallEnd(2);
160  const int domain_khi = domain.bigEnd(2);
161 
162  const GpuArray<Real, AMREX_SPACEDIM> dxInv = geom.InvCellSizeArray();
163 
164  // *****************************************************************************
165  // Take gradient of relevant quantity (p0, pres, or pert_pres = pres - p0)
166  // *****************************************************************************
167  for ( MFIter mfi(p); mfi.isValid(); ++mfi)
168  {
169  Box tbx = mfi.nodaltilebox(0);
170  Box tby = mfi.nodaltilebox(1);
171 
172  // Terrain metrics
173  const Array4<const Real>& z_cc_arr = z_phys_cc.const_array(mfi);
174 
175  const Array4<const Real>& p_arr = p.const_array(mfi);
176 
177  const Array4< Real>& gpx_arr = gradp[GpVars::gpx].array(mfi);
178  const Array4< Real>& gpy_arr = gradp[GpVars::gpy].array(mfi);
179 
180  const Array4<const Real>& mf_ux_arr = mapfac[MapFacType::u_x]->const_array(mfi);
181  const Array4<const Real>& mf_vy_arr = mapfac[MapFacType::v_y]->const_array(mfi);
182 
183  if (solverChoice.terrain_type != TerrainType::EB) {
184 
185  ParallelFor(tbx, tby,
186  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
187  {
188  //Note : mx/my == 1, so no map factor needed here
189  Real gpx = dxInv[0] * (p_arr(i,j,k) - p_arr(i-1,j,k));
190 
191  if (l_use_terrain_fitted_coords) {
192  Real met_h_xi = (z_cc_arr(i,j,k) - z_cc_arr(i-1,j,k)) * dxInv[0];
193 
194  Real dz_phys_hi, dz_phys_lo;
195  Real gpz_lo, gpz_hi;
196  if (k==domain_klo) {
197  dz_phys_hi = z_cc_arr(i ,j,k+1) - z_cc_arr(i ,j,k );
198  dz_phys_lo = z_cc_arr(i-1,j,k+1) - z_cc_arr(i-1,j,k );
199  gpz_hi = (p_arr(i ,j,k+1) - p_arr(i ,j,k )) / dz_phys_hi;
200  gpz_lo = (p_arr(i-1,j,k+1) - p_arr(i-1,j,k )) / dz_phys_lo;
201  } else if (k==domain_khi) {
202  dz_phys_hi = z_cc_arr(i ,j,k ) - z_cc_arr(i ,j,k-1);
203  dz_phys_lo = z_cc_arr(i-1,j,k ) - z_cc_arr(i-1,j,k-1);
204  gpz_hi = (p_arr(i ,j,k ) - p_arr(i ,j,k-1)) / dz_phys_hi;
205  gpz_lo = (p_arr(i-1,j,k ) - p_arr(i-1,j,k-1)) / dz_phys_lo;
206  } else {
207  dz_phys_hi = z_cc_arr(i ,j,k+1) - z_cc_arr(i ,j,k-1);
208  dz_phys_lo = z_cc_arr(i-1,j,k+1) - z_cc_arr(i-1,j,k-1);
209  gpz_hi = (p_arr(i ,j,k+1) - p_arr(i ,j,k-1)) / dz_phys_hi;
210  gpz_lo = (p_arr(i-1,j,k+1) - p_arr(i-1,j,k-1)) / dz_phys_lo;
211  }
212  Real gpx_metric = met_h_xi * myhalf * (gpz_hi + gpz_lo);
213  gpx -= gpx_metric;
214  }
215  gpx_arr(i,j,k) = gpx;
216 
217  // NOTE that the gradp array now carries the map factor!
218  gpx_arr(i,j,k) *= mf_ux_arr(i,j,0);
219  },
220  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
221  {
222  //Note : mx/my == 1, so no map factor needed here
223  Real gpy = dxInv[1] * (p_arr(i,j,k) - p_arr(i,j-1,k));
224 
225  if (l_use_terrain_fitted_coords) {
226  Real met_h_eta = (z_cc_arr(i,j,k) - z_cc_arr(i,j-1,k)) * dxInv[1];
227 
228  Real dz_phys_hi, dz_phys_lo;
229  Real gpz_lo, gpz_hi;
230  if (k==domain_klo) {
231  dz_phys_hi = z_cc_arr(i,j ,k+1) - z_cc_arr(i,j ,k );
232  dz_phys_lo = z_cc_arr(i,j-1,k+1) - z_cc_arr(i,j-1,k );
233  gpz_hi = (p_arr(i,j ,k+1) - p_arr(i,j ,k )) / dz_phys_hi;
234  gpz_lo = (p_arr(i,j-1,k+1) - p_arr(i,j-1,k )) / dz_phys_lo;
235  } else if (k==domain_khi) {
236  dz_phys_hi = z_cc_arr(i,j ,k ) - z_cc_arr(i,j ,k-1);
237  dz_phys_lo = z_cc_arr(i,j-1,k ) - z_cc_arr(i,j-1,k-1);
238  gpz_hi = (p_arr(i,j ,k ) - p_arr(i,j ,k-1)) / dz_phys_hi;
239  gpz_lo = (p_arr(i,j-1,k ) - p_arr(i,j-1,k-1)) / dz_phys_lo;
240  } else {
241  dz_phys_hi = z_cc_arr(i,j ,k+1) - z_cc_arr(i,j ,k-1);
242  dz_phys_lo = z_cc_arr(i,j-1,k+1) - z_cc_arr(i,j-1,k-1);
243  gpz_hi = (p_arr(i,j ,k+1) - p_arr(i,j ,k-1)) / dz_phys_hi;
244  gpz_lo = (p_arr(i,j-1,k+1) - p_arr(i,j-1,k-1)) / dz_phys_lo;
245  }
246  Real gpy_metric = met_h_eta * myhalf * (gpz_hi + gpz_lo);
247  gpy -= gpy_metric;
248  }
249  gpy_arr(i,j,k) = gpy;
250 
251  // NOTE that the gradp array now carries the map factor!
252  gpy_arr(i,j,k) *= mf_vy_arr(i,j,0);
253  });
254 
255  // solverChoice.terrain_type == TerrainType::EB
256  } else {
257 
258  // Pressure gradients are fitted at the centroids of cut cells, if EB and Compressible.
259  // Least-Squares Fitting: Compute slope using 3x3x3 stencil
260 
261  const bool l_fitting = false;
262 
263  const Real* dx_arr = geom.CellSize();
264  const Real dx = dx_arr[0];
265  const Real dy = dx_arr[1];
266  const Real dz = dx_arr[2];
267 
268  // EB factory
269  Array4<const EBCellFlag> cellflg = (ebfact.get_const_factory())->getMultiEBCellFlagFab()[mfi].const_array();
270 
271  // EB u-factory
272  auto const* u_factory = ebfact.get_u_const_factory();
273  Array4<const EBCellFlag> u_cellflg = u_factory->getMultiEBCellFlagFab()[mfi].const_array();
274  Array4<const Real > u_volfrac = u_factory->getVolFrac().const_array(mfi);
275  bool u_is_cut = (u_factory->getMultiEBCellFlagFab()[mfi].getType() == FabType::singlevalued);
276  Array4<const Real > u_volcent = u_is_cut ? u_factory->getCentroid().const_array(mfi) : Array4<const Real>{};
277 
278  // EB v-factory
279  auto const* v_factory = ebfact.get_v_const_factory();
280  Array4<const EBCellFlag> v_cellflg = v_factory->getMultiEBCellFlagFab()[mfi].const_array();
281  Array4<const Real > v_volfrac = v_factory->getVolFrac().const_array(mfi);
282  bool v_is_cut = (v_factory->getMultiEBCellFlagFab()[mfi].getType() == FabType::singlevalued);
283  Array4<const Real > v_volcent = v_is_cut ? v_factory->getCentroid().const_array(mfi) : Array4<const Real>{};
284 
285  if (l_fitting) {
286 
287  ParallelFor(tbx, tby,
288  [=] AMREX_GPU_DEVICE(int i, int j, int k)
289  {
290  if (u_volfrac(i,j,k) > zero) {
291 
292  if (u_cellflg(i,j,k).isSingleValued()) {
293 
294  GpuArray<Real,AMREX_SPACEDIM> slopes;
295  slopes = erf_calc_slopes_eb_staggered(Vars::xvel, Vars::cons, dx, dy, dz, i, j, k, p_arr, u_volcent, u_cellflg);
296 
297  gpx_arr(i,j,k) = slopes[0];
298 
299  } else {
300  gpx_arr(i,j,k) = dxInv[0] * (p_arr(i,j,k) - p_arr(i-1,j,k));
301  }
302 
303  } else {
304  gpx_arr(i,j,k) = zero;
305  }
306  },
307  [=] AMREX_GPU_DEVICE(int i, int j, int k)
308  {
309  if (v_volfrac(i,j,k) > zero) {
310 
311  if (v_cellflg(i,j,k).isSingleValued()) {
312 
313  GpuArray<Real,AMREX_SPACEDIM> slopes;
314  slopes = erf_calc_slopes_eb_staggered(Vars::yvel, Vars::cons, dx, dy, dz, i, j, k, p_arr, v_volcent, v_cellflg);
315 
316  gpy_arr(i,j,k) = slopes[1];
317 
318  } else {
319  gpy_arr(i,j,k) = dxInv[1] * (p_arr(i,j,k) - p_arr(i,j-1,k));
320  }
321  } else {
322  gpy_arr(i,j,k) = zero;
323  }
324  });
325 
326  } else {
327 
328  // Simple calculation: assuming pressures at cell centers
329 
330  ParallelFor(tbx, tby,
331  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
332  {
333  if (u_volfrac(i,j,k) > zero) {
334 
335  if (cellflg(i,j,k).isCovered()) {
336  // Check if all donor cells on the low side are not covered
337  if (!cellflg(i-1,j,k).isCovered() && !cellflg(i-2,j,k).isCovered() && !cellflg(i-3,j,k).isCovered()) {
338  // Use three-point stencil (second-order accurate)
339  gpx_arr(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));
340  } else {
341  // Fall back to one-sided two-point stencil (first-order)
342  gpx_arr(i,j,k) = dxInv[0] * (p_arr(i-1,j,k) - p_arr(i-2,j,k));
343  }
344  } else if (cellflg(i-1,j,k).isCovered()) {
345  // Check if all donor cells on the high side are not covered
346  if (!cellflg(i,j,k).isCovered() && !cellflg(i+1,j,k).isCovered() && !cellflg(i+2,j,k).isCovered()) {
347  // Use three-point stencil (second-order accurate)
348  gpx_arr(i,j,k) = dxInv[0] * (three*p_arr(i+1,j,k) - p_arr(i+2,j,k) - two*p_arr(i,j,k));
349  } else {
350  // Fall back to one-sided two-point stencil (first-order)
351  gpx_arr(i,j,k) = dxInv[0] * (p_arr(i+1,j,k) - p_arr(i,j,k));
352  }
353  } else {
354  gpx_arr(i,j,k) = dxInv[0] * (p_arr(i,j,k) - p_arr(i-1,j,k));
355  }
356  } else {
357  gpx_arr(i,j,k) = zero;
358  }
359  },
360  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
361  {
362  if (v_volfrac(i,j,k) > zero) {
363  if (cellflg(i,j,k).isCovered()) {
364  // Check if all donor cells on the low side are not covered
365  if (!cellflg(i,j-1,k).isCovered() && !cellflg(i,j-2,k).isCovered() && !cellflg(i,j-3,k).isCovered()) {
366  // Use three-point stencil (second-order accurate)
367  gpy_arr(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));
368  } else {
369  // Fall back to one-sided two-point stencil (first-order)
370  gpy_arr(i,j,k) = dxInv[1] * (p_arr(i,j-1,k) - p_arr(i,j-2,k));
371  }
372  } else if (cellflg(i,j-1,k).isCovered()) {
373  // Check if all donor cells on the high side are not covered
374  if (!cellflg(i,j,k).isCovered() && !cellflg(i,j+1,k).isCovered() && !cellflg(i,j+2,k).isCovered()) {
375  // Use three-point stencil (second-order accurate)
376  gpy_arr(i,j,k) = dxInv[1] * (three*p_arr(i,j+1,k) - p_arr(i,j+2,k) - two*p_arr(i,j,k));
377  } else {
378  // Fall back to one-sided two-point stencil (first-order)
379  gpy_arr(i,j,k) = dxInv[1] * (p_arr(i,j+1,k) - p_arr(i,j,k));
380  }
381  } else {
382  gpy_arr(i,j,k) = dxInv[1] * (p_arr(i,j,k) - p_arr(i,j-1,k));
383  }
384  } else {
385  gpy_arr(i,j,k) = zero;
386  }
387  });
388 
389  } // l_fitting
390 
391  } // TerrainType::EB
392 
393  } // mfi
394 }
erf_calc_slopes_eb_staggered
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > erf_calc_slopes_eb_staggered(int igrid_query, int igrid_data, amrex::Real dx, amrex::Real dy, amrex::Real dz, int i, int j, int k, amrex::Array4< amrex::Real const > const &state, amrex::Array4< amrex::Real const > const &ccent, amrex::Array4< amrex::EBCellFlag const > const &flag)
Compute centered least-squares slopes from cell data to a staggered query grid.
Definition: ERF_EBSlopes.H:318
dy
const Real dy
Definition: ERF_InitCustomPert_ABL.H:45
dx
const Real dx
Definition: ERF_InitCustomPert_ABL.H:44
three
constexpr amrex::Real three
Definition: ERF_NumericalConstants.H:32
two
constexpr amrex::Real two
Definition: ERF_NumericalConstants.H:31
zero
constexpr amrex::Real zero
Definition: ERF_NumericalConstants.H:29
eb_::get_const_factory
const std::unique_ptr< amrex::EBFArrayBoxFactory > & get_const_factory() const noexcept
Return the cell-centered EB factory.
Definition: ERF_EB.H:102
eb_::get_v_const_factory
eb_aux_ const * get_v_const_factory() const noexcept
Return the ERF auxiliary y-face EB factory.
Definition: ERF_EB.H:131
eb_::get_u_const_factory
eb_aux_ const * get_u_const_factory() const noexcept
Return the ERF auxiliary x-face EB factory.
Definition: ERF_EB.H:129
eb_aux_::getMultiEBCellFlagFab
const amrex::FabArray< amrex::EBCellFlagFab > & getMultiEBCellFlagFab() const
Return the reconstructed EB cell flags.
Definition: ERF_EBAux.cpp:1146
Vars::xvel
@ xvel
Definition: ERF_IndexDefines.H:215
Vars::cons
@ cons
Definition: ERF_IndexDefines.H:214
Vars::yvel
@ yvel
Definition: ERF_IndexDefines.H:216
WDM6SedCellScratch::dz
@ dz
Definition: ERF_AdvanceWDM6.cpp:272
SolverChoice::terrain_type
static TerrainType terrain_type
Terrain or immersed-boundary representation.
Definition: ERF_DataStruct.H:1949

Referenced by compute_gradp(), and make_gradp_pert().

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

◆ compute_gradp_z()

void compute_gradp_z ( const MultiFab &  p,
const Geometry &  geom,
const MultiFab &  z_phys_nd,
const eb_ &  ebfact,
Vector< MultiFab > &  gradp,
const SolverChoice &  solverChoice 
)

Compute the vertical component of the pressure gradient.

Parameters
[in]pPressure field.
[in]geomGeometry container.
[in]z_phys_ndPhysical height on nodes.
[in]ebfactEB factory.
[out]gradpPressure gradient components.
[in]solverChoiceSolver options.
412 {
413  const bool l_use_terrain_fitted_coords = (solverChoice.mesh_type != MeshType::ConstantDz);
414 
415  const Box domain = geom.Domain();
416  const int domain_klo = domain.smallEnd(2);
417  const int domain_khi = domain.bigEnd(2);
418 
419  const GpuArray<Real, AMREX_SPACEDIM> dxInv = geom.InvCellSizeArray();
420 
421  // *****************************************************************************
422  // Take gradient of relevant quantity (p0, pres, or pert_pres = pres - p0)
423  // *****************************************************************************
424  for ( MFIter mfi(p); mfi.isValid(); ++mfi)
425  {
426  Box tbz = mfi.nodaltilebox(2);
427 
428  // We don't compute gpz on the bottom or top domain boundary
429  if (tbz.smallEnd(2) == domain_klo) {
430  tbz.growLo(2,-1);
431  }
432  if (tbz.bigEnd(2) == domain_khi+1) {
433  tbz.growHi(2,-1);
434  }
435 
436  // Terrain metrics
437  const Array4<const Real>& z_nd_arr = z_phys_nd.const_array(mfi);
438 
439  const Array4<const Real>& p_arr = p.const_array(mfi);
440 
441  const Array4< Real>& gpz_arr = gradp[GpVars::gpz].array(mfi);
442 
443  if (solverChoice.terrain_type != TerrainType::EB) {
444 
445  ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
446  {
447  Real met_h_zeta = (l_use_terrain_fitted_coords) ? Compute_h_zeta_AtKface(i, j, k, dxInv, z_nd_arr) : 1;
448  gpz_arr(i,j,k) = dxInv[2] * ( p_arr(i,j,k)-p_arr(i,j,k-1) ) / met_h_zeta;
449  });
450 
451  } else {
452 
453  // Pressure gradients are fitted at the centroids of cut cells, if EB and Compressible.
454  // Least-Squares Fitting: Compute slope using 3x3x3 stencil
455 
456  const bool l_fitting = false;
457 
458  const Real* dx_arr = geom.CellSize();
459  const Real dx = dx_arr[0];
460  const Real dy = dx_arr[1];
461  const Real dz = dx_arr[2];
462 
463  // EB factory
464  Array4<const EBCellFlag> cellflg = (ebfact.get_const_factory())->getMultiEBCellFlagFab()[mfi].const_array();
465 
466  // EB w-factory
467  auto const* w_factory = ebfact.get_w_const_factory();
468  Array4<const EBCellFlag> w_cellflg = w_factory->getMultiEBCellFlagFab()[mfi].const_array();
469  Array4<const Real > w_volfrac = w_factory->getVolFrac().const_array(mfi);
470  bool w_is_cut = (w_factory->getMultiEBCellFlagFab()[mfi].getType() == FabType::singlevalued);
471  Array4<const Real > w_volcent = w_is_cut ? w_factory->getCentroid().const_array(mfi) : Array4<Real>{};
472 
473  if (l_fitting) {
474 
475  ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k)
476  {
477  if (w_volfrac(i,j,k) > zero) {
478 
479  if (w_cellflg(i,j,k).isSingleValued()) {
480 
481  GpuArray<Real,AMREX_SPACEDIM> slopes;
482  slopes = erf_calc_slopes_eb_staggered(Vars::zvel, Vars::cons, dx, dy, dz, i, j, k, p_arr, w_volcent, w_cellflg);
483 
484  gpz_arr(i,j,k) = slopes[2];
485 
486  } else {
487  gpz_arr(i,j,k) = dxInv[2] * (p_arr(i,j,k) - p_arr(i,j,k-1));
488  }
489  } else {
490  gpz_arr(i,j,k) = zero;
491  }
492  });
493 
494  } else {
495 
496  // Simple calculation: assuming pressures at cell centers
497 
498  ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
499  {
500  if (w_volfrac(i,j,k) > zero) {
501  if (cellflg(i,j,k).isCovered()) {
502  // Check if all donor cells on the low side are not covered
503  if (!cellflg(i,j,k-1).isCovered() && !cellflg(i,j,k-2).isCovered() && !cellflg(i,j,k-3).isCovered()) {
504  // Use three-point stencil (second-order accurate)
505  gpz_arr(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) );
506  } else {
507  // Fall back to one-sided two-point stencil (first-order)
508  gpz_arr(i,j,k) = dxInv[2] * ( p_arr(i,j,k-1) - p_arr(i,j,k-2) );
509  }
510  } else if (cellflg(i,j,k-1).isCovered()) {
511  // Check if all donor cells on the high side are not covered
512  if (!cellflg(i,j,k).isCovered() && !cellflg(i,j,k+1).isCovered() && !cellflg(i,j,k+2).isCovered()) {
513  // Use three-point stencil (second-order accurate)
514  gpz_arr(i,j,k) = dxInv[2] * ( three*p_arr(i,j,k+1) - p_arr(i,j,k+2) - two*p_arr(i,j,k) );
515  } else {
516  // Fall back to one-sided two-point stencil (first-order)
517  gpz_arr(i,j,k) = dxInv[2] * ( p_arr(i,j,k+1) - p_arr(i,j,k) );
518  }
519  } else {
520  gpz_arr(i,j,k) = dxInv[2] * ( p_arr(i,j,k)-p_arr(i,j,k-1) );
521  }
522  } else {
523  gpz_arr(i,j,k) = zero;
524  }
525  });
526 
527  } // l_fitting
528 
529  } // TerrainType::EB
530 
531  } // mfi
532 }
eb_::get_w_const_factory
eb_aux_ const * get_w_const_factory() const noexcept
Return the ERF auxiliary z-face EB factory.
Definition: ERF_EB.H:133
Vars::zvel
@ zvel
Definition: ERF_IndexDefines.H:217

Referenced by compute_gradp(), and make_gradp_pert().

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

◆ make_gradp_pert()

void make_gradp_pert ( int  level,
const SolverChoice &  solverChoice,
const Geometry &  geom,
Vector< MultiFab > &  S_data,
const MultiFab &  p0,
const MultiFab &  z_phys_nd,
const MultiFab &  z_phys_cc,
Vector< std::unique_ptr< MultiFab >> &  mapfac,
const eb_ &  ebfact,
Vector< MultiFab > &  gradp 
)

Function for computing the pressure gradient

Parameters
[in]levellevel of resolution
[in]geomgeometry container at this level
[in]S_datacurrent solution
[in]p0base ststa pressure
[in]z_phys_ndz on nodes
[in]z_phys_ccz on cell centers
[in]d_bcrec_ptrBoundary Condition Record
[in]ebfactEB factory container at this level
[out]gradppressure gradient
38 {
39  const bool l_use_moisture = (solverChoice.moisture_type != MoistureType::None);
40  const bool l_eb_terrain = (solverChoice.terrain_type == TerrainType::EB);
41 
42  const bool l_use_pert_pres = (solverChoice.use_pert_pres_gradient);
43  //
44  // Note that we only recompute gradp if compressible;
45  // if anelastic then we have computed gradp in the projection
46  // and we can reuse it, no need to recompute it
47  //
48  if (solverChoice.anelastic[level] == 0)
49  {
50  if (solverChoice.gradp_type == 1) {
51  AMREX_ASSERT_WITH_MESSAGE(solverChoice.terrain_type != TerrainType::EB,
52  "gradp_type==1 not implemented for EB");
53  }
54 
55  const int ngrow = (l_eb_terrain) ? 3 : 1;
56  MultiFab p(S_data[Vars::cons].boxArray(), S_data[Vars::cons].DistributionMap(), 1, ngrow);
57 
58  // *****************************************************************************
59  // Compute pressure
60  // *****************************************************************************
61  for ( MFIter mfi(S_data[Vars::cons]); mfi.isValid(); ++mfi)
62  {
63  Box gbx = mfi.tilebox();
64  gbx.grow(IntVect(ngrow,ngrow,ngrow));
65 
66  if (gbx.smallEnd(2) < 0) gbx.setSmall(2,0);
67  const Array4<const Real>& cell_data = S_data[Vars::cons].array(mfi);
68  const Array4< Real>& pp_arr = p.array(mfi);
69  ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
70  {
71  Real qv_for_p = (l_use_moisture) ? cell_data(i,j,k,RhoQ1_comp)/cell_data(i,j,k,Rho_comp) : zero;
72  pp_arr(i,j,k) = getPgivenRTh(cell_data(i,j,k,RhoTheta_comp),qv_for_p);
73  });
74  }
75 
76  // If we want to use the full pressure in the lateral gradients, call compute_gradp_xy here
77  if (solverChoice.gradp_type == 0 && !l_use_pert_pres) {
78  compute_gradp_xy(p,geom,z_phys_cc,mapfac,ebfact,gradp,solverChoice);
79  }
80 
81  // *****************************************************************************
82  // Compute perturbational pressure
83  // *****************************************************************************
84  for ( MFIter mfi(S_data[Vars::cons]); mfi.isValid(); ++mfi)
85  {
86  Box gbx = mfi.tilebox();
87  gbx.grow(IntVect(ngrow,ngrow,ngrow));
88 
89  if (gbx.smallEnd(2) < 0) gbx.setSmall(2,0);
90  const Array4<const Real>& p0_arr = p0.const_array(mfi);
91  const Array4< Real>& pp_arr = p.array(mfi);
92  ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
93  {
94  pp_arr(i,j,k) -= p0_arr(i,j,k);
95  });
96  }
97 
98  // If we want to use the perturbational pressure in the lateral gradients, call compute_gradp_xy here
99  if (solverChoice.gradp_type == 0 && l_use_pert_pres) {
100  compute_gradp_xy(p,geom,z_phys_cc,mapfac,ebfact,gradp,solverChoice);
101  }
102 
103  if (solverChoice.gradp_type == 0) {
104  compute_gradp_z(p,geom,z_phys_nd,ebfact,gradp,solverChoice);
105  } else {
106  compute_gradp_interpz(p,geom,z_phys_nd,z_phys_cc,mapfac,gradp,solverChoice);
107  }
108 
109  } // not anelastic
110 }
getPgivenRTh
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getPgivenRTh(const amrex::Real rhotheta, const amrex::Real qv=amrex::Real(0))
Definition: ERF_EOS.H:81
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:39
RhoTheta_comp
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:40
RhoQ1_comp
#define RhoQ1_comp
Definition: ERF_IndexDefines.H:45
compute_gradp_interpz
void compute_gradp_interpz(const MultiFab &p, const Geometry &geom, const MultiFab &z_phys_nd, const MultiFab &z_phys_cc, Vector< std::unique_ptr< MultiFab >> &mapfac, Vector< MultiFab > &gradp, const SolverChoice &solverChoice)
Compute the pressure gradient using vertical interpolation.
Definition: ERF_MakeGradP.cpp:545
AMREX_ASSERT_WITH_MESSAGE
AMREX_ASSERT_WITH_MESSAGE(wbar_cutoff_min > wbar_cutoff_max, "ERROR: wbar_cutoff_min < wbar_cutoff_max")
module_model_constants::p0
real(c_double), parameter p0
Definition: ERF_module_model_constants.F90:40
SolverChoice::moisture_type
MoistureType moisture_type
Moisture or microphysics model.
Definition: ERF_DataStruct.H:2237
SolverChoice::use_pert_pres_gradient
bool use_pert_pres_gradient
Whether momentum equations use perturbational pressure gradients.
Definition: ERF_DataStruct.H:2022
SolverChoice::anelastic
amrex::Vector< int > anelastic
Per-level flag selecting anelastic dynamics.
Definition: ERF_DataStruct.H:1981
SolverChoice::gradp_type
int gradp_type
Terrain-fitted horizontal pressure-gradient formulation.
Definition: ERF_DataStruct.H:2020

Referenced by ERF::compute_max_buoyancy_gradp_diagnostic().

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