From 76e7512116985860afa9010e23113a0b357f60d5 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Mon, 11 Dec 2023 17:50:22 -0800 Subject: [PATCH] more cleanup of moisture stuff (#1338) --- Build/cmake_with_warm_no_precip.sh | 18 - CMake/BuildERFExe.cmake | 1 - Exec/ABL/prob.cpp | 10 +- Exec/ABL_input_sounding/prob.cpp | 10 +- Exec/DevTests/MovingTerrain/prob.cpp | 11 +- Exec/DevTests/MultiBlock/prob.cpp | 8 +- Exec/DevTests/ParticlesOverWoA/prob.cpp | 10 +- Exec/Radiation/prob.cpp | 30 +- Exec/RegTests/Bubble/prob.cpp | 14 +- Exec/RegTests/CouetteFlow/prob.cpp | 10 +- Exec/RegTests/DensityCurrent/prob.cpp | 34 +- Exec/RegTests/DynamicRefinement/prob.cpp | 7 +- Exec/RegTests/EkmanSpiral_custom/prob.cpp | 10 +- Exec/RegTests/GABLS1/prob.cpp | 4 +- Exec/RegTests/IsentropicVortex/prob.cpp | 7 +- Exec/RegTests/PoiseuilleFlow/prob.cpp | 10 +- Exec/RegTests/ScalarAdvDiff/prob.cpp | 10 +- Exec/RegTests/StokesSecondProblem/prob.cpp | 10 +- Exec/RegTests/TaylorGreenVortex/prob.cpp | 10 +- Exec/RegTests/Terrain2d_Cylinder/prob.cpp | 10 +- Exec/RegTests/Terrain3d_Hemisphere/prob.cpp | 10 +- Exec/RegTests/WitchOfAgnesi/prob.cpp | 10 +- Exec/SquallLine_2D/prob.cpp | 22 +- Exec/SuperCell/prob.cpp | 22 +- Source/BoundaryConditions/ABLMost.cpp | 2 +- .../BoundaryConditions_cons.cpp | 8 +- .../BoundaryConditions_xvel.cpp | 4 +- .../BoundaryConditions_yvel.cpp | 4 +- .../BoundaryConditions_zvel.cpp | 4 +- Source/BoundaryConditions/ERF_FillPatch.cpp | 6 +- Source/BoundaryConditions/ERF_PhysBCFunct.H | 8 +- Source/DataStructs/DataStruct.H | 42 +- .../Diffusion/ComputeTurbulentViscosity.cpp | 8 - Source/Diffusion/DiffusionSrcForState_N.cpp | 34 +- Source/Diffusion/DiffusionSrcForState_T.cpp | 34 +- Source/ERF.H | 17 +- Source/ERF.cpp | 43 +- Source/ERF_make_new_level.cpp | 31 +- Source/IO/Checkpoint.cpp | 22 +- Source/IO/ERF_ReadBndryPlanes.H | 12 +- Source/IO/ERF_ReadBndryPlanes.cpp | 19 +- Source/IO/ERF_WriteBndryPlanes.cpp | 2 +- Source/IO/NCCheckpoint.cpp | 4 +- Source/IO/Plotfile.cpp | 14 +- Source/IO/Plotfile.cpp.convergence | 4 - Source/IndexDefines.H | 41 +- Source/Initialization/ERF_init_bcs.cpp | 79 ++-- Source/TimeIntegration/ERF_Advance.cpp | 20 +- Source/TimeIntegration/ERF_MRI.H | 14 +- Source/TimeIntegration/ERF_advance_dycore.cpp | 67 +-- .../ERF_advance_microphysics.cpp | 16 +- Source/TimeIntegration/ERF_fast_rhs_MT.cpp | 48 +- Source/TimeIntegration/ERF_fast_rhs_N.cpp | 50 +-- Source/TimeIntegration/ERF_fast_rhs_T.cpp | 50 +-- Source/TimeIntegration/ERF_make_buoyancy.cpp | 410 +++++++++--------- .../ERF_make_condensation_source.cpp | 77 ---- .../TimeIntegration/ERF_make_fast_coeffs.cpp | 35 +- Source/TimeIntegration/ERF_slow_rhs_inc.cpp | 4 +- Source/TimeIntegration/ERF_slow_rhs_post.cpp | 79 ++-- Source/TimeIntegration/ERF_slow_rhs_pre.cpp | 147 +++---- Source/TimeIntegration/Make.package | 4 - Source/TimeIntegration/TI_fast_rhs_fun.H | 20 +- Source/TimeIntegration/TI_headers.H | 21 +- Source/TimeIntegration/TI_slow_rhs_fun.H | 40 +- Source/TimeIntegration/TI_utils.H | 5 +- 65 files changed, 810 insertions(+), 1037 deletions(-) delete mode 100755 Build/cmake_with_warm_no_precip.sh delete mode 100644 Source/TimeIntegration/ERF_make_condensation_source.cpp 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/BoundaryConditions/ABLMost.cpp b/Source/BoundaryConditions/ABLMost.cpp index e088eeee3..cedb9202e 100644 --- a/Source/BoundaryConditions/ABLMost.cpp +++ b/Source/BoundaryConditions/ABLMost.cpp @@ -176,7 +176,7 @@ ABLMost::compute_most_bcs (const int& lev, if (var_idx == Vars::cons) { Box b2d = bx; b2d.setBig(2,zlo-1); - int n = Cons::RhoTheta; + int n = RhoTheta_comp; ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { diff --git a/Source/BoundaryConditions/BoundaryConditions_cons.cpp b/Source/BoundaryConditions/BoundaryConditions_cons.cpp index 7f07b4d27..0dab4b9ca 100644 --- a/Source/BoundaryConditions/BoundaryConditions_cons.cpp +++ b/Source/BoundaryConditions/BoundaryConditions_cons.cpp @@ -43,12 +43,12 @@ void ERFPhysBCFunct::impose_lateral_cons_bcs (const Array4& dest_arr, cons #endif const amrex::BCRec* bc_ptr = bcrs_d.data(); - GpuArray,AMREX_SPACEDIM+NVAR> l_bc_extdir_vals_d; + GpuArray,AMREX_SPACEDIM+NVAR_max> l_bc_extdir_vals_d; for (int i = 0; i < icomp+ncomp; i++) for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) l_bc_extdir_vals_d[i][ori] = m_bc_extdir_vals[bccomp+i][ori]; - GpuArray,AMREX_SPACEDIM+NVAR> l_bc_neumann_vals_d; + GpuArray,AMREX_SPACEDIM+NVAR_max> l_bc_neumann_vals_d; for (int i = 0; i < icomp+ncomp; i++) for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) l_bc_neumann_vals_d[i][ori] = m_bc_neumann_vals[bccomp+i][ori]; @@ -202,12 +202,12 @@ void ERFPhysBCFunct::impose_vertical_cons_bcs (const Array4& dest_arr, con #endif const amrex::BCRec* bc_ptr = bcrs_d.data(); - GpuArray,AMREX_SPACEDIM+NVAR> l_bc_extdir_vals_d; + GpuArray,AMREX_SPACEDIM+NVAR_max> l_bc_extdir_vals_d; for (int i = 0; i < icomp+ncomp; i++) for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) l_bc_extdir_vals_d[i][ori] = m_bc_extdir_vals[bccomp+i][ori]; - GpuArray,AMREX_SPACEDIM+NVAR> l_bc_neumann_vals_d; + GpuArray,AMREX_SPACEDIM+NVAR_max> l_bc_neumann_vals_d; for (int i = 0; i < icomp+ncomp; i++) for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) l_bc_neumann_vals_d[i][ori] = m_bc_neumann_vals[bccomp+i][ori]; diff --git a/Source/BoundaryConditions/BoundaryConditions_xvel.cpp b/Source/BoundaryConditions/BoundaryConditions_xvel.cpp index be81422f1..6bdcb5096 100644 --- a/Source/BoundaryConditions/BoundaryConditions_xvel.cpp +++ b/Source/BoundaryConditions/BoundaryConditions_xvel.cpp @@ -44,7 +44,7 @@ void ERFPhysBCFunct::impose_lateral_xvel_bcs (const Array4& dest_arr, #endif const amrex::BCRec* bc_ptr = bcrs_d.data(); - GpuArray,AMREX_SPACEDIM+NVAR> l_bc_extdir_vals_d; + GpuArray,AMREX_SPACEDIM+NVAR_max> l_bc_extdir_vals_d; for (int i = 0; i < ncomp; i++) for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) @@ -179,7 +179,7 @@ void ERFPhysBCFunct::impose_vertical_xvel_bcs (const Array4& dest_arr, #endif const amrex::BCRec* bc_ptr = bcrs_d.data(); - GpuArray,AMREX_SPACEDIM+NVAR> l_bc_extdir_vals_d; + GpuArray,AMREX_SPACEDIM+NVAR_max> l_bc_extdir_vals_d; for (int i = 0; i < ncomp; i++) for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) diff --git a/Source/BoundaryConditions/BoundaryConditions_yvel.cpp b/Source/BoundaryConditions/BoundaryConditions_yvel.cpp index ddb8cc23f..a7fd08b80 100644 --- a/Source/BoundaryConditions/BoundaryConditions_yvel.cpp +++ b/Source/BoundaryConditions/BoundaryConditions_yvel.cpp @@ -41,7 +41,7 @@ void ERFPhysBCFunct::impose_lateral_yvel_bcs (const Array4& dest_arr, #endif const amrex::BCRec* bc_ptr = bcrs_d.data(); - GpuArray, AMREX_SPACEDIM+NVAR> l_bc_extdir_vals_d; + GpuArray, AMREX_SPACEDIM+NVAR_max> l_bc_extdir_vals_d; for (int i = 0; i < ncomp; i++) for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) @@ -178,7 +178,7 @@ void ERFPhysBCFunct::impose_vertical_yvel_bcs (const Array4& dest_arr, #endif const amrex::BCRec* bc_ptr = bcrs_d.data(); - GpuArray, AMREX_SPACEDIM+NVAR> l_bc_extdir_vals_d; + GpuArray, AMREX_SPACEDIM+NVAR_max> l_bc_extdir_vals_d; for (int i = 0; i < ncomp; i++) for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) diff --git a/Source/BoundaryConditions/BoundaryConditions_zvel.cpp b/Source/BoundaryConditions/BoundaryConditions_zvel.cpp index e320c6c42..9ee5f8288 100644 --- a/Source/BoundaryConditions/BoundaryConditions_zvel.cpp +++ b/Source/BoundaryConditions/BoundaryConditions_zvel.cpp @@ -47,7 +47,7 @@ void ERFPhysBCFunct::impose_lateral_zvel_bcs (const Array4& dest_arr, #endif const amrex::BCRec* bc_ptr_w = bcrs_w_d.data(); - GpuArray,AMREX_SPACEDIM+NVAR> l_bc_extdir_vals_d; + GpuArray,AMREX_SPACEDIM+NVAR_max> l_bc_extdir_vals_d; bool l_use_terrain = (m_z_phys_nd != nullptr); @@ -200,7 +200,7 @@ void ERFPhysBCFunct::impose_vertical_zvel_bcs (const Array4& dest_arr, bool l_use_terrain = (m_z_phys_nd != nullptr); bool l_moving_terrain = (terrain_type == TerrainType::Moving); - GpuArray,AMREX_SPACEDIM+NVAR> l_bc_extdir_vals_d; + GpuArray,AMREX_SPACEDIM+NVAR_max> l_bc_extdir_vals_d; for (int i = 0; i < ncomp; i++) for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) diff --git a/Source/BoundaryConditions/ERF_FillPatch.cpp b/Source/BoundaryConditions/ERF_FillPatch.cpp index b12c33fdf..3ae8182cf 100644 --- a/Source/BoundaryConditions/ERF_FillPatch.cpp +++ b/Source/BoundaryConditions/ERF_FillPatch.cpp @@ -295,17 +295,17 @@ ERF::FillCoarsePatch (int lev, Real time, const Vector& mfs) } else if (var_idx == Vars::xvel || var_idx == Vars::xmom) { - bccomp = NVAR; + bccomp = BCVars::xvel_bc; mapper = &face_linear_interp; } else if (var_idx == Vars::yvel || var_idx == Vars::ymom) { - bccomp = NVAR+1; + bccomp = BCVars::yvel_bc; mapper = &face_linear_interp; } else if (var_idx == Vars::zvel || var_idx == Vars::zmom) { - bccomp = NVAR+2; + bccomp = BCVars::zvel_bc; mapper = &face_linear_interp; } diff --git a/Source/BoundaryConditions/ERF_PhysBCFunct.H b/Source/BoundaryConditions/ERF_PhysBCFunct.H index c7790e0f1..f660c2ecf 100644 --- a/Source/BoundaryConditions/ERF_PhysBCFunct.H +++ b/Source/BoundaryConditions/ERF_PhysBCFunct.H @@ -26,8 +26,8 @@ public: const amrex::Geometry& geom, const amrex::Vector& domain_bcs_type, const amrex::Gpu::DeviceVector& domain_bcs_type_d, const TerrainType& terrain_type, - amrex::Array,AMREX_SPACEDIM+NVAR> bc_extdir_vals, - amrex::Array,AMREX_SPACEDIM+NVAR> bc_neumann_vals, + amrex::Array,AMREX_SPACEDIM+NVAR_max> bc_extdir_vals, + amrex::Array,AMREX_SPACEDIM+NVAR_max> bc_neumann_vals, std::unique_ptr& z_phys_nd, std::unique_ptr& detJ_cc) : m_lev(lev), @@ -109,8 +109,8 @@ private: amrex::Vector m_domain_bcs_type; amrex::Gpu::DeviceVector m_domain_bcs_type_d; TerrainType m_terrain_type; - amrex::Array,AMREX_SPACEDIM+NVAR> m_bc_extdir_vals; - amrex::Array,AMREX_SPACEDIM+NVAR> m_bc_neumann_vals; + amrex::Array,AMREX_SPACEDIM+NVAR_max> m_bc_extdir_vals; + amrex::Array,AMREX_SPACEDIM+NVAR_max> m_bc_neumann_vals; amrex::MultiFab* m_z_phys_nd; amrex::MultiFab* m_detJ_cc; }; 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/Diffusion/ComputeTurbulentViscosity.cpp b/Source/Diffusion/ComputeTurbulentViscosity.cpp index 029b36c71..1a65ac245 100644 --- a/Source/Diffusion/ComputeTurbulentViscosity.cpp +++ b/Source/Diffusion/ComputeTurbulentViscosity.cpp @@ -169,16 +169,8 @@ void ComputeTurbulentViscosityLES (const amrex::MultiFab& Tau11, const amrex::Mu Real inv_Pr_t = turbChoice.Pr_t_inv; Real inv_Sc_t = turbChoice.Sc_t_inv; Real inv_sigma_k = 1.0 / turbChoice.sigma_k; -#if defined(ERF_USE_MOISTURE) // EddyDiff mapping : Theta_h Scalar_h KE_h QKE_h Q1_h Q2_h Vector Factors = {inv_Pr_t, inv_Sc_t, inv_sigma_k, inv_sigma_k, inv_Sc_t, inv_Sc_t}; // alpha = mu/Pr -#elif defined(ERF_USE_WARM_NO_PRECIP) - // EddyDiff mapping : Theta_h Scalar_h KE_h QKE_h Qv_h Qc_h - Vector Factors = {inv_Pr_t, inv_Sc_t, inv_sigma_k, inv_sigma_k, inv_Sc_t, inv_Sc_t}; // alpha = mu/Pr -#else - // EddyDiff mapping : Theta_h Scalar_h KE_h QKE_h - Vector Factors = {inv_Pr_t, inv_Sc_t, inv_sigma_k, inv_sigma_k}; // alpha = mu/Pr -#endif Gpu::AsyncVector d_Factors; d_Factors.resize(Factors.size()); Gpu::copy(Gpu::hostToDevice, Factors.begin(), Factors.end(), d_Factors.begin()); Real* fac_ptr = d_Factors.data(); diff --git a/Source/Diffusion/DiffusionSrcForState_N.cpp b/Source/Diffusion/DiffusionSrcForState_N.cpp index eb097f984..da68d7ac9 100644 --- a/Source/Diffusion/DiffusionSrcForState_N.cpp +++ b/Source/Diffusion/DiffusionSrcForState_N.cpp @@ -84,9 +84,9 @@ DiffusionSrcForState_N (const amrex::Box& bx, const amrex::Box& domain, const int qty_offset = RhoTheta_comp; // Theta, KE, QKE, Scalar - Vector alpha_eff(NUM_PRIM, 0.0); + Vector alpha_eff(NVAR_max, 0.0); if (l_consA) { - for (int i = 0; i < NUM_PRIM; ++i) { + for (int i = 0; i < NVAR_max; ++i) { switch (i) { case PrimTheta_comp: alpha_eff[PrimTheta_comp] = diffChoice.alpha_T; @@ -94,28 +94,19 @@ DiffusionSrcForState_N (const amrex::Box& bx, const amrex::Box& domain, case PrimScalar_comp: alpha_eff[PrimScalar_comp] = diffChoice.alpha_C; break; -#if defined(ERF_USE_MOISTURE) case PrimQ1_comp: alpha_eff[PrimQ1_comp] = diffChoice.alpha_C; break; case PrimQ2_comp: alpha_eff[PrimQ2_comp] = diffChoice.alpha_C; break; -#elif defined(ERF_USE_WARM_NO_PRECIP) - case PrimQv_comp: - alpha_eff[PrimQv_comp] = diffChoice.alpha_C; - break; - case PrimQc_comp: - alpha_eff[PrimQc_comp] = diffChoice.alpha_C; - break; -#endif default: alpha_eff[i] = 0.0; break; } } } else { - for (int i = 0; i < NUM_PRIM; ++i) { + for (int i = 0; i < NVAR_max; ++i) { switch (i) { case PrimTheta_comp: alpha_eff[PrimTheta_comp] = diffChoice.rhoAlpha_T; @@ -123,40 +114,21 @@ DiffusionSrcForState_N (const amrex::Box& bx, const amrex::Box& domain, case PrimScalar_comp: alpha_eff[PrimScalar_comp] = diffChoice.rhoAlpha_C; break; -#if defined(ERF_USE_MOISTURE) case PrimQ1_comp: alpha_eff[PrimQ1_comp] = diffChoice.rhoAlpha_C; break; case PrimQ2_comp: alpha_eff[PrimQ2_comp] = diffChoice.rhoAlpha_C; break; -#elif defined(ERF_USE_WARM_NO_PRECIP) - case PrimQv_comp: - alpha_eff[PrimQv_comp] = diffChoice.rhoAlpha_C; - break; - case PrimQc_comp: - alpha_eff[PrimQc_comp] = diffChoice.rhoAlpha_C; - break; -#endif default: alpha_eff[i] = 0.0; break; } } } -#if defined(ERF_USE_MOISTURE) Vector eddy_diff_idx{EddyDiff::Theta_h, EddyDiff::KE_h, EddyDiff::QKE_h, EddyDiff::Scalar_h, EddyDiff::Q1_h, EddyDiff::Q2_h}; Vector eddy_diff_idy{EddyDiff::Theta_h, EddyDiff::KE_h, EddyDiff::QKE_h, EddyDiff::Scalar_h, EddyDiff::Q1_h, EddyDiff::Q2_h}; Vector eddy_diff_idz{EddyDiff::Theta_v, EddyDiff::KE_v, EddyDiff::QKE_v, EddyDiff::Scalar_v, EddyDiff::Q1_v, EddyDiff::Q2_v}; -#elif defined(ERF_USE_WARM_NO_PRECIP) - Vector eddy_diff_idx{EddyDiff::Theta_h, EddyDiff::KE_h, EddyDiff::QKE_h, EddyDiff::Scalar_h, EddyDiff::Qv_h, EddyDiff::Qc_h}; - Vector eddy_diff_idy{EddyDiff::Theta_h, EddyDiff::KE_h, EddyDiff::QKE_h, EddyDiff::Scalar_h, EddyDiff::Qv_h, EddyDiff::Qc_h}; - Vector eddy_diff_idz{EddyDiff::Theta_v, EddyDiff::KE_v, EddyDiff::QKE_v, EddyDiff::Scalar_v, EddyDiff::Qv_v, EddyDiff::Qc_v}; -#else - Vector eddy_diff_idx{EddyDiff::Theta_h, EddyDiff::KE_h, EddyDiff::QKE_h, EddyDiff::Scalar_h}; - Vector eddy_diff_idy{EddyDiff::Theta_h, EddyDiff::KE_h, EddyDiff::QKE_h, EddyDiff::Scalar_h}; - Vector eddy_diff_idz{EddyDiff::Theta_v, EddyDiff::KE_v, EddyDiff::QKE_v, EddyDiff::Scalar_v}; -#endif // Device vectors Gpu::AsyncVector alpha_eff_d; diff --git a/Source/Diffusion/DiffusionSrcForState_T.cpp b/Source/Diffusion/DiffusionSrcForState_T.cpp index 884937916..7bbb113dd 100644 --- a/Source/Diffusion/DiffusionSrcForState_T.cpp +++ b/Source/Diffusion/DiffusionSrcForState_T.cpp @@ -91,9 +91,9 @@ DiffusionSrcForState_T (const amrex::Box& bx, const amrex::Box& domain, const int qty_offset = RhoTheta_comp; // Theta, KE, QKE, Scalar - Vector alpha_eff(NUM_PRIM, 0.0); + Vector alpha_eff(NVAR_max, 0.0); if (l_consA) { - for (int i = 0; i < NUM_PRIM; ++i) { + for (int i = 0; i < NVAR_max; ++i) { switch (i) { case PrimTheta_comp: alpha_eff[PrimTheta_comp] = diffChoice.alpha_T; @@ -101,28 +101,19 @@ DiffusionSrcForState_T (const amrex::Box& bx, const amrex::Box& domain, case PrimScalar_comp: alpha_eff[PrimScalar_comp] = diffChoice.alpha_C; break; -#if defined(ERF_USE_MOISTURE) case PrimQ1_comp: alpha_eff[PrimQ1_comp] = diffChoice.alpha_C; break; case PrimQ2_comp: alpha_eff[PrimQ2_comp] = diffChoice.alpha_C; break; -#elif defined(ERF_USE_WARM_NO_PRECIP) - case PrimQv_comp: - alpha_eff[PrimQv_comp] = diffChoice.alpha_C; - break; - case PrimQc_comp: - alpha_eff[PrimQc_comp] = diffChoice.alpha_C; - break; -#endif default: alpha_eff[i] = 0.0; break; } } } else { - for (int i = 0; i < NUM_PRIM; ++i) { + for (int i = 0; i < NVAR_max; ++i) { switch (i) { case PrimTheta_comp: alpha_eff[PrimTheta_comp] = diffChoice.rhoAlpha_T; @@ -130,21 +121,12 @@ DiffusionSrcForState_T (const amrex::Box& bx, const amrex::Box& domain, case PrimScalar_comp: alpha_eff[PrimScalar_comp] = diffChoice.rhoAlpha_C; break; -#if defined(ERF_USE_MOISTURE) case PrimQ1_comp: alpha_eff[PrimQ1_comp] = diffChoice.rhoAlpha_C; break; case PrimQ2_comp: alpha_eff[PrimQ2_comp] = diffChoice.rhoAlpha_C; break; -#elif defined(ERF_USE_WARM_NO_PRECIP) - case PrimQv_comp: - alpha_eff[PrimQv_comp] = diffChoice.rhoAlpha_C; - break; - case PrimQc_comp: - alpha_eff[PrimQc_comp] = diffChoice.rhoAlpha_C; - break; -#endif default: alpha_eff[i] = 0.0; break; @@ -152,19 +134,9 @@ DiffusionSrcForState_T (const amrex::Box& bx, const amrex::Box& domain, } } -#if defined(ERF_USE_MOISTURE) Vector eddy_diff_idx{EddyDiff::Theta_h, EddyDiff::KE_h, EddyDiff::QKE_h, EddyDiff::Scalar_h, EddyDiff::Q1_h, EddyDiff::Q2_h}; Vector eddy_diff_idy{EddyDiff::Theta_h, EddyDiff::KE_h, EddyDiff::QKE_h, EddyDiff::Scalar_h, EddyDiff::Q1_h, EddyDiff::Q2_h}; Vector eddy_diff_idz{EddyDiff::Theta_v, EddyDiff::KE_v, EddyDiff::QKE_v, EddyDiff::Scalar_v, EddyDiff::Q1_v, EddyDiff::Q2_v}; -#elif defined(ERF_USE_WARM_NO_PRECIP) - Vector eddy_diff_idx{EddyDiff::Theta_h, EddyDiff::KE_h, EddyDiff::QKE_h, EddyDiff::Scalar_h, EddyDiff::Qv_h, EddyDiff::Qc_h}; - Vector eddy_diff_idy{EddyDiff::Theta_h, EddyDiff::KE_h, EddyDiff::QKE_h, EddyDiff::Scalar_h, EddyDiff::Qv_h, EddyDiff::Qc_h}; - Vector eddy_diff_idz{EddyDiff::Theta_v, EddyDiff::KE_v, EddyDiff::QKE_v, EddyDiff::Scalar_v, EddyDiff::Qv_v, EddyDiff::Qc_v}; -#else - Vector eddy_diff_idx{EddyDiff::Theta_h, EddyDiff::KE_h, EddyDiff::QKE_h, EddyDiff::Scalar_h}; - Vector eddy_diff_idy{EddyDiff::Theta_h, EddyDiff::KE_h, EddyDiff::QKE_h, EddyDiff::Scalar_h}; - Vector eddy_diff_idz{EddyDiff::Theta_v, EddyDiff::KE_v, EddyDiff::QKE_v, EddyDiff::Scalar_v}; -#endif // Device vectors Gpu::AsyncVector alpha_eff_d; diff --git a/Source/ERF.H b/Source/ERF.H index ea9c52dee..8e3735709 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -303,13 +303,13 @@ public: void fill_from_wrfbdy (const amrex::Vector& mfs, amrex::Real time, bool cons_only = false, - int icomp_cons = 0, - int ncomp_cons = NVAR); + int icomp_cons, + int ncomp_cons); void fill_from_metgrid (const amrex::Vector& mfs, amrex::Real time, bool cons_only = false, - int icomp_cons = 0, - int ncomp_cons = NVAR); + int icomp_cons, + int ncomp_cons); #endif #ifdef ERF_USE_NETCDF @@ -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); @@ -594,10 +589,10 @@ private: amrex::Array domain_bc_type; // These hold the Dirichlet values at walls which need them ... - amrex::Array,AMREX_SPACEDIM+NVAR> m_bc_extdir_vals; + amrex::Array, AMREX_SPACEDIM+NVAR_max> m_bc_extdir_vals; // These hold the Neumann values at walls which need them ... - amrex::Array,AMREX_SPACEDIM+NVAR> m_bc_neumann_vals; + amrex::Array, AMREX_SPACEDIM+NVAR_max> m_bc_neumann_vals; // These are the "physical" boundary condition types (e.g. "inflow") amrex::GpuArray phys_bc_type; diff --git a/Source/ERF.cpp b/Source/ERF.cpp index e2b51b376..4bd5a6aff 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -319,7 +319,7 @@ ERF::post_timestep (int nstep, Real time, Real dt_lev0) if (solverChoice.coupling_type == CouplingType::TwoWay) { int src_comp_reflux = 0; - int num_comp_reflux = NVAR; + int num_comp_reflux = vars_new[0][Vars::cons].nComp(); bool use_terrain = solverChoice.use_terrain; for (int lev = finest_level-1; lev >= 0; lev--) { @@ -531,8 +531,8 @@ ERF::InitData () if (solverChoice.coupling_type == CouplingType::TwoWay) { - int src_comp_reflux = 0; - int num_comp_reflux = NVAR; + int src_comp_reflux = 0; + int num_comp_reflux = vars_new[0][Vars::cons].nComp(); AverageDown(src_comp_reflux, num_comp_reflux); } @@ -568,12 +568,13 @@ ERF::InitData () // Initialize flux registers (whether we start from scratch or restart) if (solverChoice.coupling_type == CouplingType::TwoWay) { advflux_reg[0] = nullptr; + int ncomp_reflux = vars_new[0][Vars::cons].nComp(); for (int lev = 1; lev <= finest_level; lev++) { advflux_reg[lev] = new YAFluxRegister(grids[lev], grids[lev-1], dmap[lev], dmap[lev-1], geom[lev], geom[lev-1], - ref_ratio[lev-1], lev, NVAR); + ref_ratio[lev-1], lev, ncomp_reflux); } } @@ -612,8 +613,8 @@ ERF::InitData () { amrex::IntVect ng = IntVect(0,0,0); MultiFab S(vars_new[lev][Vars::cons],make_alias,0,2); - MultiFab::Copy( *Theta_prim[lev], S, Cons::RhoTheta, 0, 1, ng); - MultiFab::Divide(*Theta_prim[lev], S, Cons::Rho , 0, 1, ng); + MultiFab::Copy( *Theta_prim[lev], S, RhoTheta_comp, 0, 1, ng); + MultiFab::Divide(*Theta_prim[lev], S, Rho_comp , 0, 1, ng); m_most->update_mac_ptrs(lev, vars_new, Theta_prim); m_most->update_fluxes(lev, t_new[lev]); } @@ -664,10 +665,12 @@ ERF::InitData () auto& lev_new = vars_new[lev]; auto& lev_old = vars_old[lev]; - MultiFab::Copy(lev_old[Vars::cons],lev_new[Vars::cons],0,0,NVAR,lev_new[Vars::cons].nGrowVect()); - MultiFab::Copy(lev_old[Vars::xvel],lev_new[Vars::xvel],0,0, 1,lev_new[Vars::xvel].nGrowVect()); - MultiFab::Copy(lev_old[Vars::yvel],lev_new[Vars::yvel],0,0, 1,lev_new[Vars::yvel].nGrowVect()); - MultiFab::Copy(lev_old[Vars::zvel],lev_new[Vars::zvel],0,0, 1,lev_new[Vars::zvel].nGrowVect()); + int ncomp = lev_new[Vars::cons].nComp(); + + MultiFab::Copy(lev_old[Vars::cons],lev_new[Vars::cons],0,0,ncomp,lev_new[Vars::cons].nGrowVect()); + MultiFab::Copy(lev_old[Vars::xvel],lev_new[Vars::xvel],0,0, 1,lev_new[Vars::xvel].nGrowVect()); + MultiFab::Copy(lev_old[Vars::yvel],lev_new[Vars::yvel],0,0, 1,lev_new[Vars::yvel].nGrowVect()); + MultiFab::Copy(lev_old[Vars::zvel],lev_new[Vars::zvel],0,0, 1,lev_new[Vars::zvel].nGrowVect()); } // Compute the minimum dz in the domain (to be used for setting the timestep) @@ -1145,7 +1148,7 @@ ERF::MakeHorizontalAverages () // First, average down all levels (if doing two-way coupling) if (solverChoice.coupling_type == CouplingType::TwoWay) { int src_comp_reflux = 0; - int num_comp_reflux = NVAR; + int num_comp_reflux = vars_new[lev][Vars::cons].nComp(); AverageDown(src_comp_reflux, num_comp_reflux); } @@ -1161,11 +1164,11 @@ ERF::MakeHorizontalAverages () auto fab_arr = mf.array(mfi); auto const cons_arr = vars_new[lev][Vars::cons].const_array(mfi); ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - Real dens = cons_arr(i, j, k, Cons::Rho); + Real dens = cons_arr(i, j, k, Rho_comp); fab_arr(i, j, k, 0) = dens; - fab_arr(i, j, k, 1) = cons_arr(i, j, k, Cons::RhoTheta) / dens; + fab_arr(i, j, k, 1) = cons_arr(i, j, k, RhoTheta_comp) / dens; if (!use_moisture) { - fab_arr(i, j, k, 2) = getPgivenRTh(cons_arr(i, j, k, Cons::RhoTheta)); + fab_arr(i, j, k, 2) = getPgivenRTh(cons_arr(i, j, k, RhoTheta_comp)); } }); } @@ -1181,10 +1184,10 @@ ERF::MakeHorizontalAverages () auto const qv_arr = qv.const_array(mfi); ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - Real dens = cons_arr(i, j, k, Cons::Rho); - fab_arr(i, j, k, 2) = getPgivenRTh(cons_arr(i, j, k, Cons::RhoTheta), qv_arr(i,j,k)); - fab_arr(i, j, k, 3) = cons_arr(i, j, k, Cons::RhoQ1) / dens; - fab_arr(i, j, k, 4) = cons_arr(i, j, k, Cons::RhoQ2) / dens; + Real dens = cons_arr(i, j, k, Rho_comp); + fab_arr(i, j, k, 2) = getPgivenRTh(cons_arr(i, j, k, RhoTheta_comp), qv_arr(i,j,k)); + fab_arr(i, j, k, 3) = cons_arr(i, j, k, RhoQ1_comp) / dens; + fab_arr(i, j, k, 4) = cons_arr(i, j, k, RhoQ2_comp) / dens; }); } @@ -1373,7 +1376,7 @@ ERF::Construct_ERFFillPatchers (int lev) auto& dm_fine = fine_new[Vars::cons].DistributionMap(); auto& dm_crse = crse_new[Vars::cons].DistributionMap(); - int ncomp = NVAR; + int ncomp = vars_new[lev][Vars::cons].nComp(); FPr_c.emplace_back(ba_fine, dm_fine, geom[lev] , ba_crse, dm_crse, geom[lev-1], @@ -1401,7 +1404,7 @@ ERF::Define_ERFFillPatchers (int lev) auto& dm_fine = fine_new[Vars::cons].DistributionMap(); auto& dm_crse = crse_new[Vars::cons].DistributionMap(); - int ncomp = NVAR; + int ncomp = fine_new[Vars::cons].nComp(); FPr_c[lev-1].Define(ba_fine, dm_fine, geom[lev] , ba_crse, dm_crse, geom[lev-1], diff --git a/Source/ERF_make_new_level.cpp b/Source/ERF_make_new_level.cpp index b48b41b92..cbd27688e 100644 --- a/Source/ERF_make_new_level.cpp +++ b/Source/ERF_make_new_level.cpp @@ -44,8 +44,10 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba, init_stuff(lev, ba, dm); - lev_new[Vars::cons].define(ba, dm, Cons::NumVars, ngrow_state); - lev_old[Vars::cons].define(ba, dm, Cons::NumVars, ngrow_state); + int ncomp_cons = (solverChoice.moisture_type == MoistureType::None) ? NVAR_max-2 : NVAR_max; + + lev_new[Vars::cons].define(ba, dm, ncomp_cons, ngrow_state); + lev_old[Vars::cons].define(ba, dm, ncomp_cons, ngrow_state); lev_new[Vars::xvel].define(convert(ba, IntVect(1,0,0)), dm, 1, ngrow_vels); lev_old[Vars::xvel].define(convert(ba, IntVect(1,0,0)), dm, 1, ngrow_vels); @@ -176,6 +178,8 @@ ERF::MakeNewLevelFromCoarse (int lev, Real time, const BoxArray& ba, auto& lev_new = vars_new[lev]; auto& lev_old = vars_old[lev]; + int ncomp = lev_new[Vars::cons].nComp(); + // ******************************************************************************************** // These are the persistent containers for the old and new data // ******************************************************************************************** @@ -228,7 +232,7 @@ ERF::MakeNewLevelFromCoarse (int lev, Real time, const BoxArray& ba, advflux_reg[lev] = new YAFluxRegister(ba, grids[lev-1], dm, dmap[lev-1], geom[lev], geom[lev-1], - ref_ratio[lev-1], lev, NVAR); + ref_ratio[lev-1], lev, ncomp); } // ***************************************************************************************************** @@ -269,11 +273,13 @@ ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapp Vector temp_lev_new(Vars::NumTypes); Vector temp_lev_old(Vars::NumTypes); + int ncomp_cons = vars_new[lev][Vars::cons].nComp(); + int ngrow_state = ComputeGhostCells(solverChoice.advChoice, solverChoice.use_NumDiff) + 1; int ngrow_vels = ComputeGhostCells(solverChoice.advChoice, solverChoice.use_NumDiff); - temp_lev_new[Vars::cons].define(ba, dm, Cons::NumVars, ngrow_state); - temp_lev_old[Vars::cons].define(ba, dm, Cons::NumVars, ngrow_state); + temp_lev_new[Vars::cons].define(ba, dm, ncomp_cons, ngrow_state); + temp_lev_old[Vars::cons].define(ba, dm, ncomp_cons, ngrow_state); temp_lev_new[Vars::xvel].define(convert(ba, IntVect(1,0,0)), dm, 1, ngrow_vels); temp_lev_old[Vars::xvel].define(convert(ba, IntVect(1,0,0)), dm, 1, ngrow_vels); @@ -292,7 +298,7 @@ ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapp advflux_reg[lev] = new YAFluxRegister(ba, grids[lev-1], dm, dmap[lev-1], geom[lev], geom[lev-1], - ref_ratio[lev-1], lev, NVAR); + ref_ratio[lev-1], lev, ncomp_cons); } //******************************************************************************************** @@ -315,10 +321,10 @@ ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapp // ******************************************************************************************** // Copy from new into old just in case // ******************************************************************************************** - MultiFab::Copy(temp_lev_old[Vars::cons],temp_lev_new[Vars::cons],0,0,NVAR,ngrow_state); - MultiFab::Copy(temp_lev_old[Vars::xvel],temp_lev_new[Vars::xvel],0,0, 1,ngrow_vels); - MultiFab::Copy(temp_lev_old[Vars::yvel],temp_lev_new[Vars::yvel],0,0, 1,ngrow_vels); - MultiFab::Copy(temp_lev_old[Vars::zvel],temp_lev_new[Vars::zvel],0,0, 1,IntVect(ngrow_vels,ngrow_vels,0)); + MultiFab::Copy(temp_lev_old[Vars::cons],temp_lev_new[Vars::cons],0,0,ncomp_cons,ngrow_state); + MultiFab::Copy(temp_lev_old[Vars::xvel],temp_lev_new[Vars::xvel],0,0, 1,ngrow_vels); + MultiFab::Copy(temp_lev_old[Vars::yvel],temp_lev_new[Vars::yvel],0,0, 1,ngrow_vels); + MultiFab::Copy(temp_lev_old[Vars::zvel],temp_lev_new[Vars::zvel],0,0, 1,IntVect(ngrow_vels,ngrow_vels,0)); // ******************************************************************************************** // Now swap the pointers @@ -464,9 +470,11 @@ ERF::initialize_integrator (int lev, MultiFab& cons_mf, MultiFab& vel_mf) const BoxArray& ba(cons_mf.boxArray()); const DistributionMapping& dm(cons_mf.DistributionMap()); + int ncomp_cons = cons_mf.nComp(); + // Initialize the integrator memory amrex::Vector int_state; // integration state data structure example - int_state.push_back(MultiFab(cons_mf, amrex::make_alias, 0, Cons::NumVars)); // cons + int_state.push_back(MultiFab(cons_mf, amrex::make_alias, 0, ncomp_cons)); // cons int_state.push_back(MultiFab(convert(ba,IntVect(1,0,0)), dm, 1, vel_mf.nGrow())); // xmom int_state.push_back(MultiFab(convert(ba,IntVect(0,1,0)), dm, 1, vel_mf.nGrow())); // ymom int_state.push_back(MultiFab(convert(ba,IntVect(0,0,1)), dm, 1, vel_mf.nGrow())); // zmom @@ -474,6 +482,7 @@ ERF::initialize_integrator (int lev, MultiFab& cons_mf, MultiFab& vel_mf) mri_integrator_mem[lev] = std::make_unique > >(int_state); mri_integrator_mem[lev]->setNoSubstepping(solverChoice.no_substepping); mri_integrator_mem[lev]->setIncompressible(solverChoice.incompressible); + mri_integrator_mem[lev]->setNcompCons(ncomp_cons); mri_integrator_mem[lev]->setForceFirstStageSingleSubstep(solverChoice.force_stage1_single_substep); } diff --git a/Source/IO/Checkpoint.cpp b/Source/IO/Checkpoint.cpp index effa60ba1..b1b587bc7 100644 --- a/Source/IO/Checkpoint.cpp +++ b/Source/IO/Checkpoint.cpp @@ -44,8 +44,10 @@ ERF::WriteCheckpointFile () const // ---- ParallelDescriptor::IOProcessor() creates the directories amrex::PreBuildDirectorHierarchy(checkpointname, "Level_", nlevels, true); + int ncomp_cons = vars_new[0][Vars::cons].nComp(); + // write Header file - if (ParallelDescriptor::IOProcessor()) { + if (ParallelDescriptor::IOProcessor()) { std::string HeaderFileName(checkpointname + "/Header"); VisMF::IO_Buffer io_buffer(VisMF::IO_Buffer_Size); @@ -70,7 +72,7 @@ ERF::WriteCheckpointFile () const // for each variable we store // conservative, cell-centered vars - HeaderFile << Cons::NumVars << "\n"; + HeaderFile << ncomp_cons << "\n"; // x-velocity on faces HeaderFile << 1 << "\n"; @@ -110,8 +112,8 @@ ERF::WriteCheckpointFile () const // Here we make copies of the MultiFab with no ghost cells for (int lev = 0; lev <= finest_level; ++lev) { - MultiFab cons(grids[lev],dmap[lev],Cons::NumVars,0); - MultiFab::Copy(cons,vars_new[lev][Vars::cons],0,0,NVAR,0); + MultiFab cons(grids[lev],dmap[lev],ncomp_cons,0); + MultiFab::Copy(cons,vars_new[lev][Vars::cons],0,0,ncomp_cons,0); VisMF::Write(cons, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "Cell")); MultiFab xvel(convert(grids[lev],IntVect(1,0,0)),dmap[lev],1,0); @@ -222,6 +224,8 @@ ERF::ReadCheckpointFile () { amrex::Print() << "Restart from checkpoint " << restart_chkfile << "\n"; + int ncomp_cons = vars_new[0][Vars::cons].nComp(); + // Header std::string File(restart_chkfile + "/Header"); @@ -251,10 +255,10 @@ ERF::ReadCheckpointFile () GotoNextLine(is); #if 0 if (solverChoice.moisture_type != MoistureType::None) { - AMREX_ASSERT(chk_ncomp == Cons::NumVars); + AMREX_ASSERT(chk_ncomp == ncomp_cons); } else { // This assumes that we are carrying RhoQ1 and RhoQ2 but not actually using them - AMREX_ASSERT(chk_ncomp == Cons::NumVars-2); + AMREX_ASSERT(chk_ncomp == ncomp_cons-2); } #endif @@ -319,14 +323,14 @@ ERF::ReadCheckpointFile () // read in the MultiFab data for (int lev = 0; lev <= finest_level; ++lev) { - MultiFab cons(grids[lev],dmap[lev],Cons::NumVars,0); + MultiFab cons(grids[lev],dmap[lev],ncomp_cons,0); VisMF::Read(cons, amrex::MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "Cell")); if (solverChoice.moisture_type != MoistureType::None) { - MultiFab::Copy(vars_new[lev][Vars::cons],cons,0,0,Cons::NumVars,0); + MultiFab::Copy(vars_new[lev][Vars::cons],cons,0,0,ncomp_cons,0); } else { // This assumes that we are carrying RhoQ1 and RhoQ2 but not actually using them - MultiFab::Copy(vars_new[lev][Vars::cons],cons,0,0,Cons::NumVars-2,0); + MultiFab::Copy(vars_new[lev][Vars::cons],cons,0,0,ncomp_cons-2,0); // We zero them here so they won't cause problems while we're still debugging vars_new[lev][Vars::cons].setVal(0.0, RhoQ1_comp, 2); diff --git a/Source/IO/ERF_ReadBndryPlanes.H b/Source/IO/ERF_ReadBndryPlanes.H index f2da53fa1..84a14139e 100644 --- a/Source/IO/ERF_ReadBndryPlanes.H +++ b/Source/IO/ERF_ReadBndryPlanes.H @@ -26,10 +26,10 @@ public: void read_time_file(); void read_input_files(amrex::Real time, amrex::Real dt, - amrex::Array,AMREX_SPACEDIM+NVAR> m_bc_extdir_vals); + amrex::Array, AMREX_SPACEDIM+NVAR_max> m_bc_extdir_vals); void read_file(int idx, amrex::Vector>& data_to_fill, - amrex::Array,AMREX_SPACEDIM+NVAR> m_bc_extdir_vals); + amrex::Array, AMREX_SPACEDIM+NVAR_max> m_bc_extdir_vals); // Return the pointer to PlaneVectors at time "time" amrex::Vector>& interp_in_time(const amrex::Real& time); @@ -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 f495e2ba6..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; } @@ -273,7 +265,7 @@ void ReadBndryPlanes::read_time_file() * @param m_bc_extdir_vals Container storing the external dirichlet boundary conditions we are reading from the input files */ void ReadBndryPlanes::read_input_files(Real time, Real dt, - Array,AMREX_SPACEDIM+NVAR> m_bc_extdir_vals) + Array,AMREX_SPACEDIM+NVAR_max> m_bc_extdir_vals) { BL_PROFILE("ERF::ReadBndryPlanes::read_input_files"); @@ -346,7 +338,7 @@ void ReadBndryPlanes::read_input_files(Real time, Real dt, * @param m_bc_extdir_vals Container storing the external dirichlet boundary conditions we are reading from the input files */ void ReadBndryPlanes::read_file(const int idx, Vector>& data_to_fill, - Array,AMREX_SPACEDIM+NVAR> m_bc_extdir_vals) + Array,AMREX_SPACEDIM+NVAR_max> m_bc_extdir_vals) { const int t_step = m_in_timesteps[idx]; const std::string chkname1 = m_filename + Concatenate("/bndry_output", t_step); @@ -358,8 +350,7 @@ void ReadBndryPlanes::read_file(const int idx, Vector, - AMREX_SPACEDIM+NVAR> l_bc_extdir_vals_d; + GpuArray, AMREX_SPACEDIM+NVAR_max> l_bc_extdir_vals_d; for (int i = 0; i < BCVars::NumTypes; i++) { @@ -413,10 +404,6 @@ void ReadBndryPlanes::read_file(const int idx, Vector // Get state variables in the same order as we define them, // since they may be in any order in the input list Vector tmp_plot_names; - for (int i = 0; i < Cons::NumVars; ++i) { + + int ncomp_cons = (solverChoice.moisture_type == MoistureType::None) ? NVAR_max-2 : NVAR_max; + + for (int i = 0; i < ncomp_cons; ++i) { if ( containerHasElement(plot_var_names, cons_names[i]) ) { tmp_plot_names.push_back(cons_names[i]); } @@ -101,6 +104,8 @@ ERF::WritePlotFile (int which, Vector plot_var_names) const Vector varnames = PlotFileVarNames(plot_var_names); const int ncomp_mf = varnames.size(); + int ncomp_cons = vars_new[0][Vars::cons].nComp(); + if (ncomp_mf == 0) return; @@ -128,12 +133,13 @@ ERF::WritePlotFile (int which, Vector plot_var_names) } } - for (int lev = 0; lev <= finest_level; ++lev) { + for (int lev = 0; lev <= finest_level; ++lev) + { int mf_comp = 0; // First, copy any of the conserved state variables into the output plotfile - AMREX_ALWAYS_ASSERT(cons_names.size() == Cons::NumVars); - for (int i = 0; i < Cons::NumVars; ++i) { + AMREX_ALWAYS_ASSERT(cons_names.size() >= ncomp_cons); + for (int i = 0; i < ncomp_cons; ++i) { if (containerHasElement(plot_var_names, cons_names[i])) { MultiFab::Copy(mf[lev],vars_new[lev][Vars::cons],i,mf_comp,1,0); mf_comp++; 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/IndexDefines.H b/Source/IndexDefines.H index 2bccf05dd..09b299a6c 100644 --- a/Source/IndexDefines.H +++ b/Source/IndexDefines.H @@ -17,7 +17,10 @@ #define RhoQ1_comp 5 #define RhoQ2_comp 6 -#define NVAR 7 +// This is the number of components if using moisture +// We use this to allocate the 1d arrays of boundary condition types, +// but not to allocate actual solution data +#define NVAR_max 7 // Cell-centered primitive variables #define PrimTheta_comp (RhoTheta_comp -1) @@ -27,8 +30,7 @@ #define PrimQ1_comp (RhoQ1_comp-1) #define PrimQ2_comp (RhoQ2_comp-1) -#define NUM_PRIM (NVAR-1) - +// NOTE: we still use this indexing even if no moisture namespace BCVars { enum { cons_bc = 0, @@ -39,9 +41,9 @@ namespace BCVars { RhoScalar_bc_comp, RhoQ1_bc_comp, RhoQ2_bc_comp, - xvel_bc = NVAR, - yvel_bc, - zvel_bc, + xvel_bc = NVAR_max, + yvel_bc = NVAR_max+1, + zvel_bc = NVAR_max+2, NumTypes }; } @@ -83,33 +85,6 @@ namespace Vars { }; } -namespace Cons { - enum { - Rho = 0, - RhoTheta, - RhoKE, - RhoQKE, - RhoScalar, - RhoQ1, - RhoQ2, - NumVars - }; -} - -#if 0 -namespace Prim { - enum { - Theta = 0, - KE, - QKE, - Scalar, - Q1, - Q2, - NumVars - }; -} -#endif - // We separate out horizontal and vertical turbulent diffusivities // These are the same for LES, but different for PBL models namespace EddyDiff { diff --git a/Source/Initialization/ERF_init_bcs.cpp b/Source/Initialization/ERF_init_bcs.cpp index 7ae16d317..4f2413ef5 100644 --- a/Source/Initialization/ERF_init_bcs.cpp +++ b/Source/Initialization/ERF_init_bcs.cpp @@ -48,7 +48,6 @@ void ERF::init_bcs () m_bc_neumann_vals[BCVars::yvel_bc][ori] = 0.0; m_bc_neumann_vals[BCVars::zvel_bc][ori] = 0.0; - ParmParse pp(bcid); std::string bc_type_in = "null"; pp.query("type", bc_type_in); @@ -111,37 +110,23 @@ void ERF::init_bcs () m_bc_extdir_vals[BCVars::RhoScalar_bc_comp][ori] = rho_in*scalar_in; } - Real qt_in = 0.; - if (input_bndry_planes && m_r2d->ingested_q1()) { - m_bc_extdir_vals[BCVars::RhoQ1_bc_comp][ori] = 0.; - } else { - if (pp.query("qt", qt_in)) - m_bc_extdir_vals[BCVars::RhoQ1_bc_comp][ori] = rho_in*qt_in; - } - Real qp_in = 0.; - if (input_bndry_planes && m_r2d->ingested_q2()) { - m_bc_extdir_vals[BCVars::RhoQ2_bc_comp][ori] = 0.; - } else { - if (pp.query("qp", qp_in)) - m_bc_extdir_vals[BCVars::RhoQ2_bc_comp][ori] = rho_in*qp_in; + if (solverChoice.moisture_type != MoistureType::None) { + Real qt_in = 0.; + if (input_bndry_planes && m_r2d->ingested_q1()) { + m_bc_extdir_vals[BCVars::RhoQ1_bc_comp][ori] = 0.; + } else { + if (pp.query("qt", qt_in)) + m_bc_extdir_vals[BCVars::RhoQ1_bc_comp][ori] = rho_in*qt_in; + } + Real qp_in = 0.; + if (input_bndry_planes && m_r2d->ingested_q2()) { + m_bc_extdir_vals[BCVars::RhoQ2_bc_comp][ori] = 0.; + } else { + if (pp.query("qp", qp_in)) + m_bc_extdir_vals[BCVars::RhoQ2_bc_comp][ori] = rho_in*qp_in; + } } -#if defined(ERF_USE_WARM_NO_PRECIP) - Real qv_in = 0.; - if (input_bndry_planes && m_r2d->ingested_qv()) { - m_bc_extdir_vals[BCVars::RhoQv_bc_comp][ori] = 0.; - } else { - if (pp.query("qv", qv_in)) - m_bc_extdir_vals[BCVars::RhoQv_bc_comp][ori] = rho_in*qv_in; - } - Real qc_in = 0.; - if (input_bndry_planes && m_r2d->ingested_qc()) { - m_bc_extdir_vals[BCVars::RhoQc_bc_comp][ori] = 0.; - } else { - if (pp.query("qc", qc_in)) - m_bc_extdir_vals[BCVars::RhoQc_bc_comp][ori] = rho_in*qc_in; - } -#endif Real KE_in = 0.; if (input_bndry_planes && m_r2d->ingested_KE()) { m_bc_extdir_vals[BCVars::RhoKE_bc_comp][ori] = 0.; @@ -267,8 +252,8 @@ void ERF::init_bcs () // // ***************************************************************************** { - domain_bcs_type.resize(AMREX_SPACEDIM+NVAR); - domain_bcs_type_d.resize(AMREX_SPACEDIM+NVAR); + domain_bcs_type.resize(AMREX_SPACEDIM+NVAR_max); + domain_bcs_type_d.resize(AMREX_SPACEDIM+NVAR_max); for (OrientationIter oit; oit; ++oit) { Orientation ori = oit(); @@ -375,34 +360,34 @@ void ERF::init_bcs () if ( bct == ERF_BC::symmetry ) { if (side == Orientation::low) { - for (int i = 0; i < NVAR; i++) + for (int i = 0; i < NVAR_max; i++) domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::reflect_even); } else { - for (int i = 0; i < NVAR; i++) + for (int i = 0; i < NVAR_max; i++) domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::reflect_even); } } else if ( bct == ERF_BC::outflow ) { if (side == Orientation::low) { - for (int i = 0; i < NVAR; i++) + for (int i = 0; i < NVAR_max; i++) domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::foextrap); } else { - for (int i = 0; i < NVAR; i++) + for (int i = 0; i < NVAR_max; i++) domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::foextrap); } } else if ( bct == ERF_BC::no_slip_wall) { if (side == Orientation::low) { - for (int i = 0; i < NVAR; i++) + for (int i = 0; i < NVAR_max; i++) domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::foextrap); if (m_bc_extdir_vals[BCVars::RhoTheta_bc_comp][ori] > 0.) domain_bcs_type[BCVars::RhoTheta_bc_comp].setLo(dir, ERFBCType::ext_dir); if (std::abs(m_bc_neumann_vals[BCVars::RhoTheta_bc_comp][ori]) > 0.) domain_bcs_type[BCVars::RhoTheta_bc_comp].setLo(dir, ERFBCType::neumann); } else { - for (int i = 0; i < NVAR; i++) + for (int i = 0; i < NVAR_max; i++) domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::foextrap); if (m_bc_extdir_vals[BCVars::RhoTheta_bc_comp][ori] > 0.) domain_bcs_type[BCVars::RhoTheta_bc_comp].setHi(dir, ERFBCType::ext_dir); @@ -413,7 +398,7 @@ void ERF::init_bcs () else if (bct == ERF_BC::slip_wall) { if (side == Orientation::low) { - for (int i = 0; i < NVAR; i++) + for (int i = 0; i < NVAR_max; i++) domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::foextrap); if (m_bc_extdir_vals[BCVars::RhoTheta_bc_comp][ori] > 0.) domain_bcs_type[BCVars::RhoTheta_bc_comp].setLo(dir, ERFBCType::ext_dir); @@ -422,7 +407,7 @@ void ERF::init_bcs () if (std::abs(m_bc_neumann_vals[BCVars::Rho_bc_comp][ori]) > 0.) domain_bcs_type[BCVars::Rho_bc_comp].setLo(dir, ERFBCType::neumann); } else { - for (int i = 0; i < NVAR; i++) + for (int i = 0; i < NVAR_max; i++) domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::foextrap); if (m_bc_extdir_vals[BCVars::RhoTheta_bc_comp][ori] > 0.) domain_bcs_type[BCVars::RhoTheta_bc_comp].setHi(dir, ERFBCType::ext_dir); @@ -435,7 +420,7 @@ void ERF::init_bcs () else if (bct == ERF_BC::inflow) { if (side == Orientation::low) { - for (int i = 0; i < NVAR; i++) { + for (int i = 0; i < NVAR_max; i++) { domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::ext_dir); if (input_bndry_planes && dir < 2 && ( ( (BCVars::cons_bc+i == BCVars::Rho_bc_comp) && m_r2d->ingested_density()) || @@ -449,7 +434,7 @@ void ERF::init_bcs () } } } else { - for (int i = 0; i < NVAR; i++) { + for (int i = 0; i < NVAR_max; i++) { domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::ext_dir); if (input_bndry_planes && dir < 2 && ( ( (BCVars::cons_bc+i == BCVars::Rho_bc_comp) && m_r2d->ingested_density()) || @@ -468,17 +453,17 @@ void ERF::init_bcs () else if (bct == ERF_BC::periodic) { if (side == Orientation::low) { - for (int i = 0; i < NVAR; i++) + for (int i = 0; i < NVAR_max; i++) domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::int_dir); } else { - for (int i = 0; i < NVAR; i++) + for (int i = 0; i < NVAR_max; i++) domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::int_dir); } } else if ( bct == ERF_BC::MOST ) { AMREX_ALWAYS_ASSERT(dir == 2 && side == Orientation::low); - for (int i = 0; i < NVAR; i++) { + for (int i = 0; i < NVAR_max; i++) { domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::foextrap); } } @@ -488,11 +473,11 @@ void ERF::init_bcs () #ifdef AMREX_USE_GPU Gpu::htod_memcpy (domain_bcs_type_d.data(), domain_bcs_type.data(), - sizeof(amrex::BCRec)*(NVAR+AMREX_SPACEDIM)); + sizeof(amrex::BCRec)*(NVAR_max+AMREX_SPACEDIM)); #else std::memcpy (domain_bcs_type_d.data(), domain_bcs_type.data(), - sizeof(amrex::BCRec)*(NVAR+AMREX_SPACEDIM)); + sizeof(amrex::BCRec)*(NVAR_max+AMREX_SPACEDIM)); #endif } diff --git a/Source/TimeIntegration/ERF_Advance.cpp b/Source/TimeIntegration/ERF_Advance.cpp index eba329f78..362464d16 100644 --- a/Source/TimeIntegration/ERF_Advance.cpp +++ b/Source/TimeIntegration/ERF_Advance.cpp @@ -37,8 +37,8 @@ ERF::Advance (int lev, Real time, Real dt_lev, int /*iteration*/, int /*ncycle*/ if (phys_bc_type[Orientation(Direction::z,Orientation::low)] == ERF_BC::MOST) { if (m_most) { amrex::IntVect ng = S_old.nGrowVect(); ng[2]=0; - MultiFab::Copy( *Theta_prim[lev], S_old, Cons::RhoTheta, 0, 1, ng); - MultiFab::Divide(*Theta_prim[lev], S_old, Cons::Rho , 0, 1, ng); + MultiFab::Copy( *Theta_prim[lev], S_old, RhoTheta_comp, 0, 1, ng); + MultiFab::Divide(*Theta_prim[lev], S_old, Rho_comp , 0, 1, ng); // NOTE: std::swap above causes the field ptrs to be out of date. // Reassign the field ptrs for MAC avg computation. m_most->update_mac_ptrs(lev, vars_old, Theta_prim); @@ -54,9 +54,10 @@ 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 defined(ERF_USE_MOISTURE) - FillPatchMoistVars(lev, *(qmoist[lev])); -#endif + + if (solverChoice.moisture_type != MoistureType::None) { + FillPatchMoistVars(lev, *qmoist[lev]); + } MultiFab* S_crse; MultiFab rU_crse, rV_crse, rW_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,10 +124,8 @@ ERF::Advance (int lev, Real time, Real dt_lev, int /*iteration*/, int /*ncycle*/ source, buoyancy, Geom(lev), dt_lev, time, &ifr); -#if defined(ERF_USE_MOISTURE) - // Update the microphysics + // Update the microphysics (moisture) advance_microphysics(lev, S_new, dt_lev); -#endif #if defined(ERF_USE_RRTMGP) // Update the radiation diff --git a/Source/TimeIntegration/ERF_MRI.H b/Source/TimeIntegration/ERF_MRI.H index 07d2d5cf7..65034e309 100644 --- a/Source/TimeIntegration/ERF_MRI.H +++ b/Source/TimeIntegration/ERF_MRI.H @@ -41,6 +41,11 @@ private: */ int incompressible; + /** + * \brief How many components in the cell-centered MultiFab + */ + int ncomp_cons; + /** * \brief Do we follow the recommendation to only perform a single substep in the first RK stage */ @@ -63,7 +68,7 @@ private: void initialize_data (const T& S_data) { // TODO: We can optimize memory by making the cell-centered part of S_sum, S_scratch - // have only 2 components, not Cons::NumVars components + // have only 2 components, not ncomp_cons components const bool include_ghost = true; amrex::IntegratorOps::CreateLike(T_store, S_data, include_ghost); S_sum = T_store[0].get(); @@ -105,6 +110,11 @@ public: // Delete the copy assignment operator MRISplitIntegrator& operator=(const MRISplitIntegrator& other) = delete; + void setNcompCons(int _ncomp_cons) + { + ncomp_cons = _ncomp_cons; + } + void setIncompressible(int _incompressible) { incompressible = _incompressible; @@ -210,7 +220,7 @@ public: /**********************************************/ /* RK3 Integration with Acoustic Sub-stepping */ /**********************************************/ - Vector num_vars = {Cons::NumVars, 1, 1, 1}; + Vector num_vars = {ncomp_cons, 1, 1, 1}; for (int i(0); i new diff --git a/Source/TimeIntegration/ERF_advance_dycore.cpp b/Source/TimeIntegration/ERF_advance_dycore.cpp index 0a1b81832..3ed147cc6 100644 --- a/Source/TimeIntegration/ERF_advance_dycore.cpp +++ b/Source/TimeIntegration/ERF_advance_dycore.cpp @@ -67,14 +67,6 @@ void ERF::advance_dycore(int level, MultiFab p_hse (base_state[level], make_alias, 1, 1); // p_0 is second component MultiFab pi_hse(base_state[level], make_alias, 2, 1); // pi_0 is second component -#if defined(ERF_USE_MOISTURE) - // TODO: Protect from grabbing non-existing data - int q_size = micro.Get_Qmoist_Size(); - MultiFab qvapor (*(qmoist[level]), make_alias, 0, 1); - MultiFab qcloud (*(qmoist[level]), make_alias, 1, 1); - MultiFab qice (*(qmoist[level]), make_alias, 2, 1); -#endif - // These pointers are used in the MRI utility functions MultiFab* r0 = &r_hse; MultiFab* p0 = &p_hse; @@ -97,7 +89,9 @@ void ERF::advance_dycore(int level, const BoxArray& ba_z = zvel_old.boxArray(); const DistributionMapping& dm = cons_old.DistributionMap(); - MultiFab S_prim (ba , dm, NUM_PRIM, cons_old.nGrowVect()); + int num_prim = cons_old.nComp() - 1; + + MultiFab S_prim (ba , dm, num_prim, cons_old.nGrowVect()); MultiFab pi_stage (ba , dm, 1, cons_old.nGrowVect()); MultiFab fast_coeffs(ba_z, dm, 5, 0); MultiFab* eddyDiffs = eddyDiffs_lev[level].get(); @@ -256,33 +250,44 @@ void ERF::advance_dycore(int level, cons_to_prim(state_old[IntVar::cons], state_old[IntVar::cons].nGrow()); } // profile -#ifdef ERF_USE_MOISTURE - 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); + int q_size = micro.Get_Qmoist_Size(); + int ncell = fine_geom.Domain().length(2); + + 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 4ab7fde22..93373efcf 100644 --- a/Source/TimeIntegration/ERF_advance_microphysics.cpp +++ b/Source/TimeIntegration/ERF_advance_microphysics.cpp @@ -2,16 +2,16 @@ using namespace amrex; -#if defined(ERF_USE_MOISTURE) void ERF::advance_microphysics (int lev, MultiFab& cons, const Real& dt_advance) { - micro.Init(cons, *(qmoist[lev]), - grids[lev], - Geom(lev), - dt_advance); - micro.Advance(); - micro.Update(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])); + } } -#endif diff --git a/Source/TimeIntegration/ERF_fast_rhs_MT.cpp b/Source/TimeIntegration/ERF_fast_rhs_MT.cpp index 28675be13..c1c689652 100644 --- a/Source/TimeIntegration/ERF_fast_rhs_MT.cpp +++ b/Source/TimeIntegration/ERF_fast_rhs_MT.cpp @@ -79,6 +79,7 @@ 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()"); @@ -258,15 +259,12 @@ 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 defined(ERF_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); -#elif defined(ERF_USE_WARM_NO_PRECIP) - Real q = 0.5 * ( prim(i,j,k,PrimQv_comp) + prim(i-1,j,k,PrimQv_comp) - +prim(i,j,k,PrimQc_comp) + prim(i-1,j,k,PrimQc_comp) ); - gpx /= (1.0 + q); -#endif + 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); + } + Real pi_c = 0.5 * (pi_stg_ca(i-1,j,k,0) + pi_stg_ca(i ,j,k,0)); Real fast_rhs_rho_u = -Gamma * R_d * pi_c * gpx; @@ -288,15 +286,12 @@ 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 defined(ERF_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); -#elif defined(ERF_USE_WARM_NO_PRECIP) - Real q = 0.5 * ( prim(i,j,k,PrimQv_comp) + prim(i,j-1,k,PrimQv_comp) - +prim(i,j,k,PrimQc_comp) + prim(i,j-1,k,PrimQc_comp) ); - gpy /= (1.0 + q); -#endif + 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); + } + Real pi_c = 0.5 * (pi_stg_ca(i,j-1,k,0) + pi_stg_ca(i,j ,k,0)); Real fast_rhs_rho_v = -Gamma * R_d * pi_c * gpy; @@ -401,17 +396,12 @@ 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 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-1,PrimQv_comp) - +prim(i,j,k,PrimQc_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_fast_rhs_N.cpp b/Source/TimeIntegration/ERF_fast_rhs_N.cpp index 729bf0d33..45f47fc80 100644 --- a/Source/TimeIntegration/ERF_fast_rhs_N.cpp +++ b/Source/TimeIntegration/ERF_fast_rhs_N.cpp @@ -58,6 +58,7 @@ 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()"); @@ -185,9 +186,7 @@ void erf_fast_rhs_N (int step, int nrk, const Array4 & stage_xmom = S_stage_data[IntVar::xmom].const_array(mfi); const Array4 & stage_ymom = S_stage_data[IntVar::ymom].const_array(mfi); -#if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP) const Array4 & prim = S_stage_prim.const_array(mfi); -#endif const Array4& slow_rhs_rho_u = S_slow_rhs[IntVar::xmom].const_array(mfi); const Array4& slow_rhs_rho_v = S_slow_rhs[IntVar::ymom].const_array(mfi); @@ -222,15 +221,12 @@ 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 defined(ERF_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); -#elif defined(ERF_USE_WARM_NO_PRECIP) - Real q = 0.5 * ( prim(i,j,k,PrimQv_comp) + prim(i-1,j,k,PrimQv_comp) - +prim(i,j,k,PrimQc_comp) + prim(i-1,j,k,PrimQc_comp) ); - gpx /= (1.0 + q); -#endif + 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); + } + Real pi_c = 0.5 * (pi_stage_ca(i-1,j,k,0) + pi_stage_ca(i,j,k,0)); Real fast_rhs_rho_u = -Gamma * R_d * pi_c * gpx; @@ -248,15 +244,12 @@ 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 defined(ERF_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); -#elif defined(ERF_USE_WARM_NO_PRECIP) - Real q = 0.5 * ( prim(i,j,k,PrimQv_comp) + prim(i,j-1,k,PrimQv_comp) - +prim(i,j,k,PrimQc_comp) + prim(i,j-1,k,PrimQc_comp) ); - gpy /= (1.0 + q); -#endif + 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); + } + Real pi_c = 0.5 * (pi_stage_ca(i,j-1,k,0) + pi_stage_ca(i,j,k,0)); Real fast_rhs_rho_v = -Gamma * R_d * pi_c * gpy; @@ -382,17 +375,12 @@ 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 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-1,PrimQv_comp) - +prim(i,j,k,PrimQc_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_fast_rhs_T.cpp b/Source/TimeIntegration/ERF_fast_rhs_T.cpp index b1ac58d6b..4efcea56d 100644 --- a/Source/TimeIntegration/ERF_fast_rhs_T.cpp +++ b/Source/TimeIntegration/ERF_fast_rhs_T.cpp @@ -65,6 +65,7 @@ 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()"); @@ -201,9 +202,7 @@ void erf_fast_rhs_T (int step, int nrk, const Array4 & stage_xmom = S_stage_data[IntVar::xmom].const_array(mfi); const Array4 & stage_ymom = S_stage_data[IntVar::ymom].const_array(mfi); -#if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP) const Array4 & prim = S_stage_prim.const_array(mfi); -#endif const Array4& old_drho_u = Delta_rho_u.array(mfi); const Array4& old_drho_v = Delta_rho_v.array(mfi); @@ -257,15 +256,12 @@ 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 defined(ERF_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); -#elif defined(ERF_USE_WARM_NO_PRECIP) - Real q = 0.5 * ( prim(i,j,k,PrimQv_comp) + prim(i-1,j,k,PrimQv_comp) - +prim(i,j,k,PrimQc_comp) + prim(i-1,j,k,PrimQc_comp) ); - gpx /= (1.0 + q); -#endif + 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); + } + Real pi_c = 0.5 * (pi_stage_ca(i-1,j,k,0) + pi_stage_ca(i ,j,k,0)); Real fast_rhs_rho_u = -Gamma * R_d * pi_c * gpx; @@ -290,15 +286,12 @@ 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 defined(ERF_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); -#elif defined(ERF_USE_WARM_NO_PRECIP) - Real q = 0.5 * ( prim(i,j,k,PrimQv_comp) + prim(i,j-1,k,PrimQv_comp) - +prim(i,j,k,PrimQc_comp) + prim(i,j-1,k,PrimQc_comp) ); - gpy /= (1.0 + q); -#endif + 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); + } + Real pi_c = 0.5 * (pi_stage_ca(i,j-1,k,0) + pi_stage_ca(i,j ,k,0)); Real fast_rhs_rho_v = -Gamma * R_d * pi_c * gpy; @@ -468,17 +461,12 @@ 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 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-1,PrimQv_comp) - +prim(i,j,k,PrimQc_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_make_buoyancy.cpp b/Source/TimeIntegration/ERF_make_buoyancy.cpp index e0e8aec42..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,101 +52,95 @@ 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.moisture_type == MoistureType::None) { + if (solverChoice.buoyancy_type == 1) { + for ( MFIter mfi(buoyancy,TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + Box tbz = mfi.tilebox(); - 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 // ****************************************************************************************** -#if defined(ERF_USE_MOISTURE) - if (solverChoice.buoyancy_type == 1) { + if (solverChoice.moisture_type == MoistureType::FastEddy) { + + AMREX_ALWAYS_ASSERT(solverChoice.buoyancy_type == 1); for ( MFIter mfi(buoyancy,TilingIfNotGPU()); mfi.isValid(); ++mfi) { @@ -166,195 +156,193 @@ 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); - amrex::ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - Real rhop_hi = cell_data(i,j,k ,Rho_comp) * (1.0 + qv_data(i,j,k ) + qc_data(i,j,k ) - + qi_data(i,j,k )) + cell_data(i,j,k,RhoQ2_comp) - r0_arr(i,j,k ); - Real rhop_lo = cell_data(i,j,k-1,Rho_comp) * (1.0 + qv_data(i,j,k-1) + qc_data(i,j,k-1) - + qi_data(i,j,k-1)) + cell_data(i,j,k-1,RhoQ2_comp) - r0_arr(i,j,k-1); + Real rhop_hi = cell_data(i,j,k ,Rho_comp) + cell_data(i,j,k ,RhoQ1_comp) + + cell_data(i,j,k ,RhoQ2_comp) - r0_arr(i,j,k ); + Real rhop_lo = cell_data(i,j,k-1,Rho_comp) + cell_data(i,j,k-1,RhoQ1_comp) + + cell_data(i,j,k-1,RhoQ2_comp) - r0_arr(i,j,k-1); buoyancy_fab(i, j, k) = grav_gpu[2] * 0.5 * ( rhop_hi + rhop_lo ); }); } // mfi - } else { + } else if (solverChoice.moisture_type != MoistureType::None) { - PlaneAverage state_ave(&(S_data[IntVar::cons]), geom, solverChoice.ave_plane); - PlaneAverage prim_ave(&S_prim , geom, solverChoice.ave_plane); + if (solverChoice.buoyancy_type == 1) { - // Compute horizontal averages of all components of each field - state_ave.compute_averages(ZDir(), state_ave.field()); - prim_ave.compute_averages(ZDir(), prim_ave.field()); + for ( MFIter mfi(buoyancy,TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + Box tbz = mfi.tilebox(); - int ncell = state_ave.ncell_line(); + // 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); - Gpu::HostVector rho_h(ncell), theta_h(ncell), qp_h(ncell); - Gpu::DeviceVector rho_d(ncell), theta_d(ncell); + const Array4 & cell_data = S_data[IntVar::cons].array(mfi); + const Array4< Real> & buoyancy_fab = buoyancy.array(mfi); - state_ave.line_average(Rho_comp, rho_h); - Gpu::copyAsync(Gpu::hostToDevice, rho_h.begin(), rho_h.end(), rho_d.begin()); + // Base state density + const Array4& r0_arr = r0->const_array(mfi); - prim_ave.line_average(PrimTheta_comp, theta_h); - Gpu::copyAsync(Gpu::hostToDevice, theta_h.begin(), theta_h.end(), theta_d.begin()); + const Array4 & qv_data = qmoist->array(mfi,0); + const Array4 & qc_data = qmoist->array(mfi,1); + const Array4 & qi_data = qmoist->array(mfi,2); - Real* rho_d_ptr = rho_d.data(); - Real* theta_d_ptr = theta_d.data(); + amrex::ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + Real rhop_hi = cell_data(i,j,k ,Rho_comp) * (1.0 + qv_data(i,j,k ) + qc_data(i,j,k ) + + qi_data(i,j,k )) + cell_data(i,j,k,RhoQ2_comp) - r0_arr(i,j,k ); + Real rhop_lo = cell_data(i,j,k-1,Rho_comp) * (1.0 + qv_data(i,j,k-1) + qc_data(i,j,k-1) + + qi_data(i,j,k-1)) + cell_data(i,j,k-1,RhoQ2_comp) - r0_arr(i,j,k-1); + buoyancy_fab(i, j, k) = grav_gpu[2] * 0.5 * ( rhop_hi + rhop_lo ); + }); + } // mfi - if (solverChoice.buoyancy_type == 2 || solverChoice.buoyancy_type == 4 ) { + } else { - prim_ave.line_average(PrimQ2_comp, qp_h); + PlaneAverage state_ave(&(S_data[IntVar::cons]), geom, solverChoice.ave_plane); + PlaneAverage prim_ave(&S_prim , geom, solverChoice.ave_plane); - Gpu::DeviceVector qp_d(ncell); + // Compute horizontal averages of all components of each field + state_ave.compute_averages(ZDir(), state_ave.field()); + prim_ave.compute_averages(ZDir(), prim_ave.field()); - // Copy data to device - Gpu::copyAsync(Gpu::hostToDevice, qp_h.begin(), qp_h.end(), qp_d.begin()); - Gpu::streamSynchronize(); + int ncell = state_ave.ncell_line(); - Real* qp_d_ptr = qp_d.data(); - Real* qv_d_ptr = qv_d.data(); - Real* qc_d_ptr = qc_d.data(); - Real* qi_d_ptr = qi_d.data(); + Gpu::HostVector rho_h(ncell), theta_h(ncell), qp_h(ncell); + Gpu::DeviceVector rho_d(ncell), theta_d(ncell); -#ifdef _OPENMP -#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) -#endif - for ( MFIter mfi(buoyancy,TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - Box tbz = mfi.tilebox(); + state_ave.line_average(Rho_comp, rho_h); + Gpu::copyAsync(Gpu::hostToDevice, rho_h.begin(), rho_h.end(), rho_d.begin()); - // 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); + prim_ave.line_average(PrimTheta_comp, theta_h); + Gpu::copyAsync(Gpu::hostToDevice, theta_h.begin(), theta_h.end(), theta_d.begin()); - const Array4< Real> & buoyancy_fab = buoyancy.array(mfi); + Real* rho_d_ptr = rho_d.data(); + Real* theta_d_ptr = theta_d.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); + if (solverChoice.buoyancy_type == 2 || solverChoice.buoyancy_type == 4 ) { - 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]); + prim_ave.line_average(PrimQ2_comp, qp_h); - 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)); + Gpu::DeviceVector qp_d(ncell); - Real qplus, qminus; + // Copy data to device + Gpu::copyAsync(Gpu::hostToDevice, qp_h.begin(), qp_h.end(), qp_d.begin()); + Gpu::streamSynchronize(); - if(solverChoice.buoyancy_type == 2){ - qplus = 0.61* ( qv_data(i,j,k)-qv_d_ptr[k]) - - (qc_data(i,j,k)-qc_d_ptr[k]+ - qi_data(i,j,k)-qi_d_ptr[k]+ - cell_prim(i,j,k,PrimQ2_comp)-qp_d_ptr[k]) - + (tempp3d-tempp1d)/tempp1d*(Real(1.0) + Real(0.61)*qv_d_ptr[k]-qc_d_ptr[k]-qi_d_ptr[k]-qp_d_ptr[k]); + Real* qp_d_ptr = qp_d.data(); + Real* qv_d_ptr = qv_d.data(); + Real* qc_d_ptr = qc_d.data(); + Real* qi_d_ptr = qi_d.data(); - qminus = 0.61 *( qv_data(i,j,k-1)-qv_d_ptr[k-1]) - - (qc_data(i,j,k-1)-qc_d_ptr[k-1]+ - qi_data(i,j,k-1)-qi_d_ptr[k-1]+ - cell_prim(i,j,k-1,PrimQ2_comp)-qp_d_ptr[k-1]) - + (tempm3d-tempm1d)/tempm1d*(Real(1.0) + Real(0.61)*qv_d_ptr[k-1]-qi_d_ptr[k-1]-qc_d_ptr[k-1]-qp_d_ptr[k-1]); - } +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(buoyancy,TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + Box tbz = mfi.tilebox(); - if(solverChoice.buoyancy_type == 4){ - qplus = 0.61* ( qv_data(i,j,k)-qv_d_ptr[k]) - - (qc_data(i,j,k)-qc_d_ptr[k]+ - qi_data(i,j,k)-qi_d_ptr[k]+ - cell_prim(i,j,k,PrimQ2_comp)-qp_d_ptr[k]) - + (cell_data(i,j,k ,RhoTheta_comp)/cell_data(i,j,k ,Rho_comp) - theta_d_ptr[k ])/theta_d_ptr[k ]; + // 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); - qminus = 0.61 *( qv_data(i,j,k-1)-qv_d_ptr[k-1]) - - (qc_data(i,j,k-1)-qc_d_ptr[k-1]+ - qi_data(i,j,k-1)-qi_d_ptr[k-1]+ - cell_prim(i,j,k-1,PrimQ2_comp)-qp_d_ptr[k-1]) - + (cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp) - theta_d_ptr[k-1])/theta_d_ptr[k-1]; - } + const Array4< Real> & buoyancy_fab = buoyancy.array(mfi); + const Array4 & cell_data = S_data[IntVar::cons].array(mfi); + const Array4 & cell_prim = S_prim.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); - Real qavg = Real(0.5) * (qplus + qminus); - Real r0avg = Real(0.5) * (rho_d_ptr[k] + rho_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]); - buoyancy_fab(i, j, k) = -qavg * r0avg * grav_gpu[2]; - }); - } // mfi + 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)); - } else if (solverChoice.buoyancy_type == 3) { + Real qplus, qminus; -#ifdef _OPENMP -#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) -#endif - for ( MFIter mfi(buoyancy,TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - Box tbz = mfi.tilebox(); + if (solverChoice.buoyancy_type == 2) { + qplus = 0.61* ( qv_data(i,j,k)-qv_d_ptr[k]) - + (qc_data(i,j,k)-qc_d_ptr[k]+ + qi_data(i,j,k)-qi_d_ptr[k]+ + cell_prim(i,j,k,PrimQ2_comp)-qp_d_ptr[k]) + + (tempp3d-tempp1d)/tempp1d*(Real(1.0) + Real(0.61)*qv_d_ptr[k]-qc_d_ptr[k]-qi_d_ptr[k]-qp_d_ptr[k]); - // 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); + qminus = 0.61 *( qv_data(i,j,k-1)-qv_d_ptr[k-1]) - + (qc_data(i,j,k-1)-qc_d_ptr[k-1]+ + qi_data(i,j,k-1)-qi_d_ptr[k-1]+ + cell_prim(i,j,k-1,PrimQ2_comp)-qp_d_ptr[k-1]) + + (tempm3d-tempm1d)/tempm1d*(Real(1.0) + Real(0.61)*qv_d_ptr[k-1]-qi_d_ptr[k-1]-qc_d_ptr[k-1]-qp_d_ptr[k-1]); - const Array4< Real> & buoyancy_fab = buoyancy.array(mfi); - - 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); + } else if (solverChoice.buoyancy_type == 4) { - 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]); + qplus = 0.61* ( qv_data(i,j,k)-qv_d_ptr[k]) - + (qc_data(i,j,k)-qc_d_ptr[k]+ + qi_data(i,j,k)-qi_d_ptr[k]+ + cell_prim(i,j,k,PrimQ2_comp)-qp_d_ptr[k]) + + (cell_data(i,j,k ,RhoTheta_comp)/cell_data(i,j,k ,Rho_comp) - theta_d_ptr[k ])/theta_d_ptr[k ]; - 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)); + qminus = 0.61 *( qv_data(i,j,k-1)-qv_d_ptr[k-1]) - + (qc_data(i,j,k-1)-qc_d_ptr[k-1]+ + qi_data(i,j,k-1)-qi_d_ptr[k-1]+ + cell_prim(i,j,k-1,PrimQ2_comp)-qp_d_ptr[k-1]) + + (cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp) - theta_d_ptr[k-1])/theta_d_ptr[k-1]; + } - Real qplus = 0.61 * qv_data(i,j,k) - (qc_data(i,j,k)+ qi_data(i,j,k)+ cell_prim(i,j,k,PrimQ2_comp)) - + (tempp3d-tempp1d)/tempp1d; + Real qavg = Real(0.5) * (qplus + qminus); + Real r0avg = Real(0.5) * (rho_d_ptr[k] + rho_d_ptr[k-1]); - Real qminus = 0.61 *qv_data(i,j,k-1) - (qc_data(i,j,k-1)+ qi_data(i,j,k-1)+ cell_prim(i,j,k-1,PrimQ2_comp)) - + (tempm3d-tempm1d)/tempm1d; + buoyancy_fab(i, j, k) = -qavg * r0avg * grav_gpu[2]; + }); + } // mfi - Real qavg = Real(0.5) * (qplus + qminus); - Real r0avg = Real(0.5) * (rho_d_ptr[k] + rho_d_ptr[k-1]); + } else if (solverChoice.buoyancy_type == 3) { - buoyancy_fab(i, j, k) = -qavg * r0avg * grav_gpu[2]; - }); - } // mfi - } // buoyancy_type - } // not buoyancy_type == 1 +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) #endif + for ( MFIter mfi(buoyancy,TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + Box tbz = mfi.tilebox(); -#if defined(ERF_USE_WARM_NO_PRECIP) + // 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); - AMREX_ALWAYS_ASSERT(solverChoice.buoyancy_type == 1); + const Array4< Real> & buoyancy_fab = buoyancy.array(mfi); - for ( MFIter mfi(buoyancy,TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - Box tbz = mfi.tilebox(); + const Array4 & cell_data = S_data[IntVar::cons].array(mfi); + const Array4 & cell_prim = S_prim.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); - // 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); + 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]); - const Array4 & cell_data = S_data[IntVar::cons].array(mfi); - const Array4< Real> & buoyancy_fab = buoyancy.array(mfi); + 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)); - // Base state density - const Array4& r0_arr = r0->const_array(mfi); + Real qplus = 0.61 * qv_data(i,j,k) - (qc_data(i,j,k)+ qi_data(i,j,k)+ cell_prim(i,j,k,PrimQ2_comp)) + + (tempp3d-tempp1d)/tempp1d; - amrex::ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - Real rhop_hi = cell_data(i,j,k ,Rho_comp) + cell_data(i,j,k ,RhoQv_comp) - + cell_data(i,j,k ,RhoQc_comp) - r0_arr(i,j,k ); - Real rhop_lo = cell_data(i,j,k-1,Rho_comp) + cell_data(i,j,k-1,RhoQv_comp) - + cell_data(i,j,k-1,RhoQc_comp) - r0_arr(i,j,k-1); - buoyancy_fab(i, j, k) = grav_gpu[2] * 0.5 * ( rhop_hi + rhop_lo ); - }); - } // mfi -#endif + Real qminus = 0.61 *qv_data(i,j,k-1) - (qc_data(i,j,k-1)+ qi_data(i,j,k-1)+ cell_prim(i,j,k-1,PrimQ2_comp)) + + (tempm3d-tempm1d)/tempm1d; + + 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 + } // not buoyancy_type == 1 + } // has moisture } 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 0871b4fb0..c18e4b8f2 100644 --- a/Source/TimeIntegration/ERF_make_fast_coeffs.cpp +++ b/Source/TimeIntegration/ERF_make_fast_coeffs.cpp @@ -32,6 +32,7 @@ void make_fast_coeffs (int /*level*/, const MultiFab& S_stage_prim, const MultiFab& pi_stage, // Exner function evaluted at least stage const amrex::Geometry geom, + bool l_use_moisture, bool l_use_terrain, Real gravity, Real c_p, std::unique_ptr& detJ_cc, @@ -135,17 +136,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) ); @@ -184,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 b8128f6d9..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, @@ -170,6 +170,8 @@ void erf_slow_rhs_post (int level, int finest_level, { std::array flux; + int ncomp = S_data[IntVar::cons].nComp(); + for ( MFIter mfi(S_data[IntVar::cons],TilingIfNotGPU()); mfi.isValid(); ++mfi) { const Box& tbx = mfi.tilebox(); @@ -178,7 +180,7 @@ void erf_slow_rhs_post (int level, int finest_level, // Define flux arrays for use in advection // ************************************************************************* for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { - flux[dir].resize(amrex::surroundingNodes(tbx,dir),NVAR); + flux[dir].resize(amrex::surroundingNodes(tbx,dir),ncomp); flux[dir].setVal(0.); } const GpuArray, AMREX_SPACEDIM> @@ -226,7 +228,7 @@ void erf_slow_rhs_post (int level, int finest_level, // ************************************************************************** // Here we fill the "current" data with "new" data because that is the result of the previous RK stage // ************************************************************************** - int nsv = Cons::NumVars-2; + int nsv = S_old[IntVar::cons].nComp() - 2; const amrex::GpuArray scomp_slow = { 2,0,0,0}; const amrex::GpuArray ncomp_slow = {nsv,0,0,0}; @@ -294,38 +296,23 @@ 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); -#ifdef ERF_USE_MOISTURE - start_comp = RhoQ1_comp; - num_comp = 2; - - AdvType moist_horiz_adv_type = ac.moistscal_horiz_adv_type; - AdvType moist_vert_adv_type = ac.moistscal_vert_adv_type; - - if (ac.use_efficient_advection){ - moist_horiz_adv_type = EfficientAdvType(nrk,ac.moistscal_horiz_adv_type); - moist_vert_adv_type = EfficientAdvType(nrk,ac.moistscal_vert_adv_type); - } - AdvectionSrcForScalars(tbx, start_comp, num_comp, avg_xmom, avg_ymom, avg_zmom, - cur_prim, cell_rhs, detJ_arr, dxInv, mf_m, - moist_horiz_adv_type, moist_vert_adv_type, l_use_terrain, flx_arr); - -#elif defined(ERF_USE_WARM_NO_PRECIP) - start_comp = RhoQv_comp; - num_comp = 2; + if (solverChoice.moisture_type != MoistureType::None) + { + start_comp = RhoQ1_comp; + num_comp = 2; - AdvType moist_horiz_adv_type = ac.moistscal_horiz_adv_type; - AdvType moist_vert_adv_type = ac.moistscal_vert_adv_type; + AdvType moist_horiz_adv_type = ac.moistscal_horiz_adv_type; + AdvType moist_vert_adv_type = ac.moistscal_vert_adv_type; - if (ac.use_efficient_advection){ - moist_horiz_adv_type = EfficientAdvType(nrk,solverChoice.moistscal_horiz_adv_type); - moist_vert_adv_type = EfficientAdvType(nrk,solverChoice.moistscal_vert_adv_type); + if (ac.use_efficient_advection){ + moist_horiz_adv_type = EfficientAdvType(nrk,ac.moistscal_horiz_adv_type); + moist_vert_adv_type = EfficientAdvType(nrk,ac.moistscal_vert_adv_type); + } + AdvectionSrcForScalars(tbx, start_comp, num_comp, avg_xmom, avg_ymom, avg_zmom, + cur_prim, cell_rhs, detJ_arr, dxInv, mf_m, + moist_horiz_adv_type, moist_vert_adv_type, l_use_terrain, flx_arr); } - AdvectionSrcForScalars(tbx, start_comp, num_comp, avg_xmom, avg_ymom, avg_zmom, - cur_prim, cell_rhs, detJ_arr, dxInv, mf_m, - moist_horiz_adv_type, moist_vert_adv_type, l_use_terrain, flx_arr); -#endif - if (l_use_diff) { Array4 diffflux_x = dflux_x->array(mfi); Array4 diffflux_y = dflux_y->array(mfi); @@ -404,21 +391,19 @@ void erf_slow_rhs_post (int level, int finest_level, new_cons, cell_rhs, mf_u, mf_v, false, false); } } -#if defined(ERF_USE_NETCDF) && (defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP)) - // Zero moist RHS in set region - if (moist_zero) { - Box bx_xlo, bx_xhi, bx_ylo, bx_yhi; - compute_interior_ghost_bxs_xy(tbx, domain, width, 0, - bx_xlo, bx_xhi, - bx_ylo, bx_yhi); - int icomp; -#if defined(ERF_USE_MOISTURE) - icomp = RhoQ1_comp; -#elif defined(ERF_USE_WARM_NO_PRECIP) - icomp = RhoQv_comp; -#endif - wrfbdy_zero_rhs_in_set_region(icomp, 1, bx_xlo, bx_xhi, bx_ylo, bx_yhi, cell_rhs); - } +#if defined(ERF_USE_NETCDF) + if (solverChoice.moisture_type != MoistureType::None) + // Zero moist RHS in set region + if (moist_zero) { + Box bx_xlo, bx_xhi, bx_ylo, bx_yhi; + compute_interior_ghost_bxs_xy(tbx, domain, width, 0, + bx_xlo, bx_xhi, + bx_ylo, bx_yhi); + int icomp; + icomp = RhoQ1_comp; + wrfbdy_zero_rhs_in_set_region(icomp, 1, bx_xlo, bx_xhi, bx_ylo, bx_yhi, cell_rhs); + } // moist_zero + } // moisture_type #endif // NOTE: Computing the RHS is done over bx (union w/ grids to evolve). @@ -549,7 +534,7 @@ void erf_slow_rhs_post (int level, int finest_level, // We only add to the flux registers in the final RK step if (l_reflux && nrk == 2) { int strt_comp_reflux = RhoTheta_comp + 1; - int num_comp_reflux = NVAR - strt_comp_reflux; + int num_comp_reflux = ncomp - strt_comp_reflux; if (level < finest_level) { fr_as_crse->CrseAdd(mfi, {{AMREX_D_DECL(&(flux[0]), &(flux[1]), &(flux[2]))}}, diff --git a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp index b406a7028..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, @@ -138,6 +136,8 @@ void erf_slow_rhs_pre (int level, int finest_level, tc.pbl_type == PBLType::MYNN25 || tc.pbl_type == PBLType::YSU ); + const bool use_moisture = (solverChoice.moisture_type != MoistureType::None); + const amrex::BCRec* bc_ptr = domain_bcs_type_d.data(); const amrex::BCRec* bc_ptr_h = domain_bcs_type.data(); @@ -261,12 +261,12 @@ void erf_slow_rhs_pre (int level, int finest_level, BL_PROFILE("slow_rhs_making_er_T"); // First create Omega using velocity (not momentum) Box gbxo = surroundingNodes(bxcc,2); - amrex::ParallelFor(gbxo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + ParallelFor(gbxo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { omega_arr(i,j,k) = (k == 0) ? 0. : OmegaFromW(i,j,k,w(i,j,k),u,v,z_nd,dxInv); }); - amrex::ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { Real met_u_h_zeta_hi = Compute_h_zeta_AtIface(i+1, j , k, dxInv, z_nd); @@ -307,7 +307,7 @@ void erf_slow_rhs_pre (int level, int finest_level, // and in the first RK stage (TKE tendencies constant for nrk>0, following WRF) if ((nrk==0) && (tc.les_type == LESType::Deardorff)) { SmnSmn_a = SmnSmn->array(mfi); - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { SmnSmn_a(i,j,k) = ComputeSmnSmn(i,j,k,s11,s22,s33,s12,s13,s23); }); @@ -349,14 +349,14 @@ void erf_slow_rhs_pre (int level, int finest_level, if (bxcc.bigEnd(1) == valid_bx.bigEnd(1)) bxcc.growHi(1, 1); // Copy from temp FABs back to tau - amrex::ParallelFor(bxcc, + ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { tau11(i,j,k) = s11(i,j,k); tau22(i,j,k) = s22(i,j,k); tau33(i,j,k) = s33(i,j,k); }); - amrex::ParallelFor(tbxxy, tbxxz, tbxyz, + ParallelFor(tbxxy, tbxxz, tbxyz, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { tau12(i,j,k) = s12(i,j,k); tau21(i,j,k) = s21(i,j,k); @@ -378,7 +378,7 @@ void erf_slow_rhs_pre (int level, int finest_level, //----------------------------------------- { BL_PROFILE("slow_rhs_making_er_N"); - amrex::ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { Real mfsq = mf_m(i,j,0)*mf_m(i,j,0); er_arr(i,j,k) = (u(i+1, j , k )/mf_u(i+1,j,0) - u(i, j, k)/mf_u(i,j,0))*dxInv[0]*mfsq + (v(i , j+1, k )/mf_v(i,j+1,0) - v(i, j, k)/mf_v(i,j,0))*dxInv[1]*mfsq + @@ -404,7 +404,7 @@ void erf_slow_rhs_pre (int level, int finest_level, // and in the first RK stage (TKE tendencies constant for nrk>0, following WRF) if ((nrk==0) && (tc.les_type == LESType::Deardorff)) { SmnSmn_a = SmnSmn->array(mfi); - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { SmnSmn_a(i,j,k) = ComputeSmnSmn(i,j,k,s11,s22,s33,s12,s13,s23); }); @@ -442,13 +442,13 @@ void erf_slow_rhs_pre (int level, int finest_level, if (bxcc.bigEnd(1) == valid_bx.bigEnd(1)) bxcc.growHi(1, 1); // Copy from temp FABs back to tau - amrex::ParallelFor(bxcc, + ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { tau11(i,j,k) = s11(i,j,k); tau22(i,j,k) = s22(i,j,k); tau33(i,j,k) = s33(i,j,k); }); - amrex::ParallelFor(tbxxy, tbxxz, tbxyz, + ParallelFor(tbxxy, tbxxz, tbxyz, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { tau12(i,j,k) = s12(i,j,k); }, @@ -542,24 +542,19 @@ 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"); - amrex::ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { //if (cell_data(i,j,k,RhoTheta_comp) < 0.) printf("BAD THETA AT %d %d %d %e %e \n", // 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.); -#if defined(ERF_USE_MOISTURE) - Real qv_for_p = qv_arr(i,j,k); -#elif defined(ERF_USE_WARM_NO_PRECIP) - Real qv_for_p = cell_data(i,j,k,RhoQv_comp) / cell_data(i,j,k,Rho_comp); -#else - Real qv_for_p = 0.; -#endif + Real qv_for_p = 0.0; + 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); }); } // end profile @@ -574,17 +569,17 @@ void erf_slow_rhs_pre (int level, int finest_level, if (l_use_terrain) { Box gbxo_lo = gbxo; gbxo_lo.setBig(2,0); - amrex::ParallelFor(gbxo_lo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + ParallelFor(gbxo_lo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { omega_arr(i,j,k) = 0.; }); Box gbxo_hi = gbxo; gbxo_hi.setSmall(2,gbxo.bigEnd(2)); - amrex::ParallelFor(gbxo_hi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + ParallelFor(gbxo_hi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { omega_arr(i,j,k) = rho_w(i,j,k); }); if (z_t) { Box gbxo_mid = gbxo; gbxo_mid.setSmall(2,1); gbxo_mid.setBig(2,gbxo.bigEnd(2)-1); - amrex::ParallelFor(gbxo_mid, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + ParallelFor(gbxo_mid, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { // We define rho on the z-face the same way as in MomentumToVelocity/VelocityToMomentum Real rho_at_face = 0.5 * (cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp)); omega_arr(i,j,k) = OmegaFromW(i,j,k,rho_w(i,j,k),rho_u,rho_v,z_nd,dxInv) - @@ -592,12 +587,12 @@ void erf_slow_rhs_pre (int level, int finest_level, }); } else { Box gbxo_mid = gbxo; gbxo_mid.setSmall(2,1); gbxo_mid.setBig(2,gbxo.bigEnd(2)-1); - amrex::ParallelFor(gbxo_mid, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + ParallelFor(gbxo_mid, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { omega_arr(i,j,k) = OmegaFromW(i,j,k,rho_w(i,j,k),rho_u,rho_v,z_nd,dxInv); }); } } else { - amrex::ParallelFor(gbxo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + ParallelFor(gbxo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { omega_arr(i,j,k) = rho_w(i,j,k); }); } @@ -689,12 +684,12 @@ void erf_slow_rhs_pre (int level, int finest_level, { auto const& src_arr = source.const_array(mfi); if (l_use_terrain && l_moving_terrain) { - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { cell_rhs(i,j,k,RhoTheta_comp) += src_arr(i,j,k,RhoTheta_comp) / detJ_arr(i,j,k); }); } else { - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { cell_rhs(i,j,k,RhoTheta_comp) += src_arr(i,j,k,RhoTheta_comp); }); @@ -706,7 +701,7 @@ void erf_slow_rhs_pre (int level, int finest_level, int n = RhoTheta_comp; int nr = Rho_comp; int np = PrimTheta_comp; - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { Real theta = cell_prim(i,j,k,np); cell_rhs(i, j, k, n) -= dptr_rayleigh_tau[k] * (theta - dptr_rayleigh_thetabar[k]) * cell_data(i,j,k,nr); @@ -715,7 +710,7 @@ void erf_slow_rhs_pre (int level, int finest_level, // Multiply the slow RHS for rho and rhotheta by detJ here so we don't have to later if (l_use_terrain && l_moving_terrain) { - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { cell_rhs(i,j,k,Rho_comp) *= detJ_arr(i,j,k); cell_rhs(i,j,k,RhoTheta_comp) *= detJ_arr(i,j,k); @@ -775,7 +770,7 @@ void erf_slow_rhs_pre (int level, int finest_level, // TERRAIN VERSION // ****************************************************************** if (l_use_terrain) { - amrex::ParallelFor(tbx, + ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // x-momentum equation @@ -803,13 +798,11 @@ void erf_slow_rhs_pre (int level, int finest_level, gpx *= mf_u(i,j,0); Real q = 0.0; -#if defined(ERF_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) ); -#elif defined(ERF_USE_WARM_NO_PRECIP) - q = 0.5 * ( cell_prim(i,j,k,PrimQv_comp) + cell_prim(i-1,j,k,PrimQv_comp) - +cell_prim(i,j,k,PrimQc_comp) + cell_prim(i-1,j,k,PrimQc_comp) ); -#endif + 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) ); + } + rho_u_rhs(i, j, k) += (-gpx - abl_pressure_grad[0]) / (1.0 + q) + rho_u_face * abl_geo_forcing[0]; @@ -838,7 +831,7 @@ void erf_slow_rhs_pre (int level, int finest_level, // ****************************************************************** // NON-TERRAIN VERSION // ****************************************************************** - amrex::ParallelFor(tbx, + ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // x-momentum equation @@ -847,13 +840,11 @@ void erf_slow_rhs_pre (int level, int finest_level, gpx *= mf_u(i,j,0); Real q = 0.0; -#if defined(ERF_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) ); -#elif defined(ERF_USE_WARM_NO_PRECIP) - q = 0.5 * ( cell_prim(i,j,k,PrimQv_comp) + cell_prim(i-1,j,k,PrimQv_comp) - +cell_prim(i,j,k,PrimQc_comp) + cell_prim(i-1,j,k,PrimQc_comp) ); -#endif + 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) ); + } + rho_u_rhs(i, j, k) += (-gpx - abl_pressure_grad[0]) / (1.0 + q) + rho_u_face * abl_geo_forcing[0]; @@ -882,7 +873,7 @@ void erf_slow_rhs_pre (int level, int finest_level, // TERRAIN VERSION // ****************************************************************** if (l_use_terrain) { - amrex::ParallelFor(tby, + ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // y-momentum equation @@ -911,13 +902,10 @@ void erf_slow_rhs_pre (int level, int finest_level, gpy *= mf_v(i,j,0); Real q = 0.0; -#if defined(ERF_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) ); -#elif defined(ERF_USE_WARM_NO_PRECIP) - q = 0.5 * ( cell_prim(i,j,k,PrimQv_comp) + cell_prim(i,j-1,k,PrimQv_comp) - +cell_prim(i,j,k,PrimQc_comp) + cell_prim(i,j-1,k,PrimQc_comp) ); -#endif + 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) ); + } rho_v_rhs(i, j, k) += (-gpy - abl_pressure_grad[1]) / (1.0_rt + q) + rho_v_face * abl_geo_forcing[1]; @@ -944,7 +932,7 @@ void erf_slow_rhs_pre (int level, int finest_level, // NON-TERRAIN VERSION // ****************************************************************** } else { - amrex::ParallelFor(tby, + ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // y-momentum equation @@ -953,13 +941,11 @@ void erf_slow_rhs_pre (int level, int finest_level, gpy *= mf_v(i,j,0); Real q = 0.0; -#if defined(ERF_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) ); -#elif defined(ERF_USE_WARM_NO_PRECIP) - q = 0.5 * ( cell_prim(i,j,k,PrimQv_comp) + cell_prim(i,j-1,k,PrimQv_comp) - +cell_prim(i,j,k,PrimQc_comp) + cell_prim(i,j-1,k,PrimQc_comp) ); -#endif + 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) ); + } + rho_v_rhs(i, j, k) += (-gpy - abl_pressure_grad[1]) / (1.0_rt + q) + rho_v_face * abl_geo_forcing[1]; @@ -986,7 +972,7 @@ void erf_slow_rhs_pre (int level, int finest_level, b2d.setSmall(2,0); b2d.setBig(2,0); // Enforce no forcing term at top and bottom boundaries - amrex::ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int) { + ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int) { rho_w_rhs(i,j, 0) = 0.; rho_w_rhs(i,j,domhi_z+1) = 0.; // TODO: generalize this }); @@ -995,25 +981,22 @@ void erf_slow_rhs_pre (int level, int finest_level, { BL_PROFILE("slow_rhs_pre_zmom"); auto rayleigh_damp_W = solverChoice.rayleigh_damp_W; + // ****************************************************************** // TERRAIN VERSION // ****************************************************************** if (l_use_terrain) { - amrex::ParallelFor(tbz, - [=] AMREX_GPU_DEVICE (int i, int j, int k) { // z-momentum equation + ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // z-momentum equation Real rho_w_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) ); Real met_h_zeta = Compute_h_zeta_AtKface(i, j, k, dxInv, z_nd); Real gpz = dxInv[2] * ( pp_arr(i,j,k)-pp_arr(i,j,k-1) ) / met_h_zeta; Real q = 0.0; -#if defined(ERF_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) ); -#elif defined(ERF_USE_WARM_NO_PRECIP) - q = 0.5 * ( cell_prim(i,j,k,PrimQv_comp) + cell_prim(i,j,k-1,PrimQv_comp) - +cell_prim(i,j,k,PrimQc_comp) + cell_prim(i,j,k-1,PrimQc_comp) ); -#endif + 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) ); + } rho_w_rhs(i, j, k) += (buoyancy_fab(i,j,k) - gpz - abl_pressure_grad[2]) / (1.0_rt + q) + rho_w_face * abl_geo_forcing[2]; @@ -1040,21 +1023,17 @@ void erf_slow_rhs_pre (int level, int finest_level, // NON-TERRAIN VERSION // ****************************************************************** } else { - amrex::ParallelFor(tbz, - [=] AMREX_GPU_DEVICE (int i, int j, int k) + ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // z-momentum equation Real rho_w_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) ); Real gpz = dxInv[2] * ( pp_arr(i,j,k)-pp_arr(i,j,k-1) ); Real q = 0.0; -#if defined(ERF_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) ); -#elif defined(ERF_USE_WARM_NO_PRECIP) - q = 0.5 * ( cell_prim(i,j,k,PrimQv_comp) + cell_prim(i,j,k-1,PrimQv_comp) - +cell_prim(i,j,k,PrimQc_comp) + cell_prim(i,j,k-1,PrimQc_comp) ); -#endif + 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) ); + } rho_w_rhs(i, j, k) += (buoyancy_fab(i,j,k) - gpz - abl_pressure_grad[2]) / (1.0_rt + q) + rho_w_face * abl_geo_forcing[2]; 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 fa990a03a..ff6289540 100644 --- a/Source/TimeIntegration/TI_fast_rhs_fun.H +++ b/Source/TimeIntegration/TI_fast_rhs_fun.H @@ -15,6 +15,8 @@ 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 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 // Per p2902 of Klemp-Skamarock-Dudhia-2007 @@ -84,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, - 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) { @@ -98,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, @@ -110,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, - 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 @@ -126,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, @@ -134,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, - 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 @@ -150,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, @@ -158,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 f4840eecf..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 @@ -202,6 +200,7 @@ void make_fast_coeffs (int level, const amrex::MultiFab& S_stage_prim, const amrex::MultiFab& pi_stage, const amrex::Geometry geom, + const bool use_moisture, const bool use_terrain, const amrex::Real gravity, const amrex::Real c_p, @@ -219,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); @@ -245,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 bbb08ea83..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 + qmoist[level], qv_d, qc_d, qi_d, 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 + qmoist[level], z_t_rk[level], Omega, source, buoyancy, Tau11, Tau22, Tau33, Tau12, Tau13, Tau21, Tau23, Tau31, Tau32, SmnSmn, eddyDiffs, Hfx3, Diss, diff --git a/Source/TimeIntegration/TI_utils.H b/Source/TimeIntegration/TI_utils.H index 88b0b587a..953b60f79 100644 --- a/Source/TimeIntegration/TI_utils.H +++ b/Source/TimeIntegration/TI_utils.H @@ -4,6 +4,9 @@ auto cons_to_prim = [&](const MultiFab& cons_state, int ng) { BL_PROFILE("cons_to_prim()"); + + int ncomp_prim = S_prim.nComp(); + #ifdef _OPENMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) #endif @@ -20,7 +23,7 @@ Real rho_theta = cons_arr(i,j,k,RhoTheta_comp); prim_arr(i,j,k,PrimTheta_comp) = rho_theta / rho; pi_stage_arr(i,j,k) = getExnergivenRTh(rho_theta, rdOcp); - for (int n = 1; n < NUM_PRIM; ++n) { + for (int n = 1; n < ncomp_prim; ++n) { prim_arr(i,j,k,PrimTheta_comp + n) = cons_arr(i,j,k,RhoTheta_comp + n) / rho; } });