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.queryAdd("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.queryAdd("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.queryAdd("dir", dir);
69 
70  Real L = Real(100.0); pp_prob.queryAdd("L" , L);
71  Real z_offset = zero; pp_prob.queryAdd("z_offset" , z_offset);
72 
73  // If hm is nonzero, then use alternate hill definition
74  Real hm = zero; pp_prob.queryAdd("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 a_mtn = Real(5000.0); pp_prob.queryAdd("a" , a_mtn);
138  Real Hm = Real(250.0); pp_prob.queryAdd("hm" , Hm);
139  Real lambda = Real(4000.0); pp_prob.queryAdd("lambda", lambda);
140  Real asq = a_mtn * a_mtn;
141 
142  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int)
143  {
144  // Clip indices for ghost-cells
145  int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);
146 
147  // Location of nodes
148  Real x = (ProbLoArr[0] + ii * dx[0] - xcen);
149 
150  Real cosx = std::cos(PI * x / lambda);
151 
152  z_arr(i,j,k0) = Hm * std::exp(-x*x/asq) * cosx * cosx;
153  });
154 
155  } else if (custom_terrain_type == "HalfCylinder") {
156 
157  Real asq = myhalf * myhalf;
158 
159  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int)
160  {
161  // Clip indices for ghost-cells
162  int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);
163 
164  // Location of nodes
165  Real x = (ProbLoArr[0] + ii * dx[0] - xcen);
166 
167  Real rsq = x*x;
168 
169  if (rsq < asq) {
170  z_arr(i,j,k0) = std::sqrt(asq - rsq);
171  } else {
172  z_arr(i,j,k0) = zero;
173  }
174  });
175 
176  } else if (custom_terrain_type == "Hemisphere") {
177 
178  Real asq = myhalf * myhalf;
179 
180  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int)
181  {
182  // Clip indices for ghost-cells
183  int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);
184  int jj = amrex::min(amrex::max(j,domlo_y),domhi_y);
185 
186  // Location of nodes
187  Real x = (ProbLoArr[0] + ii * dx[0] - xcen);
188  Real y = (ProbLoArr[1] + jj * dx[1] - ycen);
189 
190  Real rsq = x*x + y*y;
191 
192  if (rsq < asq) {
193  z_arr(i,j,k0) = std::pow(asq-rsq, myhalf);
194  } else {
195  z_arr(i,j,k0) = zero;
196  }
197  });
198 
199  } else if (custom_terrain_type == "MovingSineWave") {
200 
201  Real Ampl = zero; pp_prob.queryAdd("Ampl", Ampl);
202  Real wavelength = Real(100.); pp_prob.queryAdd("wavelength", wavelength);
203 
204  Real kp = two * PI / wavelength;
205  Real g = CONST_GRAV;
206  Real omega = std::sqrt(g * kp);
207 
208  Real time = static_cast<Real>(time_d);
209 
210  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int)
211  {
212  // Clip indices for ghost-cells
213  int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);
214 
215  // Location of nodes
216  Real x = ii * dx[0];
217 
218  // Wave height
219  Real height = Ampl * std::sin(kp * x - omega * time);
220 
221  // Populate terrain height
222  z_arr(i,j,0) = height;
223  });
224 
225  } else if (custom_terrain_type == "WindFarmTest") {
226 
227  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int)
228  {
229  // Clip indices for ghost-cells
230  int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);
231  int jj = amrex::min(amrex::max(j,domlo_y),domhi_y);
232 
233  // Location of nodes
234  Real x = (ProbLoArr[0] + ii * dx[0] - xcen);
235  Real y = (ProbLoArr[1] + jj * dx[1] - ycen);
236 
237  Real x_L = x/Real(100.0);
238  Real y_L = y/Real(100.0);
239 
240  z_arr(i,j,k0) = Real(100.0) / (one + x_L*x_L + y_L*y_L);
241  });
242 
243  } else if (custom_terrain_type == "RaisedFlat") {
244  Real z_offset = zero; pp_prob.queryAdd("z_offset" , z_offset);
245  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int)
246  {
247  z_arr(i,j,k0) = z_offset;
248  });
249  } else if (custom_terrain_type == "Cos4Hill") {
250 
251  // Get prob parameters (must be outside GPU kernel)
252  Real hm = zero; pp_prob.queryAdd("hmax", hm);
253  Real L = Real(100.0); pp_prob.queryAdd("L", L);
254  Real z_offset = zero; pp_prob.queryAdd("z_offset", z_offset);
255  Real fourL = four * L;
256 
257  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int)
258  {
259 
260  // Clip indices for ghost-cells
261  int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);
262  int jj = amrex::min(amrex::max(j,domlo_y),domhi_y);
263 
264  // Location of nodes
265  Real x = (ProbLoArr[0] + ii * dx[0] - xcen);
266  Real y = (ProbLoArr[1] + jj * dx[1] - ycen);
267  Real r = std::sqrt(x*x + y*y);
268 
269  if (r < fourL) {
270  z_arr(i,j,k0) = z_offset + hm * Real(0.0625) * std::pow(one + std::cos(PI*r/fourL), four);
271  } else {
272  z_arr(i,j,k0) = z_offset;
273  }
274  });
275 
276  } else if (custom_terrain_type == "None") {
277  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int)
278  {
279  z_arr(i,j,k0) = zero;
280  });
281  } else {
282  Abort("Don't know this custom_terrain_type");
283  }
284  }
285 }
CONST_GRAV
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:68
num
@ num
Definition: ERF_DataStruct.H:37
Coord::y
@ y
Coord::x
@ x
pp_prob
ParmParse pp_prob("prob")
dx
const Real dx
Definition: ERF_InitCustomPert_ABL.H:44
pp
ParmParse pp("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);})
four
constexpr amrex::Real four
Definition: ERF_NumericalConstants.H:40
two
constexpr amrex::Real two
Definition: ERF_NumericalConstants.H:38
one
constexpr amrex::Real one
Definition: ERF_NumericalConstants.H:37
zero
constexpr amrex::Real zero
Definition: ERF_NumericalConstants.H:36
myhalf
constexpr amrex::Real myhalf
Definition: ERF_NumericalConstants.H:41
PI
constexpr amrex::Real PI
Definition: ERF_NumericalConstants.H:46
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
MicVar_Morr::omega
@ omega
Definition: ERF_Morrison.H:55
module_model_constants::g
real(c_double), parameter g
Definition: ERF_module_model_constants.F90:19
Here is the call graph for this function: