Phase 3.1 started: added thermal_mass solver strategies (Euler/RK4) with method-dependent stability, updated interfaces/docs, and introduced regression and determinism test coverage

Author: CameronBrooks11

docs/api-reference.md

@@ -241,6 +241,7 @@ model.params["ambient_signal"] = std::string("ambient.temp");
 model.params["thermal_mass"] = 1000.0;  // J/degC
 model.params["heat_transfer_coeff"] = 10.0;  // W/degC
 model.params["initial_temp"] = 25.0;  // degC
+model.params["integration_method"] = std::string("forward_euler"); // optional: "forward_euler" (default) or "rk4"
 spec.models.push_back(model);
 ```
 
@@ -504,8 +505,12 @@ Parameters:
 - h: heat_transfer_coeff (W/degC) - Convection coefficient
 - P: power_signal (W) - Heat input
 - T_ambient: ambient_signal (degC) - Ambient temperature
+- integration_method (optional): `"forward_euler"` (default) or `"rk4"`
 
-**Stability limit:** dt < 2\*C/h
+**Stability limit:**
+
+- `forward_euler`: dt < 2\*C/h
+- `rk4`: dt < 2.785293563\*C/h (negative real-axis stability bound)
 
 Example usage:
 
@@ -518,6 +523,7 @@ model.params["ambient_signal"] = std::string("ambient");
 model.params["thermal_mass"] = 1000.0;
 model.params["heat_transfer_coeff"] = 10.0;
 model.params["initial_temp"] = 25.0;
+model.params["integration_method"] = std::string("rk4"); // optional
 ```
 
 ---

docs/numerical-methods.md

@@ -0,0 +1,57 @@
+# Numerical Methods Policy
+
+## Scope
+
+This document defines the numerical integration policy for FluxGraph models.
+It complements `docs/semantics_spec.md` and is authoritative for solver
+selection and stability interpretation.
+
+## Current Policy (Phase 3.1 baseline)
+
+1. Integration method selection is explicit per model via model parameters.
+2. If not specified, the deterministic default is `forward_euler`.
+3. Methods are selected at compile time and remain fixed at runtime.
+4. Runtime behavior must be deterministic for identical inputs and `dt`.
+
+## ThermalMassModel Integration Methods
+
+`thermal_mass` supports:
+
+1. `forward_euler` (default)
+2. `rk4` (classic fourth-order Runge-Kutta)
+
+Selection is provided via:
+
+```cpp
+model.params["integration_method"] = std::string("forward_euler");
+// or
+model.params["integration_method"] = std::string("rk4");
+```
+
+Unknown method names are compile-time errors.
+
+## Stability Limits
+
+For the linear thermal model `dT/dt = (P - h*(T - T_amb)) / C`, with
+`k = h/C`:
+
+1. `forward_euler`: `dt < 2/k = 2*C/h`
+2. `rk4`: `dt < 2.785293563/k = 2.785293563*C/h`
+
+If `h <= 0`, the model is treated as unconditionally stable for this criterion.
+
+## Validation Expectations
+
+Phase 3 validation must report:
+
+1. Error metrics (`L2`, `Linf`) versus analytical references.
+2. Convergence behavior as `dt` is refined.
+3. Determinism checks for each supported integration method.
+
+## Forward Compatibility
+
+Future methods (for example trapezoidal or implicit schemes) must define:
+
+1. Stability policy and limits.
+2. Deterministic selection and defaults.
+3. Regression and analytical validation coverage before release.

include/fluxgraph/model/interface.hpp

@@ -21,7 +21,7 @@ class IModel {
   virtual void reset() = 0;
 
   /// Compute maximum stable time step for this model
-  /// Based on Forward Euler stability criteria
+  /// Based on the model's configured integration method
   /// @return Maximum safe dt in seconds
   virtual double compute_stability_limit() const = 0;
 

include/fluxgraph/model/thermal_mass.hpp

@@ -7,6 +7,15 @@
 
 namespace fluxgraph {
 
+enum class ThermalIntegrationMethod {
+  ForwardEuler,
+  Rk4,
+};
+
+const char *to_string(ThermalIntegrationMethod method);
+ThermalIntegrationMethod
+parse_thermal_integration_method(const std::string &method_name);
+
 /// Thermal mass model: simple heat capacity with power input and ambient
 /// cooling Physics: dT/dt = (P_in - h*(T - T_amb)) / C Where:
 ///   T = temperature (degC)
@@ -32,12 +41,14 @@ class ThermalMassModel : public IModel {
                    double heat_transfer_coeff, double initial_temp,
                    const std::string &temp_signal_path,
                    const std::string &power_signal_path,
-                   const std::string &ambient_signal_path, SignalNamespace &ns);
+                   const std::string &ambient_signal_path, SignalNamespace &ns,
+                   ThermalIntegrationMethod integration_method =
+                       ThermalIntegrationMethod::ForwardEuler);
 
   void tick(double dt, SignalStore &store) override;
   void reset() override;
 
-  /// Stability limit for Forward Euler: dt < 2*C/h
+  /// Stability limit for the selected integration method
   double compute_stability_limit() const override;
 
   std::string describe() const override;
@@ -52,6 +63,10 @@ class ThermalMassModel : public IModel {
   double heat_transfer_coeff_; // h (W/K)
   double temperature_;         // Current temp (degC)
   double initial_temp_;        // Initial temp for reset (degC)
+  ThermalIntegrationMethod integration_method_;
+
+  double derivative(double temperature, double net_power,
+                    double ambient) const;
 };
 
 } // namespace fluxgraph

src/graph/compiler.cpp

@@ -293,10 +293,23 @@ void ensure_default_factories_registered_locked(FactoryRegistry &registry) {
         std::string ambient_path =
             as_string(require_param(spec.params, "ambient_signal", context),
                       context + "/ambient_signal");
+        ThermalIntegrationMethod integration_method =
+            ThermalIntegrationMethod::ForwardEuler;
+        if (auto it = spec.params.find("integration_method");
+            it != spec.params.end()) {
+          const std::string method_name =
+              as_string(it->second, context + "/integration_method");
+          try {
+            integration_method = parse_thermal_integration_method(method_name);
+          } catch (const std::invalid_argument &e) {
+            throw std::runtime_error("Invalid parameter at " + context +
+                                     "/integration_method: " + e.what());
+          }
+        }
 
         return std::make_unique<ThermalMassModel>(
             spec.id, thermal_mass, heat_transfer_coeff, initial_temp, temp_path,
-            power_path, ambient_path, ns);
+            power_path, ambient_path, ns, integration_method);
       });
 
   registry.defaults_registered = true;

src/model/thermal_mass.cpp

@@ -1,33 +1,77 @@
 #include "fluxgraph/model/thermal_mass.hpp"
 #include <limits>
+#include <stdexcept>
 #include <sstream>
 
 namespace fluxgraph {
 
+namespace {
+
+constexpr double kRk4NegativeRealAxisStabilityLimit = 2.785293563405282;
+
+} // namespace
+
+const char *to_string(ThermalIntegrationMethod method) {
+  switch (method) {
+  case ThermalIntegrationMethod::ForwardEuler:
+    return "forward_euler";
+  case ThermalIntegrationMethod::Rk4:
+    return "rk4";
+  }
+  return "forward_euler";
+}
+
+ThermalIntegrationMethod
+parse_thermal_integration_method(const std::string &method_name) {
+  if (method_name == "forward_euler") {
+    return ThermalIntegrationMethod::ForwardEuler;
+  }
+  if (method_name == "rk4") {
+    return ThermalIntegrationMethod::Rk4;
+  }
+
+  throw std::invalid_argument(
+      "Unknown thermal integration method '" + method_name +
+      "' (expected one of: forward_euler, rk4)");
+}
+
 ThermalMassModel::ThermalMassModel(const std::string &id, double thermal_mass,
                                    double heat_transfer_coeff,
                                    double initial_temp,
                                    const std::string &temp_signal_path,
                                    const std::string &power_signal_path,
                                    const std::string &ambient_signal_path,
-                                   SignalNamespace &ns)
+                                   SignalNamespace &ns,
+                                   ThermalIntegrationMethod integration_method)
     : id_(id), temp_signal_(ns.intern(temp_signal_path)),
       power_signal_(ns.intern(power_signal_path)),
       ambient_signal_(ns.intern(ambient_signal_path)),
       thermal_mass_(thermal_mass), heat_transfer_coeff_(heat_transfer_coeff),
-      temperature_(initial_temp), initial_temp_(initial_temp) {}
+      temperature_(initial_temp), initial_temp_(initial_temp),
+      integration_method_(integration_method) {}
+
+double ThermalMassModel::derivative(double temperature, double net_power,
+                                    double ambient) const {
+  const double heat_loss = heat_transfer_coeff_ * (temperature - ambient);
+  return (net_power - heat_loss) / thermal_mass_;
+}
 
 void ThermalMassModel::tick(double dt, SignalStore &store) {
   // Read inputs
   double net_power = store.read_value(power_signal_);
   double ambient = store.read_value(ambient_signal_);
 
-  // Compute heat loss to ambient
-  double heat_loss = heat_transfer_coeff_ * (temperature_ - ambient);
-
-  // Forward Euler integration: T += dT/dt * dt
-  double dT = (net_power - heat_loss) / thermal_mass_ * dt;
-  temperature_ += dT;
+  if (integration_method_ == ThermalIntegrationMethod::ForwardEuler) {
+    // Forward Euler integration: T += dT/dt * dt
+    temperature_ += derivative(temperature_, net_power, ambient) * dt;
+  } else {
+    // Classic RK4 with fixed inputs over the tick interval.
+    const double k1 = derivative(temperature_, net_power, ambient);
+    const double k2 = derivative(temperature_ + 0.5 * dt * k1, net_power, ambient);
+    const double k3 = derivative(temperature_ + 0.5 * dt * k2, net_power, ambient);
+    const double k4 = derivative(temperature_ + dt * k3, net_power, ambient);
+    temperature_ += (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4);
+  }
 
   // Write output with unit
   store.write(temp_signal_, temperature_, "degC");
@@ -37,21 +81,26 @@ void ThermalMassModel::tick(double dt, SignalStore &store) {
 void ThermalMassModel::reset() { temperature_ = initial_temp_; }
 
 double ThermalMassModel::compute_stability_limit() const {
-  // Forward Euler stability for dT/dt = -k*T: dt < 2/k
-  // For this model: k = h/C
-  // Therefore: dt < 2*C/h
+  // Stability for dT/dt = -k*T with k = h/C.
   if (heat_transfer_coeff_ <= 0.0) {
     return std::numeric_limits<double>::infinity(); // No cooling =
                                                     // unconditionally stable
   }
-  return 2.0 * thermal_mass_ / heat_transfer_coeff_;
+
+  const double tau = thermal_mass_ / heat_transfer_coeff_;
+  if (integration_method_ == ThermalIntegrationMethod::ForwardEuler) {
+    return 2.0 * tau;
+  }
+
+  return kRk4NegativeRealAxisStabilityLimit * tau;
 }
 
 std::string ThermalMassModel::describe() const {
   std::ostringstream oss;
   oss << "ThermalMass(id=" << id_ << ", C=" << thermal_mass_ << " J/K"
       << ", h=" << heat_transfer_coeff_ << " W/K"
-      << ", T0=" << initial_temp_ << " degC)";
+      << ", T0=" << initial_temp_ << " degC"
+      << ", method=" << to_string(integration_method_) << ")";
   return oss.str();
 }
 

tests/integration/determinism_test.cpp

@@ -139,6 +139,66 @@ TEST(DeterminismTest, SameInputSameOutput) {
   }
 }
 
+TEST(DeterminismTest, SameInputSameOutputRk4) {
+  auto build_graph = []() {
+    GraphSpec spec;
+
+    ModelSpec model;
+    model.id = "thermal";
+    model.type = "thermal_mass";
+    model.params["temp_signal"] = std::string("chamber.temp");
+    model.params["power_signal"] = std::string("chamber.power");
+    model.params["ambient_signal"] = std::string("ambient");
+    model.params["thermal_mass"] = 1000.0;
+    model.params["heat_transfer_coeff"] = 10.0;
+    model.params["initial_temp"] = 25.0;
+    model.params["integration_method"] = std::string("rk4");
+    spec.models.push_back(model);
+
+    return spec;
+  };
+
+  SignalNamespace ns1;
+  FunctionNamespace fn1;
+  SignalStore store1;
+  Engine engine1;
+
+  GraphCompiler compiler1;
+  auto program1 = compiler1.compile(build_graph(), ns1, fn1);
+  engine1.load(std::move(program1));
+
+  SignalNamespace ns2;
+  FunctionNamespace fn2;
+  SignalStore store2;
+  Engine engine2;
+
+  GraphCompiler compiler2;
+  auto program2 = compiler2.compile(build_graph(), ns2, fn2);
+  engine2.load(std::move(program2));
+
+  auto power_id1 = ns1.resolve("chamber.power");
+  auto ambient_id1 = ns1.resolve("ambient");
+  auto temp_id1 = ns1.resolve("chamber.temp");
+  auto power_id2 = ns2.resolve("chamber.power");
+  auto ambient_id2 = ns2.resolve("ambient");
+  auto temp_id2 = ns2.resolve("chamber.temp");
+
+  store1.write(ambient_id1, 20.0, "degC");
+  store2.write(ambient_id2, 20.0, "degC");
+
+  for (int i = 0; i < 500; ++i) {
+    const double power = (i % 200 < 100) ? 500.0 : 100.0;
+    store1.write(power_id1, power, "W");
+    store2.write(power_id2, power, "W");
+
+    engine1.tick(0.1, store1);
+    engine2.tick(0.1, store2);
+
+    EXPECT_DOUBLE_EQ(store1.read_value(temp_id1), store2.read_value(temp_id2))
+        << "Mismatch at tick " << i;
+  }
+}
+
 TEST(DeterminismTest, NoDriftOver10kTicks) {
   // Verify no floating-point drift over long simulations
 

tests/unit/compiler_test.cpp

@@ -269,6 +269,43 @@ TEST(GraphCompilerTest, ParseThermalMassModel) {
   delete model;
 }
 
+TEST(GraphCompilerTest, ParseThermalMassModelWithRk4) {
+  ModelSpec spec;
+  spec.id = "chamber_air";
+  spec.type = "thermal_mass";
+  spec.params["thermal_mass"] = 1000.0;
+  spec.params["heat_transfer_coeff"] = 10.0;
+  spec.params["initial_temp"] = 25.0;
+  spec.params["temp_signal"] = std::string("chamber_air/temperature");
+  spec.params["power_signal"] = std::string("chamber_air/power");
+  spec.params["ambient_signal"] = std::string("chamber_air/ambient");
+  spec.params["integration_method"] = std::string("rk4");
+
+  SignalNamespace ns;
+  GraphCompiler compiler;
+  std::unique_ptr<IModel> model(compiler.parse_model(spec, ns));
+
+  ASSERT_NE(model, nullptr);
+  EXPECT_NE(model->describe().find("method=rk4"), std::string::npos);
+}
+
+TEST(GraphCompilerTest, ParseThermalMassModelWithInvalidIntegrationMethodThrows) {
+  ModelSpec spec;
+  spec.id = "chamber_air";
+  spec.type = "thermal_mass";
+  spec.params["thermal_mass"] = 1000.0;
+  spec.params["heat_transfer_coeff"] = 10.0;
+  spec.params["initial_temp"] = 25.0;
+  spec.params["temp_signal"] = std::string("chamber_air/temperature");
+  spec.params["power_signal"] = std::string("chamber_air/power");
+  spec.params["ambient_signal"] = std::string("chamber_air/ambient");
+  spec.params["integration_method"] = std::string("invalid_method");
+
+  SignalNamespace ns;
+  GraphCompiler compiler;
+  EXPECT_THROW(compiler.parse_model(spec, ns), std::runtime_error);
+}
+
 TEST(GraphCompilerTest, RegisterModelFactoryRejectsInvalidInputs) {
   GraphCompiler::ModelFactory empty_factory;