SDM Refactor

Source/ERF_Constants.H

@@ -15,6 +15,7 @@ constexpr amrex::Real eightthirds  = amrex::Real(8.0/3.0);
 
 constexpr amrex::Real PI           = amrex::Real(3.14159265358979323846264338327950288);
 constexpr amrex::Real PIoTwo       = PI/two;
+constexpr amrex::Real four_thirds_pi = four/three*PI;
 
 // Physical Constants
 constexpr amrex::Real R_d          = amrex::Real(287.0);    // dry air constant for dry air [J/(kg-K)]

Source/MaterialProperties/ERF_TerminalVelocity.H

@@ -8,23 +8,38 @@
 template <typename RT /*!< real-type */>
 struct TerminalVelocity
 {
-    RT m_rho;
+    RT m_rho_w; /*!< density of water */
+
+    /*! \brief Compute viscosity coefficient - Pruppacher & Klett, 1997*/
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    auto viscCoeff( const RT a_T /*!< temperature [K] */ ) const noexcept
+    {
+        RT T_degC = a_T - RT(273.15); // [K] -> [deg Celsius]
+        RT viscosity = RT(0.0);
+        if (T_degC >= RT(0)) {
+            viscosity = (RT(1.7180) + RT(4.9e-3)*T_degC) * RT(1.0e-4);
+        } else {
+            viscosity = (  RT(1.7180) + RT(4.9e-3)*T_degC
+                         - RT(1.2e-5)*T_degC*T_degC     ) * RT(1.0e-4);
+        }
+        return viscosity;
+    }
 
     /*! \brief Terminal velocity - Rogers & Yau, 1989 */
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     RT RogersYau( const RT a_r) const noexcept
     {
-        static const RT k1 = amrex::Real(1.233e8); // m^{-1} s^{-1}
+        static const RT k1 = RT(1.233e8); // m^{-1} s^{-1}
         return (k1 * a_r * a_r);
     }
 
     /*! \brief Terminal velocity - Atlas & Ulbrich, 1977 */
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     RT AtlasUlbrich( const RT a_r) const noexcept
     {
-        static const RT a = amrex::Real(3.778);
-        static const RT b = amrex::Real(0.67);
-        RT D = 2*a_r*1000; // meter -> mm
+        static const RT a = RT(3.778);
+        static const RT b = RT(0.67);
+        RT D = RT(2)*a_r*RT(1000); // meter -> mm
         return a*(std::exp(b*std::log(D)));
     }
 
@@ -41,63 +56,56 @@ struct TerminalVelocity
                         const RT a_P,      /*!< pressure */
                         const RT a_T       /*!< temperature */) const noexcept
     {
-        RT lb4l = amrex::Real(6.62E-6);    // base length[cm] for mean free path  of air molecules
-        RT visb4l = amrex::Real(1.818E-4); // base dynamic viscosity[g/(cm*s)] for mean free path of air molecules
-        RT pb4l = amrex::Real(1013.25);    // base pressure[hPa] for mean free path of air molecules
-        RT tb4l = amrex::Real(293.15);     // base temperature[20C] for mean free path of air molecules
+        RT lb4l = RT(6.62E-6);    // base length[cm] for mean free path  of air molecules
+        RT visb4l = RT(1.818E-4); // base dynamic viscosity[g/(cm*s)] for mean free path of air molecules
+        RT pb4l = RT(1013.25);    // base pressure[hPa] for mean free path of air molecules
+        RT tb4l = RT(293.15);     // base temperature[20C] for mean free path of air molecules
 
-        RT P_hPa = a_P / amrex::Real(100.0);           // [Pa] -> [hPa]
-        RT T_degC = a_T - amrex::Real(273.15);         // [K] -> [deg Celsius]
-        RT diameter_cm = two*a_r*amrex::Real(100.0);   // [m] -> [cm]
-        RT rho_mat_cgs = m_rho/amrex::Real(1000.0);    // [kg/m^3] -> [g/cm^3]
-        RT rho_air_cgs = a_rho/amrex::Real(1000.0);    // [kg/m^3] -> [g/cm^3]
-        RT grav_cgs = CONST_GRAV * 100;   // [m/s^2] -> [cm/s^2]
+        RT P_hPa = a_P / RT(100.0);           // [Pa] -> [hPa]
+        RT diameter_cm = RT(2.0)*a_r*RT(100.0);   // [m] -> [cm]
+        RT rho_mat_cgs = m_rho_w/RT(1000.0);  // [kg/m^3] -> [g/cm^3]
+        RT rho_air_cgs = a_rho/RT(1000.0);    // [kg/m^3] -> [g/cm^3]
+        RT grav_cgs = CONST_GRAV * RT(100);   // [m/s^2] -> [cm/s^2]
 
         RT gxdrow = grav_cgs * ( rho_mat_cgs - rho_air_cgs );
 
-        RT viscosity = zero;
-        if (T_degC >= 0) {
-            viscosity = (amrex::Real(1.7180) + amrex::Real(4.9e-3)*T_degC) * amrex::Real(1.0e-4);
-        } else {
-            viscosity = (  amrex::Real(1.7180) + amrex::Real(4.9e-3)*T_degC
-                         - amrex::Real(1.2e-5)*T_degC*T_degC     ) * amrex::Real(1.0e-4);
-        }
+        RT viscosity = viscCoeff(a_T);
 
-        RT v_terminal = zero;
-        if (diameter_cm < amrex::Real(1.9e-3)) {
+        RT v_terminal = RT(0.0);
+        if (diameter_cm < RT(1.9e-3)) {
             // small cloud droplets
             RT sd_l = lb4l * (viscosity/visb4l) * (pb4l/P_hPa) * std::sqrt(a_T/tb4l);
-            RT c1  = gxdrow / (amrex::Real(18.0)*viscosity);
-            RT csc = one + amrex::Real(2.510) * (sd_l/diameter_cm);
+            RT c1  = gxdrow / (RT(18.0)*viscosity);
+            RT csc = RT(1.0) + RT(2.510) * (sd_l/diameter_cm);
             v_terminal = c1 * csc * diameter_cm * diameter_cm;
-        } else if (diameter_cm < amrex::Real(1.07e-1)) {
+        } else if (diameter_cm < RT(1.07e-1)) {
             // large cloud droplets and small raindrops
             RT sd_l = lb4l * (viscosity/visb4l) * (pb4l/P_hPa) * std::sqrt(a_T/tb4l);
-            RT c2 = rho_air_cgs * (amrex::Real(4.0)*gxdrow)/(three*viscosity*viscosity);
+            RT c2 = rho_air_cgs * (RT(4.0)*gxdrow)/(RT(three)*viscosity*viscosity);
             RT nda = c2 * diameter_cm * diameter_cm * diameter_cm;
-            RT csc = one + amrex::Real(2.510) * ( sd_l / diameter_cm );
+            RT csc = RT(1.0) + RT(2.510) * ( sd_l / diameter_cm );
             RT nre = csc * std::exp(vz_b_x(std::log(nda)));
             v_terminal = (viscosity*nre)/(rho_air_cgs*diameter_cm);
         } else {
             // large raindrops
-            RT tau = one - a_T/amrex::Real(647.0960);
-            RT sigma = amrex::Real(0.2358) * std::exp( amrex::Real(1.2560)*std::log(tau) )
-                                * ( one - amrex::Real(0.6250)*tau );
-            if( a_T<amrex::Real(267.5) ) {
-                tau = std::tanh( (a_T-amrex::Real(243.9))/amrex::Real(35.35) );
-                sigma = sigma - amrex::Real(2.854E-3) * tau + amrex::Real(1.666E-3);
+            RT tau = RT(1.0) - a_T/RT(647.0960);
+            RT sigma = RT(0.2358) * std::exp( RT(1.2560)*std::log(tau) )
+                                * ( RT(1.0) - RT(0.6250)*tau );
+            if( a_T<RT(267.5) ) {
+                tau = std::tanh( (a_T-RT(243.9))/RT(35.35) );
+                sigma = sigma - RT(2.854E-3) * tau + RT(1.666E-3);
             }
-            sigma = sigma * amrex::Real(1.E+3); // N/m => g/s^2
-            RT c3 = (amrex::Real(4.0)*gxdrow)/(three*sigma);
+            sigma = sigma * RT(1.E+3); // N/m => g/s^2
+            RT c3 = (RT(4.0)*gxdrow)/(RT(three)*sigma);
             RT bond = c3 * diameter_cm * diameter_cm;
             RT npp =    (rho_air_cgs*rho_air_cgs) * (sigma*sigma*sigma)
                         / (gxdrow*viscosity*viscosity*viscosity*viscosity);
-            npp = std::exp( std::log(npp)/amrex::Real(6.0) );
+            npp = std::exp( std::log(npp)/RT(6.0) );
             RT nre = npp * std::exp(vz_c_x(std::log(bond*npp)));
             v_terminal = (viscosity*nre)/(rho_air_cgs*diameter_cm);
         }
 
-        v_terminal /= amrex::Real(100.0); // [cm/s] -> [m/s]
+        v_terminal /= RT(100.0); // [cm/s] -> [m/s]
         return v_terminal;
     }
 
@@ -110,13 +118,13 @@ struct TerminalVelocity
         RT x1 = a_x;
         RT x2 = x1 * x1;
         RT x3 = x1 * x2;
-        RT y = - amrex::Real(0.318657E+1)
-               + amrex::Real(0.9926960)    * x1
-               - amrex::Real(0.153193E-2)  * x2
-               - amrex::Real(0.987059E-3)  * x3
-               - amrex::Real(0.578878e-3)  * x2*x2
-               + amrex::Real(0.855176E-4)  * x3*x2
-               - amrex::Real(0.327815E-5)  * x3*x3;
+        RT y = - RT(0.318657E+1)
+               + RT(0.9926960)    * x1
+               - RT(0.153193E-2)  * x2
+               - RT(0.987059E-3)  * x3
+               - RT(0.578878e-3)  * x2*x2
+               + RT(0.855176E-4)  * x3*x2
+               - RT(0.327815E-5)  * x3*x3;
         return y;
     }
 
@@ -127,12 +135,12 @@ struct TerminalVelocity
         RT x1 = a_x;
         RT x2 = x1 * x1;
         RT x3 = x1 * x2;
-        RT y = - amrex::Real(0.500015E+1)
-               + amrex::Real(0.523778E+1)  * x1
-               - amrex::Real(0.204914E+1)  * x2
-               + amrex::Real(0.47529400)   * x3
-               - amrex::Real(0.542819E-1)  * x2*x2
-               + amrex::Real(0.238449E-2)  * x3*x2;
+        RT y = - RT(0.500015E+1)
+               + RT(0.523778E+1)  * x1
+               - RT(0.204914E+1)  * x2
+               + RT(0.47529400)   * x3
+               - RT(0.542819E-1)  * x2*x2
+               + RT(0.238449E-2)  * x3*x2;
         return y;
     }
 

Source/Microphysics/SuperDropletsMoist/ERF_SuperDropletsMoist.H

@@ -313,66 +313,36 @@ class SuperDropletsMoist : public NullMoistLagrangian {
                                  amrex::MultiFab& a_mf )  const override
         {
             AMREX_ASSERT(a_mf.nComp() >= 1);
-            for (int v = 1; v < m_num_species; v++) {
-                {
-                    std::string varname = "qv_" + amrex::getEnumNameString(m_species[v]);
-                    if (a_name == varname) {
-                        amrex::MultiFab::Copy( a_mf,
-                                        *m_mic_fab_vars[s_qv_idx(v)],
-                                        0, 0, 1, amrex::IntVect::TheZeroVector() );
-                        return;
-                    }
-                }
-                {
-                    std::string varname = "qc_" + amrex::getEnumNameString(m_species[v]);
-                    if (a_name == varname) {
-                        amrex::MultiFab::Copy( a_mf,
-                                        *m_mic_fab_vars[s_qc_idx(v)],
-                                        0, 0, 1, amrex::IntVect::TheZeroVector() );
-                        return;
-                    }
-                }
-                {
-                    std::string varname = "qt_" + amrex::getEnumNameString(m_species[v]);
-                    if (a_name == varname) {
-                        amrex::MultiFab::Copy( a_mf,
-                                        *m_mic_fab_vars[s_qt_idx(v)],
-                                        0, 0, 1, amrex::IntVect::TheZeroVector() );
-                        return;
-                    }
-                }
-                {
-                    std::string varname = "sat_ratio_" + amrex::getEnumNameString(m_species[v]);
-                    if (a_name == varname) {
-                        amrex::MultiFab::Copy( a_mf,
-                                        *m_mic_fab_vars[s_sr_idx(v)],
-                                        0, 0, 1, amrex::IntVect::TheZeroVector() );
-                        return;
-                    }
+
+            // Helper to copy MultiFab if name matches
+            auto try_copy = [&](const std::string& prefix, int idx) -> bool {
+                if (a_name == prefix) {
+                    amrex::MultiFab::Copy(a_mf, *m_mic_fab_vars[idx],
+                                          0, 0, 1, amrex::IntVect::TheZeroVector());
+                    return true;
                 }
-                {
-                    std::string varname = "accum_" + amrex::getEnumNameString(m_species[v]);
-                    if (a_name == varname) {
-                        amrex::MultiFab::Copy( a_mf,
-                                        *m_mic_fab_vars[s_accum_idx(v)],
-                                        0, 0, 1, amrex::IntVect::TheZeroVector() );
-                        return;
-                    }
+                return false;
+            };
+
+            // Species variables
+            for (int v = 1; v < m_num_species; v++) {
+                std::string sp_name = amrex::getEnumNameString(m_species[v]);
+                if (try_copy("qv_" + sp_name, s_qv_idx(v)) ||
+                    try_copy("qc_" + sp_name, s_qc_idx(v)) ||
+                    try_copy("qt_" + sp_name, s_qt_idx(v)) ||
+                    try_copy("sat_ratio_" + sp_name, s_sr_idx(v)) ||
+                    try_copy("accum_" + sp_name, s_accum_idx(v))) {
+                    return;
                 }
             }
+
+            // Aerosol variables
             for (int v = 0; v < m_num_aerosols; v++) {
-                {
-                    std::string varname = "accum_" + amrex::getEnumNameString(m_aerosols[v]);
-                    if (a_name == varname) {
-                        amrex::MultiFab::Copy( a_mf,
-                                        *m_mic_fab_vars[a_accum_idx(m_num_species,v)],
-                                        0, 0, 1, amrex::IntVect::TheZeroVector() );
-                        return;
-                    }
-                }
+                std::string ae_name = amrex::getEnumNameString(m_aerosols[v]);
+                if (try_copy("accum_" + ae_name, a_accum_idx(m_num_species, v))) { return; }
             }
+
             amrex::Abort("SuperDropletsMoist::GetPlotVar() called with invalid name");
-            return;
         }
 
         AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE

Source/Microphysics/SuperDropletsMoist/ERF_SuperDropletsMoistUtils.cpp

@@ -326,14 +326,12 @@ void SuperDropletsMoist::ratioToDensity (MultiFab& a_var,
 void SuperDropletsMoist::computeQcQrWater ()
 {
     BL_PROFILE("SuperDropletsMoist::computeQcQrWater()");
-    m_super_droplets->speciesMassDensity( *(m_mic_fab_vars[MicVar_SD::q_c]),
-                                          m_idx_w,
-                                          0,
-                                          m_r_rain );
-    m_super_droplets->speciesMassDensity( *(m_mic_fab_vars[MicVar_SD::q_r]),
-                                          m_idx_w,
-                                          m_r_rain,
-                                          one );
+    m_super_droplets->cloudRainDensity( *(m_mic_fab_vars[MicVar_SD::q_c]),
+                                        0,
+                                        m_r_rain );
+    m_super_droplets->cloudRainDensity( *(m_mic_fab_vars[MicVar_SD::q_r]),
+                                        m_r_rain,
+                                        one );
 
     if (m_dimensionality == SDMSimulationDim::one_d_z) {
         for ( MFIter mfi(*m_mic_fab_vars[MicVar_SD::q_c]); mfi.isValid(); ++mfi) {

Source/Particles/ERFPC.H

@@ -333,6 +333,17 @@ class ERFPC : public amrex::ParticleContainer<  ERFParticlesRealIdxAoS::ncomps,
         // the following functions should ideally be private or protected, but need to be
         // public due to CUDA extended lambda capture rules
 
+        /*! \brief Helper for ParticleToMesh operations with non-uniform z support
+         *
+         * Generic CIC deposition from particles to a MultiFab. The caller provides
+         * a device lambda `value_func(ptd, i)` that returns the per-particle scalar
+         * to deposit (e.g. mass, number, flux). Handles both uniform and non-uniform
+         * vertical grids. Implementation in ERFPCParticleToMesh.H.
+         */
+        template<typename ValueFunc>
+        void ERFPCParticleToMesh(amrex::MultiFab& a_mf, int a_lev, int a_comp,
+                                  ValueFunc&& value_func) const;
+
         /*! Default particle initialization */
         void initializeParticlesUniformDistributionInBox (const std::unique_ptr<amrex::MultiFab>& a_ptr,
                                                           const amrex::RealBox& particle_box);