Merge links into nodes: Create link output from link solver algorithm

.pre-commit-config.yaml

@@ -15,7 +15,7 @@ repos:
         entry: uv run --frozen pytest
         language: system
         pass_filenames: false
-        types_or: [ python, c++, c ]
+        types_or: [ python, c++, c, json ]
       - id: ruff-format
         name: ruff-format
         entry: uv run --frozen ruff format

power_grid_model_c/power_grid_model/include/power_grid_model/calculation_parameters.hpp

@@ -19,7 +19,8 @@
 
 namespace power_grid_model {
 constexpr Idx disconnected = -1;
-constexpr Idx status_off = 0;
+constexpr IntS status_on = 1;
+constexpr IntS status_off = 0;
 
 // Entry of YBus, node addmittance matrix
 struct YBusElement {

power_grid_model_c/power_grid_model/include/power_grid_model/component/link.hpp

@@ -30,6 +30,18 @@ class Link final : public Edge {
     double loading(double /* max_s */, double /* max_i */) const override { return 0.0; };
     double phase_shift() const override { return 0.0; }
 
+    template <symmetry_tag sym> BranchOutput<sym> get_energized_zero_output() const {
+        BranchOutput<sym> output = get_null_output<sym>();
+        static_cast<BaseOutput&>(output) = base_output(true);
+        return output;
+    }
+
+    BranchShortCircuitOutput get_energized_zero_sc_output() const {
+        BranchShortCircuitOutput output = get_null_sc_output();
+        static_cast<BaseOutput&>(output) = base_output(true);
+        return output;
+    }
+
   private:
     double base_i_from_;
     double base_i_to_;

power_grid_model_c/power_grid_model/include/power_grid_model/container.hpp

@@ -68,6 +68,7 @@ template <class... T> class Container;
 template <class... GettableTypes, class... StorageableTypes>
 class Container<RetrievableTypes<GettableTypes...>, StorageableTypes...> {
   public:
+    using storageable_types = std::tuple<StorageableTypes...>;
     using gettable_types = std::tuple<GettableTypes...>;
 
     static constexpr size_t num_storageable = sizeof...(StorageableTypes);

power_grid_model_c/power_grid_model/include/power_grid_model/main_core/math_output_queries.hpp

@@ -21,27 +21,23 @@ namespace power_grid_model::main_core {
 
 template <typename Component, solver_output_type SolverOutputType>
 constexpr auto const& get_component_output(MathOutput<std::vector<SolverOutputType>> const& math_output,
-                                           Idx2D const& math_id) {
+                                           Idx2D math_id) {
     auto const& solver_output = math_output.solver_output[math_id.group];
 
-    auto const& component_type_output = [&solver_output]() -> auto const& {
-        // TODO(mgovers): cleanup v2: change back to std::derived_from<Component, Branch>
-        if constexpr (std::derived_from<Component, Edge> || std::derived_from<Component, Branch3>) {
-            return solver_output.branch;
-        } else if constexpr (std::same_as<Component, Source> && requires { solver_output.source; }) {
-            return solver_output.source;
-        } else if constexpr (std::same_as<Component, Shunt> && requires { solver_output.shunt; }) {
-            return solver_output.shunt;
-        } else if constexpr (std::derived_from<Component, GenericLoadGen> && requires { solver_output.load_gen; }) {
-            return solver_output.load_gen;
-        } else if constexpr (std::same_as<Component, Fault> && requires { solver_output.fault; }) {
-            return solver_output.fault;
-        } else {
-            static_assert(false, "Unsupported component type for output retrieval");
-        }
-    }();
-
-    return component_type_output[math_id.pos];
+    // TODO(mgovers): cleanup v2: change back to std::derived_from<Component, Branch>
+    if constexpr (std::derived_from<Component, Edge> || std::derived_from<Component, Branch3>) {
+        return solver_output.branch[math_id.pos];
+    } else if constexpr (std::same_as<Component, Source> && requires { solver_output.source; }) {
+        return solver_output.source[math_id.pos];
+    } else if constexpr (std::same_as<Component, Shunt> && requires { solver_output.shunt; }) {
+        return solver_output.shunt[math_id.pos];
+    } else if constexpr (std::derived_from<Component, GenericLoadGen> && requires { solver_output.load_gen; }) {
+        return solver_output.load_gen[math_id.pos];
+    } else if constexpr (std::same_as<Component, Fault> && requires { solver_output.fault; }) {
+        return solver_output.fault[math_id.pos];
+    } else {
+        static_assert(false, "Unsupported component type for output retrieval");
+    }
 }
 
 } // namespace power_grid_model::main_core

power_grid_model_c/power_grid_model/include/power_grid_model/main_core/output.hpp

@@ -128,8 +128,37 @@
     return node.get_sc_output(math_output.solver_output[math_id.group].u_bus[math_id.pos]);
 }
 
+// output link
+template <std::same_as<Link> Component, class ComponentContainer, steady_state_solver_output_type SolverOutputType>
+    requires model_component_state_c<MainModelState, ComponentContainer, Component>
+constexpr auto output_result(Component const& link, MainModelState<ComponentContainer> const& /* state */,
+                             MathOutput<std::vector<SolverOutputType>> const& math_output, Idx2D const& topo_id) {
+    using sym = decode_symmetry_v<SolverOutputType>;
+
+    if (topo_id.group == disconnected) {
+        return link.template get_null_output<sym>();
+    }
+    if (!link.edge_status()) {
+        return link.template get_energized_zero_output<sym>();
+    }
+    return link.template get_output<sym>(math_output.supernode_output[topo_id.group].link[topo_id.pos]);
+}
+template <std::same_as<Link> Component, class ComponentContainer, short_circuit_solver_output_type SolverOutputType>
+    requires model_component_state_c<MainModelState, ComponentContainer, Component>
+inline auto output_result(Component const& link, MainModelState<ComponentContainer> const& /* state */,
+                          MathOutput<std::vector<SolverOutputType>> const& math_output, Idx2D const& topo_id) {
+    if (topo_id.group == disconnected) {
+        return link.get_null_sc_output();
+    }
+    if (!link.edge_status()) {
+        return link.get_energized_zero_sc_output();
+    }
+    return link.get_sc_output(math_output.supernode_output[topo_id.group].link[topo_id.pos]);
+}
+
 // output branch
 template <std::derived_from<Edge> Component, steady_state_solver_output_type SolverOutputType>
+    requires(!std::same_as<Component, Link>) // TODO(mgovers): cleanup v2: change back to only derived_from<Branch>
 constexpr auto output_result(Component const& branch, std::vector<SolverOutputType> const& solver_output,
                              Idx2D math_id) {
     using sym = decode_symmetry_v<SolverOutputType>;
@@ -139,7 +168,9 @@
     }
     return branch.template get_output<sym>(solver_output[math_id.group].branch[math_id.pos]);
 }
+// TODO(mgovers): cleanup v2: change back to only derived_from<Branch>
 template <std::derived_from<Edge> Component, short_circuit_solver_output_type SolverOutputType>
+    requires(!std::same_as<Component, Link>) // TODO(mgovers): cleanup v2: change back to only derived_from<Branch>
 inline auto output_result(Component const& branch, std::vector<SolverOutputType> const& solver_output, Idx2D math_id) {
     if (math_id.group == disconnected) {
         return branch.get_null_sc_output();
@@ -473,6 +504,28 @@
     return voltage_regulator.get_null_sc_output();
 }
 
+template <std::same_as<Link> Component, class ComponentContainer, solver_output_type SolverOutputType,
+          non_owning_view_c ComponentOutput>
+    requires model_component_state_c<MainModelState, ComponentContainer, Component>
+constexpr void output_result(MainModelState<ComponentContainer> const& state,
+                             MathOutput<std::vector<SolverOutputType>> const& math_output, ComponentOutput output) {
+    if (auto const& link_topo_ids = state.reduced_topology->topo_node_coup.coupling.user_links_to_topo_nodes;
+        std::ranges::ssize(link_topo_ids) ==
+        get_component_size<Component>(
+            state.components)) { // TODO(mgovers): cleanup v2: this should be the only code path remaining
+        std::ranges::transform(
+            get_component_citer<Component>(state.components), link_topo_ids, std::ranges::begin(output),
+            [&state, &math_output](Component const& link, Idx2D const& topo_id) {
+                return output_result<Component, ComponentContainer>(link, state, math_output, topo_id);
+            });
+    } else {
+        detail::produce_output<Component, Idx2D>(
+            state, output, [&math_output](Component const& link, Idx2D const& math_id) {
+                return output_result<Edge>(link, math_output.solver_output, math_id);
+            });
+    }
+}
+
 // output base component
 template <std::derived_from<Base> Component, class ComponentContainer, solver_output_type SolverOutputType,
           non_owning_view_c ComponentOutput>

power_grid_model_c/power_grid_model/include/power_grid_model/main_core/topological_node_output.hpp

@@ -22,12 +22,12 @@ namespace detail {
 template <symmetry_tag sym> struct SuperNodeSolverInput {
     std::span<BranchIdx const> links;
     UserNodeValueVector<sym> node_injection;
-    UserNodeValueVector<sym> node_flow_from_branch;
+    UserNodeValueVector<sym> node_flow_through_grid;
 
     ComplexValueVector<sym> get_total_injection_per_node() const {
-        assert(node_injection.size() == node_flow_from_branch.size());
+        assert(node_injection.size() == node_flow_through_grid.size());
 
-        return std::views::zip(node_injection, node_flow_from_branch) | std::views::transform([](auto const& pair) {
+        return std::views::zip(node_injection, node_flow_through_grid) | std::views::transform([](auto const& pair) {
                    auto const& [node_inj, branch_flow] = pair;
                    return ComplexValue<sym>(node_inj + branch_flow);
                }) |
@@ -85,6 +85,9 @@ inline Idx get_node_sequence_idx(main_model_state_c auto const& state, Idx compo
         return state.comp_topo
             ->load_gen_node_idx[get_component_sequence_offset<GenericLoadGen, ComponentType>(state.components) +
                                 component_idx];
+    } else if constexpr (std::derived_from<ComponentType, Shunt>) {
+        return state.comp_topo
+            ->shunt_node_idx[get_component_sequence_offset<Shunt, ComponentType>(state.components) + component_idx];
     } else if constexpr (std::same_as<ComponentType, Fault>) {
         auto const& fault = get_component_by_sequence<Fault>(state.components, component_idx);
         return get_component_sequence_idx<Node>(state.components, fault.get_fault_object());
@@ -99,10 +102,16 @@ inline BranchIdx const& get_branch_sequence_idx(main_model_state_c auto const& s
         ->branch_node_idx[get_component_sequence_offset<Branch, ComponentType>(state.components) + component_idx];
 }
 
+template <std::derived_from<Branch3> ComponentType>
+inline Branch3Idx const& get_branch3_sequence_idx(main_model_state_c auto const& state, Idx component_idx) {
+    return state.comp_topo
+        ->branch3_node_idx[get_component_sequence_offset<Branch3, ComponentType>(state.components) + component_idx];
+}
+
 struct AddApplianceInjection {
     template <typename ComponentType, typename SolverOutputType, typename AddToTarget>
         requires flow_accumulator_c<AddToTarget, ComponentType, SolverOutputType> &&
-                 (is_in_list_c<ComponentType, Source, SymLoad, SymGenerator, AsymLoad, AsymGenerator> ||
+                 (std::derived_from<ComponentType, Appliance> ||
                   (std::same_as<ComponentType, Fault> && short_circuit_solver_output_type<SolverOutputType>))
     void operator()(main_model_state_c auto const& state, MathOutput<std::vector<SolverOutputType>> const& math_output,
                     AddToTarget accumulate_injection) const {
@@ -133,26 +142,66 @@ struct AddApplianceInjection {
             auto const& component_output = get_component_output<ComponentType>(math_output, component_math_id);
 
             auto const& branch_node_idx = get_branch_sequence_idx<ComponentType>(state, component_idx);
-            for (auto const side : {BranchSide::from, BranchSide::to}) {
-                auto const& user_node_idx = branch_node_idx[std::to_underlying(side)];
+            for (auto&& [side, user_node_idx] :
+                 std::views::zip(std::array{BranchSide::from, BranchSide::to}, branch_node_idx)) {
                 auto const& user_topo_id =
                     state.reduced_topology->topo_node_coup.coupling.user_nodes_to_topo_nodes[user_node_idx];
                 accumulate_injection.template operator()<ComponentType>(user_topo_id,
                                                                         get_injection(component_output, side));
             }
         }
     }
+
+    template <typename ComponentType, typename SolverOutputType, typename AddToTarget>
+        requires flow_accumulator_c<AddToTarget, ComponentType, SolverOutputType> &&
+                 std::derived_from<ComponentType, Branch3>
+    void operator()(main_model_state_c auto const& state, MathOutput<std::vector<SolverOutputType>> const& math_output,
+                    AddToTarget accumulate_injection) const {
+        for (auto const& [component_idx, component_math_id] : enumerate(comp_base_sequence<ComponentType>(state))) {
+            if (component_math_id.group == disconnected) {
+                continue;
+            }
+
+            auto const& branch3_node_idx = get_branch3_sequence_idx<ComponentType>(state, component_idx);
+            for (auto&& [side_pos, user_node_idx] : std::views::zip(component_math_id.pos, branch3_node_idx)) {
+                if (side_pos == disconnected) {
+                    continue;
+                }
+                auto const& component_output =
+                    get_component_output<ComponentType>(math_output, {component_math_id.group, side_pos});
+                auto const& user_topo_id =
+                    state.reduced_topology->topo_node_coup.coupling.user_nodes_to_topo_nodes[user_node_idx];
+                accumulate_injection.template operator()<ComponentType>(
+                    user_topo_id, get_injection(component_output, BranchSide::from));
+            }
+        }
+    }
 };
 
 constexpr auto add_appliance_injection = AddApplianceInjection{};
 
-template <typename InjectionComponentTypesTuple, main_model_state_c State, solver_output_type SolverOutput,
-          typename AddToTarget>
+struct user_node_contribution_tag_t {};
+struct ContributesToSteadyStateUserNodeInjection : public user_node_contribution_tag_t {
+    template <typename ComponentType>
+    static constexpr bool value = std::derived_from<ComponentType, Appliance> ||
+                                  std::derived_from<ComponentType, Branch> || std::derived_from<ComponentType, Branch3>;
+};
+struct ContributesToShortCircuitUserNodeInjection : public user_node_contribution_tag_t {
+    template <typename ComponentType>
+    static constexpr bool value = std::derived_from<ComponentType, Source> || std::derived_from<ComponentType, Shunt> ||
+                                  std::derived_from<ComponentType, Branch> ||
+                                  std::derived_from<ComponentType, Branch3> || std::derived_from<ComponentType, Fault>;
+};
+
+template <std::derived_from<user_node_contribution_tag_t> ContributesToUserNodeInjection, main_model_state_c State,
+          solver_output_type SolverOutput, typename AddToTarget>
 inline void add_flows(State const& state, MathOutput<std::vector<SolverOutput>> const& math_output,
                       AddToTarget accumulate_injection) {
-    utils::run_functor_with_tuple_return_void<InjectionComponentTypesTuple>(
+    using Container = decltype(state.components);
+
+    utils::run_functor_with_tuple_return_void<typename Container::storageable_types>(
         [&state, &math_output, &accumulate_injection]<typename ComponentType>() {
-            if constexpr (decltype(state.components)::template is_storageable_v<ComponentType>) {
+            if constexpr (ContributesToUserNodeInjection::template value<ComponentType>) {
                 add_appliance_injection.template operator()<ComponentType>(state, math_output, accumulate_injection);
             }
         });
@@ -227,26 +276,24 @@ solve_topological_nodes(LinkSolver link_solver, State const& state,
             auto const node_number = topo_node.user_nodes.size();
             return {.links = std::span{topo_node.user_links},
                     .node_injection = UserNodeValueVector<sym>(node_number),
-                    .node_flow_from_branch = UserNodeValueVector<sym>(node_number)};
+                    .node_flow_through_grid = UserNodeValueVector<sym>(node_number)};
         }) |
         std::ranges::to<std::vector>();
 
     auto const accumulate_injection = [&link_solver_input]<typename ComponentType>(Idx2D const& user_topo_id,
                                                                                    ComplexValue<sym> const& injection) {
-        if constexpr (std::derived_from<ComponentType, Branch>) {
-            link_solver_input[user_topo_id.group].node_flow_from_branch[user_topo_id.pos] += injection;
+        if constexpr (std::derived_from<ComponentType, Branch> || std::derived_from<ComponentType, Branch3> ||
+                      std::derived_from<ComponentType, Shunt>) {
+            link_solver_input[user_topo_id.group].node_flow_through_grid[user_topo_id.pos] += injection;
         } else {
             link_solver_input[user_topo_id.group].node_injection[user_topo_id.pos] += injection;
         }
     };
 
     if constexpr (steady_state_solver_output_type<SolverOutput>) {
-        using InjectionComponentTypesTuple =
-            std::tuple<Source, SymLoad, SymGenerator, AsymLoad, AsymGenerator, Line, GenericBranch, Transformer>;
-        add_flows<InjectionComponentTypesTuple>(state, math_output, accumulate_injection);
+        add_flows<ContributesToSteadyStateUserNodeInjection>(state, math_output, accumulate_injection);
     } else if constexpr (short_circuit_solver_output_type<SolverOutput>) {
-        using InjectionComponentTypesTuple = std::tuple<Source, Line, GenericBranch, Transformer, Fault>;
-        add_flows<InjectionComponentTypesTuple>(state, math_output, accumulate_injection);
+        add_flows<ContributesToShortCircuitUserNodeInjection>(state, math_output, accumulate_injection);
     }
 
     auto result =