From 0e5e4409407c7b90e696cd43c5ecc5e7ebc1c88a Mon Sep 17 00:00:00 2001 From: "Aaron M. Lattanzi" <103702284+AMLattanzi@users.noreply.github.com> Date: Tue, 28 Nov 2023 16:33:13 -0800 Subject: [PATCH] Kessler Microphysics (#1324) * Changes to EOS.H with moisture * Change all functions in EOS.H to work with moisture * Adding qp for buoyancy type 1 * Updating docs for buoyancy term * Updating docs for buoyancy term * Update docs for buoyancy * Avoiding extra divide * Update buoyancy docs * Some corrections to buoyancy docs * Some corrections to Docs in buoyancy * Correcting docs issue * Correcting docs * adding Kessler microphysics docs * adding Kessler microphysics docs * adding Kessler microphysics docs * adding Kessler microphysics docs * WIP: Adding Kessler microphysics * Guards for microphysics * WIP: Adding Kessler microphysics * WIP: kessler microphysics * Compiled and verified initial condition for Kessler * theta pert to 0 for hydrostatic case * Calling moist initial conditions * Call moist for moisture cases * Pert to density and theta * Avoiding name clashes * Complete Kessler microphysics draft * Updating squall line test case * Moving towards Newton Raphson for IC * Updated Kessler microphysics for push * Updating inputs * Adding inputs with outflow bcs * Updating periodic inputs * Removing old prob.cpp * Minor changes and removing whitespaces and tabs * Updating periodic inputs * Finalizing squall line Kessler * Updating inputs for Kessler * Updating squall line test case * enforce qv >= 0.0 * Comment out unused vars in problemBase. * Initi1D an only access use_moist_bkgd when compiled with moisture support. * Fix HIP issues. * Add Kessler to CMake compile. * Temporary patch, doesn't fix device capture. * Not ideal but compiled with nvcc. * Use PI instead of M_PI for windows. --------- Co-authored-by: Mahesh Natarajan Co-authored-by: Ann Almgren Co-authored-by: Aaron Lattanzi --- CMake/BuildERFExe.cmake | 7 +- Exec/CMakeLists.txt | 1 + Exec/Make.ERF | 5 + Exec/SquallLine_2D/CMakeLists.txt | 12 + Exec/SquallLine_2D/GNUmakefile | 35 ++ Exec/SquallLine_2D/Make.package | 2 + Exec/SquallLine_2D/README | 2 + Exec/SquallLine_2D/inputs_moisture | 84 +++ Exec/SquallLine_2D/inputs_moisture_outflow | 86 ++++ Exec/SquallLine_2D/prob.H | 89 ++++ Exec/SquallLine_2D/prob.cpp | 477 ++++++++++++++++++ Source/DataStructs/DataStruct.H | 6 +- Source/ERF.cpp | 2 + Source/Initialization/ERF_init1d.cpp | 15 +- .../Microphysics/Kessler/Diagnose_Kessler.cpp | 42 ++ Source/Microphysics/Kessler/Init_Kessler.cpp | 249 +++++++++ Source/Microphysics/Kessler/Kessler.H | 159 ++++++ Source/Microphysics/Kessler/Kessler.cpp | 238 +++++++++ Source/Microphysics/Kessler/Make.package | 5 + .../Microphysics/Kessler/Update_Kessler.cpp | 60 +++ Source/Microphysics/Microphysics.H | 1 + Source/Microphysics/Null/NullMoist.H | 10 +- Source/Microphysics/SAM/SAM.H | 6 +- Source/TimeIntegration/ERF_make_buoyancy.cpp | 28 +- Source/prob_common.H | 9 + 25 files changed, 1612 insertions(+), 18 deletions(-) create mode 100644 Exec/SquallLine_2D/CMakeLists.txt create mode 100644 Exec/SquallLine_2D/GNUmakefile create mode 100644 Exec/SquallLine_2D/Make.package create mode 100644 Exec/SquallLine_2D/README create mode 100644 Exec/SquallLine_2D/inputs_moisture create mode 100644 Exec/SquallLine_2D/inputs_moisture_outflow create mode 100644 Exec/SquallLine_2D/prob.H create mode 100644 Exec/SquallLine_2D/prob.cpp create mode 100644 Source/Microphysics/Kessler/Diagnose_Kessler.cpp create mode 100644 Source/Microphysics/Kessler/Init_Kessler.cpp create mode 100644 Source/Microphysics/Kessler/Kessler.H create mode 100644 Source/Microphysics/Kessler/Kessler.cpp create mode 100644 Source/Microphysics/Kessler/Make.package create mode 100644 Source/Microphysics/Kessler/Update_Kessler.cpp diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake index db62d9395..dbf0788ad 100644 --- a/CMake/BuildERFExe.cmake +++ b/CMake/BuildERFExe.cmake @@ -65,7 +65,11 @@ function(build_erf_lib erf_lib_name) ${SRC_DIR}/Microphysics/SAM/Precip.cpp ${SRC_DIR}/Microphysics/SAM/PrecipFall.cpp ${SRC_DIR}/Microphysics/SAM/Diagnose.cpp - ${SRC_DIR}/Microphysics/SAM/Update.cpp) + ${SRC_DIR}/Microphysics/SAM/Update.cpp + ${SRC_DIR}/Microphysics/Kessler/Init_Kessler.cpp + ${SRC_DIR}/Microphysics/Kessler/Kessler.cpp + ${SRC_DIR}/Microphysics/Kessler/Diagnose_Kessler.cpp + ${SRC_DIR}/Microphysics/Kessler/Update_Kessler.cpp) target_compile_definitions(${erf_lib_name} PUBLIC ERF_USE_MOISTURE) endif() @@ -187,6 +191,7 @@ function(build_erf_lib erf_lib_name) target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/Microphysics) target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/Microphysics/Null) target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/Microphysics/SAM) + target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/Microphysics/Kessler) endif() if(ERF_ENABLE_RRTMGP) diff --git a/Exec/CMakeLists.txt b/Exec/CMakeLists.txt index 78c93267c..101a2dfe3 100644 --- a/Exec/CMakeLists.txt +++ b/Exec/CMakeLists.txt @@ -8,6 +8,7 @@ else () add_subdirectory(ABL) add_subdirectory(SuperCell) add_subdirectory(Radiation) + add_subdirectory(SquallLine_2D) add_subdirectory(RegTests/Bubble) add_subdirectory(RegTests/CouetteFlow) add_subdirectory(RegTests/DensityCurrent) diff --git a/Exec/Make.ERF b/Exec/Make.ERF index 0505ccf8d..5d38535b7 100644 --- a/Exec/Make.ERF +++ b/Exec/Make.ERF @@ -117,6 +117,11 @@ ifeq ($(USE_MOISTURE), TRUE) include $(ERF_MOISTURE_SAM_DIR)/Make.package VPATH_LOCATIONS += $(ERF_MOISTURE_SAM_DIR) INCLUDE_LOCATIONS += $(ERF_MOISTURE_SAM_DIR) + + ERF_MOISTURE_KESSLER_DIR = $(ERF_SOURCE_DIR)/Microphysics/Kessler + include $(ERF_MOISTURE_KESSLER_DIR)/Make.package + VPATH_LOCATIONS += $(ERF_MOISTURE_KESSLER_DIR) + INCLUDE_LOCATIONS += $(ERF_MOISTURE_KESSLER_DIR) endif ifeq ($(USE_WARM_NO_PRECIP), TRUE) diff --git a/Exec/SquallLine_2D/CMakeLists.txt b/Exec/SquallLine_2D/CMakeLists.txt new file mode 100644 index 000000000..67e9bf880 --- /dev/null +++ b/Exec/SquallLine_2D/CMakeLists.txt @@ -0,0 +1,12 @@ +set(erf_exe_name squallline_2d) + +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/SquallLine_2D/GNUmakefile b/Exec/SquallLine_2D/GNUmakefile new file mode 100644 index 000000000..40f807cd1 --- /dev/null +++ b/Exec/SquallLine_2D/GNUmakefile @@ -0,0 +1,35 @@ +# AMReX +COMP = gnu +PRECISION = DOUBLE + +# Profiling +PROFILE = FALSE +TINY_PROFILE = FALSE +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 + +# Debugging +DEBUG = FALSE + +TEST = TRUE +USE_ASSERTION = TRUE + +USE_MOISTURE = TRUE +#USE_WARM_NO_PRECIP = TRUE + +# GNU Make +Bpack := ./Make.package +Blocs := . +ERF_HOME := ../.. +ERF_PROBLEM_DIR = $(ERF_HOME)/Exec/SquallLine_2D +include $(ERF_HOME)/Exec/Make.ERF diff --git a/Exec/SquallLine_2D/Make.package b/Exec/SquallLine_2D/Make.package new file mode 100644 index 000000000..aeacb72f0 --- /dev/null +++ b/Exec/SquallLine_2D/Make.package @@ -0,0 +1,2 @@ +CEXE_headers += prob.H +CEXE_sources += prob.cpp diff --git a/Exec/SquallLine_2D/README b/Exec/SquallLine_2D/README new file mode 100644 index 000000000..9242958db --- /dev/null +++ b/Exec/SquallLine_2D/README @@ -0,0 +1,2 @@ +This problem setup is the evolution of a supercell, which primarily tests the ability +of ERF to model moisture physics. diff --git a/Exec/SquallLine_2D/inputs_moisture b/Exec/SquallLine_2D/inputs_moisture new file mode 100644 index 000000000..1a6d7fb48 --- /dev/null +++ b/Exec/SquallLine_2D/inputs_moisture @@ -0,0 +1,84 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +max_step = 9000 +stop_time = 90000.0 + +amrex.fpe_trap_invalid = 1 + +fabarray.mfiter_tile_size = 2048 1024 2048 + +# PROBLEM SIZE & GEOMETRY +geometry.prob_lo = -25000. 0. 0. +geometry.prob_hi = 25000. 400. 20000. +amr.n_cell = 192 4 81 # dx=dy=dz=100 m + +# periodic in x to match WRF setup +# - as an alternative, could use symmetry at x=0 and outflow at x=25600 +geometry.is_periodic = 1 1 0 +#xlo.type = "Outflow" +#xhi.type = "Outflow" +zlo.type = "SlipWall" +zhi.type = "Outflow" + +erf.sponge_strength = 2.0 +#erf.use_zhi_sponge_damping = true +erf.zhi_sponge_start = 12000.0 + +erf.sponge_density = 1.2 +erf.sponge_x_velocity = 0.0 +erf.sponge_y_velocity = 0.0 +erf.sponge_z_velocity = 0.0 + +# TIME STEP CONTROL +erf.use_native_mri = 1 +erf.fixed_dt = 1.0 # fixed time step [s] -- Straka et al 1993 +erf.fixed_fast_dt = 0.5 # fixed time step [s] -- Straka et al 1993 +#erf.no_substepping = 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 +amr.check_file = chk # root name of checkpoint file +amr.check_int = 1000 # number of timesteps between checkpoints +#amr.restart = chk12000 + +# PLOTFILES +erf.plot_file_1 = plt # root name of plotfile +erf.plot_int_1 = 100 # number of timesteps between plotfiles +erf.plot_vars_1 = density rhotheta rhoQt rhoQp x_velocity y_velocity z_velocity pressure theta temp qt qp qv qc qi scalar pert_dens + +# SOLVER CHOICE +erf.use_gravity = true +erf.buoyancy_type = 4 +erf.use_coriolis = false +erf.use_rayleigh_damping = false + +#erf.les_type = "Smagorinsky" +erf.Cs = 0.25 +erf.les_type = "None" + +# +# diffusion coefficient from Straka, K = 75 m^2/s +# +erf.molec_diff_type = "ConstantAlpha" +#erf.molec_diff_type = "Constant" +erf.rho0_trans = 1.0 # [kg/m^3], used to convert input diffusivities +erf.dynamicViscosity = 200.0 # [kg/(m-s)] ==> nu = 75.0 m^2/s +erf.alpha_T = 00.0 # [m^2/s] +erf.alpha_C = 50.0 + +erf.moisture_model = "Kessler" +erf.use_moist_background = true + +erf.moistscal_horiz_adv_string = "Centered_2nd" +erf.moistscal_vert_adv_string = "Centered_2nd" + +# PROBLEM PARAMETERS (optional) +prob.T_0 = 300.0 +prob.U_0 = 0 +prob.T_pert = 3 diff --git a/Exec/SquallLine_2D/inputs_moisture_outflow b/Exec/SquallLine_2D/inputs_moisture_outflow new file mode 100644 index 000000000..0411708df --- /dev/null +++ b/Exec/SquallLine_2D/inputs_moisture_outflow @@ -0,0 +1,86 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +max_step = 9000 +stop_time = 90000.0 + +amrex.fpe_trap_invalid = 1 + +fabarray.mfiter_tile_size = 2048 1024 2048 + +# PROBLEM SIZE & GEOMETRY +geometry.prob_lo = -25000. 0. 0. +geometry.prob_hi = 25000. 400. 20000. +amr.n_cell = 192 4 81 # dx=dy=dz=100 m + +# periodic in x to match WRF setup +# - as an alternative, could use symmetry at x=0 and outflow at x=25600 +geometry.is_periodic = 0 1 0 +xlo.type = "Outflow" +xhi.type = "Outflow" +zlo.type = "SlipWall" +zhi.type = "Outflow" + +erf.sponge_strength = 2.0 +#erf.use_zhi_sponge_damping = true +erf.zhi_sponge_start = 12000.0 + +erf.sponge_density = 1.2 +erf.sponge_x_velocity = 0.0 +erf.sponge_y_velocity = 0.0 +erf.sponge_z_velocity = 0.0 + + + +# TIME STEP CONTROL +erf.use_native_mri = 1 +erf.fixed_dt = 1.0 # fixed time step [s] -- Straka et al 1993 +erf.fixed_fast_dt = 0.5 # fixed time step [s] -- Straka et al 1993 +#erf.no_substepping = 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 +amr.check_file = chk # root name of checkpoint file +amr.check_int = 1000 # number of timesteps between checkpoints +#amr.restart = chk09000 + +# PLOTFILES +erf.plot_file_1 = plt # root name of plotfile +erf.plot_int_1 = 100 # number of timesteps between plotfiles +erf.plot_vars_1 = density rhotheta rhoQt rhoQp x_velocity y_velocity z_velocity pressure theta temp qt qp qv qc qi scalar pert_dens + +# SOLVER CHOICE +erf.use_gravity = true +erf.buoyancy_type = 4 +erf.use_coriolis = false +erf.use_rayleigh_damping = false + +#erf.les_type = "Smagorinsky" +erf.Cs = 0.25 +erf.les_type = "None" + +# +# diffusion coefficient from Straka, K = 75 m^2/s +# +erf.molec_diff_type = "ConstantAlpha" +#erf.molec_diff_type = "Constant" +erf.rho0_trans = 1.0 # [kg/m^3], used to convert input diffusivities +erf.dynamicViscosity = 200.0 # [kg/(m-s)] ==> nu = 75.0 m^2/s +erf.alpha_T = 00.0 # [m^2/s] +erf.alpha_C = 0.0 + +erf.moisture_model = "Kessler" +erf.use_moist_background = true + +erf.moistscal_horiz_adv_string = "Centered_2nd" +erf.moistscal_vert_adv_string = "Centered_2nd" + +# PROBLEM PARAMETERS (optional) +prob.T_0 = 300.0 +prob.U_0 = 0 +prob.T_pert = 3 diff --git a/Exec/SquallLine_2D/prob.H b/Exec/SquallLine_2D/prob.H new file mode 100644 index 000000000..1d5891061 --- /dev/null +++ b/Exec/SquallLine_2D/prob.H @@ -0,0 +1,89 @@ +#ifndef _PROB_H_ +#define _PROB_H_ + +#include + +#include "AMReX_REAL.H" +#include "AMReX_ParmParse.H" +#include "AMReX_MultiFab.H" + +#include "prob_common.H" +#include "EOS.H" +#include "IndexDefines.H" +#include "TileNoZ.H" + +struct ProbParm : ProbParmDefaults { + amrex::Real T_0 = 300.0; // surface temperature == mean potential temperature + amrex::Real U_0 = 10.0; + amrex::Real V_0 = 0.0; + amrex::Real x_c = 0.0; // center of thermal perturbation + amrex::Real z_c = 3200.0; + amrex::Real x_r = 1000.0; + amrex::Real z_r = 1000.0; + amrex::Real T_pert = 5.0; // perturbation temperature + // overridden physical constants + amrex::Real C_p = 1004.0; +}; // namespace ProbParm + + +class Problem : public ProblemBase +{ +public: + Problem(); + +#include "Prob/init_density_hse_dry.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; + + + void erf_init_dens_hse_moist (amrex::MultiFab& rho_hse, + std::unique_ptr& z_phys_nd, + std::unique_ptr& z_phys_cc, + amrex::Geometry const& geom) override; + + void init_custom_terrain ( + const amrex::Geometry& geom, + amrex::MultiFab& z_phys_nd, + const amrex::Real& time) override; + + void erf_init_rayleigh ( + amrex::Vector& tau, + amrex::Vector& ubar, + amrex::Vector& vbar, + amrex::Vector& wbar, + amrex::Vector& thetabar, + amrex::Geometry const& geom) override; + +protected: + std::string name () override { return "Supercell"; } + +private: + ProbParm parms; +}; + +#endif diff --git a/Exec/SquallLine_2D/prob.cpp b/Exec/SquallLine_2D/prob.cpp new file mode 100644 index 000000000..8b7f2cc65 --- /dev/null +++ b/Exec/SquallLine_2D/prob.cpp @@ -0,0 +1,477 @@ +#include "prob.H" + +using namespace amrex; + +std::unique_ptr +amrex_probinit ( + const amrex_real* /*problo*/, + const amrex_real* /*probhi*/) +{ + return std::make_unique(); +} + +Problem::Problem () +{ + // Parse params + amrex::ParmParse pp("prob"); + pp.query("T_0", parms.T_0); + pp.query("U_0", parms.U_0); + pp.query("x_c", parms.x_c); + pp.query("z_c", parms.z_c); + pp.query("x_r", parms.x_r); + pp.query("z_r", parms.z_r); + pp.query("T_pert", parms.T_pert); +} + + +Real z_tr = 12000.0; +Real height = 1200.0; +Real R_v_by_R_d = R_v/R_d; + +Real compute_theta (Real z) +{ + Real theta_0=300.0, theta_tr=343.0, T_tr=213.0; + if(z <= z_tr) { + return theta_0 + (theta_tr - theta_0)*std::pow(z/z_tr,1.25); + } else { + return theta_tr*exp(CONST_GRAV/(Cp_d*T_tr)*(z - z_tr)); + } + +} + +AMREX_FORCE_INLINE +AMREX_GPU_HOST_DEVICE +Real compute_saturation_pressure (const Real T_b) +{ + + Real p_s = exp(34.494 - 4924.99/(T_b - 273.15 + 237.1))/std::pow(T_b - 273.15 + 105.0,1.57); + + Real T = T_b - 273.15; + + //Real p_s = 0.61121e3*exp((18.678 - T/234.5)*(T/(257.14 + T))); + + return p_s; +} + +AMREX_FORCE_INLINE +AMREX_GPU_HOST_DEVICE +Real compute_relative_humidity (Real z, Real height, Real z_tr, Real p_b, Real T_b) +{ + Real p_s = compute_saturation_pressure(T_b); + + Real q_s = 0.622*p_s/(p_b - p_s); + + if(z <= height){ + return 0.014/q_s; + }else if(z <= z_tr){ + return 1.0 - 0.75*std::pow(z/z_tr,1.25); + }else{ + return 0.25; + } +} + +AMREX_FORCE_INLINE +AMREX_GPU_HOST_DEVICE +Real compute_vapor_pressure (Real p_s, Real RH) +{ + return p_s*RH; +} + +AMREX_FORCE_INLINE +AMREX_GPU_HOST_DEVICE +Real vapor_mixing_ratio (const Real z, const Real height, 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 = 0.622*p_v/(p_b - p_v); + + if(z < height){ + return 0.014; + }else{ + return q_v; + } +} + +AMREX_FORCE_INLINE +AMREX_GPU_HOST_DEVICE +Real compute_temperature (const Real p_b, const Real theta_b) +{ + return theta_b*std::pow(p_b/p_0,R_d/Cp_d); +} + +Real compute_dewpoint_temperature (const Real z, const Real p_b, const Real T_b, const Real RH) +{ + + Real T_dp, gamma, T; + T = T_b - 273.15; + + Real a = 6.11e-3*p_0, 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; + +} + +void compute_rho (const Real z, const Real pressure, Real &theta, Real& rho, Real& q_v, Real& T_dp, Real& T_b) +{ + + theta = compute_theta(z); + T_b = compute_temperature(pressure, theta); + Real RH = compute_relative_humidity(z, height, z_tr, pressure, T_b); + q_v = vapor_mixing_ratio(z, height, pressure, T_b, RH); + rho = pressure/(R_d*T_b*(1.0 + (R_v/R_d)*q_v)); + rho = rho*(1.0 + q_v); + T_dp = compute_dewpoint_temperature(z, pressure, T_b, RH); +} + +double 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 z, 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(z, 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(z, 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; +} + +Real 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 z, const Real rho_k_minus_1) +{ + + 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, z, 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, z, rho_k_minus_1); + Real F_prime = (F_plus_dF - F)/1e-10; + Real delta_p_k = -F/F_prime; + p_k = p_k + delta_p_k; + } + + return p_k; +} + +void +init_isentropic_hse_no_terrain(Real *theta, Real* r, Real* p, Real *q_v, + const Real& dz, const Real& prob_lo_z, + const int& khi) +{ + Real z, p_b, theta_b, T_b, rho_b, RH, T_dp, rho0, theta0, T0, q_v0; + p_b = p_0; + //T_b = T0; + + // Compute the quantities at z = 0.5*dz (first cell center) + z = prob_lo_z + 0.5*dz; + p[0] = p_0; + compute_p_k(p[0], p_0, theta[0], r[0], q_v[0], T_dp, T_b, dz, z, 0.0); + + for (int k=1;k<=khi;k++){ + z = prob_lo_z + (k+0.5)*dz; + 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, z, r[k-1]); + } + + r[khi+1] = r[khi]; +} + +void +Problem::erf_init_dens_hse_moist (MultiFab& rho_hse, + std::unique_ptr& z_phys_nd, + std::unique_ptr& z_phys_cc, + Geometry const& geom) +{ + + const Real prob_lo_z = geom.ProbLo()[2]; + const Real dz = geom.CellSize()[2]; + const int khi = geom.Domain().bigEnd()[2]; + + + const Real T_sfc = 300.; + const Real rho_sfc = p_0 / (R_d*T_sfc); + const Real Thetabar = T_sfc; + + // use_terrain = 1 + if (z_phys_nd) { + + if (khi > 255) amrex::Abort("1D Arrays are hard-wired to only 256 high"); + + for ( MFIter mfi(rho_hse, TileNoZ()); mfi.isValid(); ++mfi ) + { + Array4 rho_arr = rho_hse.array(mfi); + Array4 z_cc_arr = z_phys_cc->const_array(mfi); + + // Create a flat box with same horizontal extent but only one cell in vertical + const Box& tbz = mfi.nodaltilebox(2); + Box b2d = tbz; // Copy constructor + b2d.grow(0,1); b2d.grow(1,1); // Grow by one in the lateral directions + b2d.setRange(2,0); + + ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int) { + Array1D r;; + Array1D p;; + + //init_isentropic_hse_terrain(i,j,rho_sfc,Thetabar,&(r(0)),&(p(0)),z_cc_arr,khi); + + for (int k = 0; k <= khi; k++) { + rho_arr(i,j,k) = r(k); + } + rho_arr(i,j, -1) = rho_arr(i,j,0); + rho_arr(i,j,khi+1) = rho_arr(i,j,khi); + }); + } // mfi + } 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,prob_lo_z,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(); + Real* q_v = d_q_v.data(); + Real* theta = d_t.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, + const Box& xbx, + const Box& ybx, + const Box& zbx, + Array4 const& state, + Array4 const& x_vel, + Array4 const& y_vel, + Array4 const& z_vel, + Array4 const& /*r_hse*/, + Array4 const& /*p_hse*/, + Array4 const& /*z_nd*/, + Array4 const& /*z_cc*/, +#if defined(ERF_USE_MOISTURE) + Array4 const& qv, + Array4 const& qc, + Array4 const& qi, +#elif defined(ERF_USE_WARM_NO_PRECIP) + Array4 const& , + Array4 const& , +#endif + GeometryData const& geomdata, + Array4 const& /*mf_m*/, + Array4 const& /*mf_u*/, + Array4 const& /*mf_v*/, + const SolverChoice& sc) +{ + const int khi = geomdata.Domain().bigEnd()[2]; + + AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1); + + // This is what we do at k = 0 -- note we assume p = p_0 and T = T_0 at z=0 + const amrex::Real& dz = geomdata.CellSize()[2]; + const amrex::Real& prob_lo_z = geomdata.ProbLo()[2]; + const amrex::Real& prob_hi_z = geomdata.ProbHi()[2]; + + const amrex::Real rdOcp = sc.rdOcp; + + // const amrex::Real thetatr = 343.0; + // const amrex::Real theta0 = 300.0; + const amrex::Real ztr = 12000.0; + const amrex::Real Ttr = 213.0; + const amrex::Real Ttop = 213.0; + const amrex::Real deltaz = 1000.*0.0; + const amrex::Real zs = 5000.; + const amrex::Real us = 30.; + const amrex::Real uc = 15.; + + // Call the routine to calculate the 1d background condition + + 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,prob_lo_z,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* t = d_t.data(); + Real* r = d_r.data(); + Real* q_v = d_q_v.data(); + Real* p = d_p.data(); + + Real d_z_tr = 12000.0; + Real d_height = 1200.0; + + const Real x_c = 0.0, z_c = 1.5e3, x_r = 4.0e3, z_r = 1.5e3, r_c = 1.0, theta_c = 3.0; + //const Real x_c = 0.0, z_c = 2.0e3, x_r = 10.0e3, z_r = 1.5e3, r_c = 1.0, theta_c = 3.0; + + amrex::ParallelForRNG(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) 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(); + const amrex::Real z = prob_lo[2] + (k + 0.5) * dx[2]; + const amrex::Real x = prob_lo[0] + (i + 0.5) * dx[0]; + Real rad, delta_theta, theta_total, temperature, rho, RH, scalar; + + // Introduce the warm bubble. Assume that the bubble is pressure matche with the background + + rad = std::pow(std::pow((x - x_c)/x_r,2) + std::pow((z - z_c)/z_r,2), 0.5); + + if(rad <= r_c){ + delta_theta = theta_c*std::pow(cos(PI*rad/2.0),2); + scalar = std::pow(cos(PI*rad/2.0),2); + } + else{ + delta_theta = 0.0; + scalar = 0.0; + } + + theta_total = t[k] + delta_theta; + temperature = compute_temperature(p[k], theta_total); + Real T_b = compute_temperature(p[k], t[k]); + RH = compute_relative_humidity(z, d_height, d_z_tr, p[k], T_b); + Real q_v_hot = vapor_mixing_ratio(z, d_height, p[k], T_b, RH); + rho = p[k]/(R_d*temperature*(1.0 + (R_v/R_d)*q_v_hot)); + + // Compute background quantities + Real temperature_back = compute_temperature(p[k], t[k]); + Real T_back = compute_temperature(p[k], t[k]); + Real RH_back = compute_relative_humidity(z, d_height, d_z_tr, p[k], T_back); + Real q_v_back = vapor_mixing_ratio(z, d_height, p[k], T_back, RH_back); + Real rho_back = p[k]/(R_d*temperature_back*(1.0 + (R_v/R_d)*q_v_back)); + + // This version perturbs rho but not p + state(i, j, k, RhoTheta_comp) = rho*theta_total - rho_back*t[k]*(1.0 + (R_v/R_d)*q_v_back);// rho*d_t[k]*(1.0 + R_v_by_R_d*q_v_hot); + state(i, j, k, Rho_comp) = rho - rho_back*(1.0 + q_v_back); + + // Set scalar = 0 everywhere + state(i, j, k, RhoScalar_comp) = rho*scalar; + + // mean states +#if defined(ERF_USE_MOISTURE) + state(i, j, k, RhoQt_comp) = rho*q_v_hot; + state(i, j, k, RhoQp_comp) = 0.0; + qv(i, j, k) = q_v_hot; + qc(i, j, k) = 0.0; + qi(i, j, k) = 0.0; +#elif defined(ERF_USE_WARM_NO_PRECIP) + state(i, j, k, RhoQv_comp) = 0.0;//rho*qvapor; + state(i, j, k, RhoQc_comp) = 0.0; +#endif + + }); + + // Set the x-velocity + amrex::ParallelFor(xbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + const amrex::Real z = prob_lo_z + (k+0.5) * dz; + x_vel(i,j,k) = -12.0*std::max(0.0, (2.5e3 - z)/2.5e3); + }); + + // Set the y-velocity + amrex::ParallelFor(ybx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + y_vel(i, j, k) = 0.0; + }); + + // Set the z-velocity + amrex::ParallelFor(zbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + z_vel(i, j, k) = 0.0; + }); + + amrex::Gpu::streamSynchronize(); +} + +void +Problem::init_custom_terrain (const Geometry& /*geom*/, + MultiFab& z_phys_nd, + const Real& /*time*/) +{ + // Number of ghost cells + int ngrow = z_phys_nd.nGrow(); + + // Bottom of domain + int k0 = 0; + + for ( MFIter mfi(z_phys_nd, TilingIfNotGPU()); mfi.isValid(); ++mfi ) + { + // Grown box with no z range + amrex::Box xybx = mfi.growntilebox(ngrow); + xybx.setRange(2,0); + + Array4 const& z_arr = z_phys_nd.array(mfi); + + ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) { + + // Flat terrain with z = 0 at k = 0 + z_arr(i,j,k0) = 0.0; + }); + } +} + +void +Problem::erf_init_rayleigh (amrex::Vector& tau, + amrex::Vector& ubar, + amrex::Vector& vbar, + amrex::Vector& wbar, + amrex::Vector& thetabar, + amrex::Geometry const& geom) +{ + const int khi = geom.Domain().bigEnd()[2]; + + // We just use these values to test the Rayleigh damping + for (int k = 0; k <= khi; k++) + { + tau[k] = 1.0; + ubar[k] = 2.0; + vbar[k] = 1.0; + wbar[k] = 0.0; + //thetabar[k] = parms.Theta_0; + } +} + + diff --git a/Source/DataStructs/DataStruct.H b/Source/DataStructs/DataStruct.H index ec6c285f9..970691474 100644 --- a/Source/DataStructs/DataStruct.H +++ b/Source/DataStructs/DataStruct.H @@ -64,8 +64,8 @@ struct SolverChoice { // Which expression (1,2 or 3) to use for buoyancy pp.query("buoyancy_type", buoyancy_type); - if (buoyancy_type != 1 && buoyancy_type != 2 && buoyancy_type != 3) { - amrex::Abort("buoyancy_type must be 1, 2 or 3"); + if (buoyancy_type != 1 && buoyancy_type != 2 && buoyancy_type != 3 && buoyancy_type != 4) { + amrex::Abort("buoyancy_type must be 1, 2, 3 or 4"); } // Is the terrain static or moving? @@ -162,6 +162,7 @@ struct SolverChoice { #ifdef ERF_USE_MOISTURE pp.query("mp_clouds", do_cloud); pp.query("mp_precip", do_precip); + pp.query("use_moist_background", use_moist_background); #endif // Use numerical diffusion? @@ -359,6 +360,7 @@ struct SolverChoice { // Microphysics params bool do_cloud {true}; bool do_precip {true}; + bool use_moist_background {false}; #endif }; #endif diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 34e68caf7..ab5040d55 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -965,6 +965,8 @@ ERF::ReadParameters () pp.query("moisture_model", moisture_model); if (moisture_model == "SAM") { micro.SetModel(); + } else if (moisture_model == "Kessler") { + micro.SetModel(); } else { micro.SetModel(); amrex::Print() << "WARNING: Compiled with moisture but using NullMoist model!\n"; diff --git a/Source/Initialization/ERF_init1d.cpp b/Source/Initialization/ERF_init1d.cpp index 3e3f82e6b..ce47b5f80 100644 --- a/Source/Initialization/ERF_init1d.cpp +++ b/Source/Initialization/ERF_init1d.cpp @@ -129,12 +129,21 @@ ERF::initHSE (int lev) MultiFab p_hse (base_state[lev], make_alias, 1, 1); // p_0 is second component MultiFab pi_hse(base_state[lev], make_alias, 2, 1); // pi_0 is third component +#ifdef ERF_USE_MOISTURE // Initial r_hse may or may not be in HSE -- defined in prob.cpp + if(solverChoice.use_moist_background){ + prob->erf_init_dens_hse_moist(r_hse, z_phys_nd[lev], z_phys_cc[lev], geom[lev]); + }else{ + prob->erf_init_dens_hse(r_hse, z_phys_nd[lev], z_phys_cc[lev], geom[lev]); + } +#else prob->erf_init_dens_hse(r_hse, z_phys_nd[lev], z_phys_cc[lev], geom[lev]); +#endif // This integrates up through column to update p_hse, pi_hse; // r_hse is not const b/c FillBoundary is called at the end for r_hse and p_hse erf_enforce_hse(lev, r_hse, p_hse, pi_hse, z_phys_cc[lev], z_phys_nd[lev]); + } void @@ -215,11 +224,11 @@ ERF::erf_enforce_hse (int lev, // Set value at surface from Newton iteration for rho pres_arr(i,j,k0 ) = p_0 - hz * rho_arr(i,j,k0) * l_gravity; - pi_arr(i,j,k0 ) = getExnergivenP(pres_arr(i,j,k0 ), rdOcp); + pi_arr(i,j,k0 ) = getExnergivenP(pres_arr(i,j,k0 ), rdOcp); // Set ghost cell with dz and rho at boundary pres_arr(i,j,k0-1) = p_0 + hz * rho_arr(i,j,k0) * l_gravity; - pi_arr(i,j,k0-1) = getExnergivenP(pres_arr(i,j,k0-1), rdOcp); + pi_arr(i,j,k0-1) = getExnergivenP(pres_arr(i,j,k0-1), rdOcp); Real dens_interp; if (l_use_terrain) { @@ -248,7 +257,7 @@ ERF::erf_enforce_hse (int lev, bx.setBig(0,domlo_x-1); ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { pres_arr(i,j,k) = pres_arr(domlo_x,j,k); - pi_arr(i,j,k) = getExnergivenP(pres_arr(i,j,k), rdOcp); + pi_arr(i,j,k) = getExnergivenP(pres_arr(i,j,k), rdOcp); }); } diff --git a/Source/Microphysics/Kessler/Diagnose_Kessler.cpp b/Source/Microphysics/Kessler/Diagnose_Kessler.cpp new file mode 100644 index 000000000..4b74bcf28 --- /dev/null +++ b/Source/Microphysics/Kessler/Diagnose_Kessler.cpp @@ -0,0 +1,42 @@ +#include "Microphysics.H" + +/** + * Computes diagnostic quantities like cloud ice/liquid and precipitation ice/liquid + * from the existing Microphysics variables. + */ +void Kessler::Diagnose () { + + auto qt = mic_fab_vars[MicVar_Kess::qt]; + auto qp = mic_fab_vars[MicVar_Kess::qp]; + auto qv = mic_fab_vars[MicVar_Kess::qv]; + auto qn = mic_fab_vars[MicVar_Kess::qn]; + auto qcl = mic_fab_vars[MicVar_Kess::qcl]; + auto qci = mic_fab_vars[MicVar_Kess::qci]; + auto qpl = mic_fab_vars[MicVar_Kess::qpl]; + auto qpi = mic_fab_vars[MicVar_Kess::qpi]; + auto tabs = mic_fab_vars[MicVar_Kess::tabs]; + + for ( amrex::MFIter mfi(*tabs, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { + auto qt_array = qt->array(mfi); + auto qp_array = qp->array(mfi); + auto qv_array = qv->array(mfi); + auto qn_array = qn->array(mfi); + auto qcl_array = qcl->array(mfi); + auto qci_array = qci->array(mfi); + auto qpl_array = qpl->array(mfi); + auto qpi_array = qpi->array(mfi); + auto tabs_array = tabs->array(mfi); + + const auto& box3d = mfi.tilebox(); + + amrex::ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + qv_array(i,j,k) = std::max(0.0,qt_array(i,j,k) - qn_array(i,j,k)); + amrex::Real omn = 1.0; + qcl_array(i,j,k) = qn_array(i,j,k)*omn; + qci_array(i,j,k) = qn_array(i,j,k)*(1.0-omn); + amrex::Real omp = 1.0;; + qpl_array(i,j,k) = qp_array(i,j,k)*omp; + qpi_array(i,j,k) = qp_array(i,j,k)*(1.0-omp); + }); + } +} diff --git a/Source/Microphysics/Kessler/Init_Kessler.cpp b/Source/Microphysics/Kessler/Init_Kessler.cpp new file mode 100644 index 000000000..0b97652cf --- /dev/null +++ b/Source/Microphysics/Kessler/Init_Kessler.cpp @@ -0,0 +1,249 @@ + +#include +#include "Microphysics.H" +#include "IndexDefines.H" +#include "PlaneAverage.H" +#include "EOS.H" +#include "TileNoZ.H" + +using namespace amrex; + +/** + * Initializes the Microphysics module. + * + * @param[in] cons_in Conserved variables input + * @param[in] qc_in Cloud variables input + * @param[in,out] qv_in Vapor variables input + * @param[in] qi_in Ice variables input + * @param[in] grids The boxes on which we will evolve the solution + * @param[in] geom Geometry associated with these MultiFabs and grids + * @param[in] dt_advance Timestep for the advance + */ +void Kessler::Init (const MultiFab& cons_in, MultiFab& qmoist, + const BoxArray& grids, + const Geometry& geom, + const Real& dt_advance) + { + m_geom = geom; + m_gtoe = grids; + + auto dz = m_geom.CellSize(2); + auto lowz = m_geom.ProbLo(2); + + dt = dt_advance; + + // initialize microphysics variables + for (auto ivar = 0; ivar < MicVar_Kess::NumVars; ++ivar) { + mic_fab_vars[ivar] = std::make_shared(cons_in.boxArray(), cons_in.DistributionMap(), 1, cons_in.nGrowVect()); + mic_fab_vars[ivar]->setVal(0.); + } + + // We must initialize to zero since we now need boundary values for the call to getP and we need all values filled + // The ghost cells of these arrays aren't filled in the boundary condition calls for the state + + //qmoist.setVal(0.); + + for ( MFIter mfi(cons_in, TileNoZ()); mfi.isValid(); ++mfi) { + + const auto& box3d = mfi.tilebox(); + + const auto& lo = amrex::lbound(box3d); + const auto& hi = amrex::ubound(box3d); + + nlev = box3d.length(2); + zlo = lo.z; + zhi = hi.z; + + // parameters + accrrc.resize({zlo}, {zhi}); + accrsi.resize({zlo}, {zhi}); + accrsc.resize({zlo}, {zhi}); + coefice.resize({zlo}, {zhi}); + evaps1.resize({zlo}, {zhi}); + evaps2.resize({zlo}, {zhi}); + accrgi.resize({zlo}, {zhi}); + accrgc.resize({zlo}, {zhi}); + evapg1.resize({zlo}, {zhi}); + evapg2.resize({zlo}, {zhi}); + evapr1.resize({zlo}, {zhi}); + evapr2.resize({zlo}, {zhi}); + + // data (input) + rho1d.resize({zlo}, {zhi}); + pres1d.resize({zlo}, {zhi}); + tabs1d.resize({zlo}, {zhi}); + gamaz.resize({zlo}, {zhi}); + zmid.resize({zlo}, {zhi}); + + // output + qifall.resize({zlo}, {zhi}); + tlatqi.resize({zlo}, {zhi}); + } + + auto accrrc_t = accrrc.table(); + auto accrsi_t = accrsi.table(); + auto accrsc_t = accrsc.table(); + auto coefice_t = coefice.table(); + auto evaps1_t = evaps1.table(); + auto evaps2_t = evaps2.table(); + auto accrgi_t = accrgi.table(); + auto accrgc_t = accrgc.table(); + auto evapg1_t = evapg1.table(); + auto evapg2_t = evapg2.table(); + auto evapr1_t = evapr1.table(); + auto evapr2_t = evapr2.table(); + + auto rho1d_t = rho1d.table(); + auto pres1d_t = pres1d.table(); + auto tabs1d_t = tabs1d.table(); + + auto gamaz_t = gamaz.table(); + auto zmid_t = zmid.table(); + + Real gam3 = erf_gammafff(3.0 ); + Real gamr1 = erf_gammafff(3.0+b_rain ); + Real gamr2 = erf_gammafff((5.0+b_rain)/2.0); + // Real gamr3 = erf_gammafff(4.0+b_rain ); + Real gams1 = erf_gammafff(3.0+b_snow ); + Real gams2 = erf_gammafff((5.0+b_snow)/2.0); + // Real gams3 = erf_gammafff(4.0+b_snow ); + Real gamg1 = erf_gammafff(3.0+b_grau ); + Real gamg2 = erf_gammafff((5.0+b_grau)/2.0); + // Real gamg3 = erf_gammafff(4.0+b_grau ); + + amrex::MultiFab qv(qmoist, amrex::make_alias, 0, 1); + amrex::MultiFab qc(qmoist, amrex::make_alias, 1, 1); + amrex::MultiFab qi(qmoist, amrex::make_alias, 2, 1); + + // Get the temperature, density, theta, qt and qp from input + for ( MFIter mfi(cons_in, false); mfi.isValid(); ++mfi) { + auto states_array = cons_in.array(mfi); + auto qv_array_from_moist = qv.array(mfi); + auto qc_array = qc.array(mfi); + auto qi_array = qi.array(mfi); + + auto qt_array = mic_fab_vars[MicVar_Kess::qt]->array(mfi); + auto qp_array = mic_fab_vars[MicVar_Kess::qp]->array(mfi); + auto qn_array = mic_fab_vars[MicVar_Kess::qn]->array(mfi); + auto qv_array = mic_fab_vars[MicVar_Kess::qv]->array(mfi); + auto rho_array = mic_fab_vars[MicVar_Kess::rho]->array(mfi); + auto theta_array = mic_fab_vars[MicVar_Kess::theta]->array(mfi); + auto temp_array = mic_fab_vars[MicVar_Kess::tabs]->array(mfi); + auto pres_array = mic_fab_vars[MicVar_Kess::pres]->array(mfi); + + const auto& box3d = mfi.tilebox(); + + // Get pressure, theta, temperature, density, and qt, qp + amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { + rho_array(i,j,k) = states_array(i,j,k,Rho_comp); + theta_array(i,j,k) = states_array(i,j,k,RhoTheta_comp)/states_array(i,j,k,Rho_comp); + qt_array(i,j,k) = states_array(i,j,k,RhoQt_comp)/states_array(i,j,k,Rho_comp); + qp_array(i,j,k) = states_array(i,j,k,RhoQp_comp)/states_array(i,j,k,Rho_comp); + qn_array(i,j,k) = qc_array(i,j,k) + qi_array(i,j,k); + qv_array(i,j,k) = qv_array_from_moist(i,j,k); + temp_array(i,j,k) = getTgivenRandRTh(states_array(i,j,k,Rho_comp),states_array(i,j,k,RhoTheta_comp), qv_array(i,j,k)); + pres_array(i,j,k) = getPgivenRTh(states_array(i,j,k,RhoTheta_comp))/100.; + }); + } + + // calculate the plane average variables + PlaneAverage cons_ave(&cons_in, m_geom, m_axis); + cons_ave.compute_averages(ZDir(), cons_ave.field()); + + // get host variable rho, and rhotheta + int ncell = cons_ave.ncell_line(); + + Gpu::HostVector rho_h(ncell), rhotheta_h(ncell); + cons_ave.line_average(Rho_comp, rho_h); + cons_ave.line_average(RhoTheta_comp, rhotheta_h); + + // copy data to device + Gpu::DeviceVector rho_d(ncell), rhotheta_d(ncell); + Gpu::copyAsync(Gpu::hostToDevice, rho_h.begin(), rho_h.end(), rho_d.begin()); + Gpu::copyAsync(Gpu::hostToDevice, rhotheta_h.begin(), rhotheta_h.end(), rhotheta_d.begin()); + Gpu::streamSynchronize(); + + Real* rho_dptr = rho_d.data(); + Real* rhotheta_dptr = rhotheta_d.data(); + + Real gOcp = m_gOcp; + + amrex::ParallelFor(nlev, [=] AMREX_GPU_DEVICE (int k) noexcept { + Real pressure = getPgivenRTh(rhotheta_dptr[k]); + rho1d_t(k) = rho_dptr[k]; + pres1d_t(k) = pressure/100.; + tabs1d_t(k) = getTgivenRandRTh(rho_dptr[k],rhotheta_dptr[k]); + zmid_t(k) = lowz + (k+0.5)*dz; + gamaz_t(k) = gOcp*zmid_t(k); + }); + + // This fills qv + //Diagnose(); + +#if 0 + amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int k, int j, int i) { + fluxbmk(l,j,i) = 0.0; + fluxtmk(l,j,i) = 0.0; + }); + + amrex::ParallelFor( box2d, [=] AMREX_GPU_DEVICE (int k, int l, int i0) { + mkwle (l,k) = 0.0; + mkwsb (l,k) = 0.0; + mkadv (l,k) = 0.0; + mkdiff(l,k) = 0.0; + }); +#endif + + if(round(gam3) != 2) { + std::cout << "cannot compute gamma-function in Microphysics::Init" << std::endl; + std::exit(-1); + } + + amrex::ParallelFor(nlev, [=] AMREX_GPU_DEVICE (int k) noexcept { + + Real pratio = sqrt(1.29 / rho1d_t(k)); + Real rrr1=393.0/(tabs1d_t(k)+120.0)*std::pow((tabs1d_t(k)/273.0),1.5); + Real rrr2=std::pow((tabs1d_t(k)/273.0),1.94)*(1000.0/pres1d_t(k)); + Real estw = 100.0*erf_esatw(tabs1d_t(k)); + Real esti = 100.0*erf_esati(tabs1d_t(k)); + + // accretion by snow: + Real coef1 = 0.25 * PI * nzeros * a_snow * gams1 * pratio/pow((PI * rhos * nzeros/rho1d_t(k) ) , ((3.0+b_snow)/4.0)); + Real coef2 = exp(0.025*(tabs1d_t(k) - 273.15)); + accrsi_t(k) = coef1 * coef2 * esicoef; + accrsc_t(k) = coef1 * esccoef; + coefice_t(k) = coef2; + + // evaporation of snow: + coef1 =(lsub/(tabs1d_t(k)*rv)-1.0)*lsub/(therco*rrr1*tabs1d_t(k)); + coef2 = rv*tabs1d_t(k)/(diffelq*rrr2*esti); + evaps1_t(k) = 0.65*4.0*nzeros/sqrt(PI*rhos*nzeros)/(coef1+coef2)/sqrt(rho1d_t(k)); + evaps2_t(k) = 0.49*4.0*nzeros*gams2*sqrt(a_snow/(muelq*rrr1))/pow((PI*rhos*nzeros) , ((5.0+b_snow)/8.0)) / + (coef1+coef2) * pow(rho1d_t(k) , ((1.0+b_snow)/8.0))*sqrt(pratio); + + // accretion by graupel: + coef1 = 0.25*PI*nzerog*a_grau*gamg1*pratio/pow((PI*rhog*nzerog/rho1d_t(k)) , ((3.0+b_grau)/4.0)); + coef2 = exp(0.025*(tabs1d_t(k) - 273.15)); + accrgi_t(k) = coef1 * coef2 * egicoef; + accrgc_t(k) = coef1 * egccoef; + + // evaporation of graupel: + coef1 =(lsub/(tabs1d_t(k)*rv)-1.0)*lsub/(therco*rrr1*tabs1d_t(k)); + coef2 = rv*tabs1d_t(k)/(diffelq*rrr2*esti); + evapg1_t(k) = 0.65*4.0*nzerog/sqrt(PI*rhog*nzerog)/(coef1+coef2)/sqrt(rho1d_t(k)); + evapg2_t(k) = 0.49*4.0*nzerog*gamg2*sqrt(a_grau/(muelq*rrr1))/pow((PI * rhog * nzerog) , ((5.0+b_grau)/8.0)) / + (coef1+coef2) * pow(rho1d_t(k) , ((1.0+b_grau)/8.0))*sqrt(pratio); + + // accretion by rain: + accrrc_t(k)= 0.25 * PI * nzeror * a_rain * gamr1 * pratio/pow((PI * rhor * nzeror / rho1d_t(k)) , ((3+b_rain)/4.))* erccoef; + + // evaporation of rain: + coef1 =(lcond/(tabs1d_t(k)*rv)-1.0)*lcond/(therco*rrr1*tabs1d_t(k)); + coef2 = rv*tabs1d_t(k)/(diffelq * rrr2 * estw); + evapr1_t(k) = 0.78 * 2.0 * PI * nzeror / + sqrt(PI * rhor * nzeror) / (coef1+coef2) / sqrt(rho1d_t(k)); + evapr2_t(k) = 0.31 * 2.0 * PI * nzeror * gamr2 * 0.89 * sqrt(a_rain/(muelq*rrr1))/ + pow((PI * rhor * nzeror) , ((5.0+b_rain)/8.0)) / + (coef1+coef2) * pow(rho1d_t(k) , ((1.0+b_rain)/8.0))*sqrt(pratio); + }); +} diff --git a/Source/Microphysics/Kessler/Kessler.H b/Source/Microphysics/Kessler/Kessler.H new file mode 100644 index 000000000..7ef58dd34 --- /dev/null +++ b/Source/Microphysics/Kessler/Kessler.H @@ -0,0 +1,159 @@ +/* + * Implementation 1-moment microphysics model + * NOTE: this model is based on the Kessler code, and the Klemp's paper + * 1): Joseph, Klemp, the simulation of three-dimensional convective storm dynamics, + * Journal of the atmospheric sciences, vol35, p1070 + * 2): Marat Khairoutdinov and David Randall, cloud resolving modeling of the ARM summer 1997 IOP: + * model formulation, results, unvertainties, and sensitivities, Journal of the atmospheric sciences, vol60, p607 + */ +#ifndef Kessler_H +#define Kessler_H + +#include +#include +#include + +#include +#include +#include +#include + +#include "ERF_Constants.H" +#include "Microphysics_Utils.H" +#include "IndexDefines.H" +#include "DataStruct.H" + +namespace MicVar_Kess { + enum { + // independent variables + qt = 0, + qp, + theta, // liquid/ice water potential temperature + tabs, // temperature + rho, // density + pres, // pressure + // derived variables + qr, // rain water + qv, // water vapor + qn, // cloud condensate (liquid+ice), initial to zero + qci, // cloud ice + qcl, // cloud water + qpl, // precip rain + qpi, // precip ice + // temporary variable + omega, + NumVars + }; +} + +// +// use MultiFab for 3D data, but table for 1D data +// +class Kessler : public NullMoist { + + using FabPtr = std::shared_ptr; + + public: + // constructor + Kessler () {} + + // destructor + virtual ~Kessler () = default; + + // cloud physics + void AdvanceKessler (); + + // diagnose + void Diagnose () override; + + // Set up for first time + void + define (SolverChoice& sc) override + { + docloud = sc.do_cloud; + doprecip = sc.do_precip; + m_fac_cond = lcond / sc.c_p; + m_fac_fus = lfus / sc.c_p; + m_fac_sub = lsub / sc.c_p; + m_gOcp = CONST_GRAV / sc.c_p; + m_axis = sc.ave_plane; + } + + // init + void + Init (const amrex::MultiFab& cons_in, + amrex::MultiFab& qmoist, + const amrex::BoxArray& grids, + const amrex::Geometry& geom, + const amrex::Real& dt_advance) override; + + // update ERF variables + void + Update (amrex::MultiFab& cons_in, + amrex::MultiFab& qmoist) override; + + // wrapper to do all the updating + void + Advance ( ) override + { + this->AdvanceKessler(); + this->Diagnose(); + } + + private: + // geometry + amrex::Geometry m_geom; + // valid boxes on which to evolve the solution + amrex::BoxArray m_gtoe; + + // timestep + amrex::Real dt; + + // number of vertical levels + int nlev, zlo, zhi; + + // plane average axis + int m_axis; + + // model options + bool docloud, doprecip; + + // constants + amrex::Real m_fac_cond; + amrex::Real m_fac_fus; + amrex::Real m_fac_sub; + amrex::Real m_gOcp; + + // microphysics parameters/coefficients + amrex::TableData accrrc; + amrex::TableData accrsi; + amrex::TableData accrsc; + amrex::TableData coefice; + amrex::TableData evaps1; + amrex::TableData evaps2; + amrex::TableData accrgi; + amrex::TableData accrgc; + amrex::TableData evapg1; + amrex::TableData evapg2; + amrex::TableData evapr1; + amrex::TableData evapr2; + + // vertical plane average data + amrex::TableData rho1d; + amrex::TableData pres1d; + amrex::TableData tabs1d; + amrex::TableData qt1d; + amrex::TableData qv1d; + amrex::TableData qn1d; + + // independent variables + amrex::Array mic_fab_vars; + + amrex::TableData gamaz; + amrex::TableData zmid; // mid value of vertical coordinate in physical domain + + // data (output) + amrex::TableData qifall; + amrex::TableData tlatqi; +}; +#endif diff --git a/Source/Microphysics/Kessler/Kessler.cpp b/Source/Microphysics/Kessler/Kessler.cpp new file mode 100644 index 000000000..d12cd566f --- /dev/null +++ b/Source/Microphysics/Kessler/Kessler.cpp @@ -0,0 +1,238 @@ +/* + * this file is modified from precip_proc from samxx + */ +#include "Microphysics.H" + +#include +#include + +using namespace amrex; + +/** + * Compute Precipitation-related Microphysics quantities. + */ +void Kessler::AdvanceKessler() { + + Real powr1 = (3.0 + b_rain) / 4.0; + Real powr2 = (5.0 + b_rain) / 8.0; + Real pows1 = (3.0 + b_snow) / 4.0; + Real pows2 = (5.0 + b_snow) / 8.0; + Real powg1 = (3.0 + b_grau) / 4.0; + Real powg2 = (5.0 + b_grau) / 8.0; + + auto accrrc_t = accrrc.table(); + auto accrsc_t = accrsc.table(); + auto accrsi_t = accrsi.table(); + auto accrgc_t = accrgc.table(); + auto accrgi_t = accrgi.table(); + auto coefice_t = coefice.table(); + auto evapr1_t = evapr1.table(); + auto evapr2_t = evapr2.table(); + auto evaps1_t = evaps1.table(); + auto evaps2_t = evaps2.table(); + auto evapg1_t = evapg1.table(); + auto evapg2_t = evapg2.table(); + auto pres1d_t = pres1d.table(); + + auto qt = mic_fab_vars[MicVar_Kess::qt]; + auto qp = mic_fab_vars[MicVar_Kess::qp]; + auto qn = mic_fab_vars[MicVar_Kess::qn]; + auto tabs = mic_fab_vars[MicVar_Kess::tabs]; + + auto qcl = mic_fab_vars[MicVar_Kess::qcl]; + auto theta = mic_fab_vars[MicVar_Kess::theta]; + auto qv = mic_fab_vars[MicVar_Kess::qv]; + auto rho = mic_fab_vars[MicVar_Kess::rho]; + + auto dz = m_geom.CellSize(2); + auto domain = m_geom.Domain(); + int k_lo = domain.smallEnd(2); + int k_hi = domain.bigEnd(2); + + MultiFab fz; + auto ba = tabs->boxArray(); + auto dm = tabs->DistributionMap(); + auto ngrow = tabs->nGrowVect(); + fz.define(convert(ba, IntVect(0,0,1)), dm, 1, 0); // No ghost cells + + for ( MFIter mfi(fz, TilingIfNotGPU()); mfi.isValid(); ++mfi ){ + auto rho_array = mic_fab_vars[MicVar_Kess::rho]->array(mfi); + auto qp_array = mic_fab_vars[MicVar_Kess::qp]->array(mfi); + auto fz_array = fz.array(mfi); + const Box& tbz = mfi.tilebox(); + + ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + + Real rho_avg, qp_avg; + + if (k==k_lo) { + rho_avg = rho_array(i,j,k); + qp_avg = qp_array(i,j,k); + } else if (k==k_hi+1) { + rho_avg = rho_array(i,j,k-1); + qp_avg = qp_array(i,j,k-1); + } else { + rho_avg = 0.5*(rho_array(i,j,k-1) + rho_array(i,j,k)); // Convert to g/cm^3 + qp_avg = 0.5*(qp_array(i,j,k-1) + qp_array(i,j,k)); + } + + qp_avg = std::max(0.0, qp_avg); + + Real V_terminal = 36.34*std::pow(rho_avg*0.001*qp_avg, 0.1346)*std::pow(rho_avg/1.16, -0.5); // in m/s + + fz_array(i,j,k) = rho_avg*V_terminal*qp_avg; + + /*if(k==0){ + fz_array(i,j,k) = 0; + }*/ + }); + } + + Real dtn = dt; + + /*for ( MFIter mfi(*tabs,TilingIfNotGPU()); mfi.isValid(); ++mfi) { + + auto qn_array = mic_fab_vars[MicVar_Kess::qn]->array(mfi); + auto qt_array = mic_fab_vars[MicVar_Kess::qt]->array(mfi); + auto qp_array = mic_fab_vars[MicVar_Kess::qp]->array(mfi); + auto qv_array = mic_fab_vars[MicVar_Kess::qv]->array(mfi); + + const auto& box3d = mfi.tilebox(); + + ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + + qt_array(i,j,k) = std::max(0.0, qt_array(i,j,k)); + qp_array(i,j,k) = std::max(0.0, qp_array(i,j,k)); + qn_array(i,j,k) = std::max(0.0, qn_array(i,j,k)); + + if(qt_array(i,j,k) == 0.0){ + qv_array(i,j,k) = 0.0; + qn_array(i,j,k) = 0.0; + } + }); + } + + + for ( MFIter mfi(*tabs,TilingIfNotGPU()); mfi.isValid(); ++mfi) { + auto qp_array = mic_fab_vars[MicVar_Kess::qp]->array(mfi); + auto rho_array = mic_fab_vars[MicVar_Kess::rho]->array(mfi); + + const auto& box3d = mfi.tilebox(); + auto fz_array = fz.array(mfi); + + ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + Real dq_sed = 1.0/rho_array(i,j,k)*(fz_array(i,j,k+1) - fz_array(i,j,k))/dz*dtn; + qp_array(i,j,k) = qp_array(i,j,k) + dq_sed; + }); + }*/ + + + + + // get the temperature, dentisy, theta, qt and qp from input + for ( MFIter mfi(*tabs,TilingIfNotGPU()); mfi.isValid(); ++mfi) { + auto tabs_array = mic_fab_vars[MicVar_Kess::tabs]->array(mfi); + auto qn_array = mic_fab_vars[MicVar_Kess::qn]->array(mfi); + auto qt_array = mic_fab_vars[MicVar_Kess::qt]->array(mfi); + auto qp_array = mic_fab_vars[MicVar_Kess::qp]->array(mfi); + + auto qcl_array = mic_fab_vars[MicVar_Kess::qcl]->array(mfi); + auto theta_array = theta->array(mfi); + auto qv_array = qv->array(mfi); + auto rho_array = mic_fab_vars[MicVar_Kess::rho]->array(mfi); + + const auto& box3d = mfi.tilebox(); + + auto fz_array = fz.array(mfi); + + // Expose for GPU + Real d_fac_cond = m_fac_cond; + + ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + + qt_array(i,j,k) = std::max(0.0, qt_array(i,j,k)); + qp_array(i,j,k) = std::max(0.0, qp_array(i,j,k)); + qn_array(i,j,k) = std::max(0.0, qn_array(i,j,k)); + + if(qt_array(i,j,k) == 0.0){ + qv_array(i,j,k) = 0.0; + qn_array(i,j,k) = 0.0; + } + + //------- Autoconversion/accretion + Real omn, omp, omg, qcc, qii, autor, autos, accrr, qrr, accrcs, accris, + qss, accrcg, accrig, tmp, qgg, dq_clwater_to_rain, dq_rain_to_vapor, dq_clwater_to_vapor, dq_vapor_to_clwater, qsat; + + Real pressure = getPgivenRTh(rho_array(i,j,k)*theta_array(i,j,k),qv_array(i,j,k))/100.0; + erf_qsatw(tabs_array(i,j,k), pressure, qsat); + + // If there is precipitating water (i.e. rain), and the cell is not saturated + // then the rain water can evaporate leading to extraction of latent heat, hence + // reducing temperature and creating negative buoyancy + + dq_clwater_to_rain = 0.0; + dq_rain_to_vapor = 0.0; + dq_vapor_to_clwater = 0.0; + dq_clwater_to_vapor = 0.0; + + + //Real fac = qsat*4093.0*L_v/(Cp_d*std::pow(tabs_array(i,j,k)-36.0,2)); + Real fac = qsat*L_v*L_v/(Cp_d*R_v*tabs_array(i,j,k)*tabs_array(i,j,k)); + + // If water vapor content exceeds saturation value, then vapor condenses to waterm and latent heat is released, increasing temperature + if(qv_array(i,j,k) > qsat){ + dq_vapor_to_clwater = std::min(qv_array(i,j,k), (qv_array(i,j,k)-qsat)/(1.0 + fac)); + } + + + // If water vapor is less than the satruated value, then the cloud water can evaporate, leading to evaporative cooling and + // reducing temperature + if(qv_array(i,j,k) < qsat and qn_array(i,j,k) > 0.0){ + dq_clwater_to_vapor = std::min(qn_array(i,j,k), (qsat - qv_array(i,j,k))/(1.0 + fac)); + } + + if(qp_array(i,j,k) > 0.0 && qv_array(i,j,k) < qsat) { + Real C = 1.6 + 124.9*std::pow(0.001*rho_array(i,j,k)*qp_array(i,j,k),0.2046); + dq_rain_to_vapor = 1.0/(0.001*rho_array(i,j,k))*(1.0 - qv_array(i,j,k)/qsat)*C*std::pow(0.001*rho_array(i,j,k)*qp_array(i,j,k),0.525)/ + (5.4e5 + 2.55e6/(pressure*qsat))*dtn; + // The negative sign is to make this variable (vapor formed from evaporation) + // a poistive quantity (as qv/qs < 1) + dq_rain_to_vapor = std::min({qp_array(i,j,k), dq_rain_to_vapor}); + + // Removing latent heat due to evaporation from rain water to water vapor, reduces the (potential) temperature + } + + // If there is cloud water present then do accretion and autoconversion to rain + + if (qn_array(i,j,k) > 0.0) { + qcc = qn_array(i,j,k); + + autor = 0.0; + if (qcc > qcw0) { + autor = 0.001; + } + + accrr = 0.0; + accrr = 2.2 * std::pow(qp_array(i,j,k) , 0.875); + dq_clwater_to_rain = dtn *(accrr*qcc + autor*(qcc - 0.001)); + + // If the amount of change is more than the amount of qc present, then dq = qc + dq_clwater_to_rain = std::min(dq_clwater_to_rain, qn_array(i,j,k)); + } + + + Real dq_sed = 1.0/rho_array(i,j,k)*(fz_array(i,j,k+1) - fz_array(i,j,k))/dz*dtn; + + qt_array(i,j,k) = qt_array(i,j,k) + dq_rain_to_vapor - dq_clwater_to_rain; + qp_array(i,j,k) = qp_array(i,j,k) + dq_sed + dq_clwater_to_rain - dq_rain_to_vapor; + qn_array(i,j,k) = qn_array(i,j,k) + dq_vapor_to_clwater - dq_clwater_to_vapor - dq_clwater_to_rain; + + theta_array(i,j,k) = theta_array(i,j,k) + theta_array(i,j,k)/tabs_array(i,j,k)*d_fac_cond*(dq_vapor_to_clwater - dq_clwater_to_vapor - dq_rain_to_vapor); + + qt_array(i,j,k) = std::max(0.0, qt_array(i,j,k)); + qp_array(i,j,k) = std::max(0.0, qp_array(i,j,k)); + qn_array(i,j,k) = std::max(0.0, qn_array(i,j,k)); + + }); + } +} diff --git a/Source/Microphysics/Kessler/Make.package b/Source/Microphysics/Kessler/Make.package new file mode 100644 index 000000000..6ea2aa118 --- /dev/null +++ b/Source/Microphysics/Kessler/Make.package @@ -0,0 +1,5 @@ +CEXE_sources += Init_Kessler.cpp +CEXE_sources += Diagnose_Kessler.cpp +CEXE_sources += Update_Kessler.cpp +CEXE_sources += Kessler.cpp +CEXE_headers += Kessler.H diff --git a/Source/Microphysics/Kessler/Update_Kessler.cpp b/Source/Microphysics/Kessler/Update_Kessler.cpp new file mode 100644 index 000000000..fcef4df9b --- /dev/null +++ b/Source/Microphysics/Kessler/Update_Kessler.cpp @@ -0,0 +1,60 @@ + +#include "Microphysics.H" +#include "IndexDefines.H" +#include "TileNoZ.H" + +/** + * Updates conserved and microphysics variables in the provided MultiFabs from + * the internal MultiFabs that store Microphysics module data. + * + * @param[out] cons Conserved variables + * @param[out] qmoist: qv, qc, qi, qr, qs, qg + */ +void Kessler::Update (amrex::MultiFab& cons, + amrex::MultiFab& qmoist) +{ + // copy multifab data to qc, qv, and qi + amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar_Kess::qv], 0, 0, 1, mic_fab_vars[MicVar_Kess::qv]->nGrowVect()); // vapor + amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar_Kess::qcl], 0, 1, 1, mic_fab_vars[MicVar_Kess::qcl]->nGrowVect()); // cloud water + amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar_Kess::qci], 0, 2, 1, mic_fab_vars[MicVar_Kess::qci]->nGrowVect()); // cloud ice + amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar_Kess::qpl], 0, 3, 1, mic_fab_vars[MicVar_Kess::qpl]->nGrowVect()); // rain + amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar_Kess::qpi], 0, 4, 1, mic_fab_vars[MicVar_Kess::qpi]->nGrowVect()); // snow + + // Don't need to copy this since it is filled below + // amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar_Kess::qpi], 0, 5, 1, mic_fab_vars[MicVar_Kess::qci]->nGrowVect()); // graupel + + amrex::MultiFab qgraup_mf(qmoist, amrex::make_alias, 5, 1); + + // Get the temperature, density, theta, qt and qp from input + for ( amrex::MFIter mfi(cons,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { + auto states_arr = cons.array(mfi); + + auto rho_arr = mic_fab_vars[MicVar_Kess::rho]->array(mfi); + auto theta_arr = mic_fab_vars[MicVar_Kess::theta]->array(mfi); + auto qt_arr = mic_fab_vars[MicVar_Kess::qt]->array(mfi); + auto qp_arr = mic_fab_vars[MicVar_Kess::qp]->array(mfi); + auto qpl_arr = mic_fab_vars[MicVar_Kess::qpl]->array(mfi); + auto qpi_arr = mic_fab_vars[MicVar_Kess::qpi]->array(mfi); + + auto qgraup_arr= qgraup_mf.array(mfi); + + const auto& box3d = mfi.tilebox(); + + // get potential total density, temperature, qt, qp + amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { + states_arr(i,j,k,Rho_comp) = rho_arr(i,j,k); + states_arr(i,j,k,RhoTheta_comp) = rho_arr(i,j,k)*theta_arr(i,j,k); + states_arr(i,j,k,RhoQt_comp) = rho_arr(i,j,k)*qt_arr(i,j,k); + states_arr(i,j,k,RhoQp_comp) = rho_arr(i,j,k)*qp_arr(i,j,k); + + // Graupel == precip total - rain - snow (but must be >= 0) + qgraup_arr(i,j,k) = 0.0;// + }); + } + + // Fill interior ghost cells and periodic boundaries + cons.FillBoundary(m_geom.periodicity()); + qmoist.FillBoundary(m_geom.periodicity()); +} + + diff --git a/Source/Microphysics/Microphysics.H b/Source/Microphysics/Microphysics.H index 479ce7575..1338f7674 100644 --- a/Source/Microphysics/Microphysics.H +++ b/Source/Microphysics/Microphysics.H @@ -3,6 +3,7 @@ #include #include +#include class Microphysics { diff --git a/Source/Microphysics/Null/NullMoist.H b/Source/Microphysics/Null/NullMoist.H index 994654512..1484058fa 100644 --- a/Source/Microphysics/Null/NullMoist.H +++ b/Source/Microphysics/Null/NullMoist.H @@ -13,24 +13,24 @@ class NullMoist { virtual ~NullMoist () = default; virtual void - define (SolverChoice& /*sc*/) { }; + define (SolverChoice& /*sc*/) { } virtual void Init (const amrex::MultiFab& /*cons_in*/, amrex::MultiFab& /*qmoist*/, const amrex::BoxArray& /*grids*/, const amrex::Geometry& /*geom*/, - const amrex::Real& /*dt_advance*/) { }; + const amrex::Real& /*dt_advance*/) { } virtual void - Advance ( ) { }; + Advance ( ) { } virtual void - Diagnose ( ) { }; + Diagnose ( ) { } virtual void Update (amrex::MultiFab& /*cons_in*/, - amrex::MultiFab& /*qmoist*/) { }; + amrex::MultiFab& /*qmoist*/) { } private: diff --git a/Source/Microphysics/SAM/SAM.H b/Source/Microphysics/SAM/SAM.H index 28d1636c2..e9d263fe2 100644 --- a/Source/Microphysics/SAM/SAM.H +++ b/Source/Microphysics/SAM/SAM.H @@ -75,12 +75,12 @@ class SAM : public NullMoist { // micro interface for precip fall void MicroPrecipFall (); - // diagnose - void Diagnose (); - // process microphysics void Proc (); + // diagnose + void Diagnose () override; + // Set up for first time void define (SolverChoice& sc) override diff --git a/Source/TimeIntegration/ERF_make_buoyancy.cpp b/Source/TimeIntegration/ERF_make_buoyancy.cpp index 6e3b7f39e..b5c320621 100644 --- a/Source/TimeIntegration/ERF_make_buoyancy.cpp +++ b/Source/TimeIntegration/ERF_make_buoyancy.cpp @@ -203,7 +203,7 @@ void make_buoyancy (Vector& S_data, Real* rho_d_ptr = rho_d.data(); Real* theta_d_ptr = theta_d.data(); - if (solverChoice.buoyancy_type == 2) { + if (solverChoice.buoyancy_type == 2 || solverChoice.buoyancy_type == 4 ) { prim_ave.line_average(PrimQp_comp, qp_h); @@ -245,17 +245,37 @@ void make_buoyancy (Vector& S_data, Real tempp3d = getTgivenRandRTh(cell_data(i,j,k ,Rho_comp), cell_data(i,j,k ,RhoTheta_comp)); Real tempm3d = getTgivenRandRTh(cell_data(i,j,k-1,Rho_comp), cell_data(i,j,k-1,RhoTheta_comp)); - Real qplus = 0.61* ( qv_data(i,j,k)-qv_d_ptr[k]) - + Real qplus, qminus; + + if(solverChoice.buoyancy_type == 2){ + qplus = 0.61* ( qv_data(i,j,k)-qv_d_ptr[k]) - (qc_data(i,j,k)-qc_d_ptr[k]+ qi_data(i,j,k)-qi_d_ptr[k]+ cell_prim(i,j,k,PrimQp_comp)-qp_d_ptr[k]) + (tempp3d-tempp1d)/tempp1d*(Real(1.0) + Real(0.61)*qv_d_ptr[k]-qc_d_ptr[k]-qi_d_ptr[k]-qp_d_ptr[k]); - Real qminus = 0.61 *( qv_data(i,j,k-1)-qv_d_ptr[k-1]) - + qminus = 0.61 *( qv_data(i,j,k-1)-qv_d_ptr[k-1]) - (qc_data(i,j,k-1)-qc_d_ptr[k-1]+ qi_data(i,j,k-1)-qi_d_ptr[k-1]+ cell_prim(i,j,k-1,PrimQp_comp)-qp_d_ptr[k-1]) + (tempm3d-tempm1d)/tempm1d*(Real(1.0) + Real(0.61)*qv_d_ptr[k-1]-qi_d_ptr[k-1]-qc_d_ptr[k-1]-qp_d_ptr[k-1]); + } + + if(solverChoice.buoyancy_type == 4){ + qplus = 0.61* ( qv_data(i,j,k)-qv_d_ptr[k]) - + (qc_data(i,j,k)-qc_d_ptr[k]+ + qi_data(i,j,k)-qi_d_ptr[k]+ + cell_prim(i,j,k,PrimQp_comp)-qp_d_ptr[k]) + + (cell_data(i,j,k ,RhoTheta_comp)/cell_data(i,j,k ,Rho_comp) - theta_d_ptr[k ])/theta_d_ptr[k ]; + + qminus = 0.61 *( qv_data(i,j,k-1)-qv_d_ptr[k-1]) - + (qc_data(i,j,k-1)-qc_d_ptr[k-1]+ + qi_data(i,j,k-1)-qi_d_ptr[k-1]+ + cell_prim(i,j,k-1,PrimQp_comp)-qp_d_ptr[k-1]) + + (cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp) - theta_d_ptr[k-1])/theta_d_ptr[k-1]; + } + + Real qavg = Real(0.5) * (qplus + qminus); Real r0avg = Real(0.5) * (rho_d_ptr[k] + rho_d_ptr[k-1]); @@ -305,7 +325,7 @@ void make_buoyancy (Vector& S_data, buoyancy_fab(i, j, k) = -qavg * r0avg * grav_gpu[2]; }); } // mfi - } // buoyancy_type + } // buoyancy_type } // not buoyancy_type == 1 #endif diff --git a/Source/prob_common.H b/Source/prob_common.H index 40be9eda2..931f04479 100644 --- a/Source/prob_common.H +++ b/Source/prob_common.H @@ -49,6 +49,15 @@ public: amrex::Error("Should never call erf_init_dens_hse for "+name()+" problem"); } + virtual void + erf_init_dens_hse_moist (amrex::MultiFab& /*rho_hse*/, + std::unique_ptr& /*z_phys_nd*/, + std::unique_ptr& /*z_phys_cc*/, + amrex::Geometry const& /*geom*/) + { + + } + /** * Function to perform custom initialization of a test problem *