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 Real &time)
 

Function Documentation

◆ init_my_custom_terrain()

void init_my_custom_terrain ( const Geometry &  geom,
FArrayBox &  terrain_fab,
const Real &  time 
)
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  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int)
201  {
202  // Clip indices for ghost-cells
203  int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);
204 
205  // Location of nodes
206  Real x = ii * dx[0];
207 
208  // Wave height
209  Real height = Ampl * std::sin(kp * x - omega * time);
210 
211  // Populate terrain height
212  z_arr(i,j,0) = height;
213  });
214 
215  } else if (custom_terrain_type == "WindFarmTest") {
216 
217  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int)
218  {
219  // Clip indices for ghost-cells
220  int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);
221  int jj = amrex::min(amrex::max(j,domlo_y),domhi_y);
222 
223  // Location of nodes
224  Real x = (ProbLoArr[0] + ii * dx[0] - xcen);
225  Real y = (ProbLoArr[1] + jj * dx[1] - ycen);
226 
227  Real x_L = x/Real(100.0);
228  Real y_L = y/Real(100.0);
229 
230  z_arr(i,j,k0) = Real(100.0) / (one + x_L*x_L + y_L*y_L);
231  });
232 
233  } else if (custom_terrain_type == "RaisedFlat") {
234  Real z_offset = zero; pp_prob.query("z_offset" , z_offset);
235  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int)
236  {
237  z_arr(i,j,k0) = z_offset;
238  });
239  } else if (custom_terrain_type == "Cos4Hill") {
240 
241  // Get prob parameters (must be outside GPU kernel)
242  Real hm = zero; pp_prob.query("hmax", hm);
243  Real L = Real(100.0); pp_prob.query("L", L);
244  Real z_offset = zero; pp_prob.query("z_offset", z_offset);
245  Real fourL = four * L;
246 
247  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int)
248  {
249 
250  // Clip indices for ghost-cells
251  int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);
252  int jj = amrex::min(amrex::max(j,domlo_y),domhi_y);
253 
254  // Location of nodes
255  Real x = (ProbLoArr[0] + ii * dx[0] - xcen);
256  Real y = (ProbLoArr[1] + jj * dx[1] - ycen);
257  Real r = std::sqrt(x*x + y*y);
258 
259  if (r < fourL) {
260  z_arr(i,j,k0) = z_offset + hm * Real(0.0625) * std::pow(one + std::cos(PI*r/fourL), four);
261  } else {
262  z_arr(i,j,k0) = z_offset;
263  }
264  });
265 
266  } else if (custom_terrain_type == "None") {
267  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int)
268  {
269  z_arr(i,j,k0) = zero;
270  });
271  } else {
272  Abort("Don't know this custom_terrain_type");
273  }
274  }
275 }
four
constexpr amrex::Real four
Definition: ERF_Constants.H:10
two
constexpr amrex::Real two
Definition: ERF_Constants.H:8
one
constexpr amrex::Real one
Definition: ERF_Constants.H:7
zero
constexpr amrex::Real zero
Definition: ERF_Constants.H:6
myhalf
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:11
PI
constexpr amrex::Real PI
Definition: ERF_Constants.H:24
CONST_GRAV
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:40
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:53
module_model_constants::g
real(c_double), parameter g
Definition: ERF_module_model_constants.F90:19
Here is the call graph for this function: