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:377
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:191

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

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.
509 {
510  const bool l_use_terrain_fitted_coords = (solverChoice.mesh_type != MeshType::ConstantDz);
511 
512  const Box domain = geom.Domain();
513  const int domain_klo = domain.smallEnd(2);
514  const int domain_khi = domain.bigEnd(2);
515 
516  const GpuArray<Real, AMREX_SPACEDIM> dxInv = geom.InvCellSizeArray();
517 
518  // *****************************************************************************
519  // Take gradient of relevant quantity (p0, pres, or pert_pres = pres - p0)
520  // *****************************************************************************
521  for ( MFIter mfi(p); mfi.isValid(); ++mfi)
522  {
523  Box tbx = mfi.nodaltilebox(0);
524  Box tby = mfi.nodaltilebox(1);
525  Box tbz = mfi.nodaltilebox(2);
526 
527  // We don't compute gpz on the bottom or top domain boundary
528  if (tbz.smallEnd(2) == domain_klo) {
529  tbz.growLo(2,-1);
530  }
531  if (tbz.bigEnd(2) == domain_khi+1) {
532  tbz.growHi(2,-1);
533  }
534 
535  // Terrain metrics
536  const Array4<const Real>& z_nd_arr = z_phys_nd.const_array(mfi);
537  const Array4<const Real>& z_cc_arr = z_phys_cc.const_array(mfi);
538 
539  const Array4<const Real>& p_arr = p.const_array(mfi);
540 
541  const Array4< Real>& gpx_arr = gradp[GpVars::gpx].array(mfi);
542  const Array4< Real>& gpy_arr = gradp[GpVars::gpy].array(mfi);
543  const Array4< Real>& gpz_arr = gradp[GpVars::gpz].array(mfi);
544 
545  const Array4<const Real>& mf_ux_arr = mapfac[MapFacType::u_x]->const_array(mfi);
546  const Array4<const Real>& mf_vy_arr = mapfac[MapFacType::v_y]->const_array(mfi);
547 
548  ParallelFor(tbx, tby, tbz,
549  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
550  {
551  if (l_use_terrain_fitted_coords) {
552  Real p_lo = p_arr(i-1,j,k);
553  Real p_hi = p_arr(i,j,k);
554  Real dz_int = myhalf * (z_cc_arr(i,j,k) - z_cc_arr(i-1,j,k));
555  if (dz_int > 0) {
556  // Klemp 2011, Eqn. 16: s = 1/2
557  if (k==domain_klo) {
558  p_hi = quad_interp_1d(z_cc_arr(i,j,k) - dz_int,
559  z_cc_arr(i,j,k ), p_arr(i,j,k ),
560  z_cc_arr(i,j,k+1), p_arr(i,j,k+1),
561  z_cc_arr(i,j,k+2), p_arr(i,j,k+2));
562  } else {
563  p_hi -= dz_int * ( ( p_arr(i ,j,k ) - p_arr(i ,j,k-1))
564  / (z_cc_arr(i ,j,k ) - z_cc_arr(i ,j,k-1)) );
565  }
566  if (k==domain_khi) {
567  p_lo = quad_interp_1d(z_cc_arr(i-1,j,k) + dz_int,
568  z_cc_arr(i-1,j,k-2), p_arr(i-1,j,k-2),
569  z_cc_arr(i-1,j,k-1), p_arr(i-1,j,k-1),
570  z_cc_arr(i-1,j,k ), p_arr(i-1,j,k ));
571  } else {
572  p_lo += dz_int * ( ( p_arr(i-1,j,k+1) - p_arr(i-1,j,k ))
573  / (z_cc_arr(i-1,j,k+1) - z_cc_arr(i-1,j,k )) );
574  }
575  } else if (dz_int < 0) {
576  // Klemp 2011, Eqn. 16: s = -1/2
577  if (k==domain_khi) {
578  p_hi = quad_interp_1d(z_cc_arr(i,j,k) - dz_int,
579  z_cc_arr(i,j,k-2), p_arr(i,j,k-2),
580  z_cc_arr(i,j,k-1), p_arr(i,j,k-1),
581  z_cc_arr(i,j,k ), p_arr(i,j,k ));
582  } else {
583  p_hi -= dz_int * ( ( p_arr(i ,j,k+1) - p_arr(i ,j,k ))
584  / (z_cc_arr(i ,j,k+1) - z_cc_arr(i ,j,k )) );
585  }
586  if (k==domain_klo) {
587  p_lo = quad_interp_1d(z_cc_arr(i-1,j,k) + dz_int,
588  z_cc_arr(i-1,j,k ), p_arr(i-1,j,k ),
589  z_cc_arr(i-1,j,k+1), p_arr(i-1,j,k+1),
590  z_cc_arr(i-1,j,k+2), p_arr(i-1,j,k+2));
591  } else {
592  p_lo += dz_int * ( ( p_arr(i-1,j,k ) - p_arr(i-1,j,k-1))
593  / (z_cc_arr(i-1,j,k ) - z_cc_arr(i-1,j,k-1)) );
594  }
595  }
596  gpx_arr(i,j,k) = dxInv[0] * (p_hi - p_lo);
597  } else {
598  gpx_arr(i,j,k) = dxInv[0] * (p_arr(i,j,k) - p_arr(i-1,j,k));
599  }
600 
601  // NOTE that the gradp array now carries the map factor!
602  gpx_arr(i,j,k) *= mf_ux_arr(i,j,0);
603  },
604  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
605  {
606  if (l_use_terrain_fitted_coords) {
607  Real p_lo = p_arr(i,j-1,k);
608  Real p_hi = p_arr(i,j,k);
609  Real dz_int = myhalf * (z_cc_arr(i,j,k) - z_cc_arr(i,j-1,k));
610  if (dz_int > 0) {
611  // Klemp 2011, Eqn. 16: s = 1/2
612  if (k==domain_klo) {
613  p_hi = quad_interp_1d(z_cc_arr(i,j,k) - dz_int,
614  z_cc_arr(i,j,k ), p_arr(i,j,k ),
615  z_cc_arr(i,j,k+1), p_arr(i,j,k+1),
616  z_cc_arr(i,j,k+2), p_arr(i,j,k+2));
617  } else {
618  p_hi -= dz_int * ( ( p_arr(i,j ,k ) - p_arr(i,j ,k-1))
619  / (z_cc_arr(i,j ,k ) - z_cc_arr(i,j ,k-1)) );
620  }
621  if (k==domain_khi) {
622  p_lo = quad_interp_1d(z_cc_arr(i,j-1,k) + dz_int,
623  z_cc_arr(i,j-1,k-2), p_arr(i,j-1,k-2),
624  z_cc_arr(i,j-1,k-1), p_arr(i,j-1,k-1),
625  z_cc_arr(i,j-1,k ), p_arr(i,j-1,k ));
626  } else {
627  p_lo += dz_int * ( ( p_arr(i,j-1,k+1) - p_arr(i,j-1,k ))
628  / (z_cc_arr(i,j-1,k+1) - z_cc_arr(i,j-1,k )) );
629  }
630  } else if (dz_int < 0) {
631  // Klemp 2011, Eqn. 16: s = -1/2
632  if (k==domain_khi) {
633  p_hi = quad_interp_1d(z_cc_arr(i,j,k) - dz_int,
634  z_cc_arr(i,j,k-2), p_arr(i,j,k-2),
635  z_cc_arr(i,j,k-1), p_arr(i,j,k-1),
636  z_cc_arr(i,j,k ), p_arr(i,j,k ));
637  } else {
638  p_hi -= dz_int * ( ( p_arr(i,j ,k+1) - p_arr(i,j ,k ))
639  / (z_cc_arr(i,j ,k+1) - z_cc_arr(i,j ,k )) );
640  }
641  if (k==domain_klo) {
642  p_lo = quad_interp_1d(z_cc_arr(i,j-1,k) + dz_int,
643  z_cc_arr(i,j-1,k ), p_arr(i,j-1,k ),
644  z_cc_arr(i,j-1,k+1), p_arr(i,j-1,k+1),
645  z_cc_arr(i,j-1,k+2), p_arr(i,j-1,k+2));
646  } else {
647  p_lo += dz_int * ( ( p_arr(i,j-1,k ) - p_arr(i,j-1,k-1))
648  / (z_cc_arr(i,j-1,k ) - z_cc_arr(i,j-1,k-1)) );
649  }
650  }
651  gpy_arr(i,j,k) = dxInv[1] * (p_hi - p_lo);
652  } else {
653  gpy_arr(i,j,k) = dxInv[1] * (p_arr(i,j,k) - p_arr(i,j-1,k));
654  }
655 
656  // NOTE that the gradp array now carries the map factor!
657  gpy_arr(i,j,k) *= mf_vy_arr(i,j,0);
658  },
659  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
660  {
661  // Note: identical to gradp_type == 0
662  Real met_h_zeta = (l_use_terrain_fitted_coords) ? Compute_h_zeta_AtKface(i, j, k, dxInv, z_nd_arr) : 1;
663  gpz_arr(i,j,k) = dxInv[2] * ( p_arr(i,j,k)-p_arr(i,j,k-1) ) / met_h_zeta;
664  });
665  } // mfi
666 }
myhalf
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:13
v_y
@ v_y
Definition: ERF_DataStruct.H:29
u_x
@ u_x
Definition: ERF_DataStruct.H:28
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);})
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:398
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:752
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:1848

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  } else {
256 
257  // Pressure gradients are fitted at the centroids of cut cells, if EB and Compressible.
258  // Least-Squares Fitting: Compute slope using 3x3x3 stencil
259 
260  const bool l_fitting = false;
261 
262  const Real* dx_arr = geom.CellSize();
263  const Real dx = dx_arr[0];
264  const Real dy = dx_arr[1];
265  const Real dz = dx_arr[2];
266 
267  // EB factory
268  Array4<const EBCellFlag> cellflg = (ebfact.get_const_factory())->getMultiEBCellFlagFab()[mfi].const_array();
269 
270  // EB u-factory
271  auto const* u_factory = ebfact.get_u_const_factory();
272  Array4<const EBCellFlag> u_cellflg = u_factory->getMultiEBCellFlagFab()[mfi].const_array();
273  Array4<const Real > u_volfrac = u_factory->getVolFrac().const_array(mfi);
274  bool u_is_cut = (u_factory->getMultiEBCellFlagFab()[mfi].getType() == FabType::singlevalued);
275  Array4<const Real > u_volcent = u_is_cut ? u_factory->getCentroid().const_array(mfi) : Array4<const Real>{};
276 
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  if (cellflg(i,j,k).isCovered()) {
335  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));
336  } else if (cellflg(i-1,j,k).isCovered()) {
337  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));
338  } else {
339  gpx_arr(i,j,k) = dxInv[0] * (p_arr(i,j,k) - p_arr(i-1,j,k));
340  }
341  } else {
342  gpx_arr(i,j,k) = zero;
343  }
344  },
345  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
346  {
347  if (v_volfrac(i,j,k) > zero) {
348  if (cellflg(i,j,k).isCovered()) {
349  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));
350  } else if (cellflg(i,j-1,k).isCovered()) {
351  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));
352  } else {
353  gpy_arr(i,j,k) = dxInv[1] * (p_arr(i,j,k) - p_arr(i,j-1,k));
354  }
355  } else {
356  gpy_arr(i,j,k) = zero;
357  }
358  });
359 
360  } // l_fitting
361 
362  } // TerrainType::EB
363 
364  } // mfi
365 }
three
constexpr amrex::Real three
Definition: ERF_Constants.H:11
two
constexpr amrex::Real two
Definition: ERF_Constants.H:10
zero
constexpr amrex::Real zero
Definition: ERF_Constants.H:8
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
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:121
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:119
eb_aux_::getMultiEBCellFlagFab
const amrex::FabArray< amrex::EBCellFlagFab > & getMultiEBCellFlagFab() const
Return the reconstructed EB cell flags.
Definition: ERF_EBAux.cpp:1145
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:270
SolverChoice::terrain_type
static TerrainType terrain_type
Terrain or immersed-boundary representation.
Definition: ERF_DataStruct.H:1839

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.
383 {
384  const bool l_use_terrain_fitted_coords = (solverChoice.mesh_type != MeshType::ConstantDz);
385 
386  const Box domain = geom.Domain();
387  const int domain_klo = domain.smallEnd(2);
388  const int domain_khi = domain.bigEnd(2);
389 
390  const GpuArray<Real, AMREX_SPACEDIM> dxInv = geom.InvCellSizeArray();
391 
392  // *****************************************************************************
393  // Take gradient of relevant quantity (p0, pres, or pert_pres = pres - p0)
394  // *****************************************************************************
395  for ( MFIter mfi(p); mfi.isValid(); ++mfi)
396  {
397  Box tbz = mfi.nodaltilebox(2);
398 
399  // We don't compute gpz on the bottom or top domain boundary
400  if (tbz.smallEnd(2) == domain_klo) {
401  tbz.growLo(2,-1);
402  }
403  if (tbz.bigEnd(2) == domain_khi+1) {
404  tbz.growHi(2,-1);
405  }
406 
407  // Terrain metrics
408  const Array4<const Real>& z_nd_arr = z_phys_nd.const_array(mfi);
409 
410  const Array4<const Real>& p_arr = p.const_array(mfi);
411 
412  const Array4< Real>& gpz_arr = gradp[GpVars::gpz].array(mfi);
413 
414  if (solverChoice.terrain_type != TerrainType::EB) {
415 
416  ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
417  {
418  Real met_h_zeta = (l_use_terrain_fitted_coords) ? Compute_h_zeta_AtKface(i, j, k, dxInv, z_nd_arr) : 1;
419  gpz_arr(i,j,k) = dxInv[2] * ( p_arr(i,j,k)-p_arr(i,j,k-1) ) / met_h_zeta;
420  });
421 
422  } else {
423 
424  // Pressure gradients are fitted at the centroids of cut cells, if EB and Compressible.
425  // Least-Squares Fitting: Compute slope using 3x3x3 stencil
426 
427  const bool l_fitting = false;
428 
429  const Real* dx_arr = geom.CellSize();
430  const Real dx = dx_arr[0];
431  const Real dy = dx_arr[1];
432  const Real dz = dx_arr[2];
433 
434  // EB factory
435  Array4<const EBCellFlag> cellflg = (ebfact.get_const_factory())->getMultiEBCellFlagFab()[mfi].const_array();
436 
437  // EB w-factory
438  auto const* w_factory = ebfact.get_w_const_factory();
439  Array4<const EBCellFlag> w_cellflg = w_factory->getMultiEBCellFlagFab()[mfi].const_array();
440  Array4<const Real > w_volfrac = w_factory->getVolFrac().const_array(mfi);
441  bool w_is_cut = (w_factory->getMultiEBCellFlagFab()[mfi].getType() == FabType::singlevalued);
442  Array4<const Real > w_volcent = w_is_cut ? w_factory->getCentroid().const_array(mfi) : Array4<Real>{};
443 
444  if (l_fitting) {
445 
446  ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k)
447  {
448  if (w_volfrac(i,j,k) > zero) {
449 
450  if (w_cellflg(i,j,k).isSingleValued()) {
451 
452  GpuArray<Real,AMREX_SPACEDIM> slopes;
453  slopes = erf_calc_slopes_eb_staggered(Vars::zvel, Vars::cons, dx, dy, dz, i, j, k, p_arr, w_volcent, w_cellflg);
454 
455  gpz_arr(i,j,k) = slopes[2];
456 
457  } else {
458  gpz_arr(i,j,k) = dxInv[2] * (p_arr(i,j,k) - p_arr(i,j,k-1));
459  }
460  } else {
461  gpz_arr(i,j,k) = zero;
462  }
463  });
464 
465  } else {
466 
467  // Simple calculation: assuming pressures at cell centers
468 
469  ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
470  {
471  if (w_volfrac(i,j,k) > zero) {
472  if (cellflg(i,j,k).isCovered()) {
473  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) );
474  } else if (cellflg(i,j,k-1).isCovered()) {
475  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) );
476  } else {
477  gpz_arr(i,j,k) = dxInv[2] * ( p_arr(i,j,k)-p_arr(i,j,k-1) );
478  }
479  } else {
480  gpz_arr(i,j,k) = zero;
481  }
482  });
483 
484  } // l_fitting
485 
486  } // TerrainType::EB
487 
488  } // mfi
489 }
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:123
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:502
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:2124
SolverChoice::use_pert_pres_gradient
bool use_pert_pres_gradient
Whether momentum equations use perturbational pressure gradients.
Definition: ERF_DataStruct.H:1909
SolverChoice::anelastic
amrex::Vector< int > anelastic
Per-level flag selecting anelastic dynamics.
Definition: ERF_DataStruct.H:1870
SolverChoice::gradp_type
int gradp_type
Terrain-fitted horizontal pressure-gradient formulation.
Definition: ERF_DataStruct.H:1907

Referenced by ERF::compute_max_buoyancy_gradp_diagnostic().

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