Allow multiple runs for each simulation and some more refactoring

psim/include/psim/geometry.h

@@ -12,7 +12,6 @@
 
 namespace Geometry {
 
-inline constexpr double GEOEPS = std::numeric_limits<double>::epsilon() * 1E9;
 struct Point;
 using PointPair = std::pair<Point, Point>;
 
@@ -90,18 +89,6 @@ struct Triangle {
 };
 std::ostream& operator<<(std::ostream& os, const Triangle& triangle);// NOLINT
 
-/**
- * @brief Check if two floating-point values are approximately equal within a specified tolerance.
- *
- * @param a The first value to compare.
- * @param b The second value to compare.
- * @param epsilon The tolerance used for the comparison (default: EPSILON).
- * @return true if the values are approximately equal within the given tolerance, false otherwise.
- */
-[[nodiscard]] inline bool approxEqual(double a, double b, double epsilon = GEOEPS) noexcept {// NOLINT
-    return std::abs(a - b) < epsilon;
-}
-
 class ShapeError : public std::exception {
 public:
     [[nodiscard]] const char* what() const noexcept override {

psim/include/psim/model.h

@@ -9,6 +9,17 @@
 #include "sensorInterpreter.h"
 #include <unordered_map>
 
+struct ModelParams {
+    std::size_t num_runs{ 1 };
+    std::size_t num_cells;
+    std::size_t num_sensors;
+    std::size_t measurement_steps;
+    std::size_t num_phonons;
+    double simulation_time;
+    double t_eq;
+    bool phasor_sim{ false };
+};
+
 /**
  * The Model class primarily controls the geometrical aspects of the simulation.
  * Materials must be added first with the addMaterials function. Then the sensors must be added as each
@@ -19,7 +30,6 @@
 class Model {
 public:
     using Point = Geometry::Point;
-    using SimulationType = Sensor::SimulationType;
 
     /**
      * @param num_cells - The number of cells that the model will contain.
@@ -31,13 +41,7 @@ class Model {
      * @param t_eq - The linearization (equilibrium) temperature of the system - 0 indicates a full simulation
      * @param phasor_sim - True -> phonons have uniform direction & velocity & no scattering
      */
-    Model(std::size_t num_cells,
-        std::size_t num_sensors,
-        std::size_t measurement_steps,
-        std::size_t num_phonons,
-        double simulation_time,
-        double t_eq,
-        bool phasor_sim = false);
+    Model(const ModelParams& params);
     /**
      * @param type - The type of simulation. Default is steady state.
      * @param step_interval - For periodic and transient simulations only. The distance between steps for which
@@ -61,10 +65,8 @@ class Model {
      * @param t_init - The initial temperature of the cells linked to the sensor
      * @param type - The type of simulation - steady-state periodic or transient
      */
-    void addSensor(std::size_t ID,
-        const std::string& material_name,
-        double t_init,
-        Sensor::SimulationType type);// NOLINT
+    void addSensor(std::size_t ID, const std::string& material_name, double t_init,
+        SimulationType type);// NOLINT
     /**
      * Adds a triangular cell to the system.
      * Ensures new cells are not contained in existing cells and existing cells do not contain the new cell.
@@ -109,6 +111,7 @@ class Model {
 
 private:
     SimulationType sim_type_{ SimulationType::SteadyState };
+    std::size_t num_runs_;
     std::size_t num_cells_;
     std::size_t measurement_steps_;
     double simulation_time_;
@@ -117,6 +120,7 @@ class Model {
     bool phasor_sim_;
     std::size_t start_step_{ 0 };// 0 for SS simulations -> measurements vectors are scaled down in the sensors
 
+    // TODO: Not  a big fan of this, prefer dependency injection
     ModelSimulator simulator_;
     OutputManager outputManager_;
     SensorInterpreter interpreter_;
@@ -140,13 +144,13 @@ class Model {
     [[nodiscard]] std::pair<double, double> setTemperatureBounds() noexcept;
     void initializeMaterialTables(double low_temp, double high_temp);
     [[nodiscard]] double avgTemp() const;
-    void storeResults() noexcept;
+    void storeResults(std::size_t runId) noexcept;
     /**
      * @return - The new t_eq if the system needs to be reset and re-run - less than 90% of the model's sensor
      * temperatures are stable (90%) and std::nullopt otherwise
      */
-    [[nodiscard]] std::optional<double> resetRequired() noexcept;
-    void reset() noexcept;
+    [[nodiscard]] std::optional<double> resetRequired() const noexcept;
+    void reset(bool full_reset = false) noexcept;
 };
 
 

psim/include/psim/outputManager.h

@@ -1,7 +1,9 @@
 #ifndef PSIM_OUTPUTMANAGER_H
 #define PSIM_OUTPUTMANAGER_H
 
+#include "utils.h"
 #include <filesystem>
+#include <unordered_map>
 #include <vector>
 
 struct SensorMeasurements;
@@ -10,27 +12,31 @@ namespace fs = std::filesystem;
 class OutputManager {
 public:
     OutputManager() = default;
-    void steadyStateExport(const fs::path& filepath, double time) const;
     /**
      * Exports results from each measurement step so the evolution of the system can be visualized.
      * Step intervals of 1 will write every measurement to file. If there are 100 measurements,
-     * then there will bee 100 * num_sensors entries written. Step intervals greater than 1 will write the average of
+     * then there will be 100 * num_sensors entries written. Step intervals greater than 1 will write the average of
      * the measurements to file. Step intervals of 10 with 100 measurements -> 10 * num_sensors entries. First entry is
      * the avg of steps 1-10, second is avg of steps 11-20, etc.
      * @param filepath - path and filename where the results will be written - existing file will be overwritten
-     * @param time- The time taken to run the simulation.
+     * @param time - The time taken to run the simulation.
+     * @param num_runs - The number of runs the simulation was executed for.
+     * @param type - The type of simulation that was run.
      */
-    void periodicExport(const fs::path& filepath, double time) const;
-    void addMeasurement(SensorMeasurements&& measurement) noexcept;
-    void sortMeasurements() noexcept;
+    void exportResults(const fs::path& filepath, double time, std::size_t num_runs, SimulationType type) const;
+    void addMeasurement(std::size_t runId, SensorMeasurements&& measurement) noexcept;
+    void sortMeasurements(std::size_t runId) noexcept;
     void setStepInterval(std::size_t interval) noexcept {
         step_interval_ = interval;
     }
 
 private:
     std::size_t step_interval_{ 1 };
-    std::vector<SensorMeasurements> measurements_;
+    std::unordered_map<std::size_t, std::vector<SensorMeasurements>> measurements_;
 
+    std::vector<SensorMeasurements> calculateAverages() const;
+    void steadyStateExport(std::ofstream& output) const;
+    void periodicExport(std::ofstream& output) const;
     [[nodiscard]] static std::filesystem::path adjustPath(const fs::path& filepath, const std::string& prepend);
     [[nodiscard]] static std::string getCurrentDateTime();
 };

psim/include/psim/sensor.h

@@ -3,14 +3,14 @@
 
 #include "material.h"
 #include "sensorController.h"
+#include "utils.h"
 #include <memory>
 #include <mutex>
 
 class Phonon;
 
 class Sensor {
 public:
-    enum class SimulationType { SteadyState, Periodic, Transient };
     Sensor(std::size_t ID,// NOLINT
         const Material& material,
         SimulationType type,
@@ -60,7 +60,7 @@ class Sensor {
      * @param step - The measurement step to update
      */
     void updateHeatParams(const Phonon& p, std::size_t step) noexcept;// NOLINT
-    void reset() noexcept;
+    void reset(bool full_reset) noexcept;
     void updateTables() const {
         controller_->updateTables();
     }

psim/include/psim/sensorController.h

@@ -34,7 +34,7 @@ class SensorController {
      * @return - true if the temp of this sensor is unstable (final temp not within some percentage of initial temp)
      */
     [[nodiscard]] virtual bool resetRequired(double t_final, std::vector<double>&&) noexcept;
-    virtual void reset() noexcept = 0;
+    virtual void reset(bool full_reset = false) noexcept = 0;
 
 protected:
     const Material& material_;
@@ -69,7 +69,7 @@ class SteadyStateController : public SensorController {
         return t_steady_;
     }
 
-    void reset() noexcept override;
+    void reset(bool full_reset = false) noexcept override;
 };
 
 class PeriodicController : public SensorController {
@@ -87,7 +87,7 @@ class PeriodicController : public SensorController {
         return t_steady_;
     }
 
-    void reset() noexcept override;
+    void reset(bool full_reset = false) noexcept override;
 };
 
 // May want to average results of final 10% of runs or something like this
@@ -107,7 +107,7 @@ class TransientController : public SensorController {
     void scatterUpdate(Phonon& p) const noexcept override;// NOLINT
     [[nodiscard]] bool resetRequired([[maybe_unused]] double t_final,
         std::vector<double>&& final_temps) noexcept override;
-    void reset() noexcept override;
+    void reset(bool full_reset) noexcept override;
 };
 
 #endif// PSIM_SENSORCONTROLLER_H

psim/include/psim/utils.h

@@ -1,18 +1,22 @@
 #ifndef PSIM_UTILS_H
 #define PSIM_UTILS_H
 
+#include "geometry.h"
+
 #include <array>
 #include <random>
 
-namespace Utils {
-
+enum class SimulationType { SteadyState, Periodic, Transient };
+inline constexpr double GEOEPS = std::numeric_limits<double>::epsilon() * 1E9;
 inline constexpr double PI = 3.1415926535897932384626433832795028841971693993751058209749445923;
 
+namespace Utils {
+
 // Generates a random number from a uniform distribution over [0,1].
 inline double urand() noexcept {
     static std::random_device rd;// NOLINT
     thread_local std::mt19937 generator(rd());
-    std::uniform_real_distribution dist(0., std::nextafter(1.0, 2.0));// NOLINT
+    std::uniform_real_distribution dist(0., std::nextafter(1., 2.));// NOLINT
     return dist(generator);
 }
 
@@ -51,6 +55,19 @@ std::array<std::pair<T, U>, V> zip(const std::array<T, V>& r1, const std::array<
     while (it1 != std::cend(r1) && it2 != std::cend(r2)) { result[index++] = { *it1++, *it2++ }; }
     return result;
 }
+
+/**
+ * @brief Check if two floating-point values are approximately equal within a specified tolerance.
+ *
+ * @param a The first value to compare.
+ * @param b The second value to compare.
+ * @param epsilon The tolerance used for the comparison (default: EPSILON).
+ * @return true if the values are approximately equal within the given tolerance, false otherwise.
+ */
+[[nodiscard]] inline bool approxEqual(double a, double b, double epsilon = GEOEPS) noexcept {// NOLINT
+    return std::abs(a - b) < epsilon;
+}
+
 }// namespace Utils
 
 #endif// PSIM_UTILS_H

psim/src/geometry.cpp

@@ -1,4 +1,7 @@
 #include "psim/geometry.h"
+
+#include "psim/utils.h"
+
 #include <algorithm>// for minmax, any_of
 #include <array>// for array, array<>::value_type
 #include <cmath>// for fabs, sqrt
@@ -320,7 +323,7 @@ double getLineLength(const Line& line) noexcept {
 
 // Returns 0 for a vertical line!
 double getSlope(const Point& p1, const Point& p2) noexcept {// NOLINT
-    return (approxEqual(p1.x,p2.x)) ? 0. : (p1.y - p2.y) / (p1.x - p2.x);
+    return (Utils::approxEqual(p1.x,p2.x)) ? 0. : (p1.y - p2.y) / (p1.x - p2.x);
 }
 
 double getIntercept(double x, double y, double slope) noexcept {// NOLINT

psim/src/inputManager.cpp

@@ -10,20 +10,23 @@
 using json = nlohmann::json;
 using Point = Geometry::Point;
 using Triangle = Geometry::Triangle;
-using SimulationType = Sensor::SimulationType;
 
 std::optional<Model> InputManager::deserialize(const std::filesystem::path& filepath) {// NOLINT
 
     // ************************ Helper methods **********************************
     auto buildModel = [](std::size_t num_cells, std::size_t num_sensors, const auto& s_data) {
         const bool phasor_sim = s_data.at("phasor_sim").dump() == "true";
-        return Model(num_cells,
+        ModelParams params{ 1,
+            num_cells,
             num_sensors,
             s_data.at("num_measurements"),
             s_data.at("num_phonons"),
             s_data.at("sim_time"),
             s_data.at("t_eq"),
-            phasor_sim);
+            phasor_sim };
+        if (s_data.contains("num_runs")) { params.num_runs = s_data.at("num_runs"); }
+
+        return Model(params);
     };
 
     auto addMaterial = [](auto& model, double t_eq, const auto& m_data, std::size_t id) {// NOLINT

psim/src/material.cpp

@@ -38,13 +38,13 @@ Material::Material(std::size_t mat_id, const DispersionData& disp_data, const Re
     for (const auto& freq : frequencies_) {
         const double la_gv = getGv(freq, LA_coeffs);
         velocities_la_[index] = la_gv;
-        densities_la_[index] = pow(getK(freq, LA_coeffs), 2) / 2. / pow(Utils::PI, 2) / la_gv;// NOLINT
+        densities_la_[index] = pow(getK(freq, LA_coeffs), 2) / 2. / pow(PI, 2) / la_gv;// NOLINT
         const double ta_gv = getGv(freq, TA_coeffs);
         if (!std::isnan(ta_gv)) {
             velocities_ta_[index] = ta_gv;
             // Do not divide by 2 since the TA branch is doubly degenerate
             // Can account for this in the density of states array to save computations elsewhere
-            densities_ta_[index] = pow(getK(freq, TA_coeffs), 2) / pow(Utils::PI, 2) / ta_gv;
+            densities_ta_[index] = pow(getK(freq, TA_coeffs), 2) / pow(PI, 2) / ta_gv;
         }
         ++index;
     }