Merge branch 'development' of github.com:erf-model/ERF into explicit_…

…MOST

CHANGES

@@ -1,3 +1,22 @@
+# 24.03
+    -- AMReX submodule set to 24.03 release hash (d737886)
+
+    -- Plot rain accumulation (#1478)
+
+    -- Remove FastEddy completely (#1483)
+
+    -- Rain accumulation capability (#1474)
+
+    -- Remove the FastEddy model since it coincides with Kessler (#1473)
+
+    -- Input_sounding with moisture new HSE balance (#1469)
+
+    -- SAM microphysics re-write for ERF dycore (#1461/#1447)
+
+    -- Two-way coupled AMR (#1458)
+
+    -- Numerical diffusion bug fix (#1451)
+
 # 24.02
     -- AMReX submodule set to 24.02 release hash (52393d3)
 

Docs/sphinx_doc/BestPractices.rst

@@ -0,0 +1,63 @@
+ .. role:: cpp(code)
+    :language: c++
+
+.. _GettingStarted:
+
+Best Practices
+==============
+
+Please note this section is a work in progress and aims to offer guidelines
+but, in general, your mileage may vary.
+
+Large-Eddy Simulations
+----------------------
+
+* Advection Scheme
+
+  - A WRF-like configuration is generally robust, but may under-resolve
+    turbulence features.
+
+    .. code-block:: python
+
+       erf.dycore_horiz_adv_type  = "Upwind_5th"
+       erf.dycore_vert_adv_type   = "Upwind_3rd"
+       erf.dryscal_horiz_adv_type = "Upwind_5th"
+       erf.dryscal_vert_adv_type  = "Upwind_3rd"
+
+  - Centered difference schemes will generally give some non-physical
+    numerical noise, clearly visible in the free atmosphere, but may also
+    better resolve turbulence features. With ``Centered_2nd``, the simulation
+    may remain numerically stable but without any upwinding or numerical
+    diffusion, these results should be carefully interpreted.
+
+  - For higher-order central differencing alone (i.e., without any added
+    upwinding), at least 5% numerical diffusion should be included to stabilize
+    the solution; this was tested with ``Centered_6th``. Note that this does not
+    necessarily kill the numerical noise and is only for numerical stability.
+    These options are identical to WRF's ``diff_6th_opt`` (default: off) and
+    ``diff_6th_factor`` (default: 12%) options.
+
+    .. code-block:: python
+
+       erf.use_NumDiff  = true
+       erf.NumDiffCoeff = 0.05
+
+* Time Integration
+
+  - Split timestepping offers some computational cost savings but still does
+    not allow you to run with an incompressible time-step size.
+  - The acoustic CFL should be less than 0.5, with 4--6 substeps. E.g.,:
+
+    .. code-block:: python
+
+       erf.fixed_dt           = 0.06  # slow timestep
+
+       # These are equivalent and result in a fixed fast timestep size
+       #   if dx=10, speed of sound ~ 350 m/s
+       erf.fixed_fast_dt      = 0.01  # ==> CFL~0.35
+       #erf.fixed_mri_dt_ratio = 6
+
+       # Alternatively, let ERF chose the fast timestep
+       #erf.cfl                = 0.5
+
+

Docs/sphinx_doc/Inputs.rst

@@ -1066,7 +1066,10 @@ Examples of Usage
 Moisture
 ========
 
-ERF has several different moisture models.
+ERF has several different moisture models. The models that are currently implemented
+are Eulerian models; however, ERF has the capability for Lagrangian models when
+compiled with particles.
+
 The following run-time options control how the full moisture model is used.
 
 List of Parameters

Docs/sphinx_doc/index.rst

@@ -39,6 +39,7 @@ In addition to this documentation, there is API documentation for ERF generated
 
    GettingStarted.rst
    Inputs.rst
+   BestPractices.rst
 
 .. toctree::
    :caption: THEORY

Source/DataStructs/DataStruct.H

@@ -28,6 +28,10 @@ enum struct TerrainType {
     Static, Moving
 };
 
+enum struct MoistureModelType{
+    Eulerian, Lagrangian, Undefined
+};
+
 enum struct MoistureType {
     Kessler, SAM, Kessler_NoRain, None
 };

Source/ERF.H

@@ -43,7 +43,8 @@
 #include "ParticleData.H"
 #endif
 
-#include "Microphysics.H"
+#include "EulerianMicrophysics.H"
+#include "LagrangianMicrophysics.H"
 #include "LandSurface.H"
 
 #ifdef ERF_USE_RRTMGP
@@ -374,6 +375,9 @@ private:
 
     void initialize_integrator (int lev, amrex::MultiFab& cons_mf, amrex::MultiFab& vel_mf);
 
+    // Initialize the microphysics object
+    void initializeMicrophysics (const int&);
+
     // Compute a vector of new MultiFabs by copying from valid region and filling ghost cells
     //
     // NOTE: FillPatch takes in an empty MF, and returns cell-centered + velocities (not momenta)
@@ -552,7 +556,7 @@ private:
     amrex::Vector<amrex::MultiFab> rW_old;
     amrex::Vector<amrex::MultiFab> rW_new;
 
-    Microphysics micro;
+    std::unique_ptr<Microphysics> micro;
     amrex::Vector<amrex::Vector<amrex::MultiFab*>> qmoist; // (lev,ncomp) This has up to 8 components: qt, qv, qc, qi, qp, qr, qs, qg
 
 //#if defined(ERF_USE_WINDFARM)

Source/ERF.cpp

@@ -111,10 +111,6 @@ ERF::ERF ()
 
     int nlevs_max = max_level + 1;
 
-    // NOTE: size micro before readparams (chooses the model at all levels)
-    micro.ReSize(nlevs_max);
-    qmoist.resize(nlevs_max);
-
 #ifdef ERF_USE_WINDFARM
     if(solverChoice.windfarm_type == WindFarmType::Fitch){
         Nturb.resize(nlevs_max);
@@ -132,6 +128,8 @@ ERF::ERF ()
     lsm_flux.resize(nlevs_max);
 
     ReadParameters();
+    initializeMicrophysics(nlevs_max);
+
     const std::string& pv1 = "plot_vars_1"; setPlotVariables(pv1,plot_var_names_1);
     const std::string& pv2 = "plot_vars_2"; setPlotVariables(pv2,plot_var_names_2);
 
@@ -609,11 +607,23 @@ ERF::InitData ()
         if (solverChoice.moisture_type != MoistureType::None)
         {
             for (int lev = 0; lev <= finest_level; lev++) {
-                FillPatchMoistVars(lev, *(qmoist[lev][0])); // qv component
+                if (qmoist[lev].size() > 0) FillPatchMoistVars(lev, *(qmoist[lev][0])); // qv component
             }
         }
     }
 
+#ifdef ERF_USE_PARTICLES
+    /* If using a Lagrangian microphysics model, its particle container has now been
+       constructed and initialized (calls to micro->Init). So, add its pointer to
+       ERF::particleData and remove its name from list of unallocated particle containers. */
+    if (Microphysics::modelType(solverChoice.moisture_type) == MoistureModelType::Lagrangian) {
+        const auto& pc_name( dynamic_cast<LagrangianMicrophysics&>(*micro).getName() );
+        const auto& pc_ptr( dynamic_cast<LagrangianMicrophysics&>(*micro).getParticleContainer() );
+        particleData.pushBack(pc_name, pc_ptr);
+        particleData.getNamesUnalloc().remove(pc_name);
+    }
+#endif
+
     if (input_bndry_planes) {
         // Read the "time.dat" file to know what data is available
         m_r2d->read_time_file();
@@ -801,7 +811,7 @@ ERF::InitData ()
 
     // Update micro vars before first plot file
     if (solverChoice.moisture_type != MoistureType::None) {
-        for (int lev = 0; lev <= finest_level; ++lev) micro.Update_Micro_Vars_Lev(lev, vars_new[lev][Vars::cons]);
+        for (int lev = 0; lev <= finest_level; ++lev) micro->Update_Micro_Vars_Lev(lev, vars_new[lev][Vars::cons]);
     }
 
     // check for additional plotting variables that are available after particle containers
@@ -912,6 +922,34 @@ ERF::InitData ()
     BL_PROFILE_VAR_STOP(InitData);
 }
 
+// Initialize microphysics object
+void
+ERF::initializeMicrophysics (const int& a_nlevsmax /*!< number of AMR levels */)
+{
+    if (Microphysics::modelType(solverChoice.moisture_type) == MoistureModelType::Eulerian) {
+
+        micro = std::make_unique<EulerianMicrophysics>(a_nlevsmax, solverChoice.moisture_type);
+
+    } else if (Microphysics::modelType(solverChoice.moisture_type) == MoistureModelType::Lagrangian) {
+#ifdef ERF_USE_PARTICLES
+
+        micro = std::make_unique<LagrangianMicrophysics>(a_nlevsmax, solverChoice.moisture_type);
+        /* Lagrangian microphysics models will have a particle container; it needs to be added
+           to ERF::particleData */
+        const auto& pc_name( dynamic_cast<LagrangianMicrophysics&>(*micro).getName() );
+        /* The particle container has not yet been constructed and initialized, so just add
+           its name here for now (so that functions to set plotting variables can see it). */
+        particleData.addName( pc_name );
+
+#else
+        amrex::Abort("Lagrangian microphysics can be used when compiled with ERF_USE_PARTICLES");
+#endif
+    }
+
+    qmoist.resize(a_nlevsmax);
+    return;
+}
+
 void
 ERF::restart ()
 {
@@ -1227,22 +1265,6 @@ ERF::ReadParameters ()
 
     solverChoice.init_params(max_level);
 
-    // What type of moisture model to use
-    // NOTE: Must be checked after init_params
-    if (solverChoice.moisture_type == MoistureType::SAM) {
-        micro.SetModel<SAM>();
-        amrex::Print() << "SAM moisture model!\n";
-    } else if (solverChoice.moisture_type == MoistureType::Kessler or
-               solverChoice.moisture_type == MoistureType::Kessler_NoRain) {
-        micro.SetModel<Kessler>();
-        amrex::Print() << "Kessler moisture model!\n";
-    } else if (solverChoice.moisture_type == MoistureType::None) {
-        micro.SetModel<NullMoist>();
-        amrex::Print() << "WARNING: Compiled with moisture but using NullMoist model!\n";
-    } else {
-        amrex::Abort("Dont know this moisture_type!") ;
-    }
-
     // What type of land surface model to use
     // NOTE: Must be checked after init_params
     if (solverChoice.lsm_type == LandSurfaceType::SLM) {
@@ -1307,12 +1329,13 @@ ERF::MakeHorizontalAverages ()
             auto  fab_arr = mf.array(mfi);
             auto const cons_arr = vars_new[lev][Vars::cons].const_array(mfi);
             auto const   qv_arr = qmoist[lev][0]->const_array(mfi);
+            int ncomp = vars_new[lev][Vars::cons].nComp();
 
             ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
                 Real dens = cons_arr(i, j, k, Rho_comp);
                 fab_arr(i, j, k, 2) = getPgivenRTh(cons_arr(i, j, k, RhoTheta_comp), qv_arr(i,j,k));
-                fab_arr(i, j, k, 3) = cons_arr(i, j, k, RhoQ1_comp) / dens;
-                fab_arr(i, j, k, 4) = cons_arr(i, j, k, RhoQ2_comp) / dens;
+                fab_arr(i, j, k, 3) = (ncomp > RhoQ1_comp ? cons_arr(i, j, k, RhoQ1_comp) / dens : 0.0);
+                fab_arr(i, j, k, 4) = (ncomp > RhoQ2_comp ? cons_arr(i, j, k, RhoQ2_comp) / dens : 0.0);
             });
         }
 
@@ -1604,10 +1627,6 @@ ERF::ERF (const amrex::RealBox& rb, int max_level_in,
 
     int nlevs_max = max_level + 1;
 
-    // NOTE: size micro before readparams (chooses the model at all levels)
-    micro.ReSize(nlevs_max);
-    qmoist.resize(nlevs_max);
-
 #ifdef ERF_USE_WINDFARM
     if(solverChoice.windfarm_type == WindFarmType::Fitch){
         Nturb.resize(nlevs_max);
@@ -1625,6 +1644,8 @@ ERF::ERF (const amrex::RealBox& rb, int max_level_in,
     lsm_flux.resize(nlevs_max);
 
     ReadParameters();
+    initializeMicrophysics(nlevs_max);
+
     const std::string& pv1 = "plot_vars_1"; setPlotVariables(pv1,plot_var_names_1);
     const std::string& pv2 = "plot_vars_2"; setPlotVariables(pv2,plot_var_names_2);
 

Source/ERF_make_new_level.cpp

@@ -44,7 +44,7 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba,
 
     init_stuff(lev, ba, dm);
 
-    int n_qstate   = micro.Get_Qstate_Size();
+    int n_qstate   = micro->Get_Qstate_Size();
     int ncomp_cons = NVAR_max - (NMOIST_max - n_qstate);
 
     lev_new[Vars::cons].define(ba, dm, ncomp_cons, ngrow_state);
@@ -62,15 +62,15 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba,
     //********************************************************************************************
     // Microphysics
     // *******************************************************************************************
-    int q_size  = micro.Get_Qmoist_Size();
+    int q_size  = micro->Get_Qmoist_Size(lev);
     qmoist[lev].resize(q_size);
-    micro.Define(lev, solverChoice);
+    micro->Define(lev, solverChoice);
     if (solverChoice.moisture_type != MoistureType::None)
     {
-        micro.Init(lev, vars_new[lev][Vars::cons], grids[lev], Geom(lev), 0.0); // dummy dt value
+        micro->Init(lev, vars_new[lev][Vars::cons], grids[lev], Geom(lev), 0.0); // dummy dt value
     }
     for (int mvar(0); mvar<qmoist[lev].size(); ++mvar) {
-        qmoist[lev][mvar] = micro.Get_Qmoist_Ptr(lev,mvar);
+        qmoist[lev][mvar] = micro->Get_Qmoist_Ptr(lev,mvar);
     }
 
     //********************************************************************************************
@@ -225,15 +225,15 @@ ERF::MakeNewLevelFromCoarse (int lev, Real time, const BoxArray& ba,
     //********************************************************************************************
     // Microphysics
     // *******************************************************************************************
-    int q_size  = micro.Get_Qmoist_Size();
+    int q_size  = micro->Get_Qmoist_Size(lev);
     qmoist[lev].resize(q_size);
-    micro.Define(lev, solverChoice);
+    micro->Define(lev, solverChoice);
     if (solverChoice.moisture_type != MoistureType::None)
     {
-        micro.Init(lev, vars_new[lev][Vars::cons], grids[lev], Geom(lev), 0.0); // dummy dt value
+        micro->Init(lev, vars_new[lev][Vars::cons], grids[lev], Geom(lev), 0.0); // dummy dt value
     }
     for (int mvar(0); mvar<qmoist[lev].size(); ++mvar) {
-        qmoist[lev][mvar] = micro.Get_Qmoist_Ptr(lev,mvar);
+        qmoist[lev][mvar] = micro->Get_Qmoist_Ptr(lev,mvar);
     }
 
     init_stuff(lev, ba, dm);
@@ -340,15 +340,15 @@ ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapp
     //********************************************************************************************
     // Microphysics
     // *******************************************************************************************
-    int q_size = micro.Get_Qmoist_Size();
+    int q_size = micro->Get_Qmoist_Size(lev);
     qmoist[lev].resize(q_size);
-    micro.Define(lev, solverChoice);
+    micro->Define(lev, solverChoice);
     if (solverChoice.moisture_type != MoistureType::None)
     {
-        micro.Init(lev, vars_new[lev][Vars::cons], grids[lev], Geom(lev), 0.0); // dummy dt value
+        micro->Init(lev, vars_new[lev][Vars::cons], grids[lev], Geom(lev), 0.0); // dummy dt value
     }
     for (int mvar(0); mvar<qmoist[lev].size(); ++mvar) {
-        qmoist[lev][mvar] = micro.Get_Qmoist_Ptr(lev,mvar);
+        qmoist[lev][mvar] = micro->Get_Qmoist_Ptr(lev,mvar);
     }
 
     init_stuff(lev, ba, dm);