|
| 1 | +#include "ERF_Prob.H" |
| 2 | +#include "ERF_EOS.H" |
| 3 | +#include "ERF_TerrainMetrics.H" |
| 4 | + |
| 5 | +using namespace amrex; |
| 6 | + |
| 7 | +std::unique_ptr<ProblemBase> |
| 8 | +amrex_probinit( |
| 9 | + const amrex_real* /*problo*/, |
| 10 | + const amrex_real* /*probhi*/) |
| 11 | +{ |
| 12 | + return std::make_unique<Problem>(); |
| 13 | +} |
| 14 | + |
| 15 | +Problem::Problem() |
| 16 | +{ |
| 17 | + // Parse params |
| 18 | + ParmParse pp("prob"); |
| 19 | + pp.query("T_0", parms.T_0); |
| 20 | + pp.query("U_0", parms.U_0); |
| 21 | + pp.query("x_c", parms.x_c); |
| 22 | + pp.query("z_c", parms.z_c); |
| 23 | + pp.query("x_r", parms.x_r); |
| 24 | + pp.query("z_r", parms.z_r); |
| 25 | + pp.query("T_pert", parms.T_pert); |
| 26 | + |
| 27 | + init_base_parms(parms.rho_0, parms.T_0); |
| 28 | +} |
| 29 | + |
| 30 | +void |
| 31 | +Problem::init_custom_pert( |
| 32 | + const Box& bx, |
| 33 | + const Box& xbx, |
| 34 | + const Box& ybx, |
| 35 | + const Box& zbx, |
| 36 | + Array4<Real const> const& /*state*/, |
| 37 | + Array4<Real > const& state_pert, |
| 38 | + Array4<Real > const& x_vel_pert, |
| 39 | + Array4<Real > const& y_vel_pert, |
| 40 | + Array4<Real > const& z_vel_pert, |
| 41 | + Array4<Real > const& r_hse, |
| 42 | + Array4<Real > const& p_hse, |
| 43 | + Array4<Real const> const& z_nd, |
| 44 | + Array4<Real const> const& z_cc, |
| 45 | + GeometryData const& geomdata, |
| 46 | + Array4<Real const> const& /*mf_m*/, |
| 47 | + Array4<Real const> const& /*mf_u*/, |
| 48 | + Array4<Real const> const& /*mf_v*/, |
| 49 | + const SolverChoice& sc) |
| 50 | +{ |
| 51 | + const int klo = geomdata.Domain().smallEnd()[2]; |
| 52 | + const int khi = geomdata.Domain().bigEnd()[2]; |
| 53 | + |
| 54 | + const bool use_moisture = (sc.moisture_type != MoistureType::None); |
| 55 | + |
| 56 | + AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1); |
| 57 | + |
| 58 | + ParallelFor(bx, [=, parms_d=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept |
| 59 | + { |
| 60 | + // Geometry (note we must include these here to get the data on device) |
| 61 | + const auto prob_lo = geomdata.ProbLo(); |
| 62 | + const auto dx = geomdata.CellSize(); |
| 63 | + const Real x = prob_lo[0] + (i + 0.5) * dx[0]; |
| 64 | + const Real z = z_cc(i,j,k); |
| 65 | + |
| 66 | + // Temperature that satisfies the EOS given the hydrostatically balanced (r,p) |
| 67 | + const Real Tbar_hse = p_hse(i,j,k) / (R_d * r_hse(i,j,k)); |
| 68 | + |
| 69 | + Real L = std::sqrt( |
| 70 | + std::pow((x - parms_d.x_c)/parms_d.x_r, 2) + |
| 71 | + std::pow((z - parms_d.z_c)/parms_d.z_r, 2) |
| 72 | + ); |
| 73 | + if (L <= 1.0) { |
| 74 | + Real dT = parms_d.T_pert * (std::cos(PI*L) + 1.0)/2.0; |
| 75 | + |
| 76 | + // Note: dT is a temperature perturbation, theta_perturbed is base state + perturbation in theta |
| 77 | + Real theta_perturbed = (Tbar_hse+dT)*std::pow(p_0/p_hse(i,j,k), R_d/parms_d.C_p); |
| 78 | + |
| 79 | + // This version perturbs rho but not p |
| 80 | + state_pert(i, j, k, Rho_comp) = getRhoThetagivenP(p_hse(i,j,k)) / theta_perturbed - r_hse(i,j,k); |
| 81 | + } |
| 82 | + |
| 83 | + // Set scalar = 0 everywhere |
| 84 | + state_pert(i, j, k, RhoScalar_comp) = 0.0; |
| 85 | + |
| 86 | + if (use_moisture) { |
| 87 | + state_pert(i, j, k, RhoQ1_comp) = 0.0; |
| 88 | + state_pert(i, j, k, RhoQ2_comp) = 0.0; |
| 89 | + } |
| 90 | + }); |
| 91 | + |
| 92 | + // Set the x-velocity |
| 93 | + ParallelFor(xbx, [=, parms_d=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept |
| 94 | + { |
| 95 | + Real ztop = z_nd(i,j,khi+1); |
| 96 | + Real zht = z_nd(i,j,klo); |
| 97 | + x_vel_pert(i, j, k) = parms_d.U_0 * ztop / (ztop - zht); |
| 98 | + }); |
| 99 | + |
| 100 | + // Set the y-velocity |
| 101 | + ParallelFor(ybx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept |
| 102 | + { |
| 103 | + y_vel_pert(i, j, k) = 0.0; |
| 104 | + }); |
| 105 | + |
| 106 | + const auto dx = geomdata.CellSize(); |
| 107 | + amrex::GpuArray<Real, AMREX_SPACEDIM> dxInv; |
| 108 | + dxInv[0] = 1. / dx[0]; |
| 109 | + dxInv[1] = 1. / dx[1]; |
| 110 | + dxInv[2] = 1. / dx[2]; |
| 111 | + |
| 112 | + // Set the z-velocity from impenetrable condition |
| 113 | + ParallelFor(zbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept |
| 114 | + { |
| 115 | + z_vel_pert(i, j, k) = WFromOmega(i, j, k, 0.0, x_vel_pert, y_vel_pert, z_nd, dxInv); |
| 116 | + }); |
| 117 | + |
| 118 | + amrex::Gpu::streamSynchronize(); |
| 119 | + |
| 120 | +} |
0 commit comments