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 = 0.5;
27  } else {
28  mf_m = 1.;
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 = 0.5;
43  Real num = 8. * a * a * a;
44  Real xcen = 0.5 * (ProbLoArr[0] + ProbHiArr[0]) / mf_m;
45  Real ycen = 0.5 * (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 = 100.0; pp_prob.query("L" , L);
64  Real z_offset = 0.0; pp_prob.query("z_offset" , z_offset);
65 
66  // If hm is nonzero, then use alternate hill definition
67  Real hm = 0.0; 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 + 4.0 * a * a);
99  } else {
100  Real y_L = y / L;
101  z_arr(i,j,k0) = hm / (1.0 + 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 + 4.0 * a * a);
119  } else {
120  Real r_L = r / L;
121  z_arr(i,j,k0) = hm / (1.0 + 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 = 5000.0 * 5000.0;
131  Real Hm = 250.0;
132  Real lambda = 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 = 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 = 0.5 * 0.5;
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) = pow(asq - rsq, 0.5);
163  } else {
164  z_arr(i,j,k0) = 0.0;
165  }
166  });
167 
168  } else if (custom_terrain_type == "Hemisphere") {
169 
170  Real asq = 0.5 * 0.5;
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, 0.5);
186  } else {
187  z_arr(i,j,k0) = 0.0;
188  }
189  });
190 
191  } else if (custom_terrain_type == "MovingSineWave") {
192 
193  Real Ampl = 0.0; pp_prob.query("Ampl", Ampl);
194  Real wavelength = 100.; pp_prob.query("wavelength", wavelength);
195 
196  Real kp = 2.0 * 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/100.0;
228  Real y_L = y/100.0;
229 
230  z_arr(i,j,k0) = 100.0 / (1.0 + x_L*x_L + y_L*y_L);
231  });
232 
233  } else if (custom_terrain_type == "RaisedFlat") {
234  Real z_offset = 0.0; 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 
240  } else if (custom_terrain_type == "None") {
241  ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int)
242  {
243  z_arr(i,j,k0) = 0.0;
244  });
245  } else {
246  Abort("Don't know this custom_terrain_type");
247  }
248  }
249 }
PI
constexpr amrex::Real PI
Definition: ERF_Constants.H:6
CONST_GRAV
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:21
num
@ num
Definition: ERF_DataStruct.H:24
Coord::y
@ y
Coord::x
@ x
dx
const auto dx
Definition: ERF_InitCustomPertVels_ParticleTests.H:15
pp
ParmParse pp("prob")
ParallelFor
ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { const auto prob_lo=geomdata.ProbLo();const auto dx=geomdata.CellSize();const Real x=(prob_lo[0]+(i+0.5) *dx[0])/mf_m(i, j, 0);const Real z=z_cc(i, j, k);Real L=std::sqrt(std::pow((x - x_c)/x_r, 2)+std::pow((z - z_c)/z_r, 2));if(L<=1.0) { Real dT=T_pert *(std::cos(PI *L)+1.0)/2.0;Real Tbar_hse=p_hse(i, j, k)/(R_d *r_hse(i, j, k));Real theta_perturbed=(Tbar_hse+dT) *std::pow(p_0/p_hse(i, j, k), rdOcp);Real theta_0=(Tbar_hse) *std::pow(p_0/p_hse(i, j, k), rdOcp);if(const_rho) { state_pert(i, j, k, RhoTheta_comp)=r_hse(i, j, k) *(theta_perturbed - theta_0);} else { state_pert(i, j, k, Rho_comp)=getRhoThetagivenP(p_hse(i, j, k))/theta_perturbed - r_hse(i, j, k);} } })
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: