diff --git a/Docs/sphinx_doc/Plotfiles.rst b/Docs/sphinx_doc/Plotfiles.rst index 4626358b2..59d18c183 100644 --- a/Docs/sphinx_doc/Plotfiles.rst +++ b/Docs/sphinx_doc/Plotfiles.rst @@ -209,6 +209,22 @@ Output Options | | | | | | +-----------------------------+------------------+ +| **Kmv** | Vertical | +| | Eddy Diffusivity | +| | of Momentum | ++-----------------------------+------------------+ +| **Kmh** | Horizontal | +| | Eddy Diffusivity | +| | of Momentum | ++-----------------------------+------------------+ +| **Khv** | Vertical | +| | Eddy Diffusivity | +| | of Heat | ++-----------------------------+------------------+ +| **Khh** | Horizontal | +| | Eddy Diffusivity | +| | of Heat | ++-----------------------------+------------------+ | **qt** | Nonprecipitating | | | water (qv + qc + | | | qi) | diff --git a/Exec/CMakeLists.txt b/Exec/CMakeLists.txt index 38fe84613..6474d52a0 100644 --- a/Exec/CMakeLists.txt +++ b/Exec/CMakeLists.txt @@ -14,6 +14,7 @@ else () add_subdirectory(RegTests/EkmanSpiral_custom) add_subdirectory(RegTests/EkmanSpiral_ideal) add_subdirectory(RegTests/EkmanSpiral_input_sounding) + add_subdirectory(RegTests/GABLS1) add_subdirectory(RegTests/IsentropicVortex) add_subdirectory(RegTests/PoiseuilleFlow) add_subdirectory(RegTests/ScalarAdvDiff) diff --git a/Exec/RegTests/GABLS1/CMakeLists.txt b/Exec/RegTests/GABLS1/CMakeLists.txt new file mode 100644 index 000000000..6d0ec2792 --- /dev/null +++ b/Exec/RegTests/GABLS1/CMakeLists.txt @@ -0,0 +1,12 @@ +set(erf_exe_name erf_gabls1) + +add_executable(${erf_exe_name} "") +target_sources(${erf_exe_name} + PRIVATE + prob.cpp +) + +target_include_directories(${erf_exe_name} PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}) + +include(${CMAKE_SOURCE_DIR}/CMake/BuildERFExe.cmake) +build_erf_exe(${erf_exe_name}) diff --git a/Exec/RegTests/GABLS1/GNUmakefile b/Exec/RegTests/GABLS1/GNUmakefile new file mode 100644 index 000000000..dd835a830 --- /dev/null +++ b/Exec/RegTests/GABLS1/GNUmakefile @@ -0,0 +1,37 @@ +# AMReX +COMP = gnu +PRECISION = DOUBLE + +# Profiling +PROFILE = FALSE +TINY_PROFILE = TRUE +COMM_PROFILE = FALSE +TRACE_PROFILE = FALSE +MEM_PROFILE = FALSE +USE_GPROF = FALSE + +# Performance +USE_MPI = TRUE +USE_OMP = FALSE + +USE_CUDA = FALSE +USE_HIP = FALSE +USE_SYCL = FALSE + +# Physics +USE_MOISTURE = FALSE + +# Debugging +DEBUG = FALSE + +TEST = TRUE +USE_ASSERTION = TRUE + +#USE_POISSON_SOLVE = TRUE + +# GNU Make +Bpack := ./Make.package +Blocs := . +ERF_HOME := ../.. +ERF_PROBLEM_DIR = $(ERF_HOME)/Exec/RegTests/GABLS1 +include $(ERF_HOME)/Exec/Make.ERF diff --git a/Exec/RegTests/GABLS1/Make.package b/Exec/RegTests/GABLS1/Make.package new file mode 100644 index 000000000..aeacb72f0 --- /dev/null +++ b/Exec/RegTests/GABLS1/Make.package @@ -0,0 +1,2 @@ +CEXE_headers += prob.H +CEXE_sources += prob.cpp diff --git a/Exec/RegTests/GABLS1/input_sounding_GABLS1 b/Exec/RegTests/GABLS1/input_sounding_GABLS1 new file mode 100644 index 000000000..fc02ca975 --- /dev/null +++ b/Exec/RegTests/GABLS1/input_sounding_GABLS1 @@ -0,0 +1,4 @@ +1013.2 265.0 0.0 + 0.0 265.0 0.0 8.0 0.0 + 100.0 265.0 0.0 8.0 0.0 + 400.0 268.0 0.0 8.0 0.0 diff --git a/Exec/RegTests/GABLS1/inputs_gabls1_mynn25_smag b/Exec/RegTests/GABLS1/inputs_gabls1_mynn25_smag new file mode 100644 index 000000000..f31d16159 --- /dev/null +++ b/Exec/RegTests/GABLS1/inputs_gabls1_mynn25_smag @@ -0,0 +1,77 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +stop_time = 32400.0 # 540 min = 9 h (Cuxart et al. 2006) + +amrex.fpe_trap_invalid = 1 + +fabarray.mfiter_tile_size = 1024 1024 1024 + +# PROBLEM SIZE & GEOMETRY (Cuxart et al. 2006) +geometry.prob_extent = 25 25 400 +amr.n_cell = 4 4 64 + +geometry.is_periodic = 1 1 0 + +# MOST BOUNDARY (DEFAULT IS ADIABATIC FOR THETA) +zlo.type = "Most" +erf.most.z0 = 0.1 # from Cuxart et al. 2006 +erf.most.zref = 3.125 # == dz/2 +erf.most.surf_temp = 265.0 # initial value, should match input_sounding +erf.most.surf_heating_rate = -0.25 # [K/h] from Cuxart et al. 2006 + +zhi.type = "SlipWall" +zhi.theta_grad = 0.01 # [K/m] to match the input sounding + +# INITIALIZATION (Cuxart et al. 2006) +erf.init_type = "input_sounding" +erf.init_sounding_ideal = 1 +erf.input_sounding_file = "input_sounding_GABLS1" + +# TIME STEP CONTROL +#erf.fixed_dt = 10.0 # fixed time step (Cuxart et al. 2006) +#erf.no_substepping = 1 +erf.fixed_dt = 1.0 # largest stable low Mach dt +erf.fixed_mri_dt_ratio = 6 + +# 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 = 600 # number of timesteps between checkpoints + +# PLOTFILES +erf.plot_file_1 = plt # prefix of plotfile name +erf.plot_int_1 = 600 # number of timesteps between plotfiles +erf.plot_vars_1 = density x_velocity y_velocity z_velocity pressure theta rhoQKE Kmv Khv + + +# SOLVER CHOICE +erf.dycore_vert_adv_type = "Upwind_3rd" + +erf.molec_diff_type = "None" +#erf.alpha_T = 0.0 +#erf.alpha_C = 1.0 + +erf.use_gravity = true + +# Coriolis parameter f = 1.39e-4 s^-1 (Cuxart et al. 2006) +erf.use_coriolis = true +erf.latitude = 73.0 +erf.rotational_time_period = 86455.2516813368 + +# Geostrophic wind (Cuxart et al. 2006) +erf.abl_driver_type = "GeostrophicWind" +erf.abl_geo_wind = 8.0 0.0 0.0 + +# Turbulence closure +erf.les_type = "Smagorinsky" +erf.Cs = 0.1 +erf.rho0_trans = 1.3223 # from Cuxart et al. 2006 +erf.theta_ref = 263.5 # from Cuxart et al. 2006 + +erf.pbl_type = "MYNN2.5" diff --git a/Exec/RegTests/GABLS1/prob.H b/Exec/RegTests/GABLS1/prob.H new file mode 100644 index 000000000..e4a8ee5be --- /dev/null +++ b/Exec/RegTests/GABLS1/prob.H @@ -0,0 +1,86 @@ +#ifndef _PROB_H_ +#define _PROB_H_ + +#include + +#include "AMReX_REAL.H" + +#include "prob_common.H" + +struct ProbParm : ProbParmDefaults { + amrex::Real rho_0 = 0.0; + amrex::Real T_0 = 0.0; + amrex::Real A_0 = 1.0; + amrex::Real QKE_0 = 1e-12; // min value + + amrex::Real U_0 = 0.0; + amrex::Real V_0 = 0.0; + amrex::Real W_0 = 0.0; + + // random initial perturbations (legacy code) + amrex::Real U_0_Pert_Mag = 0.0; + amrex::Real V_0_Pert_Mag = 0.0; + amrex::Real W_0_Pert_Mag = 0.0; + amrex::Real T_0_Pert_Mag = 0.0; // perturbation to rho*Theta + + // divergence-free initial perturbations + amrex::Real pert_deltaU = 0.0; + amrex::Real pert_deltaV = 0.0; + amrex::Real pert_periods_U = 5.0; + amrex::Real pert_periods_V = 5.0; + amrex::Real pert_ref_height = 100.0; + + // rayleigh damping + amrex::Real dampcoef = 0.2; // inverse time scale [1/s] + amrex::Real zdamp = 500.0; // damping depth [m] from model top + + // helper vars + amrex::Real aval; + amrex::Real bval; + amrex::Real ufac; + amrex::Real vfac; +}; // namespace ProbParm + +class Problem : public ProblemBase +{ +public: + Problem(const amrex::Real* problo, const amrex::Real* probhi); + +#include "Prob/init_constant_density_hse.H" +#include "Prob/init_rayleigh_damping.H" + + void init_custom_pert ( + const amrex::Box& bx, + const amrex::Box& xbx, + const amrex::Box& ybx, + const amrex::Box& zbx, + amrex::Array4 const& state, + amrex::Array4 const& x_vel, + amrex::Array4 const& y_vel, + amrex::Array4 const& z_vel, + amrex::Array4 const& r_hse, + amrex::Array4 const& p_hse, + amrex::Array4 const& z_nd, + amrex::Array4 const& z_cc, +#if defined(ERF_USE_MOISTURE) + amrex::Array4 const& qv, + amrex::Array4 const& qc, + amrex::Array4 const& qi, +#elif defined(ERF_USE_WARM_NO_PRECIP) + amrex::Array4 const& qv, + amrex::Array4 const& qc, +#endif + amrex::GeometryData const& geomdata, + amrex::Array4 const& mf_m, + amrex::Array4 const& mf_u, + amrex::Array4 const& mf_v, + const SolverChoice& sc) override; + +protected: + std::string name() override { return "ABL"; } + +private: + ProbParm parms; +}; + +#endif diff --git a/Exec/RegTests/GABLS1/prob.cpp b/Exec/RegTests/GABLS1/prob.cpp new file mode 100644 index 000000000..c9bffb30f --- /dev/null +++ b/Exec/RegTests/GABLS1/prob.cpp @@ -0,0 +1,181 @@ +#include "prob.H" +#include "AMReX_Random.H" + +using namespace amrex; + +std::unique_ptr +amrex_probinit(const amrex_real* problo, const amrex_real* probhi) +{ + return std::make_unique(problo, probhi); +} + +Problem::Problem(const amrex::Real* problo, const amrex::Real* probhi) +{ + // Parse params + ParmParse pp("prob"); + pp.query("rho_0", parms.rho_0); + pp.query("T_0", parms.T_0); + pp.query("A_0", parms.A_0); + pp.query("QKE_0", parms.QKE_0); + + pp.query("U_0", parms.U_0); + pp.query("V_0", parms.V_0); + pp.query("W_0", parms.W_0); + pp.query("U_0_Pert_Mag", parms.U_0_Pert_Mag); + pp.query("V_0_Pert_Mag", parms.V_0_Pert_Mag); + pp.query("W_0_Pert_Mag", parms.W_0_Pert_Mag); + pp.query("T_0_Pert_Mag", parms.T_0_Pert_Mag); + + pp.query("pert_deltaU", parms.pert_deltaU); + pp.query("pert_deltaV", parms.pert_deltaV); + pp.query("pert_periods_U", parms.pert_periods_U); + pp.query("pert_periods_V", parms.pert_periods_V); + pp.query("pert_ref_height", parms.pert_ref_height); + parms.aval = parms.pert_periods_U * 2.0 * PI / (probhi[1] - problo[1]); + parms.bval = parms.pert_periods_V * 2.0 * PI / (probhi[0] - problo[0]); + parms.ufac = parms.pert_deltaU * std::exp(0.5) / parms.pert_ref_height; + parms.vfac = parms.pert_deltaV * std::exp(0.5) / parms.pert_ref_height; + + pp.query("dampcoef", parms.dampcoef); + pp.query("zdamp", parms.zdamp); + + init_base_parms(parms.rho_0, parms.T_0); +} + +void +Problem::init_custom_pert( + const amrex::Box& bx, + const amrex::Box& xbx, + const amrex::Box& ybx, + const amrex::Box& zbx, + amrex::Array4 const& state, + amrex::Array4 const& x_vel, + amrex::Array4 const& y_vel, + amrex::Array4 const& z_vel, + amrex::Array4 const& r_hse, + amrex::Array4 const& /*p_hse*/, + amrex::Array4 const& /*z_nd*/, + amrex::Array4 const& /*z_cc*/, +#if defined(ERF_USE_MOISTURE) + amrex::Array4 const& /*qv*/, + amrex::Array4 const& /*qc*/, + amrex::Array4 const& /*qi*/, +#elif defined(ERF_USE_WARM_NO_PRECIP) + amrex::Array4 const& /*qv*/, + amrex::Array4 const& /*qc*/, +#endif + amrex::GeometryData const& geomdata, + amrex::Array4 const& /*mf_m*/, + amrex::Array4 const& /*mf_u*/, + amrex::Array4 const& /*mf_v*/, + const SolverChoice& /*sc*/) +{ + ParallelForRNG(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept { + // Geometry + const Real* prob_lo = geomdata.ProbLo(); + const Real* prob_hi = geomdata.ProbHi(); + const Real* dx = geomdata.CellSize(); + 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]; + + // Define a point (xc,yc,zc) at the center of the domain + const Real xc = 0.5 * (prob_lo[0] + prob_hi[0]); + const Real yc = 0.5 * (prob_lo[1] + prob_hi[1]); + const Real zc = 0.5 * (prob_lo[2] + prob_hi[2]); + + const Real r = std::sqrt((x-xc)*(x-xc) + (y-yc)*(y-yc) + (z-zc)*(z-zc)); + + // Add temperature perturbations + if ((z <= parms.pert_ref_height) && (parms.T_0_Pert_Mag != 0.0)) { + Real rand_double = amrex::Random(engine); // Between 0.0 and 1.0 + state(i, j, k, RhoTheta_comp) = (rand_double*2.0 - 1.0)*parms.T_0_Pert_Mag; + } + + // Set scalar = A_0*exp(-10r^2), where r is distance from center of domain + state(i, j, k, RhoScalar_comp) = parms.A_0 * exp(-10.*r*r); + + // Set an initial value for QKE = 2*e + if (z < 250) { + // Following Cuxart et al. 2006 + state(i, j, k, RhoQKE_comp) = r_hse(i,j,k) * 0.8 * std::pow(1 - z/250, 3); + } else { + state(i, j, k, RhoQKE_comp) = r_hse(i,j,k) * parms.QKE_0; + } + +#if defined(ERF_USE_MOISTURE) + state(i, j, k, RhoQt_comp) = 0.0; + state(i, j, k, RhoQp_comp) = 0.0; +#elif defined(ERF_USE_WARM_NO_PRECIP) + state(i, j, k, RhoQv_comp) = 0.0; + state(i, j, k, RhoQc_comp) = 0.0; +#endif + }); + + // Set the x-velocity + ParallelForRNG(xbx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept { + const Real* prob_lo = geomdata.ProbLo(); + const Real* dx = geomdata.CellSize(); + const Real y = prob_lo[1] + (j + 0.5) * dx[1]; + const Real z = prob_lo[2] + (k + 0.5) * dx[2]; + + // Set the x-velocity + x_vel(i, j, k) = parms.U_0; + if ((z <= parms.pert_ref_height) && (parms.U_0_Pert_Mag != 0.0)) + { + Real rand_double = amrex::Random(engine); // Between 0.0 and 1.0 + Real x_vel_prime = (rand_double*2.0 - 1.0)*parms.U_0_Pert_Mag; + x_vel(i, j, k) += x_vel_prime; + } + if (parms.pert_deltaU != 0.0) + { + const amrex::Real yl = y - prob_lo[1]; + const amrex::Real zl = z / parms.pert_ref_height; + const amrex::Real damp = std::exp(-0.5 * zl * zl); + x_vel(i, j, k) += parms.ufac * damp * z * std::cos(parms.aval * yl); + } + }); + + // Set the y-velocity + ParallelForRNG(ybx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept { + const Real* prob_lo = geomdata.ProbLo(); + const Real* dx = geomdata.CellSize(); + const Real x = prob_lo[0] + (i + 0.5) * dx[0]; + const Real z = prob_lo[2] + (k + 0.5) * dx[2]; + + // Set the y-velocity + y_vel(i, j, k) = parms.V_0; + if ((z <= parms.pert_ref_height) && (parms.V_0_Pert_Mag != 0.0)) + { + Real rand_double = amrex::Random(engine); // Between 0.0 and 1.0 + Real y_vel_prime = (rand_double*2.0 - 1.0)*parms.V_0_Pert_Mag; + y_vel(i, j, k) += y_vel_prime; + } + if (parms.pert_deltaV != 0.0) + { + const amrex::Real xl = x - prob_lo[0]; + const amrex::Real zl = z / parms.pert_ref_height; + const amrex::Real damp = std::exp(-0.5 * zl * zl); + y_vel(i, j, k) += parms.vfac * damp * z * std::cos(parms.bval * xl); + } + }); + + // Set the z-velocity + ParallelForRNG(zbx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept { + const int dom_lo_z = geomdata.Domain().smallEnd()[2]; + const int dom_hi_z = geomdata.Domain().bigEnd()[2]; + + // Set the z-velocity + if (k == dom_lo_z || k == dom_hi_z+1) + { + z_vel(i, j, k) = 0.0; + } + else if (parms.W_0_Pert_Mag != 0.0) + { + Real rand_double = amrex::Random(engine); // Between 0.0 and 1.0 + Real z_vel_prime = (rand_double*2.0 - 1.0)*parms.W_0_Pert_Mag; + z_vel(i, j, k) = parms.W_0 + z_vel_prime; + } + }); +} + diff --git a/Source/BoundaryConditions/ABLMost.H b/Source/BoundaryConditions/ABLMost.H index 73d3f26ea..9d6fd6bac 100644 --- a/Source/BoundaryConditions/ABLMost.H +++ b/Source/BoundaryConditions/ABLMost.H @@ -54,8 +54,16 @@ public: auto erf_st = pp.query("most.surf_temp", surf_temp); if (erf_st) { theta_type = ThetaCalcType::SURFACE_TEMPERATURE; + pp.query("most.surf_heating_rate", surf_heating_rate); // [K/h] + surf_heating_rate = surf_heating_rate / 3600.0; // [K/s] + if (pp.query("most.surf_temp_flux", surf_temp_flux)) { + amrex::Abort("Can only specify one of surf_temp_flux or surf_heating_rate"); + } } else { pp.query("most.surf_temp_flux", surf_temp_flux); + if (pp.query("most.surf_heating_rate", surf_heating_rate)) { + amrex::Abort("Can only specify one of surf_temp_flux or surf_heating_rate"); + } if (std::abs(surf_temp_flux) > std::numeric_limits::epsilon()) { theta_type = ThetaCalcType::HEAT_FLUX; } else { @@ -129,6 +137,22 @@ public: }// lev } + void + update_surf_temp (amrex::Real cur_time) + { + if (surf_heating_rate != 0) { + int nlevs = m_geom.size(); + for (int lev = 0; lev < nlevs; lev++) + { + t_surf[lev]->setVal(surf_temp + surf_heating_rate*cur_time); + amrex::Print() << "Surface temp at t=" << cur_time + << ": " + << surf_temp + surf_heating_rate*cur_time + << std::endl; + } + } + } + void impose_most_bcs (const int lev, const amrex::Vector& mfs, @@ -145,7 +169,9 @@ public: const FluxCalc& flux_comp); void - update_fluxes (int lev, int max_iters = 25); + update_fluxes (int lev, + amrex::Real cur_time, + int max_iters = 25); template void @@ -195,6 +221,7 @@ public: private: amrex::Real z0_const; amrex::Real surf_temp; + amrex::Real surf_heating_rate{0}; amrex::Real surf_temp_flux{0}; amrex::Real cnk_a{0.0185}; amrex::Real depth{30.0}; diff --git a/Source/BoundaryConditions/ABLMost.cpp b/Source/BoundaryConditions/ABLMost.cpp index 96bbadc01..c7daf5bbe 100644 --- a/Source/BoundaryConditions/ABLMost.cpp +++ b/Source/BoundaryConditions/ABLMost.cpp @@ -11,6 +11,7 @@ using namespace amrex; */ void ABLMost::update_fluxes (int lev, + Real cur_time, int max_iters) { // Compute plane averages for all vars (regardless of flux type) @@ -30,6 +31,7 @@ ABLMost::update_fluxes (int lev, compute_fluxes(lev, max_iters, most_flux); } } else if (theta_type == ThetaCalcType::SURFACE_TEMPERATURE) { + update_surf_temp(cur_time); if (rough_type == RoughCalcType::CONSTANT) { surface_temp most_flux(m_ma.get_zref(), surf_temp_flux); compute_fluxes(lev, max_iters, most_flux); diff --git a/Source/BoundaryConditions/MOSTStress.H b/Source/BoundaryConditions/MOSTStress.H index d6661705c..1e4204192 100644 --- a/Source/BoundaryConditions/MOSTStress.H +++ b/Source/BoundaryConditions/MOSTStress.H @@ -649,7 +649,8 @@ struct moeng_flux amrex::Real wsp = sqrt(velx*velx+vely*vely); amrex::Real num1 = (theta-theta_mean)*wsp_mean; amrex::Real num2 = (theta_mean-theta_surf)*wsp; - amrex::Real moflux = tstar*ustar*(num1+num2)/((theta_mean-theta_surf)*wsp_mean); + amrex::Real moflux = (std::abs(tstar) > EPS) ? + tstar*ustar*(num1+num2)/((theta_mean-theta_surf)*wsp_mean) : 0.0; amrex::Real deltaz = dz * (zlo - k); dest_arr(i,j,k,icomp+n) = rho*(theta - moflux*rho/eta*deltaz); @@ -783,6 +784,7 @@ struct moeng_flux private: int zlo; + static constexpr amrex::Real EPS = 1e-16; }; diff --git a/Source/DataStructs/DataStruct.H b/Source/DataStructs/DataStruct.H index 6a2abcb1d..e8941c11e 100644 --- a/Source/DataStructs/DataStruct.H +++ b/Source/DataStructs/DataStruct.H @@ -124,8 +124,14 @@ struct SolverChoice { for(int i = 0; i < AMREX_SPACEDIM; ++i) abl_pressure_grad[i] = abl_pressure_grad_in[i]; amrex::Vector abl_geo_forcing_in = {0.0, 0.0, 0.0}; - pp.queryarr("abl_geo_forcing",abl_geo_forcing_in); - for(int i = 0; i < AMREX_SPACEDIM; ++i) abl_geo_forcing[i] = abl_geo_forcing_in[i]; + if(pp.queryarr("abl_geo_forcing",abl_geo_forcing_in)) { + amrex::Print() << "Specified abl_geo_forcing: ("; + for (int i = 0; i < AMREX_SPACEDIM; ++i) { + abl_geo_forcing[i] = abl_geo_forcing_in[i]; + amrex::Print() << abl_geo_forcing[i] << " "; + } + amrex::Print() << ")" << std::endl; + } if (use_coriolis) { diff --git a/Source/Diffusion/PBLModels.cpp b/Source/Diffusion/PBLModels.cpp index 1e7ac7b6a..c890a59d4 100644 --- a/Source/Diffusion/PBLModels.cpp +++ b/Source/Diffusion/PBLModels.cpp @@ -227,7 +227,7 @@ ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel, amrex::Real SM = (R2*E2 - R1*E4)/(E2*E3 - E1*E4); amrex::Real SH = (R1*E3 - R2*E1)/(E2*E3 - E1*E4); - amrex::Real SQ = 3.0 * SM; + amrex::Real SQ = 3.0 * SM; // Nakanishi & Niino 2009 // Finally, compute the eddy viscosity/diffusivities const amrex::Real rho = cell_data(i,j,k,Rho_comp); diff --git a/Source/ERF.H b/Source/ERF.H index 456a7a8f5..55a3a11b6 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -653,7 +653,12 @@ private: // Note that the order of variable names here must match the order in Derive.cpp const amrex::Vector derived_names {"pressure", "soundspeed", "temp", "theta", "KE", "QKE", "scalar", "pres_hse", "dens_hse", "pert_pres", "pert_dens", - "dpdx", "dpdy", "pres_hse_x", "pres_hse_y", "z_phys", "detJ" , "mapfac" + "dpdx", "dpdy", "pres_hse_x", "pres_hse_y", + "z_phys", "detJ" , "mapfac", + // eddy viscosity + "Kmv","Kmh", + // eddy diffusivity of heat + "Khv","Khh" #if defined(ERF_USE_MOISTURE) ,"qt", "qp", "qv", "qc", "qi", "qrain", "qsnow", "qgraup" #elif defined(ERF_USE_WARM_NO_PRECIP) diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 7dc8b4440..ac2641b92 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -467,6 +467,24 @@ ERF::InitData () } // lev } + Real QKE_0; + if (pp.query("QKE_0", QKE_0)) { + for (int lev = 0; lev <= finest_level; lev++) { + auto& lev_new = vars_new[lev]; + for (MFIter mfi(lev_new[Vars::cons], TilingIfNotGPU()); mfi.isValid(); ++mfi) { + const Box &bx = mfi.tilebox(); + const auto &cons_arr = lev_new[Vars::cons].array(mfi); + // We want to set the lateral BC values, too + Box gbx = bx; // Copy constructor + gbx.grow(0,1); gbx.grow(1,1); // Grow by one in the lateral directions + amrex::ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + cons_arr(i,j,k,RhoQKE_comp) = cons_arr(i,j,k,Rho_comp) * QKE_0; + }); + } // mfi + } + } + + if (coupling_type == "TwoWay") { AverageDown(); } @@ -548,7 +566,7 @@ ERF::InitData () MultiFab::Copy( *Theta_prim[lev], S, Cons::RhoTheta, 0, 1, ng); MultiFab::Divide(*Theta_prim[lev], S, Cons::Rho , 0, 1, ng); m_most->update_mac_ptrs(lev, vars_new, Theta_prim); - m_most->update_fluxes(lev); + m_most->update_fluxes(lev,t_new[lev]); } } diff --git a/Source/IO/Plotfile.cpp b/Source/IO/Plotfile.cpp index 144599879..5763a585a 100644 --- a/Source/IO/Plotfile.cpp +++ b/Source/IO/Plotfile.cpp @@ -544,6 +544,23 @@ ERF::WritePlotFile (int which, Vector plot_var_names) mf_comp ++; } + if (containerHasElement(plot_var_names, "Kmv")) { + MultiFab::Copy(mf[lev],*eddyDiffs_lev[lev],EddyDiff::Mom_v,mf_comp,1,0); + mf_comp ++; + } + if (containerHasElement(plot_var_names, "Kmh")) { + MultiFab::Copy(mf[lev],*eddyDiffs_lev[lev],EddyDiff::Mom_h,mf_comp,1,0); + mf_comp ++; + } + if (containerHasElement(plot_var_names, "Khv")) { + MultiFab::Copy(mf[lev],*eddyDiffs_lev[lev],EddyDiff::Theta_v,mf_comp,1,0); + mf_comp ++; + } + if (containerHasElement(plot_var_names, "Khh")) { + MultiFab::Copy(mf[lev],*eddyDiffs_lev[lev],EddyDiff::Theta_h,mf_comp,1,0); + mf_comp ++; + } + #if defined(ERF_USE_MOISTURE) calculate_derived("qt", derived::erf_derQt); calculate_derived("qp", derived::erf_derQp); diff --git a/Source/Initialization/InputSoundingData.H b/Source/Initialization/InputSoundingData.H index 046567b4a..3e292409f 100644 --- a/Source/Initialization/InputSoundingData.H +++ b/Source/Initialization/InputSoundingData.H @@ -129,10 +129,9 @@ public: * density and pressure, then back down to get the dry density and * pressure. * - * This deviates from that code slightly because we calculate the - * correct moist theta, equal to virtual potential temperature, and - * also remove the factor (1+qv) for consistency with the the surface - * density. + * This deviates from that code slightly because we calculate a moist + * theta that is equal to virtual potential temperature, and also + * remove the factor of (1+qv) for consistency with the surface density. */ const int maxiter = 10; const int Ninp = size(); @@ -158,6 +157,18 @@ public: // integrate from surface to domain top amrex::Real qvf, dz; +#if 0 + // In this absence of moisture, this moist profile will match the + // following dry profile + amrex::Print() << "z p_m rho_m theta U V" << std::endl; + amrex::Print() << z_inp_sound[0] + << " " << pm_integ[0] + << " " << rhom_integ[0] + << " " << theta_inp_sound[0] + << " " << U_inp_sound[0] + << " " << V_inp_sound[0] + << std::endl; +#endif for (int k=1; k < size(); ++k) { qvf = 1 + (R_v/R_d-1) * qv_inp_sound[k]; @@ -171,9 +182,18 @@ public: rhom_integ[k] = 1 / ( R_d/p_0 * theta_inp_sound[k]*qvf * std::pow(pm_integ[k]/p_0, -iGamma)); } +#if 0 + amrex::Print() << z_inp_sound[k] + << " " << pm_integ[k] + << " " << rhom_integ[k] + << " " << theta_inp_sound[k] + << " " << U_inp_sound[k] + << " " << V_inp_sound[k] + << std::endl; +#endif } - // now, integrate from the top of the sounding (where it's dry) back + // now, integrate from the top of the sounding (where it's assumed dry) back // down to get the dry air column properties pd_integ[Ninp-1] = pm_integ[Ninp-1]; rhod_integ[Ninp-1] = 1 / ( diff --git a/Source/Prob/init_rayleigh_damping.H b/Source/Prob/init_rayleigh_damping.H index f0c550d51..b10d4271c 100644 --- a/Source/Prob/init_rayleigh_damping.H +++ b/Source/Prob/init_rayleigh_damping.H @@ -11,8 +11,9 @@ erf_init_rayleigh (amrex::Vector& tau, amrex::Geometry const& geom) override { const auto ztop = geom.ProbHi()[2]; - amrex::Print() << "Rayleigh damping (coef="<update_mac_ptrs(lev, vars_old, Theta_prim); - m_most->update_fluxes(lev); + m_most->update_fluxes(lev, time); } } diff --git a/Source/TimeIntegration/ERF_slow_rhs_post.cpp b/Source/TimeIntegration/ERF_slow_rhs_post.cpp index 62f7e2020..32579dbfb 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_post.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_post.cpp @@ -441,7 +441,8 @@ void erf_slow_rhs_post (int level, // NOTE: we don't include additional source terms when terrain is moving Real temp_val = detJ_arr(i,j,k) * old_cons(i,j,k,n) + dt * detJ_arr(i,j,k) * cell_rhs(i,j,k,n); cur_cons(i,j,k,n) = temp_val / detJ_new_arr(i,j,k); -#ifndef AMREX_USE_GPU + cur_cons(i,j,k,n) = amrex::max(cur_cons(i,j,k,n), 1e-12); +#if 0 if (cur_cons(i,j,k,n) < Real(0.)) { amrex::AllPrint() << "MAKING NEGATIVE QKE " << IntVect(i,j,k) << " NEW / OLD " << cur_cons(i,j,k,n) << " " << old_cons(i,j,k,n) << std::endl; @@ -482,7 +483,8 @@ void erf_slow_rhs_post (int level, const int n = start_comp + nn; cell_rhs(i,j,k,n) += src_arr(i,j,k,n); cur_cons(i,j,k,n) = old_cons(i,j,k,n) + dt * cell_rhs(i,j,k,n); -#ifndef AMREX_USE_GPU + cur_cons(i,j,k,n) = amrex::max(cur_cons(i,j,k,n), 1e-12); +#if 0 if (cur_cons(i,j,k,n) < Real(0.)) { amrex::AllPrint() << "MAKING NEGATIVE QKE " << IntVect(i,j,k) << " NEW / OLD " << cur_cons(i,j,k,n) << " " << old_cons(i,j,k,n) << std::endl;