Skip to content
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,6 @@
#include <vector>

namespace power_grid_model::main_core {
constexpr Idx isolated_component{-1};
constexpr Idx not_connected{-1};

namespace detail {
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@
#include "../common/common.hpp"

namespace power_grid_model::main_core {
constexpr Idx isolated_component{-1};

namespace detail {
template <std::derived_from<Branch> ComponentType, typename ComponentContainer>
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@
#include "calculation_parameters.hpp"
#include "container.hpp"
#include "main_model_fwd.hpp"
#include "prepare_calculate.hpp"
#include "topology.hpp"

// common
Expand All @@ -33,11 +34,9 @@
#include "main_core/calculation_input_preparation.hpp"
#include "main_core/input.hpp"
#include "main_core/main_model_type.hpp"
#include "main_core/math_state.hpp"
#include "main_core/output.hpp"
#include "main_core/topology.hpp"
#include "main_core/update.hpp"
#include "main_core/y_bus.hpp"

// stl library
#include <memory>
Expand Down Expand Up @@ -122,7 +121,6 @@ class MainModelImpl {

private:
// internal type traits
using ComponentContainer = typename ModelType::ComponentContainer;
using MainModelState = typename ModelType::MainModelState;

using SequenceIdx = typename ModelType::SequenceIdx;
Expand All @@ -131,22 +129,17 @@ class MainModelImpl {
using OwnedUpdateDataset = typename ModelType::OwnedUpdateDataset;
using ComponentFlags = typename ModelType::ComponentFlags;

static constexpr Idx isolated_component{main_core::isolated_component};
static constexpr Idx not_connected{main_core::not_connected};
static constexpr Idx sequential{main_core::utils::sequential};

public:
using ImplType = ModelType;
using Options = MainModelOptions;
using MathState = main_core::MathState;
using MetaData = meta_data::MetaData;

// constructor with 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} {
solver_preparation_context_{.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();
Expand All @@ -157,7 +150,7 @@ class MainModelImpl {
MathSolverDispatcher const& math_solver_dispatcher)
: system_frequency_{system_frequency},
meta_data_{&meta_data},
math_solver_dispatcher_{&math_solver_dispatcher} {}
solver_preparation_context_{.math_solver_dispatcher = &math_solver_dispatcher} {}

// helper function to get what components are present in the update data
ComponentFlags get_components_to_update(ConstDataset const& update_data) const {
Expand Down Expand Up @@ -206,7 +199,7 @@ class MainModelImpl {

UpdateChange const changed = main_core::update::update_component<CompType>(
state_.components, std::forward<Updates>(updates),
std::back_inserter(std::get<comp_index>(parameter_changed_components_)), sequence_idx);
std::back_inserter(std::get<comp_index>(state_status_context_.parameter_changed_components)), sequence_idx);

// update, get changed variable
update_state(changed);
Expand Down Expand Up @@ -266,18 +259,6 @@ class MainModelImpl {
std::make_shared<ComponentTopology const>(main_core::construct_topology<ModelType>(state_.components));
}

void reset_solvers() {
assert(construction_complete_);
is_topology_up_to_date_ = false;
is_sym_parameter_up_to_date_ = false;
is_asym_parameter_up_to_date_ = false;
n_math_solvers_ = 0;
main_core::clear(math_state_);
state_.math_topology.clear();
state_.topo_comp_coup.reset();
state_.comp_coup = {};
}

public:
/*
the the sequence indexer given an input array of ID's for a given component type
Expand Down Expand Up @@ -307,9 +288,11 @@ class MainModelImpl {
void update_state(UpdateChange const& changes) {
// if topology changed, everything is not up to date
// if only param changed, set param to not up to date
is_topology_up_to_date_ = is_topology_up_to_date_ && !changes.topo;
is_sym_parameter_up_to_date_ = is_sym_parameter_up_to_date_ && !changes.topo && !changes.param;
is_asym_parameter_up_to_date_ = is_asym_parameter_up_to_date_ && !changes.topo && !changes.param;
state_status_context_.is_topology_up_to_date = state_status_context_.is_topology_up_to_date && !changes.topo;
state_status_context_.is_parameter_up_to_date.sym =
state_status_context_.is_parameter_up_to_date.sym && !changes.topo && !changes.param;
state_status_context_.is_parameter_up_to_date.asym =
state_status_context_.is_parameter_up_to_date.asym && !changes.topo && !changes.param;
}

template <typename CompType> void restore_component(SequenceIdxView const& sequence_idx) {
Expand Down Expand Up @@ -363,18 +346,20 @@ class MainModelImpl {
// prepare
auto const& input = [this, &logger, prepare_input_ = std::forward<PrepareInputFn>(prepare_input)] {
Timer const timer{logger, LogEvent::prepare};
prepare_solvers<sym>();
assert(is_topology_up_to_date_ && is_parameter_up_to_date<sym>());
return prepare_input_(n_math_solvers_);
assert(construction_complete_);
prepare_solvers<sym>(state_, solver_preparation_context_, state_status_context_);
assert((state_status_context_.is_topology_up_to_date &&
is_parameter_up_to_date<sym, ImplType>(state_status_context_.is_parameter_up_to_date)));
return prepare_input_(solver_preparation_context_.n_math_solvers);
}();
// calculate
return [this, &logger, &input, solve_ = std::forward<SolveFn>(solve)] {
Timer const timer{logger, LogEvent::math_calculation};
auto& solvers = main_core::get_solvers<sym>(math_state_);
auto& y_bus_vec = main_core::get_y_bus<sym>(math_state_);
auto& solvers = main_core::get_solvers<sym>(solver_preparation_context_.math_state);
auto& y_bus_vec = main_core::get_y_bus<sym>(solver_preparation_context_.math_state);
std::vector<SolverOutputType> solver_output;
solver_output.reserve(n_math_solvers_);
for (Idx i = 0; i != n_math_solvers_; ++i) {
solver_output.reserve(solver_preparation_context_.n_math_solvers);
for (Idx i = 0; i != solver_preparation_context_.n_math_solvers; ++i) {
solver_output.emplace_back(solve_(solvers[i], y_bus_vec[i], input[i]));
}
return solver_output;
Expand Down Expand Up @@ -424,7 +409,8 @@ class MainModelImpl {

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>());
assert((state_status_context_.is_topology_up_to_date &&
is_parameter_up_to_date<sym, ImplType>(state_status_context_.is_parameter_up_to_date)));
return main_core::prepare_short_circuit_input<sym>(state_, state_.comp_coup, n_math_solvers,
voltage_scaling);
},
Expand Down Expand Up @@ -542,89 +528,27 @@ class MainModelImpl {

double system_frequency_;
MetaData const* meta_data_;
MathSolverDispatcher const* math_solver_dispatcher_;

MainModelState state_;
// math model
MathState math_state_;
Idx n_math_solvers_{0};
bool is_topology_up_to_date_{false};
bool is_sym_parameter_up_to_date_{false};
bool is_asym_parameter_up_to_date_{false};
bool is_accumulated_component_updated_{true};
bool last_updated_calculation_symmetry_mode_{false};

SolverPreparationContext solver_preparation_context_;
// math_state
// n_math_solvers
// math_solver_dispatcher

StatusCheckingContext<ImplType> state_status_context_{};
// is_topology_up_to_date
// is_sym_parameter_up_to_date
// is_asym_parameter_up_to_date
// last_updated_calculation_symmetry_mode
// parameter_changed_components

OwnedUpdateDataset cached_inverse_update_{};
UpdateChange cached_state_changes_{};
SequenceIdx parameter_changed_components_{};
#ifndef NDEBUG
// construction_complete is used for debug assertions only
bool construction_complete_{false};
#endif // !NDEBUG

template <symmetry_tag sym> bool& is_parameter_up_to_date() {
if constexpr (is_symmetric_v<sym>) {
return is_sym_parameter_up_to_date_;
} else {
return is_asym_parameter_up_to_date_;
}
}

void rebuild_topology() {
assert(construction_complete_);
// clear old solvers
reset_solvers();
ComponentConnections const comp_conn =
main_core::construct_components_connections<ModelType>(state_.components);
// re build
Topology topology{*state_.comp_topo, comp_conn};
std::tie(state_.math_topology, state_.topo_comp_coup) = topology.build_topology();
n_math_solvers_ = static_cast<Idx>(state_.math_topology.size());
is_topology_up_to_date_ = true;
is_sym_parameter_up_to_date_ = false;
is_asym_parameter_up_to_date_ = false;
}

template <symmetry_tag sym> void prepare_solvers() {
std::vector<MathSolverProxy<sym>>& solvers = main_core::get_solvers<sym>(math_state_);
// rebuild topology if needed
if (!is_topology_up_to_date_) {
rebuild_topology();
}
main_core::prepare_y_bus<sym, ModelType>(state_, n_math_solvers_, math_state_);

if (n_math_solvers_ != static_cast<Idx>(solvers.size())) {
assert(solvers.empty());
assert(n_math_solvers_ == static_cast<Idx>(state_.math_topology.size()));
assert(n_math_solvers_ == static_cast<Idx>(main_core::get_y_bus<sym>(math_state_).size()));

solvers.clear();
solvers.reserve(n_math_solvers_);
std::ranges::transform(state_.math_topology, std::back_inserter(solvers), [this](auto const& math_topo) {
return MathSolverProxy<sym>{math_solver_dispatcher_, math_topo};
});
for (Idx idx = 0; idx < n_math_solvers_; ++idx) {
main_core::get_y_bus<sym>(math_state_)[idx].register_parameters_changed_callback(
[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 =
main_core::get_math_param<sym>(state_, n_math_solvers_);
std::vector<MathModelParamIncrement> const math_param_increments =
main_core::get_math_param_increment<ModelType>(state_, n_math_solvers_, parameter_changed_components_);
if (last_updated_calculation_symmetry_mode_ == is_symmetric_v<sym>) {
main_core::update_y_bus(math_state_, math_params, math_param_increments);
} else {
main_core::update_y_bus(math_state_, math_params);
}
}
// else do nothing, set everything up to date
is_parameter_up_to_date<sym>() = true;
std::ranges::for_each(parameter_changed_components_, [](auto& comps) { comps.clear(); });
last_updated_calculation_symmetry_mode_ = is_symmetric_v<sym>;
}
};

} // namespace power_grid_model
Original file line number Diff line number Diff line change
@@ -0,0 +1,120 @@
// SPDX-FileCopyrightText: Contributors to the Power Grid Model project <[email protected]>
//
// SPDX-License-Identifier: MPL-2.0

#pragma once

#include "topology.hpp"

#include "common/common.hpp"

#include "math_solver/math_solver_dispatch.hpp"

#include "main_core/main_model_type.hpp"
#include "main_core/math_state.hpp"
#include "main_core/topology.hpp"
#include "main_core/y_bus.hpp"

namespace power_grid_model {
struct SolverPreparationContext {
main_core::MathState math_state;
Idx n_math_solvers{0};
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is it required to keep n_math_solvers?
Is it ever not equivalent to math_state.topology.size()?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I removed it (in 4f0d589) in favor of a helper function that does the calling from state.math_topology.size() indeed. It should be fine and

assert(n_math_solvers == static_cast<Idx>(main_core::get_y_bus<sym>(solver_context.math_state).size()));|

still remains as a sanity check.

Let me know what you think.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i think this is indeed legacy from the era when we would re-use a math solver when our amount of connected components in the grid remained the same, even when our topology changed (i.e. we would reuse preallocated math solver memory, but it doesn't make sense because the elephant in the room is not the math solver but the y bus structure/topology. therefore, we changed that)

as long as our cache invalidation works correctly, with the current state of the code base, that should work as intended

MathSolverDispatcher const* math_solver_dispatcher;
};

template <class ModelType>
requires(main_core::is_main_model_type_v<ModelType>)
struct StatusCheckingContext {
bool is_topology_up_to_date{false};
bool last_updated_calculation_symmetry_mode{false};
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i don't like this name with the fact that it's a bool. Maybe you can make it an enum? or otherwise rename to something along the lines of last_updated_was_symmetric?

this is basically an invariant, namely a mapping of the previous calculation type to the next. maybe you can solve it using a setter (instead of doing the checking at the user side), e.g.

class StatusCheckingContext {
    template <symmetry_tag sym> set_most_recent_calculation_symmetry() {
        last_updated_calculation_was_symmetric_ = is_symmetric_v<sym>;
    }
    template <symmetry_tag sym> is_same_calculation_symmetry() {
        return last_updated_calculation_was_symmetric_ == is_symmetric_v<sym>;
    }
};

You can also use IsParameterUpToDateHelper for caching (although slightly less memory efficient, but if we really want to go down that direction, we also don't need a separate bool for topology; we could just use bitflags or enums

note that using std::variant would be overkill.

typename ModelType::SequenceIdx parameter_changed_components{};
Copy link
Member

@mgovers mgovers Nov 17, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

it feels a little bit confusing to me that the status and the cache are part of the same struct called StatusCheckingContext.

Maybe the name is just not clear enough. Maybe something along the lines of CacheStatus or CacheValidity?

That said, maybe the confusing part is not the name of the type but the fact that one thing tells something about the cache validity while the other about the context in which it is valid? maybe that means that they shouldn't be part of the same class in the first place?

struct IsParameterUpToDateHelper {
bool sym{false};
bool asym{false};
} is_parameter_up_to_date{};
};

namespace detail {
template <class ModelType>
void reset_solvers(typename ModelType::MainModelState& state, SolverPreparationContext& solver_context,
StatusCheckingContext<ModelType>& status_context) {
status_context.is_topology_up_to_date = false;
status_context.is_parameter_up_to_date.sym = false;
status_context.is_parameter_up_to_date.asym = false;
solver_context.n_math_solvers = 0;
main_core::clear(solver_context.math_state);
state.math_topology.clear();
state.topo_comp_coup.reset();
state.comp_coup = {};
}

template <class ModelType>
void rebuild_topology(typename ModelType::MainModelState& state, SolverPreparationContext& solver_context,
StatusCheckingContext<ModelType>& status_context) {
// clear old solvers
reset_solvers(state, solver_context, status_context);
ComponentConnections const comp_conn = main_core::construct_components_connections<ModelType>(state.components);
// re build
Topology topology{*state.comp_topo, comp_conn};
std::tie(state.math_topology, state.topo_comp_coup) = topology.build_topology();
solver_context.n_math_solvers = static_cast<Idx>(state.math_topology.size());
status_context.is_topology_up_to_date = true;
status_context.is_parameter_up_to_date.sym = false;
status_context.is_parameter_up_to_date.asym = false;
}
} // namespace detail

template <symmetry_tag sym, class ModelType>
bool& is_parameter_up_to_date(
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i prefer to make this a member function

typename StatusCheckingContext<ModelType>::IsParameterUpToDateHelper& is_parameter_up_to_date) {
if constexpr (is_symmetric_v<sym>) {
return is_parameter_up_to_date.sym;
} else {
return is_parameter_up_to_date.asym;
}
}

template <symmetry_tag sym, class ModelType>
void prepare_solvers(typename ModelType::MainModelState& state, SolverPreparationContext& solver_context,
StatusCheckingContext<ModelType>& status_context) {
std::vector<MathSolverProxy<sym>>& solvers = main_core::get_solvers<sym>(solver_context.math_state);
// rebuild topology if needed
if (!status_context.is_topology_up_to_date) {
detail::rebuild_topology(state, solver_context, status_context);
}
main_core::prepare_y_bus<sym, ModelType>(state, solver_context.n_math_solvers, solver_context.math_state);

if (solver_context.n_math_solvers != static_cast<Idx>(solvers.size())) {
assert(solvers.empty());
assert(solver_context.n_math_solvers == static_cast<Idx>(state.math_topology.size()));
assert(solver_context.n_math_solvers ==
static_cast<Idx>(main_core::get_y_bus<sym>(solver_context.math_state).size()));

solvers.clear();
solvers.reserve(solver_context.n_math_solvers);
std::ranges::transform(state.math_topology, std::back_inserter(solvers),
[&solver_context](auto const& math_topo) {
return MathSolverProxy<sym>{solver_context.math_solver_dispatcher, math_topo};
});
for (Idx idx = 0; idx < solver_context.n_math_solvers; ++idx) {
main_core::get_y_bus<sym>(solver_context.math_state)[idx].register_parameters_changed_callback(
[solver = std::ref(solvers[idx])](bool changed) { solver.get().get().parameters_changed(changed); });
}
} else if (!is_parameter_up_to_date<sym, ModelType>(status_context.is_parameter_up_to_date)) {
std::vector<MathModelParam<sym>> const math_params =
main_core::get_math_param<sym>(state, solver_context.n_math_solvers);
std::vector<MathModelParamIncrement> const math_param_increments =
main_core::get_math_param_increment<ModelType>(state, solver_context.n_math_solvers,
status_context.parameter_changed_components);
if (status_context.last_updated_calculation_symmetry_mode == is_symmetric_v<sym>) {
main_core::update_y_bus(solver_context.math_state, math_params, math_param_increments);
} else {
main_core::update_y_bus(solver_context.math_state, math_params);
}
}
// else do nothing, set everything up to date
is_parameter_up_to_date<sym, ModelType>(status_context.is_parameter_up_to_date) = true;
std::ranges::for_each(status_context.parameter_changed_components, [](auto& comps) { comps.clear(); });
status_context.last_updated_calculation_symmetry_mode = is_symmetric_v<sym>;
}
} // namespace power_grid_model
Loading