ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
ERF_ProbCommon.H
Go to the documentation of this file.
1 #ifndef ERF_PROBCOMMON_H_
2 #define ERF_PROBCOMMON_H_
3 
4 #include <AMReX_ParmParse.H>
5 #include <AMReX_Geometry.H>
6 #include <AMReX_FArrayBox.H>
7 #include <AMReX_MultiFab.H>
8 
9 #include "ERF_DataStruct.H"
10 #include "ERF_EOS.H"
11 #include "ERF_HSEUtils.H"
12 #include "ERF_TileNoZ.H"
13 
15  amrex::Real rho_0 = 1.0;
16  amrex::Real T_0 = 300.0;
17 };
18 
19 template <typename T>
20 T erf_get_single_value (std::istream& is, int n)
21 {
22  T r = 0;
23  std::string line;
24  if (std::getline(is, line)) {
25  std::stringstream ss(line);
26  ss >> r;
27  if (ss.fail()) {
28  amrex::Abort("Failed to read");
29  }
30  if (ss >> r) {
31  amrex::Print() << "Trying to read line " << n << " in this file " << std::endl;
32  amrex::Abort("Wrong format: more than one number in this line");
33  }
34  } else {
35  amrex::Abort("Unable to read this line");
36  }
37  return r;
38 }
39 
40 /**
41  * Class to hold problem-specific routines
42 */
44 {
45 public:
46 
47  /**
48  * Virtual destructor to avoid data leakage with derived class
49  */
50  virtual ~ProblemBase () = default;
51 
52  /**
53  * Function to initialize the hydrostatic reference density
54  *
55  * @param[out] rho_hse hydrostatic reference density
56  * @param[in] z_phys_nd height coordinate at nodes
57  * @param[in] z_phys_cc height coordinate at cell centers
58  * @param[in] geom container for geometric information
59  */
60  virtual void
61  erf_init_dens_hse (amrex::MultiFab& /*rho_hse*/,
62  std::unique_ptr<amrex::MultiFab>& /*z_phys_nd*/,
63  std::unique_ptr<amrex::MultiFab>& /*z_phys_cc*/,
64  amrex::Geometry const& /*geom*/)
65  {
66  amrex::Print() << "Hydrostatically balanced density was NOT set"
67  << " -- an appropriate init_type should probably have been specified"
68  << " (e.g., input_sounding, WRFInput, or Metgrid)"
69  << std::endl;
70  amrex::Error("Should never call erf_init_dens_hse for "+name()+" problem");
71  }
72 
73  virtual void
74  erf_init_dens_hse_moist (amrex::MultiFab& /*rho_hse*/,
75  std::unique_ptr<amrex::MultiFab>& /*z_phys_nd*/,
76  amrex::Geometry const& /*geom*/)
77  {
78 
79  }
80 
81  /**
82  * Function to perform custom initialization of a test problem
83  *
84  * @param[in] bx cell-centered box on which to initialize scalars
85  * @param[in] xbx box on which to initialize x-component of velocity
86  * @param[in] ybx box on which to initialize y-component of velocity
87  * @param[in] zbx box on which to initialize z-component of velocity
88  * @param[out] state cell-centered variables to be filled in this routine
89  * @param[out] x_velocity x-component of velocity to be filled in this routine
90  * @param[out] y_velocity y-component of velocity to be filled in this routine
91  * @param[out] z_velocity z-component of velocity to be filled in this routine
92  * @param[out] r_hse hydrostatic reference density
93  * @param[out] p_hse hydrostatic reference pressure
94  * @param[in] z_nd height coordinate at nodes
95  * @param[in] z_cc height coordinate at cell centers
96  * @param[in] qv water vapor
97  * @param[in] qc cloud water
98  * @param[in] qi cloud ice
99  * @param[in] mf_m map factor on cell centers
100  * @param[in] mf_u map factor on x-faces
101  * @param[in] mf_v map factor on y-faces
102  * @param[in] sc SolverChoice structure that carries parameters
103  */
104  virtual void
105  init_custom_pert (const amrex::Box& /*bx*/,
106  const amrex::Box& /*xbx*/,
107  const amrex::Box& /*ybx*/,
108  const amrex::Box& /*zbx*/,
109  amrex::Array4<amrex::Real const> const& /*state*/,
110  amrex::Array4<amrex::Real > const& /*state_pert*/,
111  amrex::Array4<amrex::Real > const& /*x_vel_pert*/,
112  amrex::Array4<amrex::Real > const& /*y_vel_pert*/,
113  amrex::Array4<amrex::Real > const& /*z_vel_pert*/,
114  amrex::Array4<amrex::Real > const& /*r_hse*/,
115  amrex::Array4<amrex::Real > const& /*p_hse*/,
116  amrex::Array4<amrex::Real const> const& /*z_nd*/,
117  amrex::Array4<amrex::Real const> const& /*z_cc*/,
118  amrex::GeometryData const& /*geomdata*/,
119  amrex::Array4<amrex::Real const> const& /*mf_m*/,
120  amrex::Array4<amrex::Real const> const& /*mf_u*/,
121  amrex::Array4<amrex::Real const> const& /*mf_v*/,
122  const SolverChoice& /*sc*/)
123  {
124  amrex::Print() << "No perturbation to background fields supplied for "
125  << name() << " problem" << std::endl;
126  }
127 
128  /**
129  * Function to update user-specified temperature source terms that can
130  * vary with time and height.
131  *
132  * @param[in] time current time
133  * @param[out] rhotheta_source forcing profile
134  * @param[in] geom container for geometric information
135  * @param[in] z_phys_cc height coordinate at cell centers
136  */
137  virtual void
138  update_rhotheta_sources (const amrex::Real& /*time*/,
139  amrex::Vector<amrex::Real>& src,
140  amrex::Gpu::DeviceVector<amrex::Real>& d_src,
141  const amrex::Geometry& geom,
142  std::unique_ptr<amrex::MultiFab>& /*z_phys_cc*/)
143  {
144  if (src.empty()) return;
145 
146  amrex::Warning("Temperature forcing not defined for "+name()+" problem");
147 
148  const int khi = geom.Domain().bigEnd()[2];
149  // const amrex::Real* prob_lo = geom.ProbLo();
150  // const auto dx = geom.CellSize();
151  for (int k = 0; k <= khi; k++)
152  {
153  // const amrex::Real z_cc = prob_lo[2] + (k+0.5)* dx[2];
154  // set RHS term of RhoTheta equation based on time, z_cc here
155  src[k] = 0.0;
156  }
157 
158  // Copy from host version to device version
159  amrex::Gpu::copy(amrex::Gpu::hostToDevice, src.begin(), src.end(), d_src.begin());
160  }
161 
162  /**
163  * Function to update user-specified moisture source terms that can
164  * vary with time and height.
165  *
166  * @param[in] time current time
167  * @param[out] rhoqt_source moisture forcing profile
168  * @param[in] geom container for geometric information
169  * @param[in] z_phys_cc height coordinate at cell centers
170  */
171  virtual void
172  update_rhoqt_sources (const amrex::Real& /*time*/,
173  amrex::Vector<amrex::Real>& qsrc,
174  amrex::Gpu::DeviceVector<amrex::Real>& d_qsrc,
175  const amrex::Geometry& geom,
176  std::unique_ptr<amrex::MultiFab>& /*z_phys_cc*/)
177  {
178  if (qsrc.empty()) return;
179 
180  amrex::Warning("Moisture forcing not defined for "+name()+" problem");
181 
182  const int khi = geom.Domain().bigEnd()[2];
183  // const amrex::Real* prob_lo = geom.ProbLo();
184  // const auto dx = geom.CellSize();
185  for (int k = 0; k <= khi; k++)
186  {
187  // const amrex::Real z_cc = prob_lo[2] + (k+0.5)* dx[2];
188  // set RHS term of RhoQ1 equation based on time, z_cc here
189  qsrc[k] = 0.0;
190  }
191 
192  // Copy from host version to device version
193  amrex::Gpu::copy(amrex::Gpu::hostToDevice, qsrc.begin(), qsrc.end(), d_qsrc.begin());
194  }
195 
196  /**
197  * Function to update the vertical velocity profile, used to add subsidence
198  * source terms for x-mom, y-mom, rho*theta, rho*Q1, and rho*Q2.
199  *
200  * TODO: Currently, this is only called by InitData, so there is no time
201  * dependence.
202  *
203  * @param[in] time current time
204  * @param[out] wbar w vel forcing profile
205  * @param[in] geom container for geometric information
206  * @param[in] z_phys_cc height coordinate at cell centers
207  */
208  virtual void
209  update_w_subsidence (const amrex::Real& /*time*/,
210  amrex::Vector<amrex::Real>& wbar,
211  amrex::Gpu::DeviceVector<amrex::Real>& d_wbar,
212  const amrex::Geometry& geom,
213  std::unique_ptr<amrex::MultiFab>& /*z_phys_cc*/)
214  {
215  if (wbar.empty()) return;
216 
217  amrex::Warning("Moisture forcing not defined for "+name()+" problem");
218 
219  const int khi = geom.Domain().bigEnd()[2];
220  // const amrex::Real* prob_lo = geom.ProbLo();
221  // const auto dx = geom.CellSize();
222  for (int k = 0; k <= khi; k++)
223  {
224  // const amrex::Real z_cc = prob_lo[2] + (k+0.5)* dx[2];
225  // set vertical velocity profile based on time, z_cc here
226  wbar[k] = 0.0;
227  }
228 
229  // Copy from host version to device version
230  amrex::Gpu::copy(amrex::Gpu::hostToDevice, wbar.begin(), wbar.end(), d_wbar.begin());
231  }
232 
233  /**
234  * Function to update user-specified geostrophic wind profile.
235  *
236  * @param[in] time current time
237  * @param[out] u_geos geostrophic wind profile
238  * @param[out] v_geos geostrophic wind profile
239  * @param[in] geom container for geometric information
240  * @param[in] z_phys_cc height coordinate at cell centers
241  */
242  virtual void
243  update_geostrophic_profile (const amrex::Real& /*time*/,
244  amrex::Vector<amrex::Real>& u_geos,
245  amrex::Gpu::DeviceVector<amrex::Real>& d_u_geos,
246  amrex::Vector<amrex::Real>& v_geos,
247  amrex::Gpu::DeviceVector<amrex::Real>& d_v_geos,
248  const amrex::Geometry& geom,
249  std::unique_ptr<amrex::MultiFab>& /*z_phys_cc*/)
250  {
251  if (u_geos.empty()) return;
252 
253  amrex::Warning("Geostrophic wind profile not defined for "+name()+" problem");
254 
255  const int khi = geom.Domain().bigEnd()[2];
256  // const amrex::Real* prob_lo = geom.ProbLo();
257  // const auto dx = geom.CellSize();
258  for (int k = 0; k <= khi; k++)
259  {
260  // const amrex::Real z_cc = prob_lo[2] + (k+0.5)* dx[2];
261  // set RHS term of RhoTheta equation based on time, z_cc here
262  u_geos[k] = 0.0;
263  v_geos[k] = 0.0;
264  }
265 
266  // Copy from host version to device version
267  amrex::Gpu::copy(amrex::Gpu::hostToDevice, u_geos.begin(), u_geos.end(), d_u_geos.begin());
268  amrex::Gpu::copy(amrex::Gpu::hostToDevice, v_geos.begin(), v_geos.end(), d_v_geos.begin());
269  }
270 
271  /**
272  * Function to perform custom initialization of terrain
273  *
274  * @param[in] geom container for geometric information
275  * @param[out] z_phys_nd height coordinate at nodes
276  * @param[in] time current time
277  */
278  void
279  init_terrain_surface (const amrex::Geometry& geom,
280  amrex::FArrayBox& terrain_fab,
281  const amrex::Real& time)
282  {
283  // Check if a valid text file exists for the terrain
284  std::string fname, fname_usgs;
285  amrex::ParmParse pp("erf");
286  auto valid_fname = pp.query("terrain_file_name",fname);
287  auto valid_fname_USGS = pp.query("terrain_file_name_USGS",fname_usgs);
288 
289  bool is_usgs;
290 
291  if (valid_fname) {
292  is_usgs = false;
293  read_custom_terrain(fname,is_usgs,geom,terrain_fab,time);
294 
295  } else if (valid_fname_USGS) {
296  is_usgs = true;
297  read_custom_terrain(fname_usgs,is_usgs,geom,terrain_fab,time);
298 
299  } else {
300  init_custom_terrain (geom, terrain_fab, time);
301  }
302  }
303 
304  void
305  read_custom_terrain(const std::string& fname,
306  const bool is_usgs,
307  const amrex::Geometry& geom,
308  amrex::FArrayBox& terrain_fab,
309  const amrex::Real& /*time*/)
310  {
311  amrex::Vector<amrex::Real> m_xterrain,m_yterrain,m_zterrain;
312 
313  int nx = 0; int ny = 0;
314 
315  if (amrex::ParallelDescriptor::IOProcessor()) {
316  // Read terrain file
317  amrex::Print()<<"Reading terrain file: "<< fname<< std::endl;
318  std::ifstream file(fname);
319 
320  if (!file.is_open()) {
321  amrex::Abort("Error: Could not open the file " + fname+ "\n");
322  }
323 
324  // Check if file is empty
325  if (file.peek() == std::ifstream::traits_type::eof()) {
326  amrex::Abort("Error: The file " + fname+ " is empty.\n");
327  }
328 
329  amrex::Real value1,value2,value3;
330 
331  if (is_usgs) {
332  amrex::Real lat_min, lon_min;
333 
334  file >> lon_min >> lat_min;
335  if(std::fabs(lon_min) > 180.0) {
336  amrex::Error("The value of longitude for entry in the first line in " + fname
337  + " should not exceed 180. It is " + std::to_string(lon_min));
338  }
339  if(std::fabs(lat_min) > 90.0) {
340  amrex::Error("The value of latitude for entry in the first line in " + fname
341  + " should not exceed 90. It is " + std::to_string(lat_min));
342  }
343 
344  file >> nx >> ny;
345 
346  int counter = 0;
347  while (file >> value1 >> value2 >> value3) {
348  m_xterrain.push_back(value1);
349  if(counter%nx==0) {
350  m_yterrain.push_back(value2);
351  }
352  m_zterrain.push_back(value3);
353  counter += 1;
354  }
355  AMREX_ASSERT(m_xterrain.size() == static_cast<long int>(nx*ny));
356  AMREX_ASSERT(m_yterrain.size() == static_cast<long int>(ny));
357  AMREX_ASSERT(m_zterrain.size() == static_cast<long int>(nx*ny));
358 
359  } else {
360  int cnt = 1;
361  nx = erf_get_single_value<int>(file,cnt); cnt++;
362  ny = erf_get_single_value<int>(file,cnt); cnt++;
363  amrex::Print()<<"Expecting " << nx << " values of x, " <<
364  ny << " values of y, and " <<
365  nx*ny << " values of z" << std::endl;
366  AMREX_ALWAYS_ASSERT(nx > 0);
367  AMREX_ALWAYS_ASSERT(ny > 0);
368  m_xterrain.resize(nx);
369  m_yterrain.resize(ny);
370  m_zterrain.resize(nx * ny);
371  for (int n = 0; n < nx; n++) {
372  m_xterrain[n] = erf_get_single_value<amrex::Real>(file,cnt);
373  cnt++;
374  }
375  for (int n = 0; n < ny; n++) {
376  m_yterrain[n] = erf_get_single_value<amrex::Real>(file,cnt);
377  cnt++;
378  }
379  for (int n = 0; n < nx * ny; n++) {
380  m_zterrain[n] = erf_get_single_value<amrex::Real>(file,cnt);
381  cnt++;
382  }
383  }
384 
385  // Close the file after reading
386  file.close();
387  }
388 
389  amrex::ParallelDescriptor::Bcast(&nx,1,amrex::ParallelDescriptor::IOProcessorNumber());
390  amrex::ParallelDescriptor::Bcast(&ny,1,amrex::ParallelDescriptor::IOProcessorNumber());
391 
392  int nz = nx * ny;
393 
394  m_xterrain.resize(nx);
395  m_yterrain.resize(ny);
396  m_zterrain.resize(nz);
397 
398  amrex::ParallelDescriptor::Bcast(m_xterrain.data(),nx,amrex::ParallelDescriptor::IOProcessorNumber());
399  amrex::ParallelDescriptor::Bcast(m_yterrain.data(),ny,amrex::ParallelDescriptor::IOProcessorNumber());
400  amrex::ParallelDescriptor::Bcast(m_zterrain.data(),nz,amrex::ParallelDescriptor::IOProcessorNumber());
401 
402  // Copy data to the GPU
403  amrex::Gpu::DeviceVector<amrex::Real> d_xterrain(nx),d_yterrain(ny),d_zterrain(nz);
404  amrex::Gpu::copy(amrex::Gpu::hostToDevice, m_xterrain.begin(), m_xterrain.end(), d_xterrain.begin());
405  amrex::Gpu::copy(amrex::Gpu::hostToDevice, m_yterrain.begin(), m_yterrain.end(), d_yterrain.begin());
406  amrex::Gpu::copy(amrex::Gpu::hostToDevice, m_zterrain.begin(), m_zterrain.end(), d_zterrain.begin());
407 
408  amrex::Real* d_xt = d_xterrain.data();
409  amrex::Real* d_yt = d_yterrain.data();
410  amrex::Real* d_zt = d_zterrain.data();
411 
412  auto dx = geom.CellSizeArray();
413  auto ProbLoArr = geom.ProbLoArray();
414 
415  int ilo = geom.Domain().smallEnd(0);
416  int jlo = geom.Domain().smallEnd(1);
417  int klo = geom.Domain().smallEnd(2);
418  int ihi = geom.Domain().bigEnd(0) + 1;
419  int jhi = geom.Domain().bigEnd(1) + 1;
420 
421  amrex::Box zbx = terrain_fab.box();
422  amrex::Array4<amrex::Real> const& z_arr = terrain_fab.array();
423 
424  amrex::ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/)
425  {
426  // Clip indices for ghost cells
427  int ii = amrex::min(amrex::max(i,ilo),ihi);
428  int jj = amrex::min(amrex::max(j,jlo),jhi);
429 
430  // Location of nodes
431  amrex::Real x = ProbLoArr[0] + ii * dx[0] + 1e-3;
432  amrex::Real y = ProbLoArr[1] + jj * dx[1] + 1e-3;
433 
434  int ind11, ind12, ind21, ind22;
435  amrex::Real x1, x2, y1, y2;
436 
437  int iindex_terrain=-1;
438  int jindex_terrain=-1;
439  // *******************************************************************
440  // NOTE: usgs-format is contiguous in x
441  // *******************************************************************
442  if (is_usgs) {
443 
444  for (int it = 0; it < ny && jindex_terrain == -1; it++) {
445  if (d_yt[it] > y) {
446  jindex_terrain = it-1;
447  }
448  }
449  if (jindex_terrain == -1) {
450  jindex_terrain = ny-1;
451  }
452 
453  int gstart = (jindex_terrain )*nx;
454  int gend = (jindex_terrain+1)*nx-1;
455  for (int it = gstart; it <= gend && iindex_terrain == -1; it++) {
456  if (d_xt[it] > x) {
457  iindex_terrain = it-gstart;
458  }
459  }
460 
461  // Define the four values to interpolate between
462  ind11 = jindex_terrain*nx + iindex_terrain; // (x1,y1)
463  ind12 = ind11+nx; // (x1,y2)
464  ind21 = ind11+1; // (x2,y1)
465  ind22 = ind12+1; // (x2,y2)
466 
467  x1 = d_xt[ind11];
468  x2 = d_xt[ind21];
469  y1 = d_yt[jindex_terrain];
470  y2 = d_yt[jindex_terrain+1];
471 
472  amrex::Real denom = (x2-x1)*(y2-y1);
473  amrex::Real w_11 = (x2-x)*(y2-y)/denom; // (x1,y1)
474  amrex::Real w_12 = (x2-x)*(y-y1)/denom; // (x1,y2)
475  amrex::Real w_21 = (x-x1)*(y2-y)/denom; // (x2,y1)
476  amrex::Real w_22 = (x-x1)*(y-y1)/denom; // (x2,y2)
477 
478  // Do bilinear interpolation
479  z_arr(i,j,klo) = w_11*d_zt[ind11] + w_12*d_zt[ind12] + w_21*d_zt[ind21] + w_22*d_zt[ind22];
480 
481  } else {
482 
483  for (int it = 0; it < ny && jindex_terrain == -1; it++) {
484  if (d_yt[it] > y) {
485  jindex_terrain = it-1;
486  }
487  }
488  if (jindex_terrain == -1) {
489  jindex_terrain = ny-1;
490  }
491 
492  for (int it = 0; it < nx && iindex_terrain == -1; it++) {
493  if (d_xt[it] > x) {
494  iindex_terrain = it-1;
495  }
496  }
497  if (iindex_terrain == -1) {
498  iindex_terrain = nx-1;
499  }
500 
501  // Define the four values to interpolate between
502  x1 = d_xt[iindex_terrain];
503  x2 = d_xt[iindex_terrain+1];
504  y1 = d_yt[jindex_terrain];
505  y2 = d_yt[jindex_terrain+1];
506 
507 #if 1
508  // *******************************************************************
509  // NOTE: this format is contiguous in y to match the AMR-Wind read
510  // *******************************************************************
511  ind11 = iindex_terrain * ny + jindex_terrain; // (x1,y1)
512  ind21 = std::min(iindex_terrain+1,nx-1) * ny + jindex_terrain; // (x2,y1)
513 
514  ind12 = iindex_terrain * ny + std::min(jindex_terrain+1,ny-1); // (x1,y2)
515  ind22 = std::min(iindex_terrain+1,nx-1) * ny + std::min(jindex_terrain+1,ny-1); // (x1,y2)
516 #else
517  // *******************************************************************
518  // NOTE: this format is contiguous in x as an alternative
519  // *******************************************************************
520 
521  ind11 = jindex_terrain * nx + iindex_terrain; // (x1,y1)
522  ind12 = std::min(jindex_terrain+1,ny-1) * nx + iindex_terrain; // (x1,y2)
523 
524  ind21 = jindex_terrain * nx + std::min(iindex_terrain+1,nx-1); // (x2,y1)
525  ind22 = std::min(jindex_terrain+1,ny-1) * nx + std::min(iindex_terrain+1,nx-1); // (x2,y2)
526 #endif
527 
528  if (iindex_terrain == nx-1 && jindex_terrain == ny-1)
529  {
530  z_arr(i,j,klo) = d_zt[ind11];
531  }
532  else if (iindex_terrain != nx-1 && jindex_terrain == ny-1)
533  {
534  amrex::Real w_11 = (x2-x); // (x1,y1)
535  amrex::Real w_21 = (x-x1); // (x2,y1)
536  amrex::Real denom = (x2-x1);
537  z_arr(i,j,klo) = (w_11*d_zt[ind11] + w_21*d_zt[ind21])/denom;
538  }
539  else if (iindex_terrain == nx-1 && jindex_terrain != ny-1)
540  {
541  amrex::Real w_11 = (y2-y); // (x1,y1)
542  amrex::Real w_12 = (y-y1); // (x1,y2)
543  amrex::Real denom = (y2-y1);
544  z_arr(i,j,klo) = (w_11*d_zt[ind11] + w_12*d_zt[ind12])/denom;
545  }
546  else
547  {
548  amrex::Real w_11 = (x2-x)*(y2-y); // (x1,y1)
549  amrex::Real w_21 = (x-x1)*(y2-y); // (x2,y1)
550  amrex::Real w_12 = (x2-x)*(y-y1); // (x1,y2)
551  amrex::Real w_22 = (x-x1)*(y-y1); // (x2,y2)
552  amrex::Real denom = (x2-x1)*(y2-y1);
553  z_arr(i,j,klo) = (w_11*d_zt[ind11] + w_12*d_zt[ind12] + w_21*d_zt[ind21] + w_22*d_zt[ind22]) / denom;
554  }
555  } // usgs?
556  });
557  }
558 
559  /**
560  * Function to perform custom initialization of terrain
561  *
562  * Note: Terrain functionality can also be used to provide grid stretching.
563  *
564  * @param[in] geom container for geometric information
565  * @param[out] terrain_mf height coordinate at nodes
566  * @param[in] time current time
567  */
568  virtual void
569  init_custom_terrain (const amrex::Geometry& /*geom*/,
570  amrex::FArrayBox& terrain_fab,
571  const amrex::Real& /*time*/)
572  {
573  amrex::Print() << "Initializing flat terrain" << std::endl;
574 
575  if (SolverChoice::mesh_type == MeshType::VariableDz) {
576  SolverChoice::set_mesh_type(MeshType::StretchedDz);
577  amrex::Print() << "Resetting mesh type to StretchedDz" << std::endl;
578  }
579 
580  terrain_fab.template setVal<amrex::RunOn::Device>(0.0);
581  }
582 
583 #ifdef ERF_USE_TERRAIN_VELOCITY
584  virtual amrex::Real compute_terrain_velocity(const amrex::Real /*time*/)
585  {
586  amrex::Error("Should never call compute_terrain_velocity for "+name()+" problem");
587  }
588 #endif
589 
590  /**
591  * Function to define the quantities needed to impose Rayleigh damping
592  *
593  * @param[out] rayleigh_ptrs = {strength of Rayleigh damping, reference values for xvel/yvel/zvel/theta used to define Rayleigh damping}
594  * @param[in] geom container for geometric information
595  * @param[in] z_phys_cc height coordinate at cell centers
596  */
597  virtual void
598  erf_init_rayleigh (amrex::Vector<amrex::Vector<amrex::Real> >& /*rayleigh_ptrs*/,
599  amrex::Geometry const& /*geom*/,
600  std::unique_ptr<amrex::MultiFab>& /*z_phys_nd*/,
601  amrex::Real /*zdamp*/)
602  {
603  amrex::Error("Should never call erf_init_rayleigh for "+name()+" problem");
604  }
605 
606  /**
607  * Function to set uniform background density and potential temperature fields
608  */
609  void
610  init_uniform (const amrex::Box& bx, amrex::Array4<amrex::Real> const& state)
611  {
612  amrex::Real rho_0 = base_parms.rho_0;
613  amrex::Real T_0 = base_parms.T_0;
614  amrex::Print() << "Initializing uniform fields"
615  << " rho=" << rho_0 << " theta=" << T_0
616  << std::endl;
617 
618  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
619  {
620  state(i, j, k, Rho_comp) = rho_0;
621  state(i, j, k, RhoTheta_comp) = rho_0 * T_0;
622  });
623  }
624 
625 protected:
626  // Struct to store problem parameters
628 
629  /**
630  * Function to update default base parameters, currently only used for
631  * init_type == InitType::Uniform
632  */
633  void init_base_parms (amrex::Real rho_0, amrex::Real T_0) {
634  base_parms.rho_0 = rho_0;
635  base_parms.T_0 = T_0;
636  }
637 
638  // descriptor for problem definition
639  virtual std::string name () = 0;
640 };
641 
642 
643 /**
644  * Function to init the physical bounds of the domain
645  * and instantiate a Problem derived from ProblemBase
646 */
647 extern std::unique_ptr<ProblemBase> amrex_probinit (const amrex_real* problo,
648  const amrex_real* probhi) AMREX_ATTRIBUTE_WEAK;
649 #endif
@ wbar
Definition: ERF_DataStruct.H:91
const Box zbx
Definition: ERF_DiffSetup.H:23
#define Rho_comp
Definition: ERF_IndexDefines.H:36
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:37
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real pp(amrex::Real y)
Definition: ERF_MicrophysicsUtils.H:230
T erf_get_single_value(std::istream &is, int n)
Definition: ERF_ProbCommon.H:20
std::unique_ptr< ProblemBase > amrex_probinit(const amrex_real *problo, const amrex_real *probhi) AMREX_ATTRIBUTE_WEAK
Definition: ERF_ProbCommon.H:44
virtual void init_custom_pert(const amrex::Box &, const amrex::Box &, const amrex::Box &, const amrex::Box &, amrex::Array4< amrex::Real const > const &, amrex::Array4< amrex::Real > const &, amrex::Array4< amrex::Real > const &, amrex::Array4< amrex::Real > const &, amrex::Array4< amrex::Real > const &, amrex::Array4< amrex::Real > const &, amrex::Array4< amrex::Real > const &, amrex::Array4< amrex::Real const > const &, amrex::Array4< amrex::Real const > const &, amrex::GeometryData const &, amrex::Array4< amrex::Real const > const &, amrex::Array4< amrex::Real const > const &, amrex::Array4< amrex::Real const > const &, const SolverChoice &)
Definition: ERF_ProbCommon.H:105
virtual std::string name()=0
void init_uniform(const amrex::Box &bx, amrex::Array4< amrex::Real > const &state)
Definition: ERF_ProbCommon.H:610
void init_terrain_surface(const amrex::Geometry &geom, amrex::FArrayBox &terrain_fab, const amrex::Real &time)
Definition: ERF_ProbCommon.H:279
void read_custom_terrain(const std::string &fname, const bool is_usgs, const amrex::Geometry &geom, amrex::FArrayBox &terrain_fab, const amrex::Real &)
Definition: ERF_ProbCommon.H:305
virtual void update_w_subsidence(const amrex::Real &, amrex::Vector< amrex::Real > &wbar, amrex::Gpu::DeviceVector< amrex::Real > &d_wbar, const amrex::Geometry &geom, std::unique_ptr< amrex::MultiFab > &)
Definition: ERF_ProbCommon.H:209
virtual void erf_init_rayleigh(amrex::Vector< amrex::Vector< amrex::Real > > &, amrex::Geometry const &, std::unique_ptr< amrex::MultiFab > &, amrex::Real)
Definition: ERF_ProbCommon.H:598
void init_base_parms(amrex::Real rho_0, amrex::Real T_0)
Definition: ERF_ProbCommon.H:633
virtual void update_rhotheta_sources(const amrex::Real &, amrex::Vector< amrex::Real > &src, amrex::Gpu::DeviceVector< amrex::Real > &d_src, const amrex::Geometry &geom, std::unique_ptr< amrex::MultiFab > &)
Definition: ERF_ProbCommon.H:138
virtual void init_custom_terrain(const amrex::Geometry &, amrex::FArrayBox &terrain_fab, const amrex::Real &)
Definition: ERF_ProbCommon.H:569
virtual ~ProblemBase()=default
virtual void erf_init_dens_hse_moist(amrex::MultiFab &, std::unique_ptr< amrex::MultiFab > &, amrex::Geometry const &)
Definition: ERF_ProbCommon.H:74
ProbParmDefaults base_parms
Definition: ERF_ProbCommon.H:627
virtual void erf_init_dens_hse(amrex::MultiFab &, std::unique_ptr< amrex::MultiFab > &, std::unique_ptr< amrex::MultiFab > &, amrex::Geometry const &)
Definition: ERF_ProbCommon.H:61
virtual void update_geostrophic_profile(const amrex::Real &, amrex::Vector< amrex::Real > &u_geos, amrex::Gpu::DeviceVector< amrex::Real > &d_u_geos, amrex::Vector< amrex::Real > &v_geos, amrex::Gpu::DeviceVector< amrex::Real > &d_v_geos, const amrex::Geometry &geom, std::unique_ptr< amrex::MultiFab > &)
Definition: ERF_ProbCommon.H:243
virtual void update_rhoqt_sources(const amrex::Real &, amrex::Vector< amrex::Real > &qsrc, amrex::Gpu::DeviceVector< amrex::Real > &d_qsrc, const amrex::Geometry &geom, std::unique_ptr< amrex::MultiFab > &)
Definition: ERF_ProbCommon.H:172
@ T
Definition: ERF_IndexDefines.H:110
Definition: ERF_ProbCommon.H:14
amrex::Real T_0
Definition: ERF_ProbCommon.H:16
amrex::Real rho_0
Definition: ERF_ProbCommon.H:15
Definition: ERF_DataStruct.H:123
static MeshType mesh_type
Definition: ERF_DataStruct.H:755
static void set_mesh_type(MeshType new_mesh_type)
Definition: ERF_DataStruct.H:758