Initial mu configurability

GridKit/CommonMath.hpp

@@ -10,16 +10,22 @@ namespace GridKit
 {
   namespace Math
   {
+    template <typename RealT>
+    inline constexpr RealT DEFAULT_MU = 240.0;
+
     /**
      * @brief Smoothing scale shared by CommonMath primitives
      *
      * Used by @ref sigmoid, @ref ramp, and functions composed from them to set
      * the width of smooth transitions.
      *
+     * @warning Process-wide setting; configure before constructing models or
+     * launching workers.
+     *
      * @tparam RealT - real data type
      */
     template <typename RealT>
-    inline constexpr RealT MU = 240.0;
+    inline RealT MU = DEFAULT_MU<RealT>;
 
     /**
      * @brief Scaled sigmoid activation function

GridKit/Model/PhasorDynamics/Controller/REECB/Reecb.hpp

@@ -138,11 +138,6 @@ namespace GridKit
             ScalarT*       f);
 
       private:
-        /// Hinge width chosen so closed-circle leakage stays below the
-        /// initialization tolerance [p.u. current squared].
-        static constexpr RealT CURRENT_CIRCLE_KNEE =
-            INITIALIZATION_TOLERANCE / Math::MU<RealT>;
-
         struct InitialPoint;
 
         struct InitialCurrentLimit

GridKit/Model/PhasorDynamics/Controller/REECB/ReecbImpl.hpp

@@ -1081,7 +1081,8 @@ namespace GridKit
       Reecb<scalar_type, index_type>::sqrtramp(ValueT x)
       {
         const RealT root_width = ONE<RealT> / Math::MU<RealT>;
-        const RealT knee       = CURRENT_CIRCLE_KNEE;
+        // Keep closed-circle leakage below the initialization tolerance.
+        const RealT knee       = INITIALIZATION_TOLERANCE / Math::MU<RealT>;
 
         const ValueT absolute    = std::abs(x);
         const ValueT normalizer  = absolute + knee;
@@ -1117,12 +1118,13 @@ namespace GridKit
         }
 
         const RealT mu       = Math::MU<RealT>;
+        const RealT knee     = INITIALIZATION_TOLERANCE / mu;
         const RealT scaled_y = HALF<RealT> * mu * y;
         const RealT hinged   = y / mu
                              * (scaled_y + std::hypot(scaled_y, ONE<RealT>));
         const RealT square = hinged
-                             - QUARTER<RealT> * CURRENT_CIRCLE_KNEE
-                                   * (CURRENT_CIRCLE_KNEE / hinged);
+                             - QUARTER<RealT> * knee
+                                   * (knee / hinged);
         return std::max(ZERO<RealT>, square);
       }
 

application/PhasorDynamics/AnalysisUtilities.hpp

@@ -1,16 +1,19 @@
 #pragma once
 
+#include <cmath>
 #include <cstddef>
 #include <filesystem>
 #include <format>
 #include <fstream>
 #include <iostream>
+#include <stdexcept>
 #include <string>
 #include <vector>
 
 #include <magic_enum/magic_enum.hpp>
 #include <nlohmann/json.hpp>
 
+#include <GridKit/CommonMath.hpp>
 #include <GridKit/Model/PhasorDynamics/SystemModelData.hpp>
 #include <GridKit/Solver/Dynamic/Ida.hpp>
 #include <GridKit/Testing/TestHelpers.hpp>
@@ -60,6 +63,8 @@ namespace GridKit
       double                                         rel_tol;
       /// absolute tolerance for the solver
       double                                         abs_tol;
+      /// CommonMath smoothing scale
+      double                                         mu{Math::DEFAULT_MU<double>};
       /// fixed solver time step size, or 0 for adaptive stepping
       double                                         dt_fixed;
       /// maximum number of solver time steps, or 0 for the IDA default
@@ -82,6 +87,18 @@ namespace GridKit
       SystemModelData<>                              model_data;
     };
 
+    /**
+     * @brief Apply study-wide CommonMath configuration
+     *
+     * Must be called before constructing a system model or launching study
+     * workers.
+     */
+    template <typename RealT>
+    inline void configureCommonMath(const StudyData& study)
+    {
+      Math::MU<RealT> = static_cast<RealT>(study.mu);
+    }
+
     using json = ::nlohmann::json;
     using Log  = ::GridKit::Utilities::Logger;
 
@@ -99,8 +116,13 @@ namespace GridKit
       j.at("system_model_file").get_to(c.system_model_file);
       c.dt_monitor = j.value("dt_monitor", 0.0);
       j.at("tmax").get_to(c.tmax);
-      c.rel_tol            = j.value("rel_tol", DEFAULT_SOLVER_REL_TOL);
-      c.abs_tol            = j.value("abs_tol", DEFAULT_SOLVER_ABS_TOL);
+      c.rel_tol = j.value("rel_tol", DEFAULT_SOLVER_REL_TOL);
+      c.abs_tol = j.value("abs_tol", DEFAULT_SOLVER_ABS_TOL);
+      c.mu      = j.value("mu", Math::DEFAULT_MU<double>);
+      if (!std::isfinite(c.mu) || c.mu <= 0.0)
+      {
+        throw std::invalid_argument("\"mu\" must be a positive finite number");
+      }
       c.dt_fixed           = j.value("dt_fixed", 0.0);
       c.max_steps          = j.value("max_steps", std::size_t{0});
       c.consistent_ic_type = AnalysisManager::Sundials::IdaConsistentICType::YA_YDP;

application/PhasorDynamics/ContingencyAnalysis.cpp

@@ -154,6 +154,8 @@ int main(int argc, const char* argv[])
   checkCommandLine(argc, "ContingencyAnalysis");
   auto study_data = parseStudyData(argv[1]);
 
+  configureCommonMath<real_type>(study_data);
+
   const auto start = Clock::now();
 
   auto faults   = study_data.model_data.bus_fault;

application/PhasorDynamics/DynamicSimulation.cpp

@@ -23,6 +23,8 @@ int main(int argc, const char* argv[])
   checkCommandLine(argc, "DynamicSimulation");
   auto study = parseStudyData(argv[1]);
 
+  configureCommonMath<real_type>(study);
+
   // Instantiate system
   SystemModel<scalar_type, index_type> sys(study.model_data);
   sys.allocate();

application/PhasorDynamics/README.md

@@ -9,6 +9,7 @@
   `tmax`               | A floating-point value for max time
   `rel_tol`            | Relative solver tolerance (default: 1.0e-7)
   `abs_tol`            | Absolute solver tolerance override (default: 1.0e-9)
+  `mu`                 | CommonMath smoothing scale; must be positive and finite (default: 240.0)
   `dt_fixed`           | Fixed solver time step size, or 0 for adaptive stepping (default: 0)
   `max_steps`          | Maximum number of solver time steps, 0 for the IDA default, or a negative number for unlimited steps (default: 0)
   `consistent_ic_type` | IDA consistent initial condition calculation type; one of { "y", "ya_ydp" } (default: "ya_ydp")

tests/UnitTests/Math/SmoothnessIndicatorTests.hpp

@@ -225,6 +225,14 @@ namespace GridKit
         success *= std::isfinite(Math::ramp(far_below));
         success *= (Math::ramp(far_below) < roundoff);
 
+        // Verify ramp uses the configured runtime smoothing scale.
+        const RealT original_mu   = Math::MU<RealT>;
+        const RealT configured_mu = 500.0;
+
+        Math::MU<RealT>  = configured_mu;
+        success         *= within(Math::ramp(scalar(0.0)), std::log(scalar(2.0)) / scalar(configured_mu), roundoff);
+        Math::MU<RealT>  = original_mu;
+
         success *= (smooth_clip > lower);
         success *= (smooth_clip < upper);
         success *= std::isfinite(Math::clamp(far_above, lower, upper));

tests/UnitTests/PhasorDynamics/ControllerReecbTests.hpp

@@ -314,7 +314,7 @@ namespace GridKit
         success *= stateMatches(fixture.reecb,
                                 {{Vars::ILCAP, 2.0}},
                                 "initial current-circle capacity",
-                                kCircleTol);
+                                circleTolerance());
 
         const auto* initial_values  = fixture.reecb.y().getData();
         success                    *= scalarMatches(initial_values[index(Vars::IQMAX)],
@@ -381,7 +381,7 @@ namespace GridKit
         success                    *= stateMatches(system_base.reecb,
                                                    {{Vars::ILCAP, 2.0}},
                                 "omitted-rating current-circle capacity",
-                                kCircleTol);
+                                circleTolerance());
         success                    *= allResidualsWithinInitTolerance(system_base.reecb);
 
         return success.report(__func__);
@@ -710,7 +710,7 @@ namespace GridKit
           success *= stateMatches(boundary.reecb,
                                   {{Vars::ILCAP, test_case.capacity}},
                                   test_case.label,
-                                  kCircleTol);
+                                  circleTolerance());
           success *= allResidualsWithinInitTolerance(boundary.reecb);
         }
 
@@ -784,7 +784,7 @@ namespace GridKit
         success *= stateMatches(exhausted.reecb,
                                 {{Vars::ILCAP, 0.0}},
                                 "injection does not expand current circle",
-                                kCircleTol);
+                                circleTolerance());
         success *= allResidualsWithinInitTolerance(exhausted.reecb);
 
         Log::setVerbosity(previous_verbosity);
@@ -846,7 +846,7 @@ namespace GridKit
           RealT tolerance = kTol;
           if (variable == Vars::ILCAP)
           {
-            tolerance = kCircleTol;
+            tolerance = circleTolerance();
           }
 
           success *= variableMatches(residuals[index(variable)],
@@ -914,7 +914,7 @@ namespace GridKit
                   success *= stateMatches(fixture.reecb,
                                           {{Vars::ILCAP, 2.0}},
                                           "selector ILCAP",
-                                          kCircleTol);
+                                          circleTolerance());
 
                   // Exactly one reactive path carries the operating point.
                   const auto* y = fixture.reecb.y().getData();
@@ -1545,7 +1545,7 @@ namespace GridKit
             success *= residualsMatch(fixture.reecb,
                                       {{Vars::ILCAP, circleLeg(square) - capacity_state}},
                                       label,
-                                      kCircleTol);
+                                      circleTolerance());
           }
         }
 
@@ -1596,7 +1596,7 @@ namespace GridKit
           success *= residualsMatch(open.reecb,
                                     {{Vars::ILCAP, open_limit}},
                                     "finite open current circle",
-                                    kCircleTol);
+                                    circleTolerance());
           success *= allResidualsFinite(open.reecb);
         }
 
@@ -1635,7 +1635,7 @@ namespace GridKit
             success *= residualsMatch(fixture.reecb,
                                       {{Vars::ILCAP, ideal_capacity}},
                                       "off-axis capacity",
-                                      kCircleTol);
+                                      circleTolerance());
 
             const RealT capacity = fixture.reecb.getResidual().getData()[index(Vars::ILCAP)];
             if (!std::isfinite(capacity) || capacity < ZERO<RealT>)
@@ -1977,8 +1977,10 @@ namespace GridKit
       }
 
       /// MU-aware tolerance for ideal current-circle comparisons.
-      static constexpr RealT kCircleTol =
-          std::max(ONE<RealT> / (Math::MU<RealT> * Math::MU<RealT>), kTol);
+      static RealT circleTolerance()
+      {
+        return std::max(ONE<RealT> / (Math::MU<RealT> * Math::MU<RealT>), kTol);
+      }
 
       static constexpr size_t kBusVrColumn = index(Vars::MAXIMUM);
 

tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp

@@ -38,8 +38,11 @@ namespace GridKit
 
       static constexpr ScalarT kTol =
           static_cast<ScalarT>(100.0) * std::numeric_limits<ScalarT>::epsilon();
-      static constexpr RealT kSmoothTol =
-          std::max(ONE<RealT> / (Math::MU<RealT> * Math::MU<RealT>), kTol);
+
+      static RealT smoothTolerance()
+      {
+        return std::max(ONE<RealT> / (Math::MU<RealT> * Math::MU<RealT>), kTol);
+      }
 
       /// Construction, the monitor, and every verify() error class: missing
       /// and invalid parameters, a null bus, and an unlinked command port.
@@ -187,8 +190,8 @@ namespace GridKit
         // The latched commands restore both displaced states at their ideal
         // interior first-order rates.
         const auto* f  = latched.regca.getResidual().getData();
-        success       *= scalarMatches(f[index(Vars::IP)], 0.5, "latched active-current rate", kSmoothTol);
-        success       *= scalarMatches(f[index(Vars::IQ)], 0.3, "latched reactive-current rate", kSmoothTol);
+        success       *= scalarMatches(f[index(Vars::IP)], 0.5, "latched active-current rate", smoothTolerance());
+        success       *= scalarMatches(f[index(Vars::IQ)], 0.3, "latched reactive-current rate", smoothTolerance());
 
         return success.report(__func__);
       }
@@ -378,7 +381,7 @@ namespace GridKit
         const auto* f = residual.getData();
         for (const auto& row : expected)
         {
-          success *= scalarMatches(f[index(row.variable)], row.value, row.name, kSmoothTol);
+          success *= scalarMatches(f[index(row.variable)], row.value, row.name, smoothTolerance());
         }
 
         return success.report(__func__);
@@ -577,7 +580,7 @@ namespace GridKit
 
         // Initialization above Vhvmax activates HVRCM and preserves Q0.
         {
-          const RealT terminal_voltage = kHvrcmVoltageLimit + kHvrcmOffset;
+          const RealT terminal_voltage = kHvrcmVoltageLimit + hvrcmOffset();
 
           auto data                   = makeData();
           data.parameters[Params::q0] = 0.1;
@@ -590,7 +593,7 @@ namespace GridKit
           const auto* y              = fixture.regca.y().getData();
           const RealT extra_current  = y[index(Vars::IQEXTRA)];
           success                   *= scalarMatches(
-              extra_current, kHvrcmGain * kHvrcmOffset, "IQEXTRA with the default Khv", kSmoothTol);
+              extra_current, kHvrcmGain * hvrcmOffset(), "IQEXTRA with the default Khv", smoothTolerance());
           success *= scalarMatches(
               y[index(Vars::IQ)] - y[index(Vars::IQEXTRA)], 0.1 / terminal_voltage, "IQ preserves Q0 after HVRCM compensation");
           success *= scalarMatches(y[index(Vars::QBR)], 0.1, "QBR");
@@ -630,12 +633,12 @@ namespace GridKit
             return fixture.regca.getResidual().getData()[index(Vars::IQEXTRA)];
           };
 
-          const RealT below = residualAt(kHvrcmVoltageLimit - kHvrcmOffset, 0.0);
+          const RealT below = residualAt(kHvrcmVoltageLimit - hvrcmOffset(), 0.0);
           const RealT at    = residualAt(kHvrcmVoltageLimit, 0.0);
-          const RealT above = residualAt(kHvrcmVoltageLimit + kHvrcmOffset, 0.0);
+          const RealT above = residualAt(kHvrcmVoltageLimit + hvrcmOffset(), 0.0);
 
           success *= scalarMatches(at, kHvrcmGain * hvrcmTransition(), "HVRCM residual at the threshold");
-          success *= scalarMatches(above - below, kHvrcmGain * kHvrcmOffset, "HVRCM symmetric residual difference");
+          success *= scalarMatches(above - below, kHvrcmGain * hvrcmOffset(), "HVRCM symmetric residual difference");
 
           const RealT shifted  = residualAt(kHvrcmVoltageLimit, 0.1);
           success             *= scalarMatches(shifted - at, -0.1, "HVRCM extra-current residual shift");
@@ -832,11 +835,13 @@ namespace GridKit
       static constexpr RealT kHvrcmGain = 0.7;
 
       /// Offset whose smooth-ramp tail decays by one significand width.
-      static constexpr RealT kHvrcmOffset =
-          static_cast<RealT>(std::numeric_limits<RealT>::digits)
-          * std::numbers::ln2_v<RealT> / Math::MU<RealT>;
+      static RealT hvrcmOffset()
+      {
+        return static_cast<RealT>(std::numeric_limits<RealT>::digits)
+               * std::numbers::ln2_v<RealT> / Math::MU<RealT>;
+      }
 
-      static constexpr RealT hvrcmTransition()
+      static RealT hvrcmTransition()
       {
         return std::numbers::ln2_v<RealT> / Math::MU<RealT>;
       }