More steps towards making moisture a run-time setting (#1335)

* More steps towards making moisture a run-time setting

* fix tabs

* fix cmake build

CMake/BuildERFExe.cmake

@@ -159,6 +159,10 @@ function(build_erf_lib erf_lib_name)
        ${SRC_DIR}/Microphysics/Kessler/Kessler.cpp
        ${SRC_DIR}/Microphysics/Kessler/Diagnose_Kessler.cpp
        ${SRC_DIR}/Microphysics/Kessler/Update_Kessler.cpp
+       ${SRC_DIR}/Microphysics/FastEddy/Init_FE.cpp
+       ${SRC_DIR}/Microphysics/FastEddy/FastEddy.cpp
+       ${SRC_DIR}/Microphysics/FastEddy/Diagnose_FE.cpp
+       ${SRC_DIR}/Microphysics/FastEddy/Update_FE.cpp
   )
 
   if(NOT "${erf_exe_name}" STREQUAL "erf_unit_tests")
@@ -203,6 +207,7 @@ function(build_erf_lib erf_lib_name)
   target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/Microphysics/Null)
   target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/Microphysics/SAM)
   target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/Microphysics/Kessler)
+  target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/Microphysics/FastEddy)
 
   if(ERF_ENABLE_RRTMGP)
      target_link_libraries(${erf_lib_name} PUBLIC yakl)

Docs/sphinx_doc/Inputs.rst

@@ -893,7 +893,8 @@ List of Parameters
 | Parameter                   | Definition               | Acceptable         | Default    |
 |                             |                          | Values             |            |
 +=============================+==========================+====================+============+
-| **erf.moisture_model**      | Name of moisture model   |  "SAM", "Kessler"  | "Null"     |
+| **erf.moisture_model**      | Name of moisture model   |  "SAM", "Kessler", | "Null"     |
+|                             |                          |  "FastEddy"        |            |
 +-----------------------------+--------------------------+--------------------+------------+
 | **erf.do_cloud**            | use basic moisture model |  true / false      | true       |
 +-----------------------------+--------------------------+--------------------+------------+

Exec/Make.ERF

@@ -112,6 +112,11 @@ include $(ERF_MOISTURE_NULL_DIR)/Make.package
 VPATH_LOCATIONS   += $(ERF_MOISTURE_NULL_DIR)
 INCLUDE_LOCATIONS += $(ERF_MOISTURE_NULL_DIR)
 
+ERF_MOISTURE_FE_DIR = $(ERF_SOURCE_DIR)/Microphysics/FastEddy
+include $(ERF_MOISTURE_FE_DIR)/Make.package
+VPATH_LOCATIONS   += $(ERF_MOISTURE_FE_DIR)
+INCLUDE_LOCATIONS += $(ERF_MOISTURE_FE_DIR)
+
 ERF_MOISTURE_SAM_DIR = $(ERF_SOURCE_DIR)/Microphysics/SAM
 include $(ERF_MOISTURE_SAM_DIR)/Make.package
 VPATH_LOCATIONS   += $(ERF_MOISTURE_SAM_DIR)

Source/DataStructs/DataStruct.H

@@ -29,7 +29,7 @@ enum struct TerrainType {
 };
 
 enum struct MoistureType {
-    Kessler, SAM, None
+    Kessler, SAM, FastEddy, None
 };
 
 enum struct Coord {
@@ -90,6 +90,8 @@ struct SolverChoice {
             moisture_type = MoistureType::SAM;
         } else if (moisture_model_string == "Kessler") {
             moisture_type = MoistureType::Kessler;
+        } else if (moisture_model_string == "FastEddy") {
+            moisture_type = MoistureType::FastEddy;
         } else {
             moisture_type = MoistureType::None;
         }

Source/ERF.cpp

@@ -905,8 +905,6 @@ ERF::init_only (int lev, Real time)
     lev_new[Vars::xvel].OverrideSync(geom[lev].periodicity());
     lev_new[Vars::yvel].OverrideSync(geom[lev].periodicity());
     lev_new[Vars::zvel].OverrideSync(geom[lev].periodicity());
-
-    print_state(vars_new[0][Vars::cons], IntVect(30,30,2));
 }
 
 // read in some parameters from inputs file
@@ -965,6 +963,8 @@ ERF::ReadParameters ()
             micro.SetModel<SAM>();
         } else if (solverChoice.moisture_type == MoistureType::Kessler) {
             micro.SetModel<Kessler>();
+        } else if (solverChoice.moisture_type == MoistureType::FastEddy) {
+            micro.SetModel<FastEddy>();
         } else if (solverChoice.moisture_type == MoistureType::None) {
             micro.SetModel<NullMoist>();
             amrex::Print() << "WARNING: Compiled with moisture but using NullMoist model!\n";

Source/IO/Checkpoint.cpp

@@ -106,34 +106,34 @@ ERF::WriteCheckpointFile () const
        }
    }
 
-   // write the MultiFab data to, e.g., chk00010/Level_0/
-   // Here we make copies of the MultiFab with no ghost cells
-   for (int lev = 0; lev <= finest_level; ++lev)
-   {
-       MultiFab cons(grids[lev],dmap[lev],Cons::NumVars,0);
-       MultiFab::Copy(cons,vars_new[lev][Vars::cons],0,0,NVAR,0);
-       VisMF::Write(cons, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "Cell"));
-
-       MultiFab xvel(convert(grids[lev],IntVect(1,0,0)),dmap[lev],1,0);
-       MultiFab::Copy(xvel,vars_new[lev][Vars::xvel],0,0,1,0);
-       VisMF::Write(xvel, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "XFace"));
-
-       MultiFab yvel(convert(grids[lev],IntVect(0,1,0)),dmap[lev],1,0);
-       MultiFab::Copy(yvel,vars_new[lev][Vars::yvel],0,0,1,0);
-       VisMF::Write(yvel, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "YFace"));
-
-       MultiFab zvel(convert(grids[lev],IntVect(0,0,1)),dmap[lev],1,0);
-       MultiFab::Copy(zvel,vars_new[lev][Vars::zvel],0,0,1,0);
-       VisMF::Write(zvel, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "ZFace"));
-
-       IntVect ng;
-#ifdef ERF_USE_MOISTURE
-       // We must read and write qmoist with ghost cells because we don't directly impose BCs on these vars
-       ng = qmoist[lev].nGrowVect();
-       MultiFab moist_vars(grids[lev],dmap[lev],qmoist[lev].nComp(),ng);
-       MultiFab::Copy(moist_vars,qmoist[lev],0,0,qmoist[lev].nComp(),ng);
-       VisMF::Write(moist_vars, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "MoistVars"));
-#endif
+    // write the MultiFab data to, e.g., chk00010/Level_0/
+    // Here we make copies of the MultiFab with no ghost cells
+    for (int lev = 0; lev <= finest_level; ++lev)
+    {
+        MultiFab cons(grids[lev],dmap[lev],Cons::NumVars,0);
+        MultiFab::Copy(cons,vars_new[lev][Vars::cons],0,0,NVAR,0);
+        VisMF::Write(cons, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "Cell"));
+
+        MultiFab xvel(convert(grids[lev],IntVect(1,0,0)),dmap[lev],1,0);
+        MultiFab::Copy(xvel,vars_new[lev][Vars::xvel],0,0,1,0);
+        VisMF::Write(xvel, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "XFace"));
+
+        MultiFab yvel(convert(grids[lev],IntVect(0,1,0)),dmap[lev],1,0);
+        MultiFab::Copy(yvel,vars_new[lev][Vars::yvel],0,0,1,0);
+        VisMF::Write(yvel, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "YFace"));
+
+        MultiFab zvel(convert(grids[lev],IntVect(0,0,1)),dmap[lev],1,0);
+        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) {
+            // We must read and write qmoist with ghost cells because we don't directly impose BCs on these vars
+            ng = qmoist[lev].nGrowVect();
+            MultiFab moist_vars(grids[lev],dmap[lev],qmoist[lev].nComp(),ng);
+            MultiFab::Copy(moist_vars,qmoist[lev],0,0,qmoist[lev].nComp(),ng);
+            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();
@@ -249,7 +249,14 @@ ERF::ReadCheckpointFile ()
     // conservative, cell-centered vars
     is >> chk_ncomp;
     GotoNextLine(is);
-    AMREX_ASSERT(chk_ncomp == Cons::NumVars);
+#if 0
+    if (solverChoice.moisture_type != MoistureType::None) {
+        AMREX_ASSERT(chk_ncomp == Cons::NumVars);
+    } else {
+        // This assumes that we are carrying RhoQ1 and RhoQ2 but not actually using them
+        AMREX_ASSERT(chk_ncomp == Cons::NumVars-2);
+    }
+#endif
 
     // x-velocity on faces
     is >> chk_ncomp;
@@ -314,7 +321,16 @@ ERF::ReadCheckpointFile ()
     {
         MultiFab cons(grids[lev],dmap[lev],Cons::NumVars,0);
         VisMF::Read(cons, amrex::MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "Cell"));
-        MultiFab::Copy(vars_new[lev][Vars::cons],cons,0,0,Cons::NumVars,0);
+
+        if (solverChoice.moisture_type != MoistureType::None) {
+            MultiFab::Copy(vars_new[lev][Vars::cons],cons,0,0,Cons::NumVars,0);
+        } else {
+            // This assumes that we are carrying RhoQ1 and RhoQ2 but not actually using them
+            MultiFab::Copy(vars_new[lev][Vars::cons],cons,0,0,Cons::NumVars-2,0);
+
+            // We zero them here so they won't cause problems while we're still debugging
+            vars_new[lev][Vars::cons].setVal(0.0, RhoQ1_comp, 2);
+        }
 
         MultiFab xvel(convert(grids[lev],IntVect(1,0,0)),dmap[lev],1,0);
         VisMF::Read(xvel, amrex::MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "XFace"));
@@ -329,12 +345,14 @@ ERF::ReadCheckpointFile ()
         MultiFab::Copy(vars_new[lev][Vars::zvel],zvel,0,0,1,0);
 
         IntVect ng;
-#ifdef ERF_USE_MOISTURE
-        ng = qmoist[lev].nGrowVect();
-        MultiFab moist_vars(grids[lev],dmap[lev],qmoist[lev].nComp(),ng);
-        VisMF::Read(moist_vars, amrex::MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "MoistVars"));
-        MultiFab::Copy(qmoist[lev],moist_vars,0,0,qmoist[lev].nComp(),ng);
-#endif
+        if (solverChoice.moisture_type == MoistureType::None) {
+            qmoist[lev].setVal(0.0);
+        } else {
+            ng = qmoist[lev].nGrowVect();
+            MultiFab moist_vars(grids[lev],dmap[lev],qmoist[lev].nComp(),ng);
+            VisMF::Read(moist_vars, amrex::MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "MoistVars"));
+            MultiFab::Copy(qmoist[lev],moist_vars,0,0,qmoist[lev].nComp(),ng);
+        }
 
         // Note that we read the ghost cells of the base state (unlike above)
         ng = base_state[lev].nGrowVect();

Source/IO/ERF_Write1DProfiles.cpp

@@ -194,9 +194,7 @@ ERF::derive_diag_profiles(Gpu::HostVector<Real>& h_avg_u   , Gpu::HostVector<Rea
 
     MultiFab p_hse (base_state[lev], make_alias, 1, 1); // p_0  is second component
 
-#if defined(ERF_USE_MOISTURE)
-    MultiFab qv(qmoist[lev], make_alias, 0, 1);
-#endif
+    bool use_moisture = (solverChoice.moisture_type != MoistureType::None);
 
     for ( MFIter mfi(mf_cons,TilingIfNotGPU()); mfi.isValid(); ++mfi)
     {
@@ -206,9 +204,6 @@ ERF::derive_diag_profiles(Gpu::HostVector<Real>& h_avg_u   , Gpu::HostVector<Rea
         const Array4<Real>& v_cc_arr =  v_cc.array(mfi);
         const Array4<Real>& w_cc_arr =  w_cc.array(mfi);
         const Array4<Real>& cons_arr = mf_cons.array(mfi);
-#if defined(ERF_USE_MOISTURE)
-        const Array4<Real>&   qv_arr = qv.array(mfi);
-#endif
         const Array4<Real>&   p0_arr = p_hse.array(mfi);
 
         ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
@@ -237,18 +232,45 @@ ERF::derive_diag_profiles(Gpu::HostVector<Real>& h_avg_u   , Gpu::HostVector<Rea
             fab_arr(i, j, k,14) = tke * u_cc_arr(i,j,k);   // ku
             fab_arr(i, j, k,15) = tke * v_cc_arr(i,j,k);   // kv
             fab_arr(i, j, k,16) = tke * w_cc_arr(i,j,k);   // kw
-#if defined(ERF_USE_MOISTURE)
-            Real p = getPgivenRTh(cons_arr(i, j, k, RhoTheta_comp), qv_arr(i,j,k));
-#else
-            Real p = getPgivenRTh(cons_arr(i, j, k, RhoTheta_comp));
-#endif
-            p -= p0_arr(i,j,k);
-            fab_arr(i, j, k,17) = p;                       // p'
-            fab_arr(i, j, k,18) = p * u_cc_arr(i,j,k);     // p'u
-            fab_arr(i, j, k,19) = p * v_cc_arr(i,j,k);     // p'v
-            fab_arr(i, j, k,20) = p * w_cc_arr(i,j,k);     // p'w
+
+            if (!use_moisture) {
+                Real p = getPgivenRTh(cons_arr(i, j, k, RhoTheta_comp));
+                p -= p0_arr(i,j,k);
+                fab_arr(i, j, k,17) = p;                       // p'
+                fab_arr(i, j, k,18) = p * u_cc_arr(i,j,k);     // p'u
+                fab_arr(i, j, k,19) = p * v_cc_arr(i,j,k);     // p'v
+                fab_arr(i, j, k,20) = p * w_cc_arr(i,j,k);     // p'w
+            }
         });
-    }
+    } // mfi
+
+    if (use_moisture)
+    {
+        MultiFab qv(qmoist[lev], make_alias, 0, 1);
+
+        for ( MFIter mfi(mf_cons,TilingIfNotGPU()); mfi.isValid(); ++mfi)
+        {
+            const Box& bx = mfi.tilebox();
+            const Array4<Real>& fab_arr = mf_out.array(mfi);
+            const Array4<Real>& cons_arr = mf_cons.array(mfi);
+            const Array4<Real>& u_cc_arr =  u_cc.array(mfi);
+            const Array4<Real>& v_cc_arr =  v_cc.array(mfi);
+            const Array4<Real>& w_cc_arr =  w_cc.array(mfi);
+            const Array4<Real>&   p0_arr = p_hse.array(mfi);
+            const Array4<Real>&   qv_arr = qv.array(mfi);
+
+            ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
+            {
+                Real p = getPgivenRTh(cons_arr(i, j, k, RhoTheta_comp), qv_arr(i,j,k));
+
+                p -= p0_arr(i,j,k);
+                fab_arr(i, j, k,17) = p;                       // p'
+                fab_arr(i, j, k,18) = p * u_cc_arr(i,j,k);     // p'u
+                fab_arr(i, j, k,19) = p * v_cc_arr(i,j,k);     // p'v
+                fab_arr(i, j, k,20) = p * w_cc_arr(i,j,k);     // p'w
+            });
+        } // mfi
+    } // use_moisture
 
     h_avg_rho  = sumToLine(mf_out, 0,1,domain,zdir);
     h_avg_th   = sumToLine(mf_out, 1,1,domain,zdir);

Source/IO/Plotfile.cpp.convergence

@@ -544,7 +544,7 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
             mf_comp ++;
         }
 
-#if defined(ERF_USE_MOISTURE)
+        if (solverChoice.moisture_type != MoistureType::None) {
         calculate_derived("qt",          derived::erf_derQt);
         calculate_derived("qp",          derived::erf_derQp);
 
@@ -590,7 +590,8 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
             MultiFab::Copy(mf[lev],qg_mf,0,mf_comp,1,0);
             mf_comp += 1;
         }
-#elif defined(ERF_USE_WARM_NO_PRECIP)
+        }
+#if defined(ERF_USE_WARM_NO_PRECIP)
         calculate_derived("qv",          derived::erf_derQv);
         calculate_derived("qc",          derived::erf_derQc);
 #endif

Source/IO/ReadFromWRFInput.cpp

@@ -1,5 +1,6 @@
 #include "NCWpsFile.H"
 #include "AMReX_FArrayBox.H"
+#include "DataStruct.H"
 
 using namespace amrex;
 
@@ -16,14 +17,12 @@ read_from_wrfinput (int lev,
                     FArrayBox& NC_MSFM_fab, FArrayBox& NC_SST_fab,
                     FArrayBox& NC_C1H_fab , FArrayBox& NC_C2H_fab,
                     FArrayBox& NC_RDNW_fab,
-#if defined(ERF_USE_MOISTURE)
                     FArrayBox& NC_QVAPOR_fab,
                     FArrayBox& NC_QCLOUD_fab,
                     FArrayBox& NC_QRAIN_fab,
-#elif defined(ERF_USE_WARM_NO_PRECIP)
-#endif
                     FArrayBox& NC_PH_fab  , FArrayBox& NC_PHB_fab,
-                    FArrayBox& NC_ALB_fab , FArrayBox& NC_PB_fab)
+                    FArrayBox& NC_ALB_fab , FArrayBox& NC_PB_fab,
+                    MoistureType moisture_type)
 {
     amrex::Print() << "Loading initial data from NetCDF file at level " << lev << std::endl;
 
@@ -51,12 +50,12 @@ read_from_wrfinput (int lev,
     NC_fabs.push_back(&NC_C1H_fab);       NC_names.push_back("C1H");  NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT);         // 15
     NC_fabs.push_back(&NC_C2H_fab);       NC_names.push_back("C2H");  NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT);         // 16
     NC_fabs.push_back(&NC_RDNW_fab);      NC_names.push_back("RDNW"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT);         // 17
-#if defined(ERF_USE_MOISTURE)
-    NC_fabs.push_back(&NC_QVAPOR_fab);    NC_names.push_back("QVAPOR"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE);         // 18
-    NC_fabs.push_back(&NC_QCLOUD_fab);    NC_names.push_back("QCLOUD"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE);         // 19
-    NC_fabs.push_back(&NC_QRAIN_fab);     NC_names.push_back("QRAIN"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE);         // 20
-# elif defined(ERF_USE_WARM_NO_PRECIP)
-#endif
+
+    if (moisture_type != MoistureType::None) {
+        NC_fabs.push_back(&NC_QVAPOR_fab);    NC_names.push_back("QVAPOR"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE);         // 18
+        NC_fabs.push_back(&NC_QCLOUD_fab);    NC_names.push_back("QCLOUD"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE);         // 19
+        NC_fabs.push_back(&NC_QRAIN_fab);     NC_names.push_back("QRAIN"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE);         // 20
+    }
 
     // Read the netcdf file and fill these FABs
     amrex::Print() << "Building initial FABS from file " << fname << std::endl;