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)
 

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

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

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  compute_gradp(p,geom,z_phys_nd,z_phys_cc,ebfact,gradp,solverChoice);
70  }
71 }
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:74
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:793
SolverChoice::moisture_type
MoistureType moisture_type
Definition: ERF_DataStruct.H:873
Here is the call graph for this function: