diff --git a/Build/cmake_with_warm_no_precip.sh b/Build/cmake_with_warm_no_precip.sh deleted file mode 100755 index cf3c33526..000000000 --- a/Build/cmake_with_warm_no_precip.sh +++ /dev/null @@ -1,18 +0,0 @@ -#!/bin/bash - -# Example CMake config script for an OSX laptop with OpenMPI - -cmake -DCMAKE_INSTALL_PREFIX:PATH=./install \ - -DCMAKE_CXX_COMPILER:STRING=mpicxx \ - -DCMAKE_C_COMPILER:STRING=mpicc \ - -DCMAKE_Fortran_COMPILER:STRING=mpifort \ - -DMPIEXEC_PREFLAGS:STRING=--oversubscribe \ - -DCMAKE_BUILD_TYPE:STRING=Release \ - -DERF_DIM:STRING=3 \ - -DERF_ENABLE_MPI:BOOL=ON \ - -DERF_ENABLE_TESTS:BOOL=ON \ - -DERF_ENABLE_FCOMPARE:BOOL=ON \ - -DERF_ENABLE_DOCUMENTATION:BOOL=OFF \ - -DERF_ENABLE_WARM_NO_PRECIP:BOOL=ON \ - -DCMAKE_EXPORT_COMPILE_COMMANDS:BOOL=ON \ - .. && make -j8 diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake index b64704f24..442225e0c 100644 --- a/CMake/BuildERFExe.cmake +++ b/CMake/BuildERFExe.cmake @@ -136,7 +136,6 @@ function(build_erf_lib erf_lib_name) ${SRC_DIR}/TimeIntegration/ERF_advance_microphysics.cpp ${SRC_DIR}/TimeIntegration/ERF_advance_radiation.cpp ${SRC_DIR}/TimeIntegration/ERF_make_buoyancy.cpp - ${SRC_DIR}/TimeIntegration/ERF_make_condensation_source.cpp ${SRC_DIR}/TimeIntegration/ERF_make_fast_coeffs.cpp ${SRC_DIR}/TimeIntegration/ERF_slow_rhs_pre.cpp ${SRC_DIR}/TimeIntegration/ERF_ApplySpongeZoneBCs.cpp diff --git a/Exec/ABL/prob.cpp b/Exec/ABL/prob.cpp index 84724b1e3..bf767fe96 100644 --- a/Exec/ABL/prob.cpp +++ b/Exec/ABL/prob.cpp @@ -63,8 +63,10 @@ Problem::init_custom_pert( amrex::Array4 const& /*mf_m*/, amrex::Array4 const& /*mf_u*/, amrex::Array4 const& /*mf_v*/, - const SolverChoice& /*sc*/) + const SolverChoice& sc) { + const bool use_moisture = (sc.moisture_type != MoistureType::None); + 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(); @@ -93,8 +95,10 @@ Problem::init_custom_pert( // Set an initial value for QKE state(i, j, k, RhoQKE_comp) = parms.QKE_0; - state(i, j, k, RhoQ1_comp) = 0.0; - state(i, j, k, RhoQ2_comp) = 0.0; + if (use_moisture) { + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; + } }); // Set the x-velocity diff --git a/Exec/ABL_input_sounding/prob.cpp b/Exec/ABL_input_sounding/prob.cpp index 09c3d0090..80911c884 100644 --- a/Exec/ABL_input_sounding/prob.cpp +++ b/Exec/ABL_input_sounding/prob.cpp @@ -62,8 +62,10 @@ Problem::init_custom_pert( amrex::Array4 const& /*mf_m*/, amrex::Array4 const& /*mf_u*/, amrex::Array4 const& /*mf_v*/, - const SolverChoice& /*sc*/) + const SolverChoice& sc) { + const bool use_moisture = (sc.moisture_type != MoistureType::None); + ParallelFor(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { // Geometry const Real* prob_lo = geomdata.ProbLo(); @@ -86,8 +88,10 @@ Problem::init_custom_pert( // Set an initial value for QKE state(i, j, k, RhoQKE_comp) = parms.QKE_0; - state(i, j, k, RhoQ1_comp) = 0.0; - state(i, j, k, RhoQ2_comp) = 0.0; + if (use_moisture) { + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; + } }); // Set the x-velocity diff --git a/Exec/DevTests/MovingTerrain/prob.cpp b/Exec/DevTests/MovingTerrain/prob.cpp index fe467cab9..d4c5181bc 100644 --- a/Exec/DevTests/MovingTerrain/prob.cpp +++ b/Exec/DevTests/MovingTerrain/prob.cpp @@ -43,10 +43,12 @@ Problem::init_custom_pert( Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, Array4 const& /*mf_v*/, - const SolverChoice&) + const SolverChoice& sc) { const int khi = geomdata.Domain().bigEnd()[2]; + const bool use_moisture = (sc.moisture_type != MoistureType::None); + AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1); const Real T_sfc = parms.T_0; @@ -83,9 +85,10 @@ Problem::init_custom_pert( // Set scalar = 0 everywhere state(i, j, k, RhoScalar_comp) = state(i,j,k,Rho_comp); - state(i, j, k, RhoQ1_comp) = 0.0; - state(i, j, k, RhoQ2_comp) = 0.0; - + if (use_moisture) { + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; + } }); // Set the x-velocity diff --git a/Exec/DevTests/MultiBlock/prob.cpp b/Exec/DevTests/MultiBlock/prob.cpp index 4607e45bf..26287abe0 100644 --- a/Exec/DevTests/MultiBlock/prob.cpp +++ b/Exec/DevTests/MultiBlock/prob.cpp @@ -51,6 +51,8 @@ Problem::init_custom_pert( { const int khi = geomdata.Domain().bigEnd()[2]; + const bool use_moisture = (sc.moisture_type != MoistureType::None); + AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1); const Real rho_sfc = p_0 / (R_d*parms.T_0); @@ -86,8 +88,10 @@ Problem::init_custom_pert( // Set scalar = 0 everywhere state(i, j, k, RhoScalar_comp) = 0.0; - state(i, j, k, RhoQ1_comp) = 0.0; - state(i, j, k, RhoQ2_comp) = 0.0; + if (use_moisture) { + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; + } }); // Set the x-velocity diff --git a/Exec/DevTests/ParticlesOverWoA/prob.cpp b/Exec/DevTests/ParticlesOverWoA/prob.cpp index 8d52ee3ca..e81c90a14 100644 --- a/Exec/DevTests/ParticlesOverWoA/prob.cpp +++ b/Exec/DevTests/ParticlesOverWoA/prob.cpp @@ -48,10 +48,12 @@ Problem::init_custom_pert( Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, Array4 const& /*mf_v*/, - const SolverChoice&) + const SolverChoice& sc) { const int khi = geomdata.Domain().bigEnd()[2]; + const bool use_moisture = (sc.moisture_type != MoistureType::None); + AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1); const Real rho_sfc = p_0 / (R_d*parms.T_0); @@ -85,8 +87,10 @@ Problem::init_custom_pert( // Set scalar = 0 everywhere state(i, j, k, RhoScalar_comp) = 0.0; - state(i, j, k, RhoQ1_comp) = 0.0; - state(i, j, k, RhoQ2_comp) = 0.0; + if (use_moisture) { + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; + } }); // Set the x-velocity diff --git a/Exec/Radiation/prob.cpp b/Exec/Radiation/prob.cpp index a1e175746..2016fa9d9 100644 --- a/Exec/Radiation/prob.cpp +++ b/Exec/Radiation/prob.cpp @@ -105,14 +105,16 @@ Problem::init_custom_pert( Array4 const& /*mf_v*/, const SolverChoice& sc) { - const int khi = geomdata.Domain().bigEnd()[2]; + const int khi = geomdata.Domain().bigEnd()[2]; - AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1); + const bool use_moisture = (sc.moisture_type != MoistureType::None); - // 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]; + 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; @@ -182,15 +184,17 @@ Problem::init_custom_pert( state(i, j, k, RhoScalar_comp) = 0.0; // mean states - state(i, j, k, RhoQ1_comp) = rho*qvapor; - state(i, j, k, RhoQ2_comp) = 0.0; - qv(i, j, k) = qvapor; - qc(i, j, k) = 0.0; - qi(i, j, k) = 0.0; + if (use_moisture) { + state(i, j, k, RhoQ1_comp) = rho*qvapor; + state(i, j, k, RhoQ2_comp) = 0.0; + qv(i, j, k) = qvapor; + qc(i, j, k) = 0.0; + qi(i, j, k) = 0.0; #if defined(ERF_USE_WARM_NO_PRECIP) - state(i, j, k, RhoQ1_comp) = rho*qvapor; - state(i, j, k, RhoQ2_comp) = 0.0; + state(i, j, k, RhoQ1_comp) = rho*qvapor; + state(i, j, k, RhoQ2_comp) = 0.0; #endif + } }); // Set the x-velocity diff --git a/Exec/RegTests/Bubble/prob.cpp b/Exec/RegTests/Bubble/prob.cpp index b01a5e963..1056f6e44 100644 --- a/Exec/RegTests/Bubble/prob.cpp +++ b/Exec/RegTests/Bubble/prob.cpp @@ -58,6 +58,8 @@ Problem::init_custom_pert( { const int khi = geomdata.Domain().bigEnd()[2]; + const bool use_moisture = (sc.moisture_type != MoistureType::None); + AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1); const Real rho_sfc = p_0 / (R_d*parms.T_0); @@ -131,8 +133,10 @@ Problem::init_custom_pert( state(i, j, k, RhoScalar_comp) = 0.0; - state(i, j, k, RhoQ1_comp) = 0.0; - state(i, j, k, RhoQ2_comp) = 0.0; + if (use_moisture) { + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; + } }); } else { @@ -177,8 +181,10 @@ Problem::init_custom_pert( state(i, j, k, RhoScalar_comp) = 0.0; - state(i, j, k, RhoQ1_comp) = 0.0; - state(i, j, k, RhoQ2_comp) = 0.0; + if (use_moisture) { + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; + } }); } diff --git a/Exec/RegTests/CouetteFlow/prob.cpp b/Exec/RegTests/CouetteFlow/prob.cpp index bd088184f..a160922ad 100644 --- a/Exec/RegTests/CouetteFlow/prob.cpp +++ b/Exec/RegTests/CouetteFlow/prob.cpp @@ -44,15 +44,19 @@ Problem::init_custom_pert( Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, Array4 const& /*mf_v*/, - const SolverChoice&) + const SolverChoice& sc) { + const bool use_moisture = (sc.moisture_type != MoistureType::None); + amrex::ParallelFor(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { // Set scalar to 0 state(i, j, k, RhoScalar_comp) = 0.0; - state(i, j, k, RhoQ1_comp) = 0.0; - state(i, j, k, RhoQ2_comp) = 0.0; + if (use_moisture) { + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; + } }); diff --git a/Exec/RegTests/DensityCurrent/prob.cpp b/Exec/RegTests/DensityCurrent/prob.cpp index 87af727b5..751b31ad1 100644 --- a/Exec/RegTests/DensityCurrent/prob.cpp +++ b/Exec/RegTests/DensityCurrent/prob.cpp @@ -49,19 +49,21 @@ Problem::init_custom_pert( Array4 const& /*mf_v*/, const SolverChoice& sc) { - const int khi = geomdata.Domain().bigEnd()[2]; + const int khi = geomdata.Domain().bigEnd()[2]; - AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1); + const bool use_moisture = (sc.moisture_type != MoistureType::None); - const Real l_x_r = parms.x_r; - //const Real l_x_r = parms.x_r * mf_u(0,0,0); //used to validate constant msf - const Real l_z_r = parms.z_r; - const Real l_x_c = parms.x_c; - const Real l_z_c = parms.z_c; - const Real l_Tpt = parms.T_pert; - const Real rdOcp = sc.rdOcp; + AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1); - if (z_cc) { + const Real l_x_r = parms.x_r; + //const Real l_x_r = parms.x_r * mf_u(0,0,0); //used to validate constant msf + const Real l_z_r = parms.z_r; + const Real l_x_c = parms.x_c; + const Real l_z_c = parms.z_c; + const Real l_Tpt = parms.T_pert; + const Real rdOcp = sc.rdOcp; + + if (z_cc) { amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { // Geometry (note we must include these here to get the data on device) @@ -88,8 +90,10 @@ Problem::init_custom_pert( // Set scalar = 0 everywhere state(i, j, k, RhoScalar_comp) = 0.0; - state(i, j, k, RhoQ1_comp) = 0.0; - state(i, j, k, RhoQ2_comp) = 0.0; + if (use_moisture) { + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; + } }); } else { amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept @@ -118,8 +122,10 @@ Problem::init_custom_pert( // Set scalar = 0 everywhere state(i, j, k, RhoScalar_comp) = 0.0; - state(i, j, k, RhoQ1_comp) = 0.0; - state(i, j, k, RhoQ2_comp) = 0.0; + if (use_moisture) { + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; + } }); } diff --git a/Exec/RegTests/DynamicRefinement/prob.cpp b/Exec/RegTests/DynamicRefinement/prob.cpp index 76cbb4bae..fa911cc0d 100644 --- a/Exec/RegTests/DynamicRefinement/prob.cpp +++ b/Exec/RegTests/DynamicRefinement/prob.cpp @@ -93,6 +93,7 @@ Problem::init_custom_pert( Array4 const& /*mf_v*/, const SolverChoice& sc) { + const bool use_moisture = (sc.moisture_type != MoistureType::None); Real xc = parms.xc; Real yc = parms.yc; Real R = parms.R ; Real beta = parms.beta; @@ -126,8 +127,10 @@ Problem::init_custom_pert( // Set scalar = 0 -- unused state(i, j, k, RhoScalar_comp) = 1.0 + Omg; - state(i, j, k, RhoQ1_comp) = 0.0; - state(i, j, k, RhoQ2_comp) = 0.0; + if (use_moisture) { + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; + } }); // Set the x-velocity diff --git a/Exec/RegTests/EkmanSpiral_custom/prob.cpp b/Exec/RegTests/EkmanSpiral_custom/prob.cpp index c25554060..1c9c8810b 100644 --- a/Exec/RegTests/EkmanSpiral_custom/prob.cpp +++ b/Exec/RegTests/EkmanSpiral_custom/prob.cpp @@ -41,8 +41,10 @@ Problem::init_custom_pert( Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, Array4 const& /*mf_v*/, - const SolverChoice&) + const SolverChoice& sc) { + const bool use_moisture = (sc.moisture_type != MoistureType::None); + ParallelFor(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { // Geometry //const Real* prob_lo = geomdata.ProbLo(); @@ -54,8 +56,10 @@ Problem::init_custom_pert( // Set scalar = 0 everywhere state(i, j, k, RhoScalar_comp) = 0.0; - state(i, j, k, RhoQ1_comp) = 0.0; - state(i, j, k, RhoQ2_comp) = 0.0; + if (use_moisture) { + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; + } }); diff --git a/Exec/RegTests/GABLS1/prob.cpp b/Exec/RegTests/GABLS1/prob.cpp index 74da64e9b..2ea9dc86e 100644 --- a/Exec/RegTests/GABLS1/prob.cpp +++ b/Exec/RegTests/GABLS1/prob.cpp @@ -63,8 +63,10 @@ Problem::init_custom_pert( amrex::Array4 const& /*mf_m*/, amrex::Array4 const& /*mf_u*/, amrex::Array4 const& /*mf_v*/, - const SolverChoice& /*sc*/) + const SolverChoice& sc) { + const bool use_moisture = (sc.moisture_type != MoistureType::None); + if (state.nComp() > RhoQKE_comp) { amrex::Print() << "Initializing QKE" << std::endl; } diff --git a/Exec/RegTests/IsentropicVortex/prob.cpp b/Exec/RegTests/IsentropicVortex/prob.cpp index e011626d6..52442eae3 100644 --- a/Exec/RegTests/IsentropicVortex/prob.cpp +++ b/Exec/RegTests/IsentropicVortex/prob.cpp @@ -95,6 +95,7 @@ Problem::init_custom_pert( Array4 const& /*mf_v*/, const SolverChoice& sc) { + const bool use_moisture = (sc.moisture_type != MoistureType::None); Real xc = parms.xc; Real yc = parms.yc; Real R = parms.R ; Real beta = parms.beta; @@ -128,8 +129,10 @@ Problem::init_custom_pert( // Set scalar = 0 -- unused state(i, j, k, RhoScalar_comp) = 0.0; - state(i, j, k, RhoQ1_comp) = 0.0; - state(i, j, k, RhoQ2_comp) = 0.0; + if (use_moisture) { + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; + } }); // Set the x-velocity diff --git a/Exec/RegTests/PoiseuilleFlow/prob.cpp b/Exec/RegTests/PoiseuilleFlow/prob.cpp index 9302e2f2f..f2b2d3d09 100644 --- a/Exec/RegTests/PoiseuilleFlow/prob.cpp +++ b/Exec/RegTests/PoiseuilleFlow/prob.cpp @@ -43,8 +43,10 @@ Problem::init_custom_pert( Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, Array4 const& /*mf_v*/, - const SolverChoice&) + const SolverChoice& sc) { + const bool use_moisture = (sc.moisture_type != MoistureType::None); + amrex::ParallelFor(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { // Arbitrarily choose the radius of the bubble to be 0.05 times the length of the domain @@ -52,8 +54,10 @@ Problem::init_custom_pert( // Set scalar = A_0*exp(-10r^2), where r is distance from center of domain state(i, j, k, RhoScalar_comp) = 0.0; - state(i, j, k, RhoQ1_comp) = 0.0; - state(i, j, k, RhoQ2_comp) = 0.0; + if (use_moisture) { + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; + } }); diff --git a/Exec/RegTests/ScalarAdvDiff/prob.cpp b/Exec/RegTests/ScalarAdvDiff/prob.cpp index b2df00229..1ea764346 100644 --- a/Exec/RegTests/ScalarAdvDiff/prob.cpp +++ b/Exec/RegTests/ScalarAdvDiff/prob.cpp @@ -77,8 +77,10 @@ Problem::init_custom_pert( Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, Array4 const& /*mf_v*/, - const SolverChoice&) + const SolverChoice& sc) { + const bool use_moisture = (sc.moisture_type != MoistureType::None); + amrex::ParallelFor(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { // Geometry @@ -129,8 +131,10 @@ Problem::init_custom_pert( state(i, j, k, RhoScalar_comp) *= parms.rho_0; - state(i, j, k, RhoQ1_comp) = 0.0; - state(i, j, k, RhoQ2_comp) = 0.0; + if (use_moisture) { + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; + } }); // Set the x-velocity diff --git a/Exec/RegTests/StokesSecondProblem/prob.cpp b/Exec/RegTests/StokesSecondProblem/prob.cpp index 80ca94bd9..11cbf0263 100644 --- a/Exec/RegTests/StokesSecondProblem/prob.cpp +++ b/Exec/RegTests/StokesSecondProblem/prob.cpp @@ -39,10 +39,12 @@ Problem::init_custom_pert( Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, Array4 const& /*mf_v*/, - const SolverChoice&) + const SolverChoice& sc) { const int khi = geomdata.Domain().bigEnd()[2]; + const bool use_moisture = (sc.moisture_type != MoistureType::None); + AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1); // Geometry (note we must include these here to get the data on device) @@ -54,8 +56,10 @@ Problem::init_custom_pert( // Set scalar = 0 everywhere state(i, j, k, RhoScalar_comp) = 0.0; - state(i, j, k, RhoQ1_comp) = 0.0; - state(i, j, k, RhoQ2_comp) = 0.0; + if (use_moisture) { + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; + } }); // Set the x-velocity diff --git a/Exec/RegTests/TaylorGreenVortex/prob.cpp b/Exec/RegTests/TaylorGreenVortex/prob.cpp index 4b55a5a23..6f244b8f2 100644 --- a/Exec/RegTests/TaylorGreenVortex/prob.cpp +++ b/Exec/RegTests/TaylorGreenVortex/prob.cpp @@ -43,8 +43,10 @@ Problem::init_custom_pert( Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, Array4 const& /*mf_v*/, - const SolverChoice&) + const SolverChoice& sc) { + const bool use_moisture = (sc.moisture_type != MoistureType::None); + ParallelFor(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { // Geometry const Real* prob_lo = geomdata.ProbLo(); @@ -64,8 +66,10 @@ Problem::init_custom_pert( // Set scalar = 0 everywhere state(i, j, k, RhoScalar_comp) = 1.0 * parms.rho_0; - state(i, j, k, RhoQ1_comp) = 0.0; - state(i, j, k, RhoQ2_comp) = 0.0; + if (use_moisture) { + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; + } }); // Set the x-velocity diff --git a/Exec/RegTests/Terrain2d_Cylinder/prob.cpp b/Exec/RegTests/Terrain2d_Cylinder/prob.cpp index 5bb9cb8ac..918a51ddd 100644 --- a/Exec/RegTests/Terrain2d_Cylinder/prob.cpp +++ b/Exec/RegTests/Terrain2d_Cylinder/prob.cpp @@ -43,10 +43,12 @@ Problem::init_custom_pert( Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, Array4 const& /*mf_v*/, - const SolverChoice&) + const SolverChoice& sc) { const int khi = geomdata.Domain().bigEnd()[2]; + const bool use_moisture = (sc.moisture_type != MoistureType::None); + AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1); // Geometry (note we must include these here to get the data on device) @@ -60,8 +62,10 @@ Problem::init_custom_pert( // Set scalar = 0 everywhere state(i, j, k, RhoScalar_comp) = 0.0; - state(i, j, k, RhoQ1_comp) = 0.0; - state(i, j, k, RhoQ2_comp) = 0.0; + if (use_moisture) { + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; + } }); // Set the x-velocity diff --git a/Exec/RegTests/Terrain3d_Hemisphere/prob.cpp b/Exec/RegTests/Terrain3d_Hemisphere/prob.cpp index 3283b9af7..ee58919eb 100644 --- a/Exec/RegTests/Terrain3d_Hemisphere/prob.cpp +++ b/Exec/RegTests/Terrain3d_Hemisphere/prob.cpp @@ -42,10 +42,12 @@ Problem::init_custom_pert( Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, Array4 const& /*mf_v*/, - const SolverChoice&) + const SolverChoice& sc) { const int khi = geomdata.Domain().bigEnd()[2]; + const bool use_moisture = (sc.moisture_type != MoistureType::None); + AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1); // Geometry (note we must include these here to get the data on device) @@ -59,8 +61,10 @@ Problem::init_custom_pert( // Set scalar = 0 everywhere state(i, j, k, RhoScalar_comp) = 0.0; - state(i, j, k, RhoQ1_comp) = 0.0; - state(i, j, k, RhoQ2_comp) = 0.0; + if (use_moisture) { + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; + } }); // Set the x-velocity diff --git a/Exec/RegTests/WitchOfAgnesi/prob.cpp b/Exec/RegTests/WitchOfAgnesi/prob.cpp index 04a2e8f6b..56c970c6f 100644 --- a/Exec/RegTests/WitchOfAgnesi/prob.cpp +++ b/Exec/RegTests/WitchOfAgnesi/prob.cpp @@ -41,10 +41,12 @@ Problem::init_custom_pert( Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, Array4 const& /*mf_v*/, - const SolverChoice&) + const SolverChoice& sc) { const int khi = geomdata.Domain().bigEnd()[2]; + const bool use_moisture = (sc.moisture_type != MoistureType::None); + AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1); amrex::ParallelFor(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept @@ -58,8 +60,10 @@ Problem::init_custom_pert( // Set scalar = 0 everywhere state(i, j, k, RhoScalar_comp) = 0.0; - state(i, j, k, RhoQ1_comp) = 0.0; - state(i, j, k, RhoQ2_comp) = 0.0; + if (use_moisture) { + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; + } }); // Set the x-velocity diff --git a/Exec/SquallLine_2D/prob.cpp b/Exec/SquallLine_2D/prob.cpp index e9e90c08d..6d36b971e 100644 --- a/Exec/SquallLine_2D/prob.cpp +++ b/Exec/SquallLine_2D/prob.cpp @@ -291,9 +291,11 @@ Problem::init_custom_pert ( Array4 const& /*mf_v*/, const SolverChoice& sc) { - const int khi = geomdata.Domain().bigEnd()[2]; + const bool use_moisture = (sc.moisture_type != MoistureType::None); - AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1); + 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]; @@ -375,15 +377,17 @@ Problem::init_custom_pert ( state(i, j, k, RhoScalar_comp) = rho*scalar; // mean states - state(i, j, k, RhoQ1_comp) = rho*q_v_hot; - state(i, j, k, RhoQ2_comp) = 0.0; - qv(i, j, k) = q_v_hot; - qc(i, j, k) = 0.0; - qi(i, j, k) = 0.0; + if (use_moisture) { + state(i, j, k, RhoQ1_comp) = rho*q_v_hot; + state(i, j, k, RhoQ2_comp) = 0.0; + qv(i, j, k) = q_v_hot; + qc(i, j, k) = 0.0; + qi(i, j, k) = 0.0; #if defined(ERF_USE_WARM_NO_PRECIP) - state(i, j, k, RhoQv_comp) = 0.0;//rho*qvapor; - state(i, j, k, RhoQc_comp) = 0.0; + state(i, j, k, RhoQv_comp) = 0.0;//rho*qvapor; + state(i, j, k, RhoQc_comp) = 0.0; #endif + } }); diff --git a/Exec/SuperCell/prob.cpp b/Exec/SuperCell/prob.cpp index a1e175746..e7186c3f6 100644 --- a/Exec/SuperCell/prob.cpp +++ b/Exec/SuperCell/prob.cpp @@ -105,9 +105,11 @@ Problem::init_custom_pert( Array4 const& /*mf_v*/, const SolverChoice& sc) { - const int khi = geomdata.Domain().bigEnd()[2]; + const bool use_moisture = (sc.moisture_type != MoistureType::None); - AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1); + 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]; @@ -182,15 +184,17 @@ Problem::init_custom_pert( state(i, j, k, RhoScalar_comp) = 0.0; // mean states - state(i, j, k, RhoQ1_comp) = rho*qvapor; - state(i, j, k, RhoQ2_comp) = 0.0; - qv(i, j, k) = qvapor; - qc(i, j, k) = 0.0; - qi(i, j, k) = 0.0; + if (use_moisture) { + state(i, j, k, RhoQ1_comp) = rho*qvapor; + state(i, j, k, RhoQ2_comp) = 0.0; + qv(i, j, k) = qvapor; + qc(i, j, k) = 0.0; + qi(i, j, k) = 0.0; #if defined(ERF_USE_WARM_NO_PRECIP) - state(i, j, k, RhoQ1_comp) = rho*qvapor; - state(i, j, k, RhoQ2_comp) = 0.0; + state(i, j, k, RhoQ1_comp) = rho*qvapor; + state(i, j, k, RhoQ2_comp) = 0.0; #endif + } }); // Set the x-velocity diff --git a/Source/DataStructs/DataStruct.H b/Source/DataStructs/DataStruct.H index 80540afbc..01276a7b8 100644 --- a/Source/DataStructs/DataStruct.H +++ b/Source/DataStructs/DataStruct.H @@ -66,12 +66,30 @@ struct SolverChoice { // Do we set map scale factors to 0.5 instead of 1 for testing? pp.query("test_mapfactor", test_mapfactor); + // What type of moisture model to use + static std::string moisture_model_string = "None"; + pp.query("moisture_model", moisture_model_string); + if (moisture_model_string == "SAM") { + moisture_type = MoistureType::SAM; + } else if (moisture_model_string == "Kessler") { + moisture_type = MoistureType::Kessler; + } else if (moisture_model_string == "FastEddy") { + moisture_type = MoistureType::FastEddy; + } else { + moisture_type = MoistureType::None; + } + // 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 && buoyancy_type != 4) { amrex::Abort("buoyancy_type must be 1, 2, 3 or 4"); } + // Set a different default for moist vs dry + if (moisture_type != MoistureType::None) { + buoyancy_type = 2; // uses Tprime + } + // Is the terrain static or moving? static std::string terrain_type_string = "Static"; pp.query("terrain_type",terrain_type_string); @@ -83,19 +101,6 @@ struct SolverChoice { amrex::Abort("terrain_type can be either Moving/moving or Static/static"); } - // What type of moisture model to use - static std::string moisture_model_string = "None"; - pp.query("moisture_model", moisture_model_string); - if (moisture_model_string == "SAM") { - moisture_type = MoistureType::SAM; - } else if (moisture_model_string == "Kessler") { - moisture_type = MoistureType::Kessler; - } else if (moisture_model_string == "FastEddy") { - moisture_type = MoistureType::FastEddy; - } else { - moisture_type = MoistureType::None; - } - // Use lagged_delta_rt in the fast integrator? pp.query("use_lagged_delta_rt", use_lagged_delta_rt); @@ -110,9 +115,6 @@ struct SolverChoice { pp.query("c_p", c_p); rdOcp = R_d / c_p; -#if defined(ERF_USE_WARM_NO_PRECIP) - pp.query("tau_cond", tau_cond); -#endif #if defined(ERF_USE_POISSON_SOLVE) // Should we project the initial velocity field to make it divergence-free? @@ -317,11 +319,7 @@ struct SolverChoice { bool test_mapfactor = false; bool use_terrain = false; -#ifdef ERF_USE_MOISTURE - int buoyancy_type = 2; // uses Tprime -#else int buoyancy_type = 1; // uses rhoprime directly -#endif // Specify what additional physics/forcing modules we use bool use_gravity = false; @@ -346,10 +344,6 @@ struct SolverChoice { amrex::Real zsurf = 0.0; amrex::Real dz0; -#if defined(ERF_USE_WARM_NO_PRECIP) - amrex::Real tau_cond = 1.0; // Default time of 1 sec -- this is somewhat arbitray -#endif - #if defined(ERF_USE_POISSON_SOLVE) int project_initial_velocity = 1; #endif diff --git a/Source/ERF.H b/Source/ERF.H index 5ca3baadf..8e3735709 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -383,11 +383,6 @@ private: // includes a recursive call for finer levels void timeStep (int lev, amrex::Real time, int iteration); -#if defined(ERF_USE_WARM_NO_PRECIP) - void condensation_source (amrex::MultiFab& source, amrex::MultiFab& S_new, - amrex::Real tau_cond, amrex::Real c_p); -#endif - // advance a single level for a single time step void Advance (int lev, amrex::Real time, amrex::Real dt_lev, int iteration, int ncycle); @@ -530,7 +525,7 @@ private: amrex::Vector rW_new; Microphysics micro; - amrex::Vector qmoist; // This has 6 components: qv, qc, qi, qr, qs, qg + amrex::Vector> qmoist; // This has up to 6 components: qv, qc, qi, qr, qs, qg #if defined(ERF_USE_RRTMGP) Radiation rad; diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 22a35aa68..4bd5a6aff 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -548,7 +548,7 @@ ERF::InitData () if (solverChoice.moisture_type != MoistureType::None) { for (int lev = 0; lev <= finest_level; lev++) { - FillPatchMoistVars(lev, qmoist[lev]); + FillPatchMoistVars(lev, *(qmoist[lev])); } } } @@ -588,9 +588,9 @@ ERF::InitData () { // If not restarting we need to fill qmoist given qt and qp. if (restart_chkfile.empty()) { - micro.Init(vars_new[lev][Vars::cons], qmoist[lev], + micro.Init(vars_new[lev][Vars::cons], *(qmoist[lev]), grids[lev], Geom(lev), 0.0); // dummy value, not needed just to diagnose - micro.Update(vars_new[lev][Vars::cons], qmoist[lev]); + micro.Update(vars_new[lev][Vars::cons], *(qmoist[lev])); } } } @@ -859,6 +859,8 @@ ERF::init_only (int lev, Real time) lev_new[Vars::yvel].setVal(0.0); lev_old[Vars::yvel].setVal(0.0); lev_new[Vars::zvel].setVal(0.0); lev_old[Vars::zvel].setVal(0.0); + qmoist[lev]->setVal(0.); + // Initialize background flow (optional) if (init_type == "input_sounding") { // The base state is initialized by integrating vertically through the @@ -892,8 +894,6 @@ ERF::init_only (int lev, Real time) init_from_hse(lev); } - qmoist[lev].setVal(0.); - // Add problem-specific flow features // // Notes: @@ -1175,7 +1175,7 @@ ERF::MakeHorizontalAverages () if (use_moisture) { - MultiFab qv(qmoist[lev], make_alias, 0, 1); + MultiFab qv(*(qmoist[lev]), make_alias, 0, 1); for (MFIter mfi(mf); mfi.isValid(); ++mfi) { const Box& bx = mfi.validbox(); diff --git a/Source/ERF_make_new_level.cpp b/Source/ERF_make_new_level.cpp index a656f87c0..cbd27688e 100644 --- a/Source/ERF_make_new_level.cpp +++ b/Source/ERF_make_new_level.cpp @@ -61,7 +61,8 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba, //******************************************************************************************** // Microphysics // ******************************************************************************************* - qmoist[lev].define(ba, dm, 6, ngrow_state); // qv, qc, qi, qr, qs, qg + int q_size = micro.Get_Qmoist_Size(); + qmoist[lev] = std::make_unique(ba, dm, q_size, ngrow_state); // qv, qc, qi, qr, qs, qg // ******************************************************************************************** // Build the data structures for calculating diffusive/turbulent terms @@ -198,7 +199,8 @@ ERF::MakeNewLevelFromCoarse (int lev, Real time, const BoxArray& ba, // Microphysics // ******************************************************************************************* int ngrow_state = ComputeGhostCells(solverChoice.advChoice, solverChoice.use_NumDiff) + 1; - qmoist[lev].define(ba, dm, 6, ngrow_state); // qv, qc, qi, qr, qs, qg + int q_size = micro.Get_Qmoist_Size(); + qmoist[lev] = std::make_unique(ba, dm, q_size, ngrow_state); init_stuff(lev, ba, dm); @@ -302,7 +304,8 @@ ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapp //******************************************************************************************** // Microphysics // ******************************************************************************************* - qmoist[lev].define(ba, dm, 6, ngrow_state); // qv, qc, qi, qr, qs, qg + int q_size = micro.Get_Qmoist_Size(); + qmoist[lev]->define(ba, dm, q_size, ngrow_state); // qv, qc, qi, qr, qs, qg init_stuff(lev,ba,dm); @@ -553,6 +556,9 @@ ERF::ClearLevel (int lev) rW_new[lev].clear(); rW_old[lev].clear(); + // Clears the qmoist memory + qmoist[lev].reset(); + // Clears the integrator memory mri_integrator_mem[lev].reset(); physbcs[lev].reset(); diff --git a/Source/IO/Checkpoint.cpp b/Source/IO/Checkpoint.cpp index 13f4cd910..b1b587bc7 100644 --- a/Source/IO/Checkpoint.cpp +++ b/Source/IO/Checkpoint.cpp @@ -131,9 +131,9 @@ ERF::WriteCheckpointFile () const IntVect ng; if (solverChoice.moisture_type != MoistureType::None) { // We must read and write qmoist with ghost cells because we don't directly impose BCs on these vars - ng = qmoist[lev].nGrowVect(); - MultiFab moist_vars(grids[lev],dmap[lev],qmoist[lev].nComp(),ng); - MultiFab::Copy(moist_vars,qmoist[lev],0,0,qmoist[lev].nComp(),ng); + ng = qmoist[lev]->nGrowVect(); + MultiFab moist_vars(grids[lev],dmap[lev],qmoist[lev]->nComp(),ng); + MultiFab::Copy(moist_vars,*(qmoist[lev]),0,0,qmoist[lev]->nComp(),ng); VisMF::Write(moist_vars, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "MoistVars")); } @@ -350,12 +350,12 @@ ERF::ReadCheckpointFile () IntVect ng; if (solverChoice.moisture_type == MoistureType::None) { - qmoist[lev].setVal(0.0); + qmoist[lev]->setVal(0.0); } else { - ng = qmoist[lev].nGrowVect(); - MultiFab moist_vars(grids[lev],dmap[lev],qmoist[lev].nComp(),ng); + ng = qmoist[lev]->nGrowVect(); + MultiFab moist_vars(grids[lev],dmap[lev],qmoist[lev]->nComp(),ng); VisMF::Read(moist_vars, amrex::MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "MoistVars")); - MultiFab::Copy(qmoist[lev],moist_vars,0,0,qmoist[lev].nComp(),ng); + MultiFab::Copy(*(qmoist[lev]),moist_vars,0,0,qmoist[lev]->nComp(),ng); } // Note that we read the ghost cells of the base state (unlike above) diff --git a/Source/IO/ERF_ReadBndryPlanes.H b/Source/IO/ERF_ReadBndryPlanes.H index 4498e3e09..84a14139e 100644 --- a/Source/IO/ERF_ReadBndryPlanes.H +++ b/Source/IO/ERF_ReadBndryPlanes.H @@ -42,10 +42,6 @@ public: [[nodiscard]] int ingested_scalar() const {return is_scalar_read;} [[nodiscard]] int ingested_q1() const {return is_q1_read;} [[nodiscard]] int ingested_q2() const {return is_q2_read;} -#if defined(ERF_USE_WARM_NO_PRECIP) - int ingested_qv() const {return is_qv_read;} - int ingested_qc() const {return is_qc_read;} -#endif [[nodiscard]] int ingested_KE() const {return is_KE_read;} [[nodiscard]] int ingested_QKE() const {return is_QKE_read;} @@ -102,10 +98,6 @@ private: int is_scalar_read; int is_q1_read; int is_q2_read; -#if defined(ERF_USE_WARM_NO_PRECIP) - int is_qv_read; - int is_qc_read; -#endif int is_KE_read; int is_QKE_read; diff --git a/Source/IO/ERF_ReadBndryPlanes.cpp b/Source/IO/ERF_ReadBndryPlanes.cpp index 9d61f0654..1e7723ac9 100644 --- a/Source/IO/ERF_ReadBndryPlanes.cpp +++ b/Source/IO/ERF_ReadBndryPlanes.cpp @@ -157,10 +157,6 @@ ReadBndryPlanes::ReadBndryPlanes(const Geometry& geom, const Real& rdOcp_in) is_scalar_read = 0; is_q1_read = 0; is_q2_read = 0; -#if defined(ERF_USE_WARM_NO_PRECIP) - is_qv_read = 0; - is_qc_read = 0; -#endif is_KE_read = 0; is_QKE_read = 0; @@ -177,10 +173,6 @@ ReadBndryPlanes::ReadBndryPlanes(const Geometry& geom, const Real& rdOcp_in) if (m_var_names[i] == "scalar") is_scalar_read = 1; if (m_var_names[i] == "qt") is_q1_read = 1; if (m_var_names[i] == "qp") is_q2_read = 1; -#if defined(ERF_USE_WARM_NO_PRECIP) - if (m_var_names[i] == "qv") is_qv_read = 1; - if (m_var_names[i] == "qc") is_qc_read = 1; -#endif if (m_var_names[i] == "KE") is_KE_read = 1; if (m_var_names[i] == "QKE") is_QKE_read = 1; } @@ -412,10 +404,6 @@ void ReadBndryPlanes::read_file(const int idx, Vector& h_avg_u , Gpu::HostVector plot_var_names) const Box& bx = mfi.tilebox(); const Array4& derdat = mf[lev].array(mfi); const Array4& S_arr = vars_new[lev][Vars::cons].const_array(mfi); - const Array4& qv_arr = qmoist[0].const_array(mfi); + const Array4& qv_arr = qmoist[0]->const_array(mfi); ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { @@ -233,7 +233,7 @@ ERF::WritePlotFile (int which, Vector plot_var_names) const Array4& derdat = mf[lev].array(mfi); const Array4& p0_arr = p_hse.const_array(mfi); const Array4& S_arr = vars_new[lev][Vars::cons].const_array(mfi); - const Array4 & qv_arr = qmoist[0].const_array(mfi); + const Array4 & qv_arr = qmoist[0]->const_array(mfi); ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { @@ -599,12 +599,14 @@ ERF::WritePlotFile (int which, Vector plot_var_names) mf_comp ++; } - MultiFab qv_mf(qmoist[lev], make_alias, 0, 1); - MultiFab qc_mf(qmoist[lev], make_alias, 1, 1); - MultiFab qi_mf(qmoist[lev], make_alias, 2, 1); - MultiFab qr_mf(qmoist[lev], make_alias, 3, 1); - MultiFab qs_mf(qmoist[lev], make_alias, 4, 1); - MultiFab qg_mf(qmoist[lev], make_alias, 5, 1); + // TODO: Protect against accessing non-existent data + int q_size = micro.Get_Qmoist_Size(); + MultiFab qv_mf(*(qmoist[lev]), make_alias, 0, 1); + MultiFab qc_mf(*(qmoist[lev]), make_alias, 1, 1); + MultiFab qi_mf(*(qmoist[lev]), make_alias, 2, 1); + MultiFab qr_mf(*(qmoist[lev]), make_alias, 3, 1); + MultiFab qs_mf(*(qmoist[lev]), make_alias, 4, 1); + MultiFab qg_mf(*(qmoist[lev]), make_alias, 5, 1); if (containerHasElement(plot_var_names, "qt")) { diff --git a/Source/IO/Plotfile.cpp.convergence b/Source/IO/Plotfile.cpp.convergence index f5dc08436..b82e8e035 100644 --- a/Source/IO/Plotfile.cpp.convergence +++ b/Source/IO/Plotfile.cpp.convergence @@ -591,10 +591,6 @@ ERF::WritePlotFile (int which, Vector plot_var_names) mf_comp += 1; } } -#if defined(ERF_USE_WARM_NO_PRECIP) - calculate_derived("qv", derived::erf_derQv); - calculate_derived("qc", derived::erf_derQc); -#endif #ifdef ERF_COMPUTE_ERROR // Next, check for error in velocities and if desired, output them -- note we output none or all, not just some diff --git a/Source/Initialization/ERF_init_custom.cpp b/Source/Initialization/ERF_init_custom.cpp index fcc158133..2e2fb0586 100644 --- a/Source/Initialization/ERF_init_custom.cpp +++ b/Source/Initialization/ERF_init_custom.cpp @@ -25,7 +25,7 @@ void ERF::init_custom (int lev) { auto& lev_new = vars_new[lev]; - auto& qmoist_new = qmoist[lev]; + auto& qmoist_new = *(qmoist[lev]); MultiFab r_hse(base_state[lev], make_alias, 0, 1); // r_0 is first component MultiFab p_hse(base_state[lev], make_alias, 1, 1); // p_0 is second component @@ -42,7 +42,7 @@ ERF::init_custom (int lev) yvel_pert.setVal(0.); zvel_pert.setVal(0.); - MultiFab qmoist_pert(qmoist[lev].boxArray(), qmoist[lev].DistributionMap(), 3, qmoist[lev].nGrow()); + MultiFab qmoist_pert(qmoist[lev]->boxArray(), qmoist[lev]->DistributionMap(), 3, qmoist[lev]->nGrow()); qmoist_pert.setVal(0.); MultiFab qv_pert(qmoist_pert, amrex::make_alias, 0, 1); diff --git a/Source/Initialization/InputSoundingData.H b/Source/Initialization/InputSoundingData.H index 0496fb2cb..1f18033da 100644 --- a/Source/Initialization/InputSoundingData.H +++ b/Source/Initialization/InputSoundingData.H @@ -79,7 +79,7 @@ public: iss_z >> z >> theta >> qv >> U >> V; if (z == 0) { AMREX_ALWAYS_ASSERT(theta == theta_inp_sound_tmp[0]); - AMREX_ALWAYS_ASSERT(qv == qv_inp_sound_tmp[0]); + AMREX_ALWAYS_ASSERT(qv*0.001 == qv_inp_sound_tmp[0]); U_inp_sound_tmp[0] = U; V_inp_sound_tmp[0] = V; } else { diff --git a/Source/Microphysics/FastEddy/Diagnose_FE.cpp b/Source/Microphysics/FastEddy/Diagnose_FE.cpp index b3c1307c1..369a5c5d5 100644 --- a/Source/Microphysics/FastEddy/Diagnose_FE.cpp +++ b/Source/Microphysics/FastEddy/Diagnose_FE.cpp @@ -5,37 +5,4 @@ * from the existing Microphysics variables. */ void FastEddy::Diagnose () { - - auto qt = mic_fab_vars[MicVar_FE::qt]; - auto qp = mic_fab_vars[MicVar_FE::qp]; - auto qv = mic_fab_vars[MicVar_FE::qv]; - auto qn = mic_fab_vars[MicVar_FE::qn]; - auto qcl = mic_fab_vars[MicVar_FE::qcl]; - auto qci = mic_fab_vars[MicVar_FE::qci]; - auto qpl = mic_fab_vars[MicVar_FE::qpl]; - auto qpi = mic_fab_vars[MicVar_FE::qpi]; - auto tabs = mic_fab_vars[MicVar_FE::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); - - 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/FastEddy/FastEddy.H b/Source/Microphysics/FastEddy/FastEddy.H index df54a1736..e35354b62 100644 --- a/Source/Microphysics/FastEddy/FastEddy.H +++ b/Source/Microphysics/FastEddy/FastEddy.H @@ -20,22 +20,11 @@ namespace MicVar_FE { enum { // independent variables - qt = 0, - qp, + qv = 0, + qc, 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 }; } @@ -45,109 +34,116 @@ namespace MicVar_FE { // class FastEddy : public NullMoist { - using FabPtr = std::shared_ptr; - - public: - // constructor - FastEddy () {} - - // destructor - virtual ~FastEddy () = default; - - // cloud physics - void AdvanceFE (); - - // 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->AdvanceFE(); - this->Diagnose(); - } + using FabPtr = std::shared_ptr; + +public: + // constructor + FastEddy () {} + + // destructor + virtual ~FastEddy () = default; + + // cloud physics + void AdvanceFE (); + + // 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->AdvanceFE(); + this->Diagnose(); + } + + int + Qmoist_Size ( ) override { return FastEddy::m_qmoist_size; } 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; + // Number of qmoist variables + int m_qmoist_size = 2; + + // 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/FastEddy/FastEddy.cpp b/Source/Microphysics/FastEddy/FastEddy.cpp index e3575ed99..a07fdee14 100644 --- a/Source/Microphysics/FastEddy/FastEddy.cpp +++ b/Source/Microphysics/FastEddy/FastEddy.cpp @@ -11,204 +11,62 @@ using namespace amrex; /** * Compute Precipitation-related Microphysics quantities. */ -void FastEddy::AdvanceFE() { - - auto qt = mic_fab_vars[MicVar_FE::qt]; - auto qp = mic_fab_vars[MicVar_FE::qp]; - auto qn = mic_fab_vars[MicVar_FE::qn]; - auto tabs = mic_fab_vars[MicVar_FE::tabs]; - - auto qcl = mic_fab_vars[MicVar_FE::qcl]; - auto theta = mic_fab_vars[MicVar_FE::theta]; - auto qv = mic_fab_vars[MicVar_FE::qv]; - auto rho = mic_fab_vars[MicVar_FE::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(); - 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_FE::rho]->array(mfi); - auto qp_array = mic_fab_vars[MicVar_FE::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_FE::qn]->array(mfi); - auto qt_array = mic_fab_vars[MicVar_FE::qt]->array(mfi); - auto qp_array = mic_fab_vars[MicVar_FE::qp]->array(mfi); - auto qv_array = mic_fab_vars[MicVar_FE::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; - } - }); - } +void FastEddy::AdvanceFE () { + auto tabs = mic_fab_vars[MicVar_FE::tabs]; + // get the temperature, dentisy, theta, qt and qc from input for ( MFIter mfi(*tabs,TilingIfNotGPU()); mfi.isValid(); ++mfi) { - auto qp_array = mic_fab_vars[MicVar_FE::qp]->array(mfi); - auto rho_array = mic_fab_vars[MicVar_FE::rho]->array(mfi); + auto qv_array = mic_fab_vars[MicVar_FE::qv]->array(mfi); + auto qc_array = mic_fab_vars[MicVar_FE::qc]->array(mfi); + auto tabs_array = mic_fab_vars[MicVar_FE::tabs]->array(mfi); + auto theta_array = mic_fab_vars[MicVar_FE::theta]->array(mfi); + auto rho_array = mic_fab_vars[MicVar_FE::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; - }); - }*/ - + // Expose for GPU + Real d_fac_cond = m_fac_cond; + ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { - // 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_FE::tabs]->array(mfi); - auto qn_array = mic_fab_vars[MicVar_FE::qn]->array(mfi); - auto qt_array = mic_fab_vars[MicVar_FE::qt]->array(mfi); - auto qp_array = mic_fab_vars[MicVar_FE::qp]->array(mfi); + qc_array(i,j,k) = std::max(0.0, qc_array(i,j,k)); - auto theta_array = theta->array(mfi); - auto qv_array = qv->array(mfi); - auto rho_array = mic_fab_vars[MicVar_FE::rho]->array(mfi); + //------- Autoconversion/accretion + Real dq_clwater_to_vapor, dq_vapor_to_clwater, qsat; - const auto& box3d = mfi.tilebox(); + 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); - 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 qcc, autor, accrr, 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; + // 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_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){ + 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(qv_array(i,j,k) < qsat and qc_array(i,j,k) > 0.0){ + dq_clwater_to_vapor = std::min(qc_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 - } + qv_array(i,j,k) = qv_array(i,j,k) - dq_vapor_to_clwater + dq_clwater_to_vapor; + qc_array(i,j,k) = qc_array(i,j,k) + dq_vapor_to_clwater - dq_clwater_to_vapor; - // If there is cloud water present then do accretion and autoconversion 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); - if (qn_array(i,j,k) > 0.0) { - qcc = qn_array(i,j,k); + qv_array(i,j,k) = std::max(0.0, qv_array(i,j,k)); + qc_array(i,j,k) = std::max(0.0, qc_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/FastEddy/Init_FE.cpp b/Source/Microphysics/FastEddy/Init_FE.cpp index 033d56136..164016fc9 100644 --- a/Source/Microphysics/FastEddy/Init_FE.cpp +++ b/Source/Microphysics/FastEddy/Init_FE.cpp @@ -23,227 +23,55 @@ void FastEddy::Init (const MultiFab& cons_in, MultiFab& qmoist, const BoxArray& grids, const Geometry& geom, const Real& dt_advance) - { - m_geom = geom; - m_gtoe = grids; +{ + m_geom = geom; + m_gtoe = grids; - auto dz = m_geom.CellSize(2); - auto lowz = m_geom.ProbLo(2); + dt = dt_advance; - dt = dt_advance; + // initialize microphysics variables + for (auto ivar = 0; ivar < MicVar_FE::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.); + } - // initialize microphysics variables - for (auto ivar = 0; ivar < MicVar_FE::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 - // 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.); - //qmoist.setVal(0.); + for ( MFIter mfi(cons_in, TileNoZ()); mfi.isValid(); ++mfi) { - for ( MFIter mfi(cons_in, TileNoZ()); mfi.isValid(); ++mfi) { + const auto& box3d = mfi.tilebox(); - const auto& box3d = mfi.tilebox(); + const auto& lo = amrex::lbound(box3d); + const auto& hi = amrex::ubound(box3d); - const auto& lo = amrex::lbound(box3d); - const auto& hi = amrex::ubound(box3d); + nlev = box3d.length(2); + zlo = lo.z; + zhi = hi.z; + } - nlev = box3d.length(2); - zlo = lo.z; - zhi = hi.z; + // 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); - // 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}); + auto qv_array = mic_fab_vars[MicVar_FE::qv]->array(mfi); + auto qc_array = mic_fab_vars[MicVar_FE::qc]->array(mfi); + auto rho_array = mic_fab_vars[MicVar_FE::rho]->array(mfi); + auto theta_array = mic_fab_vars[MicVar_FE::theta]->array(mfi); + auto temp_array = mic_fab_vars[MicVar_FE::tabs]->array(mfi); - // data (input) - rho1d.resize({zlo}, {zhi}); - pres1d.resize({zlo}, {zhi}); - tabs1d.resize({zlo}, {zhi}); - gamaz.resize({zlo}, {zhi}); - zmid.resize({zlo}, {zhi}); + const auto& box3d = mfi.tilebox(); - // 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_FE::qt]->array(mfi); - auto qp_array = mic_fab_vars[MicVar_FE::qp]->array(mfi); - auto qn_array = mic_fab_vars[MicVar_FE::qn]->array(mfi); - auto qv_array = mic_fab_vars[MicVar_FE::qv]->array(mfi); - auto rho_array = mic_fab_vars[MicVar_FE::rho]->array(mfi); - auto theta_array = mic_fab_vars[MicVar_FE::theta]->array(mfi); - auto temp_array = mic_fab_vars[MicVar_FE::tabs]->array(mfi); - auto pres_array = mic_fab_vars[MicVar_FE::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,RhoQ1_comp)/states_array(i,j,k,Rho_comp); - qp_array(i,j,k) = states_array(i,j,k,RhoQ2_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); - }); + // 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); + qv_array(i,j,k) = states_array(i,j,k,RhoQ1_comp)/states_array(i,j,k,Rho_comp); + qc_array(i,j,k) = states_array(i,j,k,RhoQ2_comp)/states_array(i,j,k,Rho_comp); + 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)); + }); + } } diff --git a/Source/Microphysics/FastEddy/Update_FE.cpp b/Source/Microphysics/FastEddy/Update_FE.cpp index a157279f7..eaedc3697 100644 --- a/Source/Microphysics/FastEddy/Update_FE.cpp +++ b/Source/Microphysics/FastEddy/Update_FE.cpp @@ -13,46 +13,29 @@ void FastEddy::Update (amrex::MultiFab& cons, amrex::MultiFab& qmoist) { - // copy multifab data to qc, qv, and qi - amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar_FE::qv], 0, 0, 1, mic_fab_vars[MicVar_FE::qv]->nGrowVect()); // vapor - amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar_FE::qcl], 0, 1, 1, mic_fab_vars[MicVar_FE::qcl]->nGrowVect()); // cloud water - amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar_FE::qci], 0, 2, 1, mic_fab_vars[MicVar_FE::qci]->nGrowVect()); // cloud ice - amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar_FE::qpl], 0, 3, 1, mic_fab_vars[MicVar_FE::qpl]->nGrowVect()); // rain - amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar_FE::qpi], 0, 4, 1, mic_fab_vars[MicVar_FE::qpi]->nGrowVect()); // snow - - // Don't need to copy this since it is filled below - // amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar_FE::qpi], 0, 5, 1, mic_fab_vars[MicVar_FE::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_FE::rho]->array(mfi); - auto theta_arr = mic_fab_vars[MicVar_FE::theta]->array(mfi); - auto qt_arr = mic_fab_vars[MicVar_FE::qt]->array(mfi); - auto qp_arr = mic_fab_vars[MicVar_FE::qp]->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,RhoQ1_comp) = rho_arr(i,j,k)*qt_arr(i,j,k); - states_arr(i,j,k,RhoQ2_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()); + + // 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_FE::rho]->array(mfi); + auto theta_arr = mic_fab_vars[MicVar_FE::theta]->array(mfi); + auto qv_arr = mic_fab_vars[MicVar_FE::qv]->array(mfi); + auto qc_arr = mic_fab_vars[MicVar_FE::qc]->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,RhoTheta_comp) = rho_arr(i,j,k)*theta_arr(i,j,k); + states_arr(i,j,k,RhoQ1_comp) = rho_arr(i,j,k)*qv_arr(i,j,k); + states_arr(i,j,k,RhoQ2_comp) = rho_arr(i,j,k)*qc_arr(i,j,k); + }); + } + + // Fill interior ghost cells and periodic boundaries + cons.FillBoundary(m_geom.periodicity()); + qmoist.FillBoundary(m_geom.periodicity()); } diff --git a/Source/Microphysics/Kessler/Diagnose_Kessler.cpp b/Source/Microphysics/Kessler/Diagnose_Kessler.cpp index 43210f501..c231053dd 100644 --- a/Source/Microphysics/Kessler/Diagnose_Kessler.cpp +++ b/Source/Microphysics/Kessler/Diagnose_Kessler.cpp @@ -6,36 +6,37 @@ */ 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]; + 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); + 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); - const auto& box3d = mfi.tilebox(); + 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); - }); - } + 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 index f8a3d835f..fcf469a2a 100644 --- a/Source/Microphysics/Kessler/Init_Kessler.cpp +++ b/Source/Microphysics/Kessler/Init_Kessler.cpp @@ -19,231 +19,236 @@ using namespace amrex; * @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) +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,RhoQ1_comp)/states_array(i,j,k,Rho_comp); - qp_array(i,j,k) = states_array(i,j,k,RhoQ2_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.; + 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,RhoQ1_comp)/states_array(i,j,k,Rho_comp); + qp_array(i,j,k) = states_array(i,j,k,RhoQ2_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); }); - } - // 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(); + // 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; - }); + 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); - }); + 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 index 7ef58dd34..f771575e0 100644 --- a/Source/Microphysics/Kessler/Kessler.H +++ b/Source/Microphysics/Kessler/Kessler.H @@ -51,109 +51,115 @@ namespace MicVar_Kess { // 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(); - } + 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(); + } + + int + Qmoist_Size ( ) override { return Kessler::m_qmoist_size; } 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; + // Number of qmoist variables + int m_qmoist_size = 6; + + // 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 index 64fe35639..580799d6c 100644 --- a/Source/Microphysics/Kessler/Kessler.cpp +++ b/Source/Microphysics/Kessler/Kessler.cpp @@ -11,36 +11,36 @@ using namespace amrex; /** * Compute Precipitation-related Microphysics quantities. */ -void Kessler::AdvanceKessler() { +void Kessler::AdvanceKessler () { - 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 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 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); + 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(); - fz.define(convert(ba, IntVect(0,0,1)), dm, 1, 0); // No ghost cells + MultiFab fz; + auto ba = tabs->boxArray(); + auto dm = tabs->DistributionMap(); + 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 { + 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) { @@ -61,12 +61,12 @@ void Kessler::AdvanceKessler() { fz_array(i,j,k) = rho_avg*V_terminal*qp_avg; /*if(k==0){ - fz_array(i,j,k) = 0; - }*/ + fz_array(i,j,k) = 0; + }*/ }); } - Real dtn = dt; + Real dtn = dt; /*for ( MFIter mfi(*tabs,TilingIfNotGPU()); mfi.isValid(); ++mfi) { @@ -107,56 +107,57 @@ void Kessler::AdvanceKessler() { - // 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); + // 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 theta_array = theta->array(mfi); - auto qv_array = qv->array(mfi); - auto rho_array = mic_fab_vars[MicVar_Kess::rho]->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(); + const auto& box3d = mfi.tilebox(); - auto fz_array = fz.array(mfi); + auto fz_array = fz.array(mfi); - // Expose for GPU - Real d_fac_cond = m_fac_cond; + // Expose for GPU + Real d_fac_cond = m_fac_cond; - ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + 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)); + 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; - } + 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 qcc, autor, accrr, dq_clwater_to_rain, dq_rain_to_vapor, dq_clwater_to_vapor, dq_vapor_to_clwater, qsat; + //------- Autoconversion/accretion + Real qcc, autor, accrr, 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); + 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 + // 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; + 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)); + //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){ + 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)); } @@ -167,52 +168,52 @@ void Kessler::AdvanceKessler() { 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) { + 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}); + 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 - } + // 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); + if (qn_array(i,j,k) > 0.0) { + qcc = qn_array(i,j,k); - autor = 0.0; - if (qcc > qcw0) { - autor = 0.001; - } + 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)); + 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)); - } + // 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)); + } - if(std::fabs(fz_array(i,j,k+1)) < 1e-14) fz_array(i,j,k+1) = 0.0; - if(std::fabs(fz_array(i,j,k)) < 1e-14) fz_array(i,j,k) = 0.0; - Real dq_sed = 1.0/rho_array(i,j,k)*(fz_array(i,j,k+1) - fz_array(i,j,k))/dz*dtn; - if(std::fabs(dq_sed) < 1e-14)dq_sed = 0.0; - //dq_sed = 0.0; + if(std::fabs(fz_array(i,j,k+1)) < 1e-14) fz_array(i,j,k+1) = 0.0; + if(std::fabs(fz_array(i,j,k)) < 1e-14) fz_array(i,j,k) = 0.0; + Real dq_sed = 1.0/rho_array(i,j,k)*(fz_array(i,j,k+1) - fz_array(i,j,k))/dz*dtn; + if(std::fabs(dq_sed) < 1e-14)dq_sed = 0.0; + //dq_sed = 0.0; - 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; + 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); + 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)); + 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/Update_Kessler.cpp b/Source/Microphysics/Kessler/Update_Kessler.cpp index 304c812a3..470b839e3 100644 --- a/Source/Microphysics/Kessler/Update_Kessler.cpp +++ b/Source/Microphysics/Kessler/Update_Kessler.cpp @@ -11,48 +11,49 @@ * @param[out] qmoist: qv, qc, qi, qr, qs, qg */ void Kessler::Update (amrex::MultiFab& cons, - amrex::MultiFab& qmoist) + 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 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,RhoQ1_comp) = rho_arr(i,j,k)*qt_arr(i,j,k); - states_arr(i,j,k,RhoQ2_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()); + // 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 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,RhoQ1_comp) = rho_arr(i,j,k)*qt_arr(i,j,k); + states_arr(i,j,k,RhoQ2_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 c384ad8a0..770b6cc80 100644 --- a/Source/Microphysics/Microphysics.H +++ b/Source/Microphysics/Microphysics.H @@ -56,6 +56,9 @@ public: m_moist_model->Update(cons_in, qmoist); } + int + Get_Qmoist_Size ( ) { return m_moist_model->Qmoist_Size(); } + private: std::unique_ptr m_moist_model; diff --git a/Source/Microphysics/Null/NullMoist.H b/Source/Microphysics/Null/NullMoist.H index 1484058fa..150f62118 100644 --- a/Source/Microphysics/Null/NullMoist.H +++ b/Source/Microphysics/Null/NullMoist.H @@ -12,7 +12,8 @@ class NullMoist { virtual ~NullMoist () = default; - virtual void + virtual + void define (SolverChoice& /*sc*/) { } virtual @@ -22,17 +23,24 @@ class NullMoist { const amrex::Geometry& /*geom*/, const amrex::Real& /*dt_advance*/) { } - virtual void + virtual + void Advance ( ) { } - virtual void + virtual + void Diagnose ( ) { } - virtual void + virtual + void Update (amrex::MultiFab& /*cons_in*/, amrex::MultiFab& /*qmoist*/) { } - private: + virtual + int + Qmoist_Size ( ) { return NullMoist::m_qmoist_size; } + private: + int m_qmoist_size = 1; }; #endif diff --git a/Source/Microphysics/SAM/SAM.H b/Source/Microphysics/SAM/SAM.H index e9d263fe2..ad6434ecf 100644 --- a/Source/Microphysics/SAM/SAM.H +++ b/Source/Microphysics/SAM/SAM.H @@ -51,127 +51,133 @@ namespace MicVar { // class SAM : public NullMoist { - using FabPtr = std::shared_ptr; - - public: - // constructor - SAM () {} - - // destructor - virtual ~SAM () = default; - - // cloud physics - void Cloud (); - - // ice physics - void IceFall (); - - // precip - void Precip (); - - // precip fall - void PrecipFall (int hydro_type); - - // micro interface for precip fall - void MicroPrecipFall (); - - // process microphysics - void Proc (); - - // 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->Cloud(); - this->Diagnose(); - this->IceFall(); - this->Precip(); - this->MicroPrecipFall(); - } + using FabPtr = std::shared_ptr; + +public: + // constructor + SAM () {} + + // destructor + virtual ~SAM () = default; + + // cloud physics + void Cloud (); + + // ice physics + void IceFall (); + + // precip + void Precip (); + + // precip fall + void PrecipFall (int hydro_type); + + // micro interface for precip fall + void MicroPrecipFall (); + + // process microphysics + void Proc (); + + // 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->Cloud(); + this->Diagnose(); + this->IceFall(); + this->Precip(); + this->MicroPrecipFall(); + } + + int + Qmoist_Size ( ) override { return SAM::m_qmoist_size; } 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; + // Number of qmoist variables + int m_qmoist_size = 6; + + // 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/TimeIntegration/ERF_Advance.cpp b/Source/TimeIntegration/ERF_Advance.cpp index d3759c6b4..362464d16 100644 --- a/Source/TimeIntegration/ERF_Advance.cpp +++ b/Source/TimeIntegration/ERF_Advance.cpp @@ -54,8 +54,9 @@ ERF::Advance (int lev, Real time, Real dt_lev, int /*iteration*/, int /*ncycle*/ FillPatch(lev, time, {&vars_old[lev][Vars::cons], &vars_old[lev][Vars::xvel], &vars_old[lev][Vars::yvel], &vars_old[lev][Vars::zvel]}); + if (solverChoice.moisture_type != MoistureType::None) { - FillPatchMoistVars(lev, qmoist[lev]); + FillPatchMoistVars(lev, *qmoist[lev]); } MultiFab* S_crse; @@ -104,11 +105,6 @@ ERF::Advance (int lev, Real time, Real dt_lev, int /*iteration*/, int /*ncycle*/ // Place-holder for source array -- for now just set to 0 MultiFab source(ba,dm,nvars,1); source.setVal(0.0); -#if defined(ERF_USE_WARM_NO_PRECIP) - Real tau_cond = solverChoice.tau_cond; - Real c_p = solverChoice.c_p; - condensation_source(source, S_new, tau_cond, c_p); -#endif // We don't need to call FillPatch on cons_mf because we have fillpatch'ed S_old above MultiFab cons_mf(ba,dm,nvars,S_old.nGrowVect()); @@ -128,7 +124,7 @@ ERF::Advance (int lev, Real time, Real dt_lev, int /*iteration*/, int /*ncycle*/ source, buoyancy, Geom(lev), dt_lev, time, &ifr); - // Update the microphysics (moisture) + // Update the microphysics (moisture) advance_microphysics(lev, S_new, dt_lev); #if defined(ERF_USE_RRTMGP) diff --git a/Source/TimeIntegration/ERF_advance_dycore.cpp b/Source/TimeIntegration/ERF_advance_dycore.cpp index 7de02c47a..3ed147cc6 100644 --- a/Source/TimeIntegration/ERF_advance_dycore.cpp +++ b/Source/TimeIntegration/ERF_advance_dycore.cpp @@ -250,37 +250,44 @@ void ERF::advance_dycore(int level, cons_to_prim(state_old[IntVar::cons], state_old[IntVar::cons].nGrow()); } // profile -#if 0 && defined(ERF_USE_MOISTURE) - MultiFab qvapor (qmoist[level], make_alias, 0, 1); - MultiFab qcloud (qmoist[level], make_alias, 1, 1); - MultiFab qice (qmoist[level], make_alias, 2, 1); + int q_size = micro.Get_Qmoist_Size(); + int ncell = fine_geom.Domain().length(2); - PlaneAverage qv_ave(&qvapor, geom[level], solverChoice.ave_plane); - PlaneAverage qc_ave(&qcloud, geom[level], solverChoice.ave_plane); - PlaneAverage qi_ave(&qice , geom[level], solverChoice.ave_plane); + Gpu::DeviceVector qv_d(ncell), qc_d(ncell), qi_d(ncell); - // Compute plane averages - qv_ave.compute_averages(ZDir(), qv_ave.field()); - qc_ave.compute_averages(ZDir(), qc_ave.field()); - qi_ave.compute_averages(ZDir(), qi_ave.field()); + if (q_size >= 1) { + MultiFab qvapor (*(qmoist[level]), make_alias, 0, 1); + PlaneAverage qv_ave(&qvapor, geom[level], solverChoice.ave_plane); + qv_ave.compute_averages(ZDir(), qv_ave.field()); - // get plane averaged data - int ncell = qv_ave.ncell_line(); + Gpu::HostVector qv_h(ncell); - Gpu::HostVector qv_h(ncell), qi_h(ncell), qc_h(ncell); - Gpu::DeviceVector qv_d(ncell), qc_d(ncell), qi_d(ncell); + qv_ave.line_average(0, qv_h); + Gpu::copyAsync(Gpu::hostToDevice, qv_h.begin(), qv_h.end(), qv_d.begin()); + } - // Fill the vectors with the line averages computed above - qv_ave.line_average(0, qv_h); - qi_ave.line_average(0, qi_h); - qc_ave.line_average(0, qc_h); + if (q_size >= 2) { + MultiFab qcloud (*(qmoist[level]), make_alias, 1, 1); + PlaneAverage qc_ave(&qcloud, geom[level], solverChoice.ave_plane); + qc_ave.compute_averages(ZDir(), qc_ave.field()); - // Copy data to device - Gpu::copyAsync(Gpu::hostToDevice, qv_h.begin(), qv_h.end(), qv_d.begin()); - Gpu::copyAsync(Gpu::hostToDevice, qi_h.begin(), qi_h.end(), qi_d.begin()); - Gpu::copyAsync(Gpu::hostToDevice, qc_h.begin(), qc_h.end(), qc_d.begin()); - Gpu::streamSynchronize(); -#endif + Gpu::HostVector qc_h(ncell); + + qc_ave.line_average(0, qc_h); + Gpu::copyAsync(Gpu::hostToDevice, qc_h.begin(), qc_h.end(), qc_d.begin()); + } + + if (q_size >= 3) { + MultiFab qice (*(qmoist[level]), make_alias, 2, 1); + PlaneAverage qi_ave(&qice , geom[level], solverChoice.ave_plane); + qi_ave.compute_averages(ZDir(), qi_ave.field()); + + Gpu::HostVector qi_h(ncell); + + qi_ave.line_average(0, qi_h); + Gpu::copyAsync(Gpu::hostToDevice, qi_h.begin(), qi_h.end(), qi_d.begin()); + Gpu::streamSynchronize(); + } #include "TI_no_substep_fun.H" #include "TI_slow_rhs_fun.H" diff --git a/Source/TimeIntegration/ERF_advance_microphysics.cpp b/Source/TimeIntegration/ERF_advance_microphysics.cpp index a9abaa4ed..93373efcf 100644 --- a/Source/TimeIntegration/ERF_advance_microphysics.cpp +++ b/Source/TimeIntegration/ERF_advance_microphysics.cpp @@ -6,12 +6,12 @@ void ERF::advance_microphysics (int lev, MultiFab& cons, const Real& dt_advance) { - if (moisture_type != MoistureType::None) { - micro.Init(cons, qmoist[lev], + if (solverChoice.moisture_type != MoistureType::None) { + micro.Init(cons, *(qmoist[lev]), grids[lev], Geom(lev), dt_advance); micro.Advance(); - micro.Update(cons, qmoist[lev]); + micro.Update(cons, *(qmoist[lev])); } } diff --git a/Source/TimeIntegration/ERF_fast_rhs_MT.cpp b/Source/TimeIntegration/ERF_fast_rhs_MT.cpp index 03fc5f3d3..c1c689652 100644 --- a/Source/TimeIntegration/ERF_fast_rhs_MT.cpp +++ b/Source/TimeIntegration/ERF_fast_rhs_MT.cpp @@ -79,12 +79,11 @@ void erf_fast_rhs_MT (int step, int nrk, std::unique_ptr& mapfac_v, YAFluxRegister* fr_as_crse, YAFluxRegister* fr_as_fine, + bool l_use_moisture, bool l_reflux) { BL_PROFILE_REGION("erf_fast_rhs_MT()"); - const bool use_moisture = (solverChoice.moisture_type != Moisture::None); - Real beta_1 = 0.5 * (1.0 - beta_s); // multiplies explicit terms Real beta_2 = 0.5 * (1.0 + beta_s); // multiplies implicit terms @@ -260,7 +259,7 @@ void erf_fast_rhs_MT (int step, int nrk, Real gpx = h_zeta_old * gp_xi - h_xi_old * gp_zeta_on_iface; gpx *= mf_u(i,j,0); - if (use_moisture) { + if (l_use_moisture) { Real q = 0.5 * ( prim(i,j,k,PrimQ1_comp) + prim(i-1,j,k,PrimQ1_comp) +prim(i,j,k,PrimQ2_comp) + prim(i-1,j,k,PrimQ2_comp) ); gpx /= (1.0 + q); @@ -287,7 +286,7 @@ void erf_fast_rhs_MT (int step, int nrk, Real gpy = h_zeta_old * gp_eta - h_eta_old * gp_zeta_on_jface; gpy *= mf_v(i,j,0); - if (use_moisture) { + if (l_use_moisture) { Real q = 0.5 * ( prim(i,j,k,PrimQ1_comp) + prim(i,j-1,k,PrimQ1_comp) +prim(i,j,k,PrimQ2_comp) + prim(i,j-1,k,PrimQ2_comp) ); gpy /= (1.0 + q); @@ -397,7 +396,7 @@ void erf_fast_rhs_MT (int step, int nrk, Real coeff_P = coeffP_a(i,j,k); Real coeff_Q = coeffQ_a(i,j,k); - if (use_moisture) { + if (l_use_moisture) { Real q = 0.5 * ( prim(i,j,k,PrimQ1_comp) + prim(i,j,k-1,PrimQ1_comp) +prim(i,j,k,PrimQ2_comp) + prim(i,j,k-1,PrimQ2_comp) ); coeff_P /= (1.0 + q); diff --git a/Source/TimeIntegration/ERF_fast_rhs_N.cpp b/Source/TimeIntegration/ERF_fast_rhs_N.cpp index 77c40685c..45f47fc80 100644 --- a/Source/TimeIntegration/ERF_fast_rhs_N.cpp +++ b/Source/TimeIntegration/ERF_fast_rhs_N.cpp @@ -58,12 +58,11 @@ void erf_fast_rhs_N (int step, int nrk, std::unique_ptr& mapfac_v, YAFluxRegister* fr_as_crse, YAFluxRegister* fr_as_fine, + bool l_use_moisture, bool l_reflux) { BL_PROFILE_REGION("erf_fast_rhs_N()"); - const bool use_moisture = (solverChoice.moisture_type != Moisture::None); - Real beta_1 = 0.5 * (1.0 - beta_s); // multiplies explicit terms Real beta_2 = 0.5 * (1.0 + beta_s); // multiplies implicit terms @@ -222,7 +221,7 @@ void erf_fast_rhs_N (int step, int nrk, Real gpx = (theta_extrap(i,j,k) - theta_extrap(i-1,j,k))*dxi; gpx *= mf_u(i,j,0); - if (use_moisture) { + if (l_use_moisture) { Real q = 0.5 * ( prim(i,j,k,PrimQ1_comp) + prim(i-1,j,k,PrimQ1_comp) +prim(i,j,k,PrimQ2_comp) + prim(i-1,j,k,PrimQ2_comp) ); gpx /= (1.0 + q); @@ -245,7 +244,7 @@ void erf_fast_rhs_N (int step, int nrk, Real gpy = (theta_extrap(i,j,k) - theta_extrap(i,j-1,k))*dyi; gpy *= mf_v(i,j,0); - if (use_moisture) { + if (l_use_moisture) { Real q = 0.5 * ( prim(i,j,k,PrimQ1_comp) + prim(i,j-1,k,PrimQ1_comp) +prim(i,j,k,PrimQ2_comp) + prim(i,j-1,k,PrimQ2_comp) ); gpy /= (1.0 + q); @@ -376,7 +375,7 @@ void erf_fast_rhs_N (int step, int nrk, Real coeff_P = coeffP_a(i,j,k); Real coeff_Q = coeffQ_a(i,j,k); - if (use_moisture) { + if (l_use_moisture) { Real q = 0.5 * ( prim(i,j,k,PrimQ1_comp) + prim(i,j,k-1,PrimQ1_comp) +prim(i,j,k,PrimQ2_comp) + prim(i,j,k-1,PrimQ2_comp) ); coeff_P /= (1.0 + q); diff --git a/Source/TimeIntegration/ERF_fast_rhs_T.cpp b/Source/TimeIntegration/ERF_fast_rhs_T.cpp index 09282dbfe..4efcea56d 100644 --- a/Source/TimeIntegration/ERF_fast_rhs_T.cpp +++ b/Source/TimeIntegration/ERF_fast_rhs_T.cpp @@ -65,12 +65,11 @@ void erf_fast_rhs_T (int step, int nrk, std::unique_ptr& mapfac_v, YAFluxRegister* fr_as_crse, YAFluxRegister* fr_as_fine, + bool l_use_moisture, bool l_reflux) { BL_PROFILE_REGION("erf_fast_rhs_T()"); - const bool use_moisture = (solverChoice.moisture_type != MoistureType::None); - Real beta_1 = 0.5 * (1.0 - beta_s); // multiplies explicit terms Real beta_2 = 0.5 * (1.0 + beta_s); // multiplies implicit terms @@ -257,7 +256,7 @@ void erf_fast_rhs_T (int step, int nrk, Real gpx = gp_xi - (met_h_xi / met_h_zeta) * gp_zeta_on_iface; gpx *= mf_u(i,j,0); - if (use_moisture) { + if (l_use_moisture) { Real q = 0.5 * ( prim(i,j,k,PrimQ1_comp) + prim(i-1,j,k,PrimQ1_comp) +prim(i,j,k,PrimQ2_comp) + prim(i-1,j,k,PrimQ2_comp) ); gpx /= (1.0 + q); @@ -287,7 +286,7 @@ void erf_fast_rhs_T (int step, int nrk, Real gpy = gp_eta - (met_h_eta / met_h_zeta) * gp_zeta_on_jface; gpy *= mf_v(i,j,0); - if (use_moisture) { + if (l_use_moisture) { Real q = 0.5 * ( prim(i,j,k,PrimQ1_comp) + prim(i,j-1,k,PrimQ1_comp) +prim(i,j,k,PrimQ2_comp) + prim(i,j-1,k,PrimQ2_comp) ); gpy /= (1.0 + q); @@ -462,7 +461,7 @@ void erf_fast_rhs_T (int step, int nrk, Real coeff_P = coeffP_a(i,j,k); Real coeff_Q = coeffQ_a(i,j,k); - if (use_moisture) { + if (l_use_moisture) { Real q = 0.5 * ( prim(i,j,k,PrimQ1_comp) + prim(i,j,k-1,PrimQ1_comp) +prim(i,j,k,PrimQ2_comp) + prim(i,j,k-1,PrimQ2_comp) ); coeff_P /= (1.0 + q); diff --git a/Source/TimeIntegration/ERF_make_buoyancy.cpp b/Source/TimeIntegration/ERF_make_buoyancy.cpp index 4851f1b19..22c540c4c 100644 --- a/Source/TimeIntegration/ERF_make_buoyancy.cpp +++ b/Source/TimeIntegration/ERF_make_buoyancy.cpp @@ -24,9 +24,7 @@ using namespace amrex; * @param[in] S_data current solution * @param[in] S_prim primitive variables (i.e. conserved variables divided by density) * @param[out] buoyancy the buoyancy term computed here - * @param[in] qvapor water vapor - * @param[in] qcloud cloud water - * @param[in] qice cloud ice + * @param[in] qmoist moisture variables (in order: qv, qc, qi, ...) * @param[in] qv_d lateral average of cloud vapor * @param[in] qc_d lateral average of cloud vapor * @param[in] qd_d lateral average of cloud vapor @@ -38,12 +36,10 @@ using namespace amrex; void make_buoyancy (Vector& S_data, const MultiFab& S_prim, MultiFab& buoyancy, -#if defined(ERF_USE_MOISTURE) - const MultiFab& qmoist, + MultiFab* qmoist, Gpu::DeviceVector qv_d, Gpu::DeviceVector qc_d, Gpu::DeviceVector qi_d, -#endif const amrex::Geometry geom, const SolverChoice& solverChoice, const MultiFab* r0) @@ -56,95 +52,88 @@ void make_buoyancy (Vector& S_data, const int klo = 0; const int khi = geom.Domain().bigEnd()[2] + 1; -#if defined(ERF_USE_MOISTURE) - MultiFab qvapor(qmoist, make_alias, 0, 1); - MultiFab qcloud(qmoist, make_alias, 1, 1); - MultiFab qice (qmoist, make_alias, 2, 1); -#endif - // ****************************************************************************************** // Dry versions of buoyancy expressions (type 1 and type 2/3 -- types 2 and 3 are equivalent) // ****************************************************************************************** -#if !defined(ERF_USE_MOISTURE) && !defined(ERF_USE_WARM_NO_PRECIP) - - if (solverChoice.buoyancy_type == 1) { - for ( MFIter mfi(buoyancy,TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - Box tbz = mfi.tilebox(); + if (solverChoice.moisture_type == MoistureType::None) { + if (solverChoice.buoyancy_type == 1) { + for ( MFIter mfi(buoyancy,TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + Box tbz = mfi.tilebox(); - // We don't compute a source term for z-momentum on the bottom or top domain boundary - if (tbz.smallEnd(2) == klo) tbz.growLo(2,-1); - if (tbz.bigEnd(2) == khi) tbz.growHi(2,-1); + // We don't compute a source term for z-momentum on the bottom or top domain boundary + if (tbz.smallEnd(2) == klo) tbz.growLo(2,-1); + if (tbz.bigEnd(2) == khi) tbz.growHi(2,-1); - const Array4 & cell_data = S_data[IntVar::cons].array(mfi); - const Array4< Real> & buoyancy_fab = buoyancy.array(mfi); + const Array4 & cell_data = S_data[IntVar::cons].array(mfi); + const Array4< Real> & buoyancy_fab = buoyancy.array(mfi); - // Base state density - const Array4& r0_arr = r0->const_array(mfi); + // Base state density + const Array4& r0_arr = r0->const_array(mfi); - amrex::ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - buoyancy_fab(i, j, k) = grav_gpu[2] * 0.5 * ( (cell_data(i,j,k ) - r0_arr(i,j,k )) - +(cell_data(i,j,k-1) - r0_arr(i,j,k-1)) ); - }); - } // mfi + amrex::ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + buoyancy_fab(i, j, k) = grav_gpu[2] * 0.5 * ( (cell_data(i,j,k ) - r0_arr(i,j,k )) + +(cell_data(i,j,k-1) - r0_arr(i,j,k-1)) ); + }); + } // mfi - } else if (solverChoice.buoyancy_type == 2 || solverChoice.buoyancy_type == 3) { + } else if (solverChoice.buoyancy_type == 2 || solverChoice.buoyancy_type == 3) { - PlaneAverage state_ave(&(S_data[IntVar::cons]), geom, solverChoice.ave_plane); - PlaneAverage prim_ave(&S_prim, geom, solverChoice.ave_plane); + PlaneAverage state_ave(&(S_data[IntVar::cons]), geom, solverChoice.ave_plane); + PlaneAverage prim_ave(&S_prim, geom, solverChoice.ave_plane); - int ncell = state_ave.ncell_line(); + int ncell = state_ave.ncell_line(); - state_ave.compute_averages(ZDir(), state_ave.field()); - prim_ave.compute_averages(ZDir(), prim_ave.field()); + state_ave.compute_averages(ZDir(), state_ave.field()); + prim_ave.compute_averages(ZDir(), prim_ave.field()); - Gpu::HostVector rho_h(ncell), theta_h(ncell); - state_ave.line_average(Rho_comp, rho_h); - prim_ave.line_average(PrimTheta_comp, theta_h); + Gpu::HostVector rho_h(ncell), theta_h(ncell); + state_ave.line_average(Rho_comp, rho_h); + prim_ave.line_average(PrimTheta_comp, theta_h); - Gpu::DeviceVector rho_d(ncell); - Gpu::DeviceVector theta_d(ncell); + Gpu::DeviceVector rho_d(ncell); + Gpu::DeviceVector theta_d(ncell); - Gpu::copyAsync(Gpu::hostToDevice, rho_h.begin(), rho_h.end(), rho_d.begin()); - Gpu::copyAsync(Gpu::hostToDevice, theta_h.begin(), theta_h.end(), theta_d.begin()); + Gpu::copyAsync(Gpu::hostToDevice, rho_h.begin(), rho_h.end(), rho_d.begin()); + Gpu::copyAsync(Gpu::hostToDevice, theta_h.begin(), theta_h.end(), theta_d.begin()); - Real* rho_d_ptr = rho_d.data(); - Real* theta_d_ptr = theta_d.data(); + Real* rho_d_ptr = rho_d.data(); + Real* theta_d_ptr = theta_d.data(); #ifdef _OPENMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) #endif - for ( MFIter mfi(buoyancy,TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - Box tbz = mfi.tilebox(); + for ( MFIter mfi(buoyancy,TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + Box tbz = mfi.tilebox(); - // We don't compute a source term for z-momentum on the bottom or top boundary - if (tbz.smallEnd(2) == klo) tbz.growLo(2,-1); - if (tbz.bigEnd(2) == khi) tbz.growHi(2,-1); + // We don't compute a source term for z-momentum on the bottom or top boundary + if (tbz.smallEnd(2) == klo) tbz.growLo(2,-1); + if (tbz.bigEnd(2) == khi) tbz.growHi(2,-1); - const Array4 & cell_data = S_data[IntVar::cons].array(mfi); - const Array4< Real> & buoyancy_fab = buoyancy.array(mfi); + const Array4 & cell_data = S_data[IntVar::cons].array(mfi); + const Array4< Real> & buoyancy_fab = buoyancy.array(mfi); - amrex::ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - Real tempp1d = getTgivenRandRTh(rho_d_ptr[k ], rho_d_ptr[k ]*theta_d_ptr[k ]); - Real tempm1d = getTgivenRandRTh(rho_d_ptr[k-1], rho_d_ptr[k-1]*theta_d_ptr[k-1]); + amrex::ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + Real tempp1d = getTgivenRandRTh(rho_d_ptr[k ], rho_d_ptr[k ]*theta_d_ptr[k ]); + Real tempm1d = getTgivenRandRTh(rho_d_ptr[k-1], rho_d_ptr[k-1]*theta_d_ptr[k-1]); - 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 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 = (tempp3d-tempp1d)/tempp1d; - Real qminus = (tempm3d-tempm1d)/tempm1d; + Real qplus = (tempp3d-tempp1d)/tempp1d; + Real qminus = (tempm3d-tempm1d)/tempm1d; - Real qavg = Real(0.5) * (qplus + qminus); - Real r0avg = Real(0.5) * (rho_d_ptr[k] + rho_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]); - buoyancy_fab(i, j, k) = -qavg * r0avg * grav_gpu[2]; - }); - } // mfi - } // buoyancy_type -#endif + buoyancy_fab(i, j, k) = -qavg * r0avg * grav_gpu[2]; + }); + } // mfi + } // buoyancy_type + } // no moisture // ****************************************************************************************** // Moist versions of buoyancy expressions @@ -195,9 +184,9 @@ void make_buoyancy (Vector& S_data, // Base state density const Array4& r0_arr = r0->const_array(mfi); - const Array4 & qv_data = qvapor.array(mfi); - const Array4 & qc_data = qcloud.array(mfi); - const Array4 & qi_data = qice.array(mfi); + const Array4 & qv_data = qmoist->array(mfi,0); + const Array4 & qc_data = qmoist->array(mfi,1); + const Array4 & qi_data = qmoist->array(mfi,2); amrex::ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { @@ -262,9 +251,10 @@ void make_buoyancy (Vector& S_data, const Array4 & cell_data = S_data[IntVar::cons].array(mfi); const Array4 & cell_prim = S_prim.array(mfi); - const Array4 & qv_data = qvapor.array(mfi); - const Array4 & qc_data = qcloud.array(mfi); - const Array4 & qi_data = qice.array(mfi); + + const Array4 & qv_data = qmoist->const_array(mfi,0); + const Array4 & qc_data = qmoist->const_array(mfi,1); + const Array4 & qi_data = qmoist->const_array(mfi,2); amrex::ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { @@ -328,9 +318,9 @@ void make_buoyancy (Vector& S_data, const Array4 & cell_data = S_data[IntVar::cons].array(mfi); const Array4 & cell_prim = S_prim.array(mfi); - const Array4 & qv_data = qvapor.array(mfi); - const Array4 & qc_data = qcloud.array(mfi); - const Array4 & qi_data = qice.array(mfi); + const Array4 & qv_data = qmoist->const_array(mfi,0); + const Array4 & qc_data = qmoist->const_array(mfi,1); + const Array4 & qi_data = qmoist->const_array(mfi,2); amrex::ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { diff --git a/Source/TimeIntegration/ERF_make_condensation_source.cpp b/Source/TimeIntegration/ERF_make_condensation_source.cpp deleted file mode 100644 index 31cf7f0d0..000000000 --- a/Source/TimeIntegration/ERF_make_condensation_source.cpp +++ /dev/null @@ -1,77 +0,0 @@ -#include -#include -#include - -using namespace amrex; - -/** - * Function for computing the source terms for cloud vapor, cloud water and heat due to condensation. - * This routine is only called when using the warm moisture only formulation. - * - * @param[out] source the source terms for cloud vapor, cloud water and heat computed here - * @param[in] S current solution - * @param[in] tau_cond condensation process parameter - * @param[in] c_p Specific heat - */ - -#if defined(ERF_USE_WARM_NO_PRECIP) -void -ERF::condensation_source (MultiFab& source, MultiFab& S, Real tau_cond, Real c_p) -{ - for ( MFIter mfi(S, TilingIfNotGPU()); mfi.isValid(); ++mfi) { - - const Box& bx = mfi.tilebox(); - auto const& s_arr = S.const_array(mfi); - auto src_arr = source.array(mfi); - - // Coefficients and formula from Flatau et al (1992) - // "Polynomial Fits to Saturation Vapor Pressure", J. Applied Meteorology, 31(12), 1507-1513, 1992. - // DOI: https://doi.org/10.1175/1520-0450(1992)031<1507:pftsvp>2.0.co;2 - // as used in Munoz-Esparza et al (2022), JAMES - // "The FastEddy Resident-GPU Acclerated Large-Eddy Simulation Framework: - // Moist Dynamics Extension, Validation and Sensitivities of Modeling Non-Precipitating Shallow Cumulus Clouds" - // DOI: https://doi.org/10.1029/2021MS002904 - - Real g0 = -0.29912729e4; - Real g1 = -0.60170128e4; - Real g2 = 0.1887643854e2; - Real g3 = -0.28354721e-1; - Real g4 = 0.17838301e-4; - Real g5 = -0.84150417e-9; - Real g6 = 0.44412543e-12; - Real g7 = 0.2858487e1; - - ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept - { - Real rho_v = s_arr(i,j,k,RhoQv_comp); - Real rho_c = s_arr(i,j,k,RhoQc_comp); - - if ( (rho_v + rho_c) > 0.) { - Real rho_d = s_arr(i,j,k,Rho_comp); - Real rhotheta_d = s_arr(i,j,k,RhoTheta_comp); - Real theta_d = rhotheta_d / rho_d; - - Real T = getTgivenRandRTh(rho_d, rhotheta_d); - Real Tsq = getTgivenRandRTh(rho_d, rhotheta_d); - - Real ln_pvs = ( g0 + (g1 + (g2 + g7 * std::log(T) + (g3 + (g4 + (g5 + g6*T) * T) * T) * T) * T) ) / Tsq; - Real p_vs = std::exp(ln_pvs); // saturation vapor pressure computed using 8th-order polynomial fitting - Real rho_vs = p_vs / (R_d * T); - Real q_vs = rho_vs / rho_d; - - Real denom = 1.0 + (L_v * L_v * q_vs)/(c_p * R_v * Tsq); - Real f_cond_star = (rho_v - rho_vs) * denom; - - Real f_lim = rho_c; - Real f_cond = std::max(f_cond_star, -f_lim) / tau_cond; - - src_arr(i,j,k,RhoQv_comp) = -f_cond; - src_arr(i,j,k,RhoQc_comp) = f_cond; - - src_arr(i,j,k,RhoTheta_comp) = theta_d * L_v / (T * c_p) * f_cond; - } // net water > 0 - }); - - } // mfi -} -#endif diff --git a/Source/TimeIntegration/ERF_make_fast_coeffs.cpp b/Source/TimeIntegration/ERF_make_fast_coeffs.cpp index 634149aa5..c18e4b8f2 100644 --- a/Source/TimeIntegration/ERF_make_fast_coeffs.cpp +++ b/Source/TimeIntegration/ERF_make_fast_coeffs.cpp @@ -180,17 +180,12 @@ void make_fast_coeffs (int /*level*/, coeffP_a(i,j,k) = coeff_P; coeffQ_a(i,j,k) = coeff_Q; -#if defined(ERF_USE_MOISTURE) - Real q = 0.5 * ( prim(i,j,k,PrimQ1_comp) + prim(i,j,k-1,PrimQ1_comp) - +prim(i,j,k,PrimQ2_comp) + prim(i,j,k-1,PrimQ2_comp) ); - coeff_P /= (1.0 + q); - coeff_Q /= (1.0 + q); -#elif defined(ERF_USE_WARM_NO_PRECIP) - Real q = 0.5 * ( prim(i,j,k ,PrimQv_comp) + prim(i,j,k ,PrimQc_comp) + - prim(i,j,k-1,PrimQv_comp) + prim(i,j,k-1,PrimQc_comp) ); - coeff_P /= (1.0 + q); - coeff_Q /= (1.0 + q); -#endif + if (l_use_moisture) { + Real q = 0.5 * ( prim(i,j,k,PrimQ1_comp) + prim(i,j,k-1,PrimQ1_comp) + +prim(i,j,k,PrimQ2_comp) + prim(i,j,k-1,PrimQ2_comp) ); + coeff_P /= (1.0 + q); + coeff_Q /= (1.0 + q); + } Real theta_t_lo = 0.5 * ( prim(i,j,k-2,PrimTheta_comp) + prim(i,j,k-1,PrimTheta_comp) ); Real theta_t_mid = 0.5 * ( prim(i,j,k-1,PrimTheta_comp) + prim(i,j,k ,PrimTheta_comp) ); diff --git a/Source/TimeIntegration/ERF_slow_rhs_inc.cpp b/Source/TimeIntegration/ERF_slow_rhs_inc.cpp index 8d922fd3f..aaf81ed97 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_inc.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_inc.cpp @@ -78,9 +78,7 @@ void erf_slow_rhs_inc (int /*level*/, int nrk, const MultiFab& xvel, const MultiFab& yvel, const MultiFab& zvel, -#ifdef ERF_USE_MOISTURE - const MultiFab& qv, -#endif + const MultiFab* qv, std::unique_ptr& z_t_mf, MultiFab& Omega, const MultiFab& source, diff --git a/Source/TimeIntegration/ERF_slow_rhs_post.cpp b/Source/TimeIntegration/ERF_slow_rhs_post.cpp index 4814dc1e1..7160cbaae 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_post.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_post.cpp @@ -74,7 +74,7 @@ void erf_slow_rhs_post (int level, int finest_level, std::unique_ptr& mapfac_m, std::unique_ptr& mapfac_u, std::unique_ptr& mapfac_v, -#if defined(ERF_USE_NETCDF) && (defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP)) +#if defined(ERF_USE_NETCDF) const bool& moist_zero, const Real& bdy_time_interval, const Real& start_bdy_time, @@ -296,7 +296,7 @@ void erf_slow_rhs_post (int level, int finest_level, cur_prim, cell_rhs, detJ_arr, dxInv, mf_m, horiz_adv_type, vert_adv_type, l_use_terrain, flx_arr); - if (solverChoice.moisture_type != MoistureType::None) + if (solverChoice.moisture_type != MoistureType::None) { start_comp = RhoQ1_comp; num_comp = 2; @@ -392,7 +392,7 @@ void erf_slow_rhs_post (int level, int finest_level, } } #if defined(ERF_USE_NETCDF) - if (solverChoice.moisture_type != MoistureType::None) + if (solverChoice.moisture_type != MoistureType::None) // Zero moist RHS in set region if (moist_zero) { Box bx_xlo, bx_xhi, bx_ylo, bx_yhi; diff --git a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp index cd3e23245..4b673ac08 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp @@ -80,9 +80,7 @@ void erf_slow_rhs_pre (int level, int finest_level, const MultiFab& xvel, const MultiFab& yvel, const MultiFab& zvel, -#ifdef ERF_USE_MOISTURE - const MultiFab& qv, -#endif + const MultiFab* qv, std::unique_ptr& z_t_mf, MultiFab& Omega, const MultiFab& source, @@ -544,10 +542,8 @@ void erf_slow_rhs_pre (int level, int finest_level, Box gbx = mfi.tilebox(); gbx.grow(IntVect(1,1,0)); FArrayBox pprime; pprime.resize(gbx,1); Elixir pp_eli = pprime.elixir(); - const Array4 & pp_arr = pprime.array(); -#ifdef ERF_USE_MOISTURE - const Array4 & qv_arr = qv.const_array(mfi); -#endif + const Array4 & pp_arr = pprime.array(); + const Array4 & qv_arr = (use_moisture) ? qv->const_array(mfi) : Array4 {}; { BL_PROFILE("slow_rhs_pre_pprime"); ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept @@ -556,7 +552,7 @@ void erf_slow_rhs_pre (int level, int finest_level, // i,j,k,cell_data(i,j,k,RhoTheta_comp),cell_data(i,j,k+1,RhoTheta_comp)); AMREX_ASSERT(cell_data(i,j,k,RhoTheta_comp) > 0.); Real qv_for_p = 0.0; - if (use_moisture) { + if (use_moisture) { qv_for_p = qv_arr(i,j,k); } pp_arr(i,j,k) = getPgivenRTh(cell_data(i,j,k,RhoTheta_comp),qv_for_p) - p0_arr(i,j,k); @@ -802,7 +798,7 @@ void erf_slow_rhs_pre (int level, int finest_level, gpx *= mf_u(i,j,0); Real q = 0.0; - if (use_moisture) { + if (use_moisture) { q = 0.5 * ( cell_prim(i,j,k,PrimQ1_comp) + cell_prim(i-1,j,k,PrimQ1_comp) +cell_prim(i,j,k,PrimQ2_comp) + cell_prim(i-1,j,k,PrimQ2_comp) ); } @@ -844,7 +840,7 @@ void erf_slow_rhs_pre (int level, int finest_level, gpx *= mf_u(i,j,0); Real q = 0.0; - if (use_moisture) { + if (use_moisture) { q = 0.5 * ( cell_prim(i,j,k,PrimQ1_comp) + cell_prim(i-1,j,k,PrimQ1_comp) +cell_prim(i,j,k,PrimQ2_comp) + cell_prim(i-1,j,k,PrimQ2_comp) ); } @@ -906,7 +902,7 @@ void erf_slow_rhs_pre (int level, int finest_level, gpy *= mf_v(i,j,0); Real q = 0.0; - if (use_moisture) { + if (use_moisture) { q = 0.5 * ( cell_prim(i,j,k,PrimQ1_comp) + cell_prim(i,j-1,k,PrimQ1_comp) +cell_prim(i,j,k,PrimQ2_comp) + cell_prim(i,j-1,k,PrimQ2_comp) ); } @@ -945,7 +941,7 @@ void erf_slow_rhs_pre (int level, int finest_level, gpy *= mf_v(i,j,0); Real q = 0.0; - if (use_moisture) { + if (use_moisture) { q = 0.5 * ( cell_prim(i,j,k,PrimQ1_comp) + cell_prim(i,j-1,k,PrimQ1_comp) +cell_prim(i,j,k,PrimQ2_comp) + cell_prim(i,j-1,k,PrimQ2_comp) ); } @@ -997,7 +993,7 @@ void erf_slow_rhs_pre (int level, int finest_level, Real gpz = dxInv[2] * ( pp_arr(i,j,k)-pp_arr(i,j,k-1) ) / met_h_zeta; Real q = 0.0; - if (use_moisture) { + if (use_moisture) { q = 0.5 * ( cell_prim(i,j,k,PrimQ1_comp) + cell_prim(i,j,k-1,PrimQ1_comp) +cell_prim(i,j,k,PrimQ2_comp) + cell_prim(i,j,k-1,PrimQ2_comp) ); } @@ -1034,7 +1030,7 @@ void erf_slow_rhs_pre (int level, int finest_level, Real gpz = dxInv[2] * ( pp_arr(i,j,k)-pp_arr(i,j,k-1) ); Real q = 0.0; - if (use_moisture) { + if (use_moisture) { q = 0.5 * ( cell_prim(i,j,k,PrimQ1_comp) + cell_prim(i,j,k-1,PrimQ1_comp) +cell_prim(i,j,k,PrimQ2_comp) + cell_prim(i,j,k-1,PrimQ2_comp) ); } diff --git a/Source/TimeIntegration/Make.package b/Source/TimeIntegration/Make.package index 2a6bc3122..162bff16b 100644 --- a/Source/TimeIntegration/Make.package +++ b/Source/TimeIntegration/Make.package @@ -12,10 +12,6 @@ CEXE_sources += ERF_fast_rhs_T.cpp CEXE_sources += ERF_fast_rhs_MT.cpp CEXE_sources += ERF_ApplySpongeZoneBCs.cpp -ifeq ($(USE_WARM_NO_PRECIP),TRUE) -CEXE_sources += ERF_make_condensation_source.cpp -endif - ifeq ($(USE_POISSON_SOLVE),TRUE) CEXE_sources += ERF_slow_rhs_inc.cpp endif diff --git a/Source/TimeIntegration/TI_fast_rhs_fun.H b/Source/TimeIntegration/TI_fast_rhs_fun.H index 1ce3ce325..ff6289540 100644 --- a/Source/TimeIntegration/TI_fast_rhs_fun.H +++ b/Source/TimeIntegration/TI_fast_rhs_fun.H @@ -15,7 +15,7 @@ auto fast_rhs_fun = [&](int fast_step, int /*n_sub*/, int nrk, BL_PROFILE("fast_rhs_fun"); if (verbose) amrex::Print() << "Calling fast rhs at level " << level << " with dt = " << dtau << std::endl; - const bool use_moisture = (solverChoice.moisture_type != MoistureType::None); + const bool l_use_moisture = (solverChoice.moisture_type != MoistureType::None); // Define beta_s here so that it is consistent between where we make the fast coefficients // and where we use them @@ -86,7 +86,7 @@ auto fast_rhs_fun = [&](int fast_step, int /*n_sub*/, int nrk, // We have to call this each step since it depends on the substep time now // Note we pass in the *old* detJ here make_fast_coeffs(level, fast_coeffs, S_stage, S_prim, pi_stage, fine_geom, - use_moisture, solverChoice.use_terrain, solverChoice.gravity, solverChoice.c_p, + l_use_moisture, solverChoice.use_terrain, solverChoice.gravity, solverChoice.c_p, detJ_cc[level], r0, pi0, dtau, beta_s); if (fast_step == 0) { @@ -100,7 +100,7 @@ auto fast_rhs_fun = [&](int fast_step, int /*n_sub*/, int nrk, detJ_cc[level], detJ_cc_new[level], detJ_cc_src[level], dtau, beta_s, inv_fac, mapfac_m[level], mapfac_u[level], mapfac_v[level], - fr_as_crse, fr_as_fine, l_reflux); + fr_as_crse, fr_as_fine, l_use_moisture, l_reflux); } else { // If this is not the first substep we pass in S_data as the previous step's solution erf_fast_rhs_MT(fast_step, nrk, level, finest_level, @@ -112,14 +112,14 @@ auto fast_rhs_fun = [&](int fast_step, int /*n_sub*/, int nrk, detJ_cc[level], detJ_cc_new[level], detJ_cc_src[level], dtau, beta_s, inv_fac, mapfac_m[level], mapfac_u[level], mapfac_v[level], - fr_as_crse, fr_as_fine, l_reflux); + fr_as_crse, fr_as_fine, l_use_moisture, l_reflux); } } else if (solverChoice.use_terrain && solverChoice.terrain_type == TerrainType::Static) { if (fast_step == 0) { // If this is the first substep we make the coefficients since they are based only on stage data make_fast_coeffs(level, fast_coeffs, S_stage, S_prim, pi_stage, fine_geom, - use_moisture, solverChoice.use_terrain, solverChoice.gravity, solverChoice.c_p, + l_use_moisture, solverChoice.use_terrain, solverChoice.gravity, solverChoice.c_p, detJ_cc[level], r0, pi0, dtau, beta_s); // If this is the first substep we pass in S_old as the previous step's solution @@ -128,7 +128,7 @@ auto fast_rhs_fun = [&](int fast_step, int /*n_sub*/, int nrk, S_data, S_scratch, fine_geom, solverChoice.gravity, Omega, z_phys_nd[level], detJ_cc[level], dtau, beta_s, inv_fac, mapfac_m[level], mapfac_u[level], mapfac_v[level], - fr_as_crse, fr_as_fine, l_reflux); + fr_as_crse, fr_as_fine, l_use_moisture, l_reflux); } else { // If this is not the first substep we pass in S_data as the previous step's solution erf_fast_rhs_T(fast_step, nrk, level, finest_level, @@ -136,14 +136,14 @@ auto fast_rhs_fun = [&](int fast_step, int /*n_sub*/, int nrk, S_data, S_scratch, fine_geom, solverChoice.gravity, Omega, z_phys_nd[level], detJ_cc[level], dtau, beta_s, inv_fac, mapfac_m[level], mapfac_u[level], mapfac_v[level], - fr_as_crse, fr_as_fine, l_reflux); + fr_as_crse, fr_as_fine, l_use_moisture, l_reflux); } } else { if (fast_step == 0) { // If this is the first substep we make the coefficients since they are based only on stage data make_fast_coeffs(level, fast_coeffs, S_stage, S_prim, pi_stage, fine_geom, - use_moisture, solverChoice.use_terrain, solverChoice.gravity, solverChoice.c_p, + l_use_moisture, solverChoice.use_terrain, solverChoice.gravity, solverChoice.c_p, detJ_cc[level], r0, pi0, dtau, beta_s); // If this is the first substep we pass in S_old as the previous step's solution @@ -152,7 +152,7 @@ auto fast_rhs_fun = [&](int fast_step, int /*n_sub*/, int nrk, S_data, S_scratch, fine_geom, solverChoice.gravity, dtau, beta_s, inv_fac, mapfac_m[level], mapfac_u[level], mapfac_v[level], - fr_as_crse, fr_as_fine, l_reflux); + fr_as_crse, fr_as_fine, l_use_moisture, l_reflux); } else { // If this is not the first substep we pass in S_data as the previous step's solution erf_fast_rhs_N(fast_step, nrk, level, finest_level, @@ -160,7 +160,7 @@ auto fast_rhs_fun = [&](int fast_step, int /*n_sub*/, int nrk, S_data, S_scratch, fine_geom, solverChoice.gravity, dtau, beta_s, inv_fac, mapfac_m[level], mapfac_u[level], mapfac_v[level], - fr_as_crse, fr_as_fine, l_reflux); + fr_as_crse, fr_as_fine, l_use_moisture, l_reflux); } } diff --git a/Source/TimeIntegration/TI_headers.H b/Source/TimeIntegration/TI_headers.H index d8e64b01d..01e72c9b4 100644 --- a/Source/TimeIntegration/TI_headers.H +++ b/Source/TimeIntegration/TI_headers.H @@ -22,9 +22,7 @@ void erf_slow_rhs_pre (int level, int finest_level, int nrk, const amrex::MultiFab& xvel, const amrex::MultiFab& yvel, const amrex::MultiFab& zvel, -#if defined(ERF_USE_MOISTURE) - const amrex::MultiFab& qmoist, -#endif + const amrex::MultiFab* qmoist, std::unique_ptr& z_t, amrex::MultiFab& Omega, const amrex::MultiFab& source, @@ -91,7 +89,7 @@ void erf_slow_rhs_post (int level, int finest_level, int nrk, std::unique_ptr& mapfac_m, std::unique_ptr& mapfac_u, std::unique_ptr& mapfac_v, -#if defined(ERF_USE_NETCDF) && (defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP)) +#if defined(ERF_USE_NETCDF) const bool& moist_zero, const amrex::Real& bdy_time_interval, const amrex::Real& start_bdy_time, @@ -129,7 +127,7 @@ void erf_fast_rhs_N (int step, int nrk, int level, int finest_level, std::unique_ptr& mapfac_v, amrex::YAFluxRegister* fr_as_crse, amrex::YAFluxRegister* fr_as_fine, - bool l_reflux); + bool l_use_moisture, bool l_reflux); /** * Function for computing the fast RHS with fixed terrain @@ -156,7 +154,7 @@ void erf_fast_rhs_T (int step, int nrk, int level, int finest_level, std::unique_ptr& mapfac_v, amrex::YAFluxRegister* fr_as_crse, amrex::YAFluxRegister* fr_as_fine, - bool l_reflux); + bool l_use_moisture, bool l_reflux); /** * Function for computing the fast RHS with moving terrain @@ -190,7 +188,7 @@ void erf_fast_rhs_MT (int step, int nrk, int level, int finest_level, std::unique_ptr& mapfac_v, amrex::YAFluxRegister* fr_as_crse, amrex::YAFluxRegister* fr_as_fine, - bool l_reflux); + bool l_use_moisture, bool l_reflux); /** * Function for computing the coefficients for the tridiagonal solver used in the fast @@ -220,12 +218,10 @@ void make_fast_coeffs (int level, void make_buoyancy (amrex::Vector< amrex::MultiFab>& S_data, const amrex::MultiFab & S_prim, amrex::MultiFab& buoyancy, -#if defined(ERF_USE_MOISTURE) - const amrex::MultiFab& qmoist, + amrex::MultiFab* qmoist, const amrex::Gpu::DeviceVector qv_d, const amrex::Gpu::DeviceVector qc_d, const amrex::Gpu::DeviceVector qi_d, -#endif const amrex::Geometry geom, const SolverChoice& solverChoice, const amrex::MultiFab* r0); @@ -246,9 +242,7 @@ void erf_slow_rhs_inc (int level, int nrk, const amrex::MultiFab& xvel, const amrex::MultiFab& yvel, const amrex::MultiFab& zvel, -#if defined(ERF_USE_MOISTURE) - const amrex::MultiFab& qvapor, -#endif + const amrex::MultiFab* qvapor, std::unique_ptr& z_t, amrex::MultiFab& Omega, const amrex::MultiFab& source, diff --git a/Source/TimeIntegration/TI_slow_rhs_fun.H b/Source/TimeIntegration/TI_slow_rhs_fun.H index 3e28792fe..5ac193aa0 100644 --- a/Source/TimeIntegration/TI_slow_rhs_fun.H +++ b/Source/TimeIntegration/TI_slow_rhs_fun.H @@ -90,16 +90,12 @@ MultiFab* p0_new = &p_hse_new; make_buoyancy(S_data, S_prim, buoyancy, -#if defined(ERF_USE_MOISTURE) - qmoist[level], qv_d, qc_d, qi_d, -#endif + qmoist[level].get(), qv_d, qc_d, qi_d, fine_geom, solverChoice, r0_new); erf_slow_rhs_pre(level, finest_level, nrk, slow_dt, S_rhs, S_data, S_prim, S_scratch, xvel_new, yvel_new, zvel_new, -#if defined(ERF_USE_MOISTURE) - qmoist[level], -#endif + qmoist[level].get(), z_t_rk[level], Omega, source, buoyancy, Tau11, Tau22, Tau33, Tau12, Tau13, Tau21, Tau23, Tau31, Tau32, SmnSmn, eddyDiffs, Hfx3, Diss, @@ -187,16 +183,12 @@ // If not moving_terrain make_buoyancy(S_data, S_prim, buoyancy, -#if defined(ERF_USE_MOISTURE) - qmoist[level], qv_d, qc_d, qi_d, -#endif + qmoist[level].get(), qv_d, qc_d, qi_d, fine_geom, solverChoice, r0); erf_slow_rhs_pre(level, finest_level, nrk, slow_dt, S_rhs, S_data, S_prim, S_scratch, xvel_new, yvel_new, zvel_new, -#if defined(ERF_USE_MOISTURE) - qmoist[level], -#endif + qmoist[level].get(), z_t_rk[level], Omega, source, buoyancy, Tau11, Tau22, Tau33, Tau12, Tau13, Tau21, Tau23, Tau31, Tau32, SmnSmn, eddyDiffs, Hfx3, Diss, @@ -273,9 +265,7 @@ const Real new_stage_time, const int nrk) { -#if (!defined(ERF_USE_MOISTURE) && !defined(ERF_USE_WARM_NO_PRECIP)) amrex::ignore_unused(nrk); -#endif if (verbose) Print() << "Making slow rhs at time " << old_stage_time << " for slow variables advancing from " << old_step_time << " to " << new_stage_time << std::endl; @@ -285,11 +275,13 @@ // will be used to advance to Real slow_dt = new_stage_time - old_step_time; -#if defined(ERF_USE_NETCDF) && (defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP)) - // Flag for moisture relaxation zone bool moist_zero = false; - if (init_type=="real" && level==0 && wrfbdy_set_width>0) moist_zero = true; - if (init_type=="metgrid" && level==0 && metgrid_bdy_set_width > 0) moist_zero = true; +#if defined(ERF_USE_NETCDF) + if (solverChoice.moisture_type != MoistureType::None && level==0) { + // Flag for moisture relaxation zone + if (init_type=="real" && wrfbdy_set_width > 0) moist_zero = true; + if (init_type=="metgrid" && metgrid_bdy_set_width > 0) moist_zero = true; + } #endif // ************************************************************************* @@ -317,7 +309,7 @@ fine_geom, solverChoice, m_most, domain_bcs_type_d, z_phys_nd_src[level], detJ_cc[level], detJ_cc_new[level], mapfac_m[level], mapfac_u[level], mapfac_v[level], -#if defined(ERF_USE_NETCDF) && (defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP)) +#if defined(ERF_USE_NETCDF) moist_zero, bdy_time_interval, start_bdy_time, new_stage_time, wrfbdy_width-1, wrfbdy_set_width, bdy_data_xlo, bdy_data_xhi, bdy_data_ylo, bdy_data_yhi, @@ -333,7 +325,7 @@ fine_geom, solverChoice, m_most, domain_bcs_type_d, z_phys_nd[level], detJ_cc[level], detJ_cc[level], mapfac_m[level], mapfac_u[level], mapfac_v[level], -#if defined(ERF_USE_NETCDF) && (defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP)) +#if defined(ERF_USE_NETCDF) moist_zero, bdy_time_interval, start_bdy_time, new_stage_time, wrfbdy_width-1, wrfbdy_set_width, bdy_data_xlo, bdy_data_xhi, bdy_data_ylo, bdy_data_yhi, @@ -361,17 +353,13 @@ // If not moving_terrain make_buoyancy(S_data, S_prim, buoyancy, -#if defined(ERF_USE_MOISTURE) qmoist[level], qv_d, qc_d, qi_d, -#endif fine_geom, solverChoice, r0); erf_slow_rhs_inc(level, nrk, slow_dt, S_rhs, S_old, S_data, S_prim, S_scratch, xvel_new, yvel_new, zvel_new, -#if defined(ERF_USE_MOISTURE) qmoist[level], -#endif z_t_rk[level], Omega, source, buoyancy, Tau11, Tau22, Tau33, Tau12, Tau13, Tau21, Tau23, Tau31, Tau32, SmnSmn, eddyDiffs, Hfx3, Diss,