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