Need to do SAM still.

AMLattanzi · Dec 13, 2023 · 178d21f · 178d21f
1 parent 6bbf92f
commit 178d21f
Show file tree

Hide file tree

Showing 14 changed files with 464 additions and 94 deletions.
diff --git a/Source/ERF.cpp b/Source/ERF.cpp
@@ -157,6 +157,7 @@ ERF::ERF ()
         vars_old[lev].resize(Vars::NumTypes);
     }
 
+    micro.ReSize(nlevs_max);
     qmoist.resize(nlevs_max);
 
     mri_integrator_mem.resize(nlevs_max);
@@ -578,20 +579,17 @@ ERF::InitData ()
         }
     }
 
-    // Initialize microphysics here
-    micro.define(solverChoice);
 
-    // Call Init which will call Diagnose to fill qmoist
-    if (solverChoice.moisture_type != MoistureType::None)
+    for (int lev = 0; lev <= finest_level; ++lev)
     {
-        for (int lev = 0; lev <= finest_level; ++lev)
+        // Initialize microphysics here
+        micro.define(lev, solverChoice);
+
+        // If not restarting we need to fill qmoist given qt and qp.
+        if (solverChoice.moisture_type != MoistureType::None)
         {
-            // If not restarting we need to fill qmoist given qt and qp.
-            if (restart_chkfile.empty()) {
-                micro.Init(vars_new[lev][Vars::cons], *(qmoist[lev]),
-                           grids[lev], Geom(lev), 0.0); // dummy value, not needed just to diagnose
-                micro.Update(vars_new[lev][Vars::cons], *(qmoist[lev]));
-            }
+            micro.Init(lev, vars_new[lev][Vars::cons], *(qmoist[lev]), grids[lev], Geom(lev), 0.0); // dummy value, not needed just to diagnose
+            micro.Update_Micro_Vars_Lev(lev, vars_new[lev][Vars::cons]);
         }
     }
 
@@ -1504,6 +1502,7 @@ ERF::ERF (const amrex::RealBox& rb, int max_level_in,
     rV_old.resize(nlevs_max);
     rW_old.resize(nlevs_max);
 
+    micro.ReSize(nlevs_max);
     qmoist.resize(nlevs_max);
 
     mri_integrator_mem.resize(nlevs_max);

diff --git a/Source/Microphysics/FastEddy/FastEddy.H b/Source/Microphysics/FastEddy/FastEddy.H
@@ -70,11 +70,32 @@ public:
           const amrex::Geometry& geom,
           const amrex::Real& dt_advance) override;
 
+    // Copy state into micro vars
+    void
+    Copy_State_to_Micro (const amrex::MultiFab& cons_in) override;
+
+    // Copy micro into state vars
+    void
+    Copy_Micro_to_State (amrex::MultiFab& cons_in) override;
+
     // update ERF variables
     void
     Update (amrex::MultiFab& cons_in,
             amrex::MultiFab& qmoist) override;
 
+    void
+    Update_Micro_Vars (amrex::MultiFab& cons_in) override
+    {
+        this->Copy_State_to_Micro(cons_in);
+        this->Diagnose();
+    }
+
+    void
+    Update_State_Vars (amrex::MultiFab& cons_in) override
+    {
+        this->Copy_Micro_to_State(cons_in);
+    }
+
     // wrapper to do all the updating
     void
     Advance ( ) override
@@ -83,6 +104,13 @@ public:
         this->Diagnose();
     }
 
+    amrex::MultiFab*
+    Qmoist_Ptr (const int& varIdx) override
+    {
+        AMREX_ALWAYS_ASSERT(varIdx <= m_qmoist_size);
+        return mic_fab_vars[varIdx].get();
+    }
+
     int
     Qmoist_Size ( ) override { return FastEddy::m_qmoist_size; }
 

diff --git a/Source/Microphysics/FastEddy/FastEddy.cpp b/Source/Microphysics/FastEddy/FastEddy.cpp
@@ -46,8 +46,8 @@ void FastEddy::AdvanceFE () {
             dq_vapor_to_clwater = 0.0;
             dq_clwater_to_vapor = 0.0;
 
-            //Real fac             = qsat*4093.0*L_v/(Cp_d*std::pow(tabs_array(i,j,k)-36.0,2));
-            Real fac             = qsat*L_v*L_v/(Cp_d*R_v*tabs_array(i,j,k)*tabs_array(i,j,k));
+            //Real fac = qsat*4093.0*L_v/(Cp_d*std::pow(tabs_array(i,j,k)-36.0,2));
+            Real fac = qsat*L_v*L_v/(Cp_d*R_v*tabs_array(i,j,k)*tabs_array(i,j,k));
 
             // If water vapor content exceeds saturation value, then vapor condenses to waterm and latent heat is released, increasing temperature
             if(qv_array(i,j,k) > qsat){

diff --git a/Source/Microphysics/FastEddy/Init_FE.cpp b/Source/Microphysics/FastEddy/Init_FE.cpp
@@ -8,6 +8,89 @@
 
 using namespace amrex;
 
+/**
+ * Initializes the Microphysics module.
+ *
+ * @param[in] cons_in Conserved variables input
+ * @param[in] qc_in Cloud variables input
+ * @param[in,out] qv_in Vapor variables input
+ * @param[in] qi_in Ice variables input
+ * @param[in] grids The boxes on which we will evolve the solution
+ * @param[in] geom Geometry associated with these MultiFabs and grids
+ * @param[in] dt_advance Timestep for the advance
+ */
+void FastEddy::Init (const MultiFab& cons_in, MultiFab& qmoist,
+                     const BoxArray& grids,
+                     const Geometry& geom,
+                     const Real& dt_advance)
+{
+    dt = dt_advance;
+    m_geom = geom;
+    m_gtoe = grids;
+
+    // initialize microphysics variables
+    for (auto ivar = 0; ivar < MicVar_FE::NumVars; ++ivar) {
+        mic_fab_vars[ivar] = std::make_shared<MultiFab>(cons_in.boxArray(), cons_in.DistributionMap(), 1, cons_in.nGrowVect());
+        mic_fab_vars[ivar]->setVal(0.);
+    }
+
+    // Set class data members
+    for ( MFIter mfi(cons_in, TileNoZ()); mfi.isValid(); ++mfi) {
+
+        const auto& box3d = mfi.tilebox();
+
+        const auto& lo = amrex::lbound(box3d);
+        const auto& hi = amrex::ubound(box3d);
+
+        nlev = box3d.length(2);
+        zlo  = lo.z;
+        zhi  = hi.z;
+    }
+}
+
+/**
+ * Initializes the Microphysics module.
+ *
+ * @param[in] cons_in Conserved variables input
+ */
+void FastEddy::Copy_State_to_Micro (const MultiFab& cons_in)
+{
+    // Get the temperature, density, theta, qt and qp from input
+    for ( MFIter mfi(cons_in); mfi.isValid(); ++mfi) {
+        const auto& box3d = mfi.tilebox();
+
+        auto states_array = cons_in.array(mfi);
+
+        auto qv_array     = mic_fab_vars[MicVar_FE::qv]->array(mfi);
+        auto qc_array     = mic_fab_vars[MicVar_FE::qc]->array(mfi);
+        auto rho_array    = mic_fab_vars[MicVar_FE::rho]->array(mfi);
+        auto theta_array  = mic_fab_vars[MicVar_FE::theta]->array(mfi);
+        auto temp_array   = mic_fab_vars[MicVar_FE::tabs]->array(mfi);
+
+        // Get pressure, theta, temperature, density, and qt, qp
+        amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
+        {
+            rho_array(i,j,k)   = states_array(i,j,k,Rho_comp);
+            theta_array(i,j,k) = states_array(i,j,k,RhoTheta_comp)/states_array(i,j,k,Rho_comp);
+            qv_array(i,j,k)    = states_array(i,j,k,RhoQ1_comp)/states_array(i,j,k,Rho_comp);
+            qc_array(i,j,k)    = states_array(i,j,k,RhoQ2_comp)/states_array(i,j,k,Rho_comp);
+            temp_array(i,j,k)  = getTgivenRandRTh(states_array(i,j,k,Rho_comp),states_array(i,j,k,RhoTheta_comp), qv_array(i,j,k));
+        });
+    }
+}
+
+
+
+#if 0
+#include <AMReX_GpuContainers.H>
+#include "Microphysics.H"
+#include "IndexDefines.H"
+#include "PlaneAverage.H"
+#include "EOS.H"
+#include "TileNoZ.H"
+
+using namespace amrex;
+
 /**
  * Initializes the Microphysics module.
  *
@@ -75,3 +158,4 @@ void FastEddy::Init (const MultiFab& cons_in, MultiFab& qmoist,
         });
     }
 }
+#endif
diff --git a/Source/Microphysics/FastEddy/Update_FE.cpp b/Source/Microphysics/FastEddy/Update_FE.cpp
@@ -38,4 +38,38 @@ void FastEddy::Update (amrex::MultiFab& cons,
     qmoist.FillBoundary(m_geom.periodicity());
 }
 
+/**
+ * Updates conserved and microphysics variables in the provided MultiFabs from
+ * the internal MultiFabs that store Microphysics module data.
+ *
+ * @param[out] cons Conserved variables
+ * @param[out] qmoist: qv, qc, qi, qr, qs, qg
+ */
+void FastEddy::Copy_Micro_to_State (amrex::MultiFab& cons)
+{
+
+    // Get the temperature, density, theta, qt and qp from input
+    for (amrex::MFIter mfi(cons,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
+        const auto& box3d = mfi.tilebox();
+
+        auto states_arr = cons.array(mfi);
+
+        auto rho_arr    = mic_fab_vars[MicVar_FE::rho]->array(mfi);
+        auto theta_arr  = mic_fab_vars[MicVar_FE::theta]->array(mfi);
+        auto qv_arr     = mic_fab_vars[MicVar_FE::qv]->array(mfi);
+        auto qc_arr     = mic_fab_vars[MicVar_FE::qc]->array(mfi);
+
+        // get potential total density, temperature, qt, qp
+        amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
+        {
+            states_arr(i,j,k,RhoTheta_comp) = rho_arr(i,j,k)*theta_arr(i,j,k);
+            states_arr(i,j,k,RhoQ1_comp)    = rho_arr(i,j,k)*qv_arr(i,j,k);
+            states_arr(i,j,k,RhoQ2_comp)    = rho_arr(i,j,k)*qc_arr(i,j,k);
+        });
+    }
+
+    // Fill interior ghost cells and periodic boundaries
+    cons.FillBoundary(m_geom.periodicity());
+}
+
 
diff --git a/Source/Microphysics/Kessler/Diagnose_Kessler.cpp b/Source/Microphysics/Kessler/Diagnose_Kessler.cpp
@@ -1,4 +1,5 @@
 #include "Microphysics.H"
+#include "EOS.H"
 
 /**
  * Computes diagnostic quantities like cloud ice/liquid and precipitation ice/liquid
@@ -40,3 +41,34 @@ void Kessler::Diagnose () {
         });
     }
 }
+
+
+/**
+ * Computes diagnostic quantities like temperature and pressure
+ * from the existing Microphysics variables.
+ */
+void Kessler::Compute_Tabs_and_Pres () {
+
+    auto rho   = mic_fab_vars[MicVar_Kess::rho];
+    auto theta = mic_fab_vars[MicVar_Kess::theta];
+    auto qv    = mic_fab_vars[MicVar_Kess::qv];
+    auto tabs  = mic_fab_vars[MicVar_Kess::tabs];
+    auto pres  = mic_fab_vars[MicVar_Kess::pres];
+
+    for ( amrex::MFIter mfi(*tabs, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
+        const auto& box3d = mfi.tilebox();
+
+        auto rho_array   = rho->array(mfi);
+        auto theta_array = theta->array(mfi);
+        auto qv_array    = qv->array(mfi);
+        auto tabs_array  = tabs->array(mfi);
+        auto pres_array  = pres->array(mfi);
+
+        amrex::ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
+        {
+            amrex::Real RhoTheta = rho_array(i,j,k)*theta_array(i,j,k);
+            tabs_array(i,j,k)  = getTgivenRandRTh(rho_array(i,j,k), RhoTheta, qv_array(i,j,k));
+            pres_array(i,j,k)  = getPgivenRTh(RhoTheta, qv_array(i,j,k))/100.;
+        });
+    }
+}