ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
ERF_TerrainMetrics.H
Go to the documentation of this file.
1 #ifndef ERF_TERRAIN_METRIC_H_
2 #define ERF_TERRAIN_METRIC_H_
3 
4 #include <AMReX.H>
5 #include <AMReX_Geometry.H>
6 #include <AMReX_MultiFab.H>
7 #include <ERF_IndexDefines.H>
8 
9 /**
10  * Utility routines for constructing terrain metric terms
11  */
12 
13 // Declare functions for ERF.cpp
14 void init_zlevels (amrex::Vector<amrex::Vector<amrex::Real>> & zlevels_stag,
15  amrex::Vector<amrex::Vector<amrex::Real>> & stretched_dz_h,
16  amrex::Vector<amrex::Gpu::DeviceVector<amrex::Real>> & stretched_dz_d,
17  amrex::Vector<amrex::Geometry> const& geom,
18  amrex::Vector<amrex::IntVect> const& ref_ratio,
19  const amrex::Real grid_stretching_ratio,
20  const amrex::Real zsurf,
21  const amrex::Real dz0);
22 
23 void init_terrain_grid (int lev, const amrex::Geometry& geom,
24  amrex::MultiFab& z_phys_nd,
25  amrex::Vector<amrex::Real> const& z_levels_h,
26  amrex::GpuArray<ERF_BC, AMREX_SPACEDIM*2>& phys_bc_type);
27 
28 void init_which_terrain_grid (int lev, const amrex::Geometry& geom,
29  amrex::MultiFab& z_phys_nd,
30  amrex::Vector<amrex::Real> const& z_levels_h);
31 
32 //*****************************************************************************************
33 // Compute terrain metric terms at cell-center
34 //*****************************************************************************************
35 // Metric is at cell center
36 AMREX_FORCE_INLINE
37 AMREX_GPU_DEVICE
38 amrex::Real
39 Compute_h_zeta_AtCellCenter (const int &i, const int &j, const int &k,
40  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
41  const amrex::Array4<const amrex::Real>& z_nd)
42 {
43  amrex::Real dzInv = cellSizeInv[2];
44  amrex::Real met_h_zeta = 0.25 * dzInv *
45  ( z_nd(i+1,j,k+1) + z_nd(i+1,j+1,k+1) + z_nd(i,j,k+1) + z_nd(i,j+1,k+1)
46  -z_nd(i+1,j,k ) - z_nd(i+1,j+1,k ) - z_nd(i,j,k ) - z_nd(i,j+1,k ) );
47  return met_h_zeta;
48 }
49 
50 // Metric is at cell center
51 AMREX_GPU_DEVICE
52 AMREX_FORCE_INLINE
53 amrex::Real
54 Compute_h_xi_AtCellCenter (const int &i, const int &j, const int &k,
55  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
56  const amrex::Array4<const amrex::Real>& z_nd)
57 {
58  amrex::Real dxInv = cellSizeInv[0];
59  amrex::Real met_h_xi = 0.25 * dxInv *
60  ( z_nd(i+1,j,k) + z_nd(i+1,j+1,k) + z_nd(i+1,j,k+1) + z_nd(i+1,j+1,k+1)
61  -z_nd(i ,j,k) - z_nd(i ,j+1,k) - z_nd(i ,j,k+1) - z_nd(i ,j+1,k+1) );
62  return met_h_xi;
63 }
64 
65 // Metric is at cell center
66 AMREX_GPU_DEVICE
67 AMREX_FORCE_INLINE
68 amrex::Real
69 Compute_h_eta_AtCellCenter (const int &i, const int &j, const int &k,
70  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
71  const amrex::Array4<const amrex::Real>& z_nd)
72 {
73  amrex::Real dyInv = cellSizeInv[1];
74  amrex::Real met_h_eta = 0.25 * dyInv *
75  ( z_nd(i,j+1,k) + z_nd(i+1,j+1,k) + z_nd(i,j+1,k+1) + z_nd(i+1,j+1,k+1)
76  -z_nd(i,j ,k) - z_nd(i+1,j ,k) - z_nd(i,j ,k+1) - z_nd(i+1,j ,k+1) );
77  return met_h_eta;
78 }
79 
80 //*****************************************************************************************
81 // Compute terrain metric terms at face-centers
82 //*****************************************************************************************
83 // Metric coincides with U location
84 AMREX_GPU_DEVICE
85 AMREX_FORCE_INLINE
86 amrex::Real
87 Compute_h_zeta_AtIface (const int &i, const int &j, const int &k,
88  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
89  const amrex::Array4<const amrex::Real>& z_nd)
90 {
91  amrex::Real dzInv = cellSizeInv[2];
92  amrex::Real met_h_zeta = 0.5 * dzInv * ( z_nd(i,j,k+1) + z_nd(i,j+1,k+1)
93  - z_nd(i,j,k ) - z_nd(i,j+1,k ) );
94  return met_h_zeta;
95 }
96 
97 // Metric coincides with U location
98 AMREX_GPU_DEVICE
99 AMREX_FORCE_INLINE
100 amrex::Real
101 Compute_h_xi_AtIface (const int &i, const int &j, const int &k,
102  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
103  const amrex::Array4<const amrex::Real>& z_nd)
104 {
105  amrex::Real dxInv = cellSizeInv[0];
106  amrex::Real met_h_xi = 0.125 * dxInv *
107  ( z_nd(i+1,j,k) + z_nd(i+1,j,k+1) + z_nd(i+1,j+1,k) + z_nd(i+1,j+1,k+1)
108  -z_nd(i-1,j,k) - z_nd(i-1,j,k+1) - z_nd(i-1,j+1,k) - z_nd(i-1,j+1,k+1) );
109  return met_h_xi;
110 }
111 
112 // Metric coincides with U location
113 AMREX_GPU_DEVICE
114 AMREX_FORCE_INLINE
115 amrex::Real
116 Compute_h_eta_AtIface (const int &i, const int &j, const int &k,
117  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
118  const amrex::Array4<const amrex::Real>& z_nd)
119 {
120  amrex::Real dyInv = cellSizeInv[1];
121  amrex::Real met_h_eta = 0.5 * dyInv * ( z_nd(i,j+1,k ) + z_nd(i,j+1,k+1)
122  - z_nd(i,j ,k ) - z_nd(i,j ,k+1) );
123  return met_h_eta;
124 }
125 
126 // Metric coincides with V location
127 AMREX_GPU_DEVICE
128 AMREX_FORCE_INLINE
129 amrex::Real
130 Compute_h_zeta_AtJface (const int &i, const int &j, const int &k,
131  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
132  const amrex::Array4<const amrex::Real>& z_nd)
133 {
134  amrex::Real dzInv = cellSizeInv[2];
135  amrex::Real met_h_zeta = 0.5 * dzInv * ( z_nd(i,j,k+1) + z_nd(i+1,j,k+1)
136  - z_nd(i,j,k ) - z_nd(i+1,j,k ) );
137  return met_h_zeta;
138 }
139 
140 // Metric coincides with V location
141 AMREX_GPU_DEVICE
142 AMREX_FORCE_INLINE
143 amrex::Real
144 Compute_h_xi_AtJface (const int &i, const int &j, const int &k,
145  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
146  const amrex::Array4<const amrex::Real>& z_nd)
147 {
148  amrex::Real dxInv = cellSizeInv[0];
149  amrex::Real met_h_xi = 0.5 * dxInv * ( z_nd(i+1,j,k) + z_nd(i+1,j,k+1)
150  -z_nd(i ,j,k) - z_nd(i ,j,k+1) );
151  return met_h_xi;
152 }
153 
154 // Metric coincides with V location
155 AMREX_GPU_DEVICE
156 AMREX_FORCE_INLINE
157 amrex::Real
158 Compute_h_eta_AtJface (const int &i, const int &j, const int &k,
159  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
160  const amrex::Array4<const amrex::Real>& z_nd)
161 {
162  amrex::Real dyInv = cellSizeInv[1];
163  amrex::Real met_h_eta = 0.125 * dyInv *
164  ( z_nd(i,j+1,k) + z_nd(i,j+1,k+1) + z_nd(i+1,j+1,k) + z_nd(i+1,j+1,k+1)
165  -z_nd(i,j-1,k) - z_nd(i,j-1,k+1) - z_nd(i+1,j-1,k) - z_nd(i+1,j-1,k+1) );
166  return met_h_eta;
167 }
168 
169 // Metric coincides with K location
170 AMREX_GPU_DEVICE
171 AMREX_FORCE_INLINE
172 amrex::Real
173 Compute_h_zeta_AtKface (const int &i, const int &j, const int &k,
174  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
175  const amrex::Array4<const amrex::Real>& z_nd)
176 {
177  amrex::Real dzInv = cellSizeInv[2];
178  amrex::Real met_h_zeta = 0.125 * dzInv *
179  ( z_nd(i,j,k+1) + z_nd(i+1,j,k+1) + z_nd(i,j+1,k+1) + z_nd(i+1,j+1,k+1)
180  -z_nd(i,j,k-1) - z_nd(i+1,j,k-1) - z_nd(i,j+1,k-1) - z_nd(i+1,j+1,k-1) );
181  return met_h_zeta;
182 }
183 
184 // Metric coincides with K location
185 AMREX_GPU_DEVICE
186 AMREX_FORCE_INLINE
187 amrex::Real
188 Compute_h_xi_AtKface (const int &i, const int &j, const int &k,
189  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
190  const amrex::Array4<const amrex::Real>& z_nd)
191 {
192  amrex::Real dxInv = cellSizeInv[0];
193  amrex::Real met_h_xi = 0.5 * dxInv * ( z_nd(i+1,j,k) + z_nd(i+1,j+1,k)
194  -z_nd(i ,j,k) - z_nd(i ,j+1,k) );
195  return met_h_xi;
196 }
197 
198 // Metric coincides with K location
199 AMREX_GPU_DEVICE
200 AMREX_FORCE_INLINE
201 amrex::Real
202 Compute_h_eta_AtKface (const int &i, const int &j, const int &k,
203  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
204  const amrex::Array4<const amrex::Real>& z_nd)
205 {
206  amrex::Real dyInv = cellSizeInv[1];
207  amrex::Real met_h_eta = 0.5 * dyInv * ( z_nd(i,j+1,k) + z_nd(i+1,j+1,k)
208  -z_nd(i,j ,k) - z_nd(i+1,j ,k) );
209  return met_h_eta;
210 }
211 
212 //*****************************************************************************************
213 // Compute terrain metric terms at edge-centers
214 //*****************************************************************************************
215 // -- EdgeCenterK --
216 
217 // Metric is at edge and center Z (red pentagon)
218 AMREX_GPU_DEVICE
219 AMREX_FORCE_INLINE
220 amrex::Real
221 Compute_h_zeta_AtEdgeCenterK (const int &i, const int &j, const int &k,
222  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
223  const amrex::Array4<const amrex::Real>& z_nd)
224 {
225  amrex::Real dzInv = cellSizeInv[2];
226  amrex::Real met_h_zeta = dzInv * (z_nd(i,j,k+1) - z_nd(i,j,k));
227  return met_h_zeta;
228 }
229 
230 // Metric is at edge and center Z (red pentagon)
231 AMREX_GPU_DEVICE
232 AMREX_FORCE_INLINE
233 amrex::Real
234 Compute_h_xi_AtEdgeCenterK (const int &i, const int &j, const int &k,
235  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
236  const amrex::Array4<const amrex::Real>& z_nd)
237 {
238  amrex::Real dxInv = cellSizeInv[0];
239  amrex::Real met_h_xi = 0.25 * dxInv *
240  ( z_nd(i+1,j,k) + z_nd(i+1,j,k+1)
241  -z_nd(i-1,j,k) - z_nd(i-1,j,k+1) );
242  return met_h_xi;
243 }
244 
245 // Metric is at edge and center Z (red pentagon)
246 AMREX_GPU_DEVICE
247 AMREX_FORCE_INLINE
248 amrex::Real
249 Compute_h_eta_AtEdgeCenterK (const int &i, const int &j, const int &k,
250  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
251  const amrex::Array4<const amrex::Real>& z_nd)
252 {
253  amrex::Real dyInv = cellSizeInv[1];
254  amrex::Real met_h_eta = 0.25 * dyInv *
255  ( z_nd(i,j+1,k) + z_nd(i,j+1,k+1)
256  -z_nd(i,j-1,k) - z_nd(i,j-1,k+1) );
257  return met_h_eta;
258 }
259 
260 // -- EdgeCenterJ --
261 
262 // Metric is at edge and center Y (magenta cross)
263 AMREX_GPU_DEVICE
264 AMREX_FORCE_INLINE
265 amrex::Real
266 Compute_h_zeta_AtEdgeCenterJ (const int &i, const int &j, const int &k,
267  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
268  const amrex::Array4<const amrex::Real>& z_nd)
269 {
270  amrex::Real dzInv = cellSizeInv[2];
271  amrex::Real met_h_zeta = 0.25 * dzInv * ( z_nd(i,j,k+1) + z_nd(i,j+1,k+1)
272  -z_nd(i,j,k-1) - z_nd(i,j+1,k-1) );
273  return met_h_zeta;
274 }
275 
276 // Metric is at edge and center Y (magenta cross)
277 AMREX_GPU_DEVICE
278 AMREX_FORCE_INLINE
279 amrex::Real
280 Compute_h_xi_AtEdgeCenterJ (const int &i, const int &j, const int &k,
281  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
282  const amrex::Array4<const amrex::Real>& z_nd)
283 {
284  amrex::Real dxInv = cellSizeInv[0];
285  amrex::Real met_h_xi = 0.25 * dxInv *
286  ( z_nd(i+1,j+1,k) + z_nd(i+1,j ,k)
287  -z_nd(i-1,j+1,k) - z_nd(i-1,j ,k) );
288  return met_h_xi;
289 }
290 
291 // Metric is at edge and center Y (magenta cross)
292 AMREX_GPU_DEVICE
293 AMREX_FORCE_INLINE
294 amrex::Real
295 Compute_h_eta_AtEdgeCenterJ (const int &i, const int &j, const int &k,
296  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
297  const amrex::Array4<const amrex::Real>& z_nd)
298 {
299  amrex::Real dyInv = cellSizeInv[1];
300  amrex::Real met_h_eta = dyInv * ( z_nd(i,j+1,k) - z_nd(i,j,k) );
301  return met_h_eta;
302 }
303 
304 // -- EdgeCenterI --
305 
306 // Metric is at edge and center Y (magenta cross)
307 AMREX_GPU_DEVICE
308 AMREX_FORCE_INLINE
309 amrex::Real
310 Compute_h_zeta_AtEdgeCenterI (const int &i, const int &j, const int &k,
311  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
312  const amrex::Array4<const amrex::Real>& z_nd)
313 {
314  amrex::Real dzInv = cellSizeInv[2];
315  amrex::Real met_h_zeta = 0.25 * dzInv * ( z_nd(i,j,k+1) + z_nd(i+1,j,k+1)
316  -z_nd(i,j,k-1) - z_nd(i+1,j,k-1) );
317  return met_h_zeta;
318 }
319 
320 // Metric is at edge and center Y (magenta cross)
321 AMREX_GPU_DEVICE
322 AMREX_FORCE_INLINE
323 amrex::Real
324 Compute_h_xi_AtEdgeCenterI (const int &i, const int &j, const int &k,
325  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
326  const amrex::Array4<const amrex::Real>& z_nd)
327 {
328  amrex::Real dxInv = cellSizeInv[0];
329  amrex::Real met_h_xi = dxInv * ( z_nd(i+1,j,k) - z_nd(i,j,k) );
330  return met_h_xi;
331 }
332 
333 // Metric is at edge and center Y (magenta cross)
334 AMREX_GPU_DEVICE
335 AMREX_FORCE_INLINE
336 amrex::Real
337 Compute_h_eta_AtEdgeCenterI (const int &i, const int &j, const int &k,
338  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
339  const amrex::Array4<const amrex::Real>& z_nd)
340 {
341  amrex::Real dyInv = cellSizeInv[1];
342  amrex::Real met_h_eta = 0.25 * dyInv *
343  ( z_nd(i+1,j+1,k) + z_nd(i,j+1,k)
344  -z_nd(i+1,j-1,k) - z_nd(i,j-1,k) );
345  return met_h_eta;
346 }
347 
348 // Relative height above terrain surface at cell center from z_nd (nodal absolute height)
349 AMREX_GPU_DEVICE
350 AMREX_FORCE_INLINE
351 amrex::Real
352 Compute_Zrel_AtCellCenter (const int &i, const int &j, const int &k,
353  const amrex::Array4<const amrex::Real>& z_nd)
354 {
355  const amrex::Real z_cc = 0.125*( z_nd(i ,j ,k ) + z_nd(i ,j ,k+1) +
356  + z_nd(i+1,j ,k ) + z_nd(i ,j ,k+1)
357  + z_nd(i ,j+1,k ) + z_nd(i ,j+1,k+1)
358  + z_nd(i+1,j+1,k ) + z_nd(i ,j+1,k+1));
359 
360  // Note: we assume the z_nd array spans from the bottom to top of the domain
361  // i.e. no domain decomposition across processors in vertical direction
362  const amrex::Real z0_cc = 0.25*( z_nd(i ,j ,0) + z_nd(i ,j+1,0)
363  + z_nd(i+1,j ,0) + z_nd(i+1,j+1,0));
364 
365  return (z_cc - z0_cc);
366 }
367 
368 /**
369  * Define omega given u,v and w
370  */
371 AMREX_GPU_DEVICE
372 AMREX_FORCE_INLINE
373 amrex::Real
374 OmegaFromW (int i, int j, int k, amrex::Real w,
375  const amrex::Array4<const amrex::Real> u_arr,
376  const amrex::Array4<const amrex::Real> v_arr,
377  const amrex::Array4<const amrex::Real> z_nd,
378  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& dxInv)
379 {
380  // This is dh/dxi at z-face (i,j,k-1/2)
381  amrex::Real met_zlo_xi = 0.5 * dxInv[0] *
382  ( z_nd(i+1,j+1,k ) + z_nd(i+1,j ,k ) // hi i, lo k
383  -z_nd(i ,j+1,k ) - z_nd(i ,j ,k ) ); // lo i, lo k
384 
385  // This is dh/deta at z-face (i,j,k-1/2)
386  amrex::Real met_zlo_eta = 0.5 * dxInv[1] *
387  ( z_nd(i+1,j+1,k ) + z_nd(i ,j+1,k ) // hi j, lo k
388  -z_nd(i+1,j ,k ) - z_nd(i ,j ,k ) ); // lo j, lo k
389 
390  // Slip BC or moving terrain
391  // Use extrapolation instead of interpolation if at the bottom boundary
392  amrex::Real u = (k == 0) ? 1.5 * (0.5*(u_arr(i,j,k)+u_arr(i+1,j,k))) - 0.5*(0.5*(u_arr(i,j,k+1)+u_arr(i+1,j,k+1))) :
393  0.25 * ( u_arr(i,j,k-1) + u_arr(i+1,j,k-1) + u_arr(i,j,k) + u_arr(i+1,j,k) );
394  amrex::Real v = (k == 0) ? 1.5 * (0.5*(v_arr(i,j,k)+v_arr(i,j+1,k))) - 0.5*(0.5*(v_arr(i,j,k+1)+v_arr(i,j+1,k+1))) :
395  0.25 * ( v_arr(i,j,k-1) + v_arr(i,j+1,k-1) + v_arr(i,j,k) + v_arr(i,j+1,k) );
396 
397  amrex::Real omega = w - met_zlo_xi * u - met_zlo_eta * v;
398  return omega;
399 }
400 
401 /**
402  * Define w given scalar u,v and omega
403  */
404 AMREX_GPU_DEVICE
405 AMREX_FORCE_INLINE
406 amrex::Real
407 WFromOmega (int i, int j, int k, amrex::Real omega,
408  amrex::Real u, amrex::Real v,
409  const amrex::Array4<const amrex::Real>& z_nd,
410  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& dxInv)
411 {
412  // This is dh/dxi at z-face (i,j,k-1/2)
413  amrex::Real met_zlo_xi = 0.5 * dxInv[0] *
414  ( z_nd(i+1,j+1,k ) + z_nd(i+1,j ,k ) // hi i, lo k
415  -z_nd(i ,j+1,k ) - z_nd(i ,j ,k ) ); // lo i, lo k
416 
417  // This is dh/deta at z-face (i,j,k-1/2)
418  amrex::Real met_zlo_eta = 0.5 * dxInv[1] *
419  ( z_nd(i+1,j+1,k ) + z_nd(i ,j+1,k ) // hi j, lo k
420  -z_nd(i+1,j ,k ) - z_nd(i ,j ,k ) ); // lo j, lo k
421 
422  amrex::Real w = omega + met_zlo_xi * u + met_zlo_eta * v;
423  return w;
424 }
425 
426 /**
427  * Define w given u and v arrays and scalar omega
428  */
429 AMREX_GPU_DEVICE
430 AMREX_FORCE_INLINE
431 amrex::Real
432 WFromOmega (int i, int j, int k, amrex::Real omega,
433  const amrex::Array4<const amrex::Real>& u_arr,
434  const amrex::Array4<const amrex::Real>& v_arr,
435  const amrex::Array4<const amrex::Real>& z_nd,
436  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& dxInv)
437 {
438  // Use extrapolation instead of interpolation if at the bottom boundary
439  amrex::Real u = (k == 0) ? 1.5 * (0.5*(u_arr(i,j,k)+u_arr(i+1,j,k))) - 0.5*(0.5*(u_arr(i,j,k+1)+u_arr(i+1,j,k+1))) :
440  0.25 * ( u_arr(i,j,k-1) + u_arr(i+1,j,k-1) + u_arr(i,j,k) + u_arr(i+1,j,k) );
441  amrex::Real v = (k == 0) ? 1.5 * (0.5*(v_arr(i,j,k)+v_arr(i,j+1,k))) - 0.5*(0.5*(v_arr(i,j,k+1)+v_arr(i,j+1,k+1))) :
442  0.25 * ( v_arr(i,j,k-1) + v_arr(i,j+1,k-1) + v_arr(i,j,k) + v_arr(i,j+1,k) );
443 
444  amrex::Real w = WFromOmega(i,j,k,omega,u,v,z_nd,dxInv);
445  return w;
446 }
447 
448 //*****************************************************************************************
449 // Rotate scalar flux vector and stress tensor for MOST
450 //*****************************************************************************************
451 AMREX_GPU_DEVICE
452 AMREX_FORCE_INLINE
453 void
454 rotate_scalar_flux (const int& i,
455  const int& j,
456  const int& klo,
457  const amrex::Real& flux,
458  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& dxInv,
459  const amrex::Array4<const amrex::Real>& zphys_arr,
460  const amrex::Array4<amrex::Real>& phi1_arr,
461  const amrex::Array4<amrex::Real>& phi2_arr,
462  const amrex::Array4<amrex::Real>& phi3_arr)
463 {
464  amrex::Real h_xi = Compute_h_xi_AtCellCenter(i, j, klo, dxInv, zphys_arr);
465  amrex::Real h_eta = Compute_h_eta_AtCellCenter(i, j, klo, dxInv, zphys_arr);
466  amrex::Real InvNorm = 1.0 / std::sqrt(1.0 + h_xi*h_xi + h_eta*h_eta);
467  phi1_arr(i,j,klo) = -h_xi * flux * InvNorm;
468  phi2_arr(i,j,klo) = -h_eta * flux * InvNorm;
469  phi3_arr(i,j,klo) = flux * InvNorm;
470 }
471 
472 AMREX_GPU_DEVICE
473 AMREX_FORCE_INLINE
474 void
475 rotate_stress_tensor (const int& i,
476  const int& j,
477  const int& klo,
478  const amrex::Real& flux,
479  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& dxInv,
480  const amrex::Array4<const amrex::Real>& zphys_arr,
481  const amrex::Array4<const amrex::Real>& u_arr,
482  const amrex::Array4<const amrex::Real>& v_arr,
483  const amrex::Array4<const amrex::Real>& w_arr,
484  const amrex::Array4<amrex::Real>& tau11_arr,
485  const amrex::Array4<amrex::Real>& tau22_arr,
486  const amrex::Array4<amrex::Real>& tau33_arr,
487  const amrex::Array4<amrex::Real>& tau12_arr,
488  const amrex::Array4<amrex::Real>& tau21_arr,
489  const amrex::Array4<amrex::Real>& tau13_arr,
490  const amrex::Array4<amrex::Real>& tau31_arr,
491  const amrex::Array4<amrex::Real>& tau23_arr,
492  const amrex::Array4<amrex::Real>& tau32_arr)
493 {
494  // Unit-normal vector
495  amrex::Array1D<amrex::Real,0,2> n_hat;
496 
497  // Unit-tangent vectors
498  amrex::Array1D<amrex::Real,0,2> t_hat_1;
499  amrex::Array1D<amrex::Real,0,2> t_hat_2;
500  amrex::Array1D<amrex::Real,0,2> u_t_hat;
501 
502  // Final basis vector for right-hand coord sys
503  amrex::Array1D<amrex::Real,0,2> a_hat;
504 
505  // Rotation matrix
506  amrex::Array2D<amrex::Real,0,2,0,2> R_mat;
507 
508  // Metric data
509  amrex::Real h_xi = Compute_h_xi_AtIface(i, j, klo, dxInv, zphys_arr);
510  amrex::Real h_eta = Compute_h_eta_AtIface(i, j, klo, dxInv, zphys_arr);
511 
512  // Populate the normal vector
513  amrex::Real Inormn = 1./std::sqrt(1.0 + h_xi*h_xi + h_eta*h_eta);
514  n_hat(0) = -Inormn*h_xi; n_hat(1) = -Inormn*h_eta; n_hat(2) = Inormn;
515 
516  // Populate the tangent vectors
517  amrex::Real Inorm1 = 1./std::sqrt(1.0 + h_xi*h_xi);
518  amrex::Real Inorm2 = 1./std::sqrt(1.0 + h_eta*h_eta);
519  t_hat_1(0) = Inorm1; t_hat_2(1) = Inorm2;
520  t_hat_1(2) = Inorm1*h_xi; t_hat_2(2) = Inorm2*h_eta;
521 
522  // Populate the u_t vector
523  amrex::Real Norm_u_t = 0.0;
524  amrex::Real mag1 = (u_arr(i,j,klo) + h_xi *w_arr(i,j,klo))*Inorm1;
525  amrex::Real mag2 = (v_arr(i,j,klo) + h_eta*w_arr(i,j,klo))*Inorm2;
526  for (int icol(0); icol<3; ++icol) {
527  u_t_hat(icol) = mag1*t_hat_1(icol) + mag2*t_hat_2(icol);
528  Norm_u_t += u_t_hat(icol)*u_t_hat(icol);
529  }
530  for (int icol(0); icol<3; ++icol) {
531  u_t_hat(icol) /= std::sqrt(Norm_u_t);
532  }
533 
534  // Populate the a_hat vector
535  a_hat(0) = n_hat(1)*u_t_hat(2) - n_hat(2)*u_t_hat(1);
536  a_hat(1) = -(n_hat(0)*u_t_hat(2) - n_hat(2)*u_t_hat(0));
537  a_hat(2) = n_hat(1)*u_t_hat(1) - n_hat(1)*u_t_hat(0);
538 
539  // Copy column vectors into R_mat
540  int jrow;
541  jrow = 0;
542  for (int icol(0); icol<3; ++icol) {
543  R_mat(icol,jrow) = u_t_hat(icol);
544  }
545  jrow = 1;
546  for (int icol(0); icol<3; ++icol) {
547  R_mat(icol,jrow) = a_hat(icol);
548  }
549  jrow = 2;
550  for (int icol(0); icol<3; ++icol) {
551  R_mat(icol,jrow) = n_hat(icol);
552  }
553 
554  // Assign the stresses
555  tau11_arr(i,j,klo) = 2.0*R_mat(0,0)*R_mat(2,0)*flux;
556  tau22_arr(i,j,klo) = 2.0*R_mat(0,1)*R_mat(2,1)*flux;
557  tau33_arr(i,j,klo) = 2.0*R_mat(0,2)*R_mat(2,2)*flux;
558 
559  tau12_arr(i,j,klo) = (R_mat(0,0)*R_mat(2,1) + R_mat(2,0)*R_mat(0,1))*flux;
560  tau21_arr(i,j,klo) = tau12_arr(i,j,klo);
561 
562  tau13_arr(i,j,klo) = (R_mat(0,0)*R_mat(2,2) + R_mat(0,2)*R_mat(2,0))*flux;
563  tau31_arr(i,j,klo) = tau13_arr(i,j,klo);
564 
565  tau23_arr(i,j,klo) = (R_mat(0,1)*R_mat(2,2) + R_mat(0,2)*R_mat(2,1))*flux;
566  tau32_arr(i,j,klo) = tau23_arr(i,j,klo);
567 }
568 #endif
ERF_IndexDefines.H
Compute_h_eta_AtEdgeCenterJ
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_eta_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:295
Compute_h_xi_AtEdgeCenterK
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_xi_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:234
Compute_h_xi_AtEdgeCenterJ
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_xi_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:280
Compute_h_xi_AtIface
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_xi_AtIface(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:101
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:266
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:221
init_zlevels
void init_zlevels(amrex::Vector< amrex::Vector< amrex::Real >> &zlevels_stag, amrex::Vector< amrex::Vector< amrex::Real >> &stretched_dz_h, amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real >> &stretched_dz_d, amrex::Vector< amrex::Geometry > const &geom, amrex::Vector< amrex::IntVect > const &ref_ratio, const amrex::Real grid_stretching_ratio, const amrex::Real zsurf, const amrex::Real dz0)
Compute_h_zeta_AtCellCenter
AMREX_FORCE_INLINE AMREX_GPU_DEVICE amrex::Real Compute_h_zeta_AtCellCenter(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:39
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:310
Compute_Zrel_AtCellCenter
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_Zrel_AtCellCenter(const int &i, const int &j, const int &k, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:352
Compute_h_eta_AtCellCenter
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_eta_AtCellCenter(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:69
Compute_h_zeta_AtKface
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_zeta_AtKface(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:173
WFromOmega
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real WFromOmega(int i, int j, int k, amrex::Real omega, amrex::Real u, amrex::Real v, const amrex::Array4< const amrex::Real > &z_nd, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &dxInv)
Definition: ERF_TerrainMetrics.H:407
rotate_scalar_flux
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void rotate_scalar_flux(const int &i, const int &j, const int &klo, const amrex::Real &flux, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &dxInv, const amrex::Array4< const amrex::Real > &zphys_arr, const amrex::Array4< amrex::Real > &phi1_arr, const amrex::Array4< amrex::Real > &phi2_arr, const amrex::Array4< amrex::Real > &phi3_arr)
Definition: ERF_TerrainMetrics.H:454
OmegaFromW
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real OmegaFromW(int i, int j, int k, amrex::Real w, const amrex::Array4< const amrex::Real > u_arr, const amrex::Array4< const amrex::Real > v_arr, const amrex::Array4< const amrex::Real > z_nd, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &dxInv)
Definition: ERF_TerrainMetrics.H:374
Compute_h_zeta_AtIface
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_zeta_AtIface(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:87
rotate_stress_tensor
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void rotate_stress_tensor(const int &i, const int &j, const int &klo, const amrex::Real &flux, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &dxInv, const amrex::Array4< const amrex::Real > &zphys_arr, const amrex::Array4< const amrex::Real > &u_arr, const amrex::Array4< const amrex::Real > &v_arr, const amrex::Array4< const amrex::Real > &w_arr, const amrex::Array4< amrex::Real > &tau11_arr, const amrex::Array4< amrex::Real > &tau22_arr, const amrex::Array4< amrex::Real > &tau33_arr, const amrex::Array4< amrex::Real > &tau12_arr, const amrex::Array4< amrex::Real > &tau21_arr, const amrex::Array4< amrex::Real > &tau13_arr, const amrex::Array4< amrex::Real > &tau31_arr, const amrex::Array4< amrex::Real > &tau23_arr, const amrex::Array4< amrex::Real > &tau32_arr)
Definition: ERF_TerrainMetrics.H:475
Compute_h_xi_AtJface
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_xi_AtJface(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:144
Compute_h_xi_AtKface
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_xi_AtKface(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:188
Compute_h_eta_AtEdgeCenterK
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_eta_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:249
init_which_terrain_grid
void init_which_terrain_grid(int lev, const amrex::Geometry &geom, amrex::MultiFab &z_phys_nd, amrex::Vector< amrex::Real > const &z_levels_h)
Compute_h_eta_AtIface
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_eta_AtIface(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:116
init_terrain_grid
void init_terrain_grid(int lev, const amrex::Geometry &geom, amrex::MultiFab &z_phys_nd, amrex::Vector< amrex::Real > const &z_levels_h, amrex::GpuArray< ERF_BC, AMREX_SPACEDIM *2 > &phys_bc_type)
Compute_h_xi_AtEdgeCenterI
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_xi_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:324
Compute_h_zeta_AtJface
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_zeta_AtJface(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:130
Compute_h_xi_AtCellCenter
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_xi_AtCellCenter(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:54
Compute_h_eta_AtJface
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_eta_AtJface(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:158
Compute_h_eta_AtEdgeCenterI
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_eta_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:337
Compute_h_eta_AtKface
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_eta_AtKface(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:202
MicVar::omega
@ omega
Definition: ERF_SAM.H:49