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)
 
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)
 
void compute_gradp_z (const MultiFab &p, const Geometry &geom, const MultiFab &z_phys_nd, const eb_ &ebfact, Vector< MultiFab > &gradp, const SolverChoice &solverChoice)
 
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)
 

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 
)
121 {
122  compute_gradp_xy(p,geom,z_phys_cc,mapfac,ebfact,gradp,solverChoice);
123  compute_gradp_z(p,geom,z_phys_nd,ebfact,gradp,solverChoice);
124 }
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)
Definition: ERF_MakeGradP.cpp:347
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)
Definition: ERF_MakeGradP.cpp:127
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 
)
469 {
470  const bool l_use_terrain_fitted_coords = (solverChoice.mesh_type != MeshType::ConstantDz);
471 
472  const Box domain = geom.Domain();
473  const int domain_klo = domain.smallEnd(2);
474  const int domain_khi = domain.bigEnd(2);
475 
476  const GpuArray<Real, AMREX_SPACEDIM> dxInv = geom.InvCellSizeArray();
477 
478  // *****************************************************************************
479  // Take gradient of relevant quantity (p0, pres, or pert_pres = pres - p0)
480  // *****************************************************************************
481  for ( MFIter mfi(p); mfi.isValid(); ++mfi)
482  {
483  Box tbx = mfi.nodaltilebox(0);
484  Box tby = mfi.nodaltilebox(1);
485  Box tbz = mfi.nodaltilebox(2);
486 
487  // We don't compute gpz on the bottom or top domain boundary
488  if (tbz.smallEnd(2) == domain_klo) {
489  tbz.growLo(2,-1);
490  }
491  if (tbz.bigEnd(2) == domain_khi+1) {
492  tbz.growHi(2,-1);
493  }
494 
495  // Terrain metrics
496  const Array4<const Real>& z_nd_arr = z_phys_nd.const_array(mfi);
497  const Array4<const Real>& z_cc_arr = z_phys_cc.const_array(mfi);
498 
499  const Array4<const Real>& p_arr = p.const_array(mfi);
500 
501  const Array4< Real>& gpx_arr = gradp[GpVars::gpx].array(mfi);
502  const Array4< Real>& gpy_arr = gradp[GpVars::gpy].array(mfi);
503  const Array4< Real>& gpz_arr = gradp[GpVars::gpz].array(mfi);
504 
505  const Array4<const Real>& mf_ux_arr = mapfac[MapFacType::u_x]->const_array(mfi);
506  const Array4<const Real>& mf_vy_arr = mapfac[MapFacType::v_y]->const_array(mfi);
507 
508  ParallelFor(tbx, tby, tbz,
509  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
510  {
511  if (l_use_terrain_fitted_coords) {
512  Real p_lo = p_arr(i-1,j,k);
513  Real p_hi = p_arr(i,j,k);
514  Real dz_int = myhalf * (z_cc_arr(i,j,k) - z_cc_arr(i-1,j,k));
515  if (dz_int > 0) {
516  // Klemp 2011, Eqn. 16: s = 1/2
517  if (k==domain_klo) {
518  p_hi = quad_interp_1d(z_cc_arr(i,j,k) - dz_int,
519  z_cc_arr(i,j,k ), p_arr(i,j,k ),
520  z_cc_arr(i,j,k+1), p_arr(i,j,k+1),
521  z_cc_arr(i,j,k+2), p_arr(i,j,k+2));
522  } else {
523  p_hi -= dz_int * ( ( p_arr(i ,j,k ) - p_arr(i ,j,k-1))
524  / (z_cc_arr(i ,j,k ) - z_cc_arr(i ,j,k-1)) );
525  }
526  if (k==domain_khi) {
527  p_lo = quad_interp_1d(z_cc_arr(i-1,j,k) + dz_int,
528  z_cc_arr(i-1,j,k-2), p_arr(i-1,j,k-2),
529  z_cc_arr(i-1,j,k-1), p_arr(i-1,j,k-1),
530  z_cc_arr(i-1,j,k ), p_arr(i-1,j,k ));
531  } else {
532  p_lo += dz_int * ( ( p_arr(i-1,j,k+1) - p_arr(i-1,j,k ))
533  / (z_cc_arr(i-1,j,k+1) - z_cc_arr(i-1,j,k )) );
534  }
535  } else if (dz_int < 0) {
536  // Klemp 2011, Eqn. 16: s = -1/2
537  if (k==domain_khi) {
538  p_hi = quad_interp_1d(z_cc_arr(i,j,k) - dz_int,
539  z_cc_arr(i,j,k-2), p_arr(i,j,k-2),
540  z_cc_arr(i,j,k-1), p_arr(i,j,k-1),
541  z_cc_arr(i,j,k ), p_arr(i,j,k ));
542  } else {
543  p_hi -= dz_int * ( ( p_arr(i ,j,k+1) - p_arr(i ,j,k ))
544  / (z_cc_arr(i ,j,k+1) - z_cc_arr(i ,j,k )) );
545  }
546  if (k==domain_klo) {
547  p_lo = quad_interp_1d(z_cc_arr(i-1,j,k) + dz_int,
548  z_cc_arr(i-1,j,k ), p_arr(i-1,j,k ),
549  z_cc_arr(i-1,j,k+1), p_arr(i-1,j,k+1),
550  z_cc_arr(i-1,j,k+2), p_arr(i-1,j,k+2));
551  } else {
552  p_lo += dz_int * ( ( p_arr(i-1,j,k ) - p_arr(i-1,j,k-1))
553  / (z_cc_arr(i-1,j,k ) - z_cc_arr(i-1,j,k-1)) );
554  }
555  }
556  gpx_arr(i,j,k) = dxInv[0] * (p_hi - p_lo);
557  } else {
558  gpx_arr(i,j,k) = dxInv[0] * (p_arr(i,j,k) - p_arr(i-1,j,k));
559  }
560 
561  // NOTE that the gradp array now carries the map factor!
562  gpx_arr(i,j,k) *= mf_ux_arr(i,j,0);
563  },
564  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
565  {
566  if (l_use_terrain_fitted_coords) {
567  Real p_lo = p_arr(i,j-1,k);
568  Real p_hi = p_arr(i,j,k);
569  Real dz_int = myhalf * (z_cc_arr(i,j,k) - z_cc_arr(i,j-1,k));
570  if (dz_int > 0) {
571  // Klemp 2011, Eqn. 16: s = 1/2
572  if (k==domain_klo) {
573  p_hi = quad_interp_1d(z_cc_arr(i,j,k) - dz_int,
574  z_cc_arr(i,j,k ), p_arr(i,j,k ),
575  z_cc_arr(i,j,k+1), p_arr(i,j,k+1),
576  z_cc_arr(i,j,k+2), p_arr(i,j,k+2));
577  } else {
578  p_hi -= dz_int * ( ( p_arr(i,j ,k ) - p_arr(i,j ,k-1))
579  / (z_cc_arr(i,j ,k ) - z_cc_arr(i,j ,k-1)) );
580  }
581  if (k==domain_khi) {
582  p_lo = quad_interp_1d(z_cc_arr(i,j-1,k) + dz_int,
583  z_cc_arr(i,j-1,k-2), p_arr(i,j-1,k-2),
584  z_cc_arr(i,j-1,k-1), p_arr(i,j-1,k-1),
585  z_cc_arr(i,j-1,k ), p_arr(i,j-1,k ));
586  } else {
587  p_lo += dz_int * ( ( p_arr(i,j-1,k+1) - p_arr(i,j-1,k ))
588  / (z_cc_arr(i,j-1,k+1) - z_cc_arr(i,j-1,k )) );
589  }
590  } else if (dz_int < 0) {
591  // Klemp 2011, Eqn. 16: s = -1/2
592  if (k==domain_khi) {
593  p_hi = quad_interp_1d(z_cc_arr(i,j,k) - dz_int,
594  z_cc_arr(i,j,k-2), p_arr(i,j,k-2),
595  z_cc_arr(i,j,k-1), p_arr(i,j,k-1),
596  z_cc_arr(i,j,k ), p_arr(i,j,k ));
597  } else {
598  p_hi -= dz_int * ( ( p_arr(i,j ,k+1) - p_arr(i,j ,k ))
599  / (z_cc_arr(i,j ,k+1) - z_cc_arr(i,j ,k )) );
600  }
601  if (k==domain_klo) {
602  p_lo = quad_interp_1d(z_cc_arr(i,j-1,k) + dz_int,
603  z_cc_arr(i,j-1,k ), p_arr(i,j-1,k ),
604  z_cc_arr(i,j-1,k+1), p_arr(i,j-1,k+1),
605  z_cc_arr(i,j-1,k+2), p_arr(i,j-1,k+2));
606  } else {
607  p_lo += dz_int * ( ( p_arr(i,j-1,k ) - p_arr(i,j-1,k-1))
608  / (z_cc_arr(i,j-1,k ) - z_cc_arr(i,j-1,k-1)) );
609  }
610  }
611  gpy_arr(i,j,k) = dxInv[1] * (p_hi - p_lo);
612  } else {
613  gpy_arr(i,j,k) = dxInv[1] * (p_arr(i,j,k) - p_arr(i,j-1,k));
614  }
615 
616  // NOTE that the gradp array now carries the map factor!
617  gpy_arr(i,j,k) *= mf_vy_arr(i,j,0);
618  },
619  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
620  {
621  // Note: identical to gradp_type == 0
622  Real met_h_zeta = (l_use_terrain_fitted_coords) ? Compute_h_zeta_AtKface(i, j, k, dxInv, z_nd_arr) : 1;
623  gpz_arr(i,j,k) = dxInv[2] * ( p_arr(i,j,k)-p_arr(i,j,k-1) ) / met_h_zeta;
624  });
625  } // mfi
626 }
myhalf
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:13
v_y
@ v_y
Definition: ERF_DataStruct.H:28
u_x
@ u_x
Definition: ERF_DataStruct.H:27
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:184
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:510
GpVars::gpz
@ gpz
Definition: ERF_IndexDefines.H:188
GpVars::gpy
@ gpy
Definition: ERF_IndexDefines.H:187
GpVars::gpx
@ gpx
Definition: ERF_IndexDefines.H:186
SolverChoice::mesh_type
static MeshType mesh_type
Vertical mesh representation.
Definition: ERF_DataStruct.H:1377

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 
)
134 {
135  const bool l_use_terrain_fitted_coords = (solverChoice.mesh_type != MeshType::ConstantDz);
136 
137  const Box domain = geom.Domain();
138  const int domain_klo = domain.smallEnd(2);
139  const int domain_khi = domain.bigEnd(2);
140 
141  const GpuArray<Real, AMREX_SPACEDIM> dxInv = geom.InvCellSizeArray();
142 
143  // *****************************************************************************
144  // Take gradient of relevant quantity (p0, pres, or pert_pres = pres - p0)
145  // *****************************************************************************
146  for ( MFIter mfi(p); mfi.isValid(); ++mfi)
147  {
148  Box tbx = mfi.nodaltilebox(0);
149  Box tby = mfi.nodaltilebox(1);
150 
151  // Terrain metrics
152  const Array4<const Real>& z_cc_arr = z_phys_cc.const_array(mfi);
153 
154  const Array4<const Real>& p_arr = p.const_array(mfi);
155 
156  const Array4< Real>& gpx_arr = gradp[GpVars::gpx].array(mfi);
157  const Array4< Real>& gpy_arr = gradp[GpVars::gpy].array(mfi);
158 
159  const Array4<const Real>& mf_ux_arr = mapfac[MapFacType::u_x]->const_array(mfi);
160  const Array4<const Real>& mf_vy_arr = mapfac[MapFacType::v_y]->const_array(mfi);
161 
162  if (solverChoice.terrain_type != TerrainType::EB) {
163 
164  ParallelFor(tbx, tby,
165  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
166  {
167  //Note : mx/my == 1, so no map factor needed here
168  Real gpx = dxInv[0] * (p_arr(i,j,k) - p_arr(i-1,j,k));
169 
170  if (l_use_terrain_fitted_coords) {
171  Real met_h_xi = (z_cc_arr(i,j,k) - z_cc_arr(i-1,j,k)) * dxInv[0];
172 
173  Real dz_phys_hi, dz_phys_lo;
174  Real gpz_lo, gpz_hi;
175  if (k==domain_klo) {
176  dz_phys_hi = z_cc_arr(i ,j,k+1) - z_cc_arr(i ,j,k );
177  dz_phys_lo = z_cc_arr(i-1,j,k+1) - z_cc_arr(i-1,j,k );
178  gpz_hi = (p_arr(i ,j,k+1) - p_arr(i ,j,k )) / dz_phys_hi;
179  gpz_lo = (p_arr(i-1,j,k+1) - p_arr(i-1,j,k )) / dz_phys_lo;
180  } else if (k==domain_khi) {
181  dz_phys_hi = z_cc_arr(i ,j,k ) - z_cc_arr(i ,j,k-1);
182  dz_phys_lo = z_cc_arr(i-1,j,k ) - z_cc_arr(i-1,j,k-1);
183  gpz_hi = (p_arr(i ,j,k ) - p_arr(i ,j,k-1)) / dz_phys_hi;
184  gpz_lo = (p_arr(i-1,j,k ) - p_arr(i-1,j,k-1)) / dz_phys_lo;
185  } else {
186  dz_phys_hi = z_cc_arr(i ,j,k+1) - z_cc_arr(i ,j,k-1);
187  dz_phys_lo = z_cc_arr(i-1,j,k+1) - z_cc_arr(i-1,j,k-1);
188  gpz_hi = (p_arr(i ,j,k+1) - p_arr(i ,j,k-1)) / dz_phys_hi;
189  gpz_lo = (p_arr(i-1,j,k+1) - p_arr(i-1,j,k-1)) / dz_phys_lo;
190  }
191  Real gpx_metric = met_h_xi * myhalf * (gpz_hi + gpz_lo);
192  gpx -= gpx_metric;
193  }
194  gpx_arr(i,j,k) = gpx;
195 
196  // NOTE that the gradp array now carries the map factor!
197  gpx_arr(i,j,k) *= mf_ux_arr(i,j,0);
198  },
199  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
200  {
201  //Note : mx/my == 1, so no map factor needed here
202  Real gpy = dxInv[1] * (p_arr(i,j,k) - p_arr(i,j-1,k));
203 
204  if (l_use_terrain_fitted_coords) {
205  Real met_h_eta = (z_cc_arr(i,j,k) - z_cc_arr(i,j-1,k)) * dxInv[1];
206 
207  Real dz_phys_hi, dz_phys_lo;
208  Real gpz_lo, gpz_hi;
209  if (k==domain_klo) {
210  dz_phys_hi = z_cc_arr(i,j ,k+1) - z_cc_arr(i,j ,k );
211  dz_phys_lo = z_cc_arr(i,j-1,k+1) - z_cc_arr(i,j-1,k );
212  gpz_hi = (p_arr(i,j ,k+1) - p_arr(i,j ,k )) / dz_phys_hi;
213  gpz_lo = (p_arr(i,j-1,k+1) - p_arr(i,j-1,k )) / dz_phys_lo;
214  } else if (k==domain_khi) {
215  dz_phys_hi = z_cc_arr(i,j ,k ) - z_cc_arr(i,j ,k-1);
216  dz_phys_lo = z_cc_arr(i,j-1,k ) - z_cc_arr(i,j-1,k-1);
217  gpz_hi = (p_arr(i,j ,k ) - p_arr(i,j ,k-1)) / dz_phys_hi;
218  gpz_lo = (p_arr(i,j-1,k ) - p_arr(i,j-1,k-1)) / dz_phys_lo;
219  } else {
220  dz_phys_hi = z_cc_arr(i,j ,k+1) - z_cc_arr(i,j ,k-1);
221  dz_phys_lo = z_cc_arr(i,j-1,k+1) - z_cc_arr(i,j-1,k-1);
222  gpz_hi = (p_arr(i,j ,k+1) - p_arr(i,j ,k-1)) / dz_phys_hi;
223  gpz_lo = (p_arr(i,j-1,k+1) - p_arr(i,j-1,k-1)) / dz_phys_lo;
224  }
225  Real gpy_metric = met_h_eta * myhalf * (gpz_hi + gpz_lo);
226  gpy -= gpy_metric;
227  }
228  gpy_arr(i,j,k) = gpy;
229 
230  // NOTE that the gradp array now carries the map factor!
231  gpy_arr(i,j,k) *= mf_vy_arr(i,j,0);
232  });
233 
234  } else {
235 
236  // Pressure gradients are fitted at the centroids of cut cells, if EB and Compressible.
237  // Least-Squares Fitting: Compute slope using 3x3x3 stencil
238 
239  const bool l_fitting = false;
240 
241  const Real* dx_arr = geom.CellSize();
242  const Real dx = dx_arr[0];
243  const Real dy = dx_arr[1];
244  const Real dz = dx_arr[2];
245 
246  // EB factory
247  Array4<const EBCellFlag> cellflg = (ebfact.get_const_factory())->getMultiEBCellFlagFab()[mfi].const_array();
248 
249  // EB u-factory
250  auto const* u_factory = ebfact.get_u_const_factory();
251  Array4<const EBCellFlag> u_cellflg = u_factory->getMultiEBCellFlagFab()[mfi].const_array();
252  Array4<const Real > u_volfrac = u_factory->getVolFrac().const_array(mfi);
253  bool u_is_cut = (u_factory->getMultiEBCellFlagFab()[mfi].getType() == FabType::singlevalued);
254  Array4<const Real > u_volcent = u_is_cut ? u_factory->getCentroid().const_array(mfi) : Array4<const Real>{};
255 
256 
257  // EB v-factory
258  auto const* v_factory = ebfact.get_v_const_factory();
259  Array4<const EBCellFlag> v_cellflg = v_factory->getMultiEBCellFlagFab()[mfi].const_array();
260  Array4<const Real > v_volfrac = v_factory->getVolFrac().const_array(mfi);
261  bool v_is_cut = (v_factory->getMultiEBCellFlagFab()[mfi].getType() == FabType::singlevalued);
262  Array4<const Real > v_volcent = v_is_cut ? v_factory->getCentroid().const_array(mfi) : Array4<const Real>{};
263 
264  if (l_fitting) {
265 
266  ParallelFor(tbx, tby,
267  [=] AMREX_GPU_DEVICE(int i, int j, int k)
268  {
269  if (u_volfrac(i,j,k) > zero) {
270 
271  if (u_cellflg(i,j,k).isSingleValued()) {
272 
273  GpuArray<Real,AMREX_SPACEDIM> slopes;
274  slopes = erf_calc_slopes_eb_staggered(Vars::xvel, Vars::cons, dx, dy, dz, i, j, k, p_arr, u_volcent, u_cellflg);
275 
276  gpx_arr(i,j,k) = slopes[0];
277 
278  } else {
279  gpx_arr(i,j,k) = dxInv[0] * (p_arr(i,j,k) - p_arr(i-1,j,k));
280  }
281 
282  } else {
283  gpx_arr(i,j,k) = zero;
284  }
285  },
286  [=] AMREX_GPU_DEVICE(int i, int j, int k)
287  {
288  if (v_volfrac(i,j,k) > zero) {
289 
290  if (v_cellflg(i,j,k).isSingleValued()) {
291 
292  GpuArray<Real,AMREX_SPACEDIM> slopes;
293  slopes = erf_calc_slopes_eb_staggered(Vars::yvel, Vars::cons, dx, dy, dz, i, j, k, p_arr, v_volcent, v_cellflg);
294 
295  gpy_arr(i,j,k) = slopes[1];
296 
297  } else {
298  gpy_arr(i,j,k) = dxInv[1] * (p_arr(i,j,k) - p_arr(i,j-1,k));
299  }
300  } else {
301  gpy_arr(i,j,k) = zero;
302  }
303  });
304 
305  } else {
306 
307  // Simple calculation: assuming pressures at cell centers
308 
309  ParallelFor(tbx, tby,
310  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
311  {
312  if (u_volfrac(i,j,k) > zero) {
313  if (cellflg(i,j,k).isCovered()) {
314  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));
315  } else if (cellflg(i-1,j,k).isCovered()) {
316  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));
317  } else {
318  gpx_arr(i,j,k) = dxInv[0] * (p_arr(i,j,k) - p_arr(i-1,j,k));
319  }
320  } else {
321  gpx_arr(i,j,k) = zero;
322  }
323  },
324  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
325  {
326  if (v_volfrac(i,j,k) > zero) {
327  if (cellflg(i,j,k).isCovered()) {
328  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));
329  } else if (cellflg(i,j-1,k).isCovered()) {
330  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));
331  } else {
332  gpy_arr(i,j,k) = dxInv[1] * (p_arr(i,j,k) - p_arr(i,j-1,k));
333  }
334  } else {
335  gpy_arr(i,j,k) = zero;
336  }
337  });
338 
339  } // l_fitting
340 
341  } // TerrainType::EB
342 
343  } // mfi
344 }
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:24
dx
const Real dx
Definition: ERF_InitCustomPert_ABL.H:23
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:1149
Vars::xvel
@ xvel
Definition: ERF_IndexDefines.H:177
Vars::cons
@ cons
Definition: ERF_IndexDefines.H:176
Vars::yvel
@ yvel
Definition: ERF_IndexDefines.H:178
WSM6SedCellScratch::dz
@ dz
Definition: ERF_AdvanceWSM6.cpp:104
SolverChoice::terrain_type
static TerrainType terrain_type
Terrain or immersed-boundary representation.
Definition: ERF_DataStruct.H:1368

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 
)
353 {
354  const bool l_use_terrain_fitted_coords = (solverChoice.mesh_type != MeshType::ConstantDz);
355 
356  const Box domain = geom.Domain();
357  const int domain_klo = domain.smallEnd(2);
358  const int domain_khi = domain.bigEnd(2);
359 
360  const GpuArray<Real, AMREX_SPACEDIM> dxInv = geom.InvCellSizeArray();
361 
362  // *****************************************************************************
363  // Take gradient of relevant quantity (p0, pres, or pert_pres = pres - p0)
364  // *****************************************************************************
365  for ( MFIter mfi(p); mfi.isValid(); ++mfi)
366  {
367  Box tbz = mfi.nodaltilebox(2);
368 
369  // We don't compute gpz on the bottom or top domain boundary
370  if (tbz.smallEnd(2) == domain_klo) {
371  tbz.growLo(2,-1);
372  }
373  if (tbz.bigEnd(2) == domain_khi+1) {
374  tbz.growHi(2,-1);
375  }
376 
377  // Terrain metrics
378  const Array4<const Real>& z_nd_arr = z_phys_nd.const_array(mfi);
379 
380  const Array4<const Real>& p_arr = p.const_array(mfi);
381 
382  const Array4< Real>& gpz_arr = gradp[GpVars::gpz].array(mfi);
383 
384  if (solverChoice.terrain_type != TerrainType::EB) {
385 
386  ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
387  {
388  Real met_h_zeta = (l_use_terrain_fitted_coords) ? Compute_h_zeta_AtKface(i, j, k, dxInv, z_nd_arr) : 1;
389  gpz_arr(i,j,k) = dxInv[2] * ( p_arr(i,j,k)-p_arr(i,j,k-1) ) / met_h_zeta;
390  });
391 
392  } else {
393 
394  // Pressure gradients are fitted at the centroids of cut cells, if EB and Compressible.
395  // Least-Squares Fitting: Compute slope using 3x3x3 stencil
396 
397  const bool l_fitting = false;
398 
399  const Real* dx_arr = geom.CellSize();
400  const Real dx = dx_arr[0];
401  const Real dy = dx_arr[1];
402  const Real dz = dx_arr[2];
403 
404  // EB factory
405  Array4<const EBCellFlag> cellflg = (ebfact.get_const_factory())->getMultiEBCellFlagFab()[mfi].const_array();
406 
407  // EB w-factory
408  auto const* w_factory = ebfact.get_w_const_factory();
409  Array4<const EBCellFlag> w_cellflg = w_factory->getMultiEBCellFlagFab()[mfi].const_array();
410  Array4<const Real > w_volfrac = w_factory->getVolFrac().const_array(mfi);
411  bool w_is_cut = (w_factory->getMultiEBCellFlagFab()[mfi].getType() == FabType::singlevalued);
412  Array4<const Real > w_volcent = w_is_cut ? w_factory->getCentroid().const_array(mfi) : Array4<Real>{};
413 
414  if (l_fitting) {
415 
416  ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k)
417  {
418  if (w_volfrac(i,j,k) > zero) {
419 
420  if (w_cellflg(i,j,k).isSingleValued()) {
421 
422  GpuArray<Real,AMREX_SPACEDIM> slopes;
423  slopes = erf_calc_slopes_eb_staggered(Vars::zvel, Vars::cons, dx, dy, dz, i, j, k, p_arr, w_volcent, w_cellflg);
424 
425  gpz_arr(i,j,k) = slopes[2];
426 
427  } else {
428  gpz_arr(i,j,k) = dxInv[2] * (p_arr(i,j,k) - p_arr(i,j,k-1));
429  }
430  } else {
431  gpz_arr(i,j,k) = zero;
432  }
433  });
434 
435  } else {
436 
437  // Simple calculation: assuming pressures at cell centers
438 
439  ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
440  {
441  if (w_volfrac(i,j,k) > zero) {
442  if (cellflg(i,j,k).isCovered()) {
443  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) );
444  } else if (cellflg(i,j,k-1).isCovered()) {
445  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) );
446  } else {
447  gpz_arr(i,j,k) = dxInv[2] * ( p_arr(i,j,k)-p_arr(i,j,k-1) );
448  }
449  } else {
450  gpz_arr(i,j,k) = zero;
451  }
452  });
453 
454  } // l_fitting
455 
456  } // TerrainType::EB
457 
458  } // mfi
459 }
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:179

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:36
RhoTheta_comp
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:37
RhoQ1_comp
#define RhoQ1_comp
Definition: ERF_IndexDefines.H:42
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)
Definition: ERF_MakeGradP.cpp:462
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:1604
SolverChoice::use_pert_pres_gradient
bool use_pert_pres_gradient
Whether momentum equations use perturbational pressure gradients.
Definition: ERF_DataStruct.H:1438
SolverChoice::anelastic
amrex::Vector< int > anelastic
Per-level flag selecting anelastic dynamics.
Definition: ERF_DataStruct.H:1399
SolverChoice::gradp_type
int gradp_type
Terrain-fitted horizontal pressure-gradient formulation.
Definition: ERF_DataStruct.H:1436
Here is the call graph for this function: