Use math solver dispatch to separate math solver from main model

power_grid_model_c/power_grid_model/include/power_grid_model/main_core/math_state.hpp

@@ -4,15 +4,16 @@
 
 #pragma once
 
+#include "../math_solver/math_solver_dispatch.hpp"
 #include "../math_solver/y_bus.hpp"
 
 namespace power_grid_model::main_core {
 
 struct MathState {
     std::vector<YBus<symmetric_t>> y_bus_vec_sym;
     std::vector<YBus<asymmetric_t>> y_bus_vec_asym;
-    std::vector<MathSolver<symmetric_t>> math_solvers_sym;
-    std::vector<MathSolver<asymmetric_t>> math_solvers_asym;
+    std::vector<MathSolverProxy<symmetric_t>> math_solvers_sym;
+    std::vector<MathSolverProxy<asymmetric_t>> math_solvers_asym;
 };
 
 inline void clear(MathState& math_state) {

power_grid_model_c/power_grid_model/include/power_grid_model/main_model.hpp

@@ -22,10 +22,12 @@ class MainModel {
   public:
     using Options = MainModelOptions;
 
-    explicit MainModel(double system_frequency, ConstDataset const& input_data, Idx pos = 0)
-        : impl_{std::make_unique<Impl>(system_frequency, input_data, pos)} {}
-    explicit MainModel(double system_frequency, meta_data::MetaData const& meta_data)
-        : impl_{std::make_unique<Impl>(system_frequency, meta_data)} {};
+    explicit MainModel(double system_frequency, ConstDataset const& input_data,
+                       MathSolverDispatcher const& math_solver_dispatcher, Idx pos = 0)
+        : impl_{std::make_unique<Impl>(system_frequency, input_data, math_solver_dispatcher, pos)} {}
+    explicit MainModel(double system_frequency, meta_data::MetaData const& meta_data,
+                       MathSolverDispatcher const& math_solver_dispatcher)
+        : impl_{std::make_unique<Impl>(system_frequency, meta_data, math_solver_dispatcher)} {};
 
     // deep copy
     MainModel(MainModel const& other) {

power_grid_model_c/power_grid_model/include/power_grid_model/main_model_impl.hpp

@@ -25,7 +25,7 @@
 #include "auxiliary/output.hpp"
 
 // math model include
-#include "math_solver/math_solver.hpp"
+#include "math_solver/math_solver_dispatch.hpp"
 
 #include "optimizer/optimizer.hpp"
 
@@ -161,16 +161,22 @@ class MainModelImpl<ExtraRetrievableTypes<ExtraRetrievableType...>, ComponentLis
     using Options = MainModelOptions;
 
     // constructor with data
-    explicit MainModelImpl(double system_frequency, ConstDataset const& input_data, Idx pos = 0)
-        : system_frequency_{system_frequency}, meta_data_{&input_data.meta_data()} {
+    explicit MainModelImpl(double system_frequency, ConstDataset const& input_data,
+                           MathSolverDispatcher const& math_solver_dispatcher, Idx pos = 0)
+        : system_frequency_{system_frequency},
+          meta_data_{&input_data.meta_data()},
+          math_solver_dispatcher_{&math_solver_dispatcher} {
         assert(input_data.get_description().dataset->name == std::string_view("input"));
         add_components(input_data, pos);
         set_construction_complete();
     }
 
     // constructor with only frequency
-    explicit MainModelImpl(double system_frequency, meta_data::MetaData const& meta_data)
-        : system_frequency_{system_frequency}, meta_data_{&meta_data} {}
+    explicit MainModelImpl(double system_frequency, meta_data::MetaData const& meta_data,
+                           MathSolverDispatcher const& math_solver_dispatcher)
+        : system_frequency_{system_frequency},
+          meta_data_{&meta_data},
+          math_solver_dispatcher_{&math_solver_dispatcher} {}
 
   private:
     // helper function to get what components are present in the update data
@@ -441,25 +447,25 @@ class MainModelImpl<ExtraRetrievableTypes<ExtraRetrievableType...>, ComponentLis
     template <symmetry_tag sym> auto calculate_power_flow_(double err_tol, Idx max_iter) {
         return [this, err_tol, max_iter](MainModelState const& state,
                                          CalculationMethod calculation_method) -> std::vector<SolverOutput<sym>> {
-            return calculate_<SolverOutput<sym>, MathSolver<sym>, YBus<sym>, PowerFlowInput<sym>>(
+            return calculate_<SolverOutput<sym>, MathSolverProxy<sym>, YBus<sym>, PowerFlowInput<sym>>(
                 [&state](Idx n_math_solvers) { return prepare_power_flow_input<sym>(state, n_math_solvers); },
-                [this, err_tol, max_iter, calculation_method](MathSolver<sym>& solver, YBus<sym> const& y_bus,
+                [this, err_tol, max_iter, calculation_method](MathSolverProxy<sym>& solver, YBus<sym> const& y_bus,
                                                               PowerFlowInput<sym> const& input) {
-                    return solver.run_power_flow(input, err_tol, max_iter, calculation_info_, calculation_method,
-                                                 y_bus);
+                    return solver.get().run_power_flow(input, err_tol, max_iter, calculation_info_, calculation_method,
+                                                       y_bus);
                 });
         };
     }
 
     template <symmetry_tag sym> auto calculate_state_estimation_(double err_tol, Idx max_iter) {
         return [this, err_tol, max_iter](MainModelState const& state,
                                          CalculationMethod calculation_method) -> std::vector<SolverOutput<sym>> {
-            return calculate_<SolverOutput<sym>, MathSolver<sym>, YBus<sym>, StateEstimationInput<sym>>(
+            return calculate_<SolverOutput<sym>, MathSolverProxy<sym>, YBus<sym>, StateEstimationInput<sym>>(
                 [&state](Idx n_math_solvers) { return prepare_state_estimation_input<sym>(state, n_math_solvers); },
-                [this, err_tol, max_iter, calculation_method](MathSolver<sym>& solver, YBus<sym> const& y_bus,
+                [this, err_tol, max_iter, calculation_method](MathSolverProxy<sym>& solver, YBus<sym> const& y_bus,
                                                               StateEstimationInput<sym> const& input) {
-                    return solver.run_state_estimation(input, err_tol, max_iter, calculation_info_, calculation_method,
-                                                       y_bus);
+                    return solver.get().run_state_estimation(input, err_tol, max_iter, calculation_info_,
+                                                             calculation_method, y_bus);
                 });
         };
     }
@@ -468,14 +474,14 @@ class MainModelImpl<ExtraRetrievableTypes<ExtraRetrievableType...>, ComponentLis
         return [this,
                 voltage_scaling](MainModelState const& /*state*/,
                                  CalculationMethod calculation_method) -> std::vector<ShortCircuitSolverOutput<sym>> {
-            return calculate_<ShortCircuitSolverOutput<sym>, MathSolver<sym>, YBus<sym>, ShortCircuitInput>(
+            return calculate_<ShortCircuitSolverOutput<sym>, MathSolverProxy<sym>, YBus<sym>, ShortCircuitInput>(
                 [this, voltage_scaling](Idx /* n_math_solvers */) {
                     assert(is_topology_up_to_date_ && is_parameter_up_to_date<sym>());
                     return prepare_short_circuit_input<sym>(voltage_scaling);
                 },
-                [this, calculation_method](MathSolver<sym>& solver, YBus<sym> const& y_bus,
+                [this, calculation_method](MathSolverProxy<sym>& solver, YBus<sym> const& y_bus,
                                            ShortCircuitInput const& input) {
-                    return solver.run_short_circuit(input, calculation_info_, calculation_method, y_bus);
+                    return solver.get().run_short_circuit(input, calculation_info_, calculation_method, y_bus);
                 });
         };
     }
@@ -832,6 +838,7 @@ class MainModelImpl<ExtraRetrievableTypes<ExtraRetrievableType...>, ComponentLis
 
     double system_frequency_;
     meta_data::MetaData const* meta_data_;
+    MathSolverDispatcher const* math_solver_dispatcher_;
 
     MainModelState state_;
     // math model
@@ -859,7 +866,7 @@ class MainModelImpl<ExtraRetrievableTypes<ExtraRetrievableType...>, ComponentLis
         }
     }
 
-    template <symmetry_tag sym> std::vector<MathSolver<sym>>& get_solvers() {
+    template <symmetry_tag sym> std::vector<MathSolverProxy<sym>>& get_solvers() {
         if constexpr (is_symmetric_v<sym>) {
             return math_state_.math_solvers_sym;
         } else {
@@ -1291,7 +1298,7 @@ class MainModelImpl<ExtraRetrievableTypes<ExtraRetrievableType...>, ComponentLis
     }
 
     template <symmetry_tag sym> void prepare_solvers() {
-        std::vector<MathSolver<sym>>& solvers = get_solvers<sym>();
+        std::vector<MathSolverProxy<sym>>& solvers = get_solvers<sym>();
         // rebuild topology if needed
         if (!is_topology_up_to_date_) {
             rebuild_topology();
@@ -1305,11 +1312,14 @@ class MainModelImpl<ExtraRetrievableTypes<ExtraRetrievableType...>, ComponentLis
 
             solvers.clear();
             solvers.reserve(n_math_solvers_);
-            std::ranges::transform(state_.math_topology, std::back_inserter(solvers),
-                                   [](auto math_topo) { return MathSolver<sym>{std::move(math_topo)}; });
+            std::ranges::transform(state_.math_topology, std::back_inserter(solvers), [this](auto math_topo) {
+                return MathSolverProxy<sym>{math_solver_dispatcher_, std::move(math_topo)};
+            });
             for (Idx idx = 0; idx < n_math_solvers_; ++idx) {
                 get_y_bus<sym>()[idx].register_parameters_changed_callback(
-                    [solver = std::ref(solvers[idx])](bool changed) { solver.get().parameters_changed(changed); });
+                    [solver = std::ref(solvers[idx])](bool changed) {
+                        solver.get().get().parameters_changed(changed);
+                    });
             }
         } else if (!is_parameter_up_to_date<sym>()) {
             std::vector<MathModelParam<sym>> const math_params = get_math_param<sym>();

power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_current_pf_solver.hpp

@@ -217,9 +217,6 @@ class IterativeCurrentPFSolver : public IterativePFSolver<sym_type, IterativeCur
     }
 };
 
-template class IterativeCurrentPFSolver<symmetric_t>;
-template class IterativeCurrentPFSolver<asymmetric_t>;
-
 } // namespace iterative_current_pf
 
 using iterative_current_pf::IterativeCurrentPFSolver;

power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp

@@ -359,9 +359,6 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
     }
 };
 
-template class IterativeLinearSESolver<symmetric_t>;
-template class IterativeLinearSESolver<asymmetric_t>;
-
 } // namespace iterative_linear_se
 
 using iterative_linear_se::IterativeLinearSESolver;

power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/linear_pf_solver.hpp

@@ -109,8 +109,6 @@ template <symmetry_tag sym_type> class LinearPFSolver {
     }
 };
 
-template class LinearPFSolver<symmetric_t>;
-template class LinearPFSolver<asymmetric_t>;
 } // namespace linear_pf
 
 using linear_pf::LinearPFSolver;

power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/math_solver.hpp

@@ -7,6 +7,7 @@
 #include "iterative_current_pf_solver.hpp"
 #include "iterative_linear_se_solver.hpp"
 #include "linear_pf_solver.hpp"
+#include "math_solver_dispatch.hpp"
 #include "newton_raphson_pf_solver.hpp"
 #include "newton_raphson_se_solver.hpp"
 #include "short_circuit_solver.hpp"
@@ -24,16 +25,18 @@ namespace power_grid_model {
 
 namespace math_solver {
 
-template <symmetry_tag sym> class MathSolver {
+template <symmetry_tag sym> class MathSolver : public MathSolverBase<sym> {
   public:
     explicit MathSolver(std::shared_ptr<MathModelTopology const> const& topo_ptr)
         : topo_ptr_{topo_ptr},
           all_const_y_{std::all_of(topo_ptr->load_gen_type.cbegin(), topo_ptr->load_gen_type.cend(),
                                    [](LoadGenType x) { return x == LoadGenType::const_y; })} {}
 
+    MathSolver<sym>* clone() const final { return new MathSolver<sym>(*this); }
+
     SolverOutput<sym> run_power_flow(PowerFlowInput<sym> const& input, double err_tol, Idx max_iter,
                                      CalculationInfo& calculation_info, CalculationMethod calculation_method,
-                                     YBus<sym> const& y_bus) {
+                                     YBus<sym> const& y_bus) final {
         using enum CalculationMethod;
 
         // set method to always linear if all load_gens have const_y
@@ -57,7 +60,7 @@ template <symmetry_tag sym> class MathSolver {
 
     SolverOutput<sym> run_state_estimation(StateEstimationInput<sym> const& input, double err_tol, Idx max_iter,
                                            CalculationInfo& calculation_info, CalculationMethod calculation_method,
-                                           YBus<sym> const& y_bus) {
+                                           YBus<sym> const& y_bus) final {
         using enum CalculationMethod;
 
         switch (calculation_method) {
@@ -73,7 +76,8 @@ template <symmetry_tag sym> class MathSolver {
     }
 
     ShortCircuitSolverOutput<sym> run_short_circuit(ShortCircuitInput const& input, CalculationInfo& calculation_info,
-                                                    CalculationMethod calculation_method, YBus<sym> const& y_bus) {
+                                                    CalculationMethod calculation_method,
+                                                    YBus<sym> const& y_bus) final {
         if (calculation_method != CalculationMethod::default_method &&
             calculation_method != CalculationMethod::iec60909) {
             throw InvalidCalculationMethod{};
@@ -89,14 +93,14 @@ template <symmetry_tag sym> class MathSolver {
         return iec60909_sc_solver_.value().run_short_circuit(y_bus, input);
     }
 
-    void clear_solver() {
+    void clear_solver() final {
         newton_raphson_pf_solver_.reset();
         linear_pf_solver_.reset();
         iterative_current_pf_solver_.reset();
         iterative_linear_se_solver_.reset();
     }
 
-    void parameters_changed(bool changed) {
+    void parameters_changed(bool changed) final {
         if (iterative_current_pf_solver_.has_value()) {
             iterative_current_pf_solver_->parameters_changed(changed);
         }
@@ -175,8 +179,6 @@ template <symmetry_tag sym> class MathSolver {
     }
 };
 
-template class MathSolver<symmetric_t>;
-template class MathSolver<asymmetric_t>;
 } // namespace math_solver
 
 using math_solver::MathSolver;

power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/math_solver_dispatch.hpp

@@ -0,0 +1,116 @@
+// SPDX-FileCopyrightText: Contributors to the Power Grid Model project <powergridmodel@lfenergy.org>
+//
+// SPDX-License-Identifier: MPL-2.0
+
+// runtime dispatch for math solver
+// this can separate math solver into a different translation unit
+
+#pragma once
+
+#include "../calculation_parameters.hpp"
+#include "../common/common.hpp"
+#include "../common/exception.hpp"
+#include "../common/three_phase_tensor.hpp"
+#include "../common/timer.hpp"
+
+#include <memory>
+
+namespace power_grid_model {
+
+namespace math_solver {
+
+// forward declare YBus
+template <symmetry_tag sym> class YBus;
+
+// abstract base class
+template <symmetry_tag sym> class MathSolverBase {
+  public:
+    virtual ~MathSolverBase() = default;
+
+    virtual MathSolverBase<sym>* clone() const = 0;
+
+    virtual SolverOutput<sym> run_power_flow(PowerFlowInput<sym> const& input, double err_tol, Idx max_iter,
+                                             CalculationInfo& calculation_info, CalculationMethod calculation_method,
+                                             YBus<sym> const& y_bus) = 0;
+    virtual SolverOutput<sym> run_state_estimation(StateEstimationInput<sym> const& input, double err_tol, Idx max_iter,
+                                                   CalculationInfo& calculation_info,
+                                                   CalculationMethod calculation_method, YBus<sym> const& y_bus) = 0;
+    virtual ShortCircuitSolverOutput<sym> run_short_circuit(ShortCircuitInput const& input,
+                                                            CalculationInfo& calculation_info,
+                                                            CalculationMethod calculation_method,
+                                                            YBus<sym> const& y_bus) = 0;
+    virtual void clear_solver() = 0;
+    virtual void parameters_changed(bool changed) = 0;
+
+  protected:
+    MathSolverBase() = default;
+    MathSolverBase(MathSolverBase const&) = default;
+    MathSolverBase& operator=(MathSolverBase const&) = default;
+    MathSolverBase(MathSolverBase&&) noexcept = default;
+    MathSolverBase& operator=(MathSolverBase&&) noexcept = default;
+};
+
+// tag of math solver concrete types
+template <template <symmetry_tag> class MathSolverType> struct math_solver_tag {};
+
+class MathSolverDispatcher {
+  public:
+    template <symmetry_tag sym> struct Config {
+        template <template <class> class MathSolverType>
+        constexpr Config(math_solver_tag<MathSolverType> /* unused */)
+            : create{[](std::shared_ptr<MathModelTopology const> const& topo_ptr) -> MathSolverBase<sym>* {
+                  return new MathSolverType<sym>{topo_ptr};
+              }} {}
+
+        std::add_pointer_t<MathSolverBase<sym>*(std::shared_ptr<MathModelTopology const> const&)> create;
+    };
+
+    template <template <class> class MathSolverType>
+    constexpr MathSolverDispatcher(math_solver_tag<MathSolverType> tag) : sym_config_{tag}, asym_config_{tag} {}
+
+    template <symmetry_tag sym> Config<sym> const& get_config() const {
+        if constexpr (is_symmetric_v<sym>) {
+            return sym_config_;
+        } else {
+            return asym_config_;
+        }
+    }
+
+  private:
+    Config<symmetric_t> sym_config_;
+    Config<asymmetric_t> asym_config_;
+};
+
+template <symmetry_tag sym> class MathSolverProxy {
+  public:
+    explicit MathSolverProxy(MathSolverDispatcher const* dispatcher,
+                             std::shared_ptr<MathModelTopology const> const& topo_ptr)
+        : dispatcher_{dispatcher}, solver_{dispatcher_->get_config<sym>().create(topo_ptr)} {}
+    MathSolverProxy(MathSolverProxy const& other) : dispatcher_{other.dispatcher_}, solver_{other.get().clone()} {}
+    MathSolverProxy& operator=(MathSolverProxy const& other) {
+        if (this != &other) {
+            solver_.reset();
+            dispatcher_ = other.dispatcher_;
+            solver_.reset(other.get().clone());
+        }
+        return *this;
+    }
+    MathSolverProxy(MathSolverProxy&& other) noexcept = default;
+    MathSolverProxy& operator=(MathSolverProxy&& other) noexcept = default;
+    ~MathSolverProxy() = default;
+
+    MathSolverBase<sym>& get() { return *solver_; }
+    MathSolverBase<sym> const& get() const { return *solver_; }
+
+  private:
+    MathSolverDispatcher const* dispatcher_{};
+    std::unique_ptr<MathSolverBase<sym>> solver_;
+};
+
+} // namespace math_solver
+
+template <symmetry_tag sym> using MathSolverProxy = math_solver::MathSolverProxy<sym>;
+
+using MathSolverDispatcher = math_solver::MathSolverDispatcher;
+
+} // namespace power_grid_model

power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/newton_raphson_pf_solver.hpp

@@ -441,9 +441,6 @@ class NewtonRaphsonPFSolver : public IterativePFSolver<sym_type, NewtonRaphsonPF
     }
 };
 
-template class NewtonRaphsonPFSolver<symmetric_t>;
-template class NewtonRaphsonPFSolver<asymmetric_t>;
-
 } // namespace newton_raphson_pf
 
 using newton_raphson_pf::NewtonRaphsonPFSolver;