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

Function Documentation

◆ compute_gradp()

void compute_gradp ( const MultiFab &  p,
const Geometry &  geom,
MultiFab &  z_phys_nd,
MultiFab &  z_phys_cc,
const eb_ &  ebfact,
Vector< MultiFab > &  gradp,
const SolverChoice &  solverChoice 
)
87 {
88  const bool l_use_terrain_fitted_coords = (solverChoice.mesh_type != MeshType::ConstantDz);
89 
90  const Box domain = geom.Domain();
91  const int domain_klo = domain.smallEnd(2);
92  const int domain_khi = domain.bigEnd(2);
93 
94  const GpuArray<Real, AMREX_SPACEDIM> dxInv = geom.InvCellSizeArray();
95 
96  // *****************************************************************************
97  // Take gradient of relevant quantity (p0, pres, or pert_pres = pres - p0)
98  // *****************************************************************************
99  for ( MFIter mfi(p); mfi.isValid(); ++mfi)
100  {
101  Box tbx = mfi.nodaltilebox(0);
102  Box tby = mfi.nodaltilebox(1);
103  Box tbz = mfi.nodaltilebox(2);
104 
105  // We don't compute gpz on the bottom or top domain boundary
106  if (tbz.smallEnd(2) == domain_klo) {
107  tbz.growLo(2,-1);
108  }
109  if (tbz.bigEnd(2) == domain_khi+1) {
110  tbz.growHi(2,-1);
111  }
112 
113  // Terrain metrics
114  const Array4<const Real>& z_nd_arr = z_phys_nd.const_array(mfi);
115  const Array4<const Real>& z_cc_arr = z_phys_cc.const_array(mfi);
116 
117  const Array4<const Real>& p_arr = p.const_array(mfi);
118 
119  const Array4< Real>& gpx_arr = gradp[GpVars::gpx].array(mfi);
120  const Array4< Real>& gpy_arr = gradp[GpVars::gpy].array(mfi);
121  const Array4< Real>& gpz_arr = gradp[GpVars::gpz].array(mfi);
122 
123  if (solverChoice.terrain_type != TerrainType::EB) {
124 
125  ParallelFor(tbx, tby, tbz,
126  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
127  {
128  //Note : mx/my == 1, so no map factor needed here
129  Real gpx = dxInv[0] * (p_arr(i,j,k) - p_arr(i-1,j,k));
130 
131  if (l_use_terrain_fitted_coords) {
132  Real met_h_xi = (z_cc_arr(i,j,k) - z_cc_arr(i-1,j,k)) * dxInv[0];
133 
134  Real dz_phys_hi, dz_phys_lo;
135  Real gpz_lo, gpz_hi;
136  if (k==domain_klo) {
137  dz_phys_hi = z_cc_arr(i ,j,k+1) - z_cc_arr(i ,j,k );
138  dz_phys_lo = z_cc_arr(i-1,j,k+1) - z_cc_arr(i-1,j,k );
139  gpz_hi = (p_arr(i ,j,k+1) - p_arr(i ,j,k )) / dz_phys_hi;
140  gpz_lo = (p_arr(i-1,j,k+1) - p_arr(i-1,j,k )) / dz_phys_lo;
141  } else if (k==domain_khi) {
142  dz_phys_hi = z_cc_arr(i ,j,k ) - z_cc_arr(i ,j,k-1);
143  dz_phys_lo = z_cc_arr(i-1,j,k ) - z_cc_arr(i-1,j,k-1);
144  gpz_hi = (p_arr(i ,j,k ) - p_arr(i ,j,k-1)) / dz_phys_hi;
145  gpz_lo = (p_arr(i-1,j,k ) - p_arr(i-1,j,k-1)) / dz_phys_lo;
146  } else {
147  dz_phys_hi = z_cc_arr(i ,j,k+1) - z_cc_arr(i ,j,k-1);
148  dz_phys_lo = z_cc_arr(i-1,j,k+1) - z_cc_arr(i-1,j,k-1);
149  gpz_hi = (p_arr(i ,j,k+1) - p_arr(i ,j,k-1)) / dz_phys_hi;
150  gpz_lo = (p_arr(i-1,j,k+1) - p_arr(i-1,j,k-1)) / dz_phys_lo;
151  }
152  Real gpx_metric = met_h_xi * 0.5 * (gpz_hi + gpz_lo);
153  gpx -= gpx_metric;
154  }
155  gpx_arr(i,j,k) = gpx;
156  },
157  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
158  {
159  //Note : mx/my == 1, so no map factor needed here
160  Real gpy = dxInv[1] * (p_arr(i,j,k) - p_arr(i,j-1,k));
161 
162  if (l_use_terrain_fitted_coords) {
163  Real met_h_eta = (z_cc_arr(i,j,k) - z_cc_arr(i,j-1,k)) * dxInv[1];
164 
165  Real dz_phys_hi, dz_phys_lo;
166  Real gpz_lo, gpz_hi;
167  if (k==domain_klo) {
168  dz_phys_hi = z_cc_arr(i,j ,k+1) - z_cc_arr(i,j ,k );
169  dz_phys_lo = z_cc_arr(i,j-1,k+1) - z_cc_arr(i,j-1,k );
170  gpz_hi = (p_arr(i,j ,k+1) - p_arr(i,j ,k )) / dz_phys_hi;
171  gpz_lo = (p_arr(i,j-1,k+1) - p_arr(i,j-1,k )) / dz_phys_lo;
172  } else if (k==domain_khi) {
173  dz_phys_hi = z_cc_arr(i,j ,k ) - z_cc_arr(i,j ,k-1);
174  dz_phys_lo = z_cc_arr(i,j-1,k ) - z_cc_arr(i,j-1,k-1);
175  gpz_hi = (p_arr(i,j ,k ) - p_arr(i,j ,k-1)) / dz_phys_hi;
176  gpz_lo = (p_arr(i,j-1,k ) - p_arr(i,j-1,k-1)) / dz_phys_lo;
177  } else {
178  dz_phys_hi = z_cc_arr(i,j ,k+1) - z_cc_arr(i,j ,k-1);
179  dz_phys_lo = z_cc_arr(i,j-1,k+1) - z_cc_arr(i,j-1,k-1);
180  gpz_hi = (p_arr(i,j ,k+1) - p_arr(i,j ,k-1)) / dz_phys_hi;
181  gpz_lo = (p_arr(i,j-1,k+1) - p_arr(i,j-1,k-1)) / dz_phys_lo;
182  }
183  Real gpy_metric = met_h_eta * 0.5 * (gpz_hi + gpz_lo);
184  gpy -= gpy_metric;
185  }
186  gpy_arr(i,j,k) = gpy;
187  },
188  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
189  {
190  Real met_h_zeta = (l_use_terrain_fitted_coords) ? Compute_h_zeta_AtKface(i, j, k, dxInv, z_nd_arr) : 1;
191  gpz_arr(i,j,k) = dxInv[2] * ( p_arr(i,j,k)-p_arr(i,j,k-1) ) / met_h_zeta;
192  });
193 
194  } else {
195 
196  // Pressure gradients are fitted at the centroids of cut cells, if EB and Compressible.
197  // Least-Squares Fitting: Compute slope using 3x3x3 stencil
198 
199  const bool l_fitting = false;
200 
201  const Real* dx_arr = geom.CellSize();
202  const Real dx = dx_arr[0];
203  const Real dy = dx_arr[1];
204  const Real dz = dx_arr[2];
205 
206  // EB factory
207  Array4<const EBCellFlag> cellflg = (ebfact.get_const_factory())->getMultiEBCellFlagFab()[mfi].const_array();
208 
209  // EB u-factory
210  Array4<const EBCellFlag> u_cellflg = (ebfact.get_u_const_factory())->getMultiEBCellFlagFab()[mfi].const_array();
211  Array4<const Real > u_volfrac = (ebfact.get_u_const_factory())->getVolFrac().const_array(mfi);
212  Array4<const Real > u_volcent = (ebfact.get_u_const_factory())->getCentroid().const_array(mfi);
213 
214  // EB v-factory
215  Array4<const EBCellFlag> v_cellflg = (ebfact.get_v_const_factory())->getMultiEBCellFlagFab()[mfi].const_array();
216  Array4<const Real > v_volfrac = (ebfact.get_v_const_factory())->getVolFrac().const_array(mfi);
217  Array4<const Real > v_volcent = (ebfact.get_v_const_factory())->getCentroid().const_array(mfi);
218 
219  // EB w-factory
220  Array4<const EBCellFlag> w_cellflg = (ebfact.get_w_const_factory())->getMultiEBCellFlagFab()[mfi].const_array();
221  Array4<const Real > w_volfrac = (ebfact.get_w_const_factory())->getVolFrac().const_array(mfi);
222  Array4<const Real > w_volcent = (ebfact.get_w_const_factory())->getCentroid().const_array(mfi);
223 
224  if (l_fitting) {
225 
226  ParallelFor(tbx, tby, tbz,
227  [=] AMREX_GPU_DEVICE(int i, int j, int k)
228  {
229  if (u_volfrac(i,j,k) > 0.0) {
230 
231  if (u_cellflg(i,j,k).isSingleValued()) {
232 
233  GpuArray<Real,AMREX_SPACEDIM> slopes;
234  slopes = erf_calc_slopes_eb_staggered(Vars::xvel, Vars::cons, dx, dy, dz, i, j, k, p_arr, u_volcent, u_cellflg);
235 
236  gpx_arr(i,j,k) = slopes[0];
237 
238  } else {
239  gpx_arr(i,j,k) = dxInv[0] * (p_arr(i,j,k) - p_arr(i-1,j,k));
240  }
241 
242  } else {
243  gpx_arr(i,j,k) = 0.0;
244  }
245  },
246  [=] AMREX_GPU_DEVICE(int i, int j, int k)
247  {
248  if (v_volfrac(i,j,k) > 0.0) {
249 
250  if (v_cellflg(i,j,k).isSingleValued()) {
251 
252  GpuArray<Real,AMREX_SPACEDIM> slopes;
253  slopes = erf_calc_slopes_eb_staggered(Vars::yvel, Vars::cons, dx, dy, dz, i, j, k, p_arr, v_volcent, v_cellflg);
254 
255  gpy_arr(i,j,k) = slopes[1];
256 
257  } else {
258  gpy_arr(i,j,k) = dxInv[1] * (p_arr(i,j,k) - p_arr(i,j-1,k));
259  }
260  } else {
261  gpy_arr(i,j,k) = 0.0;
262  }
263  },
264  [=] AMREX_GPU_DEVICE(int i, int j, int k)
265  {
266  if (w_volfrac(i,j,k) > 0.0) {
267 
268  if (w_cellflg(i,j,k).isSingleValued()) {
269 
270  GpuArray<Real,AMREX_SPACEDIM> slopes;
271  slopes = erf_calc_slopes_eb_staggered(Vars::zvel, Vars::cons, dx, dy, dz, i, j, k, p_arr, w_volcent, w_cellflg);
272 
273  gpz_arr(i,j,k) = slopes[2];
274 
275  } else {
276  gpz_arr(i,j,k) = dxInv[2] * (p_arr(i,j,k) - p_arr(i,j,k-1));
277  }
278  } else {
279  gpz_arr(i,j,k) = 0.0;
280  }
281  });
282 
283  } else {
284 
285  // Simple calculation: assuming pressures at cell centers
286 
287  ParallelFor(tbx, tby, tbz,
288  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
289  {
290  if (u_volfrac(i,j,k) > 0.0) {
291  if (!cellflg(i-1,j,k).isCovered()) {
292  gpx_arr(i,j,k) = dxInv[0] * (p_arr(i,j,k) - p_arr(i-1,j,k));
293  } else {
294  if (!cellflg(i+1,j,k).isCovered()) {
295  gpx_arr(i,j,k) = dxInv[0] * (p_arr(i+1,j,k) - p_arr(i,j,k));
296  } else {
297  Abort("MakeGradP: both neighbors in x-direction are covered");
298  }
299  }
300  } else {
301  gpx_arr(i,j,k) = 0.0;
302  }
303  },
304  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
305  {
306  if (v_volfrac(i,j,k) > 0.0) {
307  if (!cellflg(i,j-1,k).isCovered()) {
308  gpy_arr(i,j,k) = dxInv[1] * (p_arr(i,j,k) - p_arr(i,j-1,k));
309  } else {
310  if (!cellflg(i,j+1,k).isCovered()) {
311  gpy_arr(i,j,k) = dxInv[1] * (p_arr(i,j+1,k) - p_arr(i,j,k));
312  } else {
313  Abort("MakeGradP: both neighbors in y-direction are covered");
314  }
315  }
316  } else {
317  gpy_arr(i,j,k) = 0.0;
318  }
319  },
320  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
321  {
322  if (w_volfrac(i,j,k) > 0.0) {
323  if (!cellflg(i,j,k-1).isCovered()) {
324  gpz_arr(i,j,k) = dxInv[2] * ( p_arr(i,j,k)-p_arr(i,j,k-1) );
325  } else {
326  if (!cellflg(i,j,k+1).isCovered()) {
327  gpz_arr(i,j,k) = dxInv[2] * ( p_arr(i,j,k+1)-p_arr(i,j,k) );
328  } else {
329  Abort("MakeGradP: both neighbors in z-direction are covered");
330  }
331  }
332  } else {
333  gpz_arr(i,j,k) = 0.0;
334  }
335  });
336 
337  } // l_fitting
338 
339  } // TerrainType::EB
340 
341  } // mfi
342 }
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:300
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:176
eb_::get_w_const_factory
eb_aux_ const * get_w_const_factory() const noexcept
Definition: ERF_EB.H:56
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:55
eb_::get_u_const_factory
eb_aux_ const * get_u_const_factory() const noexcept
Definition: ERF_EB.H:54
GpVars::gpz
@ gpz
Definition: ERF_IndexDefines.H:152
GpVars::gpy
@ gpy
Definition: ERF_IndexDefines.H:151
GpVars::gpx
@ gpx
Definition: ERF_IndexDefines.H:150
Vars::xvel
@ xvel
Definition: ERF_IndexDefines.H:141
Vars::cons
@ cons
Definition: ERF_IndexDefines.H:140
Vars::zvel
@ zvel
Definition: ERF_IndexDefines.H:143
Vars::yvel
@ yvel
Definition: ERF_IndexDefines.H:142
SolverChoice::mesh_type
static MeshType mesh_type
Definition: ERF_DataStruct.H:852
SolverChoice::terrain_type
static TerrainType terrain_type
Definition: ERF_DataStruct.H:843

Referenced by make_gradp_pert(), 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,
MultiFab &  z_phys_nd,
MultiFab &  z_phys_cc,
Vector< MultiFab > &  gradp,
const SolverChoice &  solverChoice 
)
351 {
352  const bool l_use_terrain_fitted_coords = (solverChoice.mesh_type != MeshType::ConstantDz);
353 
354  const Box domain = geom.Domain();
355  const int domain_klo = domain.smallEnd(2);
356  const int domain_khi = domain.bigEnd(2);
357 
358  const GpuArray<Real, AMREX_SPACEDIM> dxInv = geom.InvCellSizeArray();
359 
360  // *****************************************************************************
361  // Take gradient of relevant quantity (p0, pres, or pert_pres = pres - p0)
362  // *****************************************************************************
363  for ( MFIter mfi(p); mfi.isValid(); ++mfi)
364  {
365  Box tbx = mfi.nodaltilebox(0);
366  Box tby = mfi.nodaltilebox(1);
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  const Array4<const Real>& z_cc_arr = z_phys_cc.const_array(mfi);
380 
381  const Array4<const Real>& p_arr = p.const_array(mfi);
382 
383  const Array4< Real>& gpx_arr = gradp[GpVars::gpx].array(mfi);
384  const Array4< Real>& gpy_arr = gradp[GpVars::gpy].array(mfi);
385  const Array4< Real>& gpz_arr = gradp[GpVars::gpz].array(mfi);
386 
387  ParallelFor(tbx, tby, tbz,
388  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
389  {
390  if (l_use_terrain_fitted_coords) {
391  Real p_lo = p_arr(i-1,j,k);
392  Real p_hi = p_arr(i,j,k);
393  Real dz_int = 0.5 * (z_cc_arr(i,j,k) - z_cc_arr(i-1,j,k));
394  if (dz_int > 0) {
395  // Klemp 2011, Eqn. 16: s = 1/2
396  if (k==domain_klo) {
397  p_hi = quad_interp_1d(z_cc_arr(i,j,k) - dz_int,
398  z_cc_arr(i,j,k ), p_arr(i,j,k ),
399  z_cc_arr(i,j,k+1), p_arr(i,j,k+1),
400  z_cc_arr(i,j,k+2), p_arr(i,j,k+2));
401  } else {
402  p_hi -= dz_int * ( ( p_arr(i ,j,k ) - p_arr(i ,j,k-1))
403  / (z_cc_arr(i ,j,k ) - z_cc_arr(i ,j,k-1)) );
404  }
405  if (k==domain_khi) {
406  p_lo = quad_interp_1d(z_cc_arr(i-1,j,k) + dz_int,
407  z_cc_arr(i-1,j,k-2), p_arr(i-1,j,k-2),
408  z_cc_arr(i-1,j,k-1), p_arr(i-1,j,k-1),
409  z_cc_arr(i-1,j,k ), p_arr(i-1,j,k ));
410  } else {
411  p_lo += dz_int * ( ( p_arr(i-1,j,k+1) - p_arr(i-1,j,k ))
412  / (z_cc_arr(i-1,j,k+1) - z_cc_arr(i-1,j,k )) );
413  }
414  } else if (dz_int < 0) {
415  // Klemp 2011, Eqn. 16: s = -1/2
416  if (k==domain_khi) {
417  p_hi = quad_interp_1d(z_cc_arr(i,j,k) - dz_int,
418  z_cc_arr(i,j,k-2), p_arr(i,j,k-2),
419  z_cc_arr(i,j,k-1), p_arr(i,j,k-1),
420  z_cc_arr(i,j,k ), p_arr(i,j,k ));
421  } else {
422  p_hi -= dz_int * ( ( p_arr(i ,j,k+1) - p_arr(i ,j,k ))
423  / (z_cc_arr(i ,j,k+1) - z_cc_arr(i ,j,k )) );
424  }
425  if (k==domain_klo) {
426  p_lo = quad_interp_1d(z_cc_arr(i-1,j,k) + dz_int,
427  z_cc_arr(i-1,j,k ), p_arr(i-1,j,k ),
428  z_cc_arr(i-1,j,k+1), p_arr(i-1,j,k+1),
429  z_cc_arr(i-1,j,k+2), p_arr(i-1,j,k+2));
430  } else {
431  p_lo += dz_int * ( ( p_arr(i-1,j,k ) - p_arr(i-1,j,k-1))
432  / (z_cc_arr(i-1,j,k ) - z_cc_arr(i-1,j,k-1)) );
433  }
434  }
435  gpx_arr(i,j,k) = dxInv[0] * (p_hi - p_lo);
436  } else {
437  gpx_arr(i,j,k) = dxInv[0] * (p_arr(i,j,k) - p_arr(i-1,j,k));
438  }
439  },
440  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
441  {
442  if (l_use_terrain_fitted_coords) {
443  Real p_lo = p_arr(i,j-1,k);
444  Real p_hi = p_arr(i,j,k);
445  Real dz_int = 0.5 * (z_cc_arr(i,j,k) - z_cc_arr(i,j-1,k));
446  if (dz_int > 0) {
447  // Klemp 2011, Eqn. 16: s = 1/2
448  if (k==domain_klo) {
449  p_hi = quad_interp_1d(z_cc_arr(i,j,k) - dz_int,
450  z_cc_arr(i,j,k ), p_arr(i,j,k ),
451  z_cc_arr(i,j,k+1), p_arr(i,j,k+1),
452  z_cc_arr(i,j,k+2), p_arr(i,j,k+2));
453  } else {
454  p_hi -= dz_int * ( ( p_arr(i,j ,k ) - p_arr(i,j ,k-1))
455  / (z_cc_arr(i,j ,k ) - z_cc_arr(i,j ,k-1)) );
456  }
457  if (k==domain_khi) {
458  p_lo = quad_interp_1d(z_cc_arr(i,j-1,k) + dz_int,
459  z_cc_arr(i,j-1,k-2), p_arr(i,j-1,k-2),
460  z_cc_arr(i,j-1,k-1), p_arr(i,j-1,k-1),
461  z_cc_arr(i,j-1,k ), p_arr(i,j-1,k ));
462  } else {
463  p_lo += dz_int * ( ( p_arr(i,j-1,k+1) - p_arr(i,j-1,k ))
464  / (z_cc_arr(i,j-1,k+1) - z_cc_arr(i,j-1,k )) );
465  }
466  } else if (dz_int < 0) {
467  // Klemp 2011, Eqn. 16: s = -1/2
468  if (k==domain_khi) {
469  p_hi = quad_interp_1d(z_cc_arr(i,j,k) - dz_int,
470  z_cc_arr(i,j,k-2), p_arr(i,j,k-2),
471  z_cc_arr(i,j,k-1), p_arr(i,j,k-1),
472  z_cc_arr(i,j,k ), p_arr(i,j,k ));
473  } else {
474  p_hi -= dz_int * ( ( p_arr(i,j ,k+1) - p_arr(i,j ,k ))
475  / (z_cc_arr(i,j ,k+1) - z_cc_arr(i,j ,k )) );
476  }
477  if (k==domain_klo) {
478  p_lo = quad_interp_1d(z_cc_arr(i,j-1,k) + dz_int,
479  z_cc_arr(i,j-1,k ), p_arr(i,j-1,k ),
480  z_cc_arr(i,j-1,k+1), p_arr(i,j-1,k+1),
481  z_cc_arr(i,j-1,k+2), p_arr(i,j-1,k+2));
482  } else {
483  p_lo += dz_int * ( ( p_arr(i,j-1,k ) - p_arr(i,j-1,k-1))
484  / (z_cc_arr(i,j-1,k ) - z_cc_arr(i,j-1,k-1)) );
485  }
486  }
487  gpx_arr(i,j,k) = dxInv[1] * (p_hi - p_lo);
488  } else {
489  gpy_arr(i,j,k) = dxInv[1] * (p_arr(i,j,k) - p_arr(i,j-1,k));
490  }
491  },
492  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
493  {
494  // Note: identical to gradp_type == 0
495  Real met_h_zeta = (l_use_terrain_fitted_coords) ? Compute_h_zeta_AtKface(i, j, k, dxInv, z_nd_arr) : 1;
496  gpz_arr(i,j,k) = dxInv[2] * ( p_arr(i,j,k)-p_arr(i,j,k-1) ) / met_h_zeta;
497  });
498  } // mfi
499 }
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:454

Referenced by 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,
MultiFab &  p0,
MultiFab &  z_phys_nd,
MultiFab &  z_phys_cc,
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
37 {
38  const bool l_use_moisture = (solverChoice.moisture_type != MoistureType::None);
39  const bool l_eb_terrain = (solverChoice.terrain_type == TerrainType::EB);
40  //
41  // Note that we only recompute gradp if compressible;
42  // if anelastic then we have computed gradp in the projection
43  // and we can reuse it, no need to recompute it
44  //
45  if (solverChoice.anelastic[level] == 0)
46  {
47  const int ngrow = (l_eb_terrain) ? 2 : 1;
48  MultiFab p(S_data[Vars::cons].boxArray(), S_data[Vars::cons].DistributionMap(), 1, ngrow);
49 
50  // *****************************************************************************
51  // Compute pressure or perturbataional pressure
52  // *****************************************************************************
53  for ( MFIter mfi(S_data[Vars::cons]); mfi.isValid(); ++mfi)
54  {
55  Box gbx = mfi.tilebox();
56  gbx.grow(IntVect(ngrow,ngrow,ngrow));
57 
58  if (gbx.smallEnd(2) < 0) gbx.setSmall(2,0);
59  const Array4<const Real>& cell_data = S_data[Vars::cons].array(mfi);
60  const Array4<const Real>& p0_arr = p0.const_array(mfi);
61  const Array4< Real>& pptemp_arr = p.array(mfi);
62  ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
63  {
64  Real qv_for_p = (l_use_moisture) ? cell_data(i,j,k,RhoQ1_comp)/cell_data(i,j,k,Rho_comp) : 0.0;
65  pptemp_arr(i,j,k) = getPgivenRTh(cell_data(i,j,k,RhoTheta_comp),qv_for_p) - p0_arr(i,j,k);
66  });
67  }
68 
69  if (solverChoice.gradp_type == 0) {
70  compute_gradp(p,geom,z_phys_nd,z_phys_cc,ebfact,gradp,solverChoice);
71  } else if (solverChoice.gradp_type == 1) {
72  AMREX_ASSERT_WITH_MESSAGE(solverChoice.terrain_type != TerrainType::EB,
73  "gradp_type==1 not implemented for EB");
74  compute_gradp_interpz(p,geom,z_phys_nd,z_phys_cc,gradp,solverChoice);
75  }
76  }
77 }
getPgivenRTh
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getPgivenRTh(const amrex::Real rhotheta, const amrex::Real qv=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
void compute_gradp(const MultiFab &p, const Geometry &geom, MultiFab &z_phys_nd, MultiFab &z_phys_cc, const eb_ &ebfact, Vector< MultiFab > &gradp, const SolverChoice &solverChoice)
Definition: ERF_MakeGradP.cpp:80
compute_gradp_interpz
void compute_gradp_interpz(const MultiFab &p, const Geometry &geom, MultiFab &z_phys_nd, MultiFab &z_phys_cc, Vector< MultiFab > &gradp, const SolverChoice &solverChoice)
Definition: ERF_MakeGradP.cpp:345
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:868
SolverChoice::moisture_type
MoistureType moisture_type
Definition: ERF_DataStruct.H:958
SolverChoice::gradp_type
int gradp_type
Definition: ERF_DataStruct.H:881
Here is the call graph for this function: