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

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,
const MultiFab &  z_phys_nd,
const MultiFab &  z_phys_cc,
Vector< std::unique_ptr< MultiFab >> &  mapfac,
Vector< MultiFab > &  gradp,
const SolverChoice &  solverChoice 
)
357 {
358  const bool l_use_terrain_fitted_coords = (solverChoice.mesh_type != MeshType::ConstantDz);
359 
360  const Box domain = geom.Domain();
361  const int domain_klo = domain.smallEnd(2);
362  const int domain_khi = domain.bigEnd(2);
363 
364  const GpuArray<Real, AMREX_SPACEDIM> dxInv = geom.InvCellSizeArray();
365 
366  // *****************************************************************************
367  // Take gradient of relevant quantity (p0, pres, or pert_pres = pres - p0)
368  // *****************************************************************************
369  for ( MFIter mfi(p); mfi.isValid(); ++mfi)
370  {
371  Box tbx = mfi.nodaltilebox(0);
372  Box tby = mfi.nodaltilebox(1);
373  Box tbz = mfi.nodaltilebox(2);
374 
375  // We don't compute gpz on the bottom or top domain boundary
376  if (tbz.smallEnd(2) == domain_klo) {
377  tbz.growLo(2,-1);
378  }
379  if (tbz.bigEnd(2) == domain_khi+1) {
380  tbz.growHi(2,-1);
381  }
382 
383  // Terrain metrics
384  const Array4<const Real>& z_nd_arr = z_phys_nd.const_array(mfi);
385  const Array4<const Real>& z_cc_arr = z_phys_cc.const_array(mfi);
386 
387  const Array4<const Real>& p_arr = p.const_array(mfi);
388 
389  const Array4< Real>& gpx_arr = gradp[GpVars::gpx].array(mfi);
390  const Array4< Real>& gpy_arr = gradp[GpVars::gpy].array(mfi);
391  const Array4< Real>& gpz_arr = gradp[GpVars::gpz].array(mfi);
392 
393  const Array4<const Real>& mf_ux_arr = mapfac[MapFacType::u_x]->const_array(mfi);
394  const Array4<const Real>& mf_vy_arr = mapfac[MapFacType::v_y]->const_array(mfi);
395 
396  ParallelFor(tbx, tby, tbz,
397  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
398  {
399  if (l_use_terrain_fitted_coords) {
400  Real p_lo = p_arr(i-1,j,k);
401  Real p_hi = p_arr(i,j,k);
402  Real dz_int = 0.5 * (z_cc_arr(i,j,k) - z_cc_arr(i-1,j,k));
403  if (dz_int > 0) {
404  // Klemp 2011, Eqn. 16: s = 1/2
405  if (k==domain_klo) {
406  p_hi = quad_interp_1d(z_cc_arr(i,j,k) - dz_int,
407  z_cc_arr(i,j,k ), p_arr(i,j,k ),
408  z_cc_arr(i,j,k+1), p_arr(i,j,k+1),
409  z_cc_arr(i,j,k+2), p_arr(i,j,k+2));
410  } else {
411  p_hi -= dz_int * ( ( p_arr(i ,j,k ) - p_arr(i ,j,k-1))
412  / (z_cc_arr(i ,j,k ) - z_cc_arr(i ,j,k-1)) );
413  }
414  if (k==domain_khi) {
415  p_lo = quad_interp_1d(z_cc_arr(i-1,j,k) + dz_int,
416  z_cc_arr(i-1,j,k-2), p_arr(i-1,j,k-2),
417  z_cc_arr(i-1,j,k-1), p_arr(i-1,j,k-1),
418  z_cc_arr(i-1,j,k ), p_arr(i-1,j,k ));
419  } else {
420  p_lo += dz_int * ( ( p_arr(i-1,j,k+1) - p_arr(i-1,j,k ))
421  / (z_cc_arr(i-1,j,k+1) - z_cc_arr(i-1,j,k )) );
422  }
423  } else if (dz_int < 0) {
424  // Klemp 2011, Eqn. 16: s = -1/2
425  if (k==domain_khi) {
426  p_hi = quad_interp_1d(z_cc_arr(i,j,k) - dz_int,
427  z_cc_arr(i,j,k-2), p_arr(i,j,k-2),
428  z_cc_arr(i,j,k-1), p_arr(i,j,k-1),
429  z_cc_arr(i,j,k ), p_arr(i,j,k ));
430  } else {
431  p_hi -= dz_int * ( ( p_arr(i ,j,k+1) - p_arr(i ,j,k ))
432  / (z_cc_arr(i ,j,k+1) - z_cc_arr(i ,j,k )) );
433  }
434  if (k==domain_klo) {
435  p_lo = quad_interp_1d(z_cc_arr(i-1,j,k) + dz_int,
436  z_cc_arr(i-1,j,k ), p_arr(i-1,j,k ),
437  z_cc_arr(i-1,j,k+1), p_arr(i-1,j,k+1),
438  z_cc_arr(i-1,j,k+2), p_arr(i-1,j,k+2));
439  } else {
440  p_lo += dz_int * ( ( p_arr(i-1,j,k ) - p_arr(i-1,j,k-1))
441  / (z_cc_arr(i-1,j,k ) - z_cc_arr(i-1,j,k-1)) );
442  }
443  }
444  gpx_arr(i,j,k) = dxInv[0] * (p_hi - p_lo);
445  } else {
446  gpx_arr(i,j,k) = dxInv[0] * (p_arr(i,j,k) - p_arr(i-1,j,k));
447  }
448 
449  // NOTE that the gradp array now carries the map factor!
450  gpx_arr(i,j,k) *= mf_ux_arr(i,j,0);
451  },
452  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
453  {
454  if (l_use_terrain_fitted_coords) {
455  Real p_lo = p_arr(i,j-1,k);
456  Real p_hi = p_arr(i,j,k);
457  Real dz_int = 0.5 * (z_cc_arr(i,j,k) - z_cc_arr(i,j-1,k));
458  if (dz_int > 0) {
459  // Klemp 2011, Eqn. 16: s = 1/2
460  if (k==domain_klo) {
461  p_hi = quad_interp_1d(z_cc_arr(i,j,k) - dz_int,
462  z_cc_arr(i,j,k ), p_arr(i,j,k ),
463  z_cc_arr(i,j,k+1), p_arr(i,j,k+1),
464  z_cc_arr(i,j,k+2), p_arr(i,j,k+2));
465  } else {
466  p_hi -= dz_int * ( ( p_arr(i,j ,k ) - p_arr(i,j ,k-1))
467  / (z_cc_arr(i,j ,k ) - z_cc_arr(i,j ,k-1)) );
468  }
469  if (k==domain_khi) {
470  p_lo = quad_interp_1d(z_cc_arr(i,j-1,k) + dz_int,
471  z_cc_arr(i,j-1,k-2), p_arr(i,j-1,k-2),
472  z_cc_arr(i,j-1,k-1), p_arr(i,j-1,k-1),
473  z_cc_arr(i,j-1,k ), p_arr(i,j-1,k ));
474  } else {
475  p_lo += dz_int * ( ( p_arr(i,j-1,k+1) - p_arr(i,j-1,k ))
476  / (z_cc_arr(i,j-1,k+1) - z_cc_arr(i,j-1,k )) );
477  }
478  } else if (dz_int < 0) {
479  // Klemp 2011, Eqn. 16: s = -1/2
480  if (k==domain_khi) {
481  p_hi = quad_interp_1d(z_cc_arr(i,j,k) - dz_int,
482  z_cc_arr(i,j,k-2), p_arr(i,j,k-2),
483  z_cc_arr(i,j,k-1), p_arr(i,j,k-1),
484  z_cc_arr(i,j,k ), p_arr(i,j,k ));
485  } else {
486  p_hi -= dz_int * ( ( p_arr(i,j ,k+1) - p_arr(i,j ,k ))
487  / (z_cc_arr(i,j ,k+1) - z_cc_arr(i,j ,k )) );
488  }
489  if (k==domain_klo) {
490  p_lo = quad_interp_1d(z_cc_arr(i,j-1,k) + dz_int,
491  z_cc_arr(i,j-1,k ), p_arr(i,j-1,k ),
492  z_cc_arr(i,j-1,k+1), p_arr(i,j-1,k+1),
493  z_cc_arr(i,j-1,k+2), p_arr(i,j-1,k+2));
494  } else {
495  p_lo += dz_int * ( ( p_arr(i,j-1,k ) - p_arr(i,j-1,k-1))
496  / (z_cc_arr(i,j-1,k ) - z_cc_arr(i,j-1,k-1)) );
497  }
498  }
499  gpx_arr(i,j,k) = dxInv[1] * (p_hi - p_lo);
500  } else {
501  gpy_arr(i,j,k) = dxInv[1] * (p_arr(i,j,k) - p_arr(i,j-1,k));
502  }
503 
504  // NOTE that the gradp array now carries the map factor!
505  gpy_arr(i,j,k) *= mf_vy_arr(i,j,0);
506  },
507  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
508  {
509  // Note: identical to gradp_type == 0
510  Real met_h_zeta = (l_use_terrain_fitted_coords) ? Compute_h_zeta_AtKface(i, j, k, dxInv, z_nd_arr) : 1;
511  gpz_arr(i,j,k) = dxInv[2] * ( p_arr(i,j,k)-p_arr(i,j,k-1) ) / met_h_zeta;
512  });
513  } // mfi
514 }
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:463

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