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 
)

Initialize a problem-specific terrain profile on the nodal terrain FAB.

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