ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
ERF_AdvectionSrcForMom_T.H
Go to the documentation of this file.
1 #include <ERF_IndexDefines.H>
2 #include <ERF_TerrainMetrics.H>
3 #include <ERF_Interpolation.H>
4 
5 /**
6  * Function for computing the advective tendency for the x-component of momentum
7  * with metric terms and for spatial order > 2
8  *
9  * @param[in] i,j,k indices of x-face at which to create tendency
10  * @param[in] rho_u x-component of momentum
11  * @param[in] rho_v y-component of momentum
12  * @param[in] Omega component of the momentum normal to the z-coordinate surface
13  * @param[in] z_nd height coordinate at nodes
14  * @param[in] ax Area fractions on x-faces
15  * @param[in] ay Area fractions on y-faces
16  * @param[in] az Area fractions on z-faces
17  * @param[in] detJ Jacobian of the metric transformation (= 1 if use_terrain is false)
18  * @param[in] cellSizeInv inverse of the mesh spacing
19  * @param[in] mf_u map factor on x-faces
20  * @param[in] mf_v map factor on y-faces
21  */
22 template<typename InterpType_H, typename InterpType_V>
23 AMREX_GPU_DEVICE
24 AMREX_FORCE_INLINE
25 amrex::Real
26 AdvectionSrcForXMom (int i, int j, int k,
27  const amrex::Array4<const amrex::Real>& rho_u,
28  const amrex::Array4<const amrex::Real>& rho_v,
29  const amrex::Array4<const amrex::Real>& Omega,
30  const amrex::Array4<const amrex::Real>& z_nd,
31  const amrex::Array4<const amrex::Real>& ax,
32  const amrex::Array4<const amrex::Real>& /*ay*/,
33  const amrex::Array4<const amrex::Real>& az,
34  const amrex::Array4<const amrex::Real>& detJ,
35  InterpType_H interp_u_h,
36  InterpType_V interp_u_v,
37  const amrex::Real upw_frac_h,
38  const amrex::Real upw_frac_v,
39  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
40  const amrex::Array4<const amrex::Real>& mf_u_inv,
41  const amrex::Array4<const amrex::Real>& mf_v_inv)
42 {
43  amrex::Real advectionSrc;
44  auto dxInv = cellSizeInv[0], dyInv = cellSizeInv[1], dzInv = cellSizeInv[2];
45 
46  amrex::Real rho_u_avg_lo, rho_u_avg_hi;
47  amrex::Real rho_v_avg_lo, rho_v_avg_hi;
48  amrex::Real Omega_avg_lo, Omega_avg_hi;
49 
50  amrex::Real interp_hi(0.), interp_lo(0.);
51 
52  // ****************************************************************************************
53  // X-fluxes (at cell centers)
54  // ****************************************************************************************
55 
56  rho_u_avg_hi = 0.5 * (rho_u(i+1, j, k) * mf_u_inv(i+1,j ,0) + rho_u(i, j, k) * mf_u_inv(i ,j ,0));
57  interp_u_h.InterpolateInX(i+1,j,k,0,interp_hi,rho_u_avg_hi,upw_frac_h);
58 
59  rho_u_avg_lo = 0.5 * (rho_u(i-1, j, k) * mf_u_inv(i-1,j ,0) + rho_u(i, j, k) * mf_u_inv(i ,j ,0));
60  interp_u_h.InterpolateInX(i,j,k,0,interp_lo,rho_u_avg_lo,upw_frac_h);
61 
62  amrex::Real centFluxXXNext = rho_u_avg_hi * interp_hi * 0.5 * (ax(i,j,k) + ax(i+1,j,k));
63  amrex::Real centFluxXXPrev = rho_u_avg_lo * interp_lo * 0.5 * (ax(i,j,k) + ax(i-1,j,k));
64 
65  // ****************************************************************************************
66  // Y-fluxes (at edges in k-direction)
67  // ****************************************************************************************
68  rho_v_avg_hi = 0.5 * (rho_v(i, j+1, k) * mf_v_inv(i ,j+1,0) + rho_v(i-1, j+1, k) * mf_v_inv(i-1,j+1,0));
69  interp_u_h.InterpolateInY(i,j+1,k,0,interp_hi,rho_v_avg_hi,upw_frac_h);
70 
71  rho_v_avg_lo = 0.5 * (rho_v(i, j , k) * mf_v_inv(i ,j ,0) + rho_v(i-1, j , k) * mf_v_inv(i-1,j ,0));
72  interp_u_h.InterpolateInY(i,j,k,0,interp_lo,rho_v_avg_lo,upw_frac_h);
73 
74  amrex::Real edgeFluxXYNext = rho_v_avg_hi * interp_hi * Compute_h_zeta_AtEdgeCenterK(i,j+1,k,cellSizeInv,z_nd);
75  amrex::Real edgeFluxXYPrev = rho_v_avg_lo * interp_lo * Compute_h_zeta_AtEdgeCenterK(i,j ,k,cellSizeInv,z_nd);
76 
77  // ****************************************************************************************
78  // Z-fluxes (at edges in j-direction)
79  // ****************************************************************************************
80  Omega_avg_hi = 0.5 * (Omega(i, j, k+1) + Omega(i-1, j, k+1)) * 0.5 * (az(i,j,k+1) + az(i-1,j,k+1));
81  Omega_avg_lo = 0.5 * (Omega(i, j, k ) + Omega(i-1, j, k )) * 0.5 * (az(i,j,k ) + az(i-1,j,k ));
82 
83  interp_u_v.InterpolateInZ(i,j,k+1,0,interp_hi,Omega_avg_hi,upw_frac_v);
84  interp_u_v.InterpolateInZ(i,j,k ,0,interp_lo,Omega_avg_lo,upw_frac_v);
85 
86  amrex::Real edgeFluxXZNext = Omega_avg_hi * interp_hi;
87  amrex::Real edgeFluxXZPrev = Omega_avg_lo * interp_lo;
88 
89  // ****************************************************************************************
90 
91  amrex::Real mfsq = 1 / (mf_u_inv(i,j,0) * mf_u_inv(i,j,0));
92 
93  advectionSrc = (centFluxXXNext - centFluxXXPrev) * dxInv * mfsq
94  + (edgeFluxXYNext - edgeFluxXYPrev) * dyInv * mfsq
95  + (edgeFluxXZNext - edgeFluxXZPrev) * dzInv;
96  advectionSrc /= 0.5*(detJ(i,j,k) + detJ(i-1,j,k));
97 
98  return advectionSrc;
99 }
100 
101 /**
102  * Function for computing the advective tendency for the y-component of momentum
103  * with metric terms and for spatial order > 2
104  *
105  * @param[in] i,j,k indices of y-face at which to create tendency
106  * @param[in] rho_u x-component of momentum
107  * @param[in] rho_v y-component of momentum
108  * @param[in] Omega component of the momentum normal to the z-coordinate surface
109  * @param[in] z_nd height coordinate at nodes
110  * @param[in] ax Area fractions on x-faces
111  * @param[in] ay Area fractions on y-faces
112  * @param[in] az Area fractions on z-faces
113  * @param[in] detJ Jacobian of the metric transformation (= 1 if use_terrain is false)
114  * @param[in] cellSizeInv inverse of the mesh spacing
115  * @param[in] mf_u map factor on x-faces
116  * @param[in] mf_v map factor on y-faces
117  */
118 template<typename InterpType_H, typename InterpType_V>
119 AMREX_GPU_DEVICE
120 AMREX_FORCE_INLINE
121 amrex::Real
122 AdvectionSrcForYMom (int i, int j, int k,
123  const amrex::Array4<const amrex::Real>& rho_u,
124  const amrex::Array4<const amrex::Real>& rho_v,
125  const amrex::Array4<const amrex::Real>& Omega,
126  const amrex::Array4<const amrex::Real>& z_nd,
127  const amrex::Array4<const amrex::Real>& /*ax*/,
128  const amrex::Array4<const amrex::Real>& ay,
129  const amrex::Array4<const amrex::Real>& az,
130  const amrex::Array4<const amrex::Real>& detJ,
131  InterpType_H interp_v_h,
132  InterpType_V interp_v_v,
133  const amrex::Real upw_frac_h,
134  const amrex::Real upw_frac_v,
135  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
136  const amrex::Array4<const amrex::Real>& mf_u_inv,
137  const amrex::Array4<const amrex::Real>& mf_v_inv)
138 {
139  amrex::Real advectionSrc;
140  auto dxInv = cellSizeInv[0], dyInv = cellSizeInv[1], dzInv = cellSizeInv[2];
141 
142  amrex::Real rho_u_avg_lo, rho_u_avg_hi;
143  amrex::Real rho_v_avg_lo, rho_v_avg_hi;
144  amrex::Real Omega_avg_lo, Omega_avg_hi;
145 
146  amrex::Real interp_hi(0.), interp_lo(0.);
147 
148  // ****************************************************************************************
149  // x-fluxes (at edges in k-direction)
150  // ****************************************************************************************
151  rho_u_avg_hi = 0.5 * (rho_u(i+1,j,k) * mf_u_inv(i+1,j,0) + rho_u(i+1,j-1,k) * mf_u_inv(i+1,j-1,0));
152  rho_u_avg_lo = 0.5 * (rho_u(i ,j,k) * mf_u_inv(i ,j,0) + rho_u(i ,j-1,k) * mf_u_inv(i ,j-1,0));
153 
154  interp_v_h.InterpolateInX(i+1,j,k,0,interp_hi,rho_u_avg_hi,upw_frac_h);
155  interp_v_h.InterpolateInX(i ,j,k,0,interp_lo,rho_u_avg_lo,upw_frac_h);
156 
157  amrex::Real edgeFluxYXNext = rho_u_avg_hi * interp_hi * Compute_h_zeta_AtEdgeCenterK(i+1,j,k,cellSizeInv,z_nd);
158  amrex::Real edgeFluxYXPrev = rho_u_avg_lo * interp_lo * Compute_h_zeta_AtEdgeCenterK(i ,j,k,cellSizeInv,z_nd);
159 
160  // ****************************************************************************************
161  // y-fluxes (at cell centers)
162  // ****************************************************************************************
163  rho_v_avg_hi = 0.5 * (rho_v(i,j,k) * mf_v_inv(i,j,0) + rho_v(i,j+1,k) * mf_v_inv(i,j+1,0));
164  rho_v_avg_lo = 0.5 * (rho_v(i,j,k) * mf_v_inv(i,j,0) + rho_v(i,j-1,k) * mf_v_inv(i,j-1,0));
165 
166  interp_v_h.InterpolateInY(i,j+1,k,0,interp_hi,rho_v_avg_hi,upw_frac_h);
167  interp_v_h.InterpolateInY(i,j ,k,0,interp_lo,rho_v_avg_lo,upw_frac_h);
168 
169  amrex::Real centFluxYYNext = rho_v_avg_hi * 0.5 * (ay(i,j,k) + ay(i,j+1,k)) * interp_hi;
170  amrex::Real centFluxYYPrev = rho_v_avg_lo * 0.5 * (ay(i,j,k) + ay(i,j-1,k)) * interp_lo;
171 
172  // ****************************************************************************************
173  // Z-fluxes (at edges in j-direction)
174  // ****************************************************************************************
175  Omega_avg_hi = 0.5 * (Omega(i, j, k+1) + Omega(i, j-1, k+1)) * 0.5 * (az(i,j,k+1) + az(i,j-1,k+1));
176  Omega_avg_lo = 0.5 * (Omega(i, j, k ) + Omega(i, j-1, k )) * 0.5 * (az(i,j,k ) + az(i,j-1,k ));
177 
178  interp_v_v.InterpolateInZ(i,j,k+1,0,interp_hi,Omega_avg_hi,upw_frac_v);
179  interp_v_v.InterpolateInZ(i,j,k ,0,interp_lo,Omega_avg_lo,upw_frac_v);
180 
181  amrex::Real edgeFluxYZNext = Omega_avg_hi * interp_hi;
182  amrex::Real edgeFluxYZPrev = Omega_avg_lo * interp_lo;
183 
184  // ****************************************************************************************
185 
186  amrex::Real mfsq = 1 / (mf_v_inv(i,j,0) * mf_v_inv(i,j,0));
187 
188  advectionSrc = (edgeFluxYXNext - edgeFluxYXPrev) * dxInv * mfsq
189  + (centFluxYYNext - centFluxYYPrev) * dyInv * mfsq
190  + (edgeFluxYZNext - edgeFluxYZPrev) * dzInv;
191  advectionSrc /= 0.5*(detJ(i,j,k) + detJ(i,j-1,k));
192 
193  return advectionSrc;
194 }
195 
196 /**
197  * Function for computing the advective tendency for the z-component of momentum
198  * with metric terms and for spatial order > 2
199  *
200  * @param[in] i,j,k indices of z-face at which to create tendency
201  * @param[in] rho_u x-component of momentum
202  * @param[in] rho_v y-component of momentum
203  * @param[in] Omega component of the momentum normal to the z-coordinate surface
204  * @param[in] w z-component of velocity
205  * @param[in] z_nd height coordinate at nodes
206  * @param[in] ax Area fractions on x-faces
207  * @param[in] ay Area fractions on y-faces
208  * @param[in] az Area fractions on z-faces
209  * @param[in] detJ Jacobian of the metric transformation (= 1 if use_terrain is false)
210  * @param[in] cellSizeInv inverse of the mesh spacing
211  * @param[in] mf_m map factor on cell centers
212  * @param[in] mf_u map factor on x-faces
213  * @param[in] mf_v map factor on y-faces
214  * @param[in] vert_adv_type int that defines advection stencil
215  * @param[in] lo_z_face minimum k value (z-face-centered)_in the domain at this level
216  * @param[in] hi_z_face maximum k value (z-face-centered) in the domain at this level
217  */
218 template<typename InterpType_H, typename InterpType_V, typename WallInterpType>
219 AMREX_GPU_DEVICE
220 AMREX_FORCE_INLINE
221 amrex::Real
222 AdvectionSrcForZMom (int i, int j, int k,
223  const amrex::Array4<const amrex::Real>& rho_u,
224  const amrex::Array4<const amrex::Real>& rho_v,
225  const amrex::Array4<const amrex::Real>& Omega,
226  const amrex::Array4<const amrex::Real>& w,
227  const amrex::Array4<const amrex::Real>& z_nd,
228  const amrex::Array4<const amrex::Real>& /*ax*/,
229  const amrex::Array4<const amrex::Real>& /*ay*/,
230  const amrex::Array4<const amrex::Real>& az,
231  const amrex::Array4<const amrex::Real>& detJ,
232  InterpType_H interp_omega_h,
233  InterpType_V interp_omega_v,
234  WallInterpType interp_omega_wall,
235  const amrex::Real upw_frac_h,
236  const amrex::Real upw_frac_v,
237  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
238  const amrex::Array4<const amrex::Real>& mf_m,
239  const amrex::Array4<const amrex::Real>& mf_u_inv,
240  const amrex::Array4<const amrex::Real>& mf_v_inv,
241  const AdvType vert_adv_type,
242  const int lo_z_face, const int hi_z_face)
243 {
244  amrex::Real advectionSrc;
245  auto dxInv = cellSizeInv[0], dyInv = cellSizeInv[1], dzInv = cellSizeInv[2];
246 
247  amrex::Real rho_u_avg_lo, rho_u_avg_hi;
248  amrex::Real rho_v_avg_lo, rho_v_avg_hi;
249  amrex::Real Omega_avg_lo, Omega_avg_hi;
250 
251  amrex::Real interp_hi(0.), interp_lo(0.);
252 
253  // ****************************************************************************************
254  // x-fluxes (at edges in j-direction)
255  // ****************************************************************************************
256  rho_u_avg_hi = 0.5 * (rho_u(i+1, j, k) + rho_u(i+1, j, k-1)) * mf_u_inv(i+1,j ,0);
257  rho_u_avg_lo = 0.5 * (rho_u(i , j, k) + rho_u(i , j, k-1)) * mf_u_inv(i ,j ,0);
258 
259  interp_omega_h.InterpolateInX(i+1,j,k,0,interp_hi,rho_u_avg_hi,upw_frac_h);
260  interp_omega_h.InterpolateInX(i,j,k,0,interp_lo,rho_u_avg_lo,upw_frac_h);
261 
262  amrex::Real edgeFluxZXNext = rho_u_avg_hi * interp_hi * Compute_h_zeta_AtEdgeCenterJ(i+1,j,k,cellSizeInv,z_nd);
263  amrex::Real edgeFluxZXPrev = rho_u_avg_lo * interp_lo * Compute_h_zeta_AtEdgeCenterJ(i ,j,k,cellSizeInv,z_nd);
264 
265  // ****************************************************************************************
266  // y-fluxes (at edges in i-direction)
267  // ****************************************************************************************
268  rho_v_avg_hi = 0.5 * (rho_v(i, j+1, k) + rho_v(i, j+1, k-1)) * mf_v_inv(i ,j+1,0);
269  interp_omega_h.InterpolateInY(i,j+1,k,0,interp_hi,rho_v_avg_hi,upw_frac_h);
270 
271  rho_v_avg_lo = 0.5 * (rho_v(i, j , k) + rho_v(i, j , k-1)) * mf_v_inv(i ,j ,0);
272  interp_omega_h.InterpolateInY(i,j,k,0,interp_lo,rho_v_avg_lo,upw_frac_h);
273 
274  amrex::Real edgeFluxZYNext = rho_v_avg_hi * interp_hi * Compute_h_zeta_AtEdgeCenterI(i,j+1,k,cellSizeInv,z_nd);
275  amrex::Real edgeFluxZYPrev = rho_v_avg_lo * interp_lo * Compute_h_zeta_AtEdgeCenterI(i,j ,k,cellSizeInv,z_nd);
276 
277  // ****************************************************************************************
278  // z-fluxes (at cell centers)
279  // ****************************************************************************************
280 
281  Omega_avg_hi = (k == hi_z_face) ? Omega(i,j,k) * az(i,j,k) :
282  0.25 * (Omega(i,j,k) + Omega(i,j,k+1)) * (az(i,j,k) + az(i,j,k+1));
283  amrex::Real centFluxZZNext = Omega_avg_hi;
284 
285  // int l_spatial_order_hi = std::min(std::min(vert_spatial_order, 2*(hi_z_face-k)), 2*(k+1));
286  // If k == hi_z_face-1, l_spatial_order_hi = 2
287  // If k == hi_z_face-2, l_spatial_order_hi = std::min(vert_spatial_order, 4);
288  // If k == lo_z_face+11 , l_spatial_order_hi = std::min(vert_spatial_order, 4);
289  if (k == hi_z_face) {
290  centFluxZZNext *= w(i,j,k);
291  } else {
292  if (k == hi_z_face-1)
293  {
294  interp_omega_wall.InterpolateInZ(i,j,k+1,0,interp_hi,Omega_avg_hi,0,AdvType::Centered_2nd);
295  } else if (k == hi_z_face-2 || k == lo_z_face+1) {
296  if (vert_adv_type != AdvType::Centered_2nd && vert_adv_type != AdvType::Upwind_3rd) {
297  interp_omega_wall.InterpolateInZ(i,j,k+1,0,interp_hi,Omega_avg_hi,0,AdvType::Centered_4th);
298  } else {
299  interp_omega_wall.InterpolateInZ(i,j,k+1,0,interp_hi,Omega_avg_hi,upw_frac_v,vert_adv_type);
300  }
301  } else {
302  interp_omega_v.InterpolateInZ(i,j,k+1,0,interp_hi,Omega_avg_hi,upw_frac_v);
303  }
304  centFluxZZNext *= interp_hi;
305  }
306 
307  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
308 
309  Omega_avg_lo = (k == 0) ? Omega(i,j,k) * az(i,j,k) :
310  0.25 * (Omega(i,j,k) + Omega(i,j,k-1)) * (az(i,j,k) + az(i,j,k-1));
311  amrex::Real centFluxZZPrev = Omega_avg_lo;
312 
313  // int l_spatial_order_lo = std::min(std::min(vert_spatial_order, 2*(hi_z_face+1-k)), 2*k);
314  // If k == hi_z_face-1, l_spatial_order_hi = 2
315  // If k == hi_z_face-2, l_spatial_order_hi = std::min(vert_spatial_order, 4);
316  // If k == lo_z_face+1, l_spatial_order_hi = std::min(vert_spatial_order, 4);
317  if (k == 0) {
318  centFluxZZPrev *= w(i,j,k);
319  } else {
320  if (k == lo_z_face+1) {
321  interp_omega_wall.InterpolateInZ(i,j,k,0,interp_lo,Omega_avg_lo,0,AdvType::Centered_2nd);
322  } else if (k == lo_z_face+2 || k == hi_z_face-1) {
323  if (vert_adv_type != AdvType::Centered_2nd && vert_adv_type != AdvType::Upwind_3rd) {
324  interp_omega_wall.InterpolateInZ(i,j,k,0,interp_lo,Omega_avg_lo,0,AdvType::Centered_4th);
325  } else {
326  interp_omega_wall.InterpolateInZ(i,j,k,0,interp_lo,Omega_avg_lo,upw_frac_v,vert_adv_type);
327  }
328  } else {
329  interp_omega_v.InterpolateInZ(i,j,k,0,interp_lo,Omega_avg_lo,upw_frac_v);
330  }
331  centFluxZZPrev *= interp_lo;
332  }
333 
334  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
335 
336  amrex::Real mfsq = mf_m(i,j,0) * mf_m(i,j,0);
337 
338  advectionSrc = (edgeFluxZXNext - edgeFluxZXPrev) * dxInv * mfsq
339  + (edgeFluxZYNext - edgeFluxZYPrev) * dyInv * mfsq
340  + (centFluxZZNext - centFluxZZPrev) * dzInv;
341 
342  amrex::Real denom = 0.5*(detJ(i,j,k) + detJ(i,j,k-1));
343  advectionSrc /= denom;
344 
345  return advectionSrc;
346 }
347 
348 /**
349  * Wrapper function for computing the advective tendency w/ spatial order > 2.
350  */
351 template<typename InterpType_H, typename InterpType_V, typename WallInterpType>
352 void
353 AdvectionSrcForMomWrapper (const amrex::Box& bxx, const amrex::Box& bxy, const amrex::Box& bxz,
354  const amrex::Array4<amrex::Real>& rho_u_rhs,
355  const amrex::Array4<amrex::Real>& rho_v_rhs,
356  const amrex::Array4<amrex::Real>& rho_w_rhs,
357  const amrex::Array4<const amrex::Real>& rho_u,
358  const amrex::Array4<const amrex::Real>& rho_v,
359  const amrex::Array4<const amrex::Real>& Omega,
360  const amrex::Array4<const amrex::Real>& u,
361  const amrex::Array4<const amrex::Real>& v,
362  const amrex::Array4<const amrex::Real>& w,
363  const amrex::Array4<const amrex::Real>& z_nd,
364  const amrex::Array4<const amrex::Real>& ax,
365  const amrex::Array4<const amrex::Real>& ay,
366  const amrex::Array4<const amrex::Real>& az,
367  const amrex::Array4<const amrex::Real>& detJ,
368  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
369  const amrex::Array4<const amrex::Real>& mf_m,
370  const amrex::Array4<const amrex::Real>& mf_u_inv,
371  const amrex::Array4<const amrex::Real>& mf_v_inv,
372  const amrex::Real upw_frac_h,
373  const amrex::Real upw_frac_v,
374  const AdvType vert_adv_type,
375  const int lo_z_face, const int hi_z_face)
376 {
377  // Instantiate the appropriate structs
378  InterpType_H interp_u_h(u); InterpType_V interp_u_v(u); // X-MOM
379  InterpType_H interp_v_h(v); InterpType_V interp_v_v(v); // Y-MOM
380  InterpType_H interp_w_h(w); InterpType_V interp_w_v(w); // Z-MOM
381  WallInterpType interp_w_wall(w); // Z-MOM @ wall
382 
383  amrex::ParallelFor(bxx,
384  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
385  {
386  rho_u_rhs(i, j, k) = -AdvectionSrcForXMom(i, j, k, rho_u, rho_v, Omega, z_nd, ax, ay, az, detJ,
387  interp_u_h, interp_u_v,
388  upw_frac_h, upw_frac_v,
389  cellSizeInv, mf_u_inv, mf_v_inv);
390  });
391 
392  amrex::ParallelFor(bxy,
393  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
394  {
395  rho_v_rhs(i, j, k) = -AdvectionSrcForYMom(i, j, k, rho_u, rho_v, Omega, z_nd, ax, ay, az, detJ,
396  interp_v_h, interp_v_v,
397  upw_frac_h, upw_frac_v,
398  cellSizeInv, mf_u_inv, mf_v_inv);
399  });
400 
401  amrex::ParallelFor(bxz,
402  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
403  {
404  rho_w_rhs(i, j, k) = -AdvectionSrcForZMom(i, j, k, rho_u, rho_v, Omega, w, z_nd, ax, ay, az, detJ,
405  interp_w_h, interp_w_v, interp_w_wall,
406  upw_frac_h, upw_frac_v,
407  cellSizeInv, mf_m, mf_u_inv, mf_v_inv,
408  vert_adv_type, lo_z_face, hi_z_face);
409  });
410 }
411 
412 /**
413  * Wrapper function for computing the advective tendency w/ spatial order > 2.
414  */
415 template<typename InterpType_H>
416 void
417 AdvectionSrcForMomVert (const amrex::Box& bxx, const amrex::Box& bxy, const amrex::Box& bxz,
418  const amrex::Array4<amrex::Real>& rho_u_rhs,
419  const amrex::Array4<amrex::Real>& rho_v_rhs,
420  const amrex::Array4<amrex::Real>& rho_w_rhs,
421  const amrex::Array4<const amrex::Real>& rho_u,
422  const amrex::Array4<const amrex::Real>& rho_v,
423  const amrex::Array4<const amrex::Real>& Omega,
424  const amrex::Array4<const amrex::Real>& u,
425  const amrex::Array4<const amrex::Real>& v,
426  const amrex::Array4<const amrex::Real>& w,
427  const amrex::Array4<const amrex::Real>& z_nd,
428  const amrex::Array4<const amrex::Real>& ax,
429  const amrex::Array4<const amrex::Real>& ay,
430  const amrex::Array4<const amrex::Real>& az,
431  const amrex::Array4<const amrex::Real>& detJ,
432  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
433  const amrex::Array4<const amrex::Real>& mf_m,
434  const amrex::Array4<const amrex::Real>& mf_u_inv,
435  const amrex::Array4<const amrex::Real>& mf_v_inv,
436  const amrex::Real upw_frac_h,
437  const amrex::Real upw_frac_v,
438  const AdvType vert_adv_type,
439  const int lo_z_face, const int hi_z_face)
440 {
441  if (vert_adv_type == AdvType::Centered_2nd) {
442  AdvectionSrcForMomWrapper<InterpType_H,CENTERED2,UPWINDALL>(bxx, bxy, bxz,
443  rho_u_rhs, rho_v_rhs, rho_w_rhs,
444  rho_u, rho_v, Omega, u, v, w, z_nd, ax, ay, az, detJ,
445  cellSizeInv, mf_m,
446  mf_u_inv, mf_v_inv,
447  upw_frac_h, upw_frac_v,
448  vert_adv_type, lo_z_face, hi_z_face);
449  } else if (vert_adv_type == AdvType::Upwind_3rd) {
450  AdvectionSrcForMomWrapper<InterpType_H,UPWIND3,UPWINDALL>(bxx, bxy, bxz,
451  rho_u_rhs, rho_v_rhs, rho_w_rhs,
452  rho_u, rho_v, Omega, u, v, w, z_nd, ax, ay, az, detJ,
453  cellSizeInv, mf_m,
454  mf_u_inv, mf_v_inv,
455  upw_frac_h, upw_frac_v,
456  vert_adv_type, lo_z_face, hi_z_face);
457  } else if (vert_adv_type == AdvType::Centered_4th) {
458  AdvectionSrcForMomWrapper<InterpType_H,CENTERED4,UPWINDALL>(bxx, bxy, bxz,
459  rho_u_rhs, rho_v_rhs, rho_w_rhs,
460  rho_u, rho_v, Omega, u, v, w, z_nd, ax, ay, az, detJ,
461  cellSizeInv, mf_m,
462  mf_u_inv, mf_v_inv,
463  upw_frac_h, upw_frac_v,
464  vert_adv_type, lo_z_face, hi_z_face);
465  } else if (vert_adv_type == AdvType::Upwind_5th) {
466  AdvectionSrcForMomWrapper<InterpType_H,UPWIND5,UPWINDALL>(bxx, bxy, bxz,
467  rho_u_rhs, rho_v_rhs, rho_w_rhs,
468  rho_u, rho_v, Omega, u, v, w, z_nd, ax, ay, az, detJ,
469  cellSizeInv, mf_m,
470  mf_u_inv, mf_v_inv,
471  upw_frac_h, upw_frac_v,
472  vert_adv_type, lo_z_face, hi_z_face);
473  } else if (vert_adv_type == AdvType::Centered_6th) {
474  AdvectionSrcForMomWrapper<InterpType_H,CENTERED6,UPWINDALL>(bxx, bxy, bxz,
475  rho_u_rhs, rho_v_rhs, rho_w_rhs,
476  rho_u, rho_v, Omega, u, v, w, z_nd, ax, ay, az, detJ,
477  cellSizeInv, mf_m,
478  mf_u_inv, mf_v_inv,
479  upw_frac_h, upw_frac_v,
480  vert_adv_type, lo_z_face, hi_z_face);
481  } else {
482  AMREX_ASSERT_WITH_MESSAGE(false, "Unknown advection scheme!");
483  }
484 }
AdvectionSrcForYMom
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real AdvectionSrcForYMom(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 > &Omega, const amrex::Array4< const amrex::Real > &z_nd, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &ay, const amrex::Array4< const amrex::Real > &az, const amrex::Array4< const amrex::Real > &detJ, InterpType_H interp_v_h, InterpType_V interp_v_v, const amrex::Real upw_frac_h, const amrex::Real upw_frac_v, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &mf_u_inv, const amrex::Array4< const amrex::Real > &mf_v_inv)
Definition: ERF_AdvectionSrcForMom_T.H:122
AdvectionSrcForZMom
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real AdvectionSrcForZMom(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 > &Omega, const amrex::Array4< const amrex::Real > &w, const amrex::Array4< const amrex::Real > &z_nd, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &az, const amrex::Array4< const amrex::Real > &detJ, InterpType_H interp_omega_h, InterpType_V interp_omega_v, WallInterpType interp_omega_wall, const amrex::Real upw_frac_h, const amrex::Real upw_frac_v, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &mf_m, const amrex::Array4< const amrex::Real > &mf_u_inv, const amrex::Array4< const amrex::Real > &mf_v_inv, const AdvType vert_adv_type, const int lo_z_face, const int hi_z_face)
Definition: ERF_AdvectionSrcForMom_T.H:222
AdvectionSrcForMomWrapper
void AdvectionSrcForMomWrapper(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 > &Omega, const amrex::Array4< const amrex::Real > &u, const amrex::Array4< const amrex::Real > &v, const amrex::Array4< const amrex::Real > &w, const amrex::Array4< const amrex::Real > &z_nd, const amrex::Array4< const amrex::Real > &ax, const amrex::Array4< const amrex::Real > &ay, const amrex::Array4< const amrex::Real > &az, const amrex::Array4< const amrex::Real > &detJ, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &mf_m, const amrex::Array4< const amrex::Real > &mf_u_inv, const amrex::Array4< const amrex::Real > &mf_v_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_T.H:353
AdvectionSrcForXMom
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real AdvectionSrcForXMom(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 > &Omega, const amrex::Array4< const amrex::Real > &z_nd, const amrex::Array4< const amrex::Real > &ax, const amrex::Array4< const amrex::Real > &, const amrex::Array4< const amrex::Real > &az, const amrex::Array4< const amrex::Real > &detJ, InterpType_H interp_u_h, InterpType_V interp_u_v, const amrex::Real upw_frac_h, const amrex::Real upw_frac_v, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &mf_u_inv, const amrex::Array4< const amrex::Real > &mf_v_inv)
Definition: ERF_AdvectionSrcForMom_T.H:26
AdvectionSrcForMomVert
void AdvectionSrcForMomVert(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 > &Omega, const amrex::Array4< const amrex::Real > &u, const amrex::Array4< const amrex::Real > &v, const amrex::Array4< const amrex::Real > &w, const amrex::Array4< const amrex::Real > &z_nd, const amrex::Array4< const amrex::Real > &ax, const amrex::Array4< const amrex::Real > &ay, const amrex::Array4< const amrex::Real > &az, const amrex::Array4< const amrex::Real > &detJ, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &mf_m, const amrex::Array4< const amrex::Real > &mf_u_inv, const amrex::Array4< const amrex::Real > &mf_v_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_T.H:417
ERF_IndexDefines.H
AdvType
AdvType
Definition: ERF_IndexDefines.H:191
AdvType::Upwind_3rd
@ Upwind_3rd
AdvType::Centered_4th
@ Centered_4th
AdvType::Upwind_5th
@ Upwind_5th
AdvType::Centered_6th
@ Centered_6th
AdvType::Centered_2nd
@ Centered_2nd
ERF_Interpolation.H
ERF_TerrainMetrics.H
Compute_h_zeta_AtEdgeCenterJ
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_zeta_AtEdgeCenterJ(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:264
Compute_h_zeta_AtEdgeCenterK
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_zeta_AtEdgeCenterK(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:219
Compute_h_zeta_AtEdgeCenterI
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_zeta_AtEdgeCenterI(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:308