Feature/zip sparse grouped idx vector

power_grid_model_c/power_grid_model/include/power_grid_model/calculation_parameters.hpp

@@ -159,7 +159,7 @@ template <bool sym> struct StateEstimationInput {
 };
 
 struct ShortCircuitInput {
-    IdxVector fault_bus_indptr; // indptr of the fault
+    DenseGroupedIdxVector fault_buses;
     std::vector<FaultCalcParam> faults;
     ComplexVector source; // Complex u_ref of each source
 };

power_grid_model_c/power_grid_model/include/power_grid_model/grouped_index_vector.hpp

@@ -32,6 +32,8 @@ The input, ie. [0, 1, 3] should be strictly increasing
 
 namespace power_grid_model {
 
+using IdxRange = boost::iterator_range<IdxCount>;
+
 namespace detail {
 
 inline auto sparse_encode(IdxVector const& element_groups, Idx num_groups) {
@@ -46,7 +48,7 @@ inline auto sparse_encode(IdxVector const& element_groups, Idx num_groups) {
 
 inline auto sparse_decode(IdxVector const& indptr) {
     auto result = IdxVector(indptr.back());
-    for (Idx group{}; group < static_cast<Idx>(indptr.size()) - 1; ++group) {
+    for (Idx const group : boost::counting_range(Idx{}, static_cast<Idx>(indptr.size()) - 1)) {
         std::fill(std::begin(result) + indptr[group], std::begin(result) + indptr[group + 1], group);
     }
     return result;
@@ -116,10 +118,10 @@ constexpr auto from_dense = from_dense_t{};
 
 class SparseGroupedIdxVector {
   private:
-    class GroupIterator : public boost::iterator_facade<GroupIterator, Idx, boost::random_access_traversal_tag,
-                                                        boost::iterator_range<IdxCount>, Idx> {
+    class GroupIterator
+        : public boost::iterator_facade<GroupIterator, Idx, boost::random_access_traversal_tag, IdxRange, Idx> {
       public:
-        using iterator = boost::iterator_range<IdxCount>;
+        using iterator = IdxRange;
 
         GroupIterator() = default;
         explicit constexpr GroupIterator(IdxVector const& indptr, Idx group) : indptr_{&indptr}, group_{group} {}
@@ -177,13 +179,13 @@ class SparseGroupedIdxVector {
 
 class DenseGroupedIdxVector {
   private:
-    class GroupIterator : public boost::iterator_facade<GroupIterator, Idx, boost::random_access_traversal_tag,
-                                                        boost::iterator_range<IdxCount>, Idx> {
+    class GroupIterator
+        : public boost::iterator_facade<GroupIterator, Idx, boost::random_access_traversal_tag, IdxRange, Idx> {
       public:
-        using iterator = boost::iterator_range<IdxCount>;
+        using iterator = IdxRange;
 
         GroupIterator() = default;
-        explicit constexpr GroupIterator(IdxVector const& dense_vector, Idx const& group)
+        explicit constexpr GroupIterator(IdxVector const& dense_vector, Idx group)
             : dense_vector_{&dense_vector},
               group_{group},
               group_range_{std::equal_range(std::cbegin(*dense_vector_), std::cend(*dense_vector_), group)} {}

power_grid_model_c/power_grid_model/include/power_grid_model/main_model.hpp

@@ -1058,7 +1058,7 @@ class MainModelImpl<ExtraRetrievableTypes<ExtraRetrievableType...>, ComponentLis
                 fault_coup[topo_fault_indices[i][map.reorder[reordered_idx]]] = Idx2D{i, reordered_idx};
             }
 
-            sc_input[i].fault_bus_indptr = std::move(map.indptr);
+            sc_input[i].fault_buses = {from_sparse, std::move(map.indptr)};
             sc_input[i].faults.resize(state_.components.template size<Fault>());
             sc_input[i].source.resize(state_.math_topology[i]->n_source());
         }

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

@@ -13,9 +13,9 @@
 namespace power_grid_model::common_solver_functions {
 
 template <bool sym>
-void add_sources(grouped_idx_vector_type auto const& sources_per_bus, Idx const& bus_number, YBus<sym> const& y_bus,
+void add_sources(IdxRange const& sources, Idx /* bus_number */, YBus<sym> const& y_bus,
                  ComplexVector const& u_source_vector, ComplexTensor<sym>& diagonal_element, ComplexValue<sym>& u_bus) {
-    for (Idx const source_number : sources_per_bus.get_element_range(bus_number)) {
+    for (Idx const source_number : sources) {
         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
@@ -33,9 +33,9 @@ template <bool sym> void copy_y_bus(YBus<sym> const& y_bus, ComplexTensorVector<
 }
 
 template <bool sym>
-void calculate_source_result(Idx const& bus_number, YBus<sym> const& y_bus, PowerFlowInput<sym> const& input,
-                             MathOutput<sym>& output, grouped_idx_vector_type auto const& sources_per_bus) {
-    for (Idx const source : sources_per_bus.get_element_range(bus_number)) {
+void calculate_source_result(IdxRange const& sources, Idx bus_number, YBus<sym> const& y_bus,
+                             PowerFlowInput<sym> const& input, MathOutput<sym>& output) {
+    for (Idx const source : sources) {
         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]);
@@ -44,10 +44,11 @@ void calculate_source_result(Idx const& bus_number, YBus<sym> const& y_bus, Powe
 }
 
 template <bool sym, class LoadGenFunc>
-void calculate_load_gen_result(Idx const& bus_number, PowerFlowInput<sym> const& input, MathOutput<sym>& output,
-                               grouped_idx_vector_type auto const& load_gens_per_bus,
-                               LoadGenFunc const& load_gen_func) {
-    for (auto load_gen : load_gens_per_bus.get_element_range(bus_number)) {
+    requires std::invocable<std::remove_cvref_t<LoadGenFunc>, Idx> &&
+             std::same_as<std::invoke_result_t<LoadGenFunc, Idx>, LoadGenType>
+void calculate_load_gen_result(IdxRange const& load_gens, Idx bus_number, PowerFlowInput<sym> const& input,
+                               MathOutput<sym>& output, LoadGenFunc&& load_gen_func) {
+    for (Idx const load_gen : load_gens) {
         switch (LoadGenType const type = load_gen_func(load_gen); type) {
             using enum LoadGenType;
 
@@ -70,6 +71,33 @@ void calculate_load_gen_result(Idx const& bus_number, PowerFlowInput<sym> const&
     }
 }
 
+template <bool sym, typename LoadGenFunc>
+    requires std::invocable<std::remove_cvref_t<LoadGenFunc>, Idx> &&
+             std::same_as<std::invoke_result_t<LoadGenFunc, Idx>, LoadGenType>
+void calculate_result(YBus<sym> const& y_bus, PowerFlowInput<sym> const& input,
+                      grouped_idx_vector_type auto const& sources_per_bus,
+                      grouped_idx_vector_type auto const& load_gens_per_bus, MathOutput<sym>& output,
+                      LoadGenFunc&& load_gen_func) {
+    assert(sources_per_bus.size() == load_gens_per_bus.size());
+
+    // call y bus
+    output.branch = y_bus.template calculate_branch_flow<BranchMathOutput<sym>>(output.u);
+    output.shunt = y_bus.template calculate_shunt_flow<ApplianceMathOutput<sym>>(output.u);
+
+    // prepare source, load gen and node injection
+    output.source.resize(sources_per_bus.element_size());
+    output.load_gen.resize(load_gens_per_bus.element_size());
+    output.bus_injection.resize(sources_per_bus.size());
+
+    for (auto const& [bus_number, sources] : enumerated_zip_sequence(sources_per_bus)) {
+        common_solver_functions::calculate_source_result<sym>(sources, bus_number, y_bus, input, output);
+    }
+    for (auto const& [bus_number, load_gens] : enumerated_zip_sequence(load_gens_per_bus)) {
+        common_solver_functions::calculate_load_gen_result<sym>(load_gens, bus_number, input, output, load_gen_func);
+    }
+    output.bus_injection = y_bus.calculate_injection(output.u);
+}
+
 } // 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

@@ -93,11 +93,9 @@ template <bool sym> class IterativeCurrentPFSolver : public IterativePFSolver<sy
             ComplexTensorVector<sym> mat_data(y_bus.nnz_lu());
             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) {
+            for (auto const& [bus_number, sources] : enumerated_zip_sequence(sources_per_bus)) {
                 Idx const data_sequence = bus_entry[bus_number];
-                // loop sources
-                for (auto source_number : sources_per_bus.get_element_range(bus_number)) {
+                for (auto source_number : sources) {
                     // YBus_diag += Y_source
                     mat_data[data_sequence] += y_bus.math_model_param().source_param[source_number];
                 }
@@ -116,17 +114,17 @@ template <bool sym> class IterativeCurrentPFSolver : public IterativePFSolver<sy
     // Prepare matrix calculates injected current ie. RHS of solver for each iteration.
     void prepare_matrix_and_rhs(YBus<sym> const& y_bus, PowerFlowInput<sym> const& input,
                                 ComplexValueVector<sym> const& u) {
-        auto const& load_gens_per_bus = *this->load_gens_per_bus_;
-        auto const& sources_per_bus = *this->sources_per_bus_;
         std::vector<LoadGenType> const& load_gen_type = *this->load_gen_type_;
 
         // set rhs to zero for iteration start
         std::fill(rhs_u_.begin(), rhs_u_.end(), ComplexValue<sym>{0.0});
 
         // loop buses: i
-        for (Idx bus_number = 0; bus_number != this->n_bus_; ++bus_number) {
-            add_loads(bus_number, input, load_gens_per_bus, load_gen_type, u);
-            add_sources(bus_number, y_bus, input, sources_per_bus);
+        for (auto const& [bus_number, load_gens] : enumerated_zip_sequence(*this->load_gens_per_bus_)) {
+            add_loads(load_gens, bus_number, input, load_gen_type, u);
+        }
+        for (auto const& [bus_number, sources] : enumerated_zip_sequence(*this->sources_per_bus_)) {
+            add_sources(sources, bus_number, y_bus, input);
         }
     }
 
@@ -157,10 +155,9 @@ template <bool sym> class IterativeCurrentPFSolver : public IterativePFSolver<sy
     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,
-                   grouped_idx_vector_type auto const& load_gens_per_bus, std::vector<LoadGenType> const& load_gen_type,
-                   ComplexValueVector<sym> const& u) {
-        for (auto load_number : load_gens_per_bus.get_element_range(bus_number)) {
+    void add_loads(boost::iterator_range<IdxCount> const& load_gens, Idx bus_number, PowerFlowInput<sym> const& input,
+                   std::vector<LoadGenType> const& load_gen_type, ComplexValueVector<sym> const& u) {
+        for (Idx const load_number : load_gens) {
             // load type
             LoadGenType const type = load_gen_type[load_number];
             switch (type) {
@@ -184,9 +181,9 @@ template <bool sym> class IterativeCurrentPFSolver : public IterativePFSolver<sy
         }
     }
 
-    void add_sources(Idx const& bus_number, YBus<sym> const& y_bus, PowerFlowInput<sym> const& input,
-                     grouped_idx_vector_type auto const& sources_per_bus) {
-        for (Idx const source_number : sources_per_bus.get_element_range(bus_number)) {
+    void add_sources(boost::iterator_range<IdxCount> const& sources, Idx bus_number, YBus<sym> const& y_bus,
+                     PowerFlowInput<sym> const& input) {
+        for (Idx const source_number : sources) {
             // 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]});