ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
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, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
120  {
121  //Note : mx/my == 1, so no map factor needed here
122  Real gpx = dxInv[0] * (p_arr(i,j,k) - p_arr(i-1,j,k));
123 
124  if (l_use_terrain_fitted_coords) {
125  Real met_h_xi = (z_cc_arr(i,j,k) - z_cc_arr(i-1,j,k)) * dxInv[0];
126 
127  Real dz_phys_hi, dz_phys_lo;
128  Real gpz_lo, gpz_hi;
129  if (k==domain_klo) {
130  dz_phys_hi = z_cc_arr(i ,j,k+1) - z_cc_arr(i ,j,k );
131  dz_phys_lo = z_cc_arr(i-1,j,k+1) - z_cc_arr(i-1,j,k );
132  gpz_hi = (p_arr(i ,j,k+1) - p_arr(i ,j,k )) / dz_phys_hi;
133  gpz_lo = (p_arr(i-1,j,k+1) - p_arr(i-1,j,k )) / dz_phys_lo;
134  } else if (k==domain_khi) {
135  dz_phys_hi = z_cc_arr(i ,j,k ) - z_cc_arr(i ,j,k-1);
136  dz_phys_lo = z_cc_arr(i-1,j,k ) - z_cc_arr(i-1,j,k-1);
137  gpz_hi = (p_arr(i ,j,k ) - p_arr(i ,j,k-1)) / dz_phys_hi;
138  gpz_lo = (p_arr(i-1,j,k ) - p_arr(i-1,j,k-1)) / dz_phys_lo;
139  } else {
140  dz_phys_hi = z_cc_arr(i ,j,k+1) - z_cc_arr(i ,j,k-1);
141  dz_phys_lo = z_cc_arr(i-1,j,k+1) - z_cc_arr(i-1,j,k-1);
142  gpz_hi = (p_arr(i ,j,k+1) - p_arr(i ,j,k-1)) / dz_phys_hi;
143  gpz_lo = (p_arr(i-1,j,k+1) - p_arr(i-1,j,k-1)) / dz_phys_lo;
144  }
145  Real gpx_metric = met_h_xi * 0.5 * (gpz_hi + gpz_lo);
146  gpx -= gpx_metric;
147  }
148  gpx_arr(i,j,k) = gpx;
149  });
150 
151  ParallelFor(tby, [=] 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 
183  ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
184  {
185  Real met_h_zeta = (l_use_terrain_fitted_coords) ? Compute_h_zeta_AtKface(i, j, k, dxInv, z_nd_arr) : 1;
186  gpz_arr(i,j,k) = dxInv[2] * ( p_arr(i,j,k)-p_arr(i,j,k-1) ) / met_h_zeta;
187  });
188 
189  } else {
190 
191  // Pressure gradients are fitted at the centroids of cut cells, if EB and Compressible.
192  // Least-Squares Fitting: Compute slope using 3x3x3 stencil
193 
194  const bool l_fitting = false;
195 
196  const Real* dx_arr = geom.CellSize();
197  const Real dx = dx_arr[0];
198  const Real dy = dx_arr[1];
199  const Real dz = dx_arr[2];
200 
201  // EB factory
202  Array4<const EBCellFlag> cellflg = (ebfact.get_const_factory())->getMultiEBCellFlagFab()[mfi].const_array();
203 
204  // EB u-factory
205  Array4<const EBCellFlag> u_cellflg = (ebfact.get_u_const_factory())->getMultiEBCellFlagFab()[mfi].const_array();
206  Array4<const Real > u_volfrac = (ebfact.get_u_const_factory())->getVolFrac().const_array(mfi);
207  Array4<const Real > u_volcent = (ebfact.get_u_const_factory())->getCentroid().const_array(mfi);
208 
209  // EB v-factory
210  Array4<const EBCellFlag> v_cellflg = (ebfact.get_v_const_factory())->getMultiEBCellFlagFab()[mfi].const_array();
211  Array4<const Real > v_volfrac = (ebfact.get_v_const_factory())->getVolFrac().const_array(mfi);
212  Array4<const Real > v_volcent = (ebfact.get_v_const_factory())->getCentroid().const_array(mfi);
213 
214  // EB w-factory
215  Array4<const EBCellFlag> w_cellflg = (ebfact.get_w_const_factory())->getMultiEBCellFlagFab()[mfi].const_array();
216  Array4<const Real > w_volfrac = (ebfact.get_w_const_factory())->getVolFrac().const_array(mfi);
217  Array4<const Real > w_volcent = (ebfact.get_w_const_factory())->getCentroid().const_array(mfi);
218 
219  if (l_fitting) {
220 
221  ParallelFor(tbx, tby, tbz,
222  [=] AMREX_GPU_DEVICE(int i, int j, int k)
223  {
224  if (u_volfrac(i,j,k) > 0.0) {
225 
226  if (u_cellflg(i,j,k).isSingleValued()) {
227 
228  GpuArray<Real,AMREX_SPACEDIM> slopes;
229  slopes = erf_calc_slopes_eb_staggered(Vars::xvel, Vars::cons, dx, dy, dz, i, j, k, p_arr, u_volcent, u_cellflg);
230 
231  gpx_arr(i,j,k) = slopes[0];
232 
233  } else {
234  gpx_arr(i,j,k) = dxInv[0] * (p_arr(i,j,k) - p_arr(i-1,j,k));
235  }
236 
237  } else {
238  gpx_arr(i,j,k) = 0.0;
239  }
240  },
241  [=] AMREX_GPU_DEVICE(int i, int j, int k)
242  {
243  if (v_volfrac(i,j,k) > 0.0) {
244 
245  if (v_cellflg(i,j,k).isSingleValued()) {
246 
247  GpuArray<Real,AMREX_SPACEDIM> slopes;
248  slopes = erf_calc_slopes_eb_staggered(Vars::yvel, Vars::cons, dx, dy, dz, i, j, k, p_arr, v_volcent, v_cellflg);
249 
250  gpy_arr(i,j,k) = slopes[1];
251 
252  } else {
253  gpy_arr(i,j,k) = dxInv[1] * (p_arr(i,j,k) - p_arr(i,j-1,k));
254  }
255  } else {
256  gpy_arr(i,j,k) = 0.0;
257  }
258  },
259  [=] AMREX_GPU_DEVICE(int i, int j, int k)
260  {
261  if (w_volfrac(i,j,k) > 0.0) {
262 
263  if (w_cellflg(i,j,k).isSingleValued()) {
264 
265  GpuArray<Real,AMREX_SPACEDIM> slopes;
266  slopes = erf_calc_slopes_eb_staggered(Vars::zvel, Vars::cons, dx, dy, dz, i, j, k, p_arr, w_volcent, w_cellflg);
267 
268  gpz_arr(i,j,k) = slopes[2];
269 
270  } else {
271  gpz_arr(i,j,k) = dxInv[2] * (p_arr(i,j,k) - p_arr(i,j,k-1));
272  }
273  } else {
274  gpz_arr(i,j,k) = 0.0;
275  }
276  });
277 
278  } else {
279 
280  // Simple calculation: assuming pressures at cell centers
281 
282  ParallelFor(tbx, [=] 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 
299  ParallelFor(tby, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
300  {
301  if (v_volfrac(i,j,k) > 0.0) {
302  if (!cellflg(i,j-1,k).isCovered()) {
303  gpy_arr(i,j,k) = dxInv[1] * (p_arr(i,j,k) - p_arr(i,j-1,k));
304  } else {
305  if (!cellflg(i,j+1,k).isCovered()) {
306  gpy_arr(i,j,k) = dxInv[1] * (p_arr(i,j+1,k) - p_arr(i,j,k));
307  } else {
308  Abort("MakeGradP: both neighbors in y-direction are covered");
309  }
310  }
311  } else {
312  gpy_arr(i,j,k) = 0.0;
313  }
314  });
315 
316  ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
317  {
318  if (w_volfrac(i,j,k) > 0.0) {
319  if (!cellflg(i,j,k-1).isCovered()) {
320  gpz_arr(i,j,k) = dxInv[2] * ( p_arr(i,j,k)-p_arr(i,j,k-1) );
321  } else {
322  if (!cellflg(i,j,k+1).isCovered()) {
323  gpz_arr(i,j,k) = dxInv[2] * ( p_arr(i,j,k+1)-p_arr(i,j,k) );
324  } else {
325  Abort("MakeGradP: both neighbors in z-direction are covered");
326  }
327  }
328  } else {
329  gpz_arr(i,j,k) = 0.0;
330  }
331  });
332 
333  } // l_fitting
334 
335  } // TerrainType::EB
336 
337  } // mfi
338 }
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
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:755
SolverChoice::terrain_type
static TerrainType terrain_type
Definition: ERF_DataStruct.H:749

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

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:771
SolverChoice::moisture_type
MoistureType moisture_type
Definition: ERF_DataStruct.H:848
Here is the call graph for this function: