ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
ERF_AdvectionSrcForMom_N.H
Go to the documentation of this file.
1 #include "ERF_IndexDefines.H"
2 #include "ERF_Interpolation.H"
3 
4 /**
5  * Function for computing the advective tendency for the x-component of momentum
6  * without metric terms (constant OR stretched dz) and for higher-order stencils
7  *
8  * @tparam InterpType_H horizontal interpolation type
9  * @tparam InterpType_V vertical interpolation type
10  * @param[in] i i-index of x-face at which to create tendency
11  * @param[in] j j-index of x-face at which to create tendency
12  * @param[in] k k-index of x-face at which to create tendency
13  * @param[in] rho_u x-component of momentum
14  * @param[in] rho_v y-component of momentum
15  * @param[in] rho_w z-component of momentum
16  * @param[in] interp_u_h horizontal interpolator for x-component velocity
17  * @param[in] interp_u_v vertical interpolator for x-component velocity
18  * @param[in] mf_ux_inv inverse x map factor on x-faces
19  * @param[in] mf_uy_inv inverse y map factor on x-faces
20  * @param[in] mf_vx_inv inverse x map factor on y-faces
21  * @param[in] dxInv inverse grid spacing in x
22  * @param[in] dyInv inverse grid spacing in y
23  * @param[in] dzInv inverse grid spacing in z
24  * @return advective tendency for x-momentum
25  */
26 template<typename InterpType_H, typename InterpType_V>
27 AMREX_GPU_DEVICE
28 AMREX_FORCE_INLINE
29 amrex::Real
30 AdvectionSrcForXMom_N (int i, int j, int k,
31  const amrex::Array4<const amrex::Real>& rho_u,
32  const amrex::Array4<const amrex::Real>& rho_v,
33  const amrex::Array4<const amrex::Real>& rho_w,
34  InterpType_H interp_u_h,
35  InterpType_V interp_u_v,
36  const amrex::Array4<const amrex::Real>& mf_ux_inv,
37  const amrex::Array4<const amrex::Real>& mf_uy_inv,
38  const amrex::Array4<const amrex::Real>& mf_vx_inv,
39  const amrex::Real dxInv, const amrex::Real dyInv, const amrex::Real dzInv)
40 {
41  amrex::Real advectionSrc;
42 
43  amrex::Real rho_u_avg_lo, rho_u_avg_hi;
44  amrex::Real rho_v_avg_lo, rho_v_avg_hi;
45  amrex::Real rho_w_avg_lo, rho_w_avg_hi;
46 
47  amrex::Real xflux_hi; amrex::Real xflux_lo;
48  amrex::Real yflux_hi; amrex::Real yflux_lo;
49  amrex::Real zflux_hi; amrex::Real zflux_lo;
50 
51  amrex::Real interp_hi(zero), interp_lo(zero);
52 
53  rho_u_avg_hi = myhalf * (rho_u(i+1, j, k) * mf_uy_inv(i+1,j,0) + rho_u(i, j, k) * mf_uy_inv(i,j,0));
54  interp_u_h.InterpolateInX(i+1,j,k,0,interp_hi,rho_u_avg_hi);
55  xflux_hi = rho_u_avg_hi * interp_hi;
56 
57  rho_u_avg_lo = myhalf * (rho_u(i-1, j, k) * mf_uy_inv(i-1,j,0) + rho_u(i, j, k) * mf_uy_inv(i,j,0));
58  interp_u_h.InterpolateInX(i,j,k,0,interp_lo,rho_u_avg_lo);
59  xflux_lo = rho_u_avg_lo * interp_lo;
60 
61  rho_v_avg_hi = myhalf * (rho_v(i, j+1, k) * mf_vx_inv(i,j+1,0) + rho_v(i-1, j+1, k) * mf_vx_inv(i-1,j+1,0));
62  interp_u_h.InterpolateInY(i,j+1,k,0,interp_hi,rho_v_avg_hi);
63  yflux_hi = rho_v_avg_hi * interp_hi;
64 
65  rho_v_avg_lo = myhalf * (rho_v(i, j , k) * mf_vx_inv(i,j ,0) + rho_v(i-1, j , k) * mf_vx_inv(i-1,j ,0));
66  interp_u_h.InterpolateInY(i,j,k,0,interp_lo,rho_v_avg_lo);
67  yflux_lo = rho_v_avg_lo * interp_lo;
68 
69  rho_w_avg_hi = myhalf * (rho_w(i, j, k+1) + rho_w(i-1, j, k+1));
70  interp_u_v.InterpolateInZ(i,j,k+1,0,interp_hi,rho_w_avg_hi);
71  zflux_hi = rho_w_avg_hi * interp_hi;
72 
73  rho_w_avg_lo = myhalf * (rho_w(i, j, k ) + rho_w(i-1, j, k ));
74  interp_u_v.InterpolateInZ(i,j,k,0,interp_lo,rho_w_avg_lo);
75  zflux_lo = rho_w_avg_lo * interp_lo;
76 
77  amrex::Real mfsq = 1 / (mf_ux_inv(i,j,0) * mf_uy_inv(i,j,0));
78 
79  advectionSrc = (xflux_hi - xflux_lo) * dxInv * mfsq
80  + (yflux_hi - yflux_lo) * dyInv * mfsq
81  + (zflux_hi - zflux_lo) * dzInv;
82 
83  return advectionSrc;
84 }
85 
86 /**
87  * Function for computing the advective tendency for the y-component of momentum
88  * without metric terms and for higher-order stencils
89  *
90  * @tparam InterpType_H horizontal interpolation type
91  * @tparam InterpType_V vertical interpolation type
92  * @param[in] i i-index of y-face at which to create tendency
93  * @param[in] j j-index of y-face at which to create tendency
94  * @param[in] k k-index of y-face at which to create tendency
95  * @param[in] rho_u x-component of momentum
96  * @param[in] rho_v y-component of momentum
97  * @param[in] rho_w z-component of momentum
98  * @param[in] interp_v_h horizontal interpolator for y-component velocity
99  * @param[in] interp_v_v vertical interpolator for y-component velocity
100  * @param[in] mf_uy_inv inverse y map factor on x-faces
101  * @param[in] mf_vx_inv inverse x map factor on y-faces
102  * @param[in] mf_vy_inv inverse y map factor on y-faces
103  * @param[in] dxInv inverse grid spacing in x
104  * @param[in] dyInv inverse grid spacing in y
105  * @param[in] dzInv inverse grid spacing in z
106  * @return advective tendency for y-momentum
107  */
108 template<typename InterpType_H, typename InterpType_V>
109 AMREX_GPU_DEVICE
110 AMREX_FORCE_INLINE
111 amrex::Real
112 AdvectionSrcForYMom_N (int i, int j, int k,
113  const amrex::Array4<const amrex::Real>& rho_u,
114  const amrex::Array4<const amrex::Real>& rho_v,
115  const amrex::Array4<const amrex::Real>& rho_w,
116  InterpType_H interp_v_h,
117  InterpType_V interp_v_v,
118  const amrex::Array4<const amrex::Real>& mf_uy_inv,
119  const amrex::Array4<const amrex::Real>& mf_vx_inv,
120  const amrex::Array4<const amrex::Real>& mf_vy_inv,
121  const amrex::Real dxInv, const amrex::Real dyInv, const amrex::Real dzInv)
122 {
123  amrex::Real advectionSrc;
124 
125  amrex::Real rho_u_avg_lo, rho_u_avg_hi;
126  amrex::Real rho_v_avg_lo, rho_v_avg_hi;
127  amrex::Real rho_w_avg_lo, rho_w_avg_hi;
128 
129  amrex::Real xflux_hi; amrex::Real xflux_lo;
130  amrex::Real yflux_hi; amrex::Real yflux_lo;
131  amrex::Real zflux_hi; amrex::Real zflux_lo;
132 
133  amrex::Real interp_hi(zero), interp_lo(zero);
134 
135  rho_u_avg_hi = myhalf * (rho_u(i+1, j, k) * mf_uy_inv(i+1,j ,0) + rho_u(i+1, j-1, k) * mf_uy_inv(i+1,j-1,0));
136  interp_v_h.InterpolateInX(i+1,j,k,0,interp_hi,rho_u_avg_hi);
137  xflux_hi = rho_u_avg_hi * interp_hi;
138 
139  rho_u_avg_lo = myhalf * (rho_u(i , j, k) * mf_uy_inv(i ,j ,0) + rho_u(i , j-1, k) * mf_uy_inv(i ,j-1,0));
140  interp_v_h.InterpolateInX(i,j,k,0,interp_lo,rho_u_avg_lo);
141  xflux_lo = rho_u_avg_lo * interp_lo;
142 
143  rho_v_avg_hi = myhalf * (rho_v(i, j, k) * mf_vx_inv(i ,j ,0) + rho_v(i, j+1, k) * mf_vx_inv(i ,j+1,0));
144  interp_v_h.InterpolateInY(i,j+1,k,0,interp_hi,rho_v_avg_hi);
145  yflux_hi = rho_v_avg_hi * interp_hi;
146 
147  rho_v_avg_lo = myhalf * (rho_v(i, j, k) * mf_vx_inv(i ,j ,0) + rho_v(i, j-1, k) * mf_vx_inv(i ,j-1,0));
148  interp_v_h.InterpolateInY(i,j,k,0,interp_lo,rho_v_avg_lo);
149  yflux_lo = rho_v_avg_lo * interp_lo;
150 
151  rho_w_avg_hi = myhalf * (rho_w(i, j, k+1) + rho_w(i, j-1, k+1));
152  interp_v_v.InterpolateInZ(i,j,k+1,0,interp_hi,rho_w_avg_hi);
153  zflux_hi = rho_w_avg_hi * interp_hi;
154 
155  rho_w_avg_lo = myhalf * (rho_w(i, j, k ) + rho_w(i, j-1, k ));
156  interp_v_v.InterpolateInZ(i,j,k ,0,interp_lo,rho_w_avg_lo);
157  zflux_lo = rho_w_avg_lo * interp_lo;
158 
159  amrex::Real mfsq = 1 / (mf_vx_inv(i,j,0) * mf_vy_inv(i,j,0));
160 
161  advectionSrc = (xflux_hi - xflux_lo) * dxInv * mfsq
162  + (yflux_hi - yflux_lo) * dyInv * mfsq
163  + (zflux_hi - zflux_lo) * dzInv;
164 
165  return advectionSrc;
166 }
167 
168 /**
169  * Function for computing the advective tendency for the z-component of momentum
170  * without metric terms and for higher-order stencils
171  *
172  * @tparam InterpType_H horizontal interpolation type
173  * @tparam InterpType_V vertical interpolation type
174  * @tparam WallInterpType wall-aware vertical interpolation type
175  * @param[in] i i-index of z-face at which to create tendency
176  * @param[in] j j-index of z-face at which to create tendency
177  * @param[in] k k-index of z-face at which to create tendency
178  * @param[in] rho_u x-component of momentum
179  * @param[in] rho_v y-component of momentum
180  * @param[in] rho_w z-component of momentum
181  * @param[in] w z-component of velocity
182  * @param[in] interp_w_h horizontal interpolator for z-component velocity
183  * @param[in] interp_w_v vertical interpolator for z-component velocity
184  * @param[in] interp_w_wall wall-aware vertical interpolator for z-component velocity
185  * @param[in] mf_mx x map factor on cell centers
186  * @param[in] mf_my y map factor on cell centers
187  * @param[in] mf_uy_inv inverse y map factor on x-faces
188  * @param[in] mf_vx_inv inverse x map factor on y-faces
189  * @param[in] vert_adv_type vertical advection stencil
190  * @param[in] lo_z_face minimum z-face k-index at this level
191  * @param[in] hi_z_face maximum z-face k-index at this level
192  * @param[in] dxInv inverse grid spacing in x
193  * @param[in] dyInv inverse grid spacing in y
194  * @param[in] dzInv inverse grid spacing in z
195  * @return advective tendency for z-momentum
196  */
197 template<typename InterpType_H, typename InterpType_V, typename WallInterpType>
198 AMREX_GPU_DEVICE
199 AMREX_FORCE_INLINE
200 amrex::Real
201 AdvectionSrcForZMom_N (int i, int j, int k,
202  const amrex::Array4<const amrex::Real>& rho_u,
203  const amrex::Array4<const amrex::Real>& rho_v,
204  const amrex::Array4<const amrex::Real>& rho_w,
205  const amrex::Array4<const amrex::Real>& w,
206  InterpType_H interp_w_h,
207  InterpType_V interp_w_v,
208  WallInterpType interp_w_wall,
209  const amrex::Array4<const amrex::Real>& mf_mx,
210  const amrex::Array4<const amrex::Real>& mf_my,
211  const amrex::Array4<const amrex::Real>& mf_uy_inv,
212  const amrex::Array4<const amrex::Real>& mf_vx_inv,
213  const AdvType vert_adv_type,
214  const int lo_z_face, const int hi_z_face,
215  const amrex::Real dxInv, const amrex::Real dyInv, const amrex::Real dzInv)
216 {
217 
218  amrex::Real advectionSrc;
219 
220  amrex::Real rho_u_avg_lo, rho_u_avg_hi;
221  amrex::Real rho_v_avg_lo, rho_v_avg_hi;
222  amrex::Real rho_w_avg_lo, rho_w_avg_hi;
223 
224  amrex::Real xflux_hi; amrex::Real xflux_lo;
225  amrex::Real yflux_hi; amrex::Real yflux_lo;
226  amrex::Real zflux_hi; amrex::Real zflux_lo;
227 
228  amrex::Real interp_hi(zero), interp_lo(zero);
229 
230  rho_u_avg_hi = myhalf * (rho_u(i+1, j, k) + rho_u(i+1, j, k-1)) * mf_uy_inv(i+1,j ,0);
231  interp_w_h.InterpolateInX(i+1,j,k,0,interp_hi,rho_u_avg_hi);
232  xflux_hi = rho_u_avg_hi * interp_hi;
233 
234  rho_u_avg_lo = myhalf * (rho_u(i , j, k) + rho_u(i , j, k-1)) * mf_uy_inv(i ,j ,0);
235  interp_w_h.InterpolateInX(i,j,k,0,interp_lo,rho_u_avg_lo);
236  xflux_lo = rho_u_avg_lo * interp_lo;
237 
238  rho_v_avg_hi = myhalf * (rho_v(i, j+1, k) + rho_v(i, j+1, k-1)) * mf_vx_inv(i ,j+1,0);
239  interp_w_h.InterpolateInY(i,j+1,k,0,interp_hi,rho_v_avg_hi);
240  yflux_hi = rho_v_avg_hi * interp_hi;
241 
242  rho_v_avg_lo = myhalf * (rho_v(i, j , k) + rho_v(i, j , k-1)) * mf_vx_inv(i ,j ,0);
243  interp_w_h.InterpolateInY(i,j,k,0,interp_lo,rho_v_avg_lo);
244  yflux_lo = rho_v_avg_lo * interp_lo;
245 
246  // int l_spatial_order_hi = std::min(std::min(vert_spatial_order, 2*(hi_z_face-k)), 2*(k+1));
247  // If k == hi_z_face-1, l_spatial_order_hi = 2
248  // If k == hi_z_face-2, l_spatial_order_hi = std::min(vert_spatial_order, 4);
249  // If k == lo_z+face+1, l_spatial_order_hi = std::min(vert_spatial_order, 4);
250 
251  if (k == hi_z_face) {
252  zflux_hi = rho_w(i,j,k) * w(i,j,k);
253  } else {
254  rho_w_avg_hi = myhalf * (rho_w(i,j,k) + rho_w(i,j,k+1));
255  // The centre next to the top or bottom face is second order, seen from either side
256  if (k == hi_z_face-1 || k == lo_z_face) {
257  interp_w_wall.InterpolateInZ(i,j,k+1,0,interp_hi,rho_w_avg_hi,AdvType::Centered_2nd);
258  } else if (k == hi_z_face-2 || k == lo_z_face+1) {
259  if (vert_adv_type != AdvType::Centered_2nd && vert_adv_type != AdvType::Upwind_3rd) {
260  interp_w_wall.InterpolateInZ(i,j,k+1,0,interp_hi,rho_w_avg_hi,AdvType::Centered_4th);
261  } else {
262  interp_w_wall.InterpolateInZ(i,j,k+1,0,interp_hi,rho_w_avg_hi,vert_adv_type);
263  }
264  } else {
265  interp_w_v.InterpolateInZ(i,j,k+1,0,interp_hi,rho_w_avg_hi);
266  }
267  zflux_hi = rho_w_avg_hi * interp_hi;
268  }
269 
270  // int l_spatial_order_lo = std::min(std::min(vert_spatial_order, 2*(hi_z_face+1-k)), 2*k);
271  // If k == lo_z_face+1, l_spatial_order_hi = 2
272  // If k == lo_z_face+2, l_spatial_order_hi = std::min(vert_spatial_order, 4);
273  // If k == hi_z_face-1, l_spatial_order_hi = std::min(vert_spatial_order, 4);
274 
275  if (k == lo_z_face) {
276  zflux_lo = rho_w(i,j,k) * w(i,j,k);
277  } else {
278  rho_w_avg_lo = myhalf * (rho_w(i,j,k) + rho_w(i,j,k-1));
279  if (k == lo_z_face+1 || k == hi_z_face) {
280  interp_w_wall.InterpolateInZ(i,j,k,0,interp_lo,rho_w_avg_lo,AdvType::Centered_2nd);
281  } else if (k == lo_z_face+2 || k == hi_z_face-1) {
282  if (vert_adv_type != AdvType::Centered_2nd && vert_adv_type != AdvType::Upwind_3rd) {
283  interp_w_wall.InterpolateInZ(i,j,k,0,interp_lo,rho_w_avg_lo,AdvType::Centered_4th);
284  } else {
285  interp_w_wall.InterpolateInZ(i,j,k,0,interp_lo,rho_w_avg_lo,vert_adv_type);
286  }
287  } else {
288  interp_w_v.InterpolateInZ(i,j,k,0,interp_lo,rho_w_avg_lo);
289  }
290  zflux_lo = rho_w_avg_lo * interp_lo;
291  }
292 
293  amrex::Real mfsq = mf_mx(i,j,0) * mf_my(i,j,0);
294 
295  advectionSrc = (xflux_hi - xflux_lo) * dxInv * mfsq
296  + (yflux_hi - yflux_lo) * dyInv * mfsq
297  + (zflux_hi - zflux_lo) * dzInv;
298 
299  return advectionSrc;
300 }
301 
302 /**
303  * Wrapper function for computing the advective tendency w/ spatial order > two
304  *
305  * @tparam InterpType_H horizontal interpolation type
306  * @tparam InterpType_V vertical interpolation type
307  * @tparam WallInterpType wall-aware vertical interpolation type
308  * @param[in] bxx box over which x-momentum is updated
309  * @param[in] bxy box over which y-momentum is updated
310  * @param[in] bxz box over which z-momentum is updated
311  * @param[out] rho_u_rhs tendency for the x-momentum equation
312  * @param[out] rho_v_rhs tendency for the y-momentum equation
313  * @param[out] rho_w_rhs tendency for the z-momentum equation
314  * @param[in] rho_u x-component of momentum
315  * @param[in] rho_v y-component of momentum
316  * @param[in] rho_w z-component of momentum
317  * @param[in] u x-component of velocity
318  * @param[in] v y-component of velocity
319  * @param[in] w z-component of velocity
320  * @param[in] cellSizeInv inverse grid spacing
321  * @param[in] stretched_dz_d device vector of vertical cell spacing
322  * @param[in] mf_mx x map factor on cell centers
323  * @param[in] mf_ux_inv inverse x map factor on x-faces
324  * @param[in] mf_vx_inv inverse x map factor on y-faces
325  * @param[in] mf_my y map factor on cell centers
326  * @param[in] mf_uy_inv inverse y map factor on x-faces
327  * @param[in] mf_vy_inv inverse y map factor on y-faces
328  * @param[in] upw_frac_h horizontal upwind blending fraction
329  * @param[in] upw_frac_v vertical upwind blending fraction
330  * @param[in] vert_adv_type vertical advection stencil
331  * @param[in] lo_z_face minimum z-face k-index at this level
332  * @param[in] hi_z_face maximum z-face k-index at this level
333  */
334 template<typename InterpType_H, typename InterpType_V, typename WallInterpType>
335 void
336 AdvectionSrcForMomWrapper_N (const amrex::Box& bxx, const amrex::Box& bxy, const amrex::Box& bxz,
337  const amrex::Array4<amrex::Real>& rho_u_rhs,
338  const amrex::Array4<amrex::Real>& rho_v_rhs,
339  const amrex::Array4<amrex::Real>& rho_w_rhs,
340  const amrex::Array4<const amrex::Real>& rho_u,
341  const amrex::Array4<const amrex::Real>& rho_v,
342  const amrex::Array4<const amrex::Real>& rho_w,
343  const amrex::Array4<const amrex::Real>& u,
344  const amrex::Array4<const amrex::Real>& v,
345  const amrex::Array4<const amrex::Real>& w,
346  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
347  const amrex::Gpu::DeviceVector<amrex::Real>& stretched_dz_d,
348  const amrex::Array4<const amrex::Real>& mf_mx,
349  const amrex::Array4<const amrex::Real>& mf_ux_inv,
350  const amrex::Array4<const amrex::Real>& mf_vx_inv,
351  const amrex::Array4<const amrex::Real>& mf_my,
352  const amrex::Array4<const amrex::Real>& mf_uy_inv,
353  const amrex::Array4<const amrex::Real>& mf_vy_inv,
354  const amrex::Real upw_frac_h,
355  const amrex::Real upw_frac_v,
356  const AdvType vert_adv_type,
357  const int lo_z_face, const int hi_z_face)
358 {
359  // Instantiate the appropriate structs
360  InterpType_H interp_u_h(u, upw_frac_h); InterpType_V interp_u_v(u, upw_frac_v); // X-MOM
361  InterpType_H interp_v_h(v, upw_frac_h); InterpType_V interp_v_v(v, upw_frac_v); // Y-MOM
362  InterpType_H interp_w_h(w, upw_frac_h); InterpType_V interp_w_v(w, upw_frac_v); // Z-MOM
363  WallInterpType interp_w_wall(w, upw_frac_v); // Z-MOM @ wall
364 
365  auto dxInv = cellSizeInv[0];
366  auto dyInv = cellSizeInv[1];
367 
368  auto dz_ptr = stretched_dz_d.data();
369 
370  amrex::ParallelFor(bxx,
371  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
372  {
373  auto dzInv = one/dz_ptr[k];
374  rho_u_rhs(i, j, k) = -AdvectionSrcForXMom_N(i, j, k, rho_u, rho_v, rho_w,
375  interp_u_h, interp_u_v,
376  mf_ux_inv, mf_uy_inv, mf_vx_inv,
377  dxInv, dyInv, dzInv);
378  });
379 
380  amrex::ParallelFor(bxy,
381  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
382  {
383  auto dzInv = one/dz_ptr[k];
384  rho_v_rhs(i, j, k) = -AdvectionSrcForYMom_N(i, j, k, rho_u, rho_v, rho_w,
385  interp_v_h, interp_v_v,
386  mf_uy_inv, mf_vx_inv, mf_vy_inv,
387  dxInv, dyInv, dzInv);
388  });
389 
390  amrex::ParallelFor(bxz,
391  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
392  {
393  // dz_ptr holds cell widths, so the top face (k == hi_z_face) has only the cell below it
394  auto dzInv = (k == 0) ? one / dz_ptr[k] :
395  (k == hi_z_face) ? one / dz_ptr[k-1] : two/(dz_ptr[k] + dz_ptr[k-1]);
396  rho_w_rhs(i, j, k) = -AdvectionSrcForZMom_N(i, j, k, rho_u, rho_v, rho_w, w,
397  interp_w_h, interp_w_v, interp_w_wall,
398  mf_mx, mf_my, mf_uy_inv, mf_vx_inv,
399  vert_adv_type, lo_z_face, hi_z_face,
400  dxInv, dyInv, dzInv);
401  });
402 }
403 
404 /**
405  * Wrapper function for computing the advective tendency w/ spatial order > two
406  *
407  * @tparam InterpType_H horizontal interpolation type
408  * @param[in] bxx box over which x-momentum is updated
409  * @param[in] bxy box over which y-momentum is updated
410  * @param[in] bxz box over which z-momentum is updated
411  * @param[out] rho_u_rhs tendency for the x-momentum equation
412  * @param[out] rho_v_rhs tendency for the y-momentum equation
413  * @param[out] rho_w_rhs tendency for the z-momentum equation
414  * @param[in] rho_u x-component of momentum
415  * @param[in] rho_v y-component of momentum
416  * @param[in] rho_w z-component of momentum
417  * @param[in] u x-component of velocity
418  * @param[in] v y-component of velocity
419  * @param[in] w z-component of velocity
420  * @param[in] cellSizeInv inverse grid spacing
421  * @param[in] stretched_dz_d device vector of vertical cell spacing
422  * @param[in] mf_mx x map factor on cell centers
423  * @param[in] mf_ux_inv inverse x map factor on x-faces
424  * @param[in] mf_vx_inv inverse x map factor on y-faces
425  * @param[in] mf_my y map factor on cell centers
426  * @param[in] mf_uy_inv inverse y map factor on x-faces
427  * @param[in] mf_vy_inv inverse y map factor on y-faces
428  * @param[in] upw_frac_h horizontal upwind blending fraction
429  * @param[in] upw_frac_v vertical upwind blending fraction
430  * @param[in] vert_adv_type vertical advection stencil
431  * @param[in] lo_z_face minimum z-face k-index at this level
432  * @param[in] hi_z_face maximum z-face k-index at this level
433  */
434 template<typename InterpType_H>
435 void
436 AdvectionSrcForMomVert_N (const amrex::Box& bxx, const amrex::Box& bxy, const amrex::Box& bxz,
437  const amrex::Array4<amrex::Real>& rho_u_rhs,
438  const amrex::Array4<amrex::Real>& rho_v_rhs,
439  const amrex::Array4<amrex::Real>& rho_w_rhs,
440  const amrex::Array4<const amrex::Real>& rho_u,
441  const amrex::Array4<const amrex::Real>& rho_v,
442  const amrex::Array4<const amrex::Real>& rho_w,
443  const amrex::Array4<const amrex::Real>& u,
444  const amrex::Array4<const amrex::Real>& v,
445  const amrex::Array4<const amrex::Real>& w,
446  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
447  const amrex::Gpu::DeviceVector<amrex::Real>& stretched_dz_d,
448  const amrex::Array4<const amrex::Real>& mf_mx,
449  const amrex::Array4<const amrex::Real>& mf_ux_inv,
450  const amrex::Array4<const amrex::Real>& mf_vx_inv,
451  const amrex::Array4<const amrex::Real>& mf_my,
452  const amrex::Array4<const amrex::Real>& mf_uy_inv,
453  const amrex::Array4<const amrex::Real>& mf_vy_inv,
454  const amrex::Real upw_frac_h,
455  const amrex::Real upw_frac_v,
456  const AdvType vert_adv_type,
457  const int lo_z_face, const int hi_z_face)
458 {
459  if (vert_adv_type == AdvType::Centered_2nd) {
460  AdvectionSrcForMomWrapper_N<InterpType_H,CENTERED2,UPWINDALL>(bxx, bxy, bxz,
461  rho_u_rhs, rho_v_rhs, rho_w_rhs,
462  rho_u, rho_v, rho_w, u, v, w,
463  cellSizeInv, stretched_dz_d,
464  mf_mx, mf_ux_inv, mf_vx_inv,
465  mf_my, mf_uy_inv, mf_vy_inv,
466  upw_frac_h, upw_frac_v,
467  vert_adv_type,
468  lo_z_face, hi_z_face);
469  } else if (vert_adv_type == AdvType::Upwind_3rd) {
470  AdvectionSrcForMomWrapper_N<InterpType_H,UPWIND3,UPWINDALL>(bxx, bxy, bxz,
471  rho_u_rhs, rho_v_rhs, rho_w_rhs,
472  rho_u, rho_v, rho_w, u, v, w,
473  cellSizeInv, stretched_dz_d,
474  mf_mx, mf_ux_inv, mf_vx_inv,
475  mf_my, mf_uy_inv, mf_vy_inv,
476  upw_frac_h, upw_frac_v,
477  vert_adv_type,
478  lo_z_face, hi_z_face);
479  } else if (vert_adv_type == AdvType::Centered_4th) {
480  AdvectionSrcForMomWrapper_N<InterpType_H,CENTERED4,UPWINDALL>(bxx, bxy, bxz,
481  rho_u_rhs, rho_v_rhs, rho_w_rhs,
482  rho_u, rho_v, rho_w, u, v, w,
483  cellSizeInv, stretched_dz_d,
484  mf_mx, mf_ux_inv, mf_vx_inv,
485  mf_my, mf_uy_inv, mf_vy_inv,
486  upw_frac_h, upw_frac_v,
487  vert_adv_type,
488  lo_z_face, hi_z_face);
489  } else if (vert_adv_type == AdvType::Upwind_5th) {
490  AdvectionSrcForMomWrapper_N<InterpType_H,UPWIND5,UPWINDALL>(bxx, bxy, bxz,
491  rho_u_rhs, rho_v_rhs, rho_w_rhs,
492  rho_u, rho_v, rho_w, u, v, w,
493  cellSizeInv, stretched_dz_d,
494  mf_mx, mf_ux_inv, mf_vx_inv,
495  mf_my, mf_uy_inv, mf_vy_inv,
496  upw_frac_h, upw_frac_v,
497  vert_adv_type,
498  lo_z_face, hi_z_face);
499  } else if (vert_adv_type == AdvType::Centered_6th) {
500  AdvectionSrcForMomWrapper_N<InterpType_H,CENTERED6,UPWINDALL>(bxx, bxy, bxz,
501  rho_u_rhs, rho_v_rhs, rho_w_rhs,
502  rho_u, rho_v, rho_w, u, v, w,
503  cellSizeInv, stretched_dz_d,
504  mf_mx, mf_ux_inv, mf_vx_inv,
505  mf_my, mf_uy_inv, mf_vy_inv,
506  upw_frac_h, upw_frac_v,
507  vert_adv_type,
508  lo_z_face, hi_z_face);
509  } else if (vert_adv_type == AdvType::Weno_3) {
510  AdvectionSrcForMomWrapper_N<InterpType_H,WENO3,UPWINDALL>(bxx, bxy, bxz,
511  rho_u_rhs, rho_v_rhs, rho_w_rhs,
512  rho_u, rho_v, rho_w, u, v, w,
513  cellSizeInv, stretched_dz_d,
514  mf_mx, mf_ux_inv, mf_vx_inv,
515  mf_my, mf_uy_inv, mf_vy_inv,
516  upw_frac_h, upw_frac_v,
517  vert_adv_type,
518  lo_z_face, hi_z_face);
519  } else if (vert_adv_type == AdvType::Weno_3Z) {
520  AdvectionSrcForMomWrapper_N<InterpType_H,WENO_Z3,UPWINDALL>(bxx, bxy, bxz,
521  rho_u_rhs, rho_v_rhs, rho_w_rhs,
522  rho_u, rho_v, rho_w, u, v, w,
523  cellSizeInv, stretched_dz_d,
524  mf_mx, mf_ux_inv, mf_vx_inv,
525  mf_my, mf_uy_inv, mf_vy_inv,
526  upw_frac_h, upw_frac_v,
527  vert_adv_type,
528  lo_z_face, hi_z_face);
529  } else if (vert_adv_type == AdvType::Weno_3MZQ) {
530  AdvectionSrcForMomWrapper_N<InterpType_H,WENO_MZQ3,UPWINDALL>(bxx, bxy, bxz,
531  rho_u_rhs, rho_v_rhs, rho_w_rhs,
532  rho_u, rho_v, rho_w, u, v, w,
533  cellSizeInv, stretched_dz_d,
534  mf_mx, mf_ux_inv, mf_vx_inv,
535  mf_my, mf_uy_inv, mf_vy_inv,
536  upw_frac_h, upw_frac_v,
537  vert_adv_type,
538  lo_z_face, hi_z_face);
539  } else if (vert_adv_type == AdvType::Weno_5) {
540  AdvectionSrcForMomWrapper_N<InterpType_H,WENO5,UPWINDALL>(bxx, bxy, bxz,
541  rho_u_rhs, rho_v_rhs, rho_w_rhs,
542  rho_u, rho_v, rho_w, u, v, w,
543  cellSizeInv, stretched_dz_d,
544  mf_mx, mf_ux_inv, mf_vx_inv,
545  mf_my, mf_uy_inv, mf_vy_inv,
546  upw_frac_h, upw_frac_v,
547  vert_adv_type,
548  lo_z_face, hi_z_face);
549  } else if (vert_adv_type == AdvType::Weno_5Z) {
550  AdvectionSrcForMomWrapper_N<InterpType_H,WENO_Z5,UPWINDALL>(bxx, bxy, bxz,
551  rho_u_rhs, rho_v_rhs, rho_w_rhs,
552  rho_u, rho_v, rho_w, u, v, w,
553  cellSizeInv, stretched_dz_d,
554  mf_mx, mf_ux_inv, mf_vx_inv,
555  mf_my, mf_uy_inv, mf_vy_inv,
556  upw_frac_h, upw_frac_v,
557  vert_adv_type,
558  lo_z_face, hi_z_face);
559  } else if (vert_adv_type == AdvType::Weno_7) {
560  AdvectionSrcForMomWrapper_N<InterpType_H,WENO7,UPWINDALL>(bxx, bxy, bxz,
561  rho_u_rhs, rho_v_rhs, rho_w_rhs,
562  rho_u, rho_v, rho_w, u, v, w,
563  cellSizeInv, stretched_dz_d,
564  mf_mx, mf_ux_inv, mf_vx_inv,
565  mf_my, mf_uy_inv, mf_vy_inv,
566  upw_frac_h, upw_frac_v,
567  vert_adv_type,
568  lo_z_face, hi_z_face);
569  } else if (vert_adv_type == AdvType::Weno_7Z) {
570  AdvectionSrcForMomWrapper_N<InterpType_H,WENO_Z7,UPWINDALL>(bxx, bxy, bxz,
571  rho_u_rhs, rho_v_rhs, rho_w_rhs,
572  rho_u, rho_v, rho_w, u, v, w,
573  cellSizeInv, stretched_dz_d,
574  mf_mx, mf_ux_inv, mf_vx_inv,
575  mf_my, mf_uy_inv, mf_vy_inv,
576  upw_frac_h, upw_frac_v,
577  vert_adv_type,
578  lo_z_face, hi_z_face);
579  } else {
580  AMREX_ASSERT_WITH_MESSAGE(false, "Unknown advection scheme!");
581  }
582 }
AdvectionSrcForMomVert_N
void AdvectionSrcForMomVert_N(const amrex::Box &bxx, const amrex::Box &bxy, const amrex::Box &bxz, const amrex::Array4< amrex::Real > &rho_u_rhs, const amrex::Array4< amrex::Real > &rho_v_rhs, const amrex::Array4< amrex::Real > &rho_w_rhs, const amrex::Array4< const amrex::Real > &rho_u, const amrex::Array4< const amrex::Real > &rho_v, const amrex::Array4< const amrex::Real > &rho_w, const amrex::Array4< const amrex::Real > &u, const amrex::Array4< const amrex::Real > &v, const amrex::Array4< const amrex::Real > &w, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Gpu::DeviceVector< amrex::Real > &stretched_dz_d, const amrex::Array4< const amrex::Real > &mf_mx, const amrex::Array4< const amrex::Real > &mf_ux_inv, const amrex::Array4< const amrex::Real > &mf_vx_inv, const amrex::Array4< const amrex::Real > &mf_my, const amrex::Array4< const amrex::Real > &mf_uy_inv, const amrex::Array4< const amrex::Real > &mf_vy_inv, const amrex::Real upw_frac_h, const amrex::Real upw_frac_v, const AdvType vert_adv_type, const int lo_z_face, const int hi_z_face)
Definition: ERF_AdvectionSrcForMom_N.H:436
AdvectionSrcForYMom_N
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real AdvectionSrcForYMom_N(int i, int j, int k, const amrex::Array4< const amrex::Real > &rho_u, const amrex::Array4< const amrex::Real > &rho_v, const amrex::Array4< const amrex::Real > &rho_w, InterpType_H interp_v_h, InterpType_V interp_v_v, const amrex::Array4< const amrex::Real > &mf_uy_inv, const amrex::Array4< const amrex::Real > &mf_vx_inv, const amrex::Array4< const amrex::Real > &mf_vy_inv, const amrex::Real dxInv, const amrex::Real dyInv, const amrex::Real dzInv)
Definition: ERF_AdvectionSrcForMom_N.H:112
AdvectionSrcForMomWrapper_N
void AdvectionSrcForMomWrapper_N(const amrex::Box &bxx, const amrex::Box &bxy, const amrex::Box &bxz, const amrex::Array4< amrex::Real > &rho_u_rhs, const amrex::Array4< amrex::Real > &rho_v_rhs, const amrex::Array4< amrex::Real > &rho_w_rhs, const amrex::Array4< const amrex::Real > &rho_u, const amrex::Array4< const amrex::Real > &rho_v, const amrex::Array4< const amrex::Real > &rho_w, const amrex::Array4< const amrex::Real > &u, const amrex::Array4< const amrex::Real > &v, const amrex::Array4< const amrex::Real > &w, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Gpu::DeviceVector< amrex::Real > &stretched_dz_d, const amrex::Array4< const amrex::Real > &mf_mx, const amrex::Array4< const amrex::Real > &mf_ux_inv, const amrex::Array4< const amrex::Real > &mf_vx_inv, const amrex::Array4< const amrex::Real > &mf_my, const amrex::Array4< const amrex::Real > &mf_uy_inv, const amrex::Array4< const amrex::Real > &mf_vy_inv, const amrex::Real upw_frac_h, const amrex::Real upw_frac_v, const AdvType vert_adv_type, const int lo_z_face, const int hi_z_face)
Definition: ERF_AdvectionSrcForMom_N.H:336
AdvectionSrcForZMom_N
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real AdvectionSrcForZMom_N(int i, int j, int k, const amrex::Array4< const amrex::Real > &rho_u, const amrex::Array4< const amrex::Real > &rho_v, const amrex::Array4< const amrex::Real > &rho_w, const amrex::Array4< const amrex::Real > &w, InterpType_H interp_w_h, InterpType_V interp_w_v, WallInterpType interp_w_wall, const amrex::Array4< const amrex::Real > &mf_mx, const amrex::Array4< const amrex::Real > &mf_my, const amrex::Array4< const amrex::Real > &mf_uy_inv, const amrex::Array4< const amrex::Real > &mf_vx_inv, const AdvType vert_adv_type, const int lo_z_face, const int hi_z_face, const amrex::Real dxInv, const amrex::Real dyInv, const amrex::Real dzInv)
Definition: ERF_AdvectionSrcForMom_N.H:201
AdvectionSrcForXMom_N
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real AdvectionSrcForXMom_N(int i, int j, int k, const amrex::Array4< const amrex::Real > &rho_u, const amrex::Array4< const amrex::Real > &rho_v, const amrex::Array4< const amrex::Real > &rho_w, InterpType_H interp_u_h, InterpType_V interp_u_v, const amrex::Array4< const amrex::Real > &mf_ux_inv, const amrex::Array4< const amrex::Real > &mf_uy_inv, const amrex::Array4< const amrex::Real > &mf_vx_inv, const amrex::Real dxInv, const amrex::Real dyInv, const amrex::Real dzInv)
Definition: ERF_AdvectionSrcForMom_N.H:30
ERF_IndexDefines.H
AdvType
AdvType
Definition: ERF_IndexDefines.H:352
AdvType::Weno_7
@ Weno_7
AdvType::Upwind_3rd
@ Upwind_3rd
AdvType::Weno_7Z
@ Weno_7Z
AdvType::Centered_4th
@ Centered_4th
AdvType::Upwind_5th
@ Upwind_5th
AdvType::Centered_6th
@ Centered_6th
AdvType::Weno_3Z
@ Weno_3Z
AdvType::Weno_3
@ Weno_3
AdvType::Centered_2nd
@ Centered_2nd
AdvType::Weno_5Z
@ Weno_5Z
AdvType::Weno_5
@ Weno_5
AdvType::Weno_3MZQ
@ Weno_3MZQ
dxInv
amrex::GpuArray< Real, AMREX_SPACEDIM > dxInv
Definition: ERF_InitCustomPertVels_ParticleTests.H:17
ERF_Interpolation.H
ParallelFor
ParallelFor(fab_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);})
two
constexpr amrex::Real two
Definition: ERF_NumericalConstants.H:38
one
constexpr amrex::Real one
Definition: ERF_NumericalConstants.H:37
zero
constexpr amrex::Real zero
Definition: ERF_NumericalConstants.H:36
myhalf
constexpr amrex::Real myhalf
Definition: ERF_NumericalConstants.H:41
w
Real w
Definition: ERF_Plotfile2DInterpolator.cpp:23
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
AMREX_ASSERT_WITH_MESSAGE
AMREX_ASSERT_WITH_MESSAGE(wbar_cutoff_min > wbar_cutoff_max, "ERROR: wbar_cutoff_min < wbar_cutoff_max")