1 #ifndef ERF_PROBCOMMON_H_
2 #define ERF_PROBCOMMON_H_
4 #include <AMReX_ParmParse.H>
5 #include <AMReX_Geometry.H>
6 #include <AMReX_FArrayBox.H>
7 #include <AMReX_MultiFab.H>
16 amrex::Real
T_0 = 300.0;
24 if (std::getline(is, line)) {
25 std::stringstream ss(line);
28 amrex::Abort(
"Failed to read");
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");
35 amrex::Abort(
"Unable to read this line");
62 std::unique_ptr<amrex::MultiFab>& ,
63 std::unique_ptr<amrex::MultiFab>& ,
64 amrex::Geometry
const& )
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)"
70 amrex::Error(
"Should never call erf_init_dens_hse for "+
name()+
" problem");
75 std::unique_ptr<amrex::MultiFab>& ,
76 amrex::Geometry
const& )
109 amrex::Array4<amrex::Real const>
const& ,
110 amrex::Array4<amrex::Real >
const& ,
111 amrex::Array4<amrex::Real >
const& ,
112 amrex::Array4<amrex::Real >
const& ,
113 amrex::Array4<amrex::Real >
const& ,
114 amrex::Array4<amrex::Real >
const& ,
115 amrex::Array4<amrex::Real >
const& ,
116 amrex::Array4<amrex::Real const>
const& ,
117 amrex::Array4<amrex::Real const>
const& ,
118 amrex::GeometryData
const& ,
119 amrex::Array4<amrex::Real const>
const& ,
120 amrex::Array4<amrex::Real const>
const& ,
121 amrex::Array4<amrex::Real const>
const& ,
124 amrex::Print() <<
"No perturbation to background fields supplied for "
125 <<
name() <<
" problem" << std::endl;
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>& )
144 if (src.empty())
return;
146 amrex::Warning(
"Temperature forcing not defined for "+
name()+
" problem");
148 const int khi = geom.Domain().bigEnd()[2];
151 for (
int k = 0; k <= khi; k++)
159 amrex::Gpu::copy(amrex::Gpu::hostToDevice, src.begin(), src.end(), d_src.begin());
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>& )
178 if (qsrc.empty())
return;
180 amrex::Warning(
"Moisture forcing not defined for "+
name()+
" problem");
182 const int khi = geom.Domain().bigEnd()[2];
185 for (
int k = 0; k <= khi; k++)
193 amrex::Gpu::copy(amrex::Gpu::hostToDevice, qsrc.begin(), qsrc.end(), d_qsrc.begin());
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>& )
215 if (
wbar.empty())
return;
217 amrex::Warning(
"Moisture forcing not defined for "+
name()+
" problem");
219 const int khi = geom.Domain().bigEnd()[2];
222 for (
int k = 0; k <= khi; k++)
230 amrex::Gpu::copy(amrex::Gpu::hostToDevice,
wbar.begin(),
wbar.end(), d_wbar.begin());
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>& )
251 if (u_geos.empty())
return;
253 amrex::Warning(
"Geostrophic wind profile not defined for "+
name()+
" problem");
255 const int khi = geom.Domain().bigEnd()[2];
258 for (
int k = 0; k <= khi; k++)
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());
280 amrex::FArrayBox& terrain_fab,
281 const amrex::Real& time)
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);
295 }
else if (valid_fname_USGS) {
307 const amrex::Geometry& geom,
308 amrex::FArrayBox& terrain_fab,
311 amrex::Vector<amrex::Real> m_xterrain,m_yterrain,m_zterrain;
313 int nx = 0;
int ny = 0;
315 if (amrex::ParallelDescriptor::IOProcessor()) {
317 amrex::Print()<<
"Reading terrain file: "<< fname<< std::endl;
318 std::ifstream file(fname);
320 if (!file.is_open()) {
321 amrex::Abort(
"Error: Could not open the file " + fname+
"\n");
325 if (file.peek() == std::ifstream::traits_type::eof()) {
326 amrex::Abort(
"Error: The file " + fname+
" is empty.\n");
329 amrex::Real value1,value2,value3;
332 amrex::Real lat_min, lon_min;
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));
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));
347 while (file >> value1 >> value2 >> value3) {
348 m_xterrain.push_back(value1);
350 m_yterrain.push_back(value2);
352 m_zterrain.push_back(value3);
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));
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);
375 for (
int n = 0; n < ny; n++) {
376 m_yterrain[n] = erf_get_single_value<amrex::Real>(file,cnt);
379 for (
int n = 0; n < nx * ny; n++) {
380 m_zterrain[n] = erf_get_single_value<amrex::Real>(file,cnt);
389 amrex::ParallelDescriptor::Bcast(&nx,1,amrex::ParallelDescriptor::IOProcessorNumber());
390 amrex::ParallelDescriptor::Bcast(&ny,1,amrex::ParallelDescriptor::IOProcessorNumber());
394 m_xterrain.resize(nx);
395 m_yterrain.resize(ny);
396 m_zterrain.resize(nz);
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());
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());
408 amrex::Real* d_xt = d_xterrain.data();
409 amrex::Real* d_yt = d_yterrain.data();
410 amrex::Real* d_zt = d_zterrain.data();
412 auto dx = geom.CellSizeArray();
413 auto ProbLoArr = geom.ProbLoArray();
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;
421 amrex::Box
zbx = terrain_fab.box();
422 amrex::Array4<amrex::Real>
const& z_arr = terrain_fab.array();
424 amrex::ParallelFor(
zbx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int )
427 int ii = amrex::min(amrex::max(i,ilo),ihi);
428 int jj = amrex::min(amrex::max(j,jlo),jhi);
431 amrex::Real
x = ProbLoArr[0] + ii * dx[0] + 1e-3;
432 amrex::Real
y = ProbLoArr[1] + jj * dx[1] + 1e-3;
434 int ind11, ind12, ind21, ind22;
435 amrex::Real x1, x2, y1, y2;
437 int iindex_terrain=-1;
438 int jindex_terrain=-1;
444 for (
int it = 0; it < ny && jindex_terrain == -1; it++) {
446 jindex_terrain = it-1;
449 if (jindex_terrain == -1) {
450 jindex_terrain = ny-1;
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++) {
457 iindex_terrain = it-gstart;
462 ind11 = jindex_terrain*nx + iindex_terrain;
469 y1 = d_yt[jindex_terrain];
470 y2 = d_yt[jindex_terrain+1];
472 amrex::Real denom = (x2-x1)*(y2-y1);
473 amrex::Real w_11 = (x2-
x)*(y2-
y)/denom;
474 amrex::Real w_12 = (x2-
x)*(
y-y1)/denom;
475 amrex::Real w_21 = (
x-x1)*(y2-
y)/denom;
476 amrex::Real w_22 = (
x-x1)*(
y-y1)/denom;
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];
483 for (
int it = 0; it < ny && jindex_terrain == -1; it++) {
485 jindex_terrain = it-1;
488 if (jindex_terrain == -1) {
489 jindex_terrain = ny-1;
492 for (
int it = 0; it < nx && iindex_terrain == -1; it++) {
494 iindex_terrain = it-1;
497 if (iindex_terrain == -1) {
498 iindex_terrain = nx-1;
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];
511 ind11 = iindex_terrain * ny + jindex_terrain;
512 ind21 = std::min(iindex_terrain+1,nx-1) * ny + jindex_terrain;
514 ind12 = iindex_terrain * ny + std::min(jindex_terrain+1,ny-1);
515 ind22 = std::min(iindex_terrain+1,nx-1) * ny + std::min(jindex_terrain+1,ny-1);
521 ind11 = jindex_terrain * nx + iindex_terrain;
522 ind12 = std::min(jindex_terrain+1,ny-1) * nx + iindex_terrain;
524 ind21 = jindex_terrain * nx + std::min(iindex_terrain+1,nx-1);
525 ind22 = std::min(jindex_terrain+1,ny-1) * nx + std::min(iindex_terrain+1,nx-1);
528 if (iindex_terrain == nx-1 && jindex_terrain == ny-1)
530 z_arr(i,j,klo) = d_zt[ind11];
532 else if (iindex_terrain != nx-1 && jindex_terrain == ny-1)
534 amrex::Real w_11 = (x2-
x);
535 amrex::Real w_21 = (
x-x1);
536 amrex::Real denom = (x2-x1);
537 z_arr(i,j,klo) = (w_11*d_zt[ind11] + w_21*d_zt[ind21])/denom;
539 else if (iindex_terrain == nx-1 && jindex_terrain != ny-1)
541 amrex::Real w_11 = (y2-
y);
542 amrex::Real w_12 = (
y-y1);
543 amrex::Real denom = (y2-y1);
544 z_arr(i,j,klo) = (w_11*d_zt[ind11] + w_12*d_zt[ind12])/denom;
548 amrex::Real w_11 = (x2-
x)*(y2-
y);
549 amrex::Real w_21 = (
x-x1)*(y2-
y);
550 amrex::Real w_12 = (x2-
x)*(
y-y1);
551 amrex::Real w_22 = (
x-x1)*(
y-y1);
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;
570 amrex::FArrayBox& terrain_fab,
573 amrex::Print() <<
"Initializing flat terrain" << std::endl;
577 amrex::Print() <<
"Resetting mesh type to StretchedDz" << std::endl;
580 terrain_fab.template setVal<amrex::RunOn::Device>(0.0);
583 #ifdef ERF_USE_TERRAIN_VELOCITY
584 virtual amrex::Real compute_terrain_velocity(
const amrex::Real )
586 amrex::Error(
"Should never call compute_terrain_velocity for "+
name()+
" problem");
599 amrex::Geometry
const& ,
600 std::unique_ptr<amrex::MultiFab>& ,
603 amrex::Error(
"Should never call erf_init_rayleigh for "+
name()+
" problem");
610 init_uniform (
const amrex::Box& bx, amrex::Array4<amrex::Real>
const& state)
614 amrex::Print() <<
"Initializing uniform fields"
615 <<
" rho=" << rho_0 <<
" theta=" << T_0
618 ParallelFor(bx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept
639 virtual std::string
name () = 0;
648 const amrex_real* probhi) AMREX_ATTRIBUTE_WEAK;
@ 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