diff --git a/Exec/ABL/prob.H b/Exec/ABL/prob.H index af697ac01..b22c7d56b 100644 --- a/Exec/ABL/prob.H +++ b/Exec/ABL/prob.H @@ -62,9 +62,6 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, - amrex::Array4 const& qv, - amrex::Array4 const& qc, - amrex::Array4 const& qi, amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/ABL/prob.cpp b/Exec/ABL/prob.cpp index bf767fe96..c0adb291b 100644 --- a/Exec/ABL/prob.cpp +++ b/Exec/ABL/prob.cpp @@ -56,9 +56,6 @@ Problem::init_custom_pert( amrex::Array4 const& /*p_hse*/, amrex::Array4 const& /*z_nd*/, amrex::Array4 const& /*z_cc*/, - amrex::Array4 const& /*qv*/, - amrex::Array4 const& /*qc*/, - amrex::Array4 const& /*qi*/, amrex::GeometryData const& geomdata, amrex::Array4 const& /*mf_m*/, amrex::Array4 const& /*mf_u*/, diff --git a/Exec/ABL_input_sounding/prob.H b/Exec/ABL_input_sounding/prob.H index f3e709ece..c0a8a4e74 100644 --- a/Exec/ABL_input_sounding/prob.H +++ b/Exec/ABL_input_sounding/prob.H @@ -61,9 +61,6 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, - amrex::Array4 const& qv, - amrex::Array4 const& qc, - amrex::Array4 const& qi, amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/ABL_input_sounding/prob.cpp b/Exec/ABL_input_sounding/prob.cpp index 80911c884..80dddebe7 100644 --- a/Exec/ABL_input_sounding/prob.cpp +++ b/Exec/ABL_input_sounding/prob.cpp @@ -55,9 +55,6 @@ Problem::init_custom_pert( amrex::Array4 const& /*p_hse*/, amrex::Array4 const& /*z_nd*/, amrex::Array4 const& /*z_cc*/, - amrex::Array4 const& /*qv*/, - amrex::Array4 const& /*qc*/, - amrex::Array4 const& /*qi*/, amrex::GeometryData const& geomdata, amrex::Array4 const& /*mf_m*/, amrex::Array4 const& /*mf_u*/, diff --git a/Exec/DevTests/MiguelDev/prob.H b/Exec/DevTests/MiguelDev/prob.H index 468a61cc3..2bbf9bb93 100644 --- a/Exec/DevTests/MiguelDev/prob.H +++ b/Exec/DevTests/MiguelDev/prob.H @@ -51,9 +51,6 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, - amrex::Array4 const& qv, - amrex::Array4 const& qc, - amrex::Array4 const& qi, amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/DevTests/MiguelDev/prob.cpp b/Exec/DevTests/MiguelDev/prob.cpp index 5db590f00..035160882 100644 --- a/Exec/DevTests/MiguelDev/prob.cpp +++ b/Exec/DevTests/MiguelDev/prob.cpp @@ -47,9 +47,6 @@ Problem::init_custom_pert( Array4 const&, Array4 const&, Array4 const&, - Array4 const&, - Array4 const&, - Array4 const&, amrex::GeometryData const& geomdata, Array4 const&, Array4 const&, diff --git a/Exec/DevTests/MovingTerrain/prob.H b/Exec/DevTests/MovingTerrain/prob.H index 495a7eea7..fea7bd7f9 100644 --- a/Exec/DevTests/MovingTerrain/prob.H +++ b/Exec/DevTests/MovingTerrain/prob.H @@ -32,9 +32,6 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, - amrex::Array4 const& qv, - amrex::Array4 const& qc, - amrex::Array4 const& qi, amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/DevTests/MovingTerrain/prob.cpp b/Exec/DevTests/MovingTerrain/prob.cpp index d4c5181bc..9a25dc3fc 100644 --- a/Exec/DevTests/MovingTerrain/prob.cpp +++ b/Exec/DevTests/MovingTerrain/prob.cpp @@ -36,9 +36,6 @@ Problem::init_custom_pert( Array4 const& p_hse, Array4 const& z_nd, Array4 const& z_cc, - Array4 const&, - Array4 const&, - Array4 const&, GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, diff --git a/Exec/DevTests/MultiBlock/prob.H b/Exec/DevTests/MultiBlock/prob.H index e44d13361..4b58c8691 100644 --- a/Exec/DevTests/MultiBlock/prob.H +++ b/Exec/DevTests/MultiBlock/prob.H @@ -40,9 +40,6 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, - amrex::Array4 const& qv, - amrex::Array4 const& qc, - amrex::Array4 const& qi, amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/DevTests/MultiBlock/prob.cpp b/Exec/DevTests/MultiBlock/prob.cpp index 26287abe0..20be90ef0 100644 --- a/Exec/DevTests/MultiBlock/prob.cpp +++ b/Exec/DevTests/MultiBlock/prob.cpp @@ -40,9 +40,6 @@ Problem::init_custom_pert( Array4 const& p_hse, Array4 const&, Array4 const&, - Array4 const&, - Array4 const&, - Array4 const&, GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, diff --git a/Exec/DevTests/ParticlesOverWoA/prob.H b/Exec/DevTests/ParticlesOverWoA/prob.H index b0e9c7c05..71f5ef770 100644 --- a/Exec/DevTests/ParticlesOverWoA/prob.H +++ b/Exec/DevTests/ParticlesOverWoA/prob.H @@ -40,9 +40,6 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, - amrex::Array4 const& qv, - amrex::Array4 const& qc, - amrex::Array4 const& qi, amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/DevTests/ParticlesOverWoA/prob.cpp b/Exec/DevTests/ParticlesOverWoA/prob.cpp index e81c90a14..e5699da30 100644 --- a/Exec/DevTests/ParticlesOverWoA/prob.cpp +++ b/Exec/DevTests/ParticlesOverWoA/prob.cpp @@ -41,9 +41,6 @@ Problem::init_custom_pert( Array4 const& p_hse, Array4 const& z_nd, Array4 const& z_cc, - Array4 const&, - Array4 const&, - Array4 const&, GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, diff --git a/Exec/LLJ/prob.H b/Exec/LLJ/prob.H index e3e7ae0b8..0d460d496 100644 --- a/Exec/LLJ/prob.H +++ b/Exec/LLJ/prob.H @@ -33,9 +33,6 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, - amrex::Array4 const& qv, - amrex::Array4 const& qc, - amrex::Array4 const& qi, amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/LLJ/prob.cpp b/Exec/LLJ/prob.cpp index d2d78c631..6f2ad9f38 100644 --- a/Exec/LLJ/prob.cpp +++ b/Exec/LLJ/prob.cpp @@ -34,9 +34,6 @@ Problem::init_custom_pert( Array4 const&, Array4 const&, Array4 const&, - Array4 const&, - Array4 const&, - Array4 const&, GeometryData const&, Array4 const&, Array4 const&, diff --git a/Exec/Radiation/prob.H b/Exec/Radiation/prob.H index 01cc85b99..77933842f 100644 --- a/Exec/Radiation/prob.H +++ b/Exec/Radiation/prob.H @@ -40,9 +40,6 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, - amrex::Array4 const& qv, - amrex::Array4 const& qc, - amrex::Array4 const& qi, amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/Radiation/prob.cpp b/Exec/Radiation/prob.cpp index 2016fa9d9..23c967ad4 100644 --- a/Exec/Radiation/prob.cpp +++ b/Exec/Radiation/prob.cpp @@ -96,9 +96,6 @@ Problem::init_custom_pert( Array4 const& p_hse, Array4 const& /*z_nd*/, Array4 const& /*z_cc*/, - Array4 const& qv, - Array4 const& qc, - Array4 const& qi, GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, @@ -187,13 +184,6 @@ Problem::init_custom_pert( 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; -#endif } }); diff --git a/Exec/RegTests/Bubble/inputs_squall2d_x b/Exec/RegTests/Bubble/inputs_squall2d_x index d561668c4..f62ae2fe3 100644 --- a/Exec/RegTests/Bubble/inputs_squall2d_x +++ b/Exec/RegTests/Bubble/inputs_squall2d_x @@ -47,6 +47,7 @@ erf.use_rayleigh_damping = false # TODO: enable erf.les_type = "Deardorff" erf.pbl_type = "None" erf.moisture_model = "SAM" +erf.buoyancy_type = 1 erf.molec_diff_type = "Constant" erf.rho0_trans = 1.0 # [kg/m^3], used to convert input diffusivities diff --git a/Exec/RegTests/Bubble/prob.H b/Exec/RegTests/Bubble/prob.H index 0218aff1d..e8fe13dd7 100644 --- a/Exec/RegTests/Bubble/prob.H +++ b/Exec/RegTests/Bubble/prob.H @@ -56,9 +56,6 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, - amrex::Array4 const& qv, - amrex::Array4 const& qc, - amrex::Array4 const& qi, amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/RegTests/Bubble/prob.cpp b/Exec/RegTests/Bubble/prob.cpp index 1056f6e44..a70d87a1f 100644 --- a/Exec/RegTests/Bubble/prob.cpp +++ b/Exec/RegTests/Bubble/prob.cpp @@ -47,9 +47,6 @@ Problem::init_custom_pert( Array4 const& p_hse, Array4 const& /*z_nd*/, Array4 const& z_cc, - Array4 const&, - Array4 const&, - Array4 const&, GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, diff --git a/Exec/RegTests/CouetteFlow/prob.H b/Exec/RegTests/CouetteFlow/prob.H index 6bfe6a227..570d2ef92 100644 --- a/Exec/RegTests/CouetteFlow/prob.H +++ b/Exec/RegTests/CouetteFlow/prob.H @@ -35,9 +35,6 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, - amrex::Array4 const& qv, - amrex::Array4 const& qc, - amrex::Array4 const& qi, amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/RegTests/CouetteFlow/prob.cpp b/Exec/RegTests/CouetteFlow/prob.cpp index a160922ad..0bc9be19e 100644 --- a/Exec/RegTests/CouetteFlow/prob.cpp +++ b/Exec/RegTests/CouetteFlow/prob.cpp @@ -37,9 +37,6 @@ Problem::init_custom_pert( Array4 const&, Array4 const&, Array4 const&, - Array4 const&, - Array4 const&, - Array4 const&, amrex::GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, diff --git a/Exec/RegTests/DensityCurrent/prob.H b/Exec/RegTests/DensityCurrent/prob.H index c6ca99772..7d17e9d8a 100644 --- a/Exec/RegTests/DensityCurrent/prob.H +++ b/Exec/RegTests/DensityCurrent/prob.H @@ -38,9 +38,6 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, - amrex::Array4 const& qv, - amrex::Array4 const& qc, - amrex::Array4 const& qi, amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/RegTests/DensityCurrent/prob.cpp b/Exec/RegTests/DensityCurrent/prob.cpp index 751b31ad1..04ec271ce 100644 --- a/Exec/RegTests/DensityCurrent/prob.cpp +++ b/Exec/RegTests/DensityCurrent/prob.cpp @@ -40,9 +40,6 @@ Problem::init_custom_pert( Array4 const& p_hse, Array4 const& /*z_nd*/, Array4 const& z_cc, - Array4 const&, - Array4 const&, - Array4 const&, GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, diff --git a/Exec/RegTests/DynamicRefinement/prob.H b/Exec/RegTests/DynamicRefinement/prob.H index b53f1d829..1755d4979 100644 --- a/Exec/RegTests/DynamicRefinement/prob.H +++ b/Exec/RegTests/DynamicRefinement/prob.H @@ -45,9 +45,6 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, - amrex::Array4 const& qv, - amrex::Array4 const& qc, - amrex::Array4 const& qi, amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/RegTests/DynamicRefinement/prob.cpp b/Exec/RegTests/DynamicRefinement/prob.cpp index fa911cc0d..0b2b549ad 100644 --- a/Exec/RegTests/DynamicRefinement/prob.cpp +++ b/Exec/RegTests/DynamicRefinement/prob.cpp @@ -84,9 +84,6 @@ Problem::init_custom_pert( Array4 const&, Array4 const&, Array4 const&, - Array4 const&, - Array4 const&, - Array4 const&, amrex::GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, diff --git a/Exec/RegTests/EkmanSpiral_custom/prob.H b/Exec/RegTests/EkmanSpiral_custom/prob.H index 040c0af18..fedf5e9d9 100644 --- a/Exec/RegTests/EkmanSpiral_custom/prob.H +++ b/Exec/RegTests/EkmanSpiral_custom/prob.H @@ -33,9 +33,6 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, - amrex::Array4 const& qv, - amrex::Array4 const& qc, - amrex::Array4 const& qi, amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/RegTests/EkmanSpiral_custom/prob.cpp b/Exec/RegTests/EkmanSpiral_custom/prob.cpp index 1c9c8810b..e608032ee 100644 --- a/Exec/RegTests/EkmanSpiral_custom/prob.cpp +++ b/Exec/RegTests/EkmanSpiral_custom/prob.cpp @@ -34,9 +34,6 @@ Problem::init_custom_pert( Array4 const&, Array4 const&, Array4 const&, - Array4 const&, - Array4 const&, - Array4 const&, GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, diff --git a/Exec/RegTests/EkmanSpiral_ideal/prob.H b/Exec/RegTests/EkmanSpiral_ideal/prob.H index 9345870cd..192d96138 100644 --- a/Exec/RegTests/EkmanSpiral_ideal/prob.H +++ b/Exec/RegTests/EkmanSpiral_ideal/prob.H @@ -33,9 +33,6 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, - amrex::Array4 const& qv, - amrex::Array4 const& qc, - amrex::Array4 const& qi, amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/RegTests/EkmanSpiral_ideal/prob.cpp b/Exec/RegTests/EkmanSpiral_ideal/prob.cpp index 7cdc1e5cf..c69999a7b 100644 --- a/Exec/RegTests/EkmanSpiral_ideal/prob.cpp +++ b/Exec/RegTests/EkmanSpiral_ideal/prob.cpp @@ -34,9 +34,6 @@ Problem::init_custom_pert( Array4 const&, Array4 const&, Array4 const&, - Array4 const&, - Array4 const&, - Array4 const&, GeometryData const&, Array4 const&, Array4 const&, diff --git a/Exec/RegTests/EkmanSpiral_input_sounding/prob.H b/Exec/RegTests/EkmanSpiral_input_sounding/prob.H index 10d9a8836..cd96fa98b 100644 --- a/Exec/RegTests/EkmanSpiral_input_sounding/prob.H +++ b/Exec/RegTests/EkmanSpiral_input_sounding/prob.H @@ -33,9 +33,6 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, - amrex::Array4 const& qv, - amrex::Array4 const& qc, - amrex::Array4 const& qi, amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/RegTests/EkmanSpiral_input_sounding/prob.cpp b/Exec/RegTests/EkmanSpiral_input_sounding/prob.cpp index 0e691087e..ab357578e 100644 --- a/Exec/RegTests/EkmanSpiral_input_sounding/prob.cpp +++ b/Exec/RegTests/EkmanSpiral_input_sounding/prob.cpp @@ -34,9 +34,6 @@ Problem::init_custom_pert( Array4 const&, Array4 const&, Array4 const&, - Array4 const&, - Array4 const&, - Array4 const&, GeometryData const&, Array4 const&, Array4 const&, diff --git a/Exec/RegTests/GABLS1/prob.H b/Exec/RegTests/GABLS1/prob.H index 6017dca83..ebb4d529f 100644 --- a/Exec/RegTests/GABLS1/prob.H +++ b/Exec/RegTests/GABLS1/prob.H @@ -62,9 +62,6 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, - amrex::Array4 const& qv, - amrex::Array4 const& qc, - amrex::Array4 const& qi, amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/RegTests/GABLS1/prob.cpp b/Exec/RegTests/GABLS1/prob.cpp index 2ea9dc86e..748b65691 100644 --- a/Exec/RegTests/GABLS1/prob.cpp +++ b/Exec/RegTests/GABLS1/prob.cpp @@ -56,9 +56,6 @@ Problem::init_custom_pert( amrex::Array4 const& /*p_hse*/, amrex::Array4 const& /*z_nd*/, amrex::Array4 const& /*z_cc*/, - amrex::Array4 const& /*qv*/, - amrex::Array4 const& /*qc*/, - amrex::Array4 const& /*qi*/, amrex::GeometryData const& geomdata, amrex::Array4 const& /*mf_m*/, amrex::Array4 const& /*mf_u*/, diff --git a/Exec/RegTests/IsentropicVortex/prob.H b/Exec/RegTests/IsentropicVortex/prob.H index dfb671e71..d2b51270f 100644 --- a/Exec/RegTests/IsentropicVortex/prob.H +++ b/Exec/RegTests/IsentropicVortex/prob.H @@ -44,9 +44,6 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, - amrex::Array4 const& qv, - amrex::Array4 const& qc, - amrex::Array4 const& qi, amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/RegTests/IsentropicVortex/prob.cpp b/Exec/RegTests/IsentropicVortex/prob.cpp index 52442eae3..9c9aef46c 100644 --- a/Exec/RegTests/IsentropicVortex/prob.cpp +++ b/Exec/RegTests/IsentropicVortex/prob.cpp @@ -86,9 +86,6 @@ Problem::init_custom_pert( Array4 const& p_hse, Array4 const&, Array4 const&, - Array4 const&, - Array4 const&, - Array4 const&, amrex::GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, diff --git a/Exec/RegTests/PoiseuilleFlow/prob.H b/Exec/RegTests/PoiseuilleFlow/prob.H index 2ab2bff7b..120b23372 100644 --- a/Exec/RegTests/PoiseuilleFlow/prob.H +++ b/Exec/RegTests/PoiseuilleFlow/prob.H @@ -34,9 +34,6 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, - amrex::Array4 const& qv, - amrex::Array4 const& qc, - amrex::Array4 const& qi, amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/RegTests/PoiseuilleFlow/prob.cpp b/Exec/RegTests/PoiseuilleFlow/prob.cpp index f2b2d3d09..3defea69d 100644 --- a/Exec/RegTests/PoiseuilleFlow/prob.cpp +++ b/Exec/RegTests/PoiseuilleFlow/prob.cpp @@ -36,9 +36,6 @@ Problem::init_custom_pert( amrex::Array4 const&, amrex::Array4 const&, amrex::Array4 const&, - Array4 const&, - Array4 const&, - Array4 const&, amrex::GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, diff --git a/Exec/RegTests/ScalarAdvDiff/prob.H b/Exec/RegTests/ScalarAdvDiff/prob.H index 684050b4f..5adca8818 100644 --- a/Exec/RegTests/ScalarAdvDiff/prob.H +++ b/Exec/RegTests/ScalarAdvDiff/prob.H @@ -46,9 +46,6 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, - amrex::Array4 const& qv, - amrex::Array4 const& qc, - amrex::Array4 const& qi, amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/RegTests/ScalarAdvDiff/prob.cpp b/Exec/RegTests/ScalarAdvDiff/prob.cpp index 1ea764346..f9e931dfc 100644 --- a/Exec/RegTests/ScalarAdvDiff/prob.cpp +++ b/Exec/RegTests/ScalarAdvDiff/prob.cpp @@ -70,9 +70,6 @@ Problem::init_custom_pert( Array4 const&, Array4 const&, Array4 const&, - Array4 const&, - Array4 const&, - Array4 const&, GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, diff --git a/Exec/RegTests/StokesSecondProblem/prob.cpp b/Exec/RegTests/StokesSecondProblem/prob.cpp index 11cbf0263..cf141fa62 100644 --- a/Exec/RegTests/StokesSecondProblem/prob.cpp +++ b/Exec/RegTests/StokesSecondProblem/prob.cpp @@ -32,9 +32,6 @@ Problem::init_custom_pert( Array4 const& p_hse, Array4 const&, Array4 const&, - Array4 const&, - Array4 const&, - Array4 const&, GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, diff --git a/Exec/RegTests/TaylorGreenVortex/prob.H b/Exec/RegTests/TaylorGreenVortex/prob.H index 5fca0d226..0635cacbb 100644 --- a/Exec/RegTests/TaylorGreenVortex/prob.H +++ b/Exec/RegTests/TaylorGreenVortex/prob.H @@ -33,9 +33,6 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, - amrex::Array4 const& qv, - amrex::Array4 const& qc, - amrex::Array4 const& qi, amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/RegTests/TaylorGreenVortex/prob.cpp b/Exec/RegTests/TaylorGreenVortex/prob.cpp index 6f244b8f2..f90eee76c 100644 --- a/Exec/RegTests/TaylorGreenVortex/prob.cpp +++ b/Exec/RegTests/TaylorGreenVortex/prob.cpp @@ -36,9 +36,6 @@ Problem::init_custom_pert( Array4 const& p_hse, Array4 const&, Array4 const&, - Array4 const&, - Array4 const&, - Array4 const&, GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, diff --git a/Exec/RegTests/Terrain2d_Cylinder/prob.H b/Exec/RegTests/Terrain2d_Cylinder/prob.H index 46cf99f0c..182e314e9 100644 --- a/Exec/RegTests/Terrain2d_Cylinder/prob.H +++ b/Exec/RegTests/Terrain2d_Cylinder/prob.H @@ -34,9 +34,6 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, - amrex::Array4 const& qv, - amrex::Array4 const& qc, - amrex::Array4 const& qi, amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/RegTests/Terrain2d_Cylinder/prob.cpp b/Exec/RegTests/Terrain2d_Cylinder/prob.cpp index 918a51ddd..6b75165e5 100644 --- a/Exec/RegTests/Terrain2d_Cylinder/prob.cpp +++ b/Exec/RegTests/Terrain2d_Cylinder/prob.cpp @@ -36,9 +36,6 @@ Problem::init_custom_pert( Array4 const&, Array4 const& z_nd, Array4 const& z_cc, - Array4 const&, - Array4 const&, - Array4 const&, GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, diff --git a/Exec/RegTests/Terrain3d_Hemisphere/prob.H b/Exec/RegTests/Terrain3d_Hemisphere/prob.H index a36f54712..64e5368e3 100644 --- a/Exec/RegTests/Terrain3d_Hemisphere/prob.H +++ b/Exec/RegTests/Terrain3d_Hemisphere/prob.H @@ -34,9 +34,6 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, - amrex::Array4 const& qv, - amrex::Array4 const& qc, - amrex::Array4 const& qi, amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/RegTests/Terrain3d_Hemisphere/prob.cpp b/Exec/RegTests/Terrain3d_Hemisphere/prob.cpp index ee58919eb..0acce57b3 100644 --- a/Exec/RegTests/Terrain3d_Hemisphere/prob.cpp +++ b/Exec/RegTests/Terrain3d_Hemisphere/prob.cpp @@ -35,9 +35,6 @@ Problem::init_custom_pert( Array4 const&, Array4 const& z_nd, Array4 const& z_cc, - Array4 const&, - Array4 const&, - Array4 const&, GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, diff --git a/Exec/RegTests/WPS_Test/inputs_real_ChisholmView b/Exec/RegTests/WPS_Test/inputs_real_ChisholmView index d077f9e88..097433d67 100644 --- a/Exec/RegTests/WPS_Test/inputs_real_ChisholmView +++ b/Exec/RegTests/WPS_Test/inputs_real_ChisholmView @@ -43,7 +43,7 @@ erf.check_int = -100 # number of timesteps between checkpoints # PLOTFILES erf.plot_file_1 = plt # prefix of plotfile name erf.plot_int_1 = 1 # number of timesteps between plotfiles -erf.plot_vars_1 = density rhoQt x_velocity y_velocity z_velocity pressure temp theta pres_hse z_phys qt qp qv qc +erf.plot_vars_1 = density rhoQ1 x_velocity y_velocity z_velocity pressure temp theta pres_hse z_phys qt qp qv qc # SOLVER CHOICE erf.alpha_T = 0.0 diff --git a/Exec/RegTests/WitchOfAgnesi/prob.H b/Exec/RegTests/WitchOfAgnesi/prob.H index 1e4322ecf..922c5f19f 100644 --- a/Exec/RegTests/WitchOfAgnesi/prob.H +++ b/Exec/RegTests/WitchOfAgnesi/prob.H @@ -31,9 +31,6 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, - amrex::Array4 const& qv, - amrex::Array4 const& qc, - amrex::Array4 const& qi, amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/RegTests/WitchOfAgnesi/prob.cpp b/Exec/RegTests/WitchOfAgnesi/prob.cpp index 56c970c6f..cf14e4aae 100644 --- a/Exec/RegTests/WitchOfAgnesi/prob.cpp +++ b/Exec/RegTests/WitchOfAgnesi/prob.cpp @@ -34,9 +34,6 @@ Problem::init_custom_pert( Array4 const& p_hse, Array4 const& z_nd, Array4 const& z_cc, - Array4 const&, - Array4 const&, - Array4 const&, GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, diff --git a/Exec/SquallLine_2D/inputs_moisture b/Exec/SquallLine_2D/inputs_moisture index f82dd6612..58a80084f 100644 --- a/Exec/SquallLine_2D/inputs_moisture +++ b/Exec/SquallLine_2D/inputs_moisture @@ -50,7 +50,7 @@ amr.check_int = 1000 # number of timesteps between checkpoints # PLOTFILES erf.plot_file_1 = plt # root name of plotfile erf.plot_int_1 = 100 # number of timesteps between plotfiles -erf.plot_vars_1 = density rhotheta rhoQt rhoQp x_velocity y_velocity z_velocity pressure theta temp qt qp qv qc qi scalar pert_dens +erf.plot_vars_1 = density rhotheta rhoQ1 rhoQ2 rhoQ3 x_velocity y_velocity z_velocity pressure theta temp qt qp qv qc qi scalar pert_dens # SOLVER CHOICE erf.use_gravity = true diff --git a/Exec/SquallLine_2D/inputs_moisture_Gabersek b/Exec/SquallLine_2D/inputs_moisture_Gabersek index ad5b642b4..a497ac6c8 100644 --- a/Exec/SquallLine_2D/inputs_moisture_Gabersek +++ b/Exec/SquallLine_2D/inputs_moisture_Gabersek @@ -50,7 +50,7 @@ amr.check_int = 1000 # number of timesteps between checkpoints # PLOTFILES erf.plot_file_1 = plt # root name of plotfile erf.plot_int_1 = 100 # number of timesteps between plotfiles -erf.plot_vars_1 = density rhotheta rhoQt rhoQp x_velocity y_velocity z_velocity pressure theta temp qt qp qv qc qi scalar pert_dens +erf.plot_vars_1 = density rhotheta rhoQ1 rhoQ2 rhoQ3 x_velocity y_velocity z_velocity pressure theta temp qt qp qv qc qi scalar pert_dens # SOLVER CHOICE erf.use_gravity = true diff --git a/Exec/SquallLine_2D/prob.H b/Exec/SquallLine_2D/prob.H index f9f1e8617..28b637c59 100644 --- a/Exec/SquallLine_2D/prob.H +++ b/Exec/SquallLine_2D/prob.H @@ -50,9 +50,6 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, - amrex::Array4 const& qv, - amrex::Array4 const& qc, - amrex::Array4 const& qi, amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/SquallLine_2D/prob.cpp b/Exec/SquallLine_2D/prob.cpp index 6d36b971e..8e5445f22 100644 --- a/Exec/SquallLine_2D/prob.cpp +++ b/Exec/SquallLine_2D/prob.cpp @@ -282,9 +282,6 @@ Problem::init_custom_pert ( Array4 const& /*p_hse*/, Array4 const& /*z_nd*/, Array4 const& /*z_cc*/, - Array4 const& qv, - Array4 const& qc, - Array4 const& qi, GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, @@ -380,13 +377,6 @@ Problem::init_custom_pert ( 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; -#endif } }); diff --git a/Exec/SuperCell/inputs_moisture b/Exec/SuperCell/inputs_moisture index 56e7babe4..19213ec05 100644 --- a/Exec/SuperCell/inputs_moisture +++ b/Exec/SuperCell/inputs_moisture @@ -38,7 +38,7 @@ amr.check_int = 10000 # number of timesteps between checkpoints # PLOTFILES erf.plot_file_1 = plt # root name of plotfile erf.plot_int_1 = 1 # number of timesteps between plotfiles -erf.plot_vars_1 = density rhotheta rhoQt rhoQp x_velocity y_velocity z_velocity pressure theta temp qt qp qv qc qi +erf.plot_vars_1 = density rhotheta rhoQ1 rhoQ2 rhoQ3 x_velocity y_velocity z_velocity pressure theta temp qt qp qv qc qi # SOLVER CHOICE erf.use_gravity = true diff --git a/Exec/SuperCell/prob.H b/Exec/SuperCell/prob.H index 01cc85b99..77933842f 100644 --- a/Exec/SuperCell/prob.H +++ b/Exec/SuperCell/prob.H @@ -40,9 +40,6 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, - amrex::Array4 const& qv, - amrex::Array4 const& qc, - amrex::Array4 const& qi, amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/SuperCell/prob.cpp b/Exec/SuperCell/prob.cpp index e7186c3f6..a16792591 100644 --- a/Exec/SuperCell/prob.cpp +++ b/Exec/SuperCell/prob.cpp @@ -96,9 +96,6 @@ Problem::init_custom_pert( Array4 const& p_hse, Array4 const& /*z_nd*/, Array4 const& /*z_cc*/, - Array4 const& qv, - Array4 const& qc, - Array4 const& qi, GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, @@ -187,13 +184,6 @@ Problem::init_custom_pert( 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; -#endif } }); diff --git a/Source/BoundaryConditions/ERF_FillPatch.cpp b/Source/BoundaryConditions/ERF_FillPatch.cpp index 68270c55a..d7751995a 100644 --- a/Source/BoundaryConditions/ERF_FillPatch.cpp +++ b/Source/BoundaryConditions/ERF_FillPatch.cpp @@ -100,9 +100,11 @@ ERF::FillPatch (int lev, Real time, const Vector& mfs, bool fillset) // We call this even if init_type == real because this routine will fill the vertical bcs (*physbcs[lev])(mfs,icomp_cons,ncomp_cons,ngvect_cons,ngvect_vels,init_type,cons_only,BCVars::cons_bc,time); + /* // Update vars in the micro model if (solverChoice.moisture_type != MoistureType::None) micro.Update_Micro_Vars_Lev(lev, vars_new[lev][Vars::cons]); + */ } /* @@ -262,9 +264,11 @@ ERF::FillIntermediatePatch (int lev, Real time, if (!(cons_only && ncomp_cons == 1) && m_most && allow_most_bcs) m_most->impose_most_bcs(lev,mfs,eddyDiffs_lev[lev].get(),z_phys_nd[lev].get()); - // Update vars in the micro model - if (solverChoice.moisture_type != MoistureType::None) - micro.Update_Micro_Vars_Lev(lev, vars_new[lev][Vars::cons]); + /* + // Update vars in the micro model + if (solverChoice.moisture_type != MoistureType::None && (!cons_only && ncomp_cons > 2)) + micro.Update_Micro_Vars_Lev(lev, *mfs[Vars::cons]); + */ } /* diff --git a/Source/DataStructs/DataStruct.H b/Source/DataStructs/DataStruct.H index a3b31199a..2192f266a 100644 --- a/Source/DataStructs/DataStruct.H +++ b/Source/DataStructs/DataStruct.H @@ -85,6 +85,7 @@ struct SolverChoice { amrex::Abort("buoyancy_type must be 1, 2, 3 or 4"); } + /* // Set a different default for moist vs dry if (moisture_type != MoistureType::None) { if (moisture_type == MoistureType::FastEddy) { @@ -93,6 +94,7 @@ struct SolverChoice { buoyancy_type = 2; // uses Tprime } } + */ // Is the terrain static or moving? static std::string terrain_type_string = "Static"; diff --git a/Source/Diffusion/ComputeTurbulentViscosity.cpp b/Source/Diffusion/ComputeTurbulentViscosity.cpp index e74fa9f8f..2a0389277 100644 --- a/Source/Diffusion/ComputeTurbulentViscosity.cpp +++ b/Source/Diffusion/ComputeTurbulentViscosity.cpp @@ -170,8 +170,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; - // 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 + // EddyDiff mapping : Theta_h Scalar_h KE_h QKE_h Q1_h Q2_h Q3_h + Vector Factors = {inv_Pr_t, inv_Sc_t, inv_sigma_k, inv_sigma_k, inv_Sc_t, inv_Sc_t, inv_Sc_t}; // alpha = mu/Pr 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 b1d4065ab..cfbd7c078 100644 --- a/Source/Diffusion/DiffusionSrcForState_N.cpp +++ b/Source/Diffusion/DiffusionSrcForState_N.cpp @@ -100,6 +100,9 @@ DiffusionSrcForState_N (const amrex::Box& bx, const amrex::Box& domain, case PrimQ2_comp: alpha_eff[PrimQ2_comp] = diffChoice.alpha_C; break; + case PrimQ3_comp: + alpha_eff[PrimQ3_comp] = diffChoice.alpha_C; + break; default: alpha_eff[i] = 0.0; break; @@ -120,15 +123,18 @@ DiffusionSrcForState_N (const amrex::Box& bx, const amrex::Box& domain, case PrimQ2_comp: alpha_eff[PrimQ2_comp] = diffChoice.rhoAlpha_C; break; + case PrimQ3_comp: + alpha_eff[PrimQ3_comp] = diffChoice.rhoAlpha_C; + break; default: alpha_eff[i] = 0.0; break; } } } - 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}; + Vector eddy_diff_idx{EddyDiff::Theta_h, EddyDiff::KE_h, EddyDiff::QKE_h, EddyDiff::Scalar_h, EddyDiff::Q1_h, EddyDiff::Q2_h, EddyDiff::Q3_h}; + Vector eddy_diff_idy{EddyDiff::Theta_h, EddyDiff::KE_h, EddyDiff::QKE_h, EddyDiff::Scalar_h, EddyDiff::Q1_h, EddyDiff::Q2_h, EddyDiff::Q3_h}; + Vector eddy_diff_idz{EddyDiff::Theta_v, EddyDiff::KE_v, EddyDiff::QKE_v, EddyDiff::Scalar_v, EddyDiff::Q1_v, EddyDiff::Q2_v, EddyDiff::Q3_h}; // Device vectors Gpu::AsyncVector alpha_eff_d; @@ -138,7 +144,7 @@ DiffusionSrcForState_N (const amrex::Box& bx, const amrex::Box& domain, eddy_diff_idy_d.resize(eddy_diff_idy.size()); eddy_diff_idz_d.resize(eddy_diff_idz.size()); - Gpu::copy(Gpu::hostToDevice, alpha_eff.begin(), alpha_eff.end(), alpha_eff_d.begin()); + Gpu::copy(Gpu::hostToDevice, alpha_eff.begin() , alpha_eff.end() , alpha_eff_d.begin()); Gpu::copy(Gpu::hostToDevice, eddy_diff_idx.begin(), eddy_diff_idx.end(), eddy_diff_idx_d.begin()); Gpu::copy(Gpu::hostToDevice, eddy_diff_idy.begin(), eddy_diff_idy.end(), eddy_diff_idy_d.begin()); Gpu::copy(Gpu::hostToDevice, eddy_diff_idz.begin(), eddy_diff_idz.end(), eddy_diff_idz_d.begin()); diff --git a/Source/ERF.H b/Source/ERF.H index c6ff0ebe7..6b96dc255 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -646,7 +646,7 @@ private: amrex::Vector plot_var_names_1; amrex::Vector plot_var_names_2; const amrex::Vector cons_names {"density", "rhotheta", "rhoKE", "rhoQKE", "rhoadv_0", - "rhoQt", "rhoQp"}; + "rhoQ1", "rhoQ2", "rhoQ3"}; // Note that the order of variable names here must match the order in Derive.cpp const amrex::Vector derived_names {"soundspeed", "temp", "theta", "KE", "QKE", "scalar", diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 53a1f65c1..e368e0b4b 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -694,6 +694,12 @@ ERF::InitData () } } + // Update micro vars before first plot file + if (solverChoice.moisture_type != MoistureType::None) { + for (int lev = 0; lev <= finest_level; ++lev) micro.Update_Micro_Vars_Lev(lev, vars_new[lev][Vars::cons]); + } + + if (restart_chkfile.empty() && check_int > 0) { #ifdef ERF_USE_NETCDF diff --git a/Source/ERF_make_new_level.cpp b/Source/ERF_make_new_level.cpp index 7d9c77e7a..9ae9d1f1a 100644 --- a/Source/ERF_make_new_level.cpp +++ b/Source/ERF_make_new_level.cpp @@ -44,7 +44,8 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba, init_stuff(lev, ba, dm); - int ncomp_cons = (solverChoice.moisture_type == MoistureType::None) ? NVAR_max-2 : NVAR_max; + int n_qstate = micro.Get_Qstate_Size(); + int ncomp_cons = NVAR_max - (NMOIST_max - n_qstate); lev_new[Vars::cons].define(ba, dm, ncomp_cons, ngrow_state); lev_old[Vars::cons].define(ba, dm, ncomp_cons, ngrow_state); diff --git a/Source/IO/Checkpoint.cpp b/Source/IO/Checkpoint.cpp index dce2ee8a1..67c8886b5 100644 --- a/Source/IO/Checkpoint.cpp +++ b/Source/IO/Checkpoint.cpp @@ -128,20 +128,8 @@ ERF::WriteCheckpointFile () const MultiFab::Copy(zvel,vars_new[lev][Vars::zvel],0,0,1,0); VisMF::Write(zvel, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "ZFace")); - IntVect ng; - if (solverChoice.moisture_type != MoistureType::None) { - for (int mvar(0); mvarnGrowVect(); - int nvar = qmoist[lev][mvar]->nComp(); - MultiFab moist_vars(grids[lev],dmap[lev],nvar,ng); - MultiFab::Copy(moist_vars,*(qmoist[lev][mvar]),0,0,nvar,ng); - VisMF::Write(moist_vars, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "MoistVars")); - } - } - // Note that we write the ghost cells of the base state (unlike above) - ng = base_state[lev].nGrowVect(); + IntVect ng = base_state[lev].nGrowVect(); MultiFab base(grids[lev],dmap[lev],base_state[lev].nComp(),ng); MultiFab::Copy(base,base_state[lev],0,0,base.nComp(),ng); VisMF::Write(base, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "BaseState")); @@ -339,19 +327,8 @@ ERF::ReadCheckpointFile () VisMF::Read(zvel, amrex::MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "ZFace")); MultiFab::Copy(vars_new[lev][Vars::zvel],zvel,0,0,1,0); - IntVect ng; - if (solverChoice.moisture_type != MoistureType::None) { - for (int mvar(0); mvarnGrowVect(); - int nvar = qmoist[lev][mvar]->nComp(); - MultiFab moist_vars(grids[lev],dmap[lev],nvar,ng); - VisMF::Read(moist_vars, amrex::MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "MoistVars")); - MultiFab::Copy(*(qmoist[lev][mvar]),moist_vars,0,0,nvar,ng); - } - } - // Note that we read the ghost cells of the base state (unlike above) - ng = base_state[lev].nGrowVect(); + IntVect ng = base_state[lev].nGrowVect(); MultiFab base(grids[lev],dmap[lev],base_state[lev].nComp(),ng); VisMF::Read(base, amrex::MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "BaseState")); MultiFab::Copy(base_state[lev],base,0,0,base.nComp(),ng); diff --git a/Source/IO/Plotfile.cpp b/Source/IO/Plotfile.cpp index b237edd0e..b8327e13d 100644 --- a/Source/IO/Plotfile.cpp +++ b/Source/IO/Plotfile.cpp @@ -44,7 +44,8 @@ ERF::setPlotVariables (const std::string& pp_plot_var_names, Vector // since they may be in any order in the input list Vector tmp_plot_names; - int ncomp_cons = (solverChoice.moisture_type == MoistureType::None) ? NVAR_max-2 : NVAR_max; + int n_qstate = micro.Get_Qstate_Size(); + int ncomp_cons = NVAR_max - (NMOIST_max - n_qstate); for (int i = 0; i < ncomp_cons; ++i) { if ( containerHasElement(plot_var_names, cons_names[i]) ) { @@ -106,8 +107,7 @@ ERF::WritePlotFile (int which, Vector plot_var_names) int ncomp_cons = vars_new[0][Vars::cons].nComp(); - if (ncomp_mf == 0) - return; + if (ncomp_mf == 0) return; // We Fillpatch here because some of the derived quantities require derivatives // which require ghost cells to be filled. We do not need to call FillPatcher @@ -120,7 +120,7 @@ ERF::WritePlotFile (int which, Vector plot_var_names) } // Get qmoist pointers if using moisture - bool l_use_moisture = (solverChoice.moisture_type != MoistureType::None); + bool use_moisture = (solverChoice.moisture_type != MoistureType::None); for (int lev = 0; lev <= finest_level; ++lev) { for (int mvar(0); mvar plot_var_names) const Box& bx = mfi.tilebox(); const Array4& derdat = mf[lev].array(mfi); const Array4& S_arr = vars_new[lev][Vars::cons].const_array(mfi); - const Array4& qv_arr = (qmoist[lev][0]) ? qmoist[lev][0]->const_array(mfi) : - Array4 {}; ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - Real qv_for_p = (qv_arr) ? qv_arr(i,j,k) : 0; + Real qv_for_p = (use_moisture) ? S_arr(i,j,k,RhoQ1_comp)/S_arr(i,j,k,Rho_comp) : 0; const Real rhotheta = S_arr(i,j,k,RhoTheta_comp); derdat(i, j, k, mf_comp) = getPgivenRTh(rhotheta,qv_for_p); }); @@ -242,12 +240,10 @@ ERF::WritePlotFile (int which, Vector plot_var_names) const Array4& derdat = mf[lev].array(mfi); const Array4& p0_arr = p_hse.const_array(mfi); const Array4& S_arr = vars_new[lev][Vars::cons].const_array(mfi); - const Array4& qv_arr = (qmoist[lev][0]) ? qmoist[lev][0]->const_array(mfi) : - Array4 {}; ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - Real qv_for_p = (qv_arr) ? qv_arr(i,j,k) : 0; + Real qv_for_p = (use_moisture) ? S_arr(i,j,k,RhoQ1_comp)/S_arr(i,j,k,Rho_comp) : 0; const Real rhotheta = S_arr(i,j,k,RhoTheta_comp); derdat(i, j, k, mf_comp) = getPgivenRTh(rhotheta,qv_for_p) - p0_arr(i,j,k); }); @@ -610,68 +606,61 @@ ERF::WritePlotFile (int which, Vector plot_var_names) } // NOTE: Protect against accessing non-existent data - if (l_use_moisture) { + if (use_moisture) { int q_size = qmoist[lev].size(); - MultiFab qv_mf(*(qmoist[lev][0]), make_alias, 0, 1); - if (containerHasElement(plot_var_names, "qt") && (q_size >= 3)) + if (containerHasElement(plot_var_names, "qt") && (q_size >= 1)) { - MultiFab qc_mf(*(qmoist[lev][1]), make_alias, 0, 1); - MultiFab qi_mf(*(qmoist[lev][2]), make_alias, 0, 1); - MultiFab::Copy(mf[lev],qv_mf,0,mf_comp,1,0); - MultiFab::Add (mf[lev],qc_mf,0,mf_comp,1,0); - MultiFab::Add (mf[lev],qi_mf,0,mf_comp,1,0); + MultiFab qt_mf(*(qmoist[lev][0]), make_alias, 0, 1); + MultiFab::Copy(mf[lev],qt_mf,0,mf_comp,1,0); mf_comp += 1; } - if (containerHasElement(plot_var_names, "qp") && (q_size >= 6)) + if (containerHasElement(plot_var_names, "qv") && (q_size >= 2)) { - MultiFab qr_mf(*(qmoist[lev][3]), make_alias, 0, 1); - MultiFab qs_mf(*(qmoist[lev][4]), make_alias, 0, 1); - MultiFab qg_mf(*(qmoist[lev][5]), make_alias, 0, 1); - MultiFab::Copy(mf[lev],qr_mf,0,mf_comp,1,0); - MultiFab::Add (mf[lev],qs_mf,0,mf_comp,1,0); - MultiFab::Add (mf[lev],qg_mf,0,mf_comp,1,0); + MultiFab qv_mf(*(qmoist[lev][1]), make_alias, 0, 1); + MultiFab::Copy(mf[lev],qv_mf,0,mf_comp,1,0); mf_comp += 1; } - if (containerHasElement(plot_var_names, "qv")) + if (containerHasElement(plot_var_names, "qc") && (q_size >= 3)) { - MultiFab::Copy(mf[lev],qv_mf,0,mf_comp,1,0); + MultiFab qc_mf(*(qmoist[lev][2]), make_alias, 0, 1); + MultiFab::Copy(mf[lev],qc_mf,0,mf_comp,1,0); mf_comp += 1; } - if (containerHasElement(plot_var_names, "qc") && (q_size >= 2)) + if (containerHasElement(plot_var_names, "qi") && (q_size >= 4)) { - MultiFab qc_mf(*(qmoist[lev][1]), make_alias, 0, 1); - MultiFab::Copy(mf[lev],qc_mf,0,mf_comp,1,0); + MultiFab qi_mf(*(qmoist[lev][3]), make_alias, 0, 1); + MultiFab::Copy(mf[lev],qi_mf,0,mf_comp,1,0); mf_comp += 1; } - if (containerHasElement(plot_var_names, "qi") && (q_size >= 3)) + if (containerHasElement(plot_var_names, "qp") && (q_size >= 5)) { - MultiFab qi_mf(*(qmoist[lev][2]), make_alias, 0, 1); - MultiFab::Copy(mf[lev],qi_mf,0,mf_comp,1,0); + MultiFab qp_mf(*(qmoist[lev][4]), make_alias, 0, 1); + MultiFab::Copy(mf[lev],qp_mf,0,mf_comp,1,0); mf_comp += 1; } - if (containerHasElement(plot_var_names, "qrain") && (q_size >= 4)) + if (containerHasElement(plot_var_names, "qrain") && (q_size >= 6)) { - MultiFab qr_mf(*(qmoist[lev][3]), make_alias, 0, 1); + MultiFab qr_mf(*(qmoist[lev][5]), make_alias, 0, 1); MultiFab::Copy(mf[lev],qr_mf,0,mf_comp,1,0); mf_comp += 1; } - if (containerHasElement(plot_var_names, "qsnow") && (q_size >= 5)) + if (containerHasElement(plot_var_names, "qsnow") && (q_size >= 7)) { - MultiFab qs_mf(*(qmoist[lev][4]), make_alias, 0, 1); + MultiFab qs_mf(*(qmoist[lev][6]), make_alias, 0, 1); MultiFab::Copy(mf[lev],qs_mf,0,mf_comp,1,0); mf_comp += 1; } - if (containerHasElement(plot_var_names, "qgraup") && (q_size >= 6)) + if (containerHasElement(plot_var_names, "qgraup") && (q_size >= 8)) { - MultiFab qg_mf(*(qmoist[lev][5]), make_alias, 0, 1); + MultiFab qg_mf(*(qmoist[lev][7]), make_alias, 0, 1); MultiFab::Copy(mf[lev],qg_mf,0,mf_comp,1,0); mf_comp += 1; } diff --git a/Source/IndexDefines.H b/Source/IndexDefines.H index 09b299a6c..8ece30b94 100644 --- a/Source/IndexDefines.H +++ b/Source/IndexDefines.H @@ -16,11 +16,15 @@ #define RhoScalar_comp 4 #define RhoQ1_comp 5 #define RhoQ2_comp 6 +#define RhoQ3_comp 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 +#define NVAR_max 8 + +// This defines the maximum number of moisture vars +#define NMOIST_max 3 // Cell-centered primitive variables #define PrimTheta_comp (RhoTheta_comp -1) @@ -29,6 +33,7 @@ #define PrimScalar_comp (RhoScalar_comp-1) #define PrimQ1_comp (RhoQ1_comp-1) #define PrimQ2_comp (RhoQ2_comp-1) +#define PrimQ3_comp (RhoQ3_comp-1) // NOTE: we still use this indexing even if no moisture namespace BCVars { @@ -96,6 +101,7 @@ namespace EddyDiff { QKE_h, Q1_h, Q2_h, + Q3_h, Mom_v, Theta_v, Scalar_v, @@ -103,6 +109,7 @@ namespace EddyDiff { QKE_v, Q1_v, Q2_v, + Q3_v, PBL_lengthscale, NumDiffs }; diff --git a/Source/Initialization/ERF_init_custom.cpp b/Source/Initialization/ERF_init_custom.cpp index 94a832c7f..8758a164e 100644 --- a/Source/Initialization/ERF_init_custom.cpp +++ b/Source/Initialization/ERF_init_custom.cpp @@ -26,6 +26,8 @@ ERF::init_custom (int lev) { auto& lev_new = vars_new[lev]; + int n_qstate = micro.Get_Qstate_Size(); + MultiFab r_hse(base_state[lev], make_alias, 0, 1); // r_0 is first component MultiFab p_hse(base_state[lev], make_alias, 1, 1); // p_0 is second component @@ -41,13 +43,6 @@ ERF::init_custom (int lev) yvel_pert.setVal(0.); zvel_pert.setVal(0.); - MultiFab qv_pert(lev_new[Vars::cons].boxArray(), lev_new[Vars::cons].DistributionMap(), 1, lev_new[Vars::cons].nGrow()); - MultiFab qc_pert(lev_new[Vars::cons].boxArray(), lev_new[Vars::cons].DistributionMap(), 1, lev_new[Vars::cons].nGrow()); - MultiFab qi_pert(lev_new[Vars::cons].boxArray(), lev_new[Vars::cons].DistributionMap(), 1, lev_new[Vars::cons].nGrow()); - qv_pert.setVal(0.); - qc_pert.setVal(0.); - qi_pert.setVal(0.); - int fix_random_seed = 0; ParmParse pp("erf"); pp.query("fix_random_seed", fix_random_seed); // Note that the value of 1024UL is not significant -- the point here is just to set the @@ -82,12 +77,8 @@ ERF::init_custom (int lev) Array4 r_hse_arr = r_hse.array(mfi); Array4 p_hse_arr = p_hse.array(mfi); - const auto &qv_pert_arr = qv_pert.array(mfi); - const auto &qc_pert_arr = qc_pert.array(mfi); - const auto &qi_pert_arr = qi_pert.array(mfi); prob->init_custom_pert(bx, xbx, ybx, zbx, cons_pert_arr, xvel_pert_arr, yvel_pert_arr, zvel_pert_arr, r_hse_arr, p_hse_arr, z_nd_arr, z_cc_arr, - qv_pert_arr, qc_pert_arr, qi_pert_arr, geom[lev].data(), mf_m, mf_u, mf_v, solverChoice); } //mfi @@ -107,11 +98,7 @@ ERF::init_custom (int lev) if (solverChoice.moisture_type != MoistureType::None) { MultiFab::Add(lev_new[Vars::cons], cons_pert, RhoQ1_comp, RhoQ1_comp, 1, cons_pert.nGrow()); MultiFab::Add(lev_new[Vars::cons], cons_pert, RhoQ2_comp, RhoQ2_comp, 1, cons_pert.nGrow()); - - MultiFab::Add(*(qmoist[lev][0]) , qv_pert , 0, 0, 1, qv_pert.nGrow()); - MultiFab::Add(*(qmoist[lev][1]) , qc_pert , 0, 0, 1, qc_pert.nGrow()); - if (qmoist[lev].size() > 2) - MultiFab::Add(*(qmoist[lev][2]) , qi_pert , 0, 0, 1, qi_pert.nGrow()); + if (n_qstate >= 3) MultiFab::Add(lev_new[Vars::cons], cons_pert, RhoQ3_comp, RhoQ3_comp, 1, cons_pert.nGrow()); } MultiFab::Add(lev_new[Vars::xvel], xvel_pert, 0, 0, 1, xvel_pert.nGrowVect()); diff --git a/Source/Microphysics/FastEddy/FastEddy.H b/Source/Microphysics/FastEddy/FastEddy.H index 9e2113e51..9b8907f23 100644 --- a/Source/Microphysics/FastEddy/FastEddy.H +++ b/Source/Microphysics/FastEddy/FastEddy.H @@ -22,9 +22,11 @@ namespace MicVar_FE { // independent variables qv = 0, qc, + qt, + rho, // density theta, // liquid/ice water potential temperature tabs, // temperature - rho, // density + pres, // pressure NumVars }; } @@ -77,11 +79,6 @@ public: void Copy_Micro_to_State (amrex::MultiFab& cons_in) override; - // update ERF variables - void - Update (amrex::MultiFab& cons_in, - amrex::MultiFab& qmoist) override; - void Update_Micro_Vars (amrex::MultiFab& cons_in) override { @@ -115,9 +112,15 @@ public: int Qmoist_Size () override { return FastEddy::m_qmoist_size; } + int + Qstate_Size () { return FastEddy::m_qstate_size; } + private: - // Number of qmoist variables - int m_qmoist_size = 2; + // Number of qmoist variables (qt, qv, qc) + int m_qmoist_size = 3; + + // Number of qstate variables + int m_qstate_size = 2; // MicVar map (Qmoist indices -> MicVar enum) amrex::Vector MicVarMap; diff --git a/Source/Microphysics/FastEddy/FastEddy.cpp b/Source/Microphysics/FastEddy/FastEddy.cpp index 521c55a88..27210ba84 100644 --- a/Source/Microphysics/FastEddy/FastEddy.cpp +++ b/Source/Microphysics/FastEddy/FastEddy.cpp @@ -21,7 +21,7 @@ void FastEddy::AdvanceFE () auto qc_array = mic_fab_vars[MicVar_FE::qc]->array(mfi); auto tabs_array = mic_fab_vars[MicVar_FE::tabs]->array(mfi); auto theta_array = mic_fab_vars[MicVar_FE::theta]->array(mfi); - auto rho_array = mic_fab_vars[MicVar_FE::rho]->array(mfi); + auto pres_array = mic_fab_vars[MicVar_FE::pres]->array(mfi); const auto& box3d = mfi.tilebox(); @@ -36,7 +36,7 @@ void FastEddy::AdvanceFE () //------- Autoconversion/accretion Real dq_clwater_to_vapor, dq_vapor_to_clwater, qsat; - Real pressure = getPgivenRTh(rho_array(i,j,k)*theta_array(i,j,k),qv_array(i,j,k))/100.0; + Real pressure = pres_array(i,j,k); erf_qsatw(tabs_array(i,j,k), pressure, qsat); // If there is precipitating water (i.e. rain), and the cell is not saturated diff --git a/Source/Microphysics/FastEddy/Init_FE.cpp b/Source/Microphysics/FastEddy/Init_FE.cpp index 1e27d4e49..7ff34df97 100644 --- a/Source/Microphysics/FastEddy/Init_FE.cpp +++ b/Source/Microphysics/FastEddy/Init_FE.cpp @@ -1,4 +1,3 @@ - #include #include "Microphysics.H" #include "IndexDefines.H" @@ -29,7 +28,7 @@ void FastEddy::Init (const MultiFab& cons_in, m_gtoe = grids; MicVarMap.resize(m_qmoist_size); - MicVarMap = {MicVar_FE::qv, MicVar_FE::qc}; + MicVarMap = {MicVar_FE::qt, MicVar_FE::qv, MicVar_FE::qc}; // initialize microphysics variables for (auto ivar = 0; ivar < MicVar_FE::NumVars; ++ivar) { @@ -65,91 +64,14 @@ void FastEddy::Copy_State_to_Micro (const MultiFab& cons_in) auto states_array = cons_in.array(mfi); - auto qv_array = mic_fab_vars[MicVar_FE::qv]->array(mfi); - auto qc_array = mic_fab_vars[MicVar_FE::qc]->array(mfi); - auto rho_array = mic_fab_vars[MicVar_FE::rho]->array(mfi); - auto theta_array = mic_fab_vars[MicVar_FE::theta]->array(mfi); - auto temp_array = mic_fab_vars[MicVar_FE::tabs]->array(mfi); - - // Get pressure, theta, temperature, density, and qt, qp - amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - rho_array(i,j,k) = states_array(i,j,k,Rho_comp); - theta_array(i,j,k) = states_array(i,j,k,RhoTheta_comp)/states_array(i,j,k,Rho_comp); - qv_array(i,j,k) = states_array(i,j,k,RhoQ1_comp)/states_array(i,j,k,Rho_comp); - qc_array(i,j,k) = states_array(i,j,k,RhoQ2_comp)/states_array(i,j,k,Rho_comp); - temp_array(i,j,k) = getTgivenRandRTh(states_array(i,j,k,Rho_comp),states_array(i,j,k,RhoTheta_comp), qv_array(i,j,k)); - }); - } -} - - - -#if 0 -#include -#include "Microphysics.H" -#include "IndexDefines.H" -#include "PlaneAverage.H" -#include "EOS.H" -#include "TileNoZ.H" - -using namespace amrex; - -/** - * Initializes the Microphysics module. - * - * @param[in] cons_in Conserved variables input - * @param[in] qc_in Cloud variables input - * @param[in,out] qv_in Vapor variables input - * @param[in] qi_in Ice variables input - * @param[in] grids The boxes on which we will evolve the solution - * @param[in] geom Geometry associated with these MultiFabs and grids - * @param[in] dt_advance Timestep for the advance - */ -void FastEddy::Init (const MultiFab& cons_in, MultiFab& qmoist, - const BoxArray& grids, - const Geometry& geom, - const Real& dt_advance) -{ - m_geom = geom; - m_gtoe = grids; - - dt = dt_advance; - - // initialize microphysics variables - for (auto ivar = 0; ivar < MicVar_FE::NumVars; ++ivar) { - mic_fab_vars[ivar] = std::make_shared(cons_in.boxArray(), cons_in.DistributionMap(), 1, cons_in.nGrowVect()); - mic_fab_vars[ivar]->setVal(0.); - } - - // We must initialize to zero since we now need boundary values for the call to getP and we need all values filled - // The ghost cells of these arrays aren't filled in the boundary condition calls for the state - - qmoist.setVal(0.); - - for ( MFIter mfi(cons_in, TileNoZ()); mfi.isValid(); ++mfi) { - - const auto& box3d = mfi.tilebox(); - - const auto& lo = amrex::lbound(box3d); - const auto& hi = amrex::ubound(box3d); - - nlev = box3d.length(2); - zlo = lo.z; - zhi = hi.z; - } - - // Get the temperature, density, theta, qt and qp from input - for ( MFIter mfi(cons_in, false); mfi.isValid(); ++mfi) { - auto states_array = cons_in.array(mfi); - - auto qv_array = mic_fab_vars[MicVar_FE::qv]->array(mfi); - auto qc_array = mic_fab_vars[MicVar_FE::qc]->array(mfi); - auto rho_array = mic_fab_vars[MicVar_FE::rho]->array(mfi); - auto theta_array = mic_fab_vars[MicVar_FE::theta]->array(mfi); - auto temp_array = mic_fab_vars[MicVar_FE::tabs]->array(mfi); + auto qt_array = mic_fab_vars[MicVar_FE::qt]->array(mfi); + auto qv_array = mic_fab_vars[MicVar_FE::qv]->array(mfi); + auto qc_array = mic_fab_vars[MicVar_FE::qc]->array(mfi); - const auto& box3d = mfi.tilebox(); + auto rho_array = mic_fab_vars[MicVar_FE::rho]->array(mfi); + auto theta_array = mic_fab_vars[MicVar_FE::theta]->array(mfi); + auto tabs_array = mic_fab_vars[MicVar_FE::tabs]->array(mfi); + auto pres_array = mic_fab_vars[MicVar_FE::pres]->array(mfi); // Get pressure, theta, temperature, density, and qt, qp amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) @@ -158,8 +80,12 @@ void FastEddy::Init (const MultiFab& cons_in, MultiFab& qmoist, theta_array(i,j,k) = states_array(i,j,k,RhoTheta_comp)/states_array(i,j,k,Rho_comp); qv_array(i,j,k) = states_array(i,j,k,RhoQ1_comp)/states_array(i,j,k,Rho_comp); qc_array(i,j,k) = states_array(i,j,k,RhoQ2_comp)/states_array(i,j,k,Rho_comp); - temp_array(i,j,k) = getTgivenRandRTh(states_array(i,j,k,Rho_comp),states_array(i,j,k,RhoTheta_comp), qv_array(i,j,k)); + qt_array(i,j,k) = qv_array(i,j,k) + qc_array(i,j,k); + + tabs_array(i,j,k) = getTgivenRandRTh(states_array(i,j,k,Rho_comp), + states_array(i,j,k,RhoTheta_comp), + qv_array(i,j,k)); + pres_array(i,j,k) = getPgivenRTh(states_array(i,j,k,RhoTheta_comp), qv_array(i,j,k))/100.; }); } } -#endif diff --git a/Source/Microphysics/FastEddy/Update_FE.cpp b/Source/Microphysics/FastEddy/Update_FE.cpp index ba3d30055..7f405acf4 100644 --- a/Source/Microphysics/FastEddy/Update_FE.cpp +++ b/Source/Microphysics/FastEddy/Update_FE.cpp @@ -3,40 +3,6 @@ #include "IndexDefines.H" #include "TileNoZ.H" -/** - * Updates conserved and microphysics variables in the provided MultiFabs from - * the internal MultiFabs that store Microphysics module data. - * - * @param[out] cons Conserved variables - * @param[out] qmoist: qv, qc, qi, qr, qs, qg - */ -void FastEddy::Update (amrex::MultiFab& cons, - amrex::MultiFab& qmoist) -{ - // Get the temperature, density, theta, qt and qp from input - for ( amrex::MFIter mfi(cons,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { - auto states_arr = cons.array(mfi); - - auto rho_arr = mic_fab_vars[MicVar_FE::rho]->array(mfi); - auto theta_arr = mic_fab_vars[MicVar_FE::theta]->array(mfi); - auto qv_arr = mic_fab_vars[MicVar_FE::qv]->array(mfi); - auto qc_arr = mic_fab_vars[MicVar_FE::qc]->array(mfi); - const auto& box3d = mfi.tilebox(); - - // get potential total density, temperature, qt, qp - amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - states_arr(i,j,k,RhoTheta_comp) = rho_arr(i,j,k)*theta_arr(i,j,k); - states_arr(i,j,k,RhoQ1_comp) = rho_arr(i,j,k)*qv_arr(i,j,k); - states_arr(i,j,k,RhoQ2_comp) = rho_arr(i,j,k)*qc_arr(i,j,k); - }); - } - - // Fill interior ghost cells and periodic boundaries - cons.FillBoundary(m_geom.periodicity()); - qmoist.FillBoundary(m_geom.periodicity()); -} - /** * Updates conserved and microphysics variables in the provided MultiFabs from * the internal MultiFabs that store Microphysics module data. diff --git a/Source/Microphysics/Kessler/Diagnose_Kessler.cpp b/Source/Microphysics/Kessler/Diagnose_Kessler.cpp index 55b613ca3..e63a4a84d 100644 --- a/Source/Microphysics/Kessler/Diagnose_Kessler.cpp +++ b/Source/Microphysics/Kessler/Diagnose_Kessler.cpp @@ -31,44 +31,14 @@ void Kessler::Diagnose () amrex::ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - qv_array(i,j,k) = std::max(0.0,qt_array(i,j,k) - qn_array(i,j,k)); + qt_array(i,j,k) = qv_array(i,j,k) + qn_array(i,j,k); amrex::Real omn = 1.0; qcl_array(i,j,k) = qn_array(i,j,k)*omn; qci_array(i,j,k) = qn_array(i,j,k)*(1.0-omn); - amrex::Real omp = 1.0;; + amrex::Real omp = 1.0; qpl_array(i,j,k) = qp_array(i,j,k)*omp; qpi_array(i,j,k) = qp_array(i,j,k)*(1.0-omp); }); } } - -/** - * Computes diagnostic quantities like temperature and pressure - * from the existing Microphysics variables. - */ -void Kessler::Compute_Tabs_and_Pres () -{ - auto rho = mic_fab_vars[MicVar_Kess::rho]; - auto theta = mic_fab_vars[MicVar_Kess::theta]; - auto qv = mic_fab_vars[MicVar_Kess::qv]; - auto tabs = mic_fab_vars[MicVar_Kess::tabs]; - auto pres = mic_fab_vars[MicVar_Kess::pres]; - - for ( amrex::MFIter mfi(*tabs, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { - const auto& box3d = mfi.tilebox(); - - auto rho_array = rho->array(mfi); - auto theta_array = theta->array(mfi); - auto qv_array = qv->array(mfi); - auto tabs_array = tabs->array(mfi); - auto pres_array = pres->array(mfi); - - amrex::ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept - { - amrex::Real RhoTheta = rho_array(i,j,k)*theta_array(i,j,k); - tabs_array(i,j,k) = getTgivenRandRTh(rho_array(i,j,k), RhoTheta, qv_array(i,j,k)); - pres_array(i,j,k) = getPgivenRTh(RhoTheta, qv_array(i,j,k))/100.; - }); - } -} diff --git a/Source/Microphysics/Kessler/Init_Kessler.cpp b/Source/Microphysics/Kessler/Init_Kessler.cpp index 56f8c6be4..6f626550f 100644 --- a/Source/Microphysics/Kessler/Init_Kessler.cpp +++ b/Source/Microphysics/Kessler/Init_Kessler.cpp @@ -1,4 +1,3 @@ - #include #include "Microphysics.H" #include "IndexDefines.H" @@ -30,8 +29,8 @@ void Kessler::Init (const MultiFab& cons_in, m_gtoe = grids; MicVarMap.resize(m_qmoist_size); - MicVarMap = {MicVar_Kess::qv, MicVar_Kess::qcl, MicVar_Kess::qci, - MicVar_Kess::qr, MicVar_Kess::qpi, MicVar_Kess::qg}; + MicVarMap = {MicVar_Kess::qt, MicVar_Kess::qv, MicVar_Kess::qcl, MicVar_Kess::qci, + MicVar_Kess::qp, MicVar_Kess::qpl, MicVar_Kess::qpi}; // initialize microphysics variables for (auto ivar = 0; ivar < MicVar_Kess::NumVars; ++ivar) { @@ -66,278 +65,32 @@ void Kessler::Copy_State_to_Micro (const MultiFab& cons_in) auto states_array = cons_in.array(mfi); - auto qt_array = mic_fab_vars[MicVar_Kess::qt]->array(mfi); - auto qp_array = mic_fab_vars[MicVar_Kess::qp]->array(mfi); - auto qn_array = mic_fab_vars[MicVar_Kess::qn]->array(mfi); - auto qc_array = mic_fab_vars[MicVar_Kess::qcl]->array(mfi); - auto qi_array = mic_fab_vars[MicVar_Kess::qci]->array(mfi); - auto rho_array = mic_fab_vars[MicVar_Kess::rho]->array(mfi); - auto theta_array = mic_fab_vars[MicVar_Kess::theta]->array(mfi); + auto qv_array = mic_fab_vars[MicVar_Kess::qv]->array(mfi); + auto qc_array = mic_fab_vars[MicVar_Kess::qcl]->array(mfi); + auto qi_array = mic_fab_vars[MicVar_Kess::qci]->array(mfi); + auto qn_array = mic_fab_vars[MicVar_Kess::qn]->array(mfi); + auto qp_array = mic_fab_vars[MicVar_Kess::qp]->array(mfi); + + auto rho_array = mic_fab_vars[MicVar_Kess::rho]->array(mfi); + auto theta_array = mic_fab_vars[MicVar_Kess::theta]->array(mfi); + auto tabs_array = mic_fab_vars[MicVar::tabs]->array(mfi); + auto pres_array = mic_fab_vars[MicVar::pres]->array(mfi); // Get pressure, theta, temperature, density, and qt, qp amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { rho_array(i,j,k) = states_array(i,j,k,Rho_comp); theta_array(i,j,k) = states_array(i,j,k,RhoTheta_comp)/states_array(i,j,k,Rho_comp); - qt_array(i,j,k) = states_array(i,j,k,RhoQ1_comp)/states_array(i,j,k,Rho_comp); - qp_array(i,j,k) = states_array(i,j,k,RhoQ2_comp)/states_array(i,j,k,Rho_comp); + qv_array(i,j,k) = states_array(i,j,k,RhoQ1_comp)/states_array(i,j,k,Rho_comp); + qc_array(i,j,k) = states_array(i,j,k,RhoQ2_comp)/states_array(i,j,k,Rho_comp); + qp_array(i,j,k) = states_array(i,j,k,RhoQ3_comp)/states_array(i,j,k,Rho_comp); qn_array(i,j,k) = qc_array(i,j,k) + qi_array(i,j,k); - //qv_array(i,j,k) = qv_array_from_moist(i,j,k); - //temp_array(i,j,k) = getTgivenRandRTh(states_array(i,j,k,Rho_comp),states_array(i,j,k,RhoTheta_comp), qv_array(i,j,k)); - //pres_array(i,j,k) = getPgivenRTh(states_array(i,j,k,RhoTheta_comp))/100.; + + tabs_array(i,j,k) = getTgivenRandRTh(states_array(i,j,k,Rho_comp), + states_array(i,j,k,RhoTheta_comp), + qv_array(i,j,k)); + pres_array(i,j,k) = getPgivenRTh(states_array(i,j,k,RhoTheta_comp), qv_array(i,j,k))/100.; }); } } -#if 0 -/** - * Initializes the Microphysics module. - * - * @param[in] cons_in Conserved variables input - * @param[in] qc_in Cloud variables input - * @param[in,out] qv_in Vapor variables input - * @param[in] qi_in Ice variables input - * @param[in] grids The boxes on which we will evolve the solution - * @param[in] geom Geometry associated with these MultiFabs and grids - * @param[in] dt_advance Timestep for the advance - */ -void Kessler::Init (const MultiFab& cons_in, - MultiFab& qmoist, - const BoxArray& grids, - const Geometry& geom, - const Real& dt_advance) - { - m_geom = geom; - m_gtoe = grids; - - auto dz = m_geom.CellSize(2); - auto lowz = m_geom.ProbLo(2); - - dt = dt_advance; - - // initialize microphysics variables - for (auto ivar = 0; ivar < MicVar_Kess::NumVars; ++ivar) { - mic_fab_vars[ivar] = std::make_shared(cons_in.boxArray(), cons_in.DistributionMap(), 1, cons_in.nGrowVect()); - mic_fab_vars[ivar]->setVal(0.); - } - - // We must initialize to zero since we now need boundary values for the call to getP and we need all values filled - // The ghost cells of these arrays aren't filled in the boundary condition calls for the state - - //qmoist.setVal(0.); - - for ( MFIter mfi(cons_in, TileNoZ()); mfi.isValid(); ++mfi) { - - const auto& box3d = mfi.tilebox(); - - const auto& lo = amrex::lbound(box3d); - const auto& hi = amrex::ubound(box3d); - - nlev = box3d.length(2); - zlo = lo.z; - zhi = hi.z; - - - // parameters - accrrc.resize({zlo}, {zhi}); - accrsi.resize({zlo}, {zhi}); - accrsc.resize({zlo}, {zhi}); - coefice.resize({zlo}, {zhi}); - evaps1.resize({zlo}, {zhi}); - evaps2.resize({zlo}, {zhi}); - accrgi.resize({zlo}, {zhi}); - accrgc.resize({zlo}, {zhi}); - evapg1.resize({zlo}, {zhi}); - evapg2.resize({zlo}, {zhi}); - evapr1.resize({zlo}, {zhi}); - evapr2.resize({zlo}, {zhi}); - - // data (input) - rho1d.resize({zlo}, {zhi}); - pres1d.resize({zlo}, {zhi}); - tabs1d.resize({zlo}, {zhi}); - gamaz.resize({zlo}, {zhi}); - zmid.resize({zlo}, {zhi}); - - // output - qifall.resize({zlo}, {zhi}); - tlatqi.resize({zlo}, {zhi}); - - } - - - auto accrrc_t = accrrc.table(); - auto accrsi_t = accrsi.table(); - auto accrsc_t = accrsc.table(); - auto coefice_t = coefice.table(); - auto evaps1_t = evaps1.table(); - auto evaps2_t = evaps2.table(); - auto accrgi_t = accrgi.table(); - auto accrgc_t = accrgc.table(); - auto evapg1_t = evapg1.table(); - auto evapg2_t = evapg2.table(); - auto evapr1_t = evapr1.table(); - auto evapr2_t = evapr2.table(); - - auto rho1d_t = rho1d.table(); - auto pres1d_t = pres1d.table(); - auto tabs1d_t = tabs1d.table(); - - auto gamaz_t = gamaz.table(); - auto zmid_t = zmid.table(); - - Real gam3 = erf_gammafff(3.0 ); - Real gamr1 = erf_gammafff(3.0+b_rain ); - Real gamr2 = erf_gammafff((5.0+b_rain)/2.0); - // Real gamr3 = erf_gammafff(4.0+b_rain ); - Real gams1 = erf_gammafff(3.0+b_snow ); - Real gams2 = erf_gammafff((5.0+b_snow)/2.0); - // Real gams3 = erf_gammafff(4.0+b_snow ); - Real gamg1 = erf_gammafff(3.0+b_grau ); - Real gamg2 = erf_gammafff((5.0+b_grau)/2.0); - // Real gamg3 = erf_gammafff(4.0+b_grau ); - - - amrex::MultiFab qv(qmoist, amrex::make_alias, 0, 1); - amrex::MultiFab qc(qmoist, amrex::make_alias, 1, 1); - amrex::MultiFab qi(qmoist, amrex::make_alias, 2, 1); - - // Get the temperature, density, theta, qt and qp from input - for ( MFIter mfi(cons_in, false); mfi.isValid(); ++mfi) { - auto states_array = cons_in.array(mfi); - auto qv_array_from_moist = qv.array(mfi); - auto qc_array = qc.array(mfi); - auto qi_array = qi.array(mfi); - - auto qt_array = mic_fab_vars[MicVar_Kess::qt]->array(mfi); - auto qp_array = mic_fab_vars[MicVar_Kess::qp]->array(mfi); - auto qn_array = mic_fab_vars[MicVar_Kess::qn]->array(mfi); - auto qv_array = mic_fab_vars[MicVar_Kess::qv]->array(mfi); - auto rho_array = mic_fab_vars[MicVar_Kess::rho]->array(mfi); - auto theta_array = mic_fab_vars[MicVar_Kess::theta]->array(mfi); - auto temp_array = mic_fab_vars[MicVar_Kess::tabs]->array(mfi); - auto pres_array = mic_fab_vars[MicVar_Kess::pres]->array(mfi); - - const auto& box3d = mfi.tilebox(); - - // Get pressure, theta, temperature, density, and qt, qp - amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - rho_array(i,j,k) = states_array(i,j,k,Rho_comp); - theta_array(i,j,k) = states_array(i,j,k,RhoTheta_comp)/states_array(i,j,k,Rho_comp); - qt_array(i,j,k) = states_array(i,j,k,RhoQ1_comp)/states_array(i,j,k,Rho_comp); - qp_array(i,j,k) = states_array(i,j,k,RhoQ2_comp)/states_array(i,j,k,Rho_comp); - qn_array(i,j,k) = qc_array(i,j,k) + qi_array(i,j,k); - qv_array(i,j,k) = qv_array_from_moist(i,j,k); - temp_array(i,j,k) = getTgivenRandRTh(states_array(i,j,k,Rho_comp),states_array(i,j,k,RhoTheta_comp), qv_array(i,j,k)); - pres_array(i,j,k) = getPgivenRTh(states_array(i,j,k,RhoTheta_comp))/100.; - }); - } - - - // calculate the plane average variables - PlaneAverage cons_ave(&cons_in, m_geom, m_axis); - cons_ave.compute_averages(ZDir(), cons_ave.field()); - - // get host variable rho, and rhotheta - int ncell = cons_ave.ncell_line(); - - Gpu::HostVector rho_h(ncell), rhotheta_h(ncell); - cons_ave.line_average(Rho_comp, rho_h); - cons_ave.line_average(RhoTheta_comp, rhotheta_h); - - // copy data to device - Gpu::DeviceVector rho_d(ncell), rhotheta_d(ncell); - Gpu::copyAsync(Gpu::hostToDevice, rho_h.begin(), rho_h.end(), rho_d.begin()); - Gpu::copyAsync(Gpu::hostToDevice, rhotheta_h.begin(), rhotheta_h.end(), rhotheta_d.begin()); - Gpu::streamSynchronize(); - - Real* rho_dptr = rho_d.data(); - Real* rhotheta_dptr = rhotheta_d.data(); - - Real gOcp = m_gOcp; - - amrex::ParallelFor(nlev, [=] AMREX_GPU_DEVICE (int k) noexcept - { - Real pressure = getPgivenRTh(rhotheta_dptr[k]); - rho1d_t(k) = rho_dptr[k]; - pres1d_t(k) = pressure/100.; - tabs1d_t(k) = getTgivenRandRTh(rho_dptr[k],rhotheta_dptr[k]); - zmid_t(k) = lowz + (k+0.5)*dz; - gamaz_t(k) = gOcp*zmid_t(k); - }); - - // This fills qv - Diagnose(); - -#if 0 - amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int k, int j, int i) - { - fluxbmk(l,j,i) = 0.0; - fluxtmk(l,j,i) = 0.0; - }); - - amrex::ParallelFor( box2d, [=] AMREX_GPU_DEVICE (int k, int l, int i0) - { - mkwle (l,k) = 0.0; - mkwsb (l,k) = 0.0; - mkadv (l,k) = 0.0; - mkdiff(l,k) = 0.0; - }); -#endif - - if(round(gam3) != 2) { - std::cout << "cannot compute gamma-function in Microphysics::Init" << std::endl; - std::exit(-1); - } - - amrex::ParallelFor(nlev, [=] AMREX_GPU_DEVICE (int k) noexcept - { - Real pratio = sqrt(1.29 / rho1d_t(k)); - Real rrr1=393.0/(tabs1d_t(k)+120.0)*std::pow((tabs1d_t(k)/273.0),1.5); - Real rrr2=std::pow((tabs1d_t(k)/273.0),1.94)*(1000.0/pres1d_t(k)); - Real estw = 100.0*erf_esatw(tabs1d_t(k)); - Real esti = 100.0*erf_esati(tabs1d_t(k)); - - // accretion by snow: - Real coef1 = 0.25 * PI * nzeros * a_snow * gams1 * pratio/pow((PI * rhos * nzeros/rho1d_t(k) ) , ((3.0+b_snow)/4.0)); - Real coef2 = exp(0.025*(tabs1d_t(k) - 273.15)); - accrsi_t(k) = coef1 * coef2 * esicoef; - accrsc_t(k) = coef1 * esccoef; - coefice_t(k) = coef2; - - // evaporation of snow: - coef1 =(lsub/(tabs1d_t(k)*rv)-1.0)*lsub/(therco*rrr1*tabs1d_t(k)); - coef2 = rv*tabs1d_t(k)/(diffelq*rrr2*esti); - evaps1_t(k) = 0.65*4.0*nzeros/sqrt(PI*rhos*nzeros)/(coef1+coef2)/sqrt(rho1d_t(k)); - evaps2_t(k) = 0.49*4.0*nzeros*gams2*sqrt(a_snow/(muelq*rrr1))/pow((PI*rhos*nzeros) , ((5.0+b_snow)/8.0)) / - (coef1+coef2) * pow(rho1d_t(k) , ((1.0+b_snow)/8.0))*sqrt(pratio); - - // accretion by graupel: - coef1 = 0.25*PI*nzerog*a_grau*gamg1*pratio/pow((PI*rhog*nzerog/rho1d_t(k)) , ((3.0+b_grau)/4.0)); - coef2 = exp(0.025*(tabs1d_t(k) - 273.15)); - accrgi_t(k) = coef1 * coef2 * egicoef; - accrgc_t(k) = coef1 * egccoef; - - // evaporation of graupel: - coef1 =(lsub/(tabs1d_t(k)*rv)-1.0)*lsub/(therco*rrr1*tabs1d_t(k)); - coef2 = rv*tabs1d_t(k)/(diffelq*rrr2*esti); - evapg1_t(k) = 0.65*4.0*nzerog/sqrt(PI*rhog*nzerog)/(coef1+coef2)/sqrt(rho1d_t(k)); - evapg2_t(k) = 0.49*4.0*nzerog*gamg2*sqrt(a_grau/(muelq*rrr1))/pow((PI * rhog * nzerog) , ((5.0+b_grau)/8.0)) / - (coef1+coef2) * pow(rho1d_t(k) , ((1.0+b_grau)/8.0))*sqrt(pratio); - - // accretion by rain: - accrrc_t(k)= 0.25 * PI * nzeror * a_rain * gamr1 * pratio/pow((PI * rhor * nzeror / rho1d_t(k)) , ((3+b_rain)/4.))* erccoef; - - // evaporation of rain: - coef1 =(lcond/(tabs1d_t(k)*rv)-1.0)*lcond/(therco*rrr1*tabs1d_t(k)); - coef2 = rv*tabs1d_t(k)/(diffelq * rrr2 * estw); - evapr1_t(k) = 0.78 * 2.0 * PI * nzeror / - sqrt(PI * rhor * nzeror) / (coef1+coef2) / sqrt(rho1d_t(k)); - evapr2_t(k) = 0.31 * 2.0 * PI * nzeror * gamr2 * 0.89 * sqrt(a_rain/(muelq*rrr1))/ - pow((PI * rhor * nzeror) , ((5.0+b_rain)/8.0)) / - (coef1+coef2) * pow(rho1d_t(k) , ((1.0+b_rain)/8.0))*sqrt(pratio); - }); - -} -#endif diff --git a/Source/Microphysics/Kessler/Kessler.H b/Source/Microphysics/Kessler/Kessler.H index 00ade80ef..26c7c147d 100644 --- a/Source/Microphysics/Kessler/Kessler.H +++ b/Source/Microphysics/Kessler/Kessler.H @@ -95,18 +95,12 @@ public: void Copy_Micro_to_State (amrex::MultiFab& cons_in) override; - // update ERF variables - void - Update (amrex::MultiFab& cons_in, - amrex::MultiFab& qmoist) override; - // update micro vars void Update_Micro_Vars (amrex::MultiFab& cons_in) override { this->Copy_State_to_Micro(cons_in); this->Diagnose(); - this->Compute_Tabs_and_Pres(); } // update state vars @@ -124,7 +118,6 @@ public: this->AdvanceKessler(); this->Diagnose(); - this->Compute_Tabs_and_Pres(); } amrex::MultiFab* @@ -134,15 +127,18 @@ public: return mic_fab_vars[MicVarMap[varIdx]].get(); } - void - Compute_Tabs_and_Pres (); - int Qmoist_Size () override { return Kessler::m_qmoist_size; } + int + Qstate_Size () { return Kessler::m_qstate_size; } + private: - // Number of qmoist variables - int m_qmoist_size = 6; + // Number of qmoist variables (qt, qv, qcl, qci, qp, qpl, qpi) + int m_qmoist_size = 7; + + // Number of qstate variables + int m_qstate_size = 3; // MicVar map (Qmoist indices -> MicVar enum) amrex::Vector MicVarMap; @@ -173,40 +169,5 @@ private: // independent variables amrex::Array mic_fab_vars; - -#if 0 - // microphysics parameters/coefficients - amrex::TableData accrrc; - amrex::TableData accrsi; - amrex::TableData accrsc; - amrex::TableData coefice; - amrex::TableData evaps1; - amrex::TableData evaps2; - amrex::TableData accrgi; - amrex::TableData accrgc; - amrex::TableData evapg1; - amrex::TableData evapg2; - amrex::TableData evapr1; - amrex::TableData evapr2; - - // vertical plane average data - amrex::TableData rho1d; - amrex::TableData pres1d; - amrex::TableData tabs1d; - amrex::TableData qt1d; - amrex::TableData qv1d; - amrex::TableData qn1d; - - // independent variables - amrex::Array mic_fab_vars; - - amrex::TableData gamaz; - amrex::TableData zmid; // mid value of vertical coordinate in physical domain - - // data (output) - amrex::TableData qifall; - amrex::TableData tlatqi; -#endif - }; #endif diff --git a/Source/Microphysics/Kessler/Kessler.cpp b/Source/Microphysics/Kessler/Kessler.cpp index 4daee2c60..c7077fdd1 100644 --- a/Source/Microphysics/Kessler/Kessler.cpp +++ b/Source/Microphysics/Kessler/Kessler.cpp @@ -102,7 +102,7 @@ void Kessler::AdvanceKessler () //------- Autoconversion/accretion Real qcc, autor, accrr, dq_clwater_to_rain, dq_rain_to_vapor, dq_clwater_to_vapor, dq_vapor_to_clwater, qsat; - Real pressure = pres_array(i,j,k); //getPgivenRTh(rho_array(i,j,k)*theta_array(i,j,k),qv_array(i,j,k))/100.0; + Real pressure = pres_array(i,j,k); erf_qsatw(tabs_array(i,j,k), pressure, qsat); // If there is precipitating water (i.e. rain), and the cell is not saturated diff --git a/Source/Microphysics/Kessler/Update_Kessler.cpp b/Source/Microphysics/Kessler/Update_Kessler.cpp index cc085f9e0..990425e4e 100644 --- a/Source/Microphysics/Kessler/Update_Kessler.cpp +++ b/Source/Microphysics/Kessler/Update_Kessler.cpp @@ -3,60 +3,6 @@ #include "IndexDefines.H" #include "TileNoZ.H" -/** - * Updates conserved and microphysics variables in the provided MultiFabs from - * the internal MultiFabs that store Microphysics module data. - * - * @param[out] cons Conserved variables - * @param[out] qmoist: qv, qc, qi, qr, qs, qg - */ -void Kessler::Update (amrex::MultiFab& cons, - amrex::MultiFab& qmoist) -{ - // copy multifab data to qc, qv, and qi - amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar_Kess::qv], 0, 0, 1, mic_fab_vars[MicVar_Kess::qv]->nGrowVect()); // vapor - amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar_Kess::qcl], 0, 1, 1, mic_fab_vars[MicVar_Kess::qcl]->nGrowVect()); // cloud water - amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar_Kess::qci], 0, 2, 1, mic_fab_vars[MicVar_Kess::qci]->nGrowVect()); // cloud ice - amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar_Kess::qpl], 0, 3, 1, mic_fab_vars[MicVar_Kess::qpl]->nGrowVect()); // rain - amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar_Kess::qpi], 0, 4, 1, mic_fab_vars[MicVar_Kess::qpi]->nGrowVect()); // snow - - // Don't need to copy this since it is filled below - // amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar_Kess::qpi], 0, 5, 1, mic_fab_vars[MicVar_Kess::qci]->nGrowVect()); // graupel - - amrex::MultiFab qgraup_mf(qmoist, amrex::make_alias, 5, 1); - - // Get the temperature, density, theta, qt and qp from input - for ( amrex::MFIter mfi(cons,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { - auto states_arr = cons.array(mfi); - - auto rho_arr = mic_fab_vars[MicVar_Kess::rho]->array(mfi); - auto theta_arr = mic_fab_vars[MicVar_Kess::theta]->array(mfi); - auto qt_arr = mic_fab_vars[MicVar_Kess::qt]->array(mfi); - auto qp_arr = mic_fab_vars[MicVar_Kess::qp]->array(mfi); - - auto qgraup_arr= qgraup_mf.array(mfi); - - const auto& box3d = mfi.tilebox(); - - // get potential total density, temperature, qt, qp - amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - states_arr(i,j,k,Rho_comp) = rho_arr(i,j,k); - states_arr(i,j,k,RhoTheta_comp) = rho_arr(i,j,k)*theta_arr(i,j,k); - states_arr(i,j,k,RhoQ1_comp) = rho_arr(i,j,k)*qt_arr(i,j,k); - states_arr(i,j,k,RhoQ2_comp) = rho_arr(i,j,k)*qp_arr(i,j,k); - - // Graupel == precip total - rain - snow (but must be >= 0) - qgraup_arr(i,j,k) = 0.0;// - }); - } - - // Fill interior ghost cells and periodic boundaries - cons.FillBoundary(m_geom.periodicity()); - qmoist.FillBoundary(m_geom.periodicity()); -} - - /** * Updates conserved and microphysics variables in the provided MultiFabs from * the internal MultiFabs that store Microphysics module data. @@ -74,16 +20,17 @@ void Kessler::Copy_Micro_to_State (amrex::MultiFab& cons) auto rho_arr = mic_fab_vars[MicVar_Kess::rho]->array(mfi); auto theta_arr = mic_fab_vars[MicVar_Kess::theta]->array(mfi); - auto qt_arr = mic_fab_vars[MicVar_Kess::qt]->array(mfi); + auto qv_arr = mic_fab_vars[MicVar_Kess::qv]->array(mfi); + auto qc_arr = mic_fab_vars[MicVar_Kess::qcl]->array(mfi); auto qp_arr = mic_fab_vars[MicVar_Kess::qp]->array(mfi); // get potential total density, temperature, qt, qp amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - //states_arr(i,j,k,Rho_comp) = rho_arr(i,j,k); states_arr(i,j,k,RhoTheta_comp) = rho_arr(i,j,k)*theta_arr(i,j,k); - states_arr(i,j,k,RhoQ1_comp) = rho_arr(i,j,k)*qt_arr(i,j,k); - states_arr(i,j,k,RhoQ2_comp) = rho_arr(i,j,k)*qp_arr(i,j,k); + states_arr(i,j,k,RhoQ1_comp) = rho_arr(i,j,k)*qv_arr(i,j,k); + states_arr(i,j,k,RhoQ2_comp) = rho_arr(i,j,k)*qc_arr(i,j,k); + states_arr(i,j,k,RhoQ3_comp) = rho_arr(i,j,k)*qp_arr(i,j,k); }); } diff --git a/Source/Microphysics/Microphysics.H b/Source/Microphysics/Microphysics.H index 4113fe8db..cd9b4514b 100644 --- a/Source/Microphysics/Microphysics.H +++ b/Source/Microphysics/Microphysics.H @@ -55,14 +55,6 @@ public: m_moist_model[lev]->Diagnose(); } - void - Update (const int& lev, - amrex::MultiFab& cons_in, - amrex::MultiFab& qmoist) - { - m_moist_model[lev]->Update(cons_in, qmoist); - } - void Update_Micro_Vars_Lev (const int& lev, amrex::MultiFab& cons_in) { @@ -81,6 +73,9 @@ public: int Get_Qmoist_Size () { return m_moist_model[0]->Qmoist_Size(); } + int + Get_Qstate_Size () { return m_moist_model[0]->Qstate_Size(); } + private: amrex::Vector> m_moist_model; }; diff --git a/Source/Microphysics/Null/NullMoist.H b/Source/Microphysics/Null/NullMoist.H index 5dfb0836e..681513bbc 100644 --- a/Source/Microphysics/Null/NullMoist.H +++ b/Source/Microphysics/Null/NullMoist.H @@ -26,11 +26,6 @@ class NullMoist { void Advance (const amrex::Real& /*dt_advance*/) { } - virtual - void - Update (amrex::MultiFab& /*cons_in*/, - amrex::MultiFab& /*qmoist*/) { } - virtual void Update_Micro_Vars (amrex::MultiFab& /*cons_in*/) { } @@ -59,7 +54,12 @@ class NullMoist { int Qmoist_Size () { return NullMoist::m_qmoist_size; } + virtual + int + Qstate_Size () { return NullMoist::m_qstate_size; } + private: int m_qmoist_size = 1; + int m_qstate_size = 0; }; #endif diff --git a/Source/Microphysics/SAM/Diagnose_SAM.cpp b/Source/Microphysics/SAM/Diagnose_SAM.cpp index b671aabcc..60094eb02 100644 --- a/Source/Microphysics/SAM/Diagnose_SAM.cpp +++ b/Source/Microphysics/SAM/Diagnose_SAM.cpp @@ -34,7 +34,7 @@ void SAM::Diagnose () { amrex::ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - qv_array(i,j,k) = qt_array(i,j,k) - qn_array(i,j,k); + qt_array(i,j,k) = qv_array(i,j,k) + qn_array(i,j,k); amrex::Real omn = std::max(0.0, std::min(1.0,(tabs_array(i,j,k)-tbgmin)*a_bg)); qcl_array(i,j,k) = qn_array(i,j,k)*omn; qci_array(i,j,k) = qn_array(i,j,k)*(1.0-omn); @@ -49,33 +49,3 @@ void SAM::Diagnose () { } } - -/** - * Computes diagnostic quantities like temperature and pressure - * from the existing Microphysics variables. - */ -void SAM::Compute_Tabs_and_Pres () -{ - auto rho = mic_fab_vars[MicVar::rho]; - auto theta = mic_fab_vars[MicVar::theta]; - auto qv = mic_fab_vars[MicVar::qv]; - auto tabs = mic_fab_vars[MicVar::tabs]; - auto pres = mic_fab_vars[MicVar::pres]; - - for ( amrex::MFIter mfi(*tabs, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { - const auto& box3d = mfi.tilebox(); - - auto rho_array = rho->array(mfi); - auto theta_array = theta->array(mfi); - auto qv_array = qv->array(mfi); - auto tabs_array = tabs->array(mfi); - auto pres_array = pres->array(mfi); - - amrex::ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept - { - amrex::Real RhoTheta = rho_array(i,j,k)*theta_array(i,j,k); - tabs_array(i,j,k) = getTgivenRandRTh(rho_array(i,j,k), RhoTheta, qv_array(i,j,k)); - pres_array(i,j,k) = getPgivenRTh(RhoTheta, qv_array(i,j,k))/100.; - }); - } -} diff --git a/Source/Microphysics/SAM/Init_SAM.cpp b/Source/Microphysics/SAM/Init_SAM.cpp index 19970f4f3..7db0099a9 100644 --- a/Source/Microphysics/SAM/Init_SAM.cpp +++ b/Source/Microphysics/SAM/Init_SAM.cpp @@ -1,4 +1,3 @@ - #include #include "Microphysics.H" #include "IndexDefines.H" @@ -29,8 +28,8 @@ void SAM::Init (const MultiFab& cons_in, m_gtoe = grids; MicVarMap.resize(m_qmoist_size); - MicVarMap = {MicVar::qv, MicVar::qcl, MicVar::qci, - MicVar::qr, MicVar::qpi, MicVar::qg}; + MicVarMap = {MicVar::qt, MicVar::qv, MicVar::qcl, MicVar::qci, + MicVar::qp, MicVar::qpl, MicVar::qpi, MicVar::qg}; // initialize microphysics variables for (auto ivar = 0; ivar < MicVar::NumVars; ++ivar) { @@ -91,29 +90,36 @@ void SAM::Copy_State_to_Micro (const MultiFab& cons_in) auto states_array = cons_in.array(mfi); - auto qt_array = mic_fab_vars[MicVar::qt]->array(mfi); - auto qp_array = mic_fab_vars[MicVar::qp]->array(mfi); - auto qn_array = mic_fab_vars[MicVar::qn]->array(mfi); - auto qc_array = mic_fab_vars[MicVar::qcl]->array(mfi); - auto qi_array = mic_fab_vars[MicVar::qci]->array(mfi); - auto rho_array = mic_fab_vars[MicVar::rho]->array(mfi); - auto theta_array = mic_fab_vars[MicVar::theta]->array(mfi); + auto qv_array = mic_fab_vars[MicVar_Kess::qv]->array(mfi); + auto qc_array = mic_fab_vars[MicVar_Kess::qcl]->array(mfi); + auto qi_array = mic_fab_vars[MicVar_Kess::qci]->array(mfi); + auto qn_array = mic_fab_vars[MicVar_Kess::qn]->array(mfi); + auto qp_array = mic_fab_vars[MicVar_Kess::qp]->array(mfi); + + auto rho_array = mic_fab_vars[MicVar_Kess::rho]->array(mfi); + auto theta_array = mic_fab_vars[MicVar_Kess::theta]->array(mfi); + auto tabs_array = mic_fab_vars[MicVar::tabs]->array(mfi); + auto pres_array = mic_fab_vars[MicVar::pres]->array(mfi); // Get pressure, theta, temperature, density, and qt, qp amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { rho_array(i,j,k) = states_array(i,j,k,Rho_comp); theta_array(i,j,k) = states_array(i,j,k,RhoTheta_comp)/states_array(i,j,k,Rho_comp); - qt_array(i,j,k) = states_array(i,j,k,RhoQ1_comp)/states_array(i,j,k,Rho_comp); - qp_array(i,j,k) = states_array(i,j,k,RhoQ2_comp)/states_array(i,j,k,Rho_comp); + qv_array(i,j,k) = states_array(i,j,k,RhoQ1_comp)/states_array(i,j,k,Rho_comp); + qc_array(i,j,k) = states_array(i,j,k,RhoQ2_comp)/states_array(i,j,k,Rho_comp); + qp_array(i,j,k) = states_array(i,j,k,RhoQ3_comp)/states_array(i,j,k,Rho_comp); qn_array(i,j,k) = qc_array(i,j,k) + qi_array(i,j,k); - //qv_array(i,j,k) = qv_array_from_moist(i,j,k); - //temp_array(i,j,k) = getTgivenRandRTh(states_array(i,j,k,Rho_comp),states_array(i,j,k,RhoTheta_comp), qv_array(i,j,k)); - //pres_array(i,j,k) = getPgivenRTh(states_array(i,j,k,RhoTheta_comp))/100.; + + tabs_array(i,j,k) = getTgivenRandRTh(states_array(i,j,k,Rho_comp), + states_array(i,j,k,RhoTheta_comp), + qv_array(i,j,k)); + pres_array(i,j,k) = getPgivenRTh(states_array(i,j,k,RhoTheta_comp), qv_array(i,j,k))/100.; }); } } + void SAM::Compute_Coefficients () { auto dz = m_geom.CellSize(2); @@ -153,35 +159,39 @@ void SAM::Compute_Coefficients () // calculate the plane average variables PlaneAverage rho_ave(mic_fab_vars[MicVar::rho].get(), m_geom, m_axis); PlaneAverage theta_ave(mic_fab_vars[MicVar::theta].get(), m_geom, m_axis); + PlaneAverage qv_ave(mic_fab_vars[MicVar::qv].get(), m_geom, m_axis); rho_ave.compute_averages(ZDir(), rho_ave.field()); theta_ave.compute_averages(ZDir(), theta_ave.field()); + qv_ave.compute_averages(ZDir(), qv_ave.field()); // get host variable rho, and rhotheta int ncell = rho_ave.ncell_line(); - Gpu::HostVector rho_h(ncell), theta_h(ncell); + Gpu::HostVector rho_h(ncell), theta_h(ncell), qv_h(ncell); rho_ave.line_average(0, rho_h); theta_ave.line_average(0, theta_h); + qv_ave.line_average(0, qv_h); // copy data to device - Gpu::DeviceVector rho_d(ncell), theta_d(ncell); + Gpu::DeviceVector rho_d(ncell), theta_d(ncell), qv_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, qv_h.begin(), qv_h.end(), qv_d.begin()); Gpu::streamSynchronize(); Real* rho_dptr = rho_d.data(); Real* theta_dptr = theta_d.data(); + Real* qv_dptr = qv_d.data(); Real gOcp = m_gOcp; - // TODO: None of this is qv aware...? amrex::ParallelFor(nlev, [=] AMREX_GPU_DEVICE (int k) noexcept { Real RhoTheta = rho_dptr[k]*theta_dptr[k]; - Real pressure = getPgivenRTh(RhoTheta); + Real pressure = getPgivenRTh(RhoTheta, qv_dptr[k]); rho1d_t(k) = rho_dptr[k]; pres1d_t(k) = pressure/100.; - tabs1d_t(k) = getTgivenRandRTh(rho_dptr[k],RhoTheta); + tabs1d_t(k) = getTgivenRandRTh(rho_dptr[k], RhoTheta, qv_dptr[k]); zmid_t(k) = lowz + (k+0.5)*dz; gamaz_t(k) = gOcp*zmid_t(k); }); @@ -240,241 +250,3 @@ void SAM::Compute_Coefficients () (coef1+coef2) * pow(rho1d_t(k) , ((1.0+b_rain)/8.0))*sqrt(pratio); }); } - -#if 0 -/** - * Initializes the Microphysics module. - * - * @param[in] cons_in Conserved variables input - * @param[in] qc_in Cloud variables input - * @param[in,out] qv_in Vapor variables input - * @param[in] qi_in Ice variables input - * @param[in] grids The boxes on which we will evolve the solution - * @param[in] geom Geometry associated with these MultiFabs and grids - * @param[in] dt_advance Timestep for the advance - */ -void SAM::Init (const MultiFab& cons_in, MultiFab& qmoist, - const BoxArray& grids, - const Geometry& geom, - const Real& dt_advance) - { - m_geom = geom; - m_gtoe = grids; - - auto dz = m_geom.CellSize(2); - auto lowz = m_geom.ProbLo(2); - - dt = dt_advance; - - // initialize microphysics variables - for (auto ivar = 0; ivar < MicVar::NumVars; ++ivar) { - mic_fab_vars[ivar] = std::make_shared(cons_in.boxArray(), cons_in.DistributionMap(), 1, cons_in.nGrowVect()); - mic_fab_vars[ivar]->setVal(0.); - } - - // We must initialize to zero since we now need boundary values for the call to getP and we need all values filled - // The ghost cells of these arrays aren't filled in the boundary condition calls for the state - qmoist.setVal(0.); - - for ( MFIter mfi(cons_in, TileNoZ()); mfi.isValid(); ++mfi) { - - const auto& box3d = mfi.tilebox(); - - const auto& lo = amrex::lbound(box3d); - const auto& hi = amrex::ubound(box3d); - - nlev = box3d.length(2); - zlo = lo.z; - zhi = hi.z; - - // parameters - accrrc.resize({zlo}, {zhi}); - accrsi.resize({zlo}, {zhi}); - accrsc.resize({zlo}, {zhi}); - coefice.resize({zlo}, {zhi}); - evaps1.resize({zlo}, {zhi}); - evaps2.resize({zlo}, {zhi}); - accrgi.resize({zlo}, {zhi}); - accrgc.resize({zlo}, {zhi}); - evapg1.resize({zlo}, {zhi}); - evapg2.resize({zlo}, {zhi}); - evapr1.resize({zlo}, {zhi}); - evapr2.resize({zlo}, {zhi}); - - // data (input) - rho1d.resize({zlo}, {zhi}); - pres1d.resize({zlo}, {zhi}); - tabs1d.resize({zlo}, {zhi}); - gamaz.resize({zlo}, {zhi}); - zmid.resize({zlo}, {zhi}); - - // output - qifall.resize({zlo}, {zhi}); - tlatqi.resize({zlo}, {zhi}); - } - - auto accrrc_t = accrrc.table(); - auto accrsi_t = accrsi.table(); - auto accrsc_t = accrsc.table(); - auto coefice_t = coefice.table(); - auto evaps1_t = evaps1.table(); - auto evaps2_t = evaps2.table(); - auto accrgi_t = accrgi.table(); - auto accrgc_t = accrgc.table(); - auto evapg1_t = evapg1.table(); - auto evapg2_t = evapg2.table(); - auto evapr1_t = evapr1.table(); - auto evapr2_t = evapr2.table(); - - auto rho1d_t = rho1d.table(); - auto pres1d_t = pres1d.table(); - auto tabs1d_t = tabs1d.table(); - - auto gamaz_t = gamaz.table(); - auto zmid_t = zmid.table(); - - Real gam3 = erf_gammafff(3.0 ); - Real gamr1 = erf_gammafff(3.0+b_rain ); - Real gamr2 = erf_gammafff((5.0+b_rain)/2.0); - // Real gamr3 = erf_gammafff(4.0+b_rain ); - Real gams1 = erf_gammafff(3.0+b_snow ); - Real gams2 = erf_gammafff((5.0+b_snow)/2.0); - // Real gams3 = erf_gammafff(4.0+b_snow ); - Real gamg1 = erf_gammafff(3.0+b_grau ); - Real gamg2 = erf_gammafff((5.0+b_grau)/2.0); - // Real gamg3 = erf_gammafff(4.0+b_grau ); - - amrex::MultiFab qv(qmoist, amrex::make_alias, 0, 1); - amrex::MultiFab qc(qmoist, amrex::make_alias, 1, 1); - amrex::MultiFab qi(qmoist, amrex::make_alias, 2, 1); - - // Get the temperature, density, theta, qt and qp from input - for ( MFIter mfi(cons_in, false); mfi.isValid(); ++mfi) { - auto states_array = cons_in.array(mfi); - auto qc_array = qc.array(mfi); - auto qi_array = qi.array(mfi); - - auto qt_array = mic_fab_vars[MicVar::qt]->array(mfi); - auto qp_array = mic_fab_vars[MicVar::qp]->array(mfi); - auto qn_array = mic_fab_vars[MicVar::qn]->array(mfi); - auto rho_array = mic_fab_vars[MicVar::rho]->array(mfi); - auto theta_array = mic_fab_vars[MicVar::theta]->array(mfi); - auto temp_array = mic_fab_vars[MicVar::tabs]->array(mfi); - auto pres_array = mic_fab_vars[MicVar::pres]->array(mfi); - - const auto& box3d = mfi.tilebox(); - - // Get pressure, theta, temperature, density, and qt, qp - amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - rho_array(i,j,k) = states_array(i,j,k,Rho_comp); - theta_array(i,j,k) = states_array(i,j,k,RhoTheta_comp)/states_array(i,j,k,Rho_comp); - qt_array(i,j,k) = states_array(i,j,k,RhoQ1_comp)/states_array(i,j,k,Rho_comp); - qp_array(i,j,k) = states_array(i,j,k,RhoQ2_comp)/states_array(i,j,k,Rho_comp); - qn_array(i,j,k) = qc_array(i,j,k) + qi_array(i,j,k); - temp_array(i,j,k) = getTgivenRandRTh(states_array(i,j,k,Rho_comp),states_array(i,j,k,RhoTheta_comp)); - pres_array(i,j,k) = getPgivenRTh(states_array(i,j,k,RhoTheta_comp))/100.; - }); - } - - // calculate the plane average variables - PlaneAverage cons_ave(&cons_in, m_geom, m_axis); - cons_ave.compute_averages(ZDir(), cons_ave.field()); - - // get host variable rho, and rhotheta - int ncell = cons_ave.ncell_line(); - - Gpu::HostVector rho_h(ncell), rhotheta_h(ncell); - cons_ave.line_average(Rho_comp, rho_h); - cons_ave.line_average(RhoTheta_comp, rhotheta_h); - - // copy data to device - Gpu::DeviceVector rho_d(ncell), rhotheta_d(ncell); - Gpu::copyAsync(Gpu::hostToDevice, rho_h.begin(), rho_h.end(), rho_d.begin()); - Gpu::copyAsync(Gpu::hostToDevice, rhotheta_h.begin(), rhotheta_h.end(), rhotheta_d.begin()); - Gpu::streamSynchronize(); - - Real* rho_dptr = rho_d.data(); - Real* rhotheta_dptr = rhotheta_d.data(); - - Real gOcp = m_gOcp; - - amrex::ParallelFor(nlev, [=] AMREX_GPU_DEVICE (int k) noexcept { - Real pressure = getPgivenRTh(rhotheta_dptr[k]); - rho1d_t(k) = rho_dptr[k]; - pres1d_t(k) = pressure/100.; - tabs1d_t(k) = getTgivenRandRTh(rho_dptr[k],rhotheta_dptr[k]); - zmid_t(k) = lowz + (k+0.5)*dz; - gamaz_t(k) = gOcp*zmid_t(k); - }); - - // This fills qv - Diagnose(); - -#if 0 - amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int k, int j, int i) { - fluxbmk(l,j,i) = 0.0; - fluxtmk(l,j,i) = 0.0; - }); - - amrex::ParallelFor( box2d, [=] AMREX_GPU_DEVICE (int k, int l, int i0) { - mkwle (l,k) = 0.0; - mkwsb (l,k) = 0.0; - mkadv (l,k) = 0.0; - mkdiff(l,k) = 0.0; - }); -#endif - - if(round(gam3) != 2) { - std::cout << "cannot compute gamma-function in Microphysics::Init" << std::endl; - std::exit(-1); - } - - amrex::ParallelFor(nlev, [=] AMREX_GPU_DEVICE (int k) noexcept { - - Real pratio = sqrt(1.29 / rho1d_t(k)); - Real rrr1=393.0/(tabs1d_t(k)+120.0)*std::pow((tabs1d_t(k)/273.0),1.5); - Real rrr2=std::pow((tabs1d_t(k)/273.0),1.94)*(1000.0/pres1d_t(k)); - Real estw = 100.0*erf_esatw(tabs1d_t(k)); - Real esti = 100.0*erf_esati(tabs1d_t(k)); - - // accretion by snow: - Real coef1 = 0.25 * PI * nzeros * a_snow * gams1 * pratio/pow((PI * rhos * nzeros/rho1d_t(k) ) , ((3.0+b_snow)/4.0)); - Real coef2 = exp(0.025*(tabs1d_t(k) - 273.15)); - accrsi_t(k) = coef1 * coef2 * esicoef; - accrsc_t(k) = coef1 * esccoef; - coefice_t(k) = coef2; - - // evaporation of snow: - coef1 =(lsub/(tabs1d_t(k)*rv)-1.0)*lsub/(therco*rrr1*tabs1d_t(k)); - coef2 = rv*tabs1d_t(k)/(diffelq*rrr2*esti); - evaps1_t(k) = 0.65*4.0*nzeros/sqrt(PI*rhos*nzeros)/(coef1+coef2)/sqrt(rho1d_t(k)); - evaps2_t(k) = 0.49*4.0*nzeros*gams2*sqrt(a_snow/(muelq*rrr1))/pow((PI*rhos*nzeros) , ((5.0+b_snow)/8.0)) / - (coef1+coef2) * pow(rho1d_t(k) , ((1.0+b_snow)/8.0))*sqrt(pratio); - - // accretion by graupel: - coef1 = 0.25*PI*nzerog*a_grau*gamg1*pratio/pow((PI*rhog*nzerog/rho1d_t(k)) , ((3.0+b_grau)/4.0)); - coef2 = exp(0.025*(tabs1d_t(k) - 273.15)); - accrgi_t(k) = coef1 * coef2 * egicoef; - accrgc_t(k) = coef1 * egccoef; - - // evaporation of graupel: - coef1 =(lsub/(tabs1d_t(k)*rv)-1.0)*lsub/(therco*rrr1*tabs1d_t(k)); - coef2 = rv*tabs1d_t(k)/(diffelq*rrr2*esti); - evapg1_t(k) = 0.65*4.0*nzerog/sqrt(PI*rhog*nzerog)/(coef1+coef2)/sqrt(rho1d_t(k)); - evapg2_t(k) = 0.49*4.0*nzerog*gamg2*sqrt(a_grau/(muelq*rrr1))/pow((PI * rhog * nzerog) , ((5.0+b_grau)/8.0)) / - (coef1+coef2) * pow(rho1d_t(k) , ((1.0+b_grau)/8.0))*sqrt(pratio); - - // accretion by rain: - accrrc_t(k)= 0.25 * PI * nzeror * a_rain * gamr1 * pratio/pow((PI * rhor * nzeror / rho1d_t(k)) , ((3+b_rain)/4.))* erccoef; - - // evaporation of rain: - coef1 =(lcond/(tabs1d_t(k)*rv)-1.0)*lcond/(therco*rrr1*tabs1d_t(k)); - coef2 = rv*tabs1d_t(k)/(diffelq * rrr2 * estw); - evapr1_t(k) = 0.78 * 2.0 * PI * nzeror / - sqrt(PI * rhor * nzeror) / (coef1+coef2) / sqrt(rho1d_t(k)); - evapr2_t(k) = 0.31 * 2.0 * PI * nzeror * gamr2 * 0.89 * sqrt(a_rain/(muelq*rrr1))/ - pow((PI * rhor * nzeror) , ((5.0+b_rain)/8.0)) / - (coef1+coef2) * pow(rho1d_t(k) , ((1.0+b_rain)/8.0))*sqrt(pratio); - }); -} -#endif diff --git a/Source/Microphysics/SAM/SAM.H b/Source/Microphysics/SAM/SAM.H index 4cccf0995..41a416262 100644 --- a/Source/Microphysics/SAM/SAM.H +++ b/Source/Microphysics/SAM/SAM.H @@ -33,7 +33,6 @@ namespace MicVar { rho, // density pres, // pressure // derived variables - qr, // rain water qv, // water vapor qn, // cloud condensate (liquid+ice), initial to zero qci, // cloud ice @@ -107,17 +106,11 @@ public: void Copy_Micro_to_State (amrex::MultiFab& cons_in) override; - // update ERF variables - void - Update (amrex::MultiFab& cons_in, - amrex::MultiFab& qmoist) override; - void Update_Micro_Vars (amrex::MultiFab& cons_in) override { this->Copy_State_to_Micro(cons_in); this->Diagnose(); - this->Compute_Tabs_and_Pres(); this->Compute_Coefficients(); } @@ -148,20 +141,20 @@ public: } void - Compute_Tabs_and_Pres ( ); - - void - Compute_Coefficients ( ); + Compute_Coefficients (); int - Qmoist_Size ( ) override { return SAM::m_qmoist_size; } + Qmoist_Size () override { return SAM::m_qmoist_size; } - // independent variables - amrex::Array mic_fab_vars; + int + Qstate_Size () { return SAM::m_qstate_size; } private: + // Number of qmoist variables (qt, qv, qcl, qci, qp, qpl, qpi, qpg) + int m_qmoist_size = 8; + // Number of qmoist variables - int m_qmoist_size = 6; + int m_qstate_size = 3; // MicVar map (Qmoist indices -> MicVar enum) amrex::Vector MicVarMap; @@ -190,6 +183,9 @@ private: amrex::Real m_fac_sub; amrex::Real m_gOcp; + // independent variables + amrex::Array mic_fab_vars; + // microphysics parameters/coefficients amrex::TableData accrrc; amrex::TableData accrsi; diff --git a/Source/Microphysics/SAM/Update_SAM.cpp b/Source/Microphysics/SAM/Update_SAM.cpp index 08b7d33e6..0bb2d3f59 100644 --- a/Source/Microphysics/SAM/Update_SAM.cpp +++ b/Source/Microphysics/SAM/Update_SAM.cpp @@ -3,61 +3,6 @@ #include "IndexDefines.H" #include "TileNoZ.H" -/** - * Updates conserved and microphysics variables in the provided MultiFabs from - * the internal MultiFabs that store Microphysics module data. - * - * @param[out] cons Conserved variables - * @param[out] qmoist: qv, qc, qi, qr, qs, qg - */ -void SAM::Update (amrex::MultiFab& cons, - amrex::MultiFab& qmoist) -{ - // copy multifab data to qc, qv, and qi - amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar::qv], 0, 0, 1, mic_fab_vars[MicVar::qv]->nGrowVect()); // vapor - amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar::qcl], 0, 1, 1, mic_fab_vars[MicVar::qcl]->nGrowVect()); // cloud water - amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar::qci], 0, 2, 1, mic_fab_vars[MicVar::qci]->nGrowVect()); // cloud ice - amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar::qpl], 0, 3, 1, mic_fab_vars[MicVar::qpl]->nGrowVect()); // rain - amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar::qpi], 0, 4, 1, mic_fab_vars[MicVar::qpi]->nGrowVect()); // snow - - // Don't need to copy this since it is filled below - // amrex::MultiFab::Copy(qmoist, *mic_fab_vars[MicVar::qpi], 0, 5, 1, mic_fab_vars[MicVar::qci]->nGrowVect()); // graupel - - amrex::MultiFab qgraup_mf(qmoist, amrex::make_alias, 5, 1); - - // Get the temperature, density, theta, qt and qp from input - for ( amrex::MFIter mfi(cons,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { - auto states_arr = cons.array(mfi); - - auto rho_arr = mic_fab_vars[MicVar::rho]->array(mfi); - auto theta_arr = mic_fab_vars[MicVar::theta]->array(mfi); - auto qt_arr = mic_fab_vars[MicVar::qt]->array(mfi); - auto qp_arr = mic_fab_vars[MicVar::qp]->array(mfi); - auto qpl_arr = mic_fab_vars[MicVar::qpl]->array(mfi); - auto qpi_arr = mic_fab_vars[MicVar::qpi]->array(mfi); - - auto qgraup_arr= qgraup_mf.array(mfi); - - const auto& box3d = mfi.tilebox(); - - // get potential total density, temperature, qt, qp - amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - states_arr(i,j,k,Rho_comp) = rho_arr(i,j,k); - states_arr(i,j,k,RhoTheta_comp) = rho_arr(i,j,k)*theta_arr(i,j,k); - states_arr(i,j,k,RhoQ1_comp) = rho_arr(i,j,k)*qt_arr(i,j,k); - states_arr(i,j,k,RhoQ2_comp) = rho_arr(i,j,k)*qp_arr(i,j,k); - - // Graupel == precip total - rain - snow (but must be >= 0) - qgraup_arr(i,j,k) = std::max(0.0, qp_arr(i,j,k)-qpl_arr(i,j,k)-qpi_arr(i,j,k)); - }); - } - - // Fill interior ghost cells and periodic boundaries - cons.FillBoundary(m_geom.periodicity()); - qmoist.FillBoundary(m_geom.periodicity()); -} - /** * Updates conserved and microphysics variables in the provided MultiFabs from * the internal MultiFabs that store Microphysics module data. @@ -75,16 +20,17 @@ void SAM::Copy_Micro_to_State (amrex::MultiFab& cons) auto rho_arr = mic_fab_vars[MicVar_Kess::rho]->array(mfi); auto theta_arr = mic_fab_vars[MicVar_Kess::theta]->array(mfi); - auto qt_arr = mic_fab_vars[MicVar_Kess::qt]->array(mfi); + auto qv_arr = mic_fab_vars[MicVar_Kess::qv]->array(mfi); + auto qc_arr = mic_fab_vars[MicVar_Kess::qcl]->array(mfi); auto qp_arr = mic_fab_vars[MicVar_Kess::qp]->array(mfi); // get potential total density, temperature, qt, qp amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - //states_arr(i,j,k,Rho_comp) = rho_arr(i,j,k); states_arr(i,j,k,RhoTheta_comp) = rho_arr(i,j,k)*theta_arr(i,j,k); - states_arr(i,j,k,RhoQ1_comp) = rho_arr(i,j,k)*qt_arr(i,j,k); - states_arr(i,j,k,RhoQ2_comp) = rho_arr(i,j,k)*qp_arr(i,j,k); + states_arr(i,j,k,RhoQ1_comp) = rho_arr(i,j,k)*qv_arr(i,j,k); + states_arr(i,j,k,RhoQ2_comp) = rho_arr(i,j,k)*qc_arr(i,j,k); + states_arr(i,j,k,RhoQ3_comp) = rho_arr(i,j,k)*qp_arr(i,j,k); }); } diff --git a/Source/TimeIntegration/ERF_advance_dycore.cpp b/Source/TimeIntegration/ERF_advance_dycore.cpp index f6525cc5c..80f74eecb 100644 --- a/Source/TimeIntegration/ERF_advance_dycore.cpp +++ b/Source/TimeIntegration/ERF_advance_dycore.cpp @@ -251,47 +251,6 @@ void ERF::advance_dycore(int level, cons_to_prim(state_old[IntVar::cons], state_old[IntVar::cons].nGrow()); } // profile - 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); - - if (solverChoice.moisture_type != MoistureType::None) { - if (q_size >= 1) { - MultiFab qvapor (*(qmoist[level][0]), make_alias, 0, 1); - PlaneAverage qv_ave(&qvapor, geom[level], solverChoice.ave_plane); - qv_ave.compute_averages(ZDir(), qv_ave.field()); - - Gpu::HostVector qv_h(ncell); - - qv_ave.line_average(0, qv_h); - Gpu::copyAsync(Gpu::hostToDevice, qv_h.begin(), qv_h.end(), qv_d.begin()); - } - - if (q_size >= 2) { - MultiFab qcloud (*(qmoist[level][1]), make_alias, 0, 1); - PlaneAverage qc_ave(&qcloud, geom[level], solverChoice.ave_plane); - qc_ave.compute_averages(ZDir(), qc_ave.field()); - - 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][2]), make_alias, 0, 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" #include "TI_fast_rhs_fun.H" diff --git a/Source/TimeIntegration/ERF_make_buoyancy.cpp b/Source/TimeIntegration/ERF_make_buoyancy.cpp index c6dea21cb..ea1b67d03 100644 --- a/Source/TimeIntegration/ERF_make_buoyancy.cpp +++ b/Source/TimeIntegration/ERF_make_buoyancy.cpp @@ -36,10 +36,6 @@ using namespace amrex; void make_buoyancy (Vector& S_data, const MultiFab& S_prim, MultiFab& buoyancy, - Vector qmoist, - Gpu::DeviceVector qv_d, - Gpu::DeviceVector qc_d, - Gpu::DeviceVector qi_d, const amrex::Geometry geom, const SolverChoice& solverChoice, const MultiFab* r0) @@ -138,35 +134,10 @@ void make_buoyancy (Vector& S_data, // ****************************************************************************************** // Moist versions of buoyancy expressions // ****************************************************************************************** - if (solverChoice.moisture_type == MoistureType::FastEddy) { + if (solverChoice.moisture_type != MoistureType::None) { - AMREX_ALWAYS_ASSERT(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); - - 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); - - 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 ,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 if (solverChoice.moisture_type != MoistureType::None) { + if (solverChoice.moisture_type == MoistureType::FastEddy) + AMREX_ALWAYS_ASSERT(solverChoice.buoyancy_type == 1); if (solverChoice.buoyancy_type == 1) { @@ -184,16 +155,13 @@ void make_buoyancy (Vector& S_data, // Base state density const Array4& r0_arr = r0->const_array(mfi); - const Array4 & qv_data = qmoist[0]->array(mfi); - const Array4 & qc_data = qmoist[1]->array(mfi); - const Array4 & qi_data = qmoist[2]->array(mfi); + // TODO: A microphysics model may have more than q1 & q2 components for the + // non-precipitating phase. 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 @@ -209,7 +177,7 @@ void make_buoyancy (Vector& S_data, int ncell = state_ave.ncell_line(); - Gpu::HostVector rho_h(ncell), theta_h(ncell), qp_h(ncell); + Gpu::HostVector rho_h(ncell), theta_h(ncell); Gpu::DeviceVector rho_d(ncell), theta_d(ncell); state_ave.line_average(Rho_comp, rho_h); @@ -221,21 +189,28 @@ void make_buoyancy (Vector& S_data, Real* rho_d_ptr = rho_d.data(); Real* theta_d_ptr = theta_d.data(); - if (solverChoice.buoyancy_type == 2 || solverChoice.buoyancy_type == 4 ) { - - prim_ave.line_average(PrimQ2_comp, qp_h); - - Gpu::DeviceVector qp_d(ncell); - - // Copy data to device - Gpu::copyAsync(Gpu::hostToDevice, qp_h.begin(), qp_h.end(), qp_d.begin()); - Gpu::streamSynchronize(); - - 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(); + // Average valid moisture vars + int n_prim_max = NVAR_max - 1; + int n_moist_var = NMOIST_max - (S_prim.nComp() - n_prim_max); + Gpu::HostVector qv_h(ncell) , qc_h(ncell) , qp_h(ncell); + Gpu::DeviceVector qv_d(ncell,0.0), qc_d(ncell,0.0), qp_d(ncell,0.0); + if (n_moist_var >=1) { + prim_ave.line_average(PrimQ1_comp, qv_h); + Gpu::copyAsync(Gpu::hostToDevice, qv_h.begin(), qv_h.end(), qv_d.begin()); + } + if (n_moist_var >=2) { + prim_ave.line_average(PrimQ2_comp, qc_h); + Gpu::copyAsync(Gpu::hostToDevice, qc_h.begin(), qc_h.end(), qc_d.begin()); + } + if (n_moist_var >=3) { + prim_ave.line_average(PrimQ3_comp, qp_h); + Gpu::copyAsync(Gpu::hostToDevice, qp_h.begin(), qp_h.end(), qp_d.begin()); + } + Real* qv_d_ptr = qv_d.data(); + Real* qc_d_ptr = qc_d.data(); + Real* qp_d_ptr = qp_d.data(); + if (solverChoice.buoyancy_type == 2 || solverChoice.buoyancy_type == 4 ) { #ifdef _OPENMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) #endif @@ -252,10 +227,7 @@ void make_buoyancy (Vector& S_data, const Array4 & cell_data = S_data[IntVar::cons].array(mfi); const Array4 & cell_prim = S_prim.array(mfi); - const Array4 & qv_data = qmoist[0]->const_array(mfi); - const Array4 & qc_data = qmoist[1]->const_array(mfi); - const Array4 & qi_data = qmoist[2]->const_array(mfi); - + // TODO: ice has not been dealt with (q1=qv, q2=qv, q3=qp) 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 ]); @@ -266,32 +238,36 @@ void make_buoyancy (Vector& S_data, Real qplus, qminus; + Real qv_plus = (n_moist_var >= 1) ? cell_prim(i,j,k ,PrimQ1_comp) : 0.0; + Real qv_minus = (n_moist_var >= 1) ? cell_prim(i,j,k-1,PrimQ1_comp) : 0.0; + + Real qc_plus = (n_moist_var >= 2) ? cell_prim(i,j,k ,PrimQ2_comp) : 0.0; + Real qc_minus = (n_moist_var >= 2) ? cell_prim(i,j,k-1,PrimQ2_comp) : 0.0; + + Real qp_plus = (n_moist_var >= 3) ? cell_prim(i,j,k ,PrimQ3_comp) : 0.0; + Real qp_minus = (n_moist_var >= 3) ? cell_prim(i,j,k-1,PrimQ3_comp) : 0.0; + 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]); - - 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]); + qplus = 0.61 * ( qv_plus - qv_d_ptr[k] ) - + ( qc_plus - qc_d_ptr[k] + + qp_plus - qp_d_ptr[k] ) + + (tempp3d-tempp1d)/tempp1d*(Real(1.0) + Real(0.61)*qv_d_ptr[k]-qc_d_ptr[k]-qp_d_ptr[k]); - } else if (solverChoice.buoyancy_type == 4) { + qminus = 0.61 * ( qv_minus - qv_d_ptr[k-1] ) - + ( qc_minus - qc_d_ptr[k-1] + + qp_minus - qp_d_ptr[k-1] ) + + (tempm3d-tempm1d)/tempm1d*(Real(1.0) + Real(0.61)*qv_d_ptr[k-1]-qc_d_ptr[k-1]-qp_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]) + } else if (solverChoice.buoyancy_type == 4) { + qplus = 0.61 * ( qv_plus - qv_d_ptr[k] ) - + ( qc_plus - qc_d_ptr[k] + + qp_plus - qp_d_ptr[k] ) + (cell_data(i,j,k ,RhoTheta_comp)/cell_data(i,j,k ,Rho_comp) - theta_d_ptr[k ])/theta_d_ptr[k ]; - qminus = 0.61 *( qv_data(i,j,k-1)-qv_d_ptr[k-1]) - - (qc_data(i,j,k-1)-qc_d_ptr[k-1]+ - qi_data(i,j,k-1)-qi_d_ptr[k-1]+ - cell_prim(i,j,k-1,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]; + qminus = 0.61 * ( qv_minus - qv_d_ptr[k-1] ) - + ( qc_minus - qc_d_ptr[k-1] + + qp_minus - qp_d_ptr[k-1] ) + + (cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp) - theta_d_ptr[k-1])/theta_d_ptr[k-1]; } Real qavg = Real(0.5) * (qplus + qminus); @@ -302,7 +278,6 @@ void make_buoyancy (Vector& S_data, } // mfi } else if (solverChoice.buoyancy_type == 3) { - #ifdef _OPENMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) #endif @@ -318,9 +293,8 @@ void make_buoyancy (Vector& S_data, const Array4 & cell_data = S_data[IntVar::cons].array(mfi); const Array4 & cell_prim = S_prim.array(mfi); - const Array4 & qv_data = qmoist[0]->const_array(mfi); - const Array4 & qc_data = qmoist[1]->const_array(mfi); - const Array4 & qi_data = qmoist[2]->const_array(mfi); + + // TODO: ice has not been dealt with (q1=qv, q2=qv, q3=qp) amrex::ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { @@ -330,11 +304,18 @@ void make_buoyancy (Vector& S_data, Real tempp3d = getTgivenRandRTh(cell_data(i,j,k ,Rho_comp), cell_data(i,j,k ,RhoTheta_comp)); Real tempm3d = getTgivenRandRTh(cell_data(i,j,k-1,Rho_comp), cell_data(i,j,k-1,RhoTheta_comp)); - Real qplus = 0.61 * qv_data(i,j,k) - (qc_data(i,j,k)+ qi_data(i,j,k)+ cell_prim(i,j,k,PrimQ2_comp)) - + (tempp3d-tempp1d)/tempp1d; + Real qv_plus = (n_moist_var >= 1) ? cell_prim(i,j,k ,PrimQ1_comp) : 0.0; + Real qv_minus = (n_moist_var >= 1) ? cell_prim(i,j,k-1,PrimQ1_comp) : 0.0; + + Real qc_plus = (n_moist_var >= 2) ? cell_prim(i,j,k ,PrimQ2_comp) : 0.0; + Real qc_minus = (n_moist_var >= 2) ? cell_prim(i,j,k-1,PrimQ2_comp) : 0.0; + + Real qp_plus = (n_moist_var >= 3) ? cell_prim(i,j,k ,PrimQ3_comp) : 0.0; + Real qp_minus = (n_moist_var >= 3) ? cell_prim(i,j,k-1,PrimQ3_comp) : 0.0; + + Real qplus = 0.61 * qv_plus - (qc_plus + qp_plus) + (tempp3d-tempp1d)/tempp1d; - 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 qminus = 0.61 * qv_minus - (qc_minus + qp_minus) + (tempm3d-tempm1d)/tempm1d; Real qavg = Real(0.5) * (qplus + qminus); Real r0avg = Real(0.5) * (rho_d_ptr[k] + rho_d_ptr[k-1]); diff --git a/Source/TimeIntegration/ERF_slow_rhs_post.cpp b/Source/TimeIntegration/ERF_slow_rhs_post.cpp index bba64d6f9..3a8fafb6e 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_post.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_post.cpp @@ -170,8 +170,6 @@ void erf_slow_rhs_post (int level, int finest_level, { std::array flux; - int ncomp = S_data[IntVar::cons].nComp(); - int start_comp; int num_comp; @@ -183,7 +181,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),ncomp); + flux[dir].resize(amrex::surroundingNodes(tbx,dir),nvars); flux[dir].setVal(0.); } const GpuArray, AMREX_SPACEDIM> @@ -300,7 +298,7 @@ void erf_slow_rhs_post (int level, int finest_level, if (solverChoice.moisture_type != MoistureType::None) { start_comp = RhoQ1_comp; - num_comp = 2; + num_comp = nvars - start_comp; AdvType moist_horiz_adv_type = ac.moistscal_horiz_adv_type; AdvType moist_vert_adv_type = ac.moistscal_vert_adv_type; @@ -370,8 +368,9 @@ void erf_slow_rhs_post (int level, int finest_level, new_cons, cell_rhs, mf_u, mf_v, false, false); } } + start_comp = RhoScalar_comp; - num_comp = S_data[IntVar::cons].nComp() - start_comp; + num_comp = nvars - start_comp; if (l_use_terrain) { DiffusionSrcForState_T(tbx, domain, start_comp, num_comp, u, v, cur_cons, cur_prim, cell_rhs, @@ -460,7 +459,7 @@ void erf_slow_rhs_post (int level, int finest_level, } else { auto const& src_arr = source.const_array(mfi); start_comp = RhoScalar_comp; - num_comp = S_data[IntVar::cons].nComp() - start_comp; + num_comp = nvars - start_comp; ParallelFor(tbx, num_comp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int nn) noexcept { const int n = start_comp + nn; @@ -535,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 = ncomp - strt_comp_reflux; + int num_comp_reflux = nvars - 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 0a0e042eb..2eaa94cb0 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp @@ -80,7 +80,6 @@ void erf_slow_rhs_pre (int level, int finest_level, const MultiFab& xvel, const MultiFab& yvel, const MultiFab& zvel, - Vector qmoist, std::unique_ptr& z_t_mf, MultiFab& Omega, const MultiFab& source, @@ -543,7 +542,6 @@ void erf_slow_rhs_pre (int level, int finest_level, FArrayBox pprime; pprime.resize(gbx,1); Elixir pp_eli = pprime.elixir(); const Array4 & pp_arr = pprime.array(); - const Array4 & qv_arr = (use_moisture) ? qmoist[0]->const_array(mfi) : Array4 {}; { BL_PROFILE("slow_rhs_pre_pprime"); ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept @@ -551,10 +549,7 @@ void erf_slow_rhs_pre (int level, int finest_level, //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.); - Real qv_for_p = 0.0; - if (use_moisture) { - qv_for_p = qv_arr(i,j,k); - } + Real qv_for_p = (use_moisture) ? cell_data(i,j,k,RhoQ1_comp)/cell_data(i,j,k,Rho_comp) : 0.0; pp_arr(i,j,k) = getPgivenRTh(cell_data(i,j,k,RhoTheta_comp),qv_for_p) - p0_arr(i,j,k); }); } // end profile diff --git a/Source/TimeIntegration/TI_headers.H b/Source/TimeIntegration/TI_headers.H index 6e72b11d4..d7c676f7f 100644 --- a/Source/TimeIntegration/TI_headers.H +++ b/Source/TimeIntegration/TI_headers.H @@ -22,7 +22,6 @@ void erf_slow_rhs_pre (int level, int finest_level, int nrk, const amrex::MultiFab& xvel, const amrex::MultiFab& yvel, const amrex::MultiFab& zvel, - amrex::Vector qmoist, std::unique_ptr& z_t, amrex::MultiFab& Omega, const amrex::MultiFab& source, @@ -218,10 +217,6 @@ void make_fast_coeffs (int level, void make_buoyancy (amrex::Vector< amrex::MultiFab>& S_data, const amrex::MultiFab & S_prim, amrex::MultiFab& buoyancy, - amrex::Vector qmoist, - const amrex::Gpu::DeviceVector qv_d, - const amrex::Gpu::DeviceVector qc_d, - const amrex::Gpu::DeviceVector qi_d, const amrex::Geometry geom, const SolverChoice& solverChoice, const amrex::MultiFab* r0); diff --git a/Source/TimeIntegration/TI_slow_rhs_fun.H b/Source/TimeIntegration/TI_slow_rhs_fun.H index f783597e5..9b20fa5dd 100644 --- a/Source/TimeIntegration/TI_slow_rhs_fun.H +++ b/Source/TimeIntegration/TI_slow_rhs_fun.H @@ -90,12 +90,10 @@ MultiFab* p0_new = &p_hse_new; make_buoyancy(S_data, S_prim, buoyancy, - qmoist[level], 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, - qmoist[level], z_t_rk[level], Omega, source, buoyancy, Tau11, Tau22, Tau33, Tau12, Tau13, Tau21, Tau23, Tau31, Tau32, SmnSmn, eddyDiffs, Hfx3, Diss, @@ -183,12 +181,10 @@ // If not moving_terrain make_buoyancy(S_data, S_prim, buoyancy, - qmoist[level], 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, - qmoist[level], z_t_rk[level], Omega, source, buoyancy, Tau11, Tau22, Tau33, Tau12, Tau13, Tau21, Tau23, Tau31, Tau32, SmnSmn, eddyDiffs, Hfx3, Diss, @@ -353,13 +349,11 @@ // If not moving_terrain make_buoyancy(S_data, S_prim, buoyancy, - 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, - 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/prob_common.H b/Source/prob_common.H index a90aebee6..c2f78f5df 100644 --- a/Source/prob_common.H +++ b/Source/prob_common.H @@ -93,9 +93,6 @@ public: amrex::Array4 const& /*p_hse*/, amrex::Array4 const& /*z_nd*/, amrex::Array4 const& /*z_cc*/, - amrex::Array4 const& /*qv*/, - amrex::Array4 const& /*qc*/, - amrex::Array4 const& /*qi*/, amrex::GeometryData const& /*geomdata*/, amrex::Array4 const& /*mf_m*/, amrex::Array4 const& /*mf_u*/,