677 lines
30 KiB
C++
677 lines
30 KiB
C++
module;
|
|
|
|
#include <algorithm>
|
|
#include <chrono>
|
|
#include <concepts>
|
|
#include <cstddef>
|
|
#include <cstdint>
|
|
#include <memory>
|
|
#include <stdexcept>
|
|
#include <tuple>
|
|
#include <type_traits>
|
|
#include <utility>
|
|
|
|
#include <mfem.hpp>
|
|
|
|
export module mean_field:preconditioning.stellar_equilibrium;
|
|
|
|
export import :operators.stellar_equilibrium_problem;
|
|
export import :preconditioning.plan;
|
|
|
|
export namespace mean_field::preconditioning {
|
|
struct StellarPreconditionerLifecycleSnapshot final {
|
|
operators::StellarEquilibriumDependencyStamp discretization;
|
|
operators::StellarEquilibriumDependencyStamp geometry;
|
|
const void *equationOfStateIdentity{nullptr};
|
|
operators::StellarEquilibriumDependencies linearization;
|
|
std::uint64_t preparedOperatorGeneration{0};
|
|
|
|
constexpr bool operator==(const StellarPreconditionerLifecycleSnapshot &) const = default;
|
|
};
|
|
|
|
struct StellarPreconditionerPreparationChanges final {
|
|
bool discretization{false};
|
|
bool geometry{false};
|
|
bool equationOfState{false};
|
|
bool linearization{false};
|
|
|
|
[[nodiscard]] constexpr bool Any() const noexcept {
|
|
return discretization || geometry || equationOfState || linearization;
|
|
}
|
|
|
|
[[nodiscard]] constexpr bool Contains(const PreparationDependency dependency) const noexcept {
|
|
switch (dependency) {
|
|
case PreparationDependency::discretization:
|
|
return discretization;
|
|
case PreparationDependency::geometry:
|
|
return geometry;
|
|
case PreparationDependency::equation_of_state:
|
|
return equationOfState;
|
|
case PreparationDependency::linearization:
|
|
return linearization;
|
|
}
|
|
return false;
|
|
}
|
|
};
|
|
|
|
[[nodiscard]] constexpr StellarPreconditionerPreparationChanges preparationChanges(
|
|
const StellarPreconditionerLifecycleSnapshot &prepared,
|
|
const StellarPreconditionerLifecycleSnapshot ¤t
|
|
) noexcept {
|
|
return {
|
|
.discretization = prepared.discretization != current.discretization,
|
|
.geometry = prepared.geometry != current.geometry,
|
|
.equationOfState = prepared.equationOfStateIdentity != current.equationOfStateIdentity,
|
|
.linearization = prepared.linearization != current.linearization ||
|
|
prepared.preparedOperatorGeneration != current.preparedOperatorGeneration
|
|
};
|
|
}
|
|
|
|
struct StellarPreconditionerPreparationReport final {
|
|
StellarPreconditionerPreparationChanges changes;
|
|
std::uint64_t refreshedComponents{0};
|
|
|
|
[[nodiscard]] constexpr bool DidAnyWork() const noexcept {
|
|
return refreshedComponents != 0;
|
|
}
|
|
};
|
|
|
|
struct StellarPreconditionerStatistics final {
|
|
std::uint64_t setups{0};
|
|
std::uint64_t refreshChecks{0};
|
|
std::uint64_t refreshes{0};
|
|
std::uint64_t noOpRefreshes{0};
|
|
std::uint64_t componentSetups{0};
|
|
std::uint64_t componentRefreshes{0};
|
|
std::uint64_t operatorBindings{0};
|
|
std::uint64_t applications{0};
|
|
std::uint64_t backendApplications{0};
|
|
std::uint64_t innerIterations{0};
|
|
double setupSeconds{0.0};
|
|
double refreshSeconds{0.0};
|
|
double applicationSeconds{0.0};
|
|
double maximumApplicationSeconds{0.0};
|
|
};
|
|
|
|
template <typename Candidate> struct StellarEquilibriumProblemTraits {
|
|
static constexpr bool registered = false;
|
|
};
|
|
|
|
template <equilibrium::StellarEquilibriumModel Model, equilibrium::StellarDiscretizationType Discretization>
|
|
struct StellarEquilibriumProblemTraits<equilibrium::StellarEquilibriumProblem<Model, Discretization>> {
|
|
using Problem = equilibrium::StellarEquilibriumProblem<Model, Discretization>;
|
|
using Form = typename Problem::FormType;
|
|
using JacobianForm = typename Problem::JacobianFormType;
|
|
using Manifest = typename Problem::ManifestType;
|
|
|
|
static constexpr bool registered = true;
|
|
|
|
[[nodiscard]] static bool IsPrepared(const Problem &problem) noexcept {
|
|
return problem.IsPrepared();
|
|
}
|
|
|
|
[[nodiscard]] static int StateSize(const Problem &problem) noexcept {
|
|
return problem.StateSize();
|
|
}
|
|
|
|
[[nodiscard]] static int EquationSize(const Problem &problem) noexcept {
|
|
return problem.EquationSize();
|
|
}
|
|
|
|
[[nodiscard]] static const Manifest &ManifestOf(const Problem &problem) noexcept {
|
|
return problem.GetManifest();
|
|
}
|
|
|
|
[[nodiscard]] static const mfem::Operator &LinearizationOperator(const Problem &problem) noexcept {
|
|
return problem.GetLinearizationOperator();
|
|
}
|
|
|
|
[[nodiscard]] static StellarPreconditionerLifecycleSnapshot Snapshot(const Problem &problem) {
|
|
const operators::StellarEquilibriumDependencies &dependencies = problem.GetLinearizationDependencies();
|
|
return {
|
|
.discretization = dependencies.discretization,
|
|
.geometry = problem.GetGeometryDependency(),
|
|
.equationOfStateIdentity = std::addressof(problem.GetStellarModel().equationOfState()),
|
|
.linearization = dependencies,
|
|
.preparedOperatorGeneration = problem.GetPreparationGeneration()
|
|
};
|
|
}
|
|
};
|
|
|
|
template <typename Candidate>
|
|
concept StellarPreconditionerProblem = StellarEquilibriumProblemTraits<std::remove_cvref_t<Candidate>>::registered;
|
|
|
|
namespace detail {
|
|
template <typename Block> struct IsGeneratedStellarValueBlock : std::false_type { };
|
|
|
|
template <typename Generated>
|
|
struct IsGeneratedStellarValueBlock<utils::blocks::generated_value_block<Generated>> : std::true_type { };
|
|
|
|
template <typename Block> struct IsGeneratedStellarResidualBlock : std::false_type { };
|
|
|
|
template <typename Generated>
|
|
struct IsGeneratedStellarResidualBlock<utils::blocks::generated_residual_block<Generated>> : std::true_type { };
|
|
|
|
template <typename Coupling>
|
|
inline constexpr bool isPurePhysicalStellarCoupling =
|
|
!IsGeneratedStellarValueBlock<std::remove_cvref_t<typename Coupling::Value>>::value &&
|
|
!IsGeneratedStellarResidualBlock<std::remove_cvref_t<typename Coupling::Residual>>::value;
|
|
|
|
/* Pure structure contributions owned by a trusted backend are exact
|
|
* (core, specification, coupling) capabilities. Future cores and new
|
|
* edges start with no privilege: changing a built-in declaration must
|
|
* be accompanied by an explicit preconditioner decision. Generated-
|
|
* border terms remain the responsibility of specification-border
|
|
* machinery. */
|
|
template <typename PhysicalCore, typename Specification> struct StellarStructureBackendHandledCouplings {
|
|
using Type = utils::blocks::type_list<>;
|
|
};
|
|
|
|
template <>
|
|
struct StellarStructureBackendHandledCouplings<
|
|
operators::PreparedStellarEquilibriumOperator,
|
|
models::FixedTotalMass> {
|
|
using Type = utils::blocks::type_list<
|
|
operators::StellarEquilibriumJacobianCoupling<
|
|
utils::blocks::enthalpy::specific::residual,
|
|
utils::blocks::density::mass::value>,
|
|
operators::StellarEquilibriumJacobianCoupling<
|
|
utils::blocks::enthalpy::specific::residual,
|
|
utils::blocks::surface_deformation::parameters::value>>;
|
|
};
|
|
|
|
template <>
|
|
struct StellarStructureBackendHandledCouplings<
|
|
operators::PreparedStellarEquilibriumOperator,
|
|
models::FixedAngularMomentum> {
|
|
using Type = utils::blocks::type_list<
|
|
operators::StellarEquilibriumJacobianCoupling<
|
|
utils::blocks::surface_deformation::shape_equilibrium::residual,
|
|
utils::blocks::density::mass::value>,
|
|
operators::StellarEquilibriumJacobianCoupling<
|
|
utils::blocks::surface_deformation::shape_equilibrium::residual,
|
|
utils::blocks::surface_deformation::parameters::value>,
|
|
operators::StellarEquilibriumJacobianCoupling<
|
|
utils::blocks::enthalpy::specific::residual,
|
|
utils::blocks::density::mass::value>,
|
|
operators::StellarEquilibriumJacobianCoupling<
|
|
utils::blocks::enthalpy::specific::residual,
|
|
utils::blocks::surface_deformation::parameters::value>>;
|
|
};
|
|
|
|
template <>
|
|
struct StellarStructureBackendHandledCouplings<
|
|
operators::PreparedStellarEquilibriumOperator,
|
|
models::FixedCentralDensity> {
|
|
using Type = utils::blocks::type_list<operators::StellarEquilibriumJacobianCoupling<
|
|
utils::blocks::enthalpy::specific::residual,
|
|
utils::blocks::enthalpy::specific::value>>;
|
|
};
|
|
|
|
template <typename Coupling, typename Model, typename PhysicalCore, typename ModelSpecifications>
|
|
struct EveryCouplingContributionHandled;
|
|
|
|
template <
|
|
typename Coupling,
|
|
model::StellarModelType Model,
|
|
typename PhysicalCore,
|
|
models::ModelSpecification... Specifications>
|
|
struct EveryCouplingContributionHandled<
|
|
Coupling,
|
|
Model,
|
|
PhysicalCore,
|
|
models::detail::SpecificationSetStorage<Specifications...>>
|
|
final {
|
|
private:
|
|
template <typename Specification>
|
|
static constexpr bool handled =
|
|
!utils::blocks::contains_type_v<
|
|
Coupling,
|
|
typename operators::StellarEquilibriumSpecificationCompilation<Specification>::JacobianCouplings> ||
|
|
(operators::stellarEquilibriumBackendRuntimeAuthorized<Specification, Model> &&
|
|
utils::blocks::contains_type_v<
|
|
Coupling,
|
|
typename StellarStructureBackendHandledCouplings<PhysicalCore, Specification>::Type>) ||
|
|
operators::stellarEquilibriumSpecificationCouplingIsStructuralZero<Specification, Model, Coupling>;
|
|
|
|
public:
|
|
static constexpr bool value = (handled<Specifications> && ...);
|
|
};
|
|
|
|
template <
|
|
typename Remaining,
|
|
typename Model,
|
|
typename PhysicalCore,
|
|
typename ModelSpecifications,
|
|
typename Unsupported>
|
|
struct CollectUnsupportedStellarStructureCouplings;
|
|
|
|
template <typename Model, typename PhysicalCore, typename ModelSpecifications, typename Unsupported>
|
|
struct CollectUnsupportedStellarStructureCouplings<
|
|
utils::blocks::type_list<>,
|
|
Model,
|
|
PhysicalCore,
|
|
ModelSpecifications,
|
|
Unsupported> {
|
|
using Type = Unsupported;
|
|
};
|
|
|
|
template <
|
|
typename Head,
|
|
typename... Tail,
|
|
typename Model,
|
|
typename PhysicalCore,
|
|
typename ModelSpecifications,
|
|
typename... Unsupported>
|
|
struct CollectUnsupportedStellarStructureCouplings<
|
|
utils::blocks::type_list<Head, Tail...>,
|
|
Model,
|
|
PhysicalCore,
|
|
ModelSpecifications,
|
|
utils::blocks::type_list<Unsupported...>> {
|
|
private:
|
|
static constexpr bool supported =
|
|
!isPurePhysicalStellarCoupling<Head> ||
|
|
EveryCouplingContributionHandled<Head, Model, PhysicalCore, ModelSpecifications>::value;
|
|
using Next = std::conditional_t<
|
|
supported,
|
|
utils::blocks::type_list<Unsupported...>,
|
|
utils::blocks::type_list<Unsupported..., Head>>;
|
|
|
|
public:
|
|
using Type = typename CollectUnsupportedStellarStructureCouplings<
|
|
utils::blocks::type_list<Tail...>,
|
|
Model,
|
|
PhysicalCore,
|
|
ModelSpecifications,
|
|
Next>::Type;
|
|
};
|
|
|
|
template <typename Candidate, typename = void> struct DefaultStellarStructurePhysicalTopologyAudit {
|
|
using ContributionCouplings = utils::blocks::type_list<>;
|
|
using UnsupportedCouplings = utils::blocks::type_list<>;
|
|
|
|
static constexpr bool supported = false;
|
|
};
|
|
|
|
template <model::StellarModelType Model>
|
|
requires(
|
|
operators::StellarEquilibriumSystemCompilable<std::remove_cvref_t<Model>> &&
|
|
operators::hasStellarEquilibriumCoreRuntime<std::remove_cvref_t<Model>>
|
|
)
|
|
struct DefaultStellarStructurePhysicalTopologyAudit<
|
|
Model,
|
|
std::void_t<
|
|
typename operators::CompiledStellarEquilibriumSystem<
|
|
std::remove_cvref_t<Model>>::ContributionJacobianCouplings,
|
|
operators::StellarEquilibriumPhysicalCoreType<std::remove_cvref_t<Model>>>> {
|
|
private:
|
|
using Compilation = operators::CompiledStellarEquilibriumSystem<std::remove_cvref_t<Model>>;
|
|
using PhysicalCore = operators::StellarEquilibriumPhysicalCoreType<std::remove_cvref_t<Model>>;
|
|
|
|
public:
|
|
using ContributionCouplings = typename Compilation::ContributionJacobianCouplings;
|
|
using UnsupportedCouplings = typename CollectUnsupportedStellarStructureCouplings<
|
|
ContributionCouplings,
|
|
std::remove_cvref_t<Model>,
|
|
PhysicalCore,
|
|
typename std::remove_cvref_t<Model>::SpecificationTypes,
|
|
utils::blocks::type_list<>>::Type;
|
|
|
|
static constexpr bool supported = UnsupportedCouplings::size == 0;
|
|
};
|
|
} // namespace detail
|
|
|
|
/* A generated-border edge is owned by specification-border machinery and
|
|
* is deliberately ignored here. Every pure physical edge contributed by a
|
|
* model must be owned by the selected numerical structure backend, or
|
|
* every non-backend provider of that edge must prove StructuralZero. In
|
|
* particular, merely overlapping an existing base-Jacobian edge is not
|
|
* sufficient: a custom nonzero coefficient on that edge would otherwise
|
|
* disappear silently from the default preconditioner. This audit is
|
|
* detection-safe and therefore suitable for constraining factories. */
|
|
template <typename Candidate>
|
|
struct DefaultStellarStructurePhysicalTopologySupport
|
|
: detail::DefaultStellarStructurePhysicalTopologyAudit<std::remove_cvref_t<Candidate>> { };
|
|
|
|
template <typename Candidate>
|
|
inline constexpr bool defaultStellarStructurePhysicalTopologySupported =
|
|
DefaultStellarStructurePhysicalTopologySupport<std::remove_cvref_t<Candidate>>::supported;
|
|
|
|
template <typename Candidate>
|
|
concept DefaultStellarStructurePhysicalTopologySupportedFor =
|
|
defaultStellarStructurePhysicalTopologySupported<Candidate>;
|
|
|
|
namespace backend {
|
|
template <typename Component, typename Problem, typename Backend = typename Component::BackendType>
|
|
class PreparedComponent;
|
|
|
|
template <typename Component, StellarPreconditionerProblem Problem>
|
|
class PreparedComponent<Component, Problem, Identity> final {
|
|
public:
|
|
void Setup(
|
|
const Problem &,
|
|
const Component &
|
|
) noexcept {
|
|
}
|
|
|
|
[[nodiscard]] bool Refresh(
|
|
const Problem &,
|
|
const Component &,
|
|
const StellarPreconditionerPreparationChanges &
|
|
) noexcept {
|
|
return true;
|
|
}
|
|
};
|
|
|
|
template <typename Component, typename Problem>
|
|
concept PreparedComponentFor = requires(
|
|
PreparedComponent<Component, Problem> &prepared,
|
|
const Problem &problem,
|
|
const Component &component,
|
|
const StellarPreconditionerPreparationChanges &changes
|
|
) {
|
|
prepared.Setup(problem, component);
|
|
{ prepared.Refresh(problem, component, changes) } -> std::same_as<bool>;
|
|
};
|
|
} // namespace backend
|
|
|
|
namespace detail {
|
|
template <typename Form> struct IdentityPlanForForm;
|
|
|
|
template <typename... Values, typename... Residuals>
|
|
struct IdentityPlanForForm<
|
|
utils::blocks::block_form<utils::blocks::type_list<Values...>, utils::blocks::type_list<Residuals...>>> {
|
|
static_assert(sizeof...(Values) == sizeof...(Residuals));
|
|
using Type = PreconditionerPlan<IdentityBlock<Values, Residuals>...>;
|
|
|
|
[[nodiscard]] static constexpr Type Make() {
|
|
return Type{IdentityBlock<Values, Residuals>{}...};
|
|
}
|
|
};
|
|
|
|
template <typename ComponentList> struct UsesOnlyIdentityBackends : std::false_type { };
|
|
|
|
template <typename... Components>
|
|
struct UsesOnlyIdentityBackends<utils::blocks::type_list<Components...>>
|
|
: std::bool_constant<(std::same_as<typename Components::BackendType, backend::Identity> && ...)> { };
|
|
|
|
template <typename ComponentList, typename Problem> struct PreparedComponentTuple;
|
|
|
|
template <typename... Components, typename Problem>
|
|
struct PreparedComponentTuple<utils::blocks::type_list<Components...>, Problem> {
|
|
using Type = std::tuple<backend::PreparedComponent<Components, Problem>...>;
|
|
|
|
static constexpr bool available = (backend::PreparedComponentFor<Components, Problem> && ...);
|
|
};
|
|
|
|
template <typename Requirements>
|
|
[[nodiscard]] constexpr bool
|
|
componentRequiresRefresh(const StellarPreconditionerPreparationChanges &changes) noexcept {
|
|
return (Requirements::contains(PreparationDependency::discretization) && changes.discretization) ||
|
|
(Requirements::contains(PreparationDependency::geometry) && changes.geometry) ||
|
|
(Requirements::contains(PreparationDependency::equation_of_state) && changes.equationOfState) ||
|
|
(Requirements::contains(PreparationDependency::linearization) && changes.linearization);
|
|
}
|
|
} // namespace detail
|
|
|
|
template <typename Plan, typename Problem>
|
|
concept PreparedPreconditionerPlanFor =
|
|
StellarPreconditionerProblem<Problem> && PreconditionerPlanType<Plan> &&
|
|
CompletePreconditionerFor<Plan, typename StellarEquilibriumProblemTraits<Problem>::Form> &&
|
|
CompatiblePreconditionerFor<
|
|
Plan,
|
|
typename StellarEquilibriumProblemTraits<Problem>::Form,
|
|
typename StellarEquilibriumProblemTraits<Problem>::JacobianForm> &&
|
|
(!std::remove_cvref_t<Plan>::allowsOverlappingOwnership) &&
|
|
detail::UsesOnlyIdentityBackends<typename std::remove_cvref_t<Plan>::ComponentTypes>::value &&
|
|
detail::PreparedComponentTuple<
|
|
typename std::remove_cvref_t<Plan>::ComponentTypes,
|
|
std::remove_cvref_t<Problem>>::available;
|
|
|
|
template <StellarPreconditionerProblem Problem>
|
|
using IdentityPreconditionerPlanFor = typename detail::IdentityPlanForForm<
|
|
typename StellarEquilibriumProblemTraits<std::remove_cvref_t<Problem>>::Form>::Type;
|
|
|
|
template <StellarPreconditionerProblem Problem>
|
|
[[nodiscard]] constexpr IdentityPreconditionerPlanFor<Problem> makeIdentityPlan(const Problem &) {
|
|
using Form = typename StellarEquilibriumProblemTraits<std::remove_cvref_t<Problem>>::Form;
|
|
return detail::IdentityPlanForForm<Form>::Make();
|
|
}
|
|
|
|
template <StellarPreconditionerProblem Problem, typename Plan>
|
|
requires PreparedPreconditionerPlanFor<Plan, Problem>
|
|
class StellarEquilibriumPreconditioner final : public mfem::Solver {
|
|
private:
|
|
using ProblemType = std::remove_cvref_t<Problem>;
|
|
using PlanType = std::remove_cvref_t<Plan>;
|
|
using Traits = StellarEquilibriumProblemTraits<ProblemType>;
|
|
using Components = typename PlanType::ComponentTypes;
|
|
using PreparedComponents = typename detail::PreparedComponentTuple<Components, ProblemType>::Type;
|
|
using Clock = std::chrono::steady_clock;
|
|
|
|
public:
|
|
using FormType = typename Traits::Form;
|
|
using JacobianFormType = typename Traits::JacobianForm;
|
|
|
|
StellarEquilibriumPreconditioner(
|
|
ProblemType &problem,
|
|
PlanType plan
|
|
)
|
|
: mfem::Solver(Traits::StateSize(problem)),
|
|
m_problem(std::addressof(problem)),
|
|
m_manifest(std::addressof(Traits::ManifestOf(problem))),
|
|
m_linearization(std::addressof(Traits::LinearizationOperator(problem))),
|
|
m_plan(std::move(plan)) {
|
|
const Clock::time_point start = Clock::now();
|
|
VerifyPreparedProblem();
|
|
SetupComponents(std::make_index_sequence<std::tuple_size_v<PreparedComponents>>{});
|
|
m_snapshot = Traits::Snapshot(*m_problem);
|
|
m_statistics.setups = 1;
|
|
m_statistics.setupSeconds = std::chrono::duration<double>(Clock::now() - start).count();
|
|
}
|
|
|
|
StellarEquilibriumPreconditioner(const StellarEquilibriumPreconditioner &) = delete;
|
|
StellarEquilibriumPreconditioner &operator=(const StellarEquilibriumPreconditioner &) = delete;
|
|
StellarEquilibriumPreconditioner(StellarEquilibriumPreconditioner &&) = delete;
|
|
StellarEquilibriumPreconditioner &operator=(StellarEquilibriumPreconditioner &&) = delete;
|
|
|
|
void SetOperator(const mfem::Operator &operation) override {
|
|
if (operation.Height() != Height() || operation.Width() != Width()) {
|
|
throw std::invalid_argument(
|
|
"The stellar-equilibrium preconditioner received an operator with incompatible dimensions."
|
|
);
|
|
}
|
|
++m_statistics.operatorBindings;
|
|
}
|
|
|
|
void Mult(
|
|
const mfem::Vector &residual,
|
|
mfem::Vector &correction
|
|
) const override {
|
|
VerifyCurrent();
|
|
if (residual.Size() != Width()) {
|
|
throw std::invalid_argument(
|
|
"The stellar-equilibrium preconditioner received a residual with the wrong size."
|
|
);
|
|
}
|
|
if (correction.Size() != Height()) {
|
|
throw std::invalid_argument(
|
|
"The stellar-equilibrium preconditioner requires a preallocated correction of the correct size."
|
|
);
|
|
}
|
|
|
|
const Clock::time_point start = Clock::now();
|
|
correction = residual;
|
|
const double elapsed = std::chrono::duration<double>(Clock::now() - start).count();
|
|
|
|
++m_statistics.applications;
|
|
++m_statistics.backendApplications;
|
|
m_statistics.applicationSeconds += elapsed;
|
|
m_statistics.maximumApplicationSeconds = std::max(m_statistics.maximumApplicationSeconds, elapsed);
|
|
}
|
|
|
|
[[nodiscard]] StellarPreconditionerPreparationReport Refresh() {
|
|
const Clock::time_point start = Clock::now();
|
|
VerifyPreparedProblem();
|
|
|
|
const StellarPreconditionerLifecycleSnapshot current = Traits::Snapshot(*m_problem);
|
|
const StellarPreconditionerPreparationChanges changes = preparationChanges(m_snapshot, current);
|
|
++m_statistics.refreshChecks;
|
|
|
|
StellarPreconditionerPreparationReport report{.changes = changes};
|
|
if (!changes.Any()) {
|
|
++m_statistics.noOpRefreshes;
|
|
} else {
|
|
report.refreshedComponents =
|
|
RefreshComponents(changes, std::make_index_sequence<std::tuple_size_v<PreparedComponents>>{});
|
|
++m_statistics.refreshes;
|
|
m_statistics.componentRefreshes += report.refreshedComponents;
|
|
m_snapshot = current;
|
|
}
|
|
|
|
m_statistics.refreshSeconds += std::chrono::duration<double>(Clock::now() - start).count();
|
|
return report;
|
|
}
|
|
|
|
[[nodiscard]] bool IsCurrent() const {
|
|
return Traits::IsPrepared(*m_problem) && Traits::Snapshot(*m_problem) == m_snapshot;
|
|
}
|
|
|
|
[[nodiscard]] const ProblemType &GetProblem() const noexcept {
|
|
return *m_problem;
|
|
}
|
|
|
|
[[nodiscard]] const typename Traits::Manifest &GetManifest() const noexcept {
|
|
return *m_manifest;
|
|
}
|
|
|
|
[[nodiscard]] const mfem::Operator &GetLinearizationOperator() const noexcept {
|
|
return *m_linearization;
|
|
}
|
|
|
|
template <typename CorrectionBlock>
|
|
requires utils::blocks::contains_type_v<
|
|
CorrectionBlock,
|
|
typename FormType::value_blocks>
|
|
[[nodiscard]] mfem::Vector GetCorrectionBlock(mfem::Vector &correction) const {
|
|
if (correction.Size() != Height()) {
|
|
throw std::invalid_argument("A correction block view requires a complete correction vector.");
|
|
}
|
|
constexpr int index = utils::blocks::type_index_v<CorrectionBlock, typename FormType::value_blocks>;
|
|
return mfem::Vector(
|
|
correction.GetData() + m_manifest->layout().offset(utils::blocks::value_block<index>{}),
|
|
m_manifest->layout().size(utils::blocks::value_block<index>{})
|
|
);
|
|
}
|
|
|
|
template <typename ResidualBlock>
|
|
requires utils::blocks::contains_type_v<
|
|
ResidualBlock,
|
|
typename FormType::residual_blocks>
|
|
[[nodiscard]] mfem::Vector GetResidualBlock(const mfem::Vector &residual) const {
|
|
if (residual.Size() != Width()) {
|
|
throw std::invalid_argument("A residual block view requires a complete residual vector.");
|
|
}
|
|
constexpr int index = utils::blocks::type_index_v<ResidualBlock, typename FormType::residual_blocks>;
|
|
return mfem::Vector(
|
|
const_cast<mfem::real_t *>(residual.GetData()) +
|
|
m_manifest->layout().offset(utils::blocks::residual_block<index>{}),
|
|
m_manifest->layout().size(utils::blocks::residual_block<index>{})
|
|
);
|
|
}
|
|
|
|
[[nodiscard]] const PlanType &GetPlan() const noexcept {
|
|
return m_plan;
|
|
}
|
|
|
|
[[nodiscard]] const StellarPreconditionerLifecycleSnapshot &GetLifecycleSnapshot() const noexcept {
|
|
return m_snapshot;
|
|
}
|
|
|
|
[[nodiscard]] const StellarPreconditionerStatistics &GetStatistics() const noexcept {
|
|
return m_statistics;
|
|
}
|
|
|
|
private:
|
|
void VerifyPreparedProblem() const {
|
|
if (!Traits::IsPrepared(*m_problem)) {
|
|
throw std::logic_error(
|
|
"The stellar-equilibrium problem must be prepared before its preconditioner is prepared or "
|
|
"refreshed."
|
|
);
|
|
}
|
|
if (Traits::StateSize(*m_problem) <= 0 ||
|
|
Traits::StateSize(*m_problem) != Traits::EquationSize(*m_problem)) {
|
|
throw std::logic_error("A stellar-equilibrium preconditioner requires a positive square problem.");
|
|
}
|
|
|
|
const auto &layout = Traits::ManifestOf(*m_problem).layout();
|
|
if (layout.value_offsets().Last() != Traits::StateSize(*m_problem) ||
|
|
layout.residual_offsets().Last() != Traits::EquationSize(*m_problem)) {
|
|
throw std::logic_error(
|
|
"The stellar-equilibrium manifest and discrete problem dimensions are inconsistent."
|
|
);
|
|
}
|
|
}
|
|
|
|
void VerifyCurrent() const {
|
|
if (!IsCurrent()) {
|
|
throw std::logic_error(
|
|
"The stellar-equilibrium preconditioner is stale; call Refresh after preparing a new linearization."
|
|
);
|
|
}
|
|
}
|
|
|
|
template <std::size_t... Indices> void SetupComponents(std::index_sequence<Indices...>) {
|
|
((std::get<Indices>(m_preparedComponents).Setup(*m_problem, std::get<Indices>(m_plan.components())),
|
|
++m_statistics.componentSetups),
|
|
...);
|
|
}
|
|
|
|
template <std::size_t Index>
|
|
[[nodiscard]] std::uint64_t RefreshComponent(const StellarPreconditionerPreparationChanges &changes) {
|
|
using ComponentTuple = std::remove_cvref_t<decltype(m_plan.components())>;
|
|
using Component = std::tuple_element_t<Index, ComponentTuple>;
|
|
if (!detail::componentRequiresRefresh<typename Component::PreparationDependencies>(changes)) {
|
|
return 0;
|
|
}
|
|
return std::get<Index>(m_preparedComponents)
|
|
.Refresh(*m_problem, std::get<Index>(m_plan.components()), changes)
|
|
? 1U
|
|
: 0U;
|
|
}
|
|
|
|
template <std::size_t... Indices>
|
|
[[nodiscard]] std::uint64_t RefreshComponents(
|
|
const StellarPreconditionerPreparationChanges &changes,
|
|
std::index_sequence<Indices...>
|
|
) {
|
|
return (std::uint64_t{0} + ... + RefreshComponent<Indices>(changes));
|
|
}
|
|
|
|
ProblemType *m_problem;
|
|
const typename Traits::Manifest *m_manifest;
|
|
const mfem::Operator *m_linearization;
|
|
PlanType m_plan;
|
|
PreparedComponents m_preparedComponents;
|
|
StellarPreconditionerLifecycleSnapshot m_snapshot;
|
|
mutable StellarPreconditionerStatistics m_statistics;
|
|
};
|
|
|
|
template <
|
|
StellarPreconditionerProblem Problem,
|
|
typename Plan>
|
|
requires PreparedPreconditionerPlanFor<
|
|
std::remove_cvref_t<Plan>,
|
|
std::remove_cvref_t<Problem>>
|
|
[[nodiscard]] auto prepare(
|
|
Problem &problem,
|
|
Plan &&plan
|
|
) {
|
|
using ProblemType = std::remove_cvref_t<Problem>;
|
|
using PlanType = std::remove_cvref_t<Plan>;
|
|
return StellarEquilibriumPreconditioner<ProblemType, PlanType>{problem, std::forward<Plan>(plan)};
|
|
}
|
|
} // namespace mean_field::preconditioning
|