Main model refactoring: move calculation logic out of main model

CMakeLists.txt

@@ -13,11 +13,18 @@ project(power_grid_model VERSION ${PGM_VERSION})
 
 option(PGM_ENABLE_DEV_BUILD "Enable developer build (e.g.: tests)" OFF)
 
-set(CMAKE_CXX_STANDARD 23)
+set(PGM_C_STANDARD 11)
+set(PGM_CXX_STANDARD 23)
+
+if(NOT CMAKE_CXX_STANDARD OR CMAKE_CXX_STANDARD LESS ${PGM_CXX_STANDARD})
+    set(CMAKE_CXX_STANDARD ${PGM_CXX_STANDARD})
+endif()
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
 set(CMAKE_CXX_EXTENSIONS OFF)
 set(CMAKE_CXX_VISIBILITY_PRESET hidden)
-set(CMAKE_C_STANDARD 11)
+if(NOT CMAKE_C_STANDARD OR CMAKE_C_STANDARD LESS ${PGM_C_STANDARD})
+    set(CMAKE_C_STANDARD ${PGM_C_STANDARD})
+endif()
 set(CMAKE_C_STANDARD_REQUIRED ON)
 set(CMAKE_C_EXTENSIONS OFF)
 set(CMAKE_C_VISIBILITY_PRESET hidden)

power_grid_model_c/power_grid_model/include/power_grid_model/calculation.hpp

@@ -0,0 +1,121 @@
+// SPDX-FileCopyrightText: Contributors to the Power Grid Model project <powergridmodel@lfenergy.org>
+//
+// SPDX-License-Identifier: MPL-2.0
+
+#pragma once
+
+#include "main_model_fwd.hpp"
+
+#include "common/common.hpp"
+#include "common/enum.hpp"
+
+#include "main_core/calculation_input_preparation.hpp"
+#include "main_core/main_model_type.hpp"
+
+namespace power_grid_model {
+
+struct calculation_type_t {};
+struct power_flow_t : calculation_type_t {};
+struct state_estimation_t : calculation_type_t {};
+struct short_circuit_t : calculation_type_t {};
+
+template <typename T>
+concept calculation_type_tag = std::derived_from<T, calculation_type_t>;
+
+template <class Functor, class... Args>
+decltype(auto) calculation_symmetry_func_selector(CalculationSymmetry calculation_symmetry, Functor&& f,
+                                                  Args&&... args) {
+    using enum CalculationSymmetry;
+
+    switch (calculation_symmetry) {
+    case symmetric:
+        return std::forward<Functor>(f).template operator()<symmetric_t>(std::forward<Args>(args)...);
+    case asymmetric:
+        return std::forward<Functor>(f).template operator()<asymmetric_t>(std::forward<Args>(args)...);
+    default:
+        throw MissingCaseForEnumError{"Calculation symmetry selector", calculation_symmetry};
+    }
+}
+
+template <class Functor, class... Args>
+decltype(auto) calculation_type_func_selector(CalculationType calculation_type, Functor&& f, Args&&... args) {
+    using enum CalculationType;
+
+    switch (calculation_type) {
+    case CalculationType::power_flow:
+        return std::forward<Functor>(f).template operator()<power_flow_t>(std::forward<Args>(args)...);
+    case CalculationType::state_estimation:
+        return std::forward<Functor>(f).template operator()<state_estimation_t>(std::forward<Args>(args)...);
+    case CalculationType::short_circuit:
+        return std::forward<Functor>(f).template operator()<short_circuit_t>(std::forward<Args>(args)...);
+    default:
+        throw MissingCaseForEnumError{"CalculationType", calculation_type};
+    }
+}
+
+template <class Functor, class... Args>
+decltype(auto) calculation_type_symmetry_func_selector(CalculationType calculation_type,
+                                                       CalculationSymmetry calculation_symmetry, Functor&& f,
+                                                       Args&&... args) {
+    calculation_type_func_selector(
+        calculation_type,
+        []<calculation_type_tag calculation_type, typename Functor_, typename... Args_>(
+            CalculationSymmetry calculation_symmetry_, Functor_&& f_, Args_&&... args_) {
+            calculation_symmetry_func_selector(
+                calculation_symmetry_,
+                []<symmetry_tag sym, typename SubFunctor, typename... SubArgs>(SubFunctor&& sub_f,
+                                                                               SubArgs&&... sub_args) {
+                    std::forward<SubFunctor>(sub_f).template operator()<calculation_type, sym>(
+                        std::forward<SubArgs>(sub_args)...);
+                },
+                std::forward<Functor_>(f_), std::forward<Args_>(args_)...);
+        },
+        calculation_symmetry, std::forward<Functor>(f), std::forward<Args>(args)...);
+}
+
+template <typename T, typename sym> struct Calculator;
+template <symmetry_tag sym> struct Calculator<power_flow_t, sym> {
+    template <typename State>
+    static auto preparer(State const& state, ComponentToMathCoupling& /*comp_coup*/,
+                         MainModelOptions const& /*options*/) {
+        return [&state](Idx n_math_solvers) { return main_core::prepare_power_flow_input<sym>(state, n_math_solvers); };
+    }
+    static auto solver(CalculationMethod calculation_method, MainModelOptions const& options, Logger& logger) {
+        return [calculation_method, err_tol = options.err_tol, max_iter = options.max_iter,
+                &logger](MathSolverProxy<sym>& solver, YBus<sym> const& y_bus, PowerFlowInput<sym> const& input) {
+            return solver.get().run_power_flow(input, err_tol, max_iter, logger, calculation_method, y_bus);
+        };
+    }
+};
+template <symmetry_tag sym> struct Calculator<state_estimation_t, sym> {
+    template <typename State>
+    static auto preparer(State const& state, ComponentToMathCoupling& /*comp_coup*/,
+                         MainModelOptions const& /*options*/) {
+        return [&state](Idx n_math_solvers) {
+            return main_core::prepare_state_estimation_input<sym>(state, n_math_solvers);
+        };
+    }
+    static auto solver(CalculationMethod calculation_method, MainModelOptions const& options, Logger& logger) {
+        return [calculation_method, err_tol = options.err_tol, max_iter = options.max_iter,
+                &logger](MathSolverProxy<sym>& solver, YBus<sym> const& y_bus, StateEstimationInput<sym> const& input) {
+            return solver.get().run_state_estimation(input, err_tol, max_iter, logger, calculation_method, y_bus);
+        };
+    }
+};
+
+template <symmetry_tag sym> struct Calculator<short_circuit_t, sym> {
+    template <typename State>
+    static auto preparer(State const& state, ComponentToMathCoupling& comp_coup, MainModelOptions const& options) {
+        return [&state, &comp_coup, voltage_scaling = options.short_circuit_voltage_scaling](Idx n_math_solvers) {
+            return main_core::prepare_short_circuit_input<sym>(state, comp_coup, n_math_solvers, voltage_scaling);
+        };
+    }
+    static auto solver(CalculationMethod calculation_method, MainModelOptions const& /*options*/, Logger& logger) {
+        return [calculation_method, &logger](MathSolverProxy<sym>& solver, YBus<sym> const& y_bus,
+                                             ShortCircuitInput const& input) {
+            return solver.get().run_short_circuit(input, logger, calculation_method, y_bus);
+        };
+    }
+};
+
+} // namespace power_grid_model

power_grid_model_c/power_grid_model/include/power_grid_model/common/counting_iterator.hpp

@@ -21,4 +21,12 @@ struct IdxRange : public std::ranges::iota_view<Idx, Idx> {
 };
 using IdxCount = typename IdxRange::iterator;
 
+#if __cpp_lib_ranges_enumerate >= 20202L
+using std::views::enumerate;
+#else
+constexpr auto enumerate(std::ranges::viewable_range auto&& rng) {
+    return std::views::zip(IdxRange{std::ssize(rng)}, rng);
+}
+#endif
+
 } // namespace power_grid_model

power_grid_model_c/power_grid_model/include/power_grid_model/main_core/calculation_input_preparation.hpp

@@ -5,8 +5,10 @@
 #pragma once
 
 #include "state_queries.hpp"
+#include "y_bus.hpp"
 
 #include "../calculation_parameters.hpp"
+#include "../index_mapping.hpp"
 
 #include <concepts>
 #include <vector>
@@ -258,7 +260,7 @@ std::vector<ShortCircuitInput> prepare_short_circuit_input(main_model_state_c au
     std::vector<IdxVector> topo_fault_indices(state.math_topology.size());
     std::vector<IdxVector> topo_bus_indices(state.math_topology.size());
 
-    for (Idx fault_idx{0}; fault_idx < state.components.template size<Fault>(); ++fault_idx) {
+    for (Idx const fault_idx : IdxRange{state.components.template size<Fault>()}) {
         auto const& fault = state.components.template get_item_by_seq<Fault>(fault_idx);
         if (fault.status()) {
             auto const node_idx = state.components.template get_seq<Node>(fault.get_fault_object());
@@ -275,11 +277,11 @@ std::vector<ShortCircuitInput> prepare_short_circuit_input(main_model_state_c au
                                          Idx2D{.group = isolated_component, .pos = not_connected});
     std::vector<ShortCircuitInput> sc_input(n_math_solvers);
 
-    for (Idx i = 0; i != n_math_solvers; ++i) {
+    for (Idx const i : IdxRange{n_math_solvers}) {
         auto map = build_dense_mapping(topo_bus_indices[i], state.math_topology[i]->n_bus());
 
-        for (Idx reordered_idx{0}; reordered_idx < static_cast<Idx>(map.reorder.size()); ++reordered_idx) {
-            fault_coup[topo_fault_indices[i][map.reorder[reordered_idx]]] = Idx2D{.group = i, .pos = reordered_idx};
+        for (auto&& [reordered_idx, original_idx] : enumerate(map.reorder)) {
+            fault_coup[topo_fault_indices[i][original_idx]] = Idx2D{.group = i, .pos = reordered_idx};
         }
 
         sc_input[i].fault_buses = {from_dense, std::move(map.indvector), state.math_topology[i]->n_bus()};

power_grid_model_c/power_grid_model/include/power_grid_model/main_core/state_queries.hpp

@@ -9,7 +9,6 @@
 #include "../all_components.hpp"
 
 namespace power_grid_model::main_core {
-
 template <std::derived_from<Branch> ComponentType, class ComponentContainer>
     requires model_component_state_c<MainModelState, ComponentContainer, ComponentType>
 constexpr auto get_branch_nodes(MainModelState<ComponentContainer> const& state, Idx topology_sequence_idx) {

power_grid_model_c/power_grid_model/include/power_grid_model/main_core/y_bus.hpp

@@ -4,6 +4,8 @@
 
 #pragma once
 
+#include "math_state.hpp"
+
 #include "../common/common.hpp"
 
 namespace power_grid_model::main_core {