Enable non-zero mu_r

.gitignore

@@ -40,4 +40,7 @@
 .DS_Store
 
 # build directory
-build/
+build/
+
+# VS code settings directory
+.vscode

CMakeLists.txt

@@ -37,6 +37,7 @@ find_package(SUNDIALS 7.7.0 QUIET COMPONENTS arkode sunlinsollapackdense)
 if (NOT SUNDIALS_FOUND)
   include(FetchContent)
   set(FETCHCONTENT_QUIET OFF)
+  set(FETCHCONTENT_UPDATES_DISCONNECTED OFF)
   FetchContent_Declare(sundials
     GIT_REPOSITORY https://github.com/LLNL/sundials.git
     GIT_TAG v7.7.0
@@ -49,7 +50,7 @@ if (NOT SUNDIALS_FOUND)
   set(SUNDIALS_ENABLE_IDAS OFF)
   set(SUNDIALS_ENABLE_KINSOL OFF)
   set(SUNDIALS_ENABLE_LAPACK ON)
-  set(SUNDIALS_EXAMPLES_ENABLE_C OFF)
+  set(SUNDIALS_ENABLE_C_EXAMPLES OFF)
   set(BUILD_SHARED_LIBS OFF)
   FetchContent_MakeAvailable(sundials)
 endif()

FindNetCDF.cmake

@@ -50,6 +50,8 @@ if (netCDF_FOUND)
   endif ()
 
   FindNetCDF_get_is_parallel_aware("${NetCDF_INCLUDE_DIRS}")
+
+  target_link_libraries(NetCDF::NetCDF INTERFACE HDF5::HDF5)
   # Skip the rest of the logic in this file.
   return ()
 endif ()

rainshaft/evaporation.hpp

@@ -112,7 +112,7 @@ class Evaporation : public RainshaftProcess
     const double t2_term = 0.32 * cbrt(get_val(schmidt_num)) / sqrt_visc_over_rho;
     const double scaled_evap = t2_term * get_val(v_evap);
     const double t2 = scale_inv_lambdar + scaled_evap;
-    
+
     const double tau_inv = t1 * t2;
 
     if constexpr (WithGrad) {
@@ -307,7 +307,7 @@ class Evaporation : public RainshaftProcess
 std::array<RealOptGrad<WithGrad, 4>, 3> calc_evap(const RainshaftConstants& constants, const double t, const double q,
   const double nr, const double qr, const double p_mid, const RealOptGrad<WithGrad, 2> rho_dry, const RealOptGrad<WithGrad, 2> lambdar) const
   {
-    if (!regularize_qsat) 
+    if (!regularize_qsat)
     {
       // Skip the rest of this if no rain.
       if (qr < constants.qsmall)
@@ -337,7 +337,7 @@ std::array<RealOptGrad<WithGrad, 4>, 3> calc_evap(const RainshaftConstants& cons
       // regularization
       const RealOptGrad<WithGrad, 1> q_sat_dry = sat_form.q_sat_dry<WithGrad>(t, p_mid);
       const double epsilon_qsat = constants.epsilon_qsat_fac * get_val(q_sat_dry);
-      
+
       if (q <= get_val(q_sat_dry) - epsilon_qsat) { // this is the original value of q_evap
         const RealOptGrad<WithGrad, 1> diffusivity = calc_diffusivity<WithGrad>(t, p_mid);
         const RealOptGrad<WithGrad, 2> visc_over_rho = calc_visc_over_rho<WithGrad>(t, rho_dry);
@@ -370,21 +370,21 @@ std::array<RealOptGrad<WithGrad, 4>, 3> calc_evap(const RainshaftConstants& cons
           // the same factor
           const double dpoly_reg = 140.0 * pow(delta, 3) + 420.0 * pow(delta, 4) + 420.0 * pow(delta, 5) + 140.0 * pow(delta, 6);
           const RealOptGrad<WithGrad, 4> q_evap_reg = {get_val(q_evap) * poly_reg,
-                                                       {get_grad(q_evap)[0] * poly_reg + get_val(q_evap) * dpoly_reg * -get_grad(q_sat_dry)[0], 
-                                                        get_grad(q_evap)[1] * poly_reg + get_val(q_evap) * dpoly_reg, 
-                                                        get_grad(q_evap)[2] * poly_reg, 
+                                                       {get_grad(q_evap)[0] * poly_reg + get_val(q_evap) * dpoly_reg * -get_grad(q_sat_dry)[0],
+                                                        get_grad(q_evap)[1] * poly_reg + get_val(q_evap) * dpoly_reg,
+                                                        get_grad(q_evap)[2] * poly_reg,
                                                         get_grad(q_evap)[3] * poly_reg}};
           const double nr_over_qr = (nr + constants.qsmall * (1.e8)) / (qr + constants.qsmall);
           const RealOptGrad<WithGrad, 4> n_evap_reg = {get_val(q_evap_reg) * nr_over_qr,
-                                                       {get_grad(q_evap_reg)[0] * nr_over_qr, 
+                                                       {get_grad(q_evap_reg)[0] * nr_over_qr,
                                                         get_grad(q_evap_reg)[1] * nr_over_qr,
                                                         get_grad(q_evap_reg)[2] * nr_over_qr + get_val(q_evap_reg) / (qr + constants.qsmall),
                                                         get_grad(q_evap_reg)[3] * nr_over_qr - get_val(q_evap_reg) * (nr + constants.qsmall * (1.e8)) / pow(qr + constants.qsmall, 2)}};
           const double Lv_over_cp = -constants.latvap / constants.cp;
           const RealOptGrad<WithGrad, 4> t_evap_reg = {get_val(q_evap_reg) * Lv_over_cp,
-                                                       {get_grad(q_evap_reg)[0] * Lv_over_cp, 
-                                                        get_grad(q_evap_reg)[1] * Lv_over_cp, 
-                                                        get_grad(q_evap_reg)[2] * Lv_over_cp, 
+                                                       {get_grad(q_evap_reg)[0] * Lv_over_cp,
+                                                        get_grad(q_evap_reg)[1] * Lv_over_cp,
+                                                        get_grad(q_evap_reg)[2] * Lv_over_cp,
                                                         get_grad(q_evap_reg)[3] * Lv_over_cp}};
           return {q_evap_reg, n_evap_reg, t_evap_reg};
         } else {
@@ -445,47 +445,54 @@ std::array<RealOptGrad<WithGrad, 4>, 3> calc_evap(const RainshaftConstants& cons
     }
     // Factor converting D^3 to drop mass in grams.
     const double d3_to_gram = 1000. * constants.pi * constants.rhow / 6.;
-
-    // Integral for D <= 134.43 micron.
+    // for Vevapi (integral for D < 134.43 micron)
     const double int1_fac = sqrt(4579.5) * cbrt(d3_to_gram);
     const double lambdar_neg_2_5 = pow(lambdar, -2.5);
     const double lambdar_size1 = lambdar * 1.3443e-4;
-    const double term1 = int1_fac * tgamma_lower(3.5, lambdar_size1) * lambdar_neg_2_5;
-
-    // Integral for 134.43 micron < D <= 1511.64 micron.
+    // for Vevapi (integral for 134.43 micron < D < 1511.64 micron)
     const double int2_fac = sqrt(49.62) * pow(d3_to_gram, 1.0 / 6.0);
     const double lambdar_neg_2 = pow(lambdar, -2);
     const double lambdar_size2 = lambdar * 1.51164e-3;
-    const double term2 = int2_fac * (tgamma(3., lambdar_size1) - tgamma(3., lambdar_size2)) * lambdar_neg_2;
-
-    // Integral for 1511.64 micron < D <= 3477.84 micron.
+    // for Vevapi (integral for 1511.64 micron < D < 3477.84 micron)
     const double int3_fac = sqrt(17.32) * pow(d3_to_gram, 1.0 / 12.0);
     const double lambdar_neg_1_75 = pow(lambdar, -1.75);
     const double lambdar_size3 = lambdar * 3.47784e-3;
-    const double term3 = int3_fac * (tgamma(2.75, lambdar_size2) - tgamma(2.75, lambdar_size3)) * lambdar_neg_1_75;
-
-    // Integral for 3477.84 micron < D.
+    // for Vevapi (integral for 3477.84 micron < D)
     const double int4_fac = sqrt(9.17);
     const double lambdar_neg_1_5 = pow(lambdar, -1.5);
-    const double term4 = int4_fac * tgamma(2.5, lambdar_size3) * lambdar_neg_1_5;
 
-    const double v_evap = term1 + term2 + term3 + term4;
+    auto Vevapi = [&](const int i){
+      const double term1 = int1_fac * tgamma_lower(constants.mur+2.5+i, lambdar_size1)
+                                    * lambdar_neg_2_5;
+      const double term2 = int2_fac * (tgamma(constants.mur+2.0+i, lambdar_size1)
+                                     - tgamma(constants.mur+2.0+i, lambdar_size2))
+                                    * lambdar_neg_2;
+      const double term3 = int3_fac * (tgamma(constants.mur+1.75+i, lambdar_size2)
+                                     - tgamma(constants.mur+1.75+i, lambdar_size3))
+                                    * lambdar_neg_1_75;
+      const double term4 = int4_fac * tgamma(constants.mur+1.5+i, lambdar_size3)
+                                    * lambdar_neg_1_5;
+      return term1 + term2 + term3 + term4;
+    };
+
+    const double Vevap1 = Vevapi(1);
 
     if constexpr (WithGrad) {
-      const double term1_dl = (term1 - int1_fac * tgamma_lower(4.5, lambdar_size1) * lambdar_neg_2_5) / lambdar;
-      const double term2_dl = (term2 - int2_fac * (tgamma(4., lambdar_size1) - tgamma(4., lambdar_size2)) * lambdar_neg_2) / lambdar;
-      const double term3_dl = (term3 - int3_fac * (tgamma(3.75, lambdar_size2) - tgamma(3.75, lambdar_size3)) * lambdar_neg_1_75) / lambdar;
-      const double term4_dl = (term4 - int4_fac * tgamma(3.5, lambdar_size3) * lambdar_neg_1_5) / lambdar;
-      return {v_evap, {term1_dl + term2_dl + term3_dl + term4_dl}};
+      const double Vevap2 = Vevapi(2);
+
+      return {Vevap1, {((constants.mur+1)*Vevap1 - Vevap2)}};
     } else {
-      return v_evap;
+      return Vevap1;
     }
   }
 
   // Calculate characteristic velocity (v_evap) using Riemann sum over midpoints
   template <bool WithGrad = false>
   static RealOptGrad<WithGrad, 1> calc_v_evap_numerical(const RainshaftConstants &constants, const double lambdar)
   {
+    if (constants.mur != 0.0)
+      throw std::logic_error("Numerical integration of evaporation speed is not supported for mur != 0");
+
     constexpr double dd = 2.; // micron
     constexpr std::size_t low_k = 1;
     constexpr std::size_t high_k = 10000;

rainshaft/rainshaft_constants.hpp

@@ -24,13 +24,15 @@ struct RainshaftConstants {
   double dr_min;
   // Maximum allowed mean rain diameter (m)
   double dr_max;
+  // Rain shape parameter
+  double mur;
   // rho at top of column (kg/m^3)
   double rho_top;
   // nr at top of column (#/kg)
   double nr_top;
   // qr at top of column (kg/kg)
   double qr_top;
-  // constant factor for determining q_sat_dry regularization in the 
+  // constant factor for determining q_sat_dry regularization in the
   // interval (q_sat_dry, q_sat_dry + epsilon_qsat_dry*q_sat_dry)
   double epsilon_qsat_fac;
   // constant factor determining width of self-collection regularization

rainshaft/rainshaft_derived_vars.cpp

@@ -1,6 +1,5 @@
 #include "rainshaft_derived_vars.hpp"
 #include <cstddef>
-#include <cmath>
 #include <stdexcept>
 #include <numeric>
 
@@ -54,14 +53,13 @@ std::vector<double> calc_lambdar(const RainshaftConstants& constants,
   VarConst nr = state.get_variable("nr");
   VarConst qr = state.get_variable("qr");
   for (std::size_t il = 0; il != grid.nlev; ++il) {
-    double lambda_cubed;
+    double lambdar_value;
     if (regularize) {
-      lambda_cubed = constants.pi * constants.rhow * (nr[il] + constants.qsmall * (1.e8))
-        / (qr[il] + constants.qsmall);
+      lambdar_value = calc_lambdar(constants.pi, constants.rhow, constants.mur, nr[il], qr[il], constants.qsmall);
     } else {
-      lambda_cubed = constants.pi * constants.rhow * nr[il] / qr[il];
+      lambdar_value = calc_lambdar(constants.pi, constants.rhow, constants.mur, nr[il], qr[il]);
     }
-    lambdar[il] = std::cbrt(lambda_cubed);
+    lambdar[il] = lambdar_value;
   }
   return lambdar;
 }

rainshaft/rainshaft_derived_vars.hpp

@@ -2,6 +2,7 @@
 #define RAINSHAFT_DERIVED_VARS_HPP
 #include <vector>
 #include <array>
+#include <cmath>
 
 #include "rainshaft_constants.hpp"
 #include "rainshaft_grid.hpp"
@@ -82,9 +83,27 @@ class RainshaftDerivedVars {
     } else {
       return dz_i;
     }
-  } 
+  }
 };
 
+// standard definition of lambdar:
+//   lambdar^3 := pi rho_w / 6 (mur+1)(mur+2)(mur+3) nr/qr
+constexpr double calc_lambdar(const double pi, const double rhow, const double mur,
+  const double nr, const double qr)
+{
+  const double gamma_ratio = (mur+1)*(mur+2)*(mur+3);
+  return std::cbrt(pi * rhow / 6.0 * gamma_ratio * nr/qr);
+}
+
+// regularized definition of lambdar with regularization parameter qsmall:
+//   lambdar^3 := pi rho_w / 6  (mur+1)(mur+2)(mur+3) (nr+1e8*qsmall)/(qr+qsmall)
+constexpr double calc_lambdar(const double pi, const double rhow, const double mur,
+  const double nr, const double qr, const double qsmall)
+{
+  const double gamma_ratio = (mur+1)*(mur+2)*(mur+3);
+  return std::cbrt(pi * rhow / 6.0 * gamma_ratio * (nr + qsmall*1.e8)/(qr + qsmall));
+}
+
 // Convert cell widths to interface heights.
 std::vector<double> dz_to_z_int(const std::vector<double> dz);
 

rainshaft/rainshaft_driver.cpp

@@ -80,7 +80,7 @@ int main(int argc, char* argv[])
   // Load from command line to check for input file
   po::variables_map vm{};
   po::store(parse_command_line(argc, argv, desc), vm);
-  po::notify(vm);  
+  po::notify(vm);
 
   // Load from settings file if available
   if (vm.count("i")) {
@@ -90,13 +90,13 @@ int main(int argc, char* argv[])
     std::ifstream settings_file("settings_filtered.ini");
     vm = po::variables_map();
     po::store(po::parse_config_file(settings_file, desc), vm);
-    po::notify(vm); 
+    po::notify(vm);
     std::filesystem::remove("settings_filtered.ini");
   } else {
     std::cout << "Setting parameters from command line arguments..." << std::endl;
   }
   std::cout << "---------------------------------------------------" << std::endl;
-  
+
   // Setup dependencies
   option_dependency<std::string>(vm, "case_idx", "ic_file");    // specifying a particular case_idx requires input E3SM data file
 
@@ -125,12 +125,13 @@ int main(int argc, char* argv[])
   if (dt_partition_2 < 0.0) {
     dt_partition_2 = abs(dt_partition_2) * dt;
   }
- 
+
   // Set up model constants.
   // SPS: Choose rho_top in a more principled way?
+  const double mur = 0.0;
   RainshaftConstants constants{3.14159265358979323846,
                                287.04, 1.00464e3, 461.50, 997., 2.501e6,
-                               0.62197, qsmall, 9.80616, 1.e-5, 5.e-3,
+                               0.62197, qsmall, 9.80616, 1.e-5, 5.e-3, mur,
                                0.988919555598356, 1.e3, 1.e-4, epsilon_qsat_fac, epsilon_self_coll};
   // Approximate model top in meters.
   // (The grid maker will actually use the next higher-altitude E3SM level.)
@@ -275,7 +276,7 @@ int main(int argc, char* argv[])
 
     SumProcess all_processes{all_process_vec};
 
-    SizeLimiters size_limiters(constants, 10.e-6, 5.e-3, 0.);
+    SizeLimiters size_limiters(constants, 10.e-6, 5.e-3);
 
     // List of integrators that need to remain allocated for use by original scheme.
     std::vector<std::shared_ptr<RainshaftIntegrator>> backing_integrators;
@@ -294,8 +295,8 @@ int main(int argc, char* argv[])
       } else if (method_type == "forcing") {
         return std::make_unique<ForcingIntegrator>(constants, grid, size_limiters, &partition_1_processes, &partition_2_processes, state_descs, tend_descs, abs_tol, dt, dt_partition_1, dt_partition_2, cfl_substep, postprocess, regularize_lambdar, steps_per_output);
       } else if (method_type == "original") {
-        // Borrowed from original P3 settings, minimum diameter is 10 micron and maximum is 5 millimeter, mu = 0.
-        SizeLimiters size_limiters(constants, 10.e-6, 5.e-3, 0.);
+        // Borrowed from original P3 settings, minimum diameter is 10 micron and maximum is 5 millimeter.
+        SizeLimiters size_limiters(constants, 10.e-6, 5.e-3);
         std::shared_ptr<ExplicitIntegrator> local_intg = std::make_shared<ExplicitIntegrator>(constants, grid, size_limiters, &partition_2_processes, state_descs, tend_descs, abs_tol, dt, 1, rel_tol, false, regularize_lambdar);
         backing_integrators.emplace_back(local_intg);
         std::shared_ptr<LimitingIntegrator> local_lim_intg = std::make_shared<LimitingIntegrator>(constants, size_limiters, *local_intg);
@@ -316,7 +317,7 @@ int main(int argc, char* argv[])
       int error_flag;
       RainshaftSolution solution = intg->integrate(initial_time, final_time, initial_state, error_flag);
       auto after_sol = high_resolution_clock::now();
-  
+
       // Time taken for solution.
       duration<double, std::milli> walltime_ms = after_sol - before_sol;
       std::cout << "Case " << icase << ", Time: " << walltime_ms.count() << std::endl;

rainshaft/sedimentation.cpp

@@ -14,7 +14,7 @@ std::optional<LookupLinear> Sedimentation::create_lookup(const RainshaftConstant
       std::vector{5., 8595., 360000.},
       std::vector{0.001, 1.},
       [=] (const double x) {
-          const auto [v0, v3] = rain_fall_speeds(constants, 1.e6 / x, use_numerical_integration);
+          const auto [v0, v3] = rain_fall_speeds(constants, 1.e6 * (constants.mur+1.0)/x, use_numerical_integration);
           return create_v0 ? v0 : v3;
       });
   } else {
@@ -24,7 +24,7 @@ std::optional<LookupLinear> Sedimentation::create_lookup(const RainshaftConstant
 
 Sedimentation::Sedimentation(const RainshaftConstants &constants, bool use_v_table,
                              bool use_numerical_integration)
-    : use_numerical_integration(use_numerical_integration), 
+    : use_numerical_integration(use_numerical_integration),
       v0_table(create_lookup(constants, use_v_table, true, use_numerical_integration)),
       v3_table(create_lookup(constants, use_v_table, false, use_numerical_integration))
 {}
@@ -34,8 +34,8 @@ double Sedimentation::calc_max_step(const RainshaftConstants &constants,
                                     const RainshaftDerivedVars& dvars,
                                     double cfl) const
 {
-  double lambdar_top = std::cbrt(constants.pi * constants.rhow
-                                  * constants.nr_top / constants.qr_top);
+  double lambdar_top = calc_lambdar(constants.pi, constants.rhow, constants.mur,
+                                    constants.nr_top, constants.qr_top);
   const auto [v0, v3] = rain_fall_speeds(constants, dvars.rho_dry[0],
                                       lambdar_top);
   double spatial_frequency = v3 / dvars.dz[0];
@@ -54,7 +54,7 @@ void Sedimentation::calc_tend(const RainshaftConstants &constants,
                               const RainshaftDerivedVars &dvars,
                               Tendency& tend) const
 {
-  const double lambdar_top = cbrt(constants.pi * constants.rhow * constants.nr_top / constants.qr_top);
+  const double lambdar_top = calc_lambdar(constants.pi, constants.rhow, constants.mur, constants.nr_top, constants.qr_top);
   auto [v0_prev, v3_prev] = rain_fall_speeds(constants, constants.rho_top, lambdar_top);
   double nr_prev = constants.nr_top;
   double qr_prev = constants.qr_top;
@@ -72,7 +72,7 @@ void Sedimentation::calc_tend(const RainshaftConstants &constants,
     const auto [v0, v3] = rain_fall_speeds(constants, dvars.rho_dry[il], dvars.lambdar[il]);
     const double nr_flux = calc_nr_flux(nr[il], dvars.rho_dry[il], v0);
     const double qr_flux = calc_qr_flux(qr[il], dvars.rho_dry[il], v3);
-    
+
     nr_tend[il] += calc_nr_tend(dvars.dz[il], dvars.rho_dry[il], nr_flux, nr_flux_prev);
     qr_tend[il] += calc_qr_tend(dvars.dz[il], dvars.rho_dry[il], qr_flux, qr_flux_prev);
 
@@ -88,7 +88,7 @@ void Sedimentation::calc_tend_jac(const RainshaftConstants &constants,
                                   const RainshaftDerivedVars &dvars,
                                   Matrix jac) const
 {
-  const double lambdar_top = cbrt(constants.pi * constants.rhow * constants.nr_top / constants.qr_top);
+  const double lambdar_top = calc_lambdar(constants.pi, constants.rhow, constants.mur, constants.nr_top, constants.qr_top);
   auto [v0_prev, v3_prev] = rain_fall_speeds<true>(constants, {constants.rho_top, {0., 0.}}, {lambdar_top, {0., 0.}});
   double nr_prev = constants.nr_top;
   double qr_prev = constants.qr_top;