From 31ed2ef39782943eb4a4c325a0e43ebe043e3072 Mon Sep 17 00:00:00 2001 From: "Aaron M. Lattanzi" <103702284+AMLattanzi@users.noreply.github.com> Date: Thu, 4 Jan 2024 14:43:33 -0800 Subject: [PATCH] BF02 Bubble (#1364) * Add BF02 verif studies to bubble dir. * Fix device functions in bubble. * Results generated. --------- Co-authored-by: Aaron Lattanzi --- Exec/RegTests/Bubble/inputs_BF02_dry_bubble | 73 +++ Exec/RegTests/Bubble/inputs_BF02_moist_bubble | 82 ++++ Exec/RegTests/Bubble/prob.H | 113 +++-- Exec/RegTests/Bubble/prob.cpp | 426 ++++++++++++++---- Exec/SquallLine_2D/prob.cpp | 5 +- Source/ERF.H | 2 +- Source/ERF_Constants.H | 1 + Source/IO/Plotfile.cpp | 24 + Source/TimeIntegration/TI_slow_rhs_fun.H | 2 +- Source/Utils/Microphysics_Utils.H | 9 +- 10 files changed, 614 insertions(+), 123 deletions(-) create mode 100644 Exec/RegTests/Bubble/inputs_BF02_dry_bubble create mode 100644 Exec/RegTests/Bubble/inputs_BF02_moist_bubble diff --git a/Exec/RegTests/Bubble/inputs_BF02_dry_bubble b/Exec/RegTests/Bubble/inputs_BF02_dry_bubble new file mode 100644 index 000000000..82a8d5ac1 --- /dev/null +++ b/Exec/RegTests/Bubble/inputs_BF02_dry_bubble @@ -0,0 +1,73 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +max_step = 2000 +stop_time = 3600.0 + +amrex.fpe_trap_invalid = 1 + +fabarray.mfiter_tile_size = 1024 1024 1024 + +# PROBLEM SIZE & GEOMETRY +geometry.prob_extent = 20000.0 312.5 10000.0 +amr.n_cell = 256 4 128 +geometry.is_periodic = 0 1 0 +xlo.type = "SlipWall" +xhi.type = "SlipWall" +zlo.type = "SlipWall" +zhi.type = "SlipWall" + +# TIME STEP CONTROL +erf.fixed_dt = 0.5 +erf.fixed_mri_dt_ratio = 4 +#erf.no_substepping = 1 +#erf.fixed_dt = 0.1 + +# DIAGNOSTICS & VERBOSITY +erf.sum_interval = 1 # timesteps between computing mass +erf.v = 1 # verbosity in ERF.cpp +amr.v = 1 # verbosity in Amr.cpp + +# REFINEMENT / REGRIDDING +amr.max_level = 0 # maximum level number allowed + +# CHECKPOINT FILES +erf.check_file = chk # root name of checkpoint file +erf.check_int = 100 # number of timesteps between checkpoints + +# PLOTFILES +erf.plot_file_1 = plt # prefix of plotfile name +erf.plot_int_1 = 100 # number of timesteps between plotfiles +erf.plot_vars_1 = density rhotheta rhoadv_0 x_velocity y_velocity z_velocity pressure theta scalar temp pres_hse dens_hse pert_pres pert_dens + +# SOLVER CHOICES +erf.use_gravity = true +erf.use_coriolis = false +erf.use_rayleigh_damping = false + +erf.dycore_horiz_adv_type = "Upwind_5th" +erf.dycore_vert_adv_type = "Upwind_5th" +erf.dryscal_horiz_adv_type = "Upwind_5th" +erf.dryscal_vert_adv_type = "Upwind_5th" + +# PHYSICS OPTIONS +erf.les_type = "None" +erf.pbl_type = "None" +erf.moisture_model = "NullMoist" +erf.buoyancy_type = 1 +erf.use_moist_background = false + +erf.molec_diff_type = "ConstantAlpha" +erf.rho0_trans = 1.0 # [kg/m^3], used to convert input diffusivities +erf.dynamicViscosity = 1.0 # [kg/(m-s)] ==> nu = 75.0 m^2/s +erf.alpha_T = 1.0 # [m^2/s] +erf.alpha_C = 1.0 + +# PROBLEM PARAMETERS (optional) +# warm bubble input +prob.T_pert = 2.0 # theta pert magnitude +prob.x_c = 10000.0 +prob.z_c = 2000.0 +prob.x_r = 2000.0 +prob.z_r = 2000.0 +prob.T_0 = 300.0 +prob.do_moist_bubble = false +prob.T_pert_is_airtemp = false # Perturb theta \ No newline at end of file diff --git a/Exec/RegTests/Bubble/inputs_BF02_moist_bubble b/Exec/RegTests/Bubble/inputs_BF02_moist_bubble new file mode 100644 index 000000000..3f923dc3f --- /dev/null +++ b/Exec/RegTests/Bubble/inputs_BF02_moist_bubble @@ -0,0 +1,82 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +max_step = 2000 +stop_time = 3600.0 + +amrex.fpe_trap_invalid = 1 + +fabarray.mfiter_tile_size = 1024 1024 1024 + +# PROBLEM SIZE & GEOMETRY +geometry.prob_extent = 20000.0 312.5 10000.0 +amr.n_cell = 256 4 128 +geometry.is_periodic = 0 1 0 +xlo.type = "SlipWall" +xhi.type = "SlipWall" +zlo.type = "SlipWall" +zhi.type = "SlipWall" + +# TIME STEP CONTROL +erf.fixed_dt = 0.5 +erf.fixed_mri_dt_ratio = 4 +#erf.no_substepping = 1 +#erf.fixed_dt = 0.1 + +# DIAGNOSTICS & VERBOSITY +erf.sum_interval = 1 # timesteps between computing mass +erf.v = 1 # verbosity in ERF.cpp +amr.v = 1 # verbosity in Amr.cpp + +# REFINEMENT / REGRIDDING +amr.max_level = 0 # maximum level number allowed + +# CHECKPOINT FILES +erf.check_file = chk # root name of checkpoint file +erf.check_int = 100 # number of timesteps between checkpoints + +# PLOTFILES +erf.plot_file_1 = plt # prefix of plotfile name +erf.plot_int_1 = 100 # number of timesteps between plotfiles +erf.plot_vars_1 = density rhotheta rhoQ1 rhoQ2 rhoadv_0 x_velocity y_velocity z_velocity pressure theta scalar temp pres_hse dens_hse pert_pres pert_dens eq_pot_temp qt qv qc + +# SOLVER CHOICES +erf.use_gravity = true +erf.use_coriolis = false +erf.use_rayleigh_damping = false + +erf.dycore_horiz_adv_type = "Upwind_3rd" +erf.dycore_vert_adv_type = "Upwind_3rd" +erf.dryscal_horiz_adv_type = "Upwind_3rd" +erf.dryscal_vert_adv_type = "Upwind_3rd" +erf.moistscal_horiz_adv_type = "Upwind_3rd" +erf.moistscal_vert_adv_type = "Upwind_3rd" + +# PHYSICS OPTIONS +erf.les_type = "None" +erf.pbl_type = "None" +erf.moisture_model = "FastEddy" +erf.buoyancy_type = 1 +erf.use_moist_background = true + +erf.molec_diff_type = "ConstantAlpha" +erf.rho0_trans = 1.0 # [kg/m^3], used to convert input diffusivities +erf.dynamicViscosity = 1.0 # [kg/(m-s)] ==> nu = 75.0 m^2/s +erf.alpha_T = 1.0 # [m^2/s] +erf.alpha_C = 1.0 + +# INITIAL CONDITIONS +#erf.init_type = "input_sounding" +#erf.input_sounding_file = "BF02_moist_sounding" +#erf.init_sounding_ideal = true + +# PROBLEM PARAMETERS (optional) +# warm bubble input +prob.x_c = 10000.0 +prob.z_c = 2000.0 +prob.x_r = 2000.0 +prob.z_r = 2000.0 +prob.T_0 = 300.0 + +prob.do_moist_bubble = true +prob.theta_pert = 2.0 +prob.qt_init = 0.02 +prob.eq_pot_temp = 320.0 \ No newline at end of file diff --git a/Exec/RegTests/Bubble/prob.H b/Exec/RegTests/Bubble/prob.H index e8fe13dd7..9df0b7400 100644 --- a/Exec/RegTests/Bubble/prob.H +++ b/Exec/RegTests/Bubble/prob.H @@ -9,31 +9,37 @@ #include "EOS.H" struct ProbParm : ProbParmDefaults { - // background conditions - // if init_type != "" then these are perturbations and should be 0 - amrex::Real T_0 = 300.0; - amrex::Real U_0 = 0.0; - amrex::Real V_0 = 0.0; - amrex::Real W_0 = 0.0; - - // center of thermal perturbation - amrex::Real x_c = 0.0; - amrex::Real y_c = 0.0; - amrex::Real z_c = 0.0; - - // radial extent of thermal perturbation - amrex::Real x_r = 0.0; - amrex::Real y_r = 0.0; - amrex::Real z_r = 0.0; - - // perturbation temperature - amrex::Real T_pert = -15.0; - bool T_pert_is_airtemp = true; // T_pert input is air temperature - bool perturb_rho = true; // not rho*theta (i.e., p is constant); otherwise perturb rho*theta - - // rayleigh damping - amrex::Real dampcoef = 0.0; // inverse time scale [1/s] - amrex::Real zdamp = 5000.0; // damping depth [m] from model top + // background conditions + // if init_type != "" then these are perturbations and should be 0 + amrex::Real T_0 = 300.0; + amrex::Real U_0 = 0.0; + amrex::Real V_0 = 0.0; + amrex::Real W_0 = 0.0; + + // center of thermal perturbation + amrex::Real x_c = 0.0; + amrex::Real y_c = 0.0; + amrex::Real z_c = 0.0; + + // radial extent of thermal perturbation + amrex::Real x_r = 0.0; + amrex::Real y_r = 0.0; + amrex::Real z_r = 0.0; + + // perturbation temperature + amrex::Real T_pert = -15.0; + bool T_pert_is_airtemp = true; // T_pert input is air temperature + bool perturb_rho = true; // not rho*theta (i.e., p is constant); otherwise perturb rho*theta + + // rayleigh damping + amrex::Real dampcoef = 0.0; // inverse time scale [1/s] + amrex::Real zdamp = 5000.0; // damping depth [m] from model top + + // Moist bubble params + bool do_moist_bubble = false; + amrex::Real theta_pert = 2.0; + amrex::Real qt_init = 0.02; + amrex::Real eq_pot_temp = 320.0; }; // namespace ProbParm class Problem : public ProblemBase @@ -62,11 +68,66 @@ public: amrex::Array4 const& mf_v, const SolverChoice& sc) override; + void erf_init_dens_hse_moist (amrex::MultiFab& rho_hse, + std::unique_ptr& z_phys_nd, + amrex::Geometry const& geom) override; + void init_custom_terrain ( const amrex::Geometry& geom, amrex::MultiFab& z_phys_nd, const amrex::Real& time) override; + AMREX_FORCE_INLINE + AMREX_GPU_HOST_DEVICE + amrex::Real compute_theta (const amrex::Real T_b, const amrex::Real p_b); + + AMREX_FORCE_INLINE + AMREX_GPU_HOST_DEVICE + amrex::Real compute_temperature (const amrex::Real p_b); + + AMREX_FORCE_INLINE + AMREX_GPU_HOST_DEVICE + amrex::Real compute_F_for_temp(const amrex::Real T_b, const amrex::Real p_b); + + AMREX_FORCE_INLINE + AMREX_GPU_HOST_DEVICE + amrex::Real compute_p_k (amrex::Real& p_k, + const amrex::Real p_k_minus_1, + amrex::Real& theta_k, + amrex::Real& rho_k, + amrex::Real& q_v_k, + amrex::Real& T_dp, + amrex::Real& T_b, + const amrex::Real dz, + const amrex::Real rho_k_minus_1); + AMREX_FORCE_INLINE + AMREX_GPU_HOST_DEVICE + amrex::Real compute_F (const amrex::Real& p_k, + const amrex::Real& p_k_minus_1, + amrex::Real &theta_k, + amrex::Real& rho_k, + amrex::Real& q_v_k, + amrex::Real& T_dp, + amrex::Real& T_b, + const amrex::Real& dz, + const amrex::Real& rho_k_minus_1); + + AMREX_FORCE_INLINE + AMREX_GPU_HOST_DEVICE + void compute_rho (const amrex::Real& pressure, + amrex::Real &theta, + amrex::Real& rho, + amrex::Real& q_v, + amrex::Real& T_dp, + amrex::Real& T_b); + + void init_isentropic_hse_no_terrain(amrex::Real *theta, + amrex::Real* r, + amrex::Real* p, + amrex::Real *q_v, + const amrex::Real& dz, + const int& khi); + #include "Prob/init_rayleigh_damping.H" protected: @@ -111,7 +172,7 @@ perturb_rho_theta (const amrex::Real x, dT = 0.0; } else { - dT = pp.T_pert * (std::cos(PI*L) + 1.0)/2.0; + dT = pp.T_pert * std::pow(cos(PI*L/2.0),2); } // Temperature that satisfies the EOS given the hydrostatically balanced (r,p) diff --git a/Exec/RegTests/Bubble/prob.cpp b/Exec/RegTests/Bubble/prob.cpp index a70d87a1f..de5ff15d3 100644 --- a/Exec/RegTests/Bubble/prob.cpp +++ b/Exec/RegTests/Bubble/prob.cpp @@ -1,4 +1,6 @@ #include "prob.H" +#include "Microphysics_Utils.H" +#include "ERF_Constants.H" using namespace amrex; @@ -30,9 +32,234 @@ Problem::Problem() pp.query("dampcoef", parms.dampcoef); pp.query("zdamp", parms.zdamp); + pp.query("do_moist_bubble", parms.do_moist_bubble); + pp.query("theta_pert", parms.theta_pert); + pp.query("qt_init", parms.qt_init); + pp.query("eq_pot_temp", parms.eq_pot_temp); + init_base_parms(parms.rho_0, parms.T_0); } +AMREX_FORCE_INLINE +AMREX_GPU_HOST_DEVICE +Real compute_saturation_pressure (const Real T_b) +{ + return erf_esatw(T_b)*100.0; +} + +AMREX_FORCE_INLINE +AMREX_GPU_HOST_DEVICE +Real compute_relative_humidity () +{ + return 1.0; +} + +AMREX_FORCE_INLINE +AMREX_GPU_HOST_DEVICE +Real compute_vapor_pressure (const Real p_s, const Real RH) +{ + return p_s*RH; +} + +AMREX_FORCE_INLINE +AMREX_GPU_HOST_DEVICE +Real vapor_mixing_ratio (const Real p_b, const Real T_b, const Real RH) +{ + Real p_s = compute_saturation_pressure(T_b); + Real p_v = compute_vapor_pressure(p_s, RH); + Real q_v = Rd_on_Rv*p_v/(p_b - p_v); + return q_v; +} + +AMREX_FORCE_INLINE +AMREX_GPU_HOST_DEVICE +Real Problem::compute_F_for_temp(const Real T_b, const Real p_b) +{ + Real q_t = parms.qt_init; + Real fac = Cp_d + Cp_l*q_t; + Real RH = compute_relative_humidity(); + Real q_v = vapor_mixing_ratio(p_b, T_b, RH); + Real p_s = compute_saturation_pressure(T_b); + Real p_v = compute_vapor_pressure(p_s, RH); + return parms.eq_pot_temp - T_b*std::pow((p_b - p_v)/p_0, -R_d/fac)*std::exp(L_v*q_v/(fac*T_b)); +} + +AMREX_FORCE_INLINE +AMREX_GPU_HOST_DEVICE +Real Problem::compute_temperature (const Real p_b) +{ + Real T_b = 200.0, delta_T; // Initial guess + for(int iter=0; iter<20; iter++) + { + Real F = compute_F_for_temp(T_b, p_b); + Real F_plus_dF = compute_F_for_temp(T_b+1e-10, p_b); + Real F_prime = (F_plus_dF - F)/1e-10; + delta_T = -F/F_prime; + T_b = T_b + delta_T; + } + + if(std::fabs(delta_T) > 1e-8){ + amrex::Abort("Newton Raphson for temperature could not converge"); + } + + return T_b; +} + +AMREX_FORCE_INLINE +AMREX_GPU_HOST_DEVICE +Real compute_dewpoint_temperature (const Real T_b, const Real RH) +{ + Real T_dp, gamma, T; + T = T_b - 273.15; + + Real b = 18.678, c = 257.14, d = 234.5; + gamma = log(RH*exp((b - T/d)*T/(c + T))); + + T_dp = c*gamma/(b - gamma); + + return T_dp; +} + +AMREX_FORCE_INLINE +AMREX_GPU_HOST_DEVICE +Real Problem::compute_theta(const Real T_b, const Real p_b) +{ + return T_b*std::pow(p_0/p_b, R_d/Cp_d); +} + +AMREX_FORCE_INLINE +AMREX_GPU_HOST_DEVICE +Real compute_temperature_from_theta(const Real theta, const Real p) +{ + return theta*std::pow(p/p_0, R_d/Cp_d); +} + +AMREX_FORCE_INLINE +AMREX_GPU_HOST_DEVICE +void Problem::compute_rho (const Real& pressure, Real& theta, Real& rho, Real& q_v, Real& T_dp, Real& T_b) +{ + T_b = compute_temperature(pressure); + theta = compute_theta(T_b, pressure); + Real RH = compute_relative_humidity(); + q_v = vapor_mixing_ratio(pressure, T_b, RH); + rho = pressure/(R_d*T_b*(1.0 + (R_v/R_d)*q_v)); + rho = rho*(1.0 + parms.qt_init); // q_t = 0.02 a given constant for this test case + T_dp = compute_dewpoint_temperature(T_b, RH); +} + +AMREX_FORCE_INLINE +AMREX_GPU_HOST_DEVICE +Real Problem::compute_F (const Real& p_k, const Real& p_k_minus_1, Real &theta_k, Real& rho_k, Real& q_v_k, + Real& T_dp, Real& T_b, const Real& dz, const Real& rho_k_minus_1) +{ + Real F; + + if(rho_k_minus_1 == 0.0) // This loop is for the first point above the ground + { + compute_rho(p_k, theta_k, rho_k, q_v_k, T_dp, T_b); + F = p_k - p_k_minus_1 + rho_k*CONST_GRAV*dz/2.0; + } + else + { + compute_rho(p_k, theta_k, rho_k, q_v_k, T_dp, T_b); + F = p_k - p_k_minus_1 + rho_k_minus_1*CONST_GRAV*dz/2.0 + rho_k*CONST_GRAV*dz/2.0; + } + + return F; +} + +AMREX_FORCE_INLINE +AMREX_GPU_HOST_DEVICE +Real Problem::compute_p_k (Real &p_k, const Real p_k_minus_1, Real &theta_k, Real& rho_k, Real& q_v_k, + Real& T_dp, Real& T_b, const Real dz, const Real rho_k_minus_1) +{ + Real delta_p_k; + + for(int iter=0; iter<20; iter++) + { + Real F = compute_F(p_k, p_k_minus_1, theta_k, rho_k, q_v_k, T_dp, T_b, dz, rho_k_minus_1); + Real F_plus_dF = compute_F(p_k+1e-10, p_k_minus_1, theta_k, rho_k, q_v_k, T_dp, T_b, dz, rho_k_minus_1); + Real F_prime = (F_plus_dF - F)/1e-10; + delta_p_k = -F/F_prime; + p_k = p_k + delta_p_k; + } + + if(std::fabs(delta_p_k) > 1e-8){ + amrex::Abort("Newton Raphson for pressure could not converge"); + } + + return p_k; +} + +void +Problem::init_isentropic_hse_no_terrain(Real *theta, Real* r, Real* p, Real *q_v, + const Real& dz, + const int& khi) +{ + Real T_b, T_dp; + + // Compute the quantities at z = 0.5*dz (first cell center) + p[0] = p_0; + compute_p_k(p[0], p_0, theta[0], r[0], q_v[0], T_dp, T_b, dz, 0.0); + + for (int k=1;k<=khi;k++){ + p[k] = p[k-1]; + compute_p_k(p[k], p[k-1], theta[k], r[k], q_v[k], T_dp, T_b, dz, r[k-1]); + } + + r[khi+1] = r[khi]; +} + +void +Problem::erf_init_dens_hse_moist (MultiFab& rho_hse, + std::unique_ptr& z_phys_nd, + Geometry const& geom) +{ + const Real dz = geom.CellSize()[2]; + const int khi = geom.Domain().bigEnd()[2]; + + // use_terrain = 1 + if (z_phys_nd) { + amrex::Abort("Terrain not implemented for moist Bubble case!"); + } else { // use_terrain = 0 + + // These are at cell centers (unstaggered) + Vector h_r(khi+2); + Vector h_p(khi+2); + Vector h_t(khi+2); + Vector h_q_v(khi+2); + + amrex::Gpu::DeviceVector d_r(khi+2); + amrex::Gpu::DeviceVector d_p(khi+2); + amrex::Gpu::DeviceVector d_t(khi+2); + amrex::Gpu::DeviceVector d_q_v(khi+2); + + init_isentropic_hse_no_terrain(h_t.data(), h_r.data(),h_p.data(), h_q_v.data(), dz, khi); + + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, h_r.begin(), h_r.end(), d_r.begin()); + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, h_p.begin(), h_p.end(), d_p.begin()); + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, h_t.begin(), h_t.end(), d_t.begin()); + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, h_q_v.begin(), h_q_v.end(), d_q_v.begin()); + + Real* r = d_r.data(); + +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(rho_hse,TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.growntilebox(1); + const Array4 rho_hse_arr = rho_hse[mfi].array(); + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + int kk = std::max(k,0); + rho_hse_arr(i,j,k) = r[kk]; + }); + } // mfi + } // no terrain +} + + void Problem::init_custom_pert( const Box& bx, @@ -60,130 +287,151 @@ Problem::init_custom_pert( AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1); const Real rho_sfc = p_0 / (R_d*parms.T_0); - //const Real thetabar = parms.T_0; const Real dz = geomdata.CellSize()[2]; const Real prob_lo_z = geomdata.ProbLo()[2]; const Real rdOcp = sc.rdOcp; -#if 0 - // These are at cell centers (unstaggered) - Vector h_r(khi+1); - Vector h_p(khi+1); - - amrex::Gpu::DeviceVector d_r(khi+1); - amrex::Gpu::DeviceVector d_p(khi+1); -#endif - amrex::Print() << "Bubble delta T = " << parms.T_pert << " K" << std::endl; amrex::Print() << " centered at (" - << parms.x_c << " " << parms.y_c << " " << parms.z_c << ")" << std::endl; + << parms.x_c << " " << parms.y_c << " " << parms.z_c << ")" << std::endl; amrex::Print() << " with extent (" - << parms.x_r << " " << parms.y_r << " " << parms.z_r << ")" << std::endl; + << parms.x_r << " " << parms.y_r << " " << parms.z_r << ")" << std::endl; if (parms.T_0 <= 0) { amrex::Print() << "Ignoring parms.T_0 = " << parms.T_0 - << ", background fields should have been initialized with erf.init_type" - << std::endl; + << ", background fields should have been initialized with erf.init_type" + << std::endl; } if (z_cc) { // nonflat terrain + if (parms.do_moist_bubble) { + amrex::Abort("Terrain not implemented for moist Bubble case!"); + } else { + amrex::ParallelFor(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + // Geometry (note we must include these here to get the data on device) + const auto prob_lo = geomdata.ProbLo(); + const auto dx = geomdata.CellSize(); -#if 0 - if (parms.T_0 > 0) - { - // Create a flat box with same horizontal extent but only one cell in vertical - Box b2d = surroundingNodes(bx); // Copy constructor - b2d.setRange(2,0); + const Real x = prob_lo[0] + (i + 0.5) * dx[0]; + const Real y = prob_lo[1] + (j + 0.5) * dx[1]; + const Real z = z_cc(i,j,k); - ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int) - { - Array1D r;; - Array1D p;; + perturb_rho_theta(x, y, z, p_hse(i,j,k), r_hse(i,j,k), + parms, rdOcp, + state(i, j, k, Rho_comp), + state(i, j, k, RhoTheta_comp)); - init_isentropic_hse_terrain(i,j,rho_sfc,thetabar,&(r(0)),&(p(0)),z_cc,khi); + state(i, j, k, RhoScalar_comp) = 0.0; - for (int k = 0; k <= khi; k++) { - r_hse(i,j,k) = r(k); - p_hse(i,j,k) = p(k); + if (use_moisture) { + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; } - r_hse(i,j, -1) = r_hse(i,j,0); - r_hse(i,j,khi+1) = r_hse(i,j,khi); + }); - } -#endif + } // do_moist_bubble + } else { + if (parms.do_moist_bubble) { + Vector h_r(khi+2); + Vector h_p(khi+2); + Vector h_t(khi+2); + Vector h_q_v(khi+2); - amrex::ParallelFor(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept - { - // Geometry (note we must include these here to get the data on device) - const auto prob_lo = geomdata.ProbLo(); - const auto dx = geomdata.CellSize(); + amrex::Gpu::DeviceVector d_r(khi+2); + amrex::Gpu::DeviceVector d_p(khi+2); + amrex::Gpu::DeviceVector d_t(khi+2); + amrex::Gpu::DeviceVector d_q_v(khi+2); - const Real x = prob_lo[0] + (i + 0.5) * dx[0]; - const Real y = prob_lo[1] + (j + 0.5) * dx[1]; - const Real z = z_cc(i,j,k); + init_isentropic_hse_no_terrain(h_t.data(), h_r.data(),h_p.data(),h_q_v.data(),dz, khi); - perturb_rho_theta(x, y, z, p_hse(i,j,k), r_hse(i,j,k), - parms, rdOcp, - state(i, j, k, Rho_comp), - state(i, j, k, RhoTheta_comp)); + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, h_r.begin(), h_r.end(), d_r.begin()); + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, h_p.begin(), h_p.end(), d_p.begin()); + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, h_t.begin(), h_t.end(), d_t.begin()); + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, h_q_v.begin(), h_q_v.end(), d_q_v.begin()); - state(i, j, k, RhoScalar_comp) = 0.0; + Real* theta_back = d_t.data(); + Real* p_back = d_p.data(); + Real* q_v_back = d_q_v.data(); - if (use_moisture) { - state(i, j, k, RhoQ1_comp) = 0.0; - state(i, j, k, RhoQ2_comp) = 0.0; - } - }); - } else { + const Real x_c = parms.x_c, y_c = parms.y_c, z_c = parms.z_c; + const Real x_r = parms.x_r, y_r = parms.y_r, z_r = parms.z_r; + const Real theta_pert = parms.theta_pert; -#if 0 - if (parms.T_0 > 0) - { - init_isentropic_hse(rho_sfc,thetabar,h_r.data(),h_p.data(),dz,prob_lo_z,khi); + amrex::ParallelFor(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) + { + // Geometry (note we must include these here to get the data on device) + const auto prob_lo = geomdata.ProbLo(); + const auto dx = geomdata.CellSize(); + const amrex::Real x = prob_lo[0] + (i + 0.5) * dx[0]; + const amrex::Real y = prob_lo[1] + (j + 0.5) * dx[1]; + const amrex::Real z = prob_lo[2] + (k + 0.5) * dx[2]; + + Real rad, delta_theta, theta_total, rho, RH; + + // Introduce the warm bubble. Assume that the bubble is pressure matched with the background + rad = 0.0; + if (x_r > 0) rad += std::pow((x - x_c)/x_r, 2); + if (y_r > 0) rad += std::pow((y - y_c)/y_r, 2); + if (z_r > 0) rad += std::pow((z - z_c)/z_r, 2); + rad = std::sqrt(rad); + + if(rad <= 1.0){ + delta_theta = theta_pert*std::pow(cos(PI*rad/2.0),2); + } + else{ + delta_theta = 0.0; + } - amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, h_r.begin(), h_r.end(), d_r.begin()); - amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, h_p.begin(), h_p.end(), d_p.begin()); + theta_total = theta_back[k]*(delta_theta/300.0 + 1); + Real T = compute_temperature_from_theta(theta_total, p_back[k]); + rho = p_back[k]/(R_d*T*(1.0 + (R_v/R_d)*q_v_back[k])); + RH = compute_relative_humidity(); + Real q_v_hot = vapor_mixing_ratio(p_back[k], T, RH); - Real* r = d_r.data(); - Real* p = d_p.data(); + // Compute background quantities + Real T_back = compute_temperature_from_theta(theta_back[k], p_back[k]); + Real rho_back = p_back[k]/(R_d*T_back*(1.0 + (R_v/R_d)*q_v_back[k])); - amrex::ParallelFor(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int) noexcept - { - for (int k = 0; k <= khi; k++) { - r_hse(i,j,k) = r[k]; - p_hse(i,j,k) = p[k]; - } - r_hse(i,j, -1) = r_hse(i,j,0); - r_hse(i,j,khi+1) = r_hse(i,j,khi); - }); - } -#endif + // This version perturbs rho but not p + state(i, j, k, RhoTheta_comp) = rho*theta_total - rho_back*theta_back[k]*(1.0 + (R_v/R_d)*q_v_back[k]); + state(i, j, k, Rho_comp) = rho - rho_back*(1.0 + parms.qt_init); - amrex::ParallelFor(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept - { - // Geometry (note we must include these here to get the data on device) - const auto prob_lo = geomdata.ProbLo(); - const auto dx = geomdata.CellSize(); + // Set scalar = 0 everywhere + state(i, j, k, RhoScalar_comp) = 0.0; - const Real x = prob_lo[0] + (i + 0.5) * dx[0]; - const Real y = prob_lo[1] + (j + 0.5) * dx[1]; - const Real z = prob_lo[2] + (k + 0.5) * dx[2]; + // mean states + state(i, j, k, RhoQ1_comp) = rho*q_v_hot; + state(i, j, k, RhoQ2_comp) = rho*(parms.qt_init - q_v_hot); - perturb_rho_theta(x, y, z, p_hse(i,j,k), r_hse(i,j,k), - parms, rdOcp, - state(i, j, k, Rho_comp), - state(i, j, k, RhoTheta_comp)); + }); + } else { + amrex::ParallelFor(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + // Geometry (note we must include these here to get the data on device) + const auto prob_lo = geomdata.ProbLo(); + const auto dx = geomdata.CellSize(); - state(i, j, k, RhoScalar_comp) = 0.0; + const Real x = prob_lo[0] + (i + 0.5) * dx[0]; + const Real y = prob_lo[1] + (j + 0.5) * dx[1]; + const Real z = prob_lo[2] + (k + 0.5) * dx[2]; - if (use_moisture) { - state(i, j, k, RhoQ1_comp) = 0.0; - state(i, j, k, RhoQ2_comp) = 0.0; - } - }); - } + perturb_rho_theta(x, y, z, p_hse(i,j,k), r_hse(i,j,k), + parms, rdOcp, + state(i, j, k, Rho_comp), + state(i, j, k, RhoTheta_comp)); + + state(i, j, k, RhoScalar_comp) = 0.0; + + if (use_moisture) { + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; + } + }); + } // do_moist_bubble + } // use_terrain const Real u0 = parms.U_0; const Real v0 = parms.V_0; diff --git a/Exec/SquallLine_2D/prob.cpp b/Exec/SquallLine_2D/prob.cpp index 8e5445f22..a8c07a79f 100644 --- a/Exec/SquallLine_2D/prob.cpp +++ b/Exec/SquallLine_2D/prob.cpp @@ -1,4 +1,5 @@ #include "prob.H" +#include using namespace amrex; @@ -57,7 +58,7 @@ Real compute_relative_humidity (const Real z, const Real height, const Real z_tr { Real p_s = compute_saturation_pressure(T_b); - Real q_s = 0.622*p_s/(p_b - p_s); + Real q_s = Rd_on_Rv*p_s/(p_b - p_s); if(z <= height){ return 0.014/q_s; @@ -82,7 +83,7 @@ Real vapor_mixing_ratio (const Real z, const Real height, const Real p_b, const Real p_s = compute_saturation_pressure(T_b); Real p_v = compute_vapor_pressure(p_s, RH); - Real q_v = 0.622*p_v/(p_b - p_v); + Real q_v = Rd_on_Rv*p_v/(p_b - p_v); if(z < height){ return 0.014; diff --git a/Source/ERF.H b/Source/ERF.H index 6b96dc255..4c67f85a4 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -650,7 +650,7 @@ private: // Note that the order of variable names here must match the order in Derive.cpp const amrex::Vector derived_names {"soundspeed", "temp", "theta", "KE", "QKE", "scalar", - "pres_hse", "dens_hse", "pressure", "pert_pres", "pert_dens", + "pres_hse", "dens_hse", "pressure", "pert_pres", "pert_dens", "eq_pot_temp", "dpdx", "dpdy", "pres_hse_x", "pres_hse_y", "z_phys", "detJ" , "mapfac", // eddy viscosity diff --git a/Source/ERF_Constants.H b/Source/ERF_Constants.H index d39ea9f0a..ee00abeb2 100644 --- a/Source/ERF_Constants.H +++ b/Source/ERF_Constants.H @@ -11,6 +11,7 @@ constexpr amrex::Real R_d = 287.0; // dry air constant for dry air [ constexpr amrex::Real R_v = 461.505; // water vapor constant for water vapor [J/(kg-K)] constexpr amrex::Real Cp_d = 1004.5; // We have set this so that with qv=0 we get identically gamma = 1.4 constexpr amrex::Real Cp_v = 1859.0; +constexpr amrex::Real Cp_l = 4200.0; constexpr amrex::Real L_v = 2.5e6; // latent heat of vaporization (J / kg) diff --git a/Source/IO/Plotfile.cpp b/Source/IO/Plotfile.cpp index b8327e13d..b30333a29 100644 --- a/Source/IO/Plotfile.cpp +++ b/Source/IO/Plotfile.cpp @@ -268,6 +268,30 @@ ERF::WritePlotFile (int which, Vector plot_var_names) mf_comp ++; } + if (containerHasElement(plot_var_names, "eq_pot_temp")) + { +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + const Array4& derdat = mf[lev].array(mfi); + const Array4& S_arr = vars_new[lev][Vars::cons].const_array(mfi); + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + Real qv = (use_moisture) ? S_arr(i,j,k,RhoQ1_comp)/S_arr(i,j,k,Rho_comp) : 0.0; + Real qc = (use_moisture) ? S_arr(i,j,k,RhoQ2_comp)/S_arr(i,j,k,Rho_comp) : 0.0; + Real T = getTgivenRandRTh(S_arr(i,j,k,Rho_comp), S_arr(i,j,k,RhoTheta_comp), qv); + Real pressure = getPgivenRTh(S_arr(i,j,k,RhoTheta_comp), qv); + Real fac = Cp_d + Cp_l*(qv + qc); + Real pv = erf_esatw(T)*100.0; + + derdat(i, j, k, mf_comp) = T*std::pow((pressure - pv)/p_0, -R_d/fac)*std::exp(L_v*qv/(fac*T)) ; + }); + } + mf_comp ++; + } + int klo = geom[lev].Domain().smallEnd(2); int khi = geom[lev].Domain().bigEnd(2); diff --git a/Source/TimeIntegration/TI_slow_rhs_fun.H b/Source/TimeIntegration/TI_slow_rhs_fun.H index 9b20fa5dd..79434e28a 100644 --- a/Source/TimeIntegration/TI_slow_rhs_fun.H +++ b/Source/TimeIntegration/TI_slow_rhs_fun.H @@ -271,8 +271,8 @@ // will be used to advance to Real slow_dt = new_stage_time - old_step_time; - bool moist_zero = false; #if defined(ERF_USE_NETCDF) + bool moist_zero = false; if (solverChoice.moisture_type != MoistureType::None && level==0) { // Flag for moisture relaxation zone if (init_type=="real" && wrfbdy_set_width > 0) moist_zero = true; diff --git a/Source/Utils/Microphysics_Utils.H b/Source/Utils/Microphysics_Utils.H index 278f45323..4762de95a 100644 --- a/Source/Utils/Microphysics_Utils.H +++ b/Source/Utils/Microphysics_Utils.H @@ -9,6 +9,7 @@ #include #include #include +#include AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real erf_gammafff (amrex::Real x){ @@ -114,24 +115,24 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void erf_qsati (amrex::Real t, amrex::Real p, amrex::Real &qsati) { amrex::Real esati; esati = erf_esati(t); - qsati = 0.622*esati/std::max(esati,p-esati); + qsati = Rd_on_Rv*esati/std::max(esati,p-esati); } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void erf_qsatw (amrex::Real t, amrex::Real p, amrex::Real &qsatw) { amrex::Real esatw; esatw = erf_esatw(t); - qsatw = 0.622*esatw/std::max(esatw,p-esatw); + qsatw = Rd_on_Rv*esatw/std::max(esatw,p-esatw); } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void erf_dtqsati (amrex::Real t, amrex::Real p, amrex::Real &dtqsati) { - dtqsati = 0.622*erf_dtesati(t)/p; + dtqsati = Rd_on_Rv*erf_dtesati(t)/p; } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void erf_dtqsatw (amrex::Real t, amrex::Real p, amrex::Real &dtqsatw) { - dtqsatw = 0.622*erf_dtesatw(t)/p; + dtqsatw = Rd_on_Rv*erf_dtesatw(t)/p; } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE