368 amrex::Vector<amrex::Real> m_xterrain,m_yterrain,m_zterrain;
370 int nx = 0;
int ny = 0;
372 if (amrex::ParallelDescriptor::IOProcessor()) {
374 amrex::Print()<<
"Reading terrain file: "<< fname<< std::endl;
375 std::ifstream file(fname);
377 if (!file.is_open()) {
378 amrex::Abort(
"Error: Could not open the file " + fname+
"\n");
382 if (file.peek() == std::ifstream::traits_type::eof()) {
383 amrex::Abort(
"Error: The file " + fname+
" is empty.\n");
391 file >> lon_min >> lat_min;
392 if(std::fabs(lon_min) > 180.0) {
393 amrex::Error(
"The value of longitude for entry in the first line in " + fname
394 +
" should not exceed 180. It is " + std::to_string(lon_min));
396 if(std::fabs(lat_min) > 90.0) {
397 amrex::Error(
"The value of latitude for entry in the first line in " + fname
398 +
" should not exceed 90. It is " + std::to_string(lat_min));
404 while (file >> value1 >> value2 >> value3) {
405 m_xterrain.push_back(value1);
407 m_yterrain.push_back(value2);
409 m_zterrain.push_back(value3);
412 AMREX_ASSERT(m_xterrain.size() ==
static_cast<long int>(nx*ny));
413 AMREX_ASSERT(m_yterrain.size() ==
static_cast<long int>(ny));
414 AMREX_ASSERT(m_zterrain.size() ==
static_cast<long int>(nx*ny));
418 nx = erf_get_single_value<int>(file,cnt); cnt++;
419 ny = erf_get_single_value<int>(file,cnt); cnt++;
420 amrex::Print()<<
"Expecting " << nx <<
" values of x, " <<
421 ny <<
" values of y, and " <<
422 nx*ny <<
" values of z" << std::endl;
425 m_xterrain.resize(nx);
426 m_yterrain.resize(ny);
427 m_zterrain.resize(nx * ny);
428 for (
int n = 0; n < nx; n++) {
429 m_xterrain[n] = erf_get_single_value<amrex::Real>(file,cnt);
432 for (
int n = 0; n < ny; n++) {
433 m_yterrain[n] = erf_get_single_value<amrex::Real>(file,cnt);
436 for (
int n = 0; n < nx * ny; n++) {
437 m_zterrain[n] = erf_get_single_value<amrex::Real>(file,cnt);
446 amrex::ParallelDescriptor::Bcast(&nx,1,amrex::ParallelDescriptor::IOProcessorNumber());
447 amrex::ParallelDescriptor::Bcast(&ny,1,amrex::ParallelDescriptor::IOProcessorNumber());
451 m_xterrain.resize(nx);
452 m_yterrain.resize(ny);
453 m_zterrain.resize(nz);
455 amrex::ParallelDescriptor::Bcast(m_xterrain.data(),nx,amrex::ParallelDescriptor::IOProcessorNumber());
456 amrex::ParallelDescriptor::Bcast(m_yterrain.data(),ny,amrex::ParallelDescriptor::IOProcessorNumber());
457 amrex::ParallelDescriptor::Bcast(m_zterrain.data(),nz,amrex::ParallelDescriptor::IOProcessorNumber());
460 amrex::Gpu::DeviceVector<amrex::Real> d_xterrain(nx),d_yterrain(ny),d_zterrain(nz);
461 amrex::Gpu::copy(amrex::Gpu::hostToDevice, m_xterrain.begin(), m_xterrain.end(), d_xterrain.begin());
462 amrex::Gpu::copy(amrex::Gpu::hostToDevice, m_yterrain.begin(), m_yterrain.end(), d_yterrain.begin());
463 amrex::Gpu::copy(amrex::Gpu::hostToDevice, m_zterrain.begin(), m_zterrain.end(), d_zterrain.begin());
469 auto dx = geom.CellSizeArray();
470 auto ProbLoArr = geom.ProbLoArray();
472 int ilo = geom.Domain().smallEnd(0);
473 int jlo = geom.Domain().smallEnd(1);
474 int klo = geom.Domain().smallEnd(2);
475 int ihi = geom.Domain().bigEnd(0) + 1;
476 int jhi = geom.Domain().bigEnd(1) + 1;
478 amrex::Box
zbx = terrain_fab.box();
479 amrex::Array4<amrex::Real>
const& z_arr = terrain_fab.array();
484 int ii = amrex::min(amrex::max(i,ilo),ihi);
485 int jj = amrex::min(amrex::max(j,jlo),jhi);
491 int ind11, ind12, ind21, ind22;
494 int iindex_terrain=-1;
495 int jindex_terrain=-1;
501 for (
int it = 0; it < ny && jindex_terrain == -1; it++) {
503 jindex_terrain = it-1;
506 if (jindex_terrain == -1) {
507 jindex_terrain = ny-1;
510 int gstart = (jindex_terrain )*nx;
511 int gend = (jindex_terrain+1)*nx-1;
512 for (
int it = gstart; it <= gend && iindex_terrain == -1; it++) {
514 iindex_terrain = it-gstart;
519 ind11 = jindex_terrain*nx + iindex_terrain;
526 y1 = d_yt[jindex_terrain];
527 y2 = d_yt[jindex_terrain+1];
536 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];
540 for (
int it = 0; it < ny && jindex_terrain == -1; it++) {
542 jindex_terrain = it-1;
545 if (jindex_terrain == -1) {
546 jindex_terrain = ny-1;
549 for (
int it = 0; it < nx && iindex_terrain == -1; it++) {
551 iindex_terrain = it-1;
554 if (iindex_terrain == -1) {
555 iindex_terrain = nx-1;
559 x1 = d_xt[iindex_terrain];
560 x2 = d_xt[iindex_terrain+1];
561 y1 = d_yt[jindex_terrain];
562 y2 = d_yt[jindex_terrain+1];
568 ind11 = iindex_terrain * ny + jindex_terrain;
569 ind21 = std::min(iindex_terrain+1,nx-1) * ny + jindex_terrain;
571 ind12 = iindex_terrain * ny + std::min(jindex_terrain+1,ny-1);
572 ind22 = std::min(iindex_terrain+1,nx-1) * ny + std::min(jindex_terrain+1,ny-1);
578 ind11 = jindex_terrain * nx + iindex_terrain;
579 ind12 = std::min(jindex_terrain+1,ny-1) * nx + iindex_terrain;
581 ind21 = jindex_terrain * nx + std::min(iindex_terrain+1,nx-1);
582 ind22 = std::min(jindex_terrain+1,ny-1) * nx + std::min(iindex_terrain+1,nx-1);
585 if (iindex_terrain == nx-1 && jindex_terrain == ny-1)
587 z_arr(i,j,
klo) = d_zt[ind11];
589 else if (iindex_terrain != nx-1 && jindex_terrain == ny-1)
594 z_arr(i,j,
klo) = (w_11*d_zt[ind11] + w_21*d_zt[ind21])/denom;
596 else if (iindex_terrain == nx-1 && jindex_terrain != ny-1)
601 z_arr(i,j,
klo) = (w_11*d_zt[ind11] + w_12*d_zt[ind12])/denom;
610 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 auto dx
Definition: ERF_InitCustomPertVels_ParticleTests.H:15
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);} } })
AMREX_ALWAYS_ASSERT(bx.length()[2]==khi+1)
const Box zbx
Definition: ERF_SetupDiff.H:9
amrex::Real Real
Definition: ERF_ShocInterface.H:19