Refactor solvers

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

@@ -0,0 +1,75 @@
+// SPDX-FileCopyrightText: 2022 Contributors to the Power Grid Model project <dynamic.grid.calculation@alliander.com>
+//
+// SPDX-License-Identifier: MPL-2.0
+
+#pragma once
+#ifndef POWER_GRID_MODEL_COMMON_SOLVER_FUNCTIONS_HPP
+#define POWER_GRID_MODEL_COMMON_SOLVER_FUNCTIONS_HPP
+
+#include "y_bus.hpp"
+
+#include "../calculation_parameters.hpp"
+
+namespace power_grid_model::common_solver_functions {
+
+template <bool sym>
+void add_sources(IdxVector const& source_bus_indptr, Idx const& bus_number, YBus<sym> const& y_bus,
+                 ComplexVector const& u_source_vector, ComplexTensor<sym>& diagonal_element, ComplexValue<sym>& u_bus) {
+    for (Idx source_number = source_bus_indptr[bus_number]; source_number != source_bus_indptr[bus_number + 1];
+         ++source_number) {
+        ComplexTensor<sym> const y_source = y_bus.math_model_param().source_param[source_number];
+        diagonal_element += y_source; // add y_source to the diagonal of Ybus
+        u_bus += dot(y_source, ComplexValue<sym>{u_source_vector[source_number]}); // rhs += Y_source * U_source
+    }
+}
+
+template <bool sym> void copy_y_bus(YBus<sym> const& y_bus, ComplexTensorVector<sym>& mat_data) {
+    ComplexTensorVector<sym> const& ydata = y_bus.admittance();
+    std::transform(y_bus.map_lu_y_bus().cbegin(), y_bus.map_lu_y_bus().cend(), mat_data.begin(), [&](Idx k) {
+        if (k == -1) {
+            return ComplexTensor<sym>{};
+        }
+        return ydata[k];
+    });
+}
+
+template <bool sym>
+void calculate_source_result(Idx const& bus_number, YBus<sym> const& y_bus, PowerFlowInput<sym> const& input,
+                             MathOutput<sym>& output, IdxVector const& source_bus_indptr) {
+    for (Idx source = source_bus_indptr[bus_number]; source != source_bus_indptr[bus_number + 1]; ++source) {
+        ComplexValue<sym> const u_ref{input.source[source]};
+        ComplexTensor<sym> const y_ref = y_bus.math_model_param().source_param[source];
+        output.source[source].i = dot(y_ref, u_ref - output.u[bus_number]);
+        output.source[source].s = output.u[bus_number] * conj(output.source[source].i);
+    }
+}
+
+template <bool sym, class LoadGenFunc>
+void calculate_load_gen_result(Idx const& bus_number, PowerFlowInput<sym> const& input, MathOutput<sym>& output,
+                               IdxVector const& load_gen_bus_indptr, LoadGenFunc const& load_gen_func) {
+    for (Idx load_gen = load_gen_bus_indptr[bus_number]; load_gen != load_gen_bus_indptr[bus_number + 1]; ++load_gen) {
+        switch (LoadGenType const type = load_gen_func(load_gen); type) {
+            using enum LoadGenType;
+
+        case const_pq:
+            // always same power
+            output.load_gen[load_gen].s = input.s_injection[load_gen];
+            break;
+        case const_y:
+            // power is quadratic relation to voltage
+            output.load_gen[load_gen].s = input.s_injection[load_gen] * abs2(output.u[bus_number]);
+            break;
+        case const_i:
+            // power is linear relation to voltage
+            output.load_gen[load_gen].s = input.s_injection[load_gen] * cabs(output.u[bus_number]);
+            break;
+        default:
+            throw MissingCaseForEnumError("Power injection", type);
+        }
+        output.load_gen[load_gen].i = conj(output.load_gen[load_gen].s / output.u[bus_number]);
+    }
+}
+
+} // namespace power_grid_model::common_solver_functions
+
+#endif

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

@@ -55,6 +55,7 @@ Iterative Power Flow
 */
 
 #include "block_matrix.hpp"
+#include "common_solver_functions.hpp"
 #include "iterative_pf_solver.hpp"
 #include "sparse_lu_solver.hpp"
 #include "y_bus.hpp"
@@ -85,19 +86,12 @@ template <bool sym> class IterativeCurrentPFSolver : public IterativePFSolver<sy
     // Add source admittance to Y bus and set variable for prepared y bus to true
     void initialize_derived_solver(YBus<sym> const& y_bus, MathOutput<sym> const& /* output */) {
         IdxVector const& source_bus_indptr = *this->source_bus_indptr_;
-        ComplexTensorVector<sym> const& ydata = y_bus.admittance();
         IdxVector const& bus_entry = y_bus.lu_diag();
         // if Y bus is not up to date
         // re-build matrix and prefactorize Build y bus data with source admittance
         if (y_data_ptr_ != &y_bus.admittance()) {
             ComplexTensorVector<sym> mat_data(y_bus.nnz_lu());
-            // copy y bus data
-            std::transform(y_bus.map_lu_y_bus().cbegin(), y_bus.map_lu_y_bus().cend(), mat_data.begin(), [&](Idx k) {
-                if (k == -1) {
-                    return ComplexTensor<sym>{};
-                }
-                return ydata[k];
-            });
+            common_solver_functions::copy_y_bus<sym>(y_bus, mat_data);
 
             // loop bus
             for (Idx bus_number = 0; bus_number != this->n_bus_; ++bus_number) {
@@ -132,37 +126,8 @@ template <bool sym> class IterativeCurrentPFSolver : public IterativePFSolver<sy
 
         // loop buses: i
         for (Idx bus_number = 0; bus_number != this->n_bus_; ++bus_number) {
-            // loop loads/generation: j
-            for (Idx load_number = load_gen_bus_indptr[bus_number]; load_number != load_gen_bus_indptr[bus_number + 1];
-                 ++load_number) {
-                // load type
-                LoadGenType const type = load_gen_type[load_number];
-                switch (type) {
-                    using enum LoadGenType;
-
-                case const_pq:
-                    // I_inj_i = conj(S_inj_j/U_i) for constant PQ type
-                    rhs_u_[bus_number] += conj(input.s_injection[load_number] / u[bus_number]);
-                    break;
-                case const_y:
-                    // I_inj_i = conj((S_inj_j * abs(U_i)^2) / U_i) = conj((S_inj_j) * U_i for const impedance type
-                    rhs_u_[bus_number] += conj(input.s_injection[load_number]) * u[bus_number];
-                    break;
-                case const_i:
-                    // I_inj_i = conj(S_inj_j*abs(U_i)/U_i) for const current type
-                    rhs_u_[bus_number] += conj(input.s_injection[load_number] * cabs(u[bus_number]) / u[bus_number]);
-                    break;
-                default:
-                    throw MissingCaseForEnumError("Injection current calculation", type);
-                }
-            }
-            // loop sources: j
-            for (Idx source_number = source_bus_indptr[bus_number]; source_number != source_bus_indptr[bus_number + 1];
-                 ++source_number) {
-                // I_inj_i += Y_source_j * U_ref_j
-                rhs_u_[bus_number] += dot(y_bus.math_model_param().source_param[source_number],
-                                          ComplexValue<sym>{input.source[source_number]});
-            }
+            add_loads(bus_number, input, load_gen_bus_indptr, load_gen_type, u);
+            add_sources(bus_number, y_bus, input, source_bus_indptr);
         }
     }
 
@@ -192,6 +157,43 @@ template <bool sym> class IterativeCurrentPFSolver : public IterativePFSolver<sy
     // sparse solver
     SparseLUSolver<ComplexTensor<sym>, ComplexValue<sym>, ComplexValue<sym>> sparse_solver_;
     std::shared_ptr<BlockPermArray const> perm_;
+
+    void add_loads(Idx const& bus_number, PowerFlowInput<sym> const& input, IdxVector const& load_gen_bus_indptr,
+                   std::vector<LoadGenType> const& load_gen_type, ComplexValueVector<sym> const& u) {
+        for (Idx load_number = load_gen_bus_indptr[bus_number]; load_number != load_gen_bus_indptr[bus_number + 1];
+             ++load_number) {
+            // load type
+            LoadGenType const type = load_gen_type[load_number];
+            switch (type) {
+                using enum LoadGenType;
+
+            case const_pq:
+                // I_inj_i = conj(S_inj_j/U_i) for constant PQ type
+                rhs_u_[bus_number] += conj(input.s_injection[load_number] / u[bus_number]);
+                break;
+            case const_y:
+                // I_inj_i = conj((S_inj_j * abs(U_i)^2) / U_i) = conj((S_inj_j) * U_i for const impedance type
+                rhs_u_[bus_number] += conj(input.s_injection[load_number]) * u[bus_number];
+                break;
+            case const_i:
+                // I_inj_i = conj(S_inj_j*abs(U_i)/U_i) for const current type
+                rhs_u_[bus_number] += conj(input.s_injection[load_number] * cabs(u[bus_number]) / u[bus_number]);
+                break;
+            default:
+                throw MissingCaseForEnumError("Injection current calculation", type);
+            }
+        }
+    }
+
+    void add_sources(Idx const& bus_number, YBus<sym> const& y_bus, PowerFlowInput<sym> const& input,
+                     IdxVector const& source_bus_indptr) {
+        for (Idx source_number = source_bus_indptr[bus_number]; source_number != source_bus_indptr[bus_number + 1];
+             ++source_number) {
+            // I_inj_i += Y_source_j * U_ref_j
+            rhs_u_[bus_number] += dot(y_bus.math_model_param().source_param[source_number],
+                                      ComplexValue<sym>{input.source[source_number]});
+        }
+    }
 };
 
 template class IterativeCurrentPFSolver<true>;

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

@@ -11,6 +11,7 @@
  */
 
 // Check if all includes needed
+#include "common_solver_functions.hpp"
 #include "y_bus.hpp"
 
 #include "../calculation_parameters.hpp"
@@ -112,41 +113,11 @@ template <bool sym, typename DerivedSolver> class IterativePFSolver {
         output.load_gen.resize(load_gen_bus_indptr_->back());
         output.bus_injection.resize(n_bus_);
 
-        // loop all bus
-        for (Idx bus = 0; bus != n_bus_; ++bus) {
-            // source
-            for (Idx source = (*source_bus_indptr_)[bus]; source != (*source_bus_indptr_)[bus + 1]; ++source) {
-                ComplexValue<sym> const u_ref{input.source[source]};
-                ComplexTensor<sym> const y_ref = y_bus.math_model_param().source_param[source];
-                output.source[source].i = dot(y_ref, u_ref - output.u[bus]);
-                output.source[source].s = output.u[bus] * conj(output.source[source].i);
-            }
-
-            // load_gen
-            for (Idx load_gen = (*load_gen_bus_indptr_)[bus]; load_gen != (*load_gen_bus_indptr_)[bus + 1];
-                 ++load_gen) {
-                LoadGenType const type = (*load_gen_type_)[load_gen];
-                switch (type) {
-                    using enum LoadGenType;
-
-                case const_pq:
-                    // always same power
-                    output.load_gen[load_gen].s = input.s_injection[load_gen];
-                    break;
-                case const_y:
-                    // power is quadratic relation to voltage
-                    output.load_gen[load_gen].s =
-                        input.s_injection[load_gen] * cabs(output.u[bus]) * cabs(output.u[bus]);
-                    break;
-                case const_i:
-                    // power is linear relation to voltage
-                    output.load_gen[load_gen].s = input.s_injection[load_gen] * cabs(output.u[bus]);
-                    break;
-                default:
-                    throw MissingCaseForEnumError("Power injection", type);
-                }
-                output.load_gen[load_gen].i = conj(output.load_gen[load_gen].s / output.u[bus]);
-            }
+        for (Idx bus_number = 0; bus_number != n_bus_; ++bus_number) {
+            common_solver_functions::calculate_source_result<sym>(bus_number, y_bus, input, output,
+                                                                  *source_bus_indptr_);
+            common_solver_functions::calculate_load_gen_result<sym>(bus_number, input, output, *load_gen_bus_indptr_,
+                                                                    [this](Idx i) { return (*load_gen_type_)[i]; });
         }
         output.bus_injection = y_bus.calculate_injection(output.u);
     }

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

@@ -30,6 +30,7 @@ if there are sources
 
 */
 
+#include "common_solver_functions.hpp"
 #include "sparse_lu_solver.hpp"
 #include "y_bus.hpp"
 
@@ -57,9 +58,6 @@ template <bool sym> class LinearPFSolver {
 
     MathOutput<sym> run_power_flow(YBus<sym> const& y_bus, PowerFlowInput<sym> const& input,
                                    CalculationInfo& calculation_info) {
-        // getter
-        ComplexTensorVector<sym> const& ydata = y_bus.admittance();
-        IdxVector const& bus_entry = y_bus.lu_diag();
         // output
         MathOutput<sym> output;
         output.u.resize(n_bus_);
@@ -68,36 +66,8 @@ template <bool sym> class LinearPFSolver {
 
         // prepare matrix
         Timer sub_timer(calculation_info, 2221, "Prepare matrix");
-
-        // copy y bus data
-        std::transform(y_bus.map_lu_y_bus().cbegin(), y_bus.map_lu_y_bus().cend(), mat_data_.begin(), [&](Idx k) {
-            if (k == -1) {
-                return ComplexTensor<sym>{};
-            }
-            return ydata[k];
-        });
-
-        // loop to all loads and sources, j as load number
-        IdxVector const& load_gen_bus_idxptr = *load_gen_bus_indptr_;
-        IdxVector const& source_bus_indptr = *source_bus_indptr_;
-        for (Idx bus_number = 0; bus_number != n_bus_; ++bus_number) {
-            Idx const data_sequence = bus_entry[bus_number];
-            // loop loads
-            for (Idx load_number = load_gen_bus_idxptr[bus_number]; load_number != load_gen_bus_idxptr[bus_number + 1];
-                 ++load_number) {
-                // YBus_diag += -conj(S_base)
-                add_diag(mat_data_[data_sequence], -conj(input.s_injection[load_number]));
-            }
-            // loop sources
-            for (Idx source_number = source_bus_indptr[bus_number]; source_number != source_bus_indptr[bus_number + 1];
-                 ++source_number) {
-                // YBus_diag += Y_source
-                mat_data_[data_sequence] += y_bus.math_model_param().source_param[source_number];
-                // rhs += Y_source_j * U_ref_j
-                output.u[bus_number] += dot(y_bus.math_model_param().source_param[source_number],
-                                            ComplexValue<sym>{input.source[source_number]});
-            }
-        }
+        common_solver_functions::copy_y_bus<sym>(y_bus, mat_data_);
+        prepare_matrix_and_rhs(y_bus, input, output);
 
         // solve
         // u vector will have I_injection for slack bus for now
@@ -123,6 +93,30 @@ template <bool sym> class LinearPFSolver {
     SparseLUSolver<ComplexTensor<sym>, ComplexValue<sym>, ComplexValue<sym>> sparse_solver_;
     typename SparseLUSolver<ComplexTensor<sym>, ComplexValue<sym>, ComplexValue<sym>>::BlockPermArray perm_;
 
+    void prepare_matrix_and_rhs(YBus<sym> const& y_bus, PowerFlowInput<sym> const& input, MathOutput<sym>& output) {
+        // getter
+        IdxVector const& bus_entry = y_bus.lu_diag();
+        IdxVector const& load_gen_bus_idxptr = *load_gen_bus_indptr_;
+        IdxVector const& source_bus_indptr = *source_bus_indptr_;
+        for (Idx bus_number = 0; bus_number != n_bus_; ++bus_number) {
+            Idx const diagonal_position = bus_entry[bus_number];
+            auto& diagonal_element = mat_data_[diagonal_position];
+            auto& u_bus = output.u[bus_number];
+            add_loads(load_gen_bus_idxptr, bus_number, input, diagonal_element);
+            common_solver_functions::add_sources<sym>(source_bus_indptr, bus_number, y_bus, input.source,
+                                                      diagonal_element, u_bus);
+        }
+    }
+
+    static void add_loads(IdxVector const& load_gen_bus_idxptr, Idx const& bus_number, PowerFlowInput<sym> const& input,
+                          ComplexTensor<sym>& diagonal_element) {
+        for (Idx load_number = load_gen_bus_idxptr[bus_number]; load_number != load_gen_bus_idxptr[bus_number + 1];
+             ++load_number) {
+            // YBus_diag += -conj(S_base)
+            add_diag(diagonal_element, -conj(input.s_injection[load_number]));
+        }
+    }
+
     void calculate_result(YBus<sym> const& y_bus, PowerFlowInput<sym> const& input, MathOutput<sym>& output) {
         // call y bus
         output.branch = y_bus.template calculate_branch_flow<BranchMathOutput<sym>>(output.u);
@@ -133,23 +127,11 @@ template <bool sym> class LinearPFSolver {
         output.load_gen.resize(load_gen_bus_indptr_->back());
         output.bus_injection.resize(n_bus_);
 
-        // loop all bus
-        for (Idx bus = 0; bus != n_bus_; ++bus) {
-            // source
-            for (Idx source = (*source_bus_indptr_)[bus]; source != (*source_bus_indptr_)[bus + 1]; ++source) {
-                ComplexValue<sym> const u_ref{input.source[source]};
-                ComplexTensor<sym> const y_ref = y_bus.math_model_param().source_param[source];
-                output.source[source].i = dot(y_ref, u_ref - output.u[bus]);
-                output.source[source].s = output.u[bus] * conj(output.source[source].i);
-            }
-
-            // load_gen
-            for (Idx load_gen = (*load_gen_bus_indptr_)[bus]; load_gen != (*load_gen_bus_indptr_)[bus + 1];
-                 ++load_gen) {
-                // power is always quadratic relation to voltage for linear pf
-                output.load_gen[load_gen].s = input.s_injection[load_gen] * abs2(output.u[bus]);
-                output.load_gen[load_gen].i = conj(output.load_gen[load_gen].s / output.u[bus]);
-            }
+        for (Idx bus_number = 0; bus_number != n_bus_; ++bus_number) {
+            common_solver_functions::calculate_source_result<sym>(bus_number, y_bus, input, output,
+                                                                  *source_bus_indptr_);
+            common_solver_functions::calculate_load_gen_result<sym>(bus_number, input, output, *load_gen_bus_indptr_,
+                                                                    [](Idx /*i*/) { return LoadGenType::const_y; });
         }
         output.bus_injection = y_bus.calculate_injection(output.u);
     }