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::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
38  const amrex::Array4<const amrex::Real>& mf_ux_inv,
39  const amrex::Array4<const amrex::Real>& mf_uy_inv,
40  const amrex::Array4<const amrex::Real>& mf_vx_inv)
41 {
42  amrex::Real advectionSrc;
43  auto dxInv = cellSizeInv[0], dyInv = cellSizeInv[1], dzInv = cellSizeInv[2];
44 
45  amrex::Real rho_u_avg_lo, rho_u_avg_hi;
46  amrex::Real rho_v_avg_lo, rho_v_avg_hi;
47  amrex::Real Omega_avg_lo, Omega_avg_hi;
48 
49  amrex::Real interp_hi(0.), interp_lo(0.);
50 
51  // ****************************************************************************************
52  // X-fluxes (at cell centers)
53  // ****************************************************************************************
54  // average in the i direction, inverse mapfac 1/m_y = Δy/Δη
55  rho_u_avg_hi = 0.5 * (rho_u(i+1, j, k) * mf_uy_inv(i+1, j, 0) + rho_u(i, j, k) * mf_uy_inv(i, j, 0));
56  rho_u_avg_lo = 0.5 * (rho_u(i-1, j, k) * mf_uy_inv(i-1, j, 0) + rho_u(i, j, k) * mf_uy_inv(i, j, 0));
57 
58  interp_u_h.InterpolateInX(i+1,j,k,0,interp_hi,rho_u_avg_hi);
59  interp_u_h.InterpolateInX(i ,j,k,0,interp_lo,rho_u_avg_lo);
60 
61  // note: ax = (face height averaged over y) * dzInv = Δz/Δζ
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  // average in the i direction, inverse mapfac 1/m_x = Δx/Δξ
69  rho_v_avg_hi = 0.5 * (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));
70  rho_v_avg_lo = 0.5 * (rho_v(i, j , k) * mf_vx_inv(i, j , 0) + rho_v(i-1, j , k) * mf_vx_inv(i-1, j , 0));
71 
72  interp_u_h.InterpolateInY(i,j+1,k,0,interp_hi,rho_v_avg_hi);
73  interp_u_h.InterpolateInY(i,j ,k,0,interp_lo,rho_v_avg_lo);
74 
75  amrex::Real edgeFluxXYNext = rho_v_avg_hi * interp_hi * Compute_h_zeta_AtEdgeCenterK(i,j+1,k,cellSizeInv,z_nd);
76  amrex::Real edgeFluxXYPrev = rho_v_avg_lo * interp_lo * Compute_h_zeta_AtEdgeCenterK(i,j ,k,cellSizeInv,z_nd);
77 
78  // ****************************************************************************************
79  // Z-fluxes (at edges in j-direction)
80  // ****************************************************************************************
81  // average in the i direction
82  Omega_avg_hi = 0.5 * (Omega(i, j, k+1) + Omega(i-1, j, k+1));
83  Omega_avg_lo = 0.5 * (Omega(i, j, k ) + Omega(i-1, j, k ));
84 
85  interp_u_v.InterpolateInZ(i,j,k+1,0,interp_hi,Omega_avg_hi);
86  interp_u_v.InterpolateInZ(i,j,k ,0,interp_lo,Omega_avg_lo);
87 
88  amrex::Real edgeFluxXZNext = Omega_avg_hi * interp_hi * 0.5 * (az(i, j, k+1) + az(i-1, j, k+1));
89  amrex::Real edgeFluxXZPrev = Omega_avg_lo * interp_lo * 0.5 * (az(i, j, k ) + az(i-1, j, k ));
90 
91  // ****************************************************************************************
92 
93  amrex::Real mfsq = 1. / (mf_ux_inv(i,j,0) * mf_uy_inv(i,j,0)); // == m_x * m_y = (Δξ/Δx) * (Δη/Δy)
94 
95  advectionSrc = (centFluxXXNext - centFluxXXPrev) * dxInv * mfsq // ~ ρuu(Δy/Δη)(Δz/Δζ) * (1/Δξ) * (Δξ/Δx)*(Δη/Δy) = ρuu(Δz/Δζ) * (1/Δx)
96  + (edgeFluxXYNext - edgeFluxXYPrev) * dyInv * mfsq // ~ ρvu(Δx/Δξ)(Δz/Δζ) * (1/Δη) * (Δξ/Δx)*(Δη/Δy) = ρvu(Δz/Δζ) * (1/Δy)
97  + (edgeFluxXZNext - edgeFluxXZPrev) * dzInv; // ~ Ωu * (1/Δζ)
98  advectionSrc /= 0.5*(detJ(i,j,k) + detJ(i-1,j,k)); // ... then divide by (Δz/Δζ) to get back 1/Δx, 1/Δy, and 1/Δz for differencing
99 
100  return advectionSrc;
101 }
102 
103 /**
104  * Function for computing the advective tendency for the y-component of momentum
105  * with metric terms and for spatial order > 2
106  *
107  * @param[in] i,j,k indices of y-face at which to create tendency
108  * @param[in] rho_u x-component of momentum
109  * @param[in] rho_v y-component of momentum
110  * @param[in] Omega component of the momentum normal to the z-coordinate surface
111  * @param[in] z_nd height coordinate at nodes
112  * @param[in] ax Area fractions on x-faces
113  * @param[in] ay Area fractions on y-faces
114  * @param[in] az Area fractions on z-faces
115  * @param[in] detJ Jacobian of the metric transformation (= 1 if use_terrain is false)
116  * @param[in] cellSizeInv inverse of the mesh spacing
117  * @param[in] mf_u map factor on x-faces
118  * @param[in] mf_v map factor on y-faces
119  */
120 template<typename InterpType_H, typename InterpType_V>
121 AMREX_GPU_DEVICE
122 AMREX_FORCE_INLINE
123 amrex::Real
124 AdvectionSrcForYMom (int i, int j, int k,
125  const amrex::Array4<const amrex::Real>& rho_u,
126  const amrex::Array4<const amrex::Real>& rho_v,
127  const amrex::Array4<const amrex::Real>& Omega,
128  const amrex::Array4<const amrex::Real>& z_nd,
129  const amrex::Array4<const amrex::Real>& /*ax*/,
130  const amrex::Array4<const amrex::Real>& ay,
131  const amrex::Array4<const amrex::Real>& az,
132  const amrex::Array4<const amrex::Real>& detJ,
133  InterpType_H interp_v_h,
134  InterpType_V interp_v_v,
135  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
136  const amrex::Array4<const amrex::Real>& mf_uy_inv,
137  const amrex::Array4<const amrex::Real>& mf_vx_inv,
138  const amrex::Array4<const amrex::Real>& mf_vy_inv)
139 {
140  amrex::Real advectionSrc;
141  auto dxInv = cellSizeInv[0], dyInv = cellSizeInv[1], dzInv = cellSizeInv[2];
142 
143  amrex::Real rho_u_avg_lo, rho_u_avg_hi;
144  amrex::Real rho_v_avg_lo, rho_v_avg_hi;
145  amrex::Real Omega_avg_lo, Omega_avg_hi;
146 
147  amrex::Real interp_hi(0.), interp_lo(0.);
148 
149  // ****************************************************************************************
150  // X-fluxes (at edges in k-direction)
151  // ****************************************************************************************
152  // average in the j direction, inverse mapfac 1/m_y = Δy/Δη
153  rho_u_avg_hi = 0.5 * (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));
154  rho_u_avg_lo = 0.5 * (rho_u(i , j, k) * mf_uy_inv(i , j, 0) + rho_u(i , j-1, k) * mf_uy_inv(i , j-1, 0));
155 
156  interp_v_h.InterpolateInX(i+1,j,k,0,interp_hi,rho_u_avg_hi);
157  interp_v_h.InterpolateInX(i ,j,k,0,interp_lo,rho_u_avg_lo);
158 
159  amrex::Real edgeFluxYXNext = rho_u_avg_hi * interp_hi * Compute_h_zeta_AtEdgeCenterK(i+1,j,k,cellSizeInv,z_nd);
160  amrex::Real edgeFluxYXPrev = rho_u_avg_lo * interp_lo * Compute_h_zeta_AtEdgeCenterK(i ,j,k,cellSizeInv,z_nd);
161 
162  // ****************************************************************************************
163  // Y-fluxes (at cell centers)
164  // ****************************************************************************************
165  // average in the j direction, inverse mapfac 1/m_x = Δx/Δξ
166  rho_v_avg_hi = 0.5 * (rho_v(i, j+1, k) * mf_vx_inv(i, j+1, 0) + rho_v(i, j, k) * mf_vx_inv(i, j, 0));
167  rho_v_avg_lo = 0.5 * (rho_v(i, j-1, k) * mf_vx_inv(i, j-1, 0) + rho_v(i, j, k) * mf_vx_inv(i, j, 0));
168 
169  interp_v_h.InterpolateInY(i,j+1,k,0,interp_hi,rho_v_avg_hi);
170  interp_v_h.InterpolateInY(i,j ,k,0,interp_lo,rho_v_avg_lo);
171 
172  // note: ay = (face height averaged over x) * dzInv = Δz/Δζ
173  amrex::Real centFluxYYNext = rho_v_avg_hi * interp_hi * 0.5 * (ay(i,j,k) + ay(i,j+1,k));
174  amrex::Real centFluxYYPrev = rho_v_avg_lo * interp_lo * 0.5 * (ay(i,j,k) + ay(i,j-1,k));
175 
176  // ****************************************************************************************
177  // Z-fluxes (at edges in i-direction)
178  // ****************************************************************************************
179  // average in the j direction
180  Omega_avg_hi = 0.5 * (Omega(i, j, k+1) + Omega(i, j-1, k+1));
181  Omega_avg_lo = 0.5 * (Omega(i, j, k ) + Omega(i, j-1, k ));
182 
183  interp_v_v.InterpolateInZ(i,j,k+1,0,interp_hi,Omega_avg_hi);
184  interp_v_v.InterpolateInZ(i,j,k ,0,interp_lo,Omega_avg_lo);
185 
186  amrex::Real edgeFluxYZNext = Omega_avg_hi * interp_hi * 0.5 * (az(i, j, k+1) + az(i, j-1, k+1));
187  amrex::Real edgeFluxYZPrev = Omega_avg_lo * interp_lo * 0.5 * (az(i, j, k ) + az(i, j-1, k ));
188 
189  // ****************************************************************************************
190 
191  amrex::Real mfsq = 1. / (mf_vx_inv(i,j,0) * mf_vy_inv(i,j,0)); // == m_x * m_y = (Δξ/Δx) * (Δη/Δy)
192 
193  advectionSrc = (edgeFluxYXNext - edgeFluxYXPrev) * dxInv * mfsq // ~ ρuv(Δy/Δη)(Δz/Δζ) * (1/Δξ) * (Δξ/Δx)*(Δη/Δy) = ρuu(Δz/Δζ) * (1/Δx)
194  + (centFluxYYNext - centFluxYYPrev) * dyInv * mfsq // ~ ρvv(Δx/Δξ)(Δz/Δζ) * (1/Δη) * (Δξ/Δx)*(Δη/Δy) = ρvu(Δz/Δζ) * (1/Δy)
195  + (edgeFluxYZNext - edgeFluxYZPrev) * dzInv; // ~ Ωv * (1/Δζ)
196  advectionSrc /= 0.5*(detJ(i,j,k) + detJ(i,j-1,k)); // ... then divide by (Δz/Δζ) to get back 1/Δx, 1/Δy, and 1/Δz for differencing
197 
198  return advectionSrc;
199 }
200 
201 /**
202  * Function for computing the advective tendency for the z-component of momentum
203  * with metric terms and for spatial order > 2
204  *
205  * @param[in] i,j,k indices of z-face at which to create tendency
206  * @param[in] rho_u x-component of momentum
207  * @param[in] rho_v y-component of momentum
208  * @param[in] Omega component of the momentum normal to the z-coordinate surface
209  * @param[in] w z-component of velocity
210  * @param[in] z_nd height coordinate at nodes
211  * @param[in] ax Area fractions on x-faces
212  * @param[in] ay Area fractions on y-faces
213  * @param[in] az Area fractions on z-faces
214  * @param[in] detJ Jacobian of the metric transformation (= 1 if use_terrain is false)
215  * @param[in] cellSizeInv inverse of the mesh spacing
216  * @param[in] mf_m map factor on cell centers
217  * @param[in] mf_u map factor on x-faces
218  * @param[in] mf_v map factor on y-faces
219  * @param[in] vert_adv_type int that defines advection stencil
220  * @param[in] lo_z_face minimum k value (z-face-centered)_in the domain at this level
221  * @param[in] hi_z_face maximum k value (z-face-centered) in the domain at this level
222  */
223 template<typename InterpType_H, typename InterpType_V, typename WallInterpType>
224 AMREX_GPU_DEVICE
225 AMREX_FORCE_INLINE
226 amrex::Real
227 AdvectionSrcForZMom (int i, int j, int k,
228  const amrex::Array4<const amrex::Real>& rho_u,
229  const amrex::Array4<const amrex::Real>& rho_v,
230  const amrex::Array4<const amrex::Real>& Omega,
231  const amrex::Array4<const amrex::Real>& w,
232  const amrex::Array4<const amrex::Real>& z_nd,
233  const amrex::Array4<const amrex::Real>& /*ax*/,
234  const amrex::Array4<const amrex::Real>& /*ay*/,
235  const amrex::Array4<const amrex::Real>& az,
236  const amrex::Array4<const amrex::Real>& detJ,
237  InterpType_H interp_omega_h,
238  InterpType_V interp_omega_v,
239  WallInterpType interp_omega_wall,
240  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
241  const amrex::Array4<const amrex::Real>& mf_mx,
242  const amrex::Array4<const amrex::Real>& mf_my,
243  const amrex::Array4<const amrex::Real>& mf_uy_inv,
244  const amrex::Array4<const amrex::Real>& mf_vx_inv,
245  const AdvType vert_adv_type,
246  const int lo_z_face, const int hi_z_face)
247 {
248  amrex::Real advectionSrc;
249  auto dxInv = cellSizeInv[0], dyInv = cellSizeInv[1], dzInv = cellSizeInv[2];
250 
251  amrex::Real rho_u_avg_lo, rho_u_avg_hi;
252  amrex::Real rho_v_avg_lo, rho_v_avg_hi;
253  amrex::Real Omega_avg_lo, Omega_avg_hi;
254 
255  amrex::Real interp_hi(0.), interp_lo(0.);
256 
257  // ****************************************************************************************
258  // x-fluxes (at edges in j-direction)
259  // ****************************************************************************************
260  // average in the k direction, inverse mapfac 1/m_y = Δy/Δη
261  rho_u_avg_hi = 0.5 * (rho_u(i+1, j, k) + rho_u(i+1, j, k-1)) * mf_uy_inv(i+1,j ,0);
262  rho_u_avg_lo = 0.5 * (rho_u(i , j, k) + rho_u(i , j, k-1)) * mf_uy_inv(i ,j ,0);
263 
264  interp_omega_h.InterpolateInX(i+1,j,k,0,interp_hi,rho_u_avg_hi);
265  interp_omega_h.InterpolateInX(i ,j,k,0,interp_lo,rho_u_avg_lo);
266 
267  amrex::Real edgeFluxZXNext = rho_u_avg_hi * interp_hi * Compute_h_zeta_AtEdgeCenterJ(i+1,j,k,cellSizeInv,z_nd);
268  amrex::Real edgeFluxZXPrev = rho_u_avg_lo * interp_lo * Compute_h_zeta_AtEdgeCenterJ(i ,j,k,cellSizeInv,z_nd);
269 
270  // ****************************************************************************************
271  // y-fluxes (at edges in i-direction)
272  // ****************************************************************************************
273  // average in the k direction, inverse mapfac 1/m_x = Δx/Δξ
274  rho_v_avg_hi = 0.5 * (rho_v(i, j+1, k) + rho_v(i, j+1, k-1)) * mf_vx_inv(i ,j+1,0);
275  rho_v_avg_lo = 0.5 * (rho_v(i, j , k) + rho_v(i, j , k-1)) * mf_vx_inv(i ,j ,0);
276 
277  interp_omega_h.InterpolateInY(i,j+1,k,0,interp_hi,rho_v_avg_hi);
278  interp_omega_h.InterpolateInY(i,j ,k,0,interp_lo,rho_v_avg_lo);
279 
280  amrex::Real edgeFluxZYNext = rho_v_avg_hi * interp_hi * Compute_h_zeta_AtEdgeCenterI(i,j+1,k,cellSizeInv,z_nd);
281  amrex::Real edgeFluxZYPrev = rho_v_avg_lo * interp_lo * Compute_h_zeta_AtEdgeCenterI(i,j ,k,cellSizeInv,z_nd);
282 
283  // ****************************************************************************************
284  // z-fluxes (at cell centers)
285  // ****************************************************************************************
286 
287  Omega_avg_hi = (k == hi_z_face) ? Omega(i,j,k) * az(i,j,k) :
288  0.25 * (Omega(i,j,k) + Omega(i,j,k+1)) * (az(i,j,k) + az(i,j,k+1));
289  amrex::Real centFluxZZNext = Omega_avg_hi;
290 
291  // int l_spatial_order_hi = std::min(std::min(vert_spatial_order, 2*(hi_z_face-k)), 2*(k+1));
292  // If k == hi_z_face-1, l_spatial_order_hi = 2
293  // If k == hi_z_face-2, l_spatial_order_hi = std::min(vert_spatial_order, 4);
294  // If k == lo_z_face+11 , l_spatial_order_hi = std::min(vert_spatial_order, 4);
295  if (k == hi_z_face) {
296  centFluxZZNext *= w(i,j,k);
297  } else {
298  if (k == hi_z_face-1)
299  {
300  interp_omega_wall.InterpolateInZ(i,j,k+1,0,interp_hi,Omega_avg_hi,AdvType::Centered_2nd);
301  } else if (k == hi_z_face-2 || k == lo_z_face+1) {
302  if (vert_adv_type != AdvType::Centered_2nd && vert_adv_type != AdvType::Upwind_3rd) {
303  interp_omega_wall.InterpolateInZ(i,j,k+1,0,interp_hi,Omega_avg_hi,AdvType::Centered_4th);
304  } else {
305  interp_omega_wall.InterpolateInZ(i,j,k+1,0,interp_hi,Omega_avg_hi,vert_adv_type);
306  }
307  } else {
308  interp_omega_v.InterpolateInZ(i,j,k+1,0,interp_hi,Omega_avg_hi);
309  }
310  centFluxZZNext *= interp_hi;
311  }
312 
313  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
314 
315  Omega_avg_lo = (k == 0) ? Omega(i,j,k) * az(i,j,k) :
316  0.25 * (Omega(i,j,k) + Omega(i,j,k-1)) * (az(i,j,k) + az(i,j,k-1));
317  amrex::Real centFluxZZPrev = Omega_avg_lo;
318 
319  // int l_spatial_order_lo = std::min(std::min(vert_spatial_order, 2*(hi_z_face+1-k)), 2*k);
320  // If k == hi_z_face-1, l_spatial_order_hi = 2
321  // If k == hi_z_face-2, l_spatial_order_hi = std::min(vert_spatial_order, 4);
322  // If k == lo_z_face+1, l_spatial_order_hi = std::min(vert_spatial_order, 4);
323  if (k == 0) {
324  centFluxZZPrev *= w(i,j,k);
325  } else {
326  if (k == lo_z_face+1) {
327  interp_omega_wall.InterpolateInZ(i,j,k,0,interp_lo,Omega_avg_lo,AdvType::Centered_2nd);
328  } else if (k == lo_z_face+2 || k == hi_z_face-1) {
329  if (vert_adv_type != AdvType::Centered_2nd && vert_adv_type != AdvType::Upwind_3rd) {
330  interp_omega_wall.InterpolateInZ(i,j,k,0,interp_lo,Omega_avg_lo,AdvType::Centered_4th);
331  } else {
332  interp_omega_wall.InterpolateInZ(i,j,k,0,interp_lo,Omega_avg_lo,vert_adv_type);
333  }
334  } else {
335  interp_omega_v.InterpolateInZ(i,j,k,0,interp_lo,Omega_avg_lo);
336  }
337  centFluxZZPrev *= interp_lo;
338  }
339 
340  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
341 
342  amrex::Real mfsq = mf_mx(i,j,0) * mf_my(i,j,0);
343 
344  advectionSrc = (edgeFluxZXNext - edgeFluxZXPrev) * dxInv * mfsq
345  + (edgeFluxZYNext - edgeFluxZYPrev) * dyInv * mfsq
346  + (centFluxZZNext - centFluxZZPrev) * dzInv;
347 
348  amrex::Real denom = 0.5*(detJ(i,j,k) + detJ(i,j,k-1));
349  advectionSrc /= denom;
350 
351  return advectionSrc;
352 }
353 
354 /**
355  * Wrapper function for computing the advective tendency w/ spatial order > 2.
356  */
357 template<typename InterpType_H, typename InterpType_V, typename WallInterpType>
358 void
359 AdvectionSrcForMomWrapper (const amrex::Box& bxx, const amrex::Box& bxy, const amrex::Box& bxz,
360  const amrex::Array4<amrex::Real>& rho_u_rhs,
361  const amrex::Array4<amrex::Real>& rho_v_rhs,
362  const amrex::Array4<amrex::Real>& rho_w_rhs,
363  const amrex::Array4<const amrex::Real>& rho_u,
364  const amrex::Array4<const amrex::Real>& rho_v,
365  const amrex::Array4<const amrex::Real>& Omega,
366  const amrex::Array4<const amrex::Real>& u,
367  const amrex::Array4<const amrex::Real>& v,
368  const amrex::Array4<const amrex::Real>& w,
369  const amrex::Array4<const amrex::Real>& z_nd,
370  const amrex::Array4<const amrex::Real>& ax,
371  const amrex::Array4<const amrex::Real>& ay,
372  const amrex::Array4<const amrex::Real>& az,
373  const amrex::Array4<const amrex::Real>& detJ,
374  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
375  const amrex::Array4<const amrex::Real>& mf_mx,
376  const amrex::Array4<const amrex::Real>& mf_ux_inv,
377  const amrex::Array4<const amrex::Real>& mf_vx_inv,
378  const amrex::Array4<const amrex::Real>& mf_my,
379  const amrex::Array4<const amrex::Real>& mf_uy_inv,
380  const amrex::Array4<const amrex::Real>& mf_vy_inv,
381  const amrex::Real upw_frac_h,
382  const amrex::Real upw_frac_v,
383  const AdvType vert_adv_type,
384  const int lo_z_face, const int hi_z_face)
385 {
386  // Instantiate the appropriate structs
387  InterpType_H interp_u_h(u, upw_frac_h); InterpType_V interp_u_v(u, upw_frac_v); // X-MOM
388  InterpType_H interp_v_h(v, upw_frac_h); InterpType_V interp_v_v(v, upw_frac_v); // Y-MOM
389  InterpType_H interp_w_h(w, upw_frac_h); InterpType_V interp_w_v(w, upw_frac_v); // Z-MOM
390  WallInterpType interp_w_wall(w, upw_frac_v); // Z-MOM @ wall
391 
392  amrex::ParallelFor(bxx,
393  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
394  {
395  rho_u_rhs(i, j, k) = -AdvectionSrcForXMom(i, j, k, rho_u, rho_v, Omega, z_nd, ax, ay, az, detJ,
396  interp_u_h, interp_u_v,
397  cellSizeInv, mf_ux_inv, mf_uy_inv, mf_vx_inv);
398  });
399 
400  amrex::ParallelFor(bxy,
401  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
402  {
403  rho_v_rhs(i, j, k) = -AdvectionSrcForYMom(i, j, k, rho_u, rho_v, Omega, z_nd, ax, ay, az, detJ,
404  interp_v_h, interp_v_v,
405  cellSizeInv, mf_uy_inv, mf_vx_inv, mf_vy_inv);
406  });
407 
408  amrex::ParallelFor(bxz,
409  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
410  {
411  rho_w_rhs(i, j, k) = -AdvectionSrcForZMom(i, j, k, rho_u, rho_v, Omega, w, z_nd, ax, ay, az, detJ,
412  interp_w_h, interp_w_v, interp_w_wall,
413  cellSizeInv, mf_mx, mf_my, mf_uy_inv, mf_vx_inv,
414  vert_adv_type, lo_z_face, hi_z_face);
415  });
416 }
417 
418 /**
419  * Wrapper function for computing the advective tendency w/ spatial order > 2.
420  */
421 template<typename InterpType_H>
422 void
423 AdvectionSrcForMomVert (const amrex::Box& bxx, const amrex::Box& bxy, const amrex::Box& bxz,
424  const amrex::Array4<amrex::Real>& rho_u_rhs,
425  const amrex::Array4<amrex::Real>& rho_v_rhs,
426  const amrex::Array4<amrex::Real>& rho_w_rhs,
427  const amrex::Array4<const amrex::Real>& rho_u,
428  const amrex::Array4<const amrex::Real>& rho_v,
429  const amrex::Array4<const amrex::Real>& Omega,
430  const amrex::Array4<const amrex::Real>& u,
431  const amrex::Array4<const amrex::Real>& v,
432  const amrex::Array4<const amrex::Real>& w,
433  const amrex::Array4<const amrex::Real>& z_nd,
434  const amrex::Array4<const amrex::Real>& ax,
435  const amrex::Array4<const amrex::Real>& ay,
436  const amrex::Array4<const amrex::Real>& az,
437  const amrex::Array4<const amrex::Real>& detJ,
438  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
439  const amrex::Array4<const amrex::Real>& mf_mx,
440  const amrex::Array4<const amrex::Real>& mf_ux_inv,
441  const amrex::Array4<const amrex::Real>& mf_vx_inv,
442  const amrex::Array4<const amrex::Real>& mf_my,
443  const amrex::Array4<const amrex::Real>& mf_uy_inv,
444  const amrex::Array4<const amrex::Real>& mf_vy_inv,
445  const amrex::Real upw_frac_h,
446  const amrex::Real upw_frac_v,
447  const AdvType vert_adv_type,
448  const int lo_z_face, const int hi_z_face)
449 {
450  if (vert_adv_type == AdvType::Centered_2nd) {
451  AdvectionSrcForMomWrapper<InterpType_H,CENTERED2,UPWINDALL>(bxx, bxy, bxz,
452  rho_u_rhs, rho_v_rhs, rho_w_rhs,
453  rho_u, rho_v, Omega, u, v, w, z_nd, ax, ay, az, detJ,
454  cellSizeInv,
455  mf_mx, mf_ux_inv, mf_vx_inv,
456  mf_my, mf_uy_inv, mf_vy_inv,
457  upw_frac_h, upw_frac_v,
458  vert_adv_type, lo_z_face, hi_z_face);
459  } else if (vert_adv_type == AdvType::Upwind_3rd) {
460  AdvectionSrcForMomWrapper<InterpType_H,UPWIND3,UPWINDALL>(bxx, bxy, bxz,
461  rho_u_rhs, rho_v_rhs, rho_w_rhs,
462  rho_u, rho_v, Omega, u, v, w, z_nd, ax, ay, az, detJ,
463  cellSizeInv,
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, lo_z_face, hi_z_face);
468  } else if (vert_adv_type == AdvType::Centered_4th) {
469  AdvectionSrcForMomWrapper<InterpType_H,CENTERED4,UPWINDALL>(bxx, bxy, bxz,
470  rho_u_rhs, rho_v_rhs, rho_w_rhs,
471  rho_u, rho_v, Omega, u, v, w, z_nd, ax, ay, az, detJ,
472  cellSizeInv,
473  mf_mx, mf_ux_inv, mf_vx_inv,
474  mf_my, mf_uy_inv, mf_vy_inv,
475  upw_frac_h, upw_frac_v,
476  vert_adv_type, lo_z_face, hi_z_face);
477  } else if (vert_adv_type == AdvType::Upwind_5th) {
478  AdvectionSrcForMomWrapper<InterpType_H,UPWIND5,UPWINDALL>(bxx, bxy, bxz,
479  rho_u_rhs, rho_v_rhs, rho_w_rhs,
480  rho_u, rho_v, Omega, u, v, w, z_nd, ax, ay, az, detJ,
481  cellSizeInv,
482  mf_mx, mf_ux_inv, mf_vx_inv,
483  mf_my, mf_uy_inv, mf_vy_inv,
484  upw_frac_h, upw_frac_v,
485  vert_adv_type, lo_z_face, hi_z_face);
486  } else if (vert_adv_type == AdvType::Centered_6th) {
487  AdvectionSrcForMomWrapper<InterpType_H,CENTERED6,UPWINDALL>(bxx, bxy, bxz,
488  rho_u_rhs, rho_v_rhs, rho_w_rhs,
489  rho_u, rho_v, Omega, u, v, w, z_nd, ax, ay, az, detJ,
490  cellSizeInv,
491  mf_mx, mf_ux_inv, mf_vx_inv,
492  mf_my, mf_uy_inv, mf_vy_inv,
493  upw_frac_h, upw_frac_v,
494  vert_adv_type, lo_z_face, hi_z_face);
495 
496  } else if (vert_adv_type == AdvType::Weno_3) {
497  AdvectionSrcForMomWrapper<InterpType_H,WENO3,UPWINDALL>(bxx, bxy, bxz,
498  rho_u_rhs, rho_v_rhs, rho_w_rhs,
499  rho_u, rho_v, Omega, u, v, w, z_nd, ax, ay, az, detJ,
500  cellSizeInv,
501  mf_mx, mf_ux_inv, mf_vx_inv,
502  mf_my, mf_uy_inv, mf_vy_inv,
503  upw_frac_h, upw_frac_v,
504  vert_adv_type, lo_z_face, hi_z_face);
505  } else if (vert_adv_type == AdvType::Weno_3Z) {
506  AdvectionSrcForMomWrapper<InterpType_H,WENO_Z3,UPWINDALL>(bxx, bxy, bxz,
507  rho_u_rhs, rho_v_rhs, rho_w_rhs,
508  rho_u, rho_v, Omega, u, v, w, z_nd, ax, ay, az, detJ,
509  cellSizeInv,
510  mf_mx, mf_ux_inv, mf_vx_inv,
511  mf_my, mf_uy_inv, mf_vy_inv,
512  upw_frac_h, upw_frac_v,
513  vert_adv_type, lo_z_face, hi_z_face);
514  } else if (vert_adv_type == AdvType::Weno_3MZQ) {
515  AdvectionSrcForMomWrapper<InterpType_H,WENO_MZQ3,UPWINDALL>(bxx, bxy, bxz,
516  rho_u_rhs, rho_v_rhs, rho_w_rhs,
517  rho_u, rho_v, Omega, u, v, w, z_nd, ax, ay, az, detJ,
518  cellSizeInv,
519  mf_mx, mf_ux_inv, mf_vx_inv,
520  mf_my, mf_uy_inv, mf_vy_inv,
521  upw_frac_h, upw_frac_v,
522  vert_adv_type, lo_z_face, hi_z_face);
523  } else if (vert_adv_type == AdvType::Weno_5) {
524  AdvectionSrcForMomWrapper<InterpType_H,WENO5,UPWINDALL>(bxx, bxy, bxz,
525  rho_u_rhs, rho_v_rhs, rho_w_rhs,
526  rho_u, rho_v, Omega, u, v, w, z_nd, ax, ay, az, detJ,
527  cellSizeInv,
528  mf_mx, mf_ux_inv, mf_vx_inv,
529  mf_my, mf_uy_inv, mf_vy_inv,
530  upw_frac_h, upw_frac_v,
531  vert_adv_type, lo_z_face, hi_z_face);
532  } else if (vert_adv_type == AdvType::Weno_5Z) {
533  AdvectionSrcForMomWrapper<InterpType_H,WENO_Z5,UPWINDALL>(bxx, bxy, bxz,
534  rho_u_rhs, rho_v_rhs, rho_w_rhs,
535  rho_u, rho_v, Omega, u, v, w, z_nd, ax, ay, az, detJ,
536  cellSizeInv,
537  mf_mx, mf_ux_inv, mf_vx_inv,
538  mf_my, mf_uy_inv, mf_vy_inv,
539  upw_frac_h, upw_frac_v,
540  vert_adv_type, lo_z_face, hi_z_face);
541  } else if (vert_adv_type == AdvType::Weno_7) {
542  AdvectionSrcForMomWrapper<InterpType_H,WENO7,UPWINDALL>(bxx, bxy, bxz,
543  rho_u_rhs, rho_v_rhs, rho_w_rhs,
544  rho_u, rho_v, Omega, u, v, w, z_nd, ax, ay, az, detJ,
545  cellSizeInv,
546  mf_mx, mf_ux_inv, mf_vx_inv,
547  mf_my, mf_uy_inv, mf_vy_inv,
548  upw_frac_h, upw_frac_v,
549  vert_adv_type, lo_z_face, hi_z_face);
550  } else if (vert_adv_type == AdvType::Weno_7Z) {
551  AdvectionSrcForMomWrapper<InterpType_H,WENO_Z7,UPWINDALL>(bxx, bxy, bxz,
552  rho_u_rhs, rho_v_rhs, rho_w_rhs,
553  rho_u, rho_v, Omega, u, v, w, z_nd, ax, ay, az, detJ,
554  cellSizeInv,
555  mf_mx, mf_ux_inv, mf_vx_inv,
556  mf_my, mf_uy_inv, mf_vy_inv,
557  upw_frac_h, upw_frac_v,
558  vert_adv_type, lo_z_face, hi_z_face);
559  } else {
560  AMREX_ASSERT_WITH_MESSAGE(false, "Unknown advection scheme!");
561  }
562 }
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_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_T.H:359
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::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, 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)
Definition: ERF_AdvectionSrcForMom_T.H:227
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::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, 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)
Definition: ERF_AdvectionSrcForMom_T.H:124
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::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, 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)
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_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_T.H:423
ERF_IndexDefines.H
AdvType
AdvType
Definition: ERF_IndexDefines.H:221
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
ERF_Interpolation.H
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
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:273
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:228
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:317