ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Functions
ERF_InitCustomTerrain.cpp File Reference
#include "AMReX_ParmParse.H"
#include "ERF_Constants.H"
#include "ERF_TerrainMetrics.H"
Include dependency graph for ERF_InitCustomTerrain.cpp:

Functions

void init_my_custom_terrain (const Geometry &geom, FArrayBox &terrain_fab, const double &time_d)
 

Function Documentation

◆ init_my_custom_terrain()

void init_my_custom_terrain ( const Geometry &  geom,
FArrayBox &  terrain_fab,
const double &  time_d 
)
15 {
16  //
17  // We put this here as a convenience for testing the map factor implementation
18  // Note that these factors must match those in Source/ERF_MakeNewArrays.cpp
19  //
20  ParmParse pp("erf");
21  bool test_mapfactor = false;
22  pp.query("test_mapfactor",test_mapfactor);
23 
24  Real mf_m;
25  if (test_mapfactor) {
26  mf_m = myhalf;
27  } else {
28  mf_m = one;
29  }
30 
31  // Domain cell size and real bounds
32  auto dx = geom.CellSizeArray();
33  auto ProbLoArr = geom.ProbLoArray();
34  auto ProbHiArr = geom.ProbHiArray();
35 
36  const amrex::Box& domain = geom.Domain();
37  int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1;
38  int domlo_y = domain.smallEnd(1); int domhi_y = domain.bigEnd(1) + 1;
39  int domlo_z = domain.smallEnd(2);
40 
41  // User function parameters
42  Real a = myhalf;
43  Real num = Real(8.) * a * a * a;
44  Real xcen = myhalf * (ProbLoArr[0] + ProbHiArr[0]) / mf_m;
45  Real ycen = myhalf * (ProbLoArr[1] + ProbHiArr[1]) / mf_m;
46 
47  // Populate bottom plane
48  int k0 = domlo_z;
49 
50  std::string custom_terrain_type = "None";
51  ParmParse pp_prob("prob"); pp_prob.query("custom_terrain_type", custom_terrain_type);
52 
53  amrex::Box zbx = terrain_fab.box();
54  if (zbx.smallEnd(2) <= k0)
55  {
56  amrex::Array4<Real> const& z_arr = terrain_fab.array();
57 
58  if (custom_terrain_type == "WoA") {
59 
60  // Default to x-direction
61  int dir = 0; pp_prob.query("dir", dir);
62 
63  Real L = Real(100.0); pp_prob.query("L" , L);
64  Real z_offset = zero; pp_prob.query("z_offset" , z_offset);
65 
66  // If hm is nonzero, then use alternate hill definition
67  Real hm = zero; pp_prob.query("hmax" , hm);
68 
69  // This is a 2D hill with variation in only the x-direction
70  if (dir == 0) {
71  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int)
72  {
73  // Clip indices for ghost-cells
74  int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);
75 
76  // Location of nodes
77  Real x = (ProbLoArr[0] + ii * dx[0] - xcen) * mf_m;
78 
79  // WoA Hill in x-direction
80  if (hm==0) {
81  z_arr(i,j,k0) = num / (x*x + 4 * a * a);
82  } else {
83  Real x_L = x / L;
84  z_arr(i,j,k0) = hm / (1 + x_L*x_L) + z_offset;
85  }
86  });
87  } else if (dir == 1) {
88  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int)
89  {
90  // Clip indices for ghost-cells
91  int jj = amrex::min(amrex::max(j,domlo_y),domhi_y);
92 
93  // Location of nodes
94  Real y = (ProbLoArr[1] + jj * dx[1] - ycen) * mf_m;
95 
96  // WoA Hill in y-direction
97  if (hm==0) {
98  z_arr(i,j,k0) = num / (y*y + Real(4.0) * a * a);
99  } else {
100  Real y_L = y / L;
101  z_arr(i,j,k0) = hm / (one + y_L*y_L) + z_offset;
102  }
103  });
104  } else if (dir == 2) {
105  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int)
106  {
107  // Clip indices for ghost-cells
108  int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);
109  int jj = amrex::min(amrex::max(j,domlo_y),domhi_y);
110 
111  // Location of nodes
112  Real x = (ProbLoArr[0] + ii * dx[0] - xcen) * mf_m;
113  Real y = (ProbLoArr[1] + jj * dx[1] - ycen) * mf_m;
114  Real r = std::sqrt(x*x + y*y);
115 
116  // WoA Hill in radial direction
117  if (hm==0) {
118  z_arr(i,j,k0) = num / (r*r + Real(4.0) * a * a);
119  } else {
120  Real r_L = r / L;
121  z_arr(i,j,k0) = hm / (one + r_L*r_L) + z_offset;
122  }
123  });
124  } else {
125  amrex::Abort("Unknown dir in ERF_Prob.cpp");
126  }
127 
128  } else if (custom_terrain_type == "ScharMountain") {
129 
130  Real asq = Real(5000.0) * Real(5000.0);
131  Real Hm = Real(250.0);
132  Real lambda = Real(4000.0);
133 
134  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int)
135  {
136  // Clip indices for ghost-cells
137  int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);
138 
139  // Location of nodes
140  Real x = (ProbLoArr[0] + ii * dx[0] - xcen);
141 
142  Real cosx = std::cos(PI * x / lambda);
143 
144  z_arr(i,j,k0) = Hm * std::exp(-x*x/asq) * cosx * cosx;
145  });
146 
147  } else if (custom_terrain_type == "HalfCylinder") {
148 
149  Real asq = myhalf * myhalf;
150 
151  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int)
152  {
153  // Clip indices for ghost-cells
154  int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);
155 
156  // Location of nodes
157  Real x = (ProbLoArr[0] + ii * dx[0] - xcen);
158 
159  Real rsq = x*x;
160 
161  if (rsq < asq) {
162  z_arr(i,j,k0) = std::sqrt(asq - rsq);
163  } else {
164  z_arr(i,j,k0) = zero;
165  }
166  });
167 
168  } else if (custom_terrain_type == "Hemisphere") {
169 
170  Real asq = myhalf * myhalf;
171 
172  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int)
173  {
174  // Clip indices for ghost-cells
175  int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);
176  int jj = amrex::min(amrex::max(j,domlo_y),domhi_y);
177 
178  // Location of nodes
179  Real x = (ProbLoArr[0] + ii * dx[0] - xcen);
180  Real y = (ProbLoArr[1] + jj * dx[1] - ycen);
181 
182  Real rsq = x*x + y*y;
183 
184  if (rsq < asq) {
185  z_arr(i,j,k0) = std::pow(asq-rsq, myhalf);
186  } else {
187  z_arr(i,j,k0) = zero;
188  }
189  });
190 
191  } else if (custom_terrain_type == "MovingSineWave") {
192 
193  Real Ampl = zero; pp_prob.query("Ampl", Ampl);
194  Real wavelength = Real(100.); pp_prob.query("wavelength", wavelength);
195 
196  Real kp = two * PI / wavelength;
197  Real g = CONST_GRAV;
198  Real omega = std::sqrt(g * kp);
199 
200  Real time = static_cast<Real>(time_d);
201 
202  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int)
203  {
204  // Clip indices for ghost-cells
205  int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);
206 
207  // Location of nodes
208  Real x = ii * dx[0];
209 
210  // Wave height
211  Real height = Ampl * std::sin(kp * x - omega * time);
212 
213  // Populate terrain height
214  z_arr(i,j,0) = height;
215  });
216 
217  } else if (custom_terrain_type == "WindFarmTest") {
218 
219  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int)
220  {
221  // Clip indices for ghost-cells
222  int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);
223  int jj = amrex::min(amrex::max(j,domlo_y),domhi_y);
224 
225  // Location of nodes
226  Real x = (ProbLoArr[0] + ii * dx[0] - xcen);
227  Real y = (ProbLoArr[1] + jj * dx[1] - ycen);
228 
229  Real x_L = x/Real(100.0);
230  Real y_L = y/Real(100.0);
231 
232  z_arr(i,j,k0) = Real(100.0) / (one + x_L*x_L + y_L*y_L);
233  });
234 
235  } else if (custom_terrain_type == "RaisedFlat") {
236  Real z_offset = zero; pp_prob.query("z_offset" , z_offset);
237  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int)
238  {
239  z_arr(i,j,k0) = z_offset;
240  });
241  } else if (custom_terrain_type == "Cos4Hill") {
242 
243  // Get prob parameters (must be outside GPU kernel)
244  Real hm = zero; pp_prob.query("hmax", hm);
245  Real L = Real(100.0); pp_prob.query("L", L);
246  Real z_offset = zero; pp_prob.query("z_offset", z_offset);
247  Real fourL = four * L;
248 
249  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int)
250  {
251 
252  // Clip indices for ghost-cells
253  int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);
254  int jj = amrex::min(amrex::max(j,domlo_y),domhi_y);
255 
256  // Location of nodes
257  Real x = (ProbLoArr[0] + ii * dx[0] - xcen);
258  Real y = (ProbLoArr[1] + jj * dx[1] - ycen);
259  Real r = std::sqrt(x*x + y*y);
260 
261  if (r < fourL) {
262  z_arr(i,j,k0) = z_offset + hm * Real(0.0625) * std::pow(one + std::cos(PI*r/fourL), four);
263  } else {
264  z_arr(i,j,k0) = z_offset;
265  }
266  });
267 
268  } else if (custom_terrain_type == "None") {
269  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int)
270  {
271  z_arr(i,j,k0) = zero;
272  });
273  } else {
274  Abort("Don't know this custom_terrain_type");
275  }
276  }
277 }
four
constexpr amrex::Real four
Definition: ERF_Constants.H:12
two
constexpr amrex::Real two
Definition: ERF_Constants.H:10
one
constexpr amrex::Real one
Definition: ERF_Constants.H:9
zero
constexpr amrex::Real zero
Definition: ERF_Constants.H:8
myhalf
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:13
PI
constexpr amrex::Real PI
Definition: ERF_Constants.H:42
CONST_GRAV
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:64
num
@ num
Definition: ERF_DataStruct.H:24
Coord::y
@ y
Coord::x
@ x
pp
ParmParse pp("prob")
dx
const Real dx
Definition: ERF_InitCustomPert_ABL.H:23
pp_prob
ParmParse pp_prob("prob")
Ampl
Real Ampl
Definition: ERF_InitCustomPert_MovingTerrain.H:4
wavelength
Real wavelength
Definition: ERF_InitCustomPert_MovingTerrain.H:5
kp
Real kp
Definition: ERF_InitCustomPert_MovingTerrain.H:8
height
Real height
Definition: ERF_InitCustomPert_SquallLine.H:33
ParallelFor
ParallelFor(grown_box, [=] AMREX_GPU_DEVICE(int i, int j, int k) { qrcuten_arr(i, j, k)=Real(0);qscuten_arr(i, j, k)=Real(0);qicuten_arr(i, j, k)=Real(0);})
zbx
const Box zbx
Definition: ERF_SetupDiff.H:9
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
MicVar_Morr::omega
@ omega
Definition: ERF_Morrison.H:54
module_model_constants::g
real(c_double), parameter g
Definition: ERF_module_model_constants.F90:19
Here is the call graph for this function: