377 amrex::Vector<amrex::Real> m_xterrain,m_yterrain,m_zterrain;
379 int nx = 0;
int ny = 0;
381 if (amrex::ParallelDescriptor::IOProcessor()) {
383 amrex::Print()<<
"Reading terrain file: "<< fname<< std::endl;
384 std::ifstream file(fname);
386 if (!file.is_open()) {
387 amrex::Abort(
"Error: Could not open the file " + fname+
"\n");
391 if (file.peek() == std::ifstream::traits_type::eof()) {
392 amrex::Abort(
"Error: The file " + fname+
" is empty.\n");
400 file >> lon_min >> lat_min;
401 if(std::fabs(lon_min) > 180.0) {
402 amrex::Error(
"The value of longitude for entry in the first line in " + fname
403 +
" should not exceed 180. It is " + std::to_string(lon_min));
405 if(std::fabs(lat_min) > 90.0) {
406 amrex::Error(
"The value of latitude for entry in the first line in " + fname
407 +
" should not exceed 90. It is " + std::to_string(lat_min));
413 while (file >> value1 >> value2 >> value3) {
414 m_xterrain.push_back(value1);
416 m_yterrain.push_back(value2);
418 m_zterrain.push_back(value3);
421 AMREX_ASSERT(m_xterrain.size() ==
static_cast<long int>(nx*ny));
422 AMREX_ASSERT(m_yterrain.size() ==
static_cast<long int>(ny));
423 AMREX_ASSERT(m_zterrain.size() ==
static_cast<long int>(nx*ny));
427 nx = erf_get_single_value<int>(file,cnt); cnt++;
428 ny = erf_get_single_value<int>(file,cnt); cnt++;
429 amrex::Print()<<
"Expecting " << nx <<
" values of x, " <<
430 ny <<
" values of y, and " <<
431 nx*ny <<
" values of z" << std::endl;
434 m_xterrain.resize(nx);
435 m_yterrain.resize(ny);
436 m_zterrain.resize(nx * ny);
437 for (
int n = 0; n < nx; n++) {
438 m_xterrain[n] = erf_get_single_value<amrex::Real>(file,cnt);
441 for (
int n = 0; n < ny; n++) {
442 m_yterrain[n] = erf_get_single_value<amrex::Real>(file,cnt);
445 for (
int n = 0; n < nx * ny; n++) {
446 m_zterrain[n] = erf_get_single_value<amrex::Real>(file,cnt);
455 amrex::ParallelDescriptor::Bcast(&nx,1,amrex::ParallelDescriptor::IOProcessorNumber());
456 amrex::ParallelDescriptor::Bcast(&ny,1,amrex::ParallelDescriptor::IOProcessorNumber());
460 m_xterrain.resize(nx);
461 m_yterrain.resize(ny);
462 m_zterrain.resize(nz);
464 amrex::ParallelDescriptor::Bcast(m_xterrain.data(),nx,amrex::ParallelDescriptor::IOProcessorNumber());
465 amrex::ParallelDescriptor::Bcast(m_yterrain.data(),ny,amrex::ParallelDescriptor::IOProcessorNumber());
466 amrex::ParallelDescriptor::Bcast(m_zterrain.data(),nz,amrex::ParallelDescriptor::IOProcessorNumber());
469 amrex::Gpu::DeviceVector<amrex::Real> d_xterrain(nx),d_yterrain(ny),d_zterrain(nz);
470 amrex::Gpu::copy(amrex::Gpu::hostToDevice, m_xterrain.begin(), m_xterrain.end(), d_xterrain.begin());
471 amrex::Gpu::copy(amrex::Gpu::hostToDevice, m_yterrain.begin(), m_yterrain.end(), d_yterrain.begin());
472 amrex::Gpu::copy(amrex::Gpu::hostToDevice, m_zterrain.begin(), m_zterrain.end(), d_zterrain.begin());
478 auto dx = geom.CellSizeArray();
479 auto ProbLoArr = geom.ProbLoArray();
481 int ilo = geom.Domain().smallEnd(0);
482 int jlo = geom.Domain().smallEnd(1);
483 int klo = geom.Domain().smallEnd(2);
484 int ihi = geom.Domain().bigEnd(0) + 1;
485 int jhi = geom.Domain().bigEnd(1) + 1;
487 amrex::Box
zbx = terrain_fab.box();
488 amrex::Array4<amrex::Real>
const& z_arr = terrain_fab.array();
493 int ii = amrex::min(amrex::max(i,ilo),ihi);
494 int jj = amrex::min(amrex::max(j,jlo),jhi);
500 int ind11, ind12, ind21, ind22;
503 int iindex_terrain=-1;
504 int jindex_terrain=-1;
510 for (
int it = 0; it < ny && jindex_terrain == -1; it++) {
512 jindex_terrain = it-1;
515 if (jindex_terrain == -1) {
516 jindex_terrain = ny-1;
519 int gstart = (jindex_terrain )*nx;
520 int gend = (jindex_terrain+1)*nx-1;
521 for (
int it = gstart; it <= gend && iindex_terrain == -1; it++) {
523 iindex_terrain = it-gstart;
528 ind11 = jindex_terrain*nx + iindex_terrain;
535 y1 = d_yt[jindex_terrain];
536 y2 = d_yt[jindex_terrain+1];
545 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];
549 for (
int it = 0; it < ny && jindex_terrain == -1; it++) {
551 jindex_terrain = it-1;
554 if (jindex_terrain == -1) {
555 jindex_terrain = ny-1;
558 for (
int it = 0; it < nx && iindex_terrain == -1; it++) {
560 iindex_terrain = it-1;
563 if (iindex_terrain == -1) {
564 iindex_terrain = nx-1;
568 x1 = d_xt[iindex_terrain];
569 x2 = d_xt[iindex_terrain+1];
570 y1 = d_yt[jindex_terrain];
571 y2 = d_yt[jindex_terrain+1];
577 ind11 = iindex_terrain * ny + jindex_terrain;
578 ind21 = std::min(iindex_terrain+1,nx-1) * ny + jindex_terrain;
580 ind12 = iindex_terrain * ny + std::min(jindex_terrain+1,ny-1);
581 ind22 = std::min(iindex_terrain+1,nx-1) * ny + std::min(jindex_terrain+1,ny-1);
587 ind11 = jindex_terrain * nx + iindex_terrain;
588 ind12 = std::min(jindex_terrain+1,ny-1) * nx + iindex_terrain;
590 ind21 = jindex_terrain * nx + std::min(iindex_terrain+1,nx-1);
591 ind22 = std::min(jindex_terrain+1,ny-1) * nx + std::min(iindex_terrain+1,nx-1);
594 if (iindex_terrain == nx-1 && jindex_terrain == ny-1)
596 z_arr(i,j,
klo) = d_zt[ind11];
598 else if (iindex_terrain != nx-1 && jindex_terrain == ny-1)
603 z_arr(i,j,
klo) = (w_11*d_zt[ind11] + w_21*d_zt[ind21])/denom;
605 else if (iindex_terrain == nx-1 && jindex_terrain != ny-1)
610 z_arr(i,j,
klo) = (w_11*d_zt[ind11] + w_12*d_zt[ind12])/denom;
619 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;
const int klo
Definition: ERF_InitCustomPertVels_ParticleTests.H:4
const Real dx
Definition: ERF_InitCustomPert_ABL.H:23
AMREX_ALWAYS_ASSERT(bx.length()[2]==khi+1)
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);} } })
const Box zbx
Definition: ERF_SetupDiff.H:9
amrex::Real Real
Definition: ERF_ShocInterface.H:19