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:342
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:176

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

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  Array4<const EBCellFlag> u_cellflg = (ebfact.get_u_const_factory())->getMultiEBCellFlagFab()[mfi].const_array();
251  Array4<const Real > u_volfrac = (ebfact.get_u_const_factory())->getVolFrac().const_array(mfi);
252  Array4<const Real > u_volcent = (ebfact.get_u_const_factory())->getCentroid().const_array(mfi);
253 
254  // EB v-factory
255  Array4<const EBCellFlag> v_cellflg = (ebfact.get_v_const_factory())->getMultiEBCellFlagFab()[mfi].const_array();
256  Array4<const Real > v_volfrac = (ebfact.get_v_const_factory())->getVolFrac().const_array(mfi);
257  Array4<const Real > v_volcent = (ebfact.get_v_const_factory())->getCentroid().const_array(mfi);
258 
259  if (l_fitting) {
260 
261  ParallelFor(tbx, tby,
262  [=] AMREX_GPU_DEVICE(int i, int j, int k)
263  {
264  if (u_volfrac(i,j,k) > zero) {
265 
266  if (u_cellflg(i,j,k).isSingleValued()) {
267 
268  GpuArray<Real,AMREX_SPACEDIM> slopes;
269  slopes = erf_calc_slopes_eb_staggered(Vars::xvel, Vars::cons, dx, dy, dz, i, j, k, p_arr, u_volcent, u_cellflg);
270 
271  gpx_arr(i,j,k) = slopes[0];
272 
273  } else {
274  gpx_arr(i,j,k) = dxInv[0] * (p_arr(i,j,k) - p_arr(i-1,j,k));
275  }
276 
277  } else {
278  gpx_arr(i,j,k) = zero;
279  }
280  },
281  [=] AMREX_GPU_DEVICE(int i, int j, int k)
282  {
283  if (v_volfrac(i,j,k) > zero) {
284 
285  if (v_cellflg(i,j,k).isSingleValued()) {
286 
287  GpuArray<Real,AMREX_SPACEDIM> slopes;
288  slopes = erf_calc_slopes_eb_staggered(Vars::yvel, Vars::cons, dx, dy, dz, i, j, k, p_arr, v_volcent, v_cellflg);
289 
290  gpy_arr(i,j,k) = slopes[1];
291 
292  } else {
293  gpy_arr(i,j,k) = dxInv[1] * (p_arr(i,j,k) - p_arr(i,j-1,k));
294  }
295  } else {
296  gpy_arr(i,j,k) = zero;
297  }
298  });
299 
300  } else {
301 
302  // Simple calculation: assuming pressures at cell centers
303 
304  ParallelFor(tbx, tby,
305  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
306  {
307  if (u_volfrac(i,j,k) > zero) {
308  if (cellflg(i,j,k).isCovered()) {
309  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));
310  } else if (cellflg(i-1,j,k).isCovered()) {
311  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));
312  } else {
313  gpx_arr(i,j,k) = dxInv[0] * (p_arr(i,j,k) - p_arr(i-1,j,k));
314  }
315  } else {
316  gpx_arr(i,j,k) = zero;
317  }
318  },
319  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
320  {
321  if (v_volfrac(i,j,k) > zero) {
322  if (cellflg(i,j,k).isCovered()) {
323  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));
324  } else if (cellflg(i,j-1,k).isCovered()) {
325  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));
326  } else {
327  gpy_arr(i,j,k) = dxInv[1] * (p_arr(i,j,k) - p_arr(i,j-1,k));
328  }
329  } else {
330  gpy_arr(i,j,k) = zero;
331  }
332  });
333 
334  } // l_fitting
335 
336  } // TerrainType::EB
337 
338  } // mfi
339 }
three
constexpr amrex::Real three
Definition: ERF_Constants.H:9
two
constexpr amrex::Real two
Definition: ERF_Constants.H:8
zero
constexpr amrex::Real zero
Definition: ERF_Constants.H:6
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)
Definition: ERF_EBSlopes.H:289
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
Definition: ERF_EB.H:46
eb_::get_v_const_factory
eb_aux_ const * get_v_const_factory() const noexcept
Definition: ERF_EB.H:51
eb_::get_u_const_factory
eb_aux_ const * get_u_const_factory() const noexcept
Definition: ERF_EB.H:50
Vars::xvel
@ xvel
Definition: ERF_IndexDefines.H:175
Vars::cons
@ cons
Definition: ERF_IndexDefines.H:174
Vars::yvel
@ yvel
Definition: ERF_IndexDefines.H:176
WSM6SedCellScratch::dz
@ dz
Definition: ERF_AdvanceWSM6.cpp:104
SolverChoice::terrain_type
static TerrainType terrain_type
Definition: ERF_DataStruct.H:1125

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

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:455
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::anelastic
amrex::Vector< int > anelastic
Definition: ERF_DataStruct.H:1152
SolverChoice::use_pert_pres_gradient
bool use_pert_pres_gradient
Definition: ERF_DataStruct.H:1189
SolverChoice::moisture_type
MoistureType moisture_type
Definition: ERF_DataStruct.H:1299
SolverChoice::gradp_type
int gradp_type
Definition: ERF_DataStruct.H:1187
Here is the call graph for this function: