diff --git a/Source/DataStructs/ERF_DataStruct.H b/Source/DataStructs/ERF_DataStruct.H index dbb90e15c..b049e2894 100644 --- a/Source/DataStructs/ERF_DataStruct.H +++ b/Source/DataStructs/ERF_DataStruct.H @@ -222,6 +222,7 @@ struct SolverChoice { } pp.query("moisture_tight_coupling",moisture_tight_coupling); + pp.query("macrophysics_acoustic_coupling",macrophysics_acoustic_coupling); } // Which expression (1,2/3 or 4) to use for buoyancy @@ -1265,7 +1266,8 @@ struct SolverChoice { // Microphysics params MoistureComponentIndices moisture_indices; - bool moisture_tight_coupling {false}; + bool moisture_tight_coupling {false}; + bool macrophysics_acoustic_coupling {false}; std::string windfarm_loc_table, windfarm_spec_table, windfarm_spec_table_extra; std::string windfarm_blade_table, windfarm_airfoil_tables; diff --git a/Source/Diffusion/ERF_Diffusion.H b/Source/Diffusion/ERF_Diffusion.H index 631f8e07a..1ae5a3d28 100644 --- a/Source/Diffusion/ERF_Diffusion.H +++ b/Source/Diffusion/ERF_Diffusion.H @@ -456,9 +456,10 @@ void ComputeStrain_T (amrex::Box bxcc, amrex::Box tbxxy, void ImplicitDiffForStateLU_N (const amrex::Box& bx, const amrex::Box& domain, - const int level, - const int n, - const amrex::Real dt, + const int& level, + const int& n, + const int& qty_index, + const amrex::Real& dt, const amrex::GpuArray& bc_neumann_vals, const amrex::Array4< amrex::Real>& cell_data, const amrex::GpuArray& cellSizeInv, @@ -466,14 +467,15 @@ void ImplicitDiffForStateLU_N (const amrex::Box& bx, const amrex::Array4& mu_turb, const SolverChoice& solverChoice, const amrex::BCRec* bc_ptr, - const bool use_SurfLayer, - const amrex::Real implicit_fac); + const bool& use_SurfLayer, + const amrex::Real& implicit_fac); void ImplicitDiffForStateLU_S (const amrex::Box& bx, const amrex::Box& domain, - const int level, - const int n, - const amrex::Real dt, + const int& level, + const int& n, + const int& qty_index, + const amrex::Real& dt, const amrex::GpuArray& bc_neumann_vals, const amrex::Array4< amrex::Real>& cell_data, const amrex::Gpu::DeviceVector& stretched_dz_d, @@ -481,14 +483,15 @@ void ImplicitDiffForStateLU_S (const amrex::Box& bx, const amrex::Array4& mu_turb, const SolverChoice& solverChoice, const amrex::BCRec* bc_ptr, - const bool use_SurfLayer, - const amrex::Real implicit_fac); + const bool& use_SurfLayer, + const amrex::Real& implicit_fac); void ImplicitDiffForStateLU_T (const amrex::Box& bx, const amrex::Box& domain, - const int level, - const int n, - const amrex::Real dt, + const int& level, + const int& n, + const int& qty_index, + const amrex::Real& dt, const amrex::GpuArray& bc_neumann_vals, const amrex::Array4< amrex::Real>& cell_data, const amrex::Array4& z_nd, @@ -498,14 +501,14 @@ void ImplicitDiffForStateLU_T (const amrex::Box& bx, const amrex::Array4& mu_turb, const SolverChoice& solverChoice, const amrex::BCRec* bc_ptr, - const bool use_SurfLayer, - const amrex::Real implicit_fac); + const bool& use_SurfLayer, + const amrex::Real& implicit_fac); template void ImplicitDiffForMomLU_N (const amrex::Box& bx, const amrex::Box& domain, - const int level, - const amrex::Real dt, + const int& level, + const amrex::Real& dt, const amrex::Array4& cell_data, const amrex::Array4< amrex::Real>& face_data, const amrex::Array4& tau, @@ -514,14 +517,14 @@ void ImplicitDiffForMomLU_N (const amrex::Box& bx, const amrex::Array4& mu_turb, const SolverChoice& solverChoice, const amrex::BCRec* bc_ptr, - const bool use_SurfLayer, - const amrex::Real implicit_fac); + const bool& use_SurfLayer, + const amrex::Real& implicit_fac); template void ImplicitDiffForMomLU_S (const amrex::Box& bx, const amrex::Box& domain, - const int level, - const amrex::Real dt, + const int& level, + const amrex::Real& dt, const amrex::Array4& cell_data, const amrex::Array4< amrex::Real>& face_data, const amrex::Array4& tau, @@ -530,14 +533,14 @@ void ImplicitDiffForMomLU_S (const amrex::Box& bx, const amrex::Array4& mu_turb, const SolverChoice& solverChoice, const amrex::BCRec* bc_ptr, - const bool use_SurfLayer, - const amrex::Real implicit_fac); + const bool& use_SurfLayer, + const amrex::Real& implicit_fac); template void ImplicitDiffForMomLU_T (const amrex::Box& bx, const amrex::Box& domain, - const int level, - const amrex::Real dt, + const int& level, + const amrex::Real& dt, const amrex::Array4& cell_data, const amrex::Array4< amrex::Real>& face_data, const amrex::Array4& tau, @@ -548,6 +551,6 @@ void ImplicitDiffForMomLU_T (const amrex::Box& bx, const amrex::Array4& mu_turb, const SolverChoice& solverChoice, const amrex::BCRec* bc_ptr, - const bool use_SurfLayer, - const amrex::Real implicit_fac); + const bool& use_SurfLayer, + const amrex::Real& implicit_fac); #endif diff --git a/Source/Diffusion/ERF_ImplicitDiff_N.cpp b/Source/Diffusion/ERF_ImplicitDiff_N.cpp index af3421ab6..aea3f5321 100644 --- a/Source/Diffusion/ERF_ImplicitDiff_N.cpp +++ b/Source/Diffusion/ERF_ImplicitDiff_N.cpp @@ -28,9 +28,10 @@ using namespace amrex; void ImplicitDiffForStateLU_N (const Box& bx, const Box& domain, - const int level, - const int n, - const Real dt, + const int& level, + const int& n, + const int& qty_index, + const Real& dt, const GpuArray& bc_neumann_vals, const Array4< Real>& cell_data, const GpuArray& cellSizeInv, @@ -38,14 +39,13 @@ ImplicitDiffForStateLU_N (const Box& bx, const Array4& mu_turb, const SolverChoice &solverChoice, const BCRec* bc_ptr, - const bool use_SurfLayer, - const Real implicit_fac) + const bool& use_SurfLayer, + const Real& implicit_fac) { BL_PROFILE_VAR("ImplicitDiffForState_N()",ImplicitDiffForState_N); // setup quantities for getRhoAlpha() #include "ERF_SetupVertDiff.H" - const int qty_index = n; const int prim_index = qty_index - 1; const int prim_scal_index = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? PrimScalar_comp : prim_index; @@ -203,8 +203,8 @@ template void ImplicitDiffForMomLU_N (const Box& bx, const Box& /*domain*/, - const int level, - const Real dt, + const int& level, + const Real& dt, const Array4& cell_data, const Array4< Real>& face_data, const Array4& tau, @@ -213,8 +213,8 @@ ImplicitDiffForMomLU_N (const Box& bx, const Array4& mu_turb, const SolverChoice &solverChoice, const BCRec* bc_ptr, - const bool use_SurfLayer, - const Real implicit_fac) + const bool& use_SurfLayer, + const Real& implicit_fac) { BL_PROFILE_VAR("ImplicitDiffForMom_N()",ImplicitDiffForMom_N); @@ -430,8 +430,8 @@ ImplicitDiffForMomLU_N (const Box& bx, template void ImplicitDiffForMomLU_N ( \ const Box&, \ const Box&, \ - const int, \ - const Real, \ + const int&, \ + const Real&, \ const Array4&, \ const Array4< Real>&, \ const Array4&, \ @@ -440,8 +440,8 @@ ImplicitDiffForMomLU_N (const Box& bx, const Array4&, \ const SolverChoice&, \ const BCRec*, \ - const bool, \ - const Real); + const bool&, \ + const Real&); INSTANTIATE_IMPLICIT_DIFF_FOR_MOM_LU(0) INSTANTIATE_IMPLICIT_DIFF_FOR_MOM_LU(1) INSTANTIATE_IMPLICIT_DIFF_FOR_MOM_LU(2) diff --git a/Source/Diffusion/ERF_ImplicitDiff_S.cpp b/Source/Diffusion/ERF_ImplicitDiff_S.cpp index 6fbad2f2c..02cc35f4e 100644 --- a/Source/Diffusion/ERF_ImplicitDiff_S.cpp +++ b/Source/Diffusion/ERF_ImplicitDiff_S.cpp @@ -26,9 +26,10 @@ using namespace amrex; void ImplicitDiffForStateLU_S (const Box& bx, const Box& domain, - const int level, - const int n, - const Real dt, + const int& level, + const int& n, + const int& qty_index, + const Real& dt, const GpuArray& bc_neumann_vals, const Array4< Real>& cell_data, const Gpu::DeviceVector& stretched_dz_d, @@ -36,14 +37,13 @@ ImplicitDiffForStateLU_S (const Box& bx, const Array4& mu_turb, const SolverChoice &solverChoice, const BCRec* bc_ptr, - const bool use_SurfLayer, - const Real implicit_fac) + const bool& use_SurfLayer, + const Real& implicit_fac) { BL_PROFILE_VAR("ImplicitDiffForState_S()",ImplicitDiffForState_S); // setup quantities for getRhoAlpha() #include "ERF_SetupVertDiff.H" - const int qty_index = n; const int prim_index = qty_index - 1; const int prim_scal_index = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? PrimScalar_comp : prim_index; @@ -213,8 +213,8 @@ template void ImplicitDiffForMomLU_S (const Box& bx, const Box& /*domain*/, - const int level, - const Real dt, + const int& level, + const Real& dt, const Array4& cell_data, const Array4< Real>& face_data, const Array4& tau, @@ -223,8 +223,8 @@ ImplicitDiffForMomLU_S (const Box& bx, const Array4& mu_turb, const SolverChoice &solverChoice, const BCRec* bc_ptr, - const bool use_SurfLayer, - const Real implicit_fac) + const bool& use_SurfLayer, + const Real& implicit_fac) { BL_PROFILE_VAR("ImplicitDiffForMom_S()",ImplicitDiffForMom_S); @@ -453,8 +453,8 @@ ImplicitDiffForMomLU_S (const Box& bx, template void ImplicitDiffForMomLU_S ( \ const Box&, \ const Box&, \ - const int, \ - const Real, \ + const int&, \ + const Real&, \ const Array4&, \ const Array4< Real>&, \ const Array4&, \ @@ -463,8 +463,8 @@ ImplicitDiffForMomLU_S (const Box& bx, const Array4&, \ const SolverChoice&, \ const BCRec*, \ - const bool, \ - const Real); + const bool&, \ + const Real&); INSTANTIATE_IMPLICIT_DIFF_FOR_MOM_LU(0) INSTANTIATE_IMPLICIT_DIFF_FOR_MOM_LU(1) INSTANTIATE_IMPLICIT_DIFF_FOR_MOM_LU(2) diff --git a/Source/Diffusion/ERF_ImplicitDiff_T.cpp b/Source/Diffusion/ERF_ImplicitDiff_T.cpp index 16c937ad3..6c94b326c 100644 --- a/Source/Diffusion/ERF_ImplicitDiff_T.cpp +++ b/Source/Diffusion/ERF_ImplicitDiff_T.cpp @@ -28,9 +28,10 @@ using namespace amrex; void ImplicitDiffForStateLU_T (const Box& bx, const Box& domain, - const int level, - const int n, - const Real dt, + const int& level, + const int& n, + const int& qty_index, + const Real& dt, const GpuArray& bc_neumann_vals, const Array4< Real>& cell_data, const Array4& z_nd, @@ -40,14 +41,13 @@ ImplicitDiffForStateLU_T (const Box& bx, const Array4& mu_turb, const SolverChoice &solverChoice, const BCRec* bc_ptr, - const bool use_SurfLayer, - const Real implicit_fac) + const bool& use_SurfLayer, + const Real& implicit_fac) { BL_PROFILE_VAR("ImplicitDiffForState_T()",ImplicitDiffForState_T); // setup quantities for getRhoAlpha() #include "ERF_SetupVertDiff.H" - const int qty_index = n; const int prim_index = qty_index - 1; const int prim_scal_index = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? PrimScalar_comp : prim_index; @@ -221,8 +221,8 @@ template void ImplicitDiffForMomLU_T (const Box& bx, const Box& /*domain*/, - const int level, - const Real dt, + const int& level, + const Real& dt, const Array4& cell_data, const Array4< Real>& face_data, const Array4& tau, @@ -233,8 +233,8 @@ ImplicitDiffForMomLU_T (const Box& bx, const Array4& mu_turb, const SolverChoice &solverChoice, const BCRec* bc_ptr, - const bool use_SurfLayer, - const Real implicit_fac) + const bool& use_SurfLayer, + const Real& implicit_fac) { BL_PROFILE_VAR("ImplicitDiffForMom_T()",ImplicitDiffForMom_T); @@ -473,8 +473,8 @@ ImplicitDiffForMomLU_T (const Box& bx, template void ImplicitDiffForMomLU_T ( \ const Box&, \ const Box&, \ - const int, \ - const Real, \ + const int&, \ + const Real&, \ const Array4&, \ const Array4< Real>&, \ const Array4&, \ @@ -485,8 +485,8 @@ ImplicitDiffForMomLU_T (const Box& bx, const Array4&, \ const SolverChoice&, \ const BCRec*, \ - const bool, \ - const Real); + const bool&, \ + const Real&); INSTANTIATE_IMPLICIT_DIFF_FOR_MOM_LU(0) INSTANTIATE_IMPLICIT_DIFF_FOR_MOM_LU(1) INSTANTIATE_IMPLICIT_DIFF_FOR_MOM_LU(2) diff --git a/Source/ERF.H b/Source/ERF.H index 61289e60a..9c4cd1e06 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -24,7 +24,7 @@ #include #include -#include +#include "ERF_EB.H" #ifdef ERF_USE_FFT #include @@ -36,20 +36,20 @@ #include "ERF_ProbCommon.H" -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include +#include "ERF_IndexDefines.H" +#include "ERF_DataStruct.H" +#include "ERF_TurbPertStruct.H" +#include "ERF_InputSoundingData.H" +#include "ERF_InputSpongeData.H" +#include "ERF_SurfaceLayer.H" +#include "ERF_Derive.H" +#include "ERF_ReadBndryPlanes.H" +#include "ERF_WriteBndryPlanes.H" +#include "ERF_MRI.H" +#include "ERF_PhysBCFunct.H" +#include "ERF_FillPatcher.H" +#include "ERF_SampleData.H" +#include "ERF_ForestDrag.H" #ifdef ERF_USE_PARTICLES #include "ERF_ParticleData.H" @@ -63,7 +63,7 @@ #include "ERF_WindFarm.H" #endif -#include +#include "ERF_RadiationModels.H" #ifdef ERF_USE_SHOC #include "ERF_ShocInterface.H" @@ -880,6 +880,7 @@ private: std::unique_ptr micro; amrex::Vector> qmoist; // (lev,ncomp) rain_accum, snow_accum, graup_accum + amrex::Vector> macrophysics_src; // RhoTheta source for macro coupling // Variables for wind farm parametrization models diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 1903db87f..b76c7b978 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -2156,6 +2156,7 @@ ERF::initializeMicrophysics (const int& a_nlevsmax /*!< number of AMR levels */) } qmoist.resize(a_nlevsmax); + macrophysics_src.resize(a_nlevsmax); return; } diff --git a/Source/ERF_MakeNewArrays.cpp b/Source/ERF_MakeNewArrays.cpp index 7242e12f8..0a86cf6f3 100644 --- a/Source/ERF_MakeNewArrays.cpp +++ b/Source/ERF_MakeNewArrays.cpp @@ -43,7 +43,6 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm, // ******************************************************************************************** // Allocate terrain arrays // ******************************************************************************************** - BoxArray ba_nd(ba); ba_nd.surroundingNodes(); @@ -465,6 +464,16 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm, rad_fluxes[lev]->setVal(zero); } + //********************************************************* + // Macrophysics heating source terms + //********************************************************* + if (solverChoice.moisture_type != MoistureType::None && + solverChoice.macrophysics_acoustic_coupling) + { + macrophysics_src[lev] = std::make_unique(ba, dm, 3, 0); + macrophysics_src[lev]->setVal(0); + } + //********************************************************* // Turbulent perturbation region initialization //********************************************************* diff --git a/Source/Microphysics/Kessler/ERF_Kessler.H b/Source/Microphysics/Kessler/ERF_Kessler.H index c6ace4afc..264927bb3 100644 --- a/Source/Microphysics/Kessler/ERF_Kessler.H +++ b/Source/Microphysics/Kessler/ERF_Kessler.H @@ -56,7 +56,7 @@ public: { m_fac_cond = L_v / sc.c_p; m_axis = sc.ave_plane; - m_do_cond = (!sc.use_shoc); + m_do_cond = (!sc.use_shoc && !sc.macrophysics_acoustic_coupling); } // init diff --git a/Source/Microphysics/Morrison/ERF_Morrison.H b/Source/Microphysics/Morrison/ERF_Morrison.H index a4d22b3d0..09211084b 100644 --- a/Source/Microphysics/Morrison/ERF_Morrison.H +++ b/Source/Microphysics/Morrison/ERF_Morrison.H @@ -72,7 +72,7 @@ public: { m_axis = sc.ave_plane; m_rdOcp = sc.rdOcp; - m_do_cond = (!sc.use_shoc); + m_do_cond = (!sc.use_shoc && !sc.macrophysics_acoustic_coupling); } // init diff --git a/Source/Microphysics/SAM/ERF_SAM.H b/Source/Microphysics/SAM/ERF_SAM.H index 637f14ab1..4cbfa5054 100644 --- a/Source/Microphysics/SAM/ERF_SAM.H +++ b/Source/Microphysics/SAM/ERF_SAM.H @@ -83,7 +83,7 @@ public: m_fac_sub = lsub / sc.c_p; m_axis = sc.ave_plane; m_rdOcp = sc.rdOcp; - m_do_cond = (!sc.use_shoc); + m_do_cond = (!sc.use_shoc && !sc.macrophysics_acoustic_coupling); } // init diff --git a/Source/Microphysics/SatAdj/ERF_SatAdj.H b/Source/Microphysics/SatAdj/ERF_SatAdj.H index 707dbcc14..02393bd75 100644 --- a/Source/Microphysics/SatAdj/ERF_SatAdj.H +++ b/Source/Microphysics/SatAdj/ERF_SatAdj.H @@ -58,7 +58,7 @@ public: { m_fac_cond = lcond / sc.c_p; m_rdOcp = sc.rdOcp; - m_do_cond = (!sc.use_shoc); + m_do_cond = (!sc.use_shoc && !sc.macrophysics_acoustic_coupling); } // init diff --git a/Source/Microphysics/WSM6/ERF_WSM6.H b/Source/Microphysics/WSM6/ERF_WSM6.H index 26b0d9ce6..0478b6a5b 100644 --- a/Source/Microphysics/WSM6/ERF_WSM6.H +++ b/Source/Microphysics/WSM6/ERF_WSM6.H @@ -44,7 +44,7 @@ public: void Define(SolverChoice& sc) override { m_axis = sc.ave_plane; - m_do_cond = (!sc.use_shoc); + m_do_cond = (!sc.use_shoc && !sc.macrophysics_acoustic_coupling); } void Init(const amrex::MultiFab& cons_in, diff --git a/Source/SourceTerms/ERF_MakeSources.cpp b/Source/SourceTerms/ERF_MakeSources.cpp index 634b2dff8..baed75591 100644 --- a/Source/SourceTerms/ERF_MakeSources.cpp +++ b/Source/SourceTerms/ERF_MakeSources.cpp @@ -4,10 +4,10 @@ #include #include -#include -#include -#include -#include +#include "ERF_NumericalDiffusion.H" +#include "ERF_SrcHeaders.H" +#include "ERF_TI_slow_headers.H" +#include "ERF_MOSTStress.H" using namespace amrex; @@ -46,6 +46,7 @@ void make_sources (int level, const Geometry geom, const SolverChoice& solverChoice, Vector>& mapfac, + const MultiFab* macrophysics_src, const MultiFab* rhotheta_src, const MultiFab* rhoqt_src, const Real* dptr_wbar_sub, @@ -67,7 +68,7 @@ void make_sources (int level, source.setVal(0.0,Rho_comp,2); } - const bool l_use_ndiff = solverChoice.use_num_diff; + const bool l_use_ndiff = solverChoice.use_num_diff; TurbChoice tc = solverChoice.turbChoice[level]; const bool l_use_KE = tc.use_tke; @@ -92,6 +93,9 @@ void make_sources (int level, // flag for a moisture model bool has_moisture = (solverChoice.moisture_type != MoistureType::None); + // flag for tight macrophysics coupling + bool acoustic_macro = solverChoice.macrophysics_acoustic_coupling; + // ***************************************************************************** // Planar averages for subsidence terms // ***************************************************************************** @@ -194,18 +198,19 @@ void make_sources (int level, // ***************************************************************************** // Define source term for cell-centered conserved variables, from - // one user-defined source terms for (rho theta) and (rho q_t) - // two radiation for (rho theta) - // three Rayleigh damping for (rho theta) - // Real(4.) custom forcing for (rho theta) and (rho Q1) - // Real(5.) custom subsidence for (rho theta) and (rho Q1) - // Real(6.) numerical diffusion for (rho theta) - // Real(7.) sponging - // Real(8.) turbulent perturbation - // Real(9.) nudging towards input sounding values (only for theta) + // 1. user-defined source terms for (rho theta) and (rho q_t) + // 2. radiation for (rho theta) + // 3. Rayleigh damping for (rho theta) + // 4. custom forcing for (rho theta) and (rho Q1) + // 5. custom subsidence for (rho theta) and (rho Q1) + // 6. numerical diffusion for (rho theta) + // 7. sponging + // 8. turbulent perturbation + // 9. nudging towards input sounding values (only for theta) // 10a. Immersed forcing for terrain // 10b. Immersed forcing for buildings - // Real(11.) Four stream radiation source for (rho theta) + // 11. Four stream radiation source for (rho theta) + // 12. Macrophysics source for (rho theta & rho Q1) // ***************************************************************************** // *********************************************************************************************** @@ -223,18 +228,20 @@ void make_sources (int level, const Array4& cell_prim = S_prim.array(mfi); const Array4 & cell_src = source.array(mfi); - const Array4& r0 = r_hse.const_array(mfi); + const Array4& mac_src = (macrophysics_src) ? macrophysics_src->const_array(mfi) : + Array4{}; + + const Array4& r0 = r_hse.const_array(mfi); const Array4& th0 = th_hse.const_array(mfi); const Array4& qv0 = qv_hse.const_array(mfi); const Array4& z_cc_arr = z_phys_cc->const_array(mfi); const Array4& t_blank_arr = (terrain_blank) ? terrain_blank->const_array(mfi) : - Array4{}; - + Array4{}; // ************************************************************************************* - // two Add radiation source terms to (rho theta) + // 2. Add radiation source terms to (rho theta) // ************************************************************************************* if (solverChoice.rad_type != RadiationType::None && is_slow_step) { auto const& qheating_arr = qheating_rates->const_array(mfi); @@ -247,7 +254,7 @@ void make_sources (int level, // ************************************************************************************* - // three Add Rayleigh damping for (rho theta) + // 3. Add Rayleigh damping for (rho theta) // ************************************************************************************* Real dampcoef = solverChoice.dampingChoice.rayleigh_dampcoef; @@ -267,7 +274,7 @@ void make_sources (int level, } // ************************************************************************************* - // Real(4.) Add custom forcing for (rho theta) + // 4. Add custom forcing for (rho theta) // ************************************************************************************* if (solverChoice.custom_rhotheta_forcing && is_slow_step) { const int n = RhoTheta_comp; @@ -303,7 +310,7 @@ void make_sources (int level, } // ************************************************************************************* - // Real(4.) Add custom forcing for RhoQ1 + // 4. Add custom forcing for RhoQ1 // ************************************************************************************* if (solverChoice.custom_moisture_forcing && is_slow_step) { const int n = RhoQ1_comp; @@ -339,7 +346,7 @@ void make_sources (int level, } // ************************************************************************************* - // Real(5.) Add custom subsidence for (rho theta) + // 5. Add custom subsidence for (rho theta) // ************************************************************************************* if (solverChoice.custom_w_subsidence && is_slow_step && solverChoice.do_theta_advection) { const int n = RhoTheta_comp; @@ -366,7 +373,7 @@ void make_sources (int level, } // ************************************************************************************* - // Real(5.) Add custom subsidence for RhoQ1 and RhoQ2 + // 5. Add custom subsidence for RhoQ1 and RhoQ2 // ************************************************************************************* if (solverChoice.custom_w_subsidence && (solverChoice.moisture_type != MoistureType::None) && is_slow_step) { const int nv = RhoQ1_comp; @@ -399,7 +406,7 @@ void make_sources (int level, } // ************************************************************************************* - // Real(6.) Add numerical diffusion for rho and (rho theta) + // 6. Add numerical diffusion for rho and (rho theta) // ************************************************************************************* if (l_use_ndiff && is_slow_step) { int sc; @@ -427,7 +434,7 @@ void make_sources (int level, } // ************************************************************************************* - // Real(7.) Add sponging + // 7. Add sponging // ************************************************************************************* if (!(solverChoice.spongeChoice.sponge_type == SpongeType::Input_Sponge) && is_slow_step){ const int n_qstate = S_data[IntVars::cons].nComp() - (NDRY + NSCALARS); @@ -441,7 +448,7 @@ void make_sources (int level, // ************************************************************************************* - // Real(8.) Add perturbation + // 8. Add perturbation // ************************************************************************************* if (solverChoice.use_source_perturbation(level) && is_slow_step) { auto m_ixtype = S_data[IntVars::cons].boxArray().ixType(); // Conserved term @@ -450,7 +457,7 @@ void make_sources (int level, } // ************************************************************************************* - // Real(9.) Add nudging towards value specified in input sounding + // 9. Add nudging towards value specified in input sounding // ************************************************************************************* if (solverChoice.nudging_from_input_sounding && is_slow_step) { @@ -679,7 +686,7 @@ void make_sources (int level, } // ************************************************************************************* - // Real(11.) Add 4 stream radiation src to RhoTheta + // 11. Add 4 stream radiation src to RhoTheta // ************************************************************************************* if (solverChoice.four_stream_radiation && has_moisture && is_slow_step) { @@ -736,6 +743,18 @@ void make_sources (int level, } }); } + + // ************************************************************************************* + // 12. Add macrophysics source for RhoTheta + // ************************************************************************************* + if (is_slow_step && has_moisture && acoustic_macro) { + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + cell_src(i,j,k,RhoTheta_comp) += mac_src(i,j,k,0); + cell_src(i,j,k,RhoQ1_comp ) += mac_src(i,j,k,1); + cell_src(i,j,k,RhoQ2_comp ) += mac_src(i,j,k,2); + }); + } } // mfi } // OMP } diff --git a/Source/SourceTerms/ERF_SrcHeaders.H b/Source/SourceTerms/ERF_SrcHeaders.H index bbe5b29c0..8cd7caa09 100644 --- a/Source/SourceTerms/ERF_SrcHeaders.H +++ b/Source/SourceTerms/ERF_SrcHeaders.H @@ -69,6 +69,7 @@ void make_sources (int level, int nrk, const amrex::Geometry geom, const SolverChoice& solverChoice, amrex::Vector>& mapfac, + const amrex::MultiFab* macrophysics_src, const amrex::MultiFab* rhotheta_src, const amrex::MultiFab* rhoqt_src, const amrex::Real* dptr_wbar_sub, diff --git a/Source/TimeIntegration/ERF_Advance.cpp b/Source/TimeIntegration/ERF_Advance.cpp index 3e4857e20..c750f1be0 100644 --- a/Source/TimeIntegration/ERF_Advance.cpp +++ b/Source/TimeIntegration/ERF_Advance.cpp @@ -1,8 +1,9 @@ -#include -#include +#include "ERF.H" +#include "ERF_Utils.H" +#include "ERF_MacroPhysics.H" #ifdef ERF_USE_WINDFARM -#include +#include "ERF_WindFarm.H" #endif using namespace amrex; @@ -204,6 +205,15 @@ ERF::Advance (int lev, Real time, Real dt_lev, int iteration, int /*ncycle*/) state_new.push_back(MultiFab(rV_new[lev], amrex::make_alias, 0, 1)); // ymom state_new.push_back(MultiFab(rW_new[lev], amrex::make_alias, 0, 1)); // zmom + // ************************************************************************************** + // Macrophysics source term computed pre-step + // ************************************************************************************** + if (solverChoice.moisture_type != MoistureType::None && + solverChoice.macrophysics_acoustic_coupling) + { + GetMacroPhysicsSrc(real_width, dt_lev, Geom(lev), macrophysics_src[lev].get(), S_old); + } + // ************************************************************************************** // Tests on the reasonableness of the solution before the dycore // ************************************************************************************** diff --git a/Source/TimeIntegration/ERF_ImplicitPost.H b/Source/TimeIntegration/ERF_ImplicitPost.H index d700060fe..6ecd6b241 100644 --- a/Source/TimeIntegration/ERF_ImplicitPost.H +++ b/Source/TimeIntegration/ERF_ImplicitPost.H @@ -3,7 +3,8 @@ // ***************************************************************************** const bool l_do_implicit_moist = solverChoice.implicit_moisture_diffusion; const Real l_vert_implicit_fac = solverChoice.vert_implicit_fac[nrk]; - if ( (l_do_implicit_moist) && + if ( false && + (l_do_implicit_moist) && (l_vert_implicit_fac > zero) && (solverChoice.moisture_type != MoistureType::None) ) { @@ -34,19 +35,19 @@ const Array4& q1fx_z = Q1fx3->const_array(mfi); if (l_use_stretched_dz) { - ImplicitDiffForStateLU_S(bx, fine_geom.Domain(), level, RhoQ1_comp, + ImplicitDiffForStateLU_S(bx, fine_geom.Domain(), level, RhoQ1_comp, RhoQ1_comp, slow_dt, l_bc_neumann_vals_d, cell_data, stretched_dz_d[level], q1fx_z, mu_turb, solverChoice, bc_ptr_h, l_use_SurfLayer, l_vert_implicit_fac); } else if (l_use_terrain_fitted_coords) { - ImplicitDiffForStateLU_T(bx, fine_geom.Domain(), level, RhoQ1_comp, + ImplicitDiffForStateLU_T(bx, fine_geom.Domain(), level, RhoQ1_comp, RhoQ1_comp, slow_dt, l_bc_neumann_vals_d, cell_data, z_nd_arr, detJ_arr, dxInv, q1fx_z, mu_turb, solverChoice, bc_ptr_h, l_use_SurfLayer, l_vert_implicit_fac); } else { // no terrain - ImplicitDiffForStateLU_N(bx, fine_geom.Domain(), level, RhoQ1_comp, + ImplicitDiffForStateLU_N(bx, fine_geom.Domain(), level, RhoQ1_comp, RhoQ1_comp, slow_dt, l_bc_neumann_vals_d, cell_data, dxInv, q1fx_z, mu_turb, solverChoice, bc_ptr_h, l_use_SurfLayer, l_vert_implicit_fac); diff --git a/Source/TimeIntegration/ERF_ImplicitPre.H b/Source/TimeIntegration/ERF_ImplicitPre.H index a7a24348f..7dad86a8e 100644 --- a/Source/TimeIntegration/ERF_ImplicitPre.H +++ b/Source/TimeIntegration/ERF_ImplicitPre.H @@ -6,8 +6,10 @@ if (l_vert_implicit_fac > zero) { - const bool l_do_implicit_theta = solverChoice.implicit_thermal_diffusion; - const bool l_do_implicit_mom = solverChoice.implicit_momentum_diffusion; + const bool l_has_moisture = (solverChoice.moisture_type != MoistureType::None); + const bool l_do_implicit_theta = solverChoice.implicit_thermal_diffusion; + const bool l_do_implicit_moisture = solverChoice.implicit_moisture_diffusion; + const bool l_do_implicit_mom = solverChoice.implicit_momentum_diffusion; // If we're doing an implicit solve for momenta (u and v only), // then we use the explicit solution at this point to derive @@ -61,16 +63,24 @@ const Array4 tau13 = Tau[level][TauType::tau13]->array(mfi); const Array4 tau23 = Tau[level][TauType::tau23]->array(mfi); [[maybe_unused]] const Array4 tau33 = Tau[level][TauType::tau33]->array(mfi); - const Array4& hfx_z = Hfx3->const_array(mfi); + const Array4& hfx_z = Hfx3->const_array(mfi); + const Array4& q1fx_z = Q1fx3->const_array(mfi); if (l_use_stretched_dz) { if (l_do_implicit_theta) { - ImplicitDiffForStateLU_S(bx, fine_geom.Domain(), level, RhoTheta_comp, + ImplicitDiffForStateLU_S(bx, fine_geom.Domain(), level, 1, RhoTheta_comp, slow_dt, l_bc_neumann_vals_d, cell_data, stretched_dz_d[level], hfx_z, mu_turb, solverChoice, bc_ptr_h, l_use_SurfLayer, l_vert_implicit_fac); } + if (l_do_implicit_moisture && l_has_moisture) { + ImplicitDiffForStateLU_S(bx, fine_geom.Domain(), level, 2, RhoQ1_comp, + slow_dt, l_bc_neumann_vals_d, cell_data, + stretched_dz_d[level], q1fx_z, + mu_turb, solverChoice, + bc_ptr_h, l_use_SurfLayer, l_vert_implicit_fac); + } if (l_do_implicit_mom) { ImplicitDiffForMomLU_S<0>(tbx, fine_geom.Domain(), level, slow_dt, cell_data, rho_u, tau13, tau13_corr, @@ -93,12 +103,19 @@ } } else if (l_use_terrain_fitted_coords) { if (l_do_implicit_theta) { - ImplicitDiffForStateLU_T(bx, fine_geom.Domain(), level, RhoTheta_comp, + ImplicitDiffForStateLU_T(bx, fine_geom.Domain(), level, 1, RhoTheta_comp, slow_dt, l_bc_neumann_vals_d, cell_data, z_nd_arr, detJ_arr, dxInv, hfx_z, mu_turb, solverChoice, bc_ptr_h, l_use_SurfLayer, l_vert_implicit_fac); } + if (l_do_implicit_moisture && l_has_moisture) { + ImplicitDiffForStateLU_T(bx, fine_geom.Domain(), level, 2, RhoQ1_comp, + slow_dt, l_bc_neumann_vals_d, cell_data, + z_nd_arr, detJ_arr, dxInv, q1fx_z, + mu_turb, solverChoice, + bc_ptr_h, l_use_SurfLayer, l_vert_implicit_fac); + } if (l_do_implicit_mom) { ImplicitDiffForMomLU_T<0>(tbx, fine_geom.Domain(), level, slow_dt, cell_data, rho_u, tau13, tau13_corr, @@ -121,11 +138,17 @@ } } else { // no terrain if (l_do_implicit_theta) { - ImplicitDiffForStateLU_N(bx, fine_geom.Domain(), level, RhoTheta_comp, + ImplicitDiffForStateLU_N(bx, fine_geom.Domain(), level, 1, RhoTheta_comp, slow_dt, l_bc_neumann_vals_d, cell_data, dxInv, hfx_z, mu_turb, solverChoice, bc_ptr_h, l_use_SurfLayer, l_vert_implicit_fac); } + if (l_do_implicit_moisture && l_has_moisture) { + ImplicitDiffForStateLU_N(bx, fine_geom.Domain(), level, 2, RhoQ1_comp, + slow_dt, l_bc_neumann_vals_d, cell_data, + dxInv, q1fx_z, mu_turb, solverChoice, + bc_ptr_h, l_use_SurfLayer, l_vert_implicit_fac); + } if (l_do_implicit_mom) { ImplicitDiffForMomLU_N<0>(tbx, fine_geom.Domain(), level, slow_dt, cell_data, rho_u, tau13, tau13_corr, diff --git a/Source/TimeIntegration/ERF_MakeFastCoeffs.cpp b/Source/TimeIntegration/ERF_MakeFastCoeffs.cpp index b3e43acfd..531df69e4 100644 --- a/Source/TimeIntegration/ERF_MakeFastCoeffs.cpp +++ b/Source/TimeIntegration/ERF_MakeFastCoeffs.cpp @@ -123,16 +123,32 @@ void make_fast_coeffs (int /*level*/, Real detJ_on_kface = myhalf * (detJ(i,j,k) + detJ(i,j,k-1)); Real inv_detJ_on_kface = one / detJ_on_kface; - Real qv_p = (l_use_moisture) ? prim(i,j,k ,PrimQ1_comp) : zero; - Real qv_q = (l_use_moisture) ? prim(i,j,k-1,PrimQ1_comp) : zero; - - Real coeff_P = -Gamma * R_d * dzi * inv_detJ_on_kface * pi_c * (one + RvOverRd*qv_p) + Real thd_km2 = prim(i,j,k-2,PrimTheta_comp); + Real thd_km1 = prim(i,j,k-1,PrimTheta_comp); + Real thd_k = prim(i,j,k ,PrimTheta_comp); + Real thd_kp1 = prim(i,j,k+1,PrimTheta_comp); + + Real qv_km2 = (l_use_moisture) ? prim(i,j,k-2,PrimQ1_comp) : zero; + Real qv_km1 = (l_use_moisture) ? prim(i,j,k-1,PrimQ1_comp) : zero; + Real qv_k = (l_use_moisture) ? prim(i,j,k ,PrimQ1_comp) : zero; + Real qv_kp1 = (l_use_moisture) ? prim(i,j,k+1,PrimQ1_comp) : zero; + + Real thm_km2 = thd_km2 * (Real(1.0) + RvOverRd*qv_km2); + Real thm_km1 = thd_km1 * (Real(1.0) + RvOverRd*qv_km1); + Real thm_k = thd_k * (Real(1.0) + RvOverRd*qv_k ); + Real thm_kp1 = thd_kp1 * (Real(1.0) + RvOverRd*qv_kp1); + + Real thm_t_lo = myhalf * (thm_km2 + thm_km1); + Real thm_t_mid = myhalf * (thm_km1 + thm_k ); + Real thm_t_hi = myhalf * (thm_k + thm_kp1); + + Real coeff_P = -Gamma * R_d * dzi * inv_detJ_on_kface * pi_c + halfg * R_d * rhobar_hi * pi_stage_ca(i,j,k) / - ( c_v * pibar_hi * stage_cons(i,j,k,RhoTheta_comp) ); + ( c_v * pibar_hi * stage_cons(i,j,k,RhoTheta_comp) * (one + RvOverRd*qv_k) ); - Real coeff_Q = Gamma * R_d * dzi * inv_detJ_on_kface * pi_c * (one + RvOverRd*qv_q) + Real coeff_Q = Gamma * R_d * dzi * inv_detJ_on_kface * pi_c + halfg * R_d * rhobar_lo * pi_stage_ca(i,j,k-1) / - ( c_v * pibar_lo * stage_cons(i,j,k-1,RhoTheta_comp) ); + ( c_v * pibar_lo * stage_cons(i,j,k-1,RhoTheta_comp) * (one + RvOverRd*qv_km1) ); coeffP_a(i,j,k) = coeff_P; coeffQ_a(i,j,k) = coeff_Q; @@ -144,17 +160,14 @@ void make_fast_coeffs (int /*level*/, coeff_Q /= (one + q); } - Real theta_t_lo = myhalf * ( prim(i,j,k-2,PrimTheta_comp) + prim(i,j,k-1,PrimTheta_comp) ); - Real theta_t_mid = myhalf * ( prim(i,j,k-1,PrimTheta_comp) + prim(i,j,k ,PrimTheta_comp) ); - Real theta_t_hi = myhalf * ( prim(i,j,k ,PrimTheta_comp) + prim(i,j,k+1,PrimTheta_comp) ); - // LHS for tri-diagonal system Real D = dtau * dtau * beta_2 * beta_2 * dzi; - coeffA_a(i,j,k) = D * (one/detJ(i,j,k-1)) * ( halfg - coeff_Q * theta_t_lo ); - coeffC_a(i,j,k) = D * (one/detJ(i,j,k )) * (-halfg + coeff_P * theta_t_hi ); + coeffA_a(i,j,k) = D * (one/detJ(i,j,k-1)) * ( halfg - coeff_Q * thm_t_lo ); + coeffC_a(i,j,k) = D * (one/detJ(i,j,k )) * (-halfg + coeff_P * thm_t_hi ); - coeffB_a(i,j,k) = one + D * ( (coeff_Q/detJ(i,j,k-1) - coeff_P/detJ(i,j,k)) * theta_t_mid - + halfg * (Real(1.0)/detJ(i,j,k) - Real(1.0)/detJ(i,j,k-1)) ); + coeffB_a(i,j,k) = one + D * ( ( coeff_Q / detJ(i,j,k-1) + - coeff_P / detJ(i,j,k ) ) * thm_t_mid + + halfg * ( Real(1.0)/detJ(i,j,k) - Real(1.0)/detJ(i,j,k-1) ) ); }); } else { @@ -172,13 +185,13 @@ void make_fast_coeffs (int /*level*/, Real qv_p = (l_use_moisture) ? prim(i,j,k ,PrimQ1_comp) : zero; Real qv_q = (l_use_moisture) ? prim(i,j,k-1,PrimQ1_comp) : zero; - Real coeff_P = -Gamma * R_d * dzi * pi_c * (one + RvOverRd*qv_p) + Real coeff_P = -Gamma * R_d * dzi * pi_c + halfg * R_d * rhobar_hi * pi_stage_ca(i,j,k) / - ( c_v * pibar_hi * stage_cons(i,j,k,RhoTheta_comp) ); + ( c_v * pibar_hi * stage_cons(i,j,k,RhoTheta_comp) * (one + RvOverRd*qv_p) ); - Real coeff_Q = Gamma * R_d * dzi * pi_c * (one + RvOverRd*qv_q) + Real coeff_Q = Gamma * R_d * dzi * pi_c + halfg * R_d * rhobar_lo * pi_stage_ca(i,j,k-1) / - ( c_v * pibar_lo * stage_cons(i,j,k-1,RhoTheta_comp) ); + ( c_v * pibar_lo * stage_cons(i,j,k-1,RhoTheta_comp) * (one + RvOverRd*qv_q) ); coeffP_a(i,j,k) = coeff_P; coeffQ_a(i,j,k) = coeff_Q; @@ -196,10 +209,11 @@ void make_fast_coeffs (int /*level*/, // LHS for tri-diagonal system Real D = dtau * dtau * beta_2 * beta_2 * dzi; - coeffA_a(i,j,k) = D * ( halfg - coeff_Q * theta_t_lo ); - coeffC_a(i,j,k) = D * (-halfg + coeff_P * theta_t_hi ); + coeffA_a(i,j,k) = D * ( halfg - coeff_Q * (one + RvOverRd*qv_q) * theta_t_lo ); + coeffC_a(i,j,k) = D * (-halfg + coeff_P * (one + RvOverRd*qv_p) * theta_t_hi ); - coeffB_a(i,j,k) = one + D * (coeff_Q - coeff_P) * theta_t_mid; + coeffB_a(i,j,k) = one + D * ( coeff_Q * (one + RvOverRd*qv_q) + - coeff_P * (one + RvOverRd*qv_p) ) * theta_t_mid; }); } diff --git a/Source/TimeIntegration/ERF_SlowRhsPost.cpp b/Source/TimeIntegration/ERF_SlowRhsPost.cpp index 4d9b0bc67..2e97f8700 100644 --- a/Source/TimeIntegration/ERF_SlowRhsPost.cpp +++ b/Source/TimeIntegration/ERF_SlowRhsPost.cpp @@ -175,11 +175,11 @@ void erf_slow_rhs_post (int level, int finest_level, } // Valid vars - Vector is_valid_slow_var; is_valid_slow_var.resize(RhoQ1_comp+1,0); + Vector is_valid_slow_var; is_valid_slow_var.resize(RhoQ2_comp+1,0); if (l_use_KE) { is_valid_slow_var[ RhoKE_comp] = 1; } if (l_do_scalar) { is_valid_slow_var[RhoScalar_comp] = 1; } if (solverChoice.moisture_type != MoistureType::None) { - is_valid_slow_var[RhoQ1_comp] = 1; + is_valid_slow_var[RhoQ2_comp] = 1; } // ************************************************************************* @@ -281,12 +281,12 @@ void erf_slow_rhs_post (int level, int finest_level, const GpuArray ncomp_slow = {nsv,0,0,0}; // ************************************************************************** - // Note that here we do copy only the "slow" variables, not (rho) or (rho theta) + // Note that here we do copy only the "slow" variables, not (rho, rhoTheta, RhoQ1) // ************************************************************************** ParallelFor(tbx, ncomp_slow[IntVars::cons], [=] AMREX_GPU_DEVICE (int i, int j, int k, int nn) { const int n = scomp_slow[IntVars::cons] + nn; - cur_cons(i,j,k,n) = new_cons(i,j,k,n); + if (n != RhoQ1_comp) { cur_cons(i,j,k,n) = new_cons(i,j,k,n); } }); // We have projected the velocities stored in S_data but we will use @@ -372,14 +372,14 @@ void erf_slow_rhs_post (int level, int finest_level, // // Note that we either advect and diffuse all or none of the moisture variables // - for (int ivar(RhoKE_comp); ivar<= RhoQ1_comp; ++ivar) + for (int ivar(RhoKE_comp); ivar<= RhoQ2_comp; ++ivar) { if (is_valid_slow_var[ivar]) { start_comp = ivar; num_comp = 1; - if (ivar == RhoQ1_comp) { + if (ivar == RhoQ2_comp) { horiz_adv_type = ac.moistscal_horiz_adv_type; vert_adv_type = ac.moistscal_vert_adv_type; horiz_upw_frac = ac.moistscal_horiz_upw_frac; @@ -390,7 +390,7 @@ void erf_slow_rhs_post (int level, int finest_level, vert_adv_type = EfficientAdvType(nrk,ac.moistscal_vert_adv_type); } - num_comp = n_qstate; + num_comp = n_qstate - 1; } else { horiz_adv_type = ac.dryscal_horiz_adv_type; @@ -436,7 +436,7 @@ void erf_slow_rhs_post (int level, int finest_level, if (l_use_diff) { // Allow for implicit moisture diffusion - const Real l_vert_implicit_fac = (solverChoice.implicit_moisture_diffusion) ? + const Real l_vert_implicit_fac = (false) ? //(solverChoice.implicit_moisture_diffusion) ? solverChoice.vert_implicit_fac[nrk] : zero; const Array4 tm_arr = t_mean_mf ? t_mean_mf->const_array(mfi) : Array4{}; @@ -477,7 +477,7 @@ void erf_slow_rhs_post (int level, int finest_level, } // loop ivar #if defined(ERF_USE_NETCDF) - if (moist_set_rhs_bool) + if (false) //moist_set_rhs_bool) { Real bdy_factor = solverChoice.bdy_nudge_factor; const Array4 & new_cons_const = S_new[IntVars::cons].const_array(mfi); @@ -507,14 +507,14 @@ void erf_slow_rhs_post (int level, int finest_level, auto const& src_arr = source.const_array(mfi); - for (int ivar(RhoKE_comp); ivar<= RhoQ1_comp; ++ivar) + for (int ivar(RhoKE_comp); ivar<= RhoQ2_comp; ++ivar) { if (is_valid_slow_var[ivar]) { start_comp = ivar; num_comp = 1; - if (ivar == RhoQ1_comp) { - num_comp = nvars - RhoQ1_comp; + if (ivar == RhoQ2_comp) { + num_comp = nvars - RhoQ2_comp; } else if (ivar == RhoScalar_comp) { num_comp = NSCALARS; } diff --git a/Source/TimeIntegration/ERF_SlowRhsPre.cpp b/Source/TimeIntegration/ERF_SlowRhsPre.cpp index 9311d3a33..4d1471caf 100644 --- a/Source/TimeIntegration/ERF_SlowRhsPre.cpp +++ b/Source/TimeIntegration/ERF_SlowRhsPre.cpp @@ -131,6 +131,12 @@ void erf_slow_rhs_pre (int level, int finest_level, const AdvType l_vert_adv_type = solverChoice.advChoice.dycore_vert_adv_type; const Real l_horiz_upw_frac = solverChoice.advChoice.dycore_horiz_upw_frac; const Real l_vert_upw_frac = solverChoice.advChoice.dycore_vert_upw_frac; + + const AdvType l_moist_horiz_adv_type = solverChoice.advChoice.moistscal_horiz_adv_type; + const AdvType l_moist_vert_adv_type = solverChoice.advChoice.moistscal_vert_adv_type; + const Real l_moist_horiz_upw_frac = solverChoice.advChoice.moistscal_horiz_upw_frac; + const Real l_moist_vert_upw_frac = solverChoice.advChoice.moistscal_vert_upw_frac; + const bool l_use_stretched_dz = (solverChoice.mesh_type == MeshType::StretchedDz); const bool l_use_terrain_fitted_coords = (solverChoice.mesh_type == MeshType::VariableDz); const bool l_moving_terrain = (solverChoice.terrain_type == TerrainType::MovingFittedMesh); @@ -562,32 +568,53 @@ void erf_slow_rhs_pre (int level, int finest_level, dxInv, mf_mx, mf_my, mf_uy, mf_vx, flx_arr, l_fixed_rho); AdvectionSrcForScalars(bx, icomp, ncomp, - avg_xmom_arr, avg_ymom_arr, avg_zmom_arr, - cell_prim, cell_rhs, - detJ_arr, dxInv, mf_mx, mf_my, - l_horiz_adv_type, l_vert_adv_type, - l_horiz_upw_frac, l_vert_upw_frac, - flx_arr, domain, bc_ptr_h); + avg_xmom_arr, avg_ymom_arr, avg_zmom_arr, + cell_prim, cell_rhs, + detJ_arr, dxInv, mf_mx, mf_my, + l_horiz_adv_type, l_vert_adv_type, + l_horiz_upw_frac, l_vert_upw_frac, + flx_arr, domain, bc_ptr_h); + if (l_use_moisture) { + AdvectionSrcForScalars(bx, RhoQ1_comp, ncomp, + avg_xmom_arr, avg_ymom_arr, avg_zmom_arr, + cell_prim, cell_rhs, + detJ_arr, dxInv, mf_mx, mf_my, + l_moist_horiz_adv_type, l_moist_vert_adv_type, + l_moist_horiz_upw_frac, l_moist_vert_upw_frac, + flx_arr, domain, bc_ptr_h); + } } else { EBAdvectionSrcForRho(bx, cell_rhs, - rho_u, rho_v, omega_arr, - avg_xmom_arr, avg_ymom_arr, avg_zmom_arr, - mask_arr, cfg_arr, - ax_arr, ay_arr, az_arr, - fcx_arr, fcy_arr, fcz_arr, detJ_arr, - dxInv, mf_mx, mf_my, mf_uy, mf_vx, - flx_arr, l_fixed_rho, - already_on_centroids); + rho_u, rho_v, omega_arr, + avg_xmom_arr, avg_ymom_arr, avg_zmom_arr, + mask_arr, cfg_arr, + ax_arr, ay_arr, az_arr, + fcx_arr, fcy_arr, fcz_arr, detJ_arr, + dxInv, mf_mx, mf_my, mf_uy, mf_vx, + flx_arr, l_fixed_rho, + already_on_centroids); EBAdvectionSrcForScalars(bx, icomp, ncomp, - avg_xmom_arr, avg_ymom_arr, avg_zmom_arr, - cell_prim, cell_rhs, - mask_arr, cfg_arr, ax_arr, ay_arr, az_arr, - fcx_arr, fcy_arr, fcz_arr, - detJ_arr, dxInv, mf_mx, mf_my, - l_horiz_adv_type, l_vert_adv_type, - l_horiz_upw_frac, l_vert_upw_frac, - flx_arr, domain, bc_ptr_h, - already_on_centroids); + avg_xmom_arr, avg_ymom_arr, avg_zmom_arr, + cell_prim, cell_rhs, + mask_arr, cfg_arr, ax_arr, ay_arr, az_arr, + fcx_arr, fcy_arr, fcz_arr, + detJ_arr, dxInv, mf_mx, mf_my, + l_horiz_adv_type, l_vert_adv_type, + l_horiz_upw_frac, l_vert_upw_frac, + flx_arr, domain, bc_ptr_h, + already_on_centroids); + if (l_use_moisture) { + EBAdvectionSrcForScalars(bx, RhoQ1_comp, ncomp, + avg_xmom_arr, avg_ymom_arr, avg_zmom_arr, + cell_prim, cell_rhs, + mask_arr, cfg_arr, ax_arr, ay_arr, az_arr, + fcx_arr, fcy_arr, fcz_arr, + detJ_arr, dxInv, mf_mx, mf_my, + l_moist_horiz_adv_type, l_moist_vert_adv_type, + l_moist_horiz_upw_frac, l_moist_vert_upw_frac, + flx_arr, domain, bc_ptr_h, + already_on_centroids); + } } if (l_use_diff) { @@ -629,6 +656,17 @@ void erf_slow_rhs_pre (int level, int finest_level, hfx_z, q1fx_z, q2fx_z, diss, mu_turb, solverChoice, level, tm_arr, grav_gpu, bc_ptr_d, l_use_SurfLayer, l_vert_implicit_fac); + if (l_use_moisture) { + DiffusionSrcForState_S(bx, domain, RhoQ1_comp, 1, u, v, + cell_data, cell_prim, cell_rhs, + diffflux_x, diffflux_y, diffflux_z, + stretched_dz_d, dxInv, SmnSmn_a, + mf_mx, mf_ux, mf_vx, + mf_my, mf_uy, mf_vy, + hfx_z, q1fx_z, q2fx_z, diss, + mu_turb, solverChoice, level, + tm_arr, grav_gpu, bc_ptr_d, l_use_SurfLayer, l_vert_implicit_fac); + } } else if (l_use_terrain_fitted_coords) { DiffusionSrcForState_T(bx, domain, n_start, n_comp, l_rotate, u, v, cell_data, cell_prim, cell_rhs, @@ -640,6 +678,18 @@ void erf_slow_rhs_pre (int level, int finest_level, hfx_x, hfx_y, hfx_z, q1fx_x, q1fx_y, q1fx_z, q2fx_z, diss, mu_turb, solverChoice, level, tm_arr, grav_gpu, bc_ptr_d, l_use_SurfLayer, l_vert_implicit_fac); + if (l_use_moisture) { + DiffusionSrcForState_T(bx, domain, RhoQ1_comp, 1, l_rotate, u, v, + cell_data, cell_prim, cell_rhs, + diffflux_x, diffflux_y, diffflux_z, + z_nd, z_cc, ax_arr, ay_arr, az_arr, detJ_arr, + dxInv, SmnSmn_a, + mf_mx, mf_ux, mf_vx, + mf_my, mf_uy, mf_vy, + hfx_x, hfx_y, hfx_z, q1fx_x, q1fx_y, q1fx_z, q2fx_z, diss, + mu_turb, solverChoice, level, + tm_arr, grav_gpu, bc_ptr_d, l_use_SurfLayer, l_vert_implicit_fac); + } } else if (l_use_eb) { DiffusionSrcForState_EB(bx, domain, n_start, n_comp, u, v, cell_data, cell_prim, cell_rhs, @@ -650,6 +700,17 @@ void erf_slow_rhs_pre (int level, int finest_level, hfx_z, q1fx_z, q2fx_z, hfx_EB, mu_turb, solverChoice, level, bc_ptr_d, l_use_SurfLayer); + if (l_use_moisture) { + DiffusionSrcForState_EB(bx, domain, RhoQ1_comp, 1, u, v, + cell_data, cell_prim, cell_rhs, + diffflux_x, diffflux_y, diffflux_z, + cfg_arr, ax_arr, ay_arr, az_arr, detJ_arr, + barea_arr, bcent_arr, + dx, dxInv, + hfx_z, q1fx_z, q2fx_z, hfx_EB, + mu_turb, solverChoice, level, + bc_ptr_d, l_use_SurfLayer); + } } else { DiffusionSrcForState_N(bx, domain, n_start, n_comp, u, v, cell_data, cell_prim, cell_rhs, @@ -660,6 +721,17 @@ void erf_slow_rhs_pre (int level, int finest_level, hfx_z, q1fx_z, q2fx_z, diss, mu_turb, solverChoice, level, tm_arr, grav_gpu, bc_ptr_d, l_use_SurfLayer, l_vert_implicit_fac); + if (l_use_moisture) { + DiffusionSrcForState_N(bx, domain, RhoQ1_comp, 1, u, v, + cell_data, cell_prim, cell_rhs, + diffflux_x, diffflux_y, diffflux_z, + dxInv, SmnSmn_a, + mf_mx, mf_ux, mf_vx, + mf_my, mf_uy, mf_vy, + hfx_z, q1fx_z, q2fx_z, diss, + mu_turb, solverChoice, level, + tm_arr, grav_gpu, bc_ptr_d, l_use_SurfLayer, l_vert_implicit_fac); + } } } @@ -668,6 +740,9 @@ void erf_slow_rhs_pre (int level, int finest_level, { cell_rhs(i,j,k,Rho_comp) += source_arr(i,j,k,Rho_comp); cell_rhs(i,j,k,RhoTheta_comp) += source_arr(i,j,k,RhoTheta_comp); + if (l_use_moisture) { + cell_rhs(i,j,k,RhoQ1_comp) += source_arr(i,j,k,RhoQ1_comp); + } }); // If anelastic and in second RK stage, take average of old-time and new-time source diff --git a/Source/TimeIntegration/ERF_Substep_T.cpp b/Source/TimeIntegration/ERF_Substep_T.cpp index fd6b59676..f1e338056 100644 --- a/Source/TimeIntegration/ERF_Substep_T.cpp +++ b/Source/TimeIntegration/ERF_Substep_T.cpp @@ -112,7 +112,6 @@ void erf_substep_T (int step, int /*nrk*/, MultiFab coeff_P_mf(fast_coeffs, make_alias, 3, 1); MultiFab coeff_Q_mf(fast_coeffs, make_alias, 4, 1); - // ************************************************************************* // Set gravity as a vector const Array grav{zero, zero, -gravity}; const GpuArray grav_gpu{grav[0], grav[1], grav[2]}; @@ -122,13 +121,12 @@ void erf_substep_T (int step, int /*nrk*/, // ************************************************************************* // First set up some arrays we'll need // ************************************************************************* - #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif for ( MFIter mfi(S_stage_data[IntVars::cons],TilingIfNotGPU()); mfi.isValid(); ++mfi) { - const Array4 & cur_cons = S_data[IntVars::cons].array(mfi); + const Array4 & cur_cons = S_data[IntVars::cons].array(mfi); const Array4& prev_cons = S_prev[IntVars::cons].const_array(mfi); const Array4& stage_cons = S_stage_data[IntVars::cons].const_array(mfi); const Array4& lagged_arr = lagged_delta_rt.array(mfi); @@ -137,7 +135,7 @@ void erf_substep_T (int step, int /*nrk*/, const Array4& old_drho_u = Delta_rho_u.array(mfi); const Array4& old_drho_v = Delta_rho_v.array(mfi); const Array4& old_drho_w = Delta_rho_w.array(mfi); - const Array4& old_drho_theta = Delta_rho_theta.array(mfi); + const Array4& old_drho_thm = Delta_rho_theta.array(mfi); const Array4& prev_xmom = S_prev[IntVars::xmom].const_array(mfi); const Array4& prev_ymom = S_prev[IntVars::ymom].const_array(mfi); @@ -150,13 +148,6 @@ void erf_substep_T (int step, int /*nrk*/, Box bx = mfi.validbox(); Box gbx = mfi.tilebox(); gbx.grow(1); - if (step == 0) { - ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - cur_cons(i,j,k,Rho_comp) = prev_cons(i,j,k,Rho_comp); - cur_cons(i,j,k,RhoTheta_comp) = prev_cons(i,j,k,RhoTheta_comp); - }); - } // step = 0 - Box gtbx = mfi.nodaltilebox(0); gtbx.grow(IntVect(1,1,0)); Box gtby = mfi.nodaltilebox(1); gtby.grow(IntVect(1,1,0)); Box gtbz = mfi.nodaltilebox(2); gtbz.grow(IntVect(1,1,0)); @@ -164,6 +155,20 @@ void erf_substep_T (int step, int /*nrk*/, const auto& bx_lo = lbound(bx); const auto& bx_hi = ubound(bx); + // ************************************************************************* + // Transfer from S_old to S_data in first step(update in place in S_data) + // ************************************************************************* + if (step == 0) { + ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + cur_cons(i,j,k,Rho_comp) = prev_cons(i,j,k,Rho_comp); + cur_cons(i,j,k,RhoTheta_comp) = prev_cons(i,j,k,RhoTheta_comp); + if (l_use_moisture) { cur_cons(i,j,k,RhoQ1_comp) = prev_cons(i,j,k,RhoQ1_comp); } + }); + } + + // ************************************************************************* + // Compute perturbational lateral momenta + // ************************************************************************* ParallelFor(gtbx, gtby, gtbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { old_drho_u(i,j,k) = prev_xmom(i,j,k) - stage_xmom(i,j,k); @@ -185,43 +190,29 @@ void erf_substep_T (int step, int /*nrk*/, old_drho_w(i,j,k) = prev_zmom(i,j,k) - stage_zmom(i,j,k); }); - const Array4& theta_extrap = extrap.array(mfi); - const Array4& prim = S_stage_prim.const_array(mfi); - + // ************************************************************************* + // Compute pert density, moist potential temperature, and thm extrap + // ************************************************************************* + const Array4& thm_extrap = extrap.array(mfi); + const Array4& prim = S_stage_prim.const_array(mfi); ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - old_drho(i,j,k) = cur_cons(i,j,k,Rho_comp) - stage_cons(i,j,k,Rho_comp); - old_drho_theta(i,j,k) = cur_cons(i,j,k,RhoTheta_comp) - stage_cons(i,j,k,RhoTheta_comp); + Real qv_stg = (l_use_moisture) ? prim(i,j,k,PrimQ1_comp) : zero; + Real qv_cur = (l_use_moisture) ? cur_cons(i,j,k,RhoQ1_comp ) / cur_cons(i,j,k,Rho_comp) : zero; + old_drho(i,j,k) = cur_cons(i,j,k,Rho_comp) - stage_cons(i,j,k,Rho_comp); + old_drho_thm(i,j,k) = cur_cons(i,j,k,RhoTheta_comp) * (Real(1.0) + RvOverRd*qv_cur) + - stage_cons(i,j,k,RhoTheta_comp) * (Real(1.0) + RvOverRd*qv_stg); if (step == 0) { - theta_extrap(i,j,k) = old_drho_theta(i,j,k); + thm_extrap(i,j,k) = old_drho_thm(i,j,k); } else { - theta_extrap(i,j,k) = old_drho_theta(i,j,k) + beta_d * - ( old_drho_theta(i,j,k) - lagged_arr(i,j,k) ); + thm_extrap(i,j,k) = old_drho_thm(i,j,k) + beta_d * ( old_drho_thm(i,j,k) - lagged_arr(i,j,k) ); } - - // NOTE: qv is not changing over the fast steps so we use the stage data - Real qv = (l_use_moisture) ? prim(i,j,k,PrimQ1_comp) : zero; - theta_extrap(i,j,k) *= (one + RvOverRd*qv); + lagged_arr(i,j,k) = old_drho_thm(i,j,k); }); } // mfi -#ifdef _OPENMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for ( MFIter mfi(S_stage_data[IntVars::cons],TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - // We define lagged_delta_rt for our next step as the current delta_rt - Box gbx = mfi.tilebox(); gbx.grow(1); - const Array4& old_drho_theta = Delta_rho_theta.array(mfi); - const Array4& lagged_arr = lagged_delta_rt.array(mfi); - ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - lagged_arr(i,j,k) = old_drho_theta(i,j,k); - }); - } // mfi - - // ************************************************************************* - // Define updates in the current RK stage - // ************************************************************************* - + // ********************************************************************* + // Update lateral momenta (explicit) + // ********************************************************************* #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif @@ -250,7 +241,6 @@ void erf_substep_T (int step, int /*nrk*/, const Array4& cur_xmom = S_data[IntVars::xmom].array(mfi); const Array4& cur_ymom = S_data[IntVars::ymom].array(mfi); - // These store the advection momenta which we will use to update the slow variables const Array4& avg_xmom_arr = avg_xmom.array(mfi); const Array4& avg_ymom_arr = avg_ymom.array(mfi); @@ -258,21 +248,11 @@ void erf_substep_T (int step, int /*nrk*/, const Array4& pi_stage_ca = pi_stage.const_array(mfi); - const Array4& theta_extrap = extrap.array(mfi); + const Array4& thm_extrap = extrap.array(mfi); - // Map factors const Array4& mf_ux = mapfac[MapFacType::u_x]->const_array(mfi); const Array4& mf_vy = mapfac[MapFacType::v_y]->const_array(mfi); - // Create old_drho_u/v/w/theta = U'', V'', W'', Theta'' in the docs - // Note that we do the Copy and Subtract including one ghost cell - // so that we don't have to fill ghost cells of the new MultiFabs - // Initialize New_rho_u/v/w to Delta_rho_u/v/w so that - // the ghost cells in New_rho_u/v/w will match old_drho_u/v/w - - // ********************************************************************* - // Define updates in the RHS of {x, y, z}-momentum equations - // ********************************************************************* { BL_PROFILE("substep_xymom_T"); @@ -285,12 +265,12 @@ void erf_substep_T (int step, int /*nrk*/, // Add (negative) gradient of (rho theta) multiplied by lagged "pi" Real met_h_xi = Compute_h_xi_AtIface (i, j, k, dxInv, z_nd); Real met_h_zeta = Compute_h_zeta_AtIface(i, j, k, dxInv, z_nd); - Real gp_xi = (theta_extrap(i,j,k) - theta_extrap(i-1,j,k)) * dxi; + Real gp_xi = (thm_extrap(i,j,k) - thm_extrap(i-1,j,k)) * dxi; Real gp_zeta_on_iface = (k == 0) ? - myhalf * dzi * ( theta_extrap(i-1,j,k+1) + theta_extrap(i,j,k+1) - - theta_extrap(i-1,j,k ) - theta_extrap(i,j,k ) ) : - fourth * dzi * ( theta_extrap(i-1,j,k+1) + theta_extrap(i,j,k+1) - - theta_extrap(i-1,j,k-1) - theta_extrap(i,j,k-1) ); + myhalf * dzi * ( thm_extrap(i-1,j,k+1) + thm_extrap(i,j,k+1) + - thm_extrap(i-1,j,k ) - thm_extrap(i,j,k ) ) : + fourth * dzi * ( thm_extrap(i-1,j,k+1) + thm_extrap(i,j,k+1) + - thm_extrap(i-1,j,k-1) - thm_extrap(i,j,k-1) ); Real gpx = gp_xi - (met_h_xi / met_h_zeta) * gp_zeta_on_iface; gpx *= mf_ux(i,j,0); @@ -317,12 +297,12 @@ void erf_substep_T (int step, int /*nrk*/, // Add (negative) gradient of (rho theta) multiplied by lagged "pi" Real met_h_eta = Compute_h_eta_AtJface(i, j, k, dxInv, z_nd); Real met_h_zeta = Compute_h_zeta_AtJface(i, j, k, dxInv, z_nd); - Real gp_eta = (theta_extrap(i,j,k) -theta_extrap(i,j-1,k)) * dyi; + Real gp_eta = (thm_extrap(i,j,k) - thm_extrap(i,j-1,k)) * dyi; Real gp_zeta_on_jface = (k == 0) ? - myhalf * dzi * ( theta_extrap(i,j,k+1) + theta_extrap(i,j-1,k+1) - - theta_extrap(i,j,k ) - theta_extrap(i,j-1,k ) ) : - fourth * dzi * ( theta_extrap(i,j,k+1) + theta_extrap(i,j-1,k+1) - - theta_extrap(i,j,k-1) - theta_extrap(i,j-1,k-1) ); + myhalf * dzi * ( thm_extrap(i,j,k+1) + thm_extrap(i,j-1,k+1) + - thm_extrap(i,j,k ) - thm_extrap(i,j-1,k ) ) : + fourth * dzi * ( thm_extrap(i,j,k+1) + thm_extrap(i,j-1,k+1) + - thm_extrap(i,j,k-1) - thm_extrap(i,j-1,k-1) ); Real gpy = gp_eta - (met_h_eta / met_h_zeta) * gp_zeta_on_jface; gpy *= mf_vy(i,j,0); @@ -369,11 +349,11 @@ void erf_substep_T (int step, int /*nrk*/, const Array4 & stage_zmom = S_stage_data[IntVars::zmom].const_array(mfi); const Array4 & prim = S_stage_prim.const_array(mfi); - const Array4& old_drho_u = Delta_rho_u.array(mfi); - const Array4& old_drho_v = Delta_rho_v.array(mfi); - const Array4& old_drho_w = Delta_rho_w.array(mfi); - const Array4& old_drho = Delta_rho.array(mfi); - const Array4& old_drho_theta = Delta_rho_theta.array(mfi); + const Array4& old_drho_u = Delta_rho_u.array(mfi); + const Array4& old_drho_v = Delta_rho_v.array(mfi); + const Array4& old_drho_w = Delta_rho_w.array(mfi); + const Array4& old_drho = Delta_rho.array(mfi); + const Array4& old_drho_thm = Delta_rho_theta.array(mfi); const Array4& slow_rhs_cons = S_slow_rhs[IntVars::cons].const_array(mfi); const Array4& slow_rhs_rho_w = S_slow_rhs[IntVars::zmom].const_array(mfi); @@ -384,7 +364,6 @@ void erf_substep_T (int step, int /*nrk*/, const Array4& cur_cons = S_data[IntVars::cons].array(mfi); const Array4& cur_zmom = S_data[IntVars::zmom].array(mfi); - // These store the advection momenta which we will use to update the slow variables const Array4& avg_zmom_arr = avg_zmom.array(mfi); const Array4& z_nd = z_phys_nd->const_array(mfi); @@ -392,7 +371,6 @@ void erf_substep_T (int step, int /*nrk*/, const Array4< Real>& omega_arr = Omega.array(mfi); - // Map factors const Array4& mf_mx = mapfac[MapFacType::m_x]->const_array(mfi); const Array4& mf_my = mapfac[MapFacType::m_y]->const_array(mfi); const Array4& mf_ux = mapfac[MapFacType::u_x]->const_array(mfi); @@ -400,19 +378,13 @@ void erf_substep_T (int step, int /*nrk*/, const Array4& mf_vx = mapfac[MapFacType::v_x]->const_array(mfi); const Array4& mf_vy = mapfac[MapFacType::v_y]->const_array(mfi); - // Create old_drho_u/v/w/theta = U'', V'', W'', Theta'' in the docs - // Note that we do the Copy and Subtract including one ghost cell - // so that we don't have to fill ghost cells of the new MultiFabs - // Initialize New_rho_u/v/w to Delta_rho_u/v/w so that - // the ghost cells in New_rho_u/v/w will match old_drho_u/v/w - FArrayBox temp_rhs_fab; FArrayBox RHS_fab; FArrayBox soln_fab; RHS_fab.resize (tbz,1,The_Async_Arena()); soln_fab.resize (tbz,1,The_Async_Arena()); - temp_rhs_fab.resize(tbz,2,The_Async_Arena()); + temp_rhs_fab.resize(tbz,3,The_Async_Arena()); auto const& RHS_a = RHS_fab.array(); auto const& soln_a = soln_fab.array(); @@ -434,7 +406,37 @@ void erf_substep_T (int step, int /*nrk*/, const GpuArray, AMREX_SPACEDIM> flx_arr{{AMREX_D_DECL(flux[0].array(), flux[1].array(), flux[2].array())}}; - // ********************************************************************* + // ************************************************************************* + // Define old contravariant velocity (start of sub-step) + // ************************************************************************* + { + Box gbxo = mfi.nodaltilebox(2); + Box gbxo_mid = gbxo; + + if (gbxo.smallEnd(2) == domlo.z) { + Box gbxo_lo = gbxo; gbxo_lo.setBig(2,gbxo.smallEnd(2)); + gbxo_mid.setSmall(2,gbxo.smallEnd(2)+1); + ParallelFor(gbxo_lo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + omega_arr(i,j,k) = zero; + }); + } + if (gbxo.bigEnd(2) == domhi.z+1) { + Box gbxo_hi = gbxo; gbxo_hi.setSmall(2,gbxo.bigEnd(2)); + gbxo_mid.setBig(2,gbxo.bigEnd(2)-1); + ParallelFor(gbxo_hi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + omega_arr(i,j,k) = old_drho_w(i,j,k); + }); + } + ParallelFor(gbxo_mid, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + omega_arr(i,j,k) = OmegaFromW(i,j,k,old_drho_w(i,j,k), + old_drho_u,old_drho_v, + mf_ux,mf_vy,z_nd,dxInv); + }); + } + + // ************************************************************************* + // Define lateral update to rho & rho theta & qv^{*} + // ************************************************************************* { BL_PROFILE("fast_T_making_rho_rhs"); ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { @@ -462,12 +464,21 @@ void erf_substep_T (int step, int /*nrk*/, Real mfsq = mf_mx(i,j,0) * mf_my(i,j,0); // NOTE: we are saving the (1/J) weighting for later when we add this to rho and theta - temp_rhs_arr(i,j,k,0) = ( xflux_hi - xflux_lo ) * dxi * mfsq + - ( yflux_hi - yflux_lo ) * dyi * mfsq; - temp_rhs_arr(i,j,k,1) = (( xflux_hi * (prim(i,j,k,0) + prim(i+1,j,k,0)) - - xflux_lo * (prim(i,j,k,0) + prim(i-1,j,k,0)) ) * dxi * mfsq+ - ( yflux_hi * (prim(i,j,k,0) + prim(i,j+1,k,0)) - - yflux_lo * (prim(i,j,k,0) + prim(i,j-1,k,0)) ) * dyi * mfsq) * myhalf; + temp_rhs_arr(i,j,k,Rho_comp) = ( xflux_hi - xflux_lo ) * dxi * mfsq + + ( yflux_hi - yflux_lo ) * dyi * mfsq; + temp_rhs_arr(i,j,k,RhoTheta_comp) = (( xflux_hi * (prim(i,j,k) + prim(i+1,j,k)) - + xflux_lo * (prim(i,j,k) + prim(i-1,j,k)) ) * dxi * mfsq+ + ( yflux_hi * (prim(i,j,k) + prim(i,j+1,k)) - + yflux_lo * (prim(i,j,k) + prim(i,j-1,k)) ) * dyi * mfsq) * myhalf; + if (l_use_moisture) { + temp_rhs_arr(i,j,k,RhoTheta_comp+1) = (( xflux_hi * (prim(i,j,k,PrimQ1_comp) + prim(i+1,j,k,PrimQ1_comp)) - + xflux_lo * (prim(i,j,k,PrimQ1_comp) + prim(i-1,j,k,PrimQ1_comp)) ) * dxi * mfsq+ + ( yflux_hi * (prim(i,j,k,PrimQ1_comp) + prim(i,j+1,k,PrimQ1_comp)) - + yflux_lo * (prim(i,j,k,PrimQ1_comp) + prim(i,j-1,k,PrimQ1_comp)) ) * dyi * mfsq) * myhalf; + } else { + temp_rhs_arr(i,j,k,RhoTheta_comp+1) = zero; + } + if (l_reflux) { (flx_arr[0])(i,j,k,0) = xflux_lo; @@ -486,34 +497,7 @@ void erf_substep_T (int step, int /*nrk*/, } } }); - } // end profile - - // ********************************************************************* - { - Box gbxo = mfi.nodaltilebox(2); - Box gbxo_mid = gbxo; - - if (gbxo.smallEnd(2) == domlo.z) { - Box gbxo_lo = gbxo; gbxo_lo.setBig(2,gbxo.smallEnd(2)); - gbxo_mid.setSmall(2,gbxo.smallEnd(2)+1); - ParallelFor(gbxo_lo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - omega_arr(i,j,k) = zero; - }); - } - if (gbxo.bigEnd(2) == domhi.z+1) { - Box gbxo_hi = gbxo; gbxo_hi.setSmall(2,gbxo.bigEnd(2)); - gbxo_mid.setBig(2,gbxo.bigEnd(2)-1); - ParallelFor(gbxo_hi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - omega_arr(i,j,k) = old_drho_w(i,j,k); - }); } - ParallelFor(gbxo_mid, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - omega_arr(i,j,k) = OmegaFromW(i,j,k,old_drho_w(i,j,k), - old_drho_u,old_drho_v, - mf_ux,mf_vy,z_nd,dxInv); - }); - } // end profile - // ********************************************************************* Box bx_shrunk_in_k = bx; int klo = tbz.smallEnd(2); @@ -526,6 +510,9 @@ void erf_substep_T (int step, int /*nrk*/, // We define halfg to match the notes (which is why we take the absolute value) Real halfg = std::abs(myhalf * grav_gpu[2]); + // ************************************************************************* + // RHS for tridiagonal system + // ************************************************************************* { BL_PROFILE("fast_loop_on_shrunk_t"); //Note we don't act on the bottom or top boundaries of the domain @@ -534,19 +521,60 @@ void erf_substep_T (int step, int /*nrk*/, Real coeff_P = coeffP_a(i,j,k); Real coeff_Q = coeffQ_a(i,j,k); - Real theta_t_lo = myhalf * ( prim(i,j,k-2,PrimTheta_comp) + prim(i,j,k-1,PrimTheta_comp) ); - Real theta_t_mid = myhalf * ( prim(i,j,k-1,PrimTheta_comp) + prim(i,j,k ,PrimTheta_comp) ); - Real theta_t_hi = myhalf * ( prim(i,j,k ,PrimTheta_comp) + prim(i,j,k+1,PrimTheta_comp) ); + // Cell center vars for constructing theta_m + Real thd_km2 = prim(i,j,k-2,PrimTheta_comp); + Real thd_km1 = prim(i,j,k-1,PrimTheta_comp); + Real thd_k = prim(i,j,k ,PrimTheta_comp); + Real thd_kp1 = prim(i,j,k+1,PrimTheta_comp); + + Real qv_km2 = (l_use_moisture) ? prim(i,j,k-2,PrimQ1_comp) : zero; + Real qv_km1 = (l_use_moisture) ? prim(i,j,k-1,PrimQ1_comp) : zero; + Real qv_k = (l_use_moisture) ? prim(i,j,k ,PrimQ1_comp) : zero; + Real qv_kp1 = (l_use_moisture) ? prim(i,j,k+1,PrimQ1_comp) : zero; + + Real thm_km2 = thd_km2 * (Real(1.0) + RvOverRd*qv_km2); + Real thm_km1 = thd_km1 * (Real(1.0) + RvOverRd*qv_km1); + Real thm_k = thd_k * (Real(1.0) + RvOverRd*qv_k ); + Real thm_kp1 = thd_kp1 * (Real(1.0) + RvOverRd*qv_kp1); + + // Stage theta_m at w-faces + Real thm_t_lo = myhalf * (thm_km2 + thm_km1); + Real thm_t_mid = myhalf * (thm_km1 + thm_k ); + Real thm_t_hi = myhalf * (thm_k + thm_kp1); + + // Stage coefficients for linearization + Real A_T_k = (one + RvOverRd*qv_k ); + Real A_T_km1 = (one + RvOverRd*qv_km1); + Real A_Q_k = RvOverRd*prim(i,j,k ,PrimTheta_comp); + Real A_Q_km1 = RvOverRd*prim(i,j,k-1,PrimTheta_comp); + Real A_D_k = -RvOverRd*prim(i,j,k ,PrimTheta_comp)*qv_k; + Real A_D_km1 = -RvOverRd*prim(i,j,k-1,PrimTheta_comp)*qv_km1; + + // Sums of slow RHS + Real Q_srhs_k = (l_use_moisture) ? slow_rhs_cons(i,j,k ,RhoQ1_comp) : zero; + Real Q_srhs_km1 = (l_use_moisture) ? slow_rhs_cons(i,j,k-1,RhoQ1_comp) : zero; + Real A_sum_srhs_k = A_D_k * slow_rhs_cons(i,j,k ,Rho_comp) + + A_T_k * slow_rhs_cons(i,j,k ,RhoTheta_comp) + + A_Q_k * Q_srhs_k; + Real A_sum_srhs_km1 = A_D_km1 * slow_rhs_cons(i,j,k-1,Rho_comp) + + A_T_km1 * slow_rhs_cons(i,j,k-1,RhoTheta_comp) + + A_Q_km1 * Q_srhs_km1; + + // Sums of the temp RHS + Real A_sum_trhs_k = A_D_k * temp_rhs_arr(i,j,k ,Rho_comp) + + A_T_k * temp_rhs_arr(i,j,k ,RhoTheta_comp) + + A_Q_k * temp_rhs_arr(i,j,k ,RhoTheta_comp+1); + Real A_sum_trhs_km1 = A_D_km1 * temp_rhs_arr(i,j,k-1,Rho_comp) + + A_T_km1 * temp_rhs_arr(i,j,k-1,RhoTheta_comp) + + A_Q_km1 * temp_rhs_arr(i,j,k-1,RhoTheta_comp+1); // line 2 last two terms (order dtau) - Real R0_tmp = -halfg * old_drho(i,j,k ) + coeff_P * old_drho_theta(i,j,k ) - -halfg * old_drho(i,j,k-1) + coeff_Q * old_drho_theta(i,j,k-1); + Real R0_tmp = -halfg * old_drho(i,j,k ) + coeff_P * old_drho_thm(i,j,k ) + -halfg * old_drho(i,j,k-1) + coeff_Q * old_drho_thm(i,j,k-1); // line 3 residuals (order dtau^2) one <-> beta_2 - Real R1_tmp = -halfg * ( slow_rhs_cons(i,j,k ,Rho_comp) + slow_rhs_cons(i,j,k-1,Rho_comp) ); - - R1_tmp += coeff_P * slow_rhs_cons(i,j,k ,RhoTheta_comp) - + coeff_Q * slow_rhs_cons(i,j,k-1,RhoTheta_comp); + Real R1_tmp = -halfg * ( slow_rhs_cons(i,j,k ,Rho_comp) + slow_rhs_cons(i,j,k-1,Rho_comp) ); + R1_tmp += coeff_P * A_sum_srhs_k + coeff_Q * A_sum_srhs_km1; Real Omega_kp1 = omega_arr(i,j,k+1); Real Omega_k = omega_arr(i,j,k ); @@ -559,8 +587,8 @@ void erf_substep_T (int step, int /*nrk*/, + temp_rhs_arr(i,j,k,Rho_comp)/detJ(i,j,k) + temp_rhs_arr(i,j,k-1,Rho_comp)/detJ(i,j,k-1) ); // consolidate lines 6&7 (order dtau^2) - R1_tmp += -( coeff_P/detJ(i,j,k ) * ( beta_1 * dzi * (Omega_kp1*theta_t_hi - Omega_k*theta_t_mid) + temp_rhs_arr(i,j,k ,RhoTheta_comp) ) - + coeff_Q/detJ(i,j,k-1) * ( beta_1 * dzi * (Omega_k*theta_t_mid - Omega_km1*theta_t_lo) + temp_rhs_arr(i,j,k-1,RhoTheta_comp) ) ); + R1_tmp += -( coeff_P/detJ(i,j,k ) * ( beta_1 * dzi * (Omega_kp1*thm_t_hi - Omega_k*thm_t_mid ) + A_sum_trhs_k ) + + coeff_Q/detJ(i,j,k-1) * ( beta_1 * dzi * (Omega_k*thm_t_mid - Omega_km1*thm_t_lo) + A_sum_trhs_km1) ); // line 1 RHS_a(i,j,k) = old_drho_w(i,j,k) + dtau * (slow_rhs_rho_w(i,j,k) + zmom_src_arr(i,j,k) + R0_tmp + dtau*beta_2*R1_tmp); @@ -572,19 +600,22 @@ void erf_substep_T (int step, int /*nrk*/, }); } // end profile - Box b2d = tbz; // Copy constructor + Box b2d = tbz; b2d.setRange(2,0); auto const lo = lbound(bx); auto const hi = ubound(bx); + // ************************************************************************* + // Solve tridiagonal system + // ************************************************************************* { BL_PROFILE("substep_b2d_loop_t"); #ifdef AMREX_USE_GPU ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int) { // w_klo, w_khi given by specified Dirichlet values - RHS_a(i,j,lo.z ) = dtau * (slow_rhs_rho_w(i,j,lo.z) + zmom_src_arr(i,j,lo.z)); + RHS_a(i,j,lo.z ) = dtau * (slow_rhs_rho_w(i,j,lo.z ) + zmom_src_arr(i,j,lo.z )); RHS_a(i,j,hi.z+1) = dtau * (slow_rhs_rho_w(i,j,hi.z+1) + zmom_src_arr(i,j,hi.z+1)); // w = specified Dirichlet value at k = lo.z @@ -598,7 +629,7 @@ void erf_substep_T (int step, int /*nrk*/, cur_zmom(i,j,hi.z+1) = stage_zmom(i,j,hi.z+1) + soln_a(i,j,hi.z+1); for (int k = hi.z; k >= lo.z; k--) { - soln_a(i,j,k) -= ( coeffC_a(i,j,k) * inv_coeffB_a(i,j,k) ) *soln_a(i,j,k+1); + soln_a(i,j,k) -= ( coeffC_a(i,j,k) * inv_coeffB_a(i,j,k) ) * soln_a(i,j,k+1); } }); #else @@ -671,13 +702,13 @@ void erf_substep_T (int step, int /*nrk*/, }); // ************************************************************************** - // Define updates in the RHS of rho and (rho theta) + // Define updates in the RHS of rho, rho theta, rho qv // ************************************************************************** { BL_PROFILE("fast_rho_final_update"); ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - Real zflux_lo = beta_2 * soln_a(i,j,k ) + beta_1 * omega_arr(i,j,k); + Real zflux_lo = beta_2 * soln_a(i,j,k ) + beta_1 * omega_arr(i,j,k ); Real zflux_hi = beta_2 * soln_a(i,j,k+1) + beta_1 * omega_arr(i,j,k+1); // Note that in the solve we effectively impose new_drho_w(i,j,vbx_hi.z+1)=0 @@ -695,16 +726,21 @@ void erf_substep_T (int step, int /*nrk*/, } } - Real fast_rhs_rho = -(temp_rhs_arr(i,j,k,0) + ( zflux_hi - zflux_lo ) * dzi) / detJ(i,j,k); - + Real fast_rhs_rho = -(temp_rhs_arr(i,j,k,Rho_comp) + ( zflux_hi - zflux_lo ) * dzi) / detJ(i,j,k); cur_cons(i,j,k,0) += dtau * (slow_rhs_cons(i,j,k,0) + fast_rhs_rho); - Real fast_rhs_rhotheta = -( temp_rhs_arr(i,j,k,1) + myhalf * - ( zflux_hi * (prim(i,j,k) + prim(i,j,k+1)) - - zflux_lo * (prim(i,j,k) + prim(i,j,k-1)) ) * dzi ) / detJ(i,j,k); - + Real fast_rhs_rhotheta = -( temp_rhs_arr(i,j,k,RhoTheta_comp) + myhalf * + ( zflux_hi * (prim(i,j,k) + prim(i,j,k+1)) - + zflux_lo * (prim(i,j,k) + prim(i,j,k-1)) ) * dzi ) / detJ(i,j,k); cur_cons(i,j,k,1) += dtau * (slow_rhs_cons(i,j,k,1) + fast_rhs_rhotheta); + if (l_use_moisture) { + Real fast_rhs_rhoqv = -( temp_rhs_arr(i,j,k,RhoTheta_comp+1) + myhalf * + ( zflux_hi * (prim(i,j,k,PrimQ1_comp) + prim(i,j,k+1,PrimQ1_comp)) - + zflux_lo * (prim(i,j,k,PrimQ1_comp) + prim(i,j,k-1,PrimQ1_comp)) ) * dzi ) / detJ(i,j,k); + cur_cons(i,j,k,RhoQ1_comp) += dtau * (slow_rhs_cons(i,j,k,RhoQ1_comp) + fast_rhs_rhoqv); + } + if (l_reflux) { (flx_arr[2])(i,j,k,1) = (flx_arr[2])(i,j,k,0) * myhalf * (prim(i,j,k) + prim(i,j,k-1)); } @@ -712,6 +748,9 @@ void erf_substep_T (int step, int /*nrk*/, // add in source terms for cell-centered conserved variables cur_cons(i,j,k,Rho_comp) += dtau * cc_src_arr(i,j,k,Rho_comp); cur_cons(i,j,k,RhoTheta_comp) += dtau * cc_src_arr(i,j,k,RhoTheta_comp); + if (l_use_moisture) { + cur_cons(i,j,k,RhoQ1_comp) = amrex::max(zero, cur_cons(i,j,k,RhoQ1_comp) + dtau * cc_src_arr(i,j,k,RhoQ1_comp)); + } }); } // end profile diff --git a/Source/TimeIntegration/ERF_TI_slow_rhs_pre.H b/Source/TimeIntegration/ERF_TI_slow_rhs_pre.H index a90994c38..e98465dca 100644 --- a/Source/TimeIntegration/ERF_TI_slow_rhs_pre.H +++ b/Source/TimeIntegration/ERF_TI_slow_rhs_pre.H @@ -105,9 +105,9 @@ qheating_rates[level].get(), terrain_blank, fine_geom, solverChoice, mapfac[level], - rhotheta_src[level].get(), rhoqt_src[level].get(), + macrophysics_src[level].get(), rhotheta_src[level].get(), rhoqt_src[level].get(), dptr_wbar_sub, d_rayleigh_ptrs_at_lev, d_sinesq_at_lev, - (solverChoice.hindcast_surface_bcs? &surface_state_interp[level] : nullptr), + (solverChoice.hindcast_surface_bcs? &surface_state_interp[level] : nullptr), input_sounding_data, turbPert, true); // ***************************************************************************** @@ -169,10 +169,13 @@ fr_as_crse, fr_as_fine); if ((solverChoice.vert_implicit_fac[nrk] > zero) && solverChoice.implicit_before_substep) { - MultiFab scratch(S_data[IntVars::cons].boxArray(),S_data[IntVars::cons].DistributionMap(), 2, + MultiFab scratch(S_data[IntVars::cons].boxArray(),S_data[IntVars::cons].DistributionMap(), 3, S_data[IntVars::cons].nGrowVect()); - MultiFab::Copy(scratch, S_old[IntVars::cons], 0, 0, 2, S_data[IntVars::cons].nGrowVect()); // scratch := S_old (for rho, rhotheta) - MultiFab::Saxpy(scratch, slow_dt, S_rhs[IntVars::cons], 0, 0, 2, 0); // scratch := S_old + slow_dt*Src (for rho, rhotheta) + MultiFab::Copy(scratch, S_old[IntVars::cons], 0, 0, 2, S_data[IntVars::cons].nGrowVect()); // scratch := S_old (for rho, rhotheta) + MultiFab::Copy(scratch, S_old[IntVars::cons], RhoQ1_comp, 2, 1, S_data[IntVars::cons].nGrowVect()); // scratch := S_old (for rhoQ1) + MultiFab::Saxpy(scratch, slow_dt, S_rhs[IntVars::cons], 0, 0, 2, 0); // scratch := S_old + slow_dt*Src (for rho, rhotheta) + MultiFab::Saxpy(scratch, slow_dt, S_rhs[IntVars::cons], RhoQ1_comp, 2, 1, 0); // scratch := S_old + slow_dt*Src (for rhoQ1) + scratch.FillBoundary(geom[level].periodicity()); MultiFab scratch_xmom(S_data[IntVars::xmom].boxArray(), diff --git a/Source/TimeIntegration/ERF_TI_substep_fun.H b/Source/TimeIntegration/ERF_TI_substep_fun.H index 4496dc678..cb449bd15 100644 --- a/Source/TimeIntegration/ERF_TI_substep_fun.H +++ b/Source/TimeIntegration/ERF_TI_substep_fun.H @@ -64,7 +64,7 @@ auto acoustic_substepping_fun = [&](int fast_step, int n_sub, int nrk, qheating_rates[level].get(), terrain_blank, fine_geom, solverChoice, mapfac[level], - rhotheta_src[level].get(), rhoqt_src[level].get(), + macrophysics_src[level].get(), rhotheta_src[level].get(), rhoqt_src[level].get(), dptr_wbar_sub, d_rayleigh_ptrs_at_lev, d_sinesq_at_lev, (solverChoice.hindcast_surface_bcs? &surface_state_interp[level] : nullptr), diff --git a/Source/Utils/ERF_MacroPhysics.H b/Source/Utils/ERF_MacroPhysics.H new file mode 100644 index 000000000..3b197a97a --- /dev/null +++ b/Source/Utils/ERF_MacroPhysics.H @@ -0,0 +1,136 @@ +#ifndef ERF_MACROPHYS_H +#define ERF_MACROPHYS_H + +#include "ERF_EOS.H" +#include "ERF_Constants.H" +#include "ERF_IndexDefines.H" +#include "ERF_MicrophysicsUtils.H" + +AMREX_GPU_HOST_DEVICE +AMREX_FORCE_INLINE +static void +NewtonIterSat (const amrex::Real& fac_cond, + amrex::Real& tabs, + const amrex::Real& pres, + amrex::Real& qv) +{ + // Solution tolerance + amrex::Real tol = amrex::Real(1.0e-6); + + // Saturation moisture fractions + amrex::Real qsat, dqsat; + + // Newton iteration vars + int niter; + amrex::Real fff, dfff; + amrex::Real dtabs; + + // Initial guess for temperature + amrex::Real tabs_old = tabs; + + niter = 0; + dtabs = 1; + + //================================================== + // Newton iteration to qv=qsat (cloud phase only) + //================================================== + do { + // Saturation moisture fractions + erf_qsatw(tabs, pres, qsat); + erf_dtqsatw(tabs, pres, dqsat); + + // Function for root finding: + // 0 = -T_new + T_old + L_v/C_p * (qv - qsat) + fff = -tabs + tabs_old + fac_cond*(qv - qsat); + + // Derivative of function (T_new iterated on) + dfff = -one - fac_cond*dqsat; + + // Update the temperature + dtabs = -fff/dfff; + tabs += dtabs; + + // Update iteration + niter = niter+1; + } while(std::abs(dtabs) > tol && niter < 20); + + // Update from last iteration (dq = dq/dt * dt) + qv = qsat + dqsat*dtabs; +} + +void +GetMacroPhysicsSrc (int& real_width, + amrex::Real& dt, + amrex::Geometry& geom, + amrex::MultiFab* macro_src, + amrex::MultiFab& cons_data) +{ + amrex::Real idt = amrex::Real(1.) / dt; + amrex::Real d_rdOcp = R_d / Cp_d; + amrex::Real d_fac_cond = L_v / Cp_d; + +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (amrex::MFIter mfi(cons_data,TileNoZ()); mfi.isValid(); ++mfi) + { + amrex::Box tbx = mfi.tilebox(); + + const amrex::Array4& mac_src = macro_src->array(mfi); + const amrex::Array4& cell_data = cons_data.array(mfi); + + // No macrophysics in the nudging region + const amrex::Box& domain = geom.Domain(); + int i_lo = domain.smallEnd(0); + int i_hi = domain.bigEnd(0); + int j_lo = domain.smallEnd(1); + int j_hi = domain.bigEnd(1); + if (tbx.smallEnd(0) == i_lo) { tbx.growLo(0,-real_width); } + if (tbx.bigEnd(0) == i_hi) { tbx.growHi(0,-real_width); } + if (tbx.smallEnd(1) == j_lo) { tbx.growLo(1,-real_width); } + if (tbx.bigEnd(1) == j_hi) { tbx.growHi(1,-real_width); } + + // Theta tendency and update qv/qc in place + ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + // Conserved CC vars at start of step + amrex::Real rho_s = cell_data(i,j,k,Rho_comp); + amrex::Real rt_s = cell_data(i,j,k,RhoTheta_comp); + amrex::Real th_s = cell_data(i,j,k,RhoTheta_comp) / rho_s; + amrex::Real qv_s = cell_data(i,j,k,RhoQ1_comp ) / rho_s; + amrex::Real qc_s = cell_data(i,j,k,RhoQ2_comp ) / rho_s; + + // To be updated in place + amrex::Real qs; + amrex::Real qv_f = qv_s; + amrex::Real tabs_f = getTgivenRandRTh(rho_s, rt_s, qv_s); + + // Invariants + amrex::Real qt = qv_s + qc_s; + amrex::Real pres = getPgivenRTh(rt_s, qv_s); + + // Valid Newton iteration + erf_qsatw(tabs_f, pres*amrex::Real(0.01), qs); + if (qt > qs) { + NewtonIterSat(d_fac_cond, tabs_f, pres*amrex::Real(0.01), qv_f); + } + // Put qc into qv and double check Newton iter criteria + else { + qv_f += qc_s; + tabs_f -= d_fac_cond * qc_s; + erf_qsatw(tabs_f, pres*amrex::Real(0.01), qs); + if (qv_f > qs) { + NewtonIterSat(d_fac_cond, tabs_f, pres*amrex::Real(0.01), qv_f); + } + } + + // Populate tendency and data + amrex::Real th_f = getThgivenTandP(tabs_f, pres, d_rdOcp); + mac_src(i,j,k,0) = rho_s * idt * ( th_f - th_s ); + mac_src(i,j,k,1) = rho_s * idt * ( qv_f - qv_s ); + mac_src(i,j,k,2) = rho_s * idt * ( amrex::max(amrex::Real(0.0),qt - qv_f) - qc_s ); + }); + } // mfi +} + +#endif diff --git a/Source/Utils/Make.package b/Source/Utils/Make.package index 3960c5fe7..deebc3821 100644 --- a/Source/Utils/Make.package +++ b/Source/Utils/Make.package @@ -17,6 +17,7 @@ CEXE_headers += ERF_SatMethods.H CEXE_headers += ERF_DirectionSelector.H CEXE_headers += ERF_StormDiagnostics.H CEXE_headers += ERF_HurricaneDiagnostics.H +CEXE_headers += ERF_MacroPhysics.H CEXE_sources += ERF_AverageDown.cpp CEXE_sources += ERF_ChopGrids.cpp