module;

#include <algorithm>
#include <array>
#include <chrono>
#include <cmath>
#include <concepts>
#include <cstddef>
#include <cstdint>
#include <exception>
#include <expected>
#include <functional>
#include <limits>
#include <memory>
#include <optional>
#include <span>
#include <stdexcept>
#include <string>
#include <string_view>
#include <type_traits>
#include <utility>
#include <vector>

#include <mfem.hpp>
#include <mpi.h>

export module mean_field:solver.stellar_context;

export import :normalization.stellar_equilibrium;
export import :preconditioning.stellar_recipe;
export import :seed.stellar_equilibrium_projection;
export import :solver.linear_backend;
export import :solver.newton;
export import :solver.stellar_structure;

export namespace mean_field::solver {
    template <typename Model, typename Discretization, typename PreconditionerPrescription, typename LinearBackend>
    class StellarEquilibriumContext;

    template <typename Context, typename NewtonConfiguration, typename Observer> class StellarEquilibriumSolver;
} // namespace mean_field::solver

export namespace mean_field::solver::detail {
    // Internal, synchronous experiment access; not a stable solver API.
    struct StellarEquilibriumContextDiagnostics;
}

namespace mean_field::solver::detail {
    template <typename Model, typename Discretization>
    using StellarContextProblem =
        equilibrium::StellarEquilibriumProblem<std::remove_cvref_t<Model>, std::remove_cvref_t<Discretization>>;

    template <typename Problem>
    using StellarContextNormalizedOperator =
        normalization::NormalizedStellarEquilibriumOperator<std::remove_cvref_t<Problem>>;

    template <typename Prescription, typename Problem>
    using StellarContextPhysicalInverse =
        preconditioning::PreparedStellarInverseType<std::remove_cvref_t<Prescription>, std::remove_cvref_t<Problem>>;

    template <typename Problem, typename PhysicalInverse>
    using StellarContextNormalizedInverse = normalization::
        NormalizedStellarPreconditioner<std::remove_cvref_t<Problem>, std::remove_cvref_t<PhysicalInverse>>;

    template <typename Configuration, typename Operator, typename Preconditioner>
    using StellarContextLinearBackend = PreparedLinearBackendType<
        std::remove_cvref_t<Configuration>,
        std::remove_cvref_t<Operator>,
        std::remove_cvref_t<Preconditioner>>;

    struct StellarEquilibriumContextAssembly;

    [[nodiscard]] constexpr std::string_view
    PreparationStageName(const operators::StellarEquilibriumPreparationStage stage) noexcept {
        using Stage = operators::StellarEquilibriumPreparationStage;
        switch (stage) {
        case Stage::generated_geometry:
            return "generated geometry";
        case Stage::gravity:
            return "gravity";
        case Stage::barotropic_closure:
            return "barotropic closure";
        case Stage::hydrostatic_equilibrium:
            return "hydrostatic equilibrium";
        case Stage::displacement_residual:
            return "displacement residual";
        case Stage::pressure_force:
            return "pressure force";
        case Stage::gravity_displacement_force:
            return "gravity displacement force";
        case Stage::rotational_displacement_force:
            return "rotational displacement force";
        case Stage::displacement_composition:
            return "displacement composition";
        case Stage::mass_normalization:
            return "mass normalization";
        case Stage::model_specification:
            return "model specification";
        case Stage::unspecified:
            return "unspecified preparation";
        }
        return "unknown preparation";
    }

    [[nodiscard]] inline bool AllStellarRanksSatisfy(
        const bool localValue,
        const MPI_Comm communicator
    ) {
        int localValueAsInteger  = localValue ? 1 : 0;
        int globalValueAsInteger = 0;
        if (MPI_Allreduce(&localValueAsInteger, &globalValueAsInteger, 1, MPI_INT, MPI_MIN, communicator) !=
            MPI_SUCCESS) {
            throw std::runtime_error("The stellar-equilibrium solver could not complete a collective validity check.");
        }
        return globalValueAsInteger != 0;
    }

    [[nodiscard]] inline int MaximumStellarRankValue(
        const int localValue,
        const MPI_Comm communicator
    ) {
        int globalValue = 0;
        if (MPI_Allreduce(&localValue, &globalValue, 1, MPI_INT, MPI_MAX, communicator) != MPI_SUCCESS) {
            throw std::runtime_error("The stellar-equilibrium solver could not combine a distributed status.");
        }
        return globalValue;
    }

    inline void RequireCollectivelyIdenticalNewtonOptions(
        const nonlinear::NewtonOptions &options,
        const MPI_Comm communicator
    ) {
        const std::array<double, 9> localRealValues{
            options.relativeTolerance,
            options.absoluteTolerance,
            options.linearSolve.relativeTolerance,
            options.linearSolve.absoluteTolerance,
            options.backtracking.initialStepLength,
            options.backtracking.contractionFactor,
            options.backtracking.fractionToBoundarySafety,
            options.backtracking.sufficientDecrease,
            options.backtracking.minimumStepLength
        };
        std::array<double, 9> minimumRealValues{};
        std::array<double, 9> maximumRealValues{};
        const std::array<int, 3> localIntegerValues{
            options.maximumIterations, options.linearSolve.maximumIterations, options.backtracking.maximumTrials
        };
        std::array<int, 3> minimumIntegerValues{};
        std::array<int, 3> maximumIntegerValues{};

        if (MPI_Allreduce(
                localRealValues.data(), minimumRealValues.data(), static_cast<int>(localRealValues.size()), MPI_DOUBLE,
                MPI_MIN, communicator
            ) != MPI_SUCCESS ||
            MPI_Allreduce(
                localRealValues.data(), maximumRealValues.data(), static_cast<int>(localRealValues.size()), MPI_DOUBLE,
                MPI_MAX, communicator
            ) != MPI_SUCCESS ||
            MPI_Allreduce(
                localIntegerValues.data(), minimumIntegerValues.data(), static_cast<int>(localIntegerValues.size()),
                MPI_INT, MPI_MIN, communicator
            ) != MPI_SUCCESS ||
            MPI_Allreduce(
                localIntegerValues.data(), maximumIntegerValues.data(), static_cast<int>(localIntegerValues.size()),
                MPI_INT, MPI_MAX, communicator
            ) != MPI_SUCCESS) {
            throw std::runtime_error("The stellar-equilibrium solver could not validate its distributed options.");
        }

        if (minimumRealValues != maximumRealValues || minimumIntegerValues != maximumIntegerValues) {
            throw std::invalid_argument(
                "Newton requires identical nonlinear, line-search, and linear-solve options on every communicator "
                "rank."
            );
        }
    }
} // namespace mean_field::solver::detail

export namespace mean_field::solver {
    template <typename Candidate, typename Model>
    concept StellarEquilibriumInitialStateFor =
        equilibrium::StellarEquilibriumModel<std::remove_cvref_t<Model>> &&
        std::move_constructible<std::remove_cvref_t<Candidate>> &&
        seed::RadialProfileProjectableModel<std::remove_cvref_t<Model>> &&
        seed::RadialSeedStrategyFor<std::remove_cvref_t<Candidate>, std::remove_cvref_t<Model>>;

    template <typename Model>
    concept DefaultStellarEquilibriumInitialStateAvailableFor =
        equilibrium::StellarEquilibriumModel<std::remove_cvref_t<Model>> &&
        std::remove_cvref_t<Model>::template containsSpecification<models::FixedCentralDensity> &&
        seed::RadialProfileProjectableModel<std::remove_cvref_t<Model>> &&
        seed::RadialSeedStrategyFor<seed::LaneEmden, std::remove_cvref_t<Model>>;

    template <typename Model, typename Discretization, typename PreconditionerPrescription, typename LinearBackend>
    concept StellarEquilibriumContextConfiguration =
        equilibrium::StellarEquilibriumModel<std::remove_cvref_t<Model>> &&
        equilibrium::StellarDiscretizationType<std::remove_cvref_t<Discretization>> &&
        std::move_constructible<std::remove_cvref_t<Model>> &&
        std::move_constructible<std::remove_cvref_t<Discretization>> &&
        preconditioning::StellarPreconditionerPrescription<std::remove_cvref_t<PreconditionerPrescription>> &&
        LinearBackendConfiguration<std::remove_cvref_t<LinearBackend>> &&
        equilibrium::StellarEquilibriumModelDiscretizationCompatible<
            std::remove_cvref_t<Model>,
            std::remove_cvref_t<Discretization>> &&
        normalization::NormalizableStellarEquilibriumProblem<detail::StellarContextProblem<Model, Discretization>> &&
        preconditioning::StellarPreconditionerRuntimeAvailableFor<
            std::remove_cvref_t<PreconditionerPrescription>,
            detail::StellarContextProblem<Model, Discretization>> &&
        LinearBackendRuntimeAvailableFor<
            std::remove_cvref_t<LinearBackend>,
            detail::StellarContextNormalizedOperator<detail::StellarContextProblem<Model, Discretization>>,
            detail::StellarContextNormalizedInverse<
                detail::StellarContextProblem<Model, Discretization>,
                detail::StellarContextPhysicalInverse<
                    PreconditionerPrescription,
                    detail::StellarContextProblem<Model, Discretization>>>>;

    /*
     * The user-owned lifetime root for one discretized stellar equilibrium
     * problem. It owns the physical problem, accepted/trial state, prepared
     * preconditioner, linear backend, and all Newton work vectors. Solvers
     * borrow it exclusively and contain no numerical storage.
     */
    template <typename Model, typename Discretization, typename PreconditionerPrescription, typename LinearBackend>
    class StellarEquilibriumContext final {
    public:
        static_assert(StellarEquilibriumContextConfiguration<
                      Model,
                      Discretization,
                      PreconditionerPrescription,
                      LinearBackend>);

        using ModelType                      = std::remove_cvref_t<Model>;
        using DiscretizationType             = std::remove_cvref_t<Discretization>;
        using PreconditionerPrescriptionType = std::remove_cvref_t<PreconditionerPrescription>;
        using LinearBackendConfigurationType = std::remove_cvref_t<LinearBackend>;
        using ProblemType                    = detail::StellarContextProblem<ModelType, DiscretizationType>;
        using EvaluationReport               = StellarEquilibriumEvaluationReport<ProblemType>;

        StellarEquilibriumContext(const StellarEquilibriumContext &)            = delete;
        StellarEquilibriumContext &operator=(const StellarEquilibriumContext &) = delete;
        StellarEquilibriumContext(StellarEquilibriumContext &&)                 = delete;
        StellarEquilibriumContext &operator=(StellarEquilibriumContext &&)      = delete;
        ~StellarEquilibriumContext()                                            = default;

        [[nodiscard]] bool isReady() const {
            return m_state != nullptr && m_state->IsReady();
        }

        [[nodiscard]] bool hasActiveSolver() const noexcept {
            return m_hasActiveSolver;
        }

        [[nodiscard]] MPI_Comm communicator() const & {
            if (m_state == nullptr) {
                throw std::logic_error("A stellar-equilibrium context has no runtime state.");
            }
            return m_state->Problem().GetCommunicator();
        }

        [[nodiscard]] MPI_Comm communicator() const && = delete;

    private:
        template <typename, typename, typename> friend class StellarEquilibriumSolver;
        friend struct detail::StellarEquilibriumContextAssembly;
        friend struct detail::StellarEquilibriumContextDiagnostics;

        using NormalizedOperatorType = detail::StellarContextNormalizedOperator<ProblemType>;
        using PhysicalInverseType = detail::StellarContextPhysicalInverse<PreconditionerPrescriptionType, ProblemType>;
        using NormalizedInverseType = detail::StellarContextNormalizedInverse<ProblemType, PhysicalInverseType>;
        using LinearBackendType     = detail::
            StellarContextLinearBackend<LinearBackendConfigurationType, NormalizedOperatorType, NormalizedInverseType>;
        using Storage       = detail::StellarStructureStorage<ProblemType>;
        using InitialBundle = std::shared_ptr<Storage>;

        struct DependencyLedger final {
            std::array<std::byte, 9> identities{};
            operators::StellarEquilibriumDependencies dependencies{};

            explicit DependencyLedger(const ProblemType &problem) {
                if (problem.IsPrepared()) {
                    dependencies = problem.GetLinearizationDependencies();
                    AdvanceState();
                } else {
                    dependencies = {
                        .discretization     = Stamp(0),
                        .density            = Stamp(1),
                        .surfaceDeformation = Stamp(2),
                        .gravityGradient    = Stamp(3),
                        .gravityPotential   = Stamp(4),
                        .enthalpy           = Stamp(5),
                        .bernoulliConstant  = Stamp(6),
                        .rotation           = Stamp(7),
                        .targetMass         = Stamp(8)
                    };
                }
            }

            void AdvanceState() noexcept {
                ++dependencies.density.revision;
                ++dependencies.surfaceDeformation.revision;
                ++dependencies.gravityGradient.revision;
                ++dependencies.gravityPotential.revision;
                ++dependencies.enthalpy.revision;
                ++dependencies.bernoulliConstant.revision;
                ++dependencies.rotation.revision;
                ++dependencies.targetMass.revision;
            }

        private:
            [[nodiscard]] operators::StellarEquilibriumDependencyStamp Stamp(const std::size_t index) const noexcept {
                return {
                    .identity =
                        static_cast<std::uint64_t>(reinterpret_cast<std::uintptr_t>(std::addressof(identities[index]))),
                    .revision = 1
                };
            }
        };

        struct State final {
            explicit State(
                InitialBundle initial,
                PreconditionerPrescriptionType preconditionerPrescription,
                LinearBackendConfigurationType linearBackendConfiguration
            )
                : storage(std::move(initial)),
                  dependencyLedger(RequireProblem(storage)),
                  acceptedNormalizedState(RequireProblem(storage).StateSize()),
                  trialNormalizedState(RequireProblem(storage).StateSize()),
                  acceptedNormalizedResidual(RequireProblem(storage).EquationSize()),
                  trialNormalizedResidual(RequireProblem(storage).EquationSize()),
                  linearRightHandSide(RequireProblem(storage).EquationSize()),
                  normalizedCorrection(RequireProblem(storage).StateSize()),
                  physicalCorrection(RequireProblem(storage).StateSize()),
                  volumeDisplacementDirection(
                      RequireProblem(storage).GetPhysicalOperator().GetDomainDeformation().volumeDisplacementSize()
                  ),
                  candidatePhysicalState(RequireProblem(storage).StateSize()),
                  normalizedOperator(std::make_unique<NormalizedOperatorType>(RequireProblem(storage))) {
                ValidateInitialState();
                InitializeWorkspaces();
                PrepareInitialOperator();
                InitializeGeometryPreflightRules();

                physicalInverse   = std::unique_ptr<PhysicalInverseType>{new PhysicalInverseType(
                    PreparePhysicalInverse(std::move(preconditionerPrescription), RequireProblem(storage))
                )};
                normalizedInverse = std::unique_ptr<NormalizedInverseType>{
                    new NormalizedInverseType(normalizedOperator->MakeScaledPreconditioner(*physicalInverse))
                };
                linearBackend = std::unique_ptr<LinearBackendType>{new LinearBackendType(PrepareLinearBackend(
                    std::move(linearBackendConfiguration), *normalizedOperator, *normalizedInverse,
                    RequireProblem(storage).GetCommunicator()
                ))};

                if (!IsReady()) {
                    throw std::logic_error("The stellar-equilibrium context did not produce a complete runtime.");
                }
            }

            State(const State &)            = delete;
            State &operator=(const State &) = delete;
            State(State &&)                 = delete;
            State &operator=(State &&)      = delete;

            [[nodiscard]] bool IsReady() const {
                const ProblemType *problem = ProblemPointer(storage);
                return problem != nullptr && storage->physicalState != nullptr && normalizedOperator != nullptr &&
                       physicalInverse != nullptr && normalizedInverse != nullptr && linearBackend != nullptr &&
                       normalizedOperator->IsPrepared() && physicalInverse->IsCurrent() &&
                       normalizedInverse->IsCurrent() && linearBackend->IsReady() &&
                       std::addressof(physicalInverse->GetProblem()) == problem &&
                       std::addressof(linearBackend->GetOperator()) == normalizedOperator.get() &&
                       std::addressof(linearBackend->GetPreconditioner()) == normalizedInverse.get() &&
                       CommunicatorsCompatible(linearBackend->GetCommunicator(), problem->GetCommunicator()) &&
                       normalizedOperator->Height() == problem->EquationSize() &&
                       normalizedOperator->Width() == problem->StateSize() &&
                       std::isfinite(acceptedMinimumJacobianDeterminant) && acceptedMinimumJacobianDeterminant > 0.0 &&
                       physicalInverse->Height() == problem->StateSize() &&
                       physicalInverse->Width() == problem->EquationSize() &&
                       normalizedInverse->Height() == problem->StateSize() &&
                       normalizedInverse->Width() == problem->EquationSize() &&
                       linearBackend->RightHandSideSize() == problem->EquationSize() &&
                       linearBackend->CorrectionSize() == problem->StateSize() &&
                       acceptedNormalizedState.Size() == problem->StateSize() &&
                       trialNormalizedState.Size() == problem->StateSize() &&
                       acceptedNormalizedResidual.Size() == problem->EquationSize() &&
                       trialNormalizedResidual.Size() == problem->EquationSize() &&
                       linearRightHandSide.Size() == problem->EquationSize() &&
                       normalizedCorrection.Size() == problem->StateSize() &&
                       physicalCorrection.Size() == problem->StateSize() &&
                       volumeDisplacementDirection.Size() ==
                           problem->GetPhysicalOperator().GetDomainDeformation().volumeDisplacementSize() &&
                       candidatePhysicalState.Size() == problem->StateSize() &&
                       storage->physicalState->Size() == problem->StateSize();
            }

            [[nodiscard]] ProblemType &Problem() noexcept {
                return *storage->problem;
            }

            [[nodiscard]] const ProblemType &Problem() const noexcept {
                return *storage->problem;
            }

            [[nodiscard]] mfem::Vector &AcceptedPhysicalState() noexcept {
                return *storage->physicalState;
            }

            [[nodiscard]] const mfem::Vector &AcceptedPhysicalState() const noexcept {
                return *storage->physicalState;
            }

            void BeginEvaluation() {
                AdvanceViewGeneration();
                storage->certification = detail::StellarViewCertification::unavailable;
            }

            [[nodiscard]] std::expected<
                deformation::DomainDeformationGeometryReport,
                operators::StellarEquilibriumPreparationRejection>
            PrepareTrial() {
                dependencyLedger.AdvanceState();
                auto result = TryPrepareOperator(trialNormalizedState);
                if (!result.has_value()) {
                    return std::unexpected(result.error());
                }
                normalizedOperator->BuildResidual(trialNormalizedResidual);
                return result->physical.generatedGeometry;
            }

            void RestoreAccepted() {
                dependencyLedger.AdvanceState();
                const auto preparation             = PrepareOperator(acceptedNormalizedState);
                acceptedMinimumJacobianDeterminant = preparation.physical.generatedGeometry.minimumJacobianDeterminant;
                normalizedOperator->BuildResidual(acceptedNormalizedResidual);
                AcceptedPhysicalState() = normalizedOperator->GetPhysicalState();
                trialNormalizedState    = acceptedNormalizedState;
                trialNormalizedResidual = acceptedNormalizedResidual;
                candidatePhysicalState  = AcceptedPhysicalState();
                physicalInverse->Refresh();
                if (!IsReady()) {
                    throw std::logic_error("The stellar-equilibrium context could not restore its accepted state.");
                }
            }

            void PrepareTrialCommit() {
                physicalInverse->Refresh();
                if (!IsReady()) {
                    throw std::logic_error(
                        "The stellar-equilibrium context could not prepare to commit a Newton step."
                    );
                }
            }

            void CommitPreparedTrial(const double minimumJacobianDeterminant) noexcept {
                CopyValues(trialNormalizedState, acceptedNormalizedState);
                CopyValues(trialNormalizedResidual, acceptedNormalizedResidual);
                CopyValues(normalizedOperator->GetPhysicalState(), AcceptedPhysicalState());
                CopyValues(AcceptedPhysicalState(), candidatePhysicalState);
                acceptedMinimumJacobianDeterminant = minimumJacobianDeterminant;
                storage->certification             = detail::StellarViewCertification::unavailable;
            }

            void AdvanceCorrectionWarmStart(const double acceptedStepLength) {
                /*
                 * Under the frozen-Jacobian Newton model,
                 *
                 *   F(x + alpha p) ~= (1 - alpha) F(x),
                 *
                 * so (1 - alpha) p is the corresponding first guess for the
                 * correction at the accepted state.  This reuses only a
                 * vector, not a Krylov basis.  If an extreme but valid user
                 * step overflows the heuristic, fall back collectively to a
                 * cold start without invalidating the accepted state.
                 */
                normalizedCorrection *= 1.0 - acceptedStepLength;
                if (!detail::AllStellarRanksSatisfy(AllFinite(normalizedCorrection), Problem().GetCommunicator())) {
                    normalizedCorrection = 0.0;
                }
            }

            [[nodiscard]] deformation::LargestSafeNewtonStepSizeEstimate EstimateLargestSafeStepSize(
                const double maximumStepSize,
                const double fractionToBoundarySafety
            ) {
                normalizedOperator->DenormalizeState(normalizedCorrection, physicalCorrection);
                const auto &physicalOperator = Problem().GetPhysicalOperator();
                Problem().BuildVolumeDisplacementDirection(physicalCorrection, volumeDisplacementDirection);

                const fem::FEM &finiteElements =
                    equilibrium::detail::StellarEquilibriumProblemFactory::FiniteElementModel(Problem());
                return deformation::estimate_largest_safe_newton_step_size(
                    Problem().GetDiscretization().domainMapper(), *finiteElements.displacementFes,
                    *finiteElements.compactificationCoordinate, physicalOperator.GetGeneratedVolumeDisplacement(),
                    volumeDisplacementDirection, geometryPreflightRules,
                    {.maximumStepSize          = maximumStepSize,
                     .determinantFloor         = 0.0,
                     .fractionToBoundarySafety = fractionToBoundarySafety}
                );
            }

            std::shared_ptr<Storage> storage;
            DependencyLedger dependencyLedger;
            mfem::Vector acceptedNormalizedState;
            mfem::Vector trialNormalizedState;
            mfem::Vector acceptedNormalizedResidual;
            mfem::Vector trialNormalizedResidual;
            mfem::Vector linearRightHandSide;
            mfem::Vector normalizedCorrection;
            mfem::Vector physicalCorrection;
            mfem::Vector volumeDisplacementDirection;
            mfem::Vector candidatePhysicalState;
            std::vector<deformation::NewtonStepGeometryRule> geometryPreflightRules;
            double acceptedMinimumJacobianDeterminant{std::numeric_limits<double>::quiet_NaN()};
            std::unique_ptr<NormalizedOperatorType> normalizedOperator;
            std::unique_ptr<PhysicalInverseType> physicalInverse;
            std::unique_ptr<NormalizedInverseType> normalizedInverse;
            std::unique_ptr<LinearBackendType> linearBackend;

        private:
            [[nodiscard]] static ProblemType *ProblemPointer(const std::shared_ptr<Storage> &candidate) noexcept {
                return candidate == nullptr ? nullptr : candidate->problem.get();
            }

            [[nodiscard]] static ProblemType &RequireProblem(const std::shared_ptr<Storage> &candidate) {
                ProblemType *problem = ProblemPointer(candidate);
                if (problem == nullptr) {
                    throw std::invalid_argument("A stellar-equilibrium context requires an owned physical problem.");
                }
                return *problem;
            }

            void AdvanceViewGeneration() {
                if (storage->viewGeneration == std::numeric_limits<std::uint64_t>::max()) {
                    throw std::overflow_error("The stellar-equilibrium view generation was exhausted.");
                }
                ++storage->viewGeneration;
            }

            void ValidateInitialState() const {
                const auto &problem         = RequireProblem(storage);
                const MPI_Comm communicator = problem.GetCommunicator();
                if (communicator == MPI_COMM_NULL) {
                    throw std::invalid_argument("A stellar-equilibrium context requires a non-null communicator.");
                }
                if (!detail::AllStellarRanksSatisfy(
                        storage->physicalState != nullptr && storage->physicalState->Size() == problem.StateSize(),
                        communicator
                    )) {
                    throw std::invalid_argument(
                        "A stellar-equilibrium initial state does not match its local physical problem on every rank."
                    );
                }
                if (!detail::AllStellarRanksSatisfy(AllFinite(*storage->physicalState), communicator)) {
                    throw std::invalid_argument("The stellar-equilibrium initial state must be finite on every rank.");
                }
                if constexpr (ProblemType::generatedRotationProviderCount == 0) {
                    if (!detail::AllStellarRanksSatisfy(storage->rotation != nullptr, communicator)) {
                        throw std::invalid_argument(
                            "A model without generated rotation requires a prescribed rigid rotation on every rank."
                        );
                    }
                } else {
                    if (!detail::AllStellarRanksSatisfy(storage->rotation == nullptr, communicator)) {
                        throw std::invalid_argument(
                            "A model with generated rotation cannot accept a prescribed rigid rotation on any rank."
                        );
                    }
                }
            }

            [[nodiscard]] static bool AllFinite(const mfem::Vector &values) noexcept {
                for (int index = 0; index < values.Size(); ++index) {
                    if (!std::isfinite(values(index))) {
                        return false;
                    }
                }
                return true;
            }

            static void CopyValues(
                const mfem::Vector &source,
                mfem::Vector &destination
            ) noexcept {
                for (int index = 0; index < source.Size(); ++index) {
                    destination(index) = source(index);
                }
            }

            [[nodiscard]] static bool CommunicatorsCompatible(
                const MPI_Comm left,
                const MPI_Comm right
            ) {
                if (left == MPI_COMM_NULL || right == MPI_COMM_NULL) {
                    return false;
                }
                int comparison = MPI_UNEQUAL;
                if (MPI_Comm_compare(left, right, &comparison) != MPI_SUCCESS) {
                    return false;
                }
                return comparison == MPI_IDENT || comparison == MPI_CONGRUENT;
            }

            void InitializeWorkspaces() {
                normalizedOperator->NormalizeState(AcceptedPhysicalState(), acceptedNormalizedState);
                trialNormalizedState        = acceptedNormalizedState;
                acceptedNormalizedResidual  = 0.0;
                trialNormalizedResidual     = 0.0;
                linearRightHandSide         = 0.0;
                normalizedCorrection        = 0.0;
                physicalCorrection          = 0.0;
                volumeDisplacementDirection = 0.0;
                candidatePhysicalState      = AcceptedPhysicalState();
            }

            void PrepareInitialOperator() {
                const auto preparation             = PrepareOperator(acceptedNormalizedState);
                acceptedMinimumJacobianDeterminant = preparation.physical.generatedGeometry.minimumJacobianDeterminant;
                AcceptedPhysicalState()            = normalizedOperator->GetPhysicalState();
                candidatePhysicalState             = AcceptedPhysicalState();
                normalizedOperator->BuildResidual(acceptedNormalizedResidual);
                trialNormalizedResidual = acceptedNormalizedResidual;
            }

            void AppendGeometryPreflightRule(
                const int element,
                const mfem::IntegrationRule &integrationRule
            ) {
                geometryPreflightRules.push_back({.element = element, .integrationRule = &integrationRule});
            }

            void InitializeGeometryPreflightRules() {
                const ProblemType &problem = Problem();
                const fem::FEM &finiteElements =
                    equilibrium::detail::StellarEquilibriumProblemFactory::FiniteElementModel(problem);
                if (finiteElements.mesh == nullptr || finiteElements.displacementFes == nullptr ||
                    finiteElements.compactificationCoordinate == nullptr) {
                    throw std::logic_error(
                        "The Newton geometry preflight requires complete displacement geometry data."
                    );
                }
                if (!problem.GetPhysicalOperator().GetDomainDeformation().descriptor().linearOnReferenceGeometry) {
                    throw std::invalid_argument(
                        "The Newton geometry preflight requires a domain deformation that is linear on the "
                        "reference geometry."
                    );
                }

                geometryPreflightRules.clear();
                geometryPreflightRules.reserve(
                    static_cast<std::size_t>(finiteElements.mesh->GetNE()) * static_cast<std::size_t>(10)
                );

                const int dimension = problem.GetDiscretization().domainMapper().GetDimension();
                for (int element = 0; element < finiteElements.mesh->GetNE(); ++element) {
                    const mfem::FiniteElement *finiteElement = finiteElements.displacementFes->GetFE(element);
                    mfem::ElementTransformation *transformation =
                        finiteElements.mesh->GetElementTransformation(element);
                    if (finiteElement == nullptr || transformation == nullptr) {
                        throw std::logic_error(
                            "The Newton geometry preflight encountered incomplete element geometry data."
                        );
                    }
                    const int geometryInspectionOrder =
                        std::max(finiteElement->GetOrder() + 2, 2 * dimension * finiteElement->GetOrder());
                    AppendGeometryPreflightRule(
                        element, mfem::IntRules.Get(transformation->GetGeometryType(), geometryInspectionOrder)
                    );
                }

                const auto appendPreparedRules = [this](const auto &preparedOperator) {
                    preparedOperator.VisitMappedGeometryRules(
                        [this](const int element, const mfem::IntegrationRule &integrationRule) {
                            AppendGeometryPreflightRule(element, integrationRule);
                        }
                    );
                };

                const auto &physicalOperator = problem.GetPhysicalOperator();
                const auto &gravityGeometry  = physicalOperator.GetGravityContext().GetGeometryContext();
                appendPreparedRules(gravityGeometry.GetMassOperator());
                appendPreparedRules(gravityGeometry.GetSourceOperator());
                appendPreparedRules(physicalOperator.GetBarotropicClosureOperator());
                appendPreparedRules(physicalOperator.GetHydrostaticOperator());

                const auto &displacementOperator = physicalOperator.GetDisplacementOperator();
                appendPreparedRules(displacementOperator.GetPressureOperator());
                appendPreparedRules(displacementOperator.GetGravityOperator());
                appendPreparedRules(displacementOperator.GetRotationalOperator());
                appendPreparedRules(physicalOperator.GetMassNormalizationOperator());

                if constexpr (ProblemType::hasFixedAngularMomentum) {
                    appendPreparedRules(problem.GetPreparedOperator().GetAngularMomentumConstraint());
                }

                const auto ruleLess = [](const deformation::NewtonStepGeometryRule &left,
                                         const deformation::NewtonStepGeometryRule &right) {
                    if (left.element != right.element) {
                        return left.element < right.element;
                    }
                    return std::less<const mfem::IntegrationRule *>{}(left.integrationRule, right.integrationRule);
                };
                std::sort(geometryPreflightRules.begin(), geometryPreflightRules.end(), ruleLess);
                geometryPreflightRules.erase(
                    std::unique(
                        geometryPreflightRules.begin(), geometryPreflightRules.end(),
                        [](const deformation::NewtonStepGeometryRule &left,
                           const deformation::NewtonStepGeometryRule &right) {
                            return left.element == right.element && left.integrationRule == right.integrationRule;
                        }
                    ),
                    geometryPreflightRules.end()
                );
            }

            [[nodiscard]] auto PrepareOperator(const mfem::Vector &normalizedState) {
                if constexpr (ProblemType::generatedRotationProviderCount == 0) {
                    return normalizedOperator->Prepare(
                        normalizedState, dependencyLedger.dependencies, *storage->rotation
                    );
                } else {
                    return normalizedOperator->Prepare(normalizedState, dependencyLedger.dependencies);
                }
            }

            [[nodiscard]] auto TryPrepareOperator(const mfem::Vector &normalizedState) {
                if constexpr (ProblemType::generatedRotationProviderCount == 0) {
                    return normalizedOperator->TryPrepare(
                        normalizedState, dependencyLedger.dependencies, *storage->rotation
                    );
                } else {
                    return normalizedOperator->TryPrepare(normalizedState, dependencyLedger.dependencies);
                }
            }

            [[nodiscard]] static PhysicalInverseType PreparePhysicalInverse(
                PreconditionerPrescriptionType prescription,
                const ProblemType &physicalProblem
            ) {
                using preconditioning::prepareStellarPreconditioner;
                return prepareStellarPreconditioner(std::move(prescription), physicalProblem);
            }

            [[nodiscard]] static LinearBackendType PrepareLinearBackend(
                LinearBackendConfigurationType configuration,
                const NormalizedOperatorType &operation,
                NormalizedInverseType &preconditioner,
                MPI_Comm communicator
            ) {
                return prepareLinearBackend(std::move(configuration), operation, preconditioner, communicator);
            }
        };

        template <StellarEquilibriumInitialStateFor<ModelType> InitialState>
        explicit StellarEquilibriumContext(
            ModelType model,
            DiscretizationType discretization,
            PreconditionerPrescriptionType preconditionerPrescription,
            LinearBackendConfigurationType linearBackend,
            InitialState initialState,
            seed::StellarEquilibriumProjectionOptions projectionOptions,
            std::optional<physics::RigidRotation> prescribedRotation
        )
            : m_state(
                  std::make_unique<State>(
                      BuildInitialBundle(
                          std::move(model),
                          std::move(discretization),
                          std::move(initialState),
                          projectionOptions,
                          std::move(prescribedRotation)
                      ),
                      std::move(preconditionerPrescription),
                      std::move(linearBackend)
                  )
              ) {
        }

        template <StellarEquilibriumInitialStateFor<ModelType> InitialState>
        [[nodiscard]] static InitialBundle BuildInitialBundle(
            ModelType model,
            DiscretizationType discretization,
            InitialState initialState,
            const seed::StellarEquilibriumProjectionOptions &projectionOptions,
            std::optional<physics::RigidRotation> prescribedRotation
        ) {
            auto problem = equilibrium::detail::StellarEquilibriumProblemFactory::CreateOwned(
                std::move(model), std::move(discretization)
            );
            auto physicalState = std::make_unique<mfem::Vector>(
                seed::makeProjectedEquilibriumState(*problem, initialState, projectionOptions).values
            );
            std::unique_ptr<physics::RigidRotation> rotation;
            if (prescribedRotation.has_value()) {
                rotation = std::make_unique<physics::RigidRotation>(std::move(*prescribedRotation));
            }
            return std::make_shared<Storage>(std::move(problem), std::move(physicalState), std::move(rotation));
        }

        void AcquireSolver() {
            if (m_hasActiveSolver) {
                throw std::logic_error("A stellar-equilibrium context already has an active solver.");
            }
            if (!isReady()) {
                throw std::logic_error("A solver requires a complete, current stellar-equilibrium context.");
            }
            m_hasActiveSolver = true;
        }

        void ReleaseSolver() noexcept {
            m_hasActiveSolver = false;
        }

        std::unique_ptr<State> m_state;
        bool m_hasActiveSolver{false};
    };

    template <typename Candidate> struct IsStellarEquilibriumContext : std::false_type { };

    template <typename Model, typename Discretization, typename PreconditionerPrescription, typename LinearBackend>
    struct IsStellarEquilibriumContext<
        StellarEquilibriumContext<Model, Discretization, PreconditionerPrescription, LinearBackend>> : std::true_type {
    };

    template <typename Candidate>
    concept StellarEquilibriumContextType = IsStellarEquilibriumContext<std::remove_cvref_t<Candidate>>::value;
} // namespace mean_field::solver

export namespace mean_field::solver::detail {
    struct StellarEquilibriumContextDiagnostics final {
        // The callback must not retain references to runtime storage. It may
        // prepare trial states, but accepted vectors must remain unchanged.
        // Restore the production preparation and correction on every exit.
        template <typename Context, typename Callback>
        static void WithState(Context &context, Callback &&callback) {
            if (context.hasActiveSolver() || !context.isReady()) {
                throw std::logic_error("Diagnostics require a ready context with no active solver.");
            }
            auto &state = *context.m_state;
            mfem::Vector savedCorrection(state.normalizedCorrection);
            context.AcquireSolver();
            try {
                state.BeginEvaluation();
                std::invoke(std::forward<Callback>(callback), state,
                    equilibrium::detail::StellarEquilibriumProblemFactory::FiniteElementModel(state.Problem()));
                state.RestoreAccepted();
                state.normalizedCorrection = savedCorrection;
                context.ReleaseSolver();
            } catch (...) {
                const auto original = std::current_exception();
                try {
                    state.RestoreAccepted();
                    state.normalizedCorrection = savedCorrection;
                } catch (...) {
                    context.ReleaseSolver();
                    throw;
                }
                context.ReleaseSolver();
                std::rethrow_exception(original);
            }
        }
    };
}

namespace mean_field::solver::detail {
    [[nodiscard]] inline physics::RigidRotation ZeroRigidRotation() {
        mfem::Vector angularVelocity(3);
        mfem::Vector center(3);
        angularVelocity = 0.0;
        center          = 0.0;
        return physics::RigidRotation{angularVelocity, center};
    }

    struct StellarEquilibriumContextAssembly final {
        template <
            typename Model,
            typename Discretization,
            typename PreconditionerPrescription,
            typename LinearBackend,
            typename InitialState>
        [[nodiscard]] static auto FromSeed(
            Model model,
            Discretization discretization,
            PreconditionerPrescription preconditionerPrescription,
            LinearBackend linearBackend,
            InitialState initialState,
            seed::StellarEquilibriumProjectionOptions projectionOptions,
            std::optional<physics::RigidRotation> prescribedRotation
        ) {
            using Context = StellarEquilibriumContext<
                std::remove_cvref_t<Model>, std::remove_cvref_t<Discretization>,
                std::remove_cvref_t<PreconditionerPrescription>, std::remove_cvref_t<LinearBackend>>;
            return Context{std::move(model),
                           std::move(discretization),
                           std::move(preconditionerPrescription),
                           std::move(linearBackend),
                           std::move(initialState),
                           std::move(projectionOptions),
                           std::move(prescribedRotation)};
        }
    };
} // namespace mean_field::solver::detail

export namespace mean_field::solver {
    template <
        typename Model,
        equilibrium::StellarDiscretizationType Discretization,
        preconditioning::StellarPreconditionerPrescription PreconditionerPrescription,
        LinearBackendConfiguration LinearBackend,
        typename InitialState>
        requires StellarEquilibriumInitialStateFor<
                     InitialState,
                     Model> &&
                 StellarEquilibriumContextConfiguration<
                     Model,
                     Discretization,
                     PreconditionerPrescription,
                     LinearBackend>
    [[nodiscard]] auto makeContext(
        Model model,
        Discretization discretization,
        PreconditionerPrescription preconditionerPrescription,
        LinearBackend linearBackend,
        InitialState initialState
    ) {
        std::optional<physics::RigidRotation> prescribedRotation;
        if constexpr (detail::StellarContextProblem<Model, Discretization>::generatedRotationProviderCount == 0) {
            prescribedRotation.emplace(detail::ZeroRigidRotation());
        }
        return detail::StellarEquilibriumContextAssembly::FromSeed(
            std::move(model), std::move(discretization), std::move(preconditionerPrescription),
            std::move(linearBackend), std::move(initialState), seed::StellarEquilibriumProjectionOptions{},
            std::move(prescribedRotation)
        );
    }

    template <
        typename Model,
        equilibrium::StellarDiscretizationType Discretization,
        preconditioning::StellarPreconditionerPrescription PreconditionerPrescription,
        LinearBackendConfiguration LinearBackend,
        typename InitialState>
        requires StellarEquilibriumInitialStateFor<
                     InitialState,
                     Model> &&
                 StellarEquilibriumContextConfiguration<
                     Model,
                     Discretization,
                     PreconditionerPrescription,
                     LinearBackend>
    [[nodiscard]] auto makeContext(
        Model model,
        Discretization discretization,
        PreconditionerPrescription preconditionerPrescription,
        LinearBackend linearBackend,
        InitialState initialState,
        seed::StellarEquilibriumProjectionOptions projectionOptions
    ) {
        std::optional<physics::RigidRotation> prescribedRotation;
        if constexpr (detail::StellarContextProblem<Model, Discretization>::generatedRotationProviderCount == 0) {
            prescribedRotation.emplace(detail::ZeroRigidRotation());
        }
        return detail::StellarEquilibriumContextAssembly::FromSeed(
            std::move(model), std::move(discretization), std::move(preconditionerPrescription),
            std::move(linearBackend), std::move(initialState), std::move(projectionOptions),
            std::move(prescribedRotation)
        );
    }

    template <
        typename Model,
        equilibrium::StellarDiscretizationType Discretization,
        preconditioning::StellarPreconditionerPrescription PreconditionerPrescription,
        LinearBackendConfiguration LinearBackend,
        typename InitialState>
        requires StellarEquilibriumInitialStateFor<
                     InitialState,
                     Model> &&
                 StellarEquilibriumContextConfiguration<
                     Model,
                     Discretization,
                     PreconditionerPrescription,
                     LinearBackend> &&
                 (detail::StellarContextProblem<
                      Model,
                      Discretization>::generatedRotationProviderCount == 0)
    [[nodiscard]] auto makeContext(
        Model model,
        Discretization discretization,
        PreconditionerPrescription preconditionerPrescription,
        LinearBackend linearBackend,
        InitialState initialState,
        physics::RigidRotation prescribedRotation
    ) {
        return detail::StellarEquilibriumContextAssembly::FromSeed(
            std::move(model), std::move(discretization), std::move(preconditionerPrescription),
            std::move(linearBackend), std::move(initialState), seed::StellarEquilibriumProjectionOptions{},
            std::optional<physics::RigidRotation>{std::move(prescribedRotation)}
        );
    }

    template <
        typename Model,
        equilibrium::StellarDiscretizationType Discretization,
        preconditioning::StellarPreconditionerPrescription PreconditionerPrescription,
        LinearBackendConfiguration LinearBackend,
        typename InitialState>
        requires StellarEquilibriumInitialStateFor<
                     InitialState,
                     Model> &&
                 StellarEquilibriumContextConfiguration<
                     Model,
                     Discretization,
                     PreconditionerPrescription,
                     LinearBackend> &&
                 (detail::StellarContextProblem<
                      Model,
                      Discretization>::generatedRotationProviderCount == 0)
    [[nodiscard]] auto makeContext(
        Model model,
        Discretization discretization,
        PreconditionerPrescription preconditionerPrescription,
        LinearBackend linearBackend,
        InitialState initialState,
        seed::StellarEquilibriumProjectionOptions projectionOptions,
        physics::RigidRotation prescribedRotation
    ) {
        return detail::StellarEquilibriumContextAssembly::FromSeed(
            std::move(model), std::move(discretization), std::move(preconditionerPrescription),
            std::move(linearBackend), std::move(initialState), std::move(projectionOptions),
            std::optional<physics::RigidRotation>{std::move(prescribedRotation)}
        );
    }

    template <
        typename Model,
        equilibrium::StellarDiscretizationType Discretization,
        preconditioning::StellarPreconditionerPrescription PreconditionerPrescription,
        LinearBackendConfiguration LinearBackend>
        requires DefaultStellarEquilibriumInitialStateAvailableFor<Model> && StellarEquilibriumContextConfiguration<
                                                                                 Model,
                                                                                 Discretization,
                                                                                 PreconditionerPrescription,
                                                                                 LinearBackend>
    [[nodiscard]] auto makeContext(
        Model model,
        Discretization discretization,
        PreconditionerPrescription preconditionerPrescription,
        LinearBackend linearBackend
    ) {
        return makeContext(
            std::move(model), std::move(discretization), std::move(preconditionerPrescription),
            std::move(linearBackend), seed::LaneEmden{}
        );
    }

    template <
        typename Model,
        equilibrium::StellarDiscretizationType Discretization,
        preconditioning::StellarPreconditionerPrescription PreconditionerPrescription,
        LinearBackendConfiguration LinearBackend>
        requires DefaultStellarEquilibriumInitialStateAvailableFor<Model> && StellarEquilibriumContextConfiguration<
                                                                                 Model,
                                                                                 Discretization,
                                                                                 PreconditionerPrescription,
                                                                                 LinearBackend>
    [[nodiscard]] auto makeContext(
        Model model,
        Discretization discretization,
        PreconditionerPrescription preconditionerPrescription,
        LinearBackend linearBackend,
        seed::StellarEquilibriumProjectionOptions projectionOptions
    ) {
        return makeContext(
            std::move(model), std::move(discretization), std::move(preconditionerPrescription),
            std::move(linearBackend), seed::LaneEmden{}, std::move(projectionOptions)
        );
    }

    template <
        typename Model,
        equilibrium::StellarDiscretizationType Discretization,
        preconditioning::StellarPreconditionerPrescription PreconditionerPrescription,
        LinearBackendConfiguration LinearBackend>
        requires DefaultStellarEquilibriumInitialStateAvailableFor<Model> &&
                 StellarEquilibriumContextConfiguration<
                     Model,
                     Discretization,
                     PreconditionerPrescription,
                     LinearBackend> &&
                 (detail::StellarContextProblem<
                      Model,
                      Discretization>::generatedRotationProviderCount == 0)
    [[nodiscard]] auto makeContext(
        Model model,
        Discretization discretization,
        PreconditionerPrescription preconditionerPrescription,
        LinearBackend linearBackend,
        physics::RigidRotation prescribedRotation
    ) {
        return makeContext(
            std::move(model), std::move(discretization), std::move(preconditionerPrescription),
            std::move(linearBackend), seed::LaneEmden{}, std::move(prescribedRotation)
        );
    }

    template <
        typename Model,
        equilibrium::StellarDiscretizationType Discretization,
        preconditioning::StellarPreconditionerPrescription PreconditionerPrescription,
        LinearBackendConfiguration LinearBackend>
        requires DefaultStellarEquilibriumInitialStateAvailableFor<Model> &&
                 StellarEquilibriumContextConfiguration<
                     Model,
                     Discretization,
                     PreconditionerPrescription,
                     LinearBackend> &&
                 (detail::StellarContextProblem<
                      Model,
                      Discretization>::generatedRotationProviderCount == 0)
    [[nodiscard]] auto makeContext(
        Model model,
        Discretization discretization,
        PreconditionerPrescription preconditionerPrescription,
        LinearBackend linearBackend,
        seed::StellarEquilibriumProjectionOptions projectionOptions,
        physics::RigidRotation prescribedRotation
    ) {
        return makeContext(
            std::move(model), std::move(discretization), std::move(preconditionerPrescription),
            std::move(linearBackend), seed::LaneEmden{}, std::move(projectionOptions), std::move(prescribedRotation)
        );
    }
} // namespace mean_field::solver

export namespace mean_field::solver {
    template <typename Context, typename NewtonConfiguration, typename Observer> class StellarEquilibriumSolver final {
    private:
        using Clock = std::chrono::steady_clock;

        struct IterationTimings final {
            Clock::time_point start{};
            double geometryPreflightSeconds{0.0};
            double lineSearchSeconds{0.0};
            double trialPreparationSeconds{0.0};
            double metricEvaluationSeconds{0.0};
            double preconditionerRefreshSeconds{0.0};
            double rollbackSeconds{0.0};
            double observerSeconds{0.0};
            std::optional<deformation::LargestSafeNewtonStepSizeEstimate> geometryPreflight;
        };

    public:
        static_assert(StellarEquilibriumContextType<Context>);
        static_assert(nonlinear::NewtonConfiguration<std::remove_cvref_t<NewtonConfiguration>>);
        static_assert(nonlinear::NewtonObserver<std::remove_cvref_t<Observer>>);
        static_assert(std::move_constructible<std::remove_cvref_t<Observer>>);

        using ContextType      = std::remove_cvref_t<Context>;
        using NewtonType       = std::remove_cvref_t<NewtonConfiguration>;
        using ObserverType     = std::remove_cvref_t<Observer>;
        using ProblemType      = typename ContextType::ProblemType;
        using EvaluationReport = StellarEquilibriumEvaluationReport<ProblemType>;

        StellarEquilibriumSolver(
            ContextType &context,
            NewtonType newtonConfiguration,
            ObserverType observer
        )
            : m_context(std::addressof(context)),
              m_newton(std::move(newtonConfiguration)),
              m_observer(std::move(observer)) {
            m_newton.options().Validate();
            m_context->AcquireSolver();
        }

        StellarEquilibriumSolver(const StellarEquilibriumSolver &)            = delete;
        StellarEquilibriumSolver &operator=(const StellarEquilibriumSolver &) = delete;
        StellarEquilibriumSolver(StellarEquilibriumSolver &&)                 = delete;
        StellarEquilibriumSolver &operator=(StellarEquilibriumSolver &&)      = delete;

        ~StellarEquilibriumSolver() {
            if (m_context != nullptr) {
                m_context->ReleaseSolver();
            }
        }

        [[nodiscard]] bool isReady() const {
            return m_context != nullptr && m_context->hasActiveSolver() && m_context->isReady() && !m_evaluating;
        }

        [[nodiscard]] EvaluationReport evaluate() {
            EvaluationGuard evaluationGuard{m_evaluating};
            auto &state                 = *m_context->m_state;
            const MPI_Comm communicator = state.Problem().GetCommunicator();
            if (!detail::AllStellarRanksSatisfy(state.IsReady(), communicator)) {
                throw std::logic_error(
                    "A Newton evaluation requires a complete, current context on every communicator rank."
                );
            }

            const auto &options = m_newton.options();
            options.Validate();
            detail::RequireCollectivelyIdenticalNewtonOptions(options, communicator);
            state.BeginEvaluation();

            StellarEquilibriumEvaluationDiagnostics diagnostics{};
            nonlinear::MetricEvaluation acceptedMetric = EvaluateMetric(state.acceptedNormalizedResidual);
            diagnostics.initialResidualNorm            = acceptedMetric.residualNorm;
            diagnostics.finalResidualNorm              = acceptedMetric.residualNorm;

            if (!detail::AllStellarRanksSatisfy(ValidMetric(acceptedMetric), communicator)) {
                return Failure(
                    state, std::move(diagnostics), StellarEquilibriumFailureReason::non_finite_residual,
                    "The initial normalized residual metric is non-finite on at least one rank."
                );
            }

            const double convergenceThreshold = options.ConvergenceThreshold(acceptedMetric.residualNorm);
            if (detail::AllStellarRanksSatisfy(acceptedMetric.residualNorm <= convergenceThreshold, communicator)) {
                return Success(state, std::move(diagnostics));
            }

            double nextLineSearchStepLength = options.backtracking.initialStepLength;

            for (int iteration = 0; iteration < options.maximumIterations; ++iteration) {
                ++diagnostics.attemptedNonlinearIterations;
                NotifyBefore(iteration, diagnostics.initialResidualNorm, acceptedMetric, state);
                // Observer work is deliberately outside solver timing. Trial
                // callbacks are subtracted below because they occur inside the
                // line-search control flow.
                IterationTimings timings{.start = Clock::now()};

                state.linearRightHandSide = state.acceptedNormalizedResidual;
                state.linearRightHandSide *= -1.0;
                LinearSolveReport linearReport = state.linearBackend->Solve(
                    state.linearRightHandSide, state.normalizedCorrection, options.linearSolve
                );
                diagnostics.lastLinearSolve = linearReport;
                diagnostics.totalLinearSolveSeconds += linearReport.solveSeconds;

                const int expectedCorrectionSize = state.acceptedNormalizedState.Size();
                if (!detail::AllStellarRanksSatisfy(
                        state.normalizedCorrection.Size() == expectedCorrectionSize, communicator
                    )) {
                    state.normalizedCorrection.SetSize(expectedCorrectionSize);
                    state.normalizedCorrection = 0.0;
                    throw std::logic_error(
                        "A prepared linear backend changed the Newton correction vector's required size."
                    );
                }

                const bool correctionIsFinite =
                    detail::AllStellarRanksSatisfy(AllFinite(state.normalizedCorrection), communicator);
                const bool linearSolveConverged =
                    detail::AllStellarRanksSatisfy(linearReport.Converged(), communicator);
                if (!correctionIsFinite) {
                    state.normalizedCorrection = 0.0;
                }
                if (!linearSolveConverged || !correctionIsFinite) {
                    NotifyAfter(
                        iteration, nonlinear::IterationDisposition::linear_solve_failure, false, 0.0, 0,
                        diagnostics.initialResidualNorm, acceptedMetric, acceptedMetric, linearReport, timings, state
                    );
                    return Failure(
                        state, std::move(diagnostics), StellarEquilibriumFailureReason::linear_solve_failure,
                        "The linearized Newton step did not converge to a finite correction."
                    );
                }

                const nonlinear::MetricEvaluation previousMetric = acceptedMetric;
                const Clock::time_point geometryPreflightStart   = Clock::now();
                timings.geometryPreflight                        = state.EstimateLargestSafeStepSize(
                    nextLineSearchStepLength, options.backtracking.fractionToBoundarySafety
                );
                timings.geometryPreflightSeconds =
                    std::chrono::duration<double>(Clock::now() - geometryPreflightStart).count();
                diagnostics.totalGeometryPreflightSeconds += timings.geometryPreflightSeconds;
                diagnostics.lastGeometryPreflight = timings.geometryPreflight;
                if (timings.geometryPreflight->limitedByGeometry) {
                    ++diagnostics.geometryLimitedIterations;
                }

                if (timings.geometryPreflight->stepSize < options.backtracking.minimumStepLength) {
                    diagnostics.finalResidualNorm = acceptedMetric.residualNorm;
                    NotifyAfter(
                        iteration, nonlinear::IterationDisposition::globalization_failure, false, 0.0, 0,
                        diagnostics.initialResidualNorm, previousMetric, acceptedMetric, linearReport, timings, state
                    );
                    return Failure(
                        state, std::move(diagnostics), StellarEquilibriumFailureReason::globalization_failure,
                        "The geometry preflight found no orientation-preserving Newton step at or above the "
                        "configured minimum step length."
                    );
                }

                bool accepted             = false;
                double acceptedStepLength = 0.0;
                int lineSearchTrials      = 0;
                nonlinear::MetricEvaluation trialMetric{};
                StellarEquilibriumFailureReason rejectionReason =
                    StellarEquilibriumFailureReason::globalization_failure;
                std::string rejectionMessage = "The backtracking line search found no acceptable Newton step.";
                double stepLength            = timings.geometryPreflight->stepSize;

                const Clock::time_point lineSearchStart = Clock::now();
                try {
                    for (int trial = 0; trial < options.backtracking.maximumTrials; ++trial) {
                        ++lineSearchTrials;
                        ++diagnostics.totalLineSearchTrials;
                        nonlinear::LineSearchTrialDisposition trialDisposition =
                            nonlinear::LineSearchTrialDisposition::insufficient_decrease;
                        std::optional<nonlinear::MetricEvaluation> observedTrialMetric;
                        std::optional<double> minimumJacobianDeterminant;
                        std::string_view rejectionSource;
                        double preparationSeconds      = 0.0;
                        double metricSeconds           = 0.0;
                        bool residualAvailable         = false;
                        bool invertedGeometryRejection = false;

                        state.trialNormalizedState     = state.acceptedNormalizedState;
                        state.trialNormalizedState.Add(stepLength, state.normalizedCorrection);
                        state.normalizedOperator->DenormalizeState(
                            state.trialNormalizedState, state.candidatePhysicalState
                        );

                        if (!detail::AllStellarRanksSatisfy(
                                AllFinite(state.trialNormalizedState) && AllFinite(state.candidatePhysicalState),
                                communicator
                            )) {
                            rejectionReason = StellarEquilibriumFailureReason::non_finite_state;
                            rejectionMessage =
                                "The backtracking line search produced a non-finite candidate state on at least one "
                                "rank.";
                            trialDisposition = nonlinear::LineSearchTrialDisposition::non_finite_state;
                            rejectionSource  = "candidate state";
                            ++diagnostics.nonFiniteLineSearchTrials;
                        } else {
                            const Clock::time_point preparationStart = Clock::now();
                            auto preparationResult                   = state.PrepareTrial();
                            preparationSeconds = std::chrono::duration<double>(Clock::now() - preparationStart).count();
                            timings.trialPreparationSeconds += preparationSeconds;
                            diagnostics.totalTrialPreparationSeconds += preparationSeconds;

                            if (!preparationResult.has_value()) {
                                const auto &preparationRejection = preparationResult.error();
                                rejectionSource = detail::PreparationStageName(preparationRejection.stage);
                                invertedGeometryRejection =
                                    preparationRejection.reason ==
                                    operators::StellarEquilibriumPreparationRejectionReason::inverted_geometry;
                                switch (preparationRejection.reason) {
                                case operators::StellarEquilibriumPreparationRejectionReason::inverted_geometry:
                                case operators::StellarEquilibriumPreparationRejectionReason::thermodynamic_domain:
                                    rejectionReason = StellarEquilibriumFailureReason::inadmissible_state;
                                    rejectionMessage =
                                        "The candidate Newton state was physically inadmissible on at least one rank.";
                                    trialDisposition = nonlinear::LineSearchTrialDisposition::inadmissible_state;
                                    ++diagnostics.inadmissibleLineSearchTrials;
                                    break;
                                case operators::StellarEquilibriumPreparationRejectionReason::non_finite_geometry:
                                    rejectionReason = StellarEquilibriumFailureReason::non_finite_state;
                                    rejectionMessage =
                                        "The candidate Newton state produced non-finite mapped geometry on at least "
                                        "one rank.";
                                    trialDisposition = nonlinear::LineSearchTrialDisposition::non_finite_state;
                                    ++diagnostics.nonFiniteLineSearchTrials;
                                    break;
                                case operators::StellarEquilibriumPreparationRejectionReason::non_finite_thermodynamics:
                                    if (preparationRejection.thermodynamicErrorCode ==
                                        eos::EvaluationErrorCode::nonfinite_input) {
                                        rejectionReason = StellarEquilibriumFailureReason::non_finite_state;
                                        rejectionMessage =
                                            "The candidate Newton state produced non-finite thermodynamic input "
                                            "on at least one rank.";
                                        trialDisposition = nonlinear::LineSearchTrialDisposition::non_finite_state;
                                    } else {
                                        rejectionReason = StellarEquilibriumFailureReason::non_finite_residual;
                                        rejectionMessage =
                                            "The candidate Newton state produced non-finite thermodynamic data on "
                                            "at least one rank.";
                                        trialDisposition = nonlinear::LineSearchTrialDisposition::non_finite_residual;
                                    }
                                    ++diagnostics.nonFiniteLineSearchTrials;
                                    break;
                                case operators::StellarEquilibriumPreparationRejectionReason::inadmissible_physics:
                                    rejectionReason = StellarEquilibriumFailureReason::inadmissible_state;
                                    rejectionMessage =
                                        "The candidate Newton state violated a physical admissibility condition on "
                                        "at least one rank.";
                                    trialDisposition = nonlinear::LineSearchTrialDisposition::inadmissible_state;
                                    ++diagnostics.inadmissibleLineSearchTrials;
                                    break;
                                case operators::StellarEquilibriumPreparationRejectionReason::non_finite_physics:
                                    rejectionReason = StellarEquilibriumFailureReason::non_finite_residual;
                                    rejectionMessage =
                                        "The candidate Newton state produced non-finite physical data on at least "
                                        "one rank.";
                                    trialDisposition = nonlinear::LineSearchTrialDisposition::non_finite_residual;
                                    ++diagnostics.nonFiniteLineSearchTrials;
                                    break;
                                default:
                                    throw std::logic_error(
                                        "The prepared stellar operator returned an unknown candidate-rejection "
                                        "reason."
                                    );
                                }
                                if (preparationRejection.reason ==
                                        operators::StellarEquilibriumPreparationRejectionReason::inverted_geometry ||
                                    preparationRejection.reason ==
                                        operators::StellarEquilibriumPreparationRejectionReason::non_finite_geometry) {
                                    if (std::isfinite(preparationRejection.minimumJacobianDeterminant)) {
                                        minimumJacobianDeterminant = preparationRejection.minimumJacobianDeterminant;
                                    }
                                }
                            } else if (!detail::AllStellarRanksSatisfy(
                                           AllFinite(state.trialNormalizedResidual), communicator
                                       )) {
                                residualAvailable = true;
                                rejectionReason   = StellarEquilibriumFailureReason::non_finite_residual;
                                rejectionMessage =
                                    "The backtracking line search produced a non-finite normalized residual on at "
                                    "least one rank.";
                                trialDisposition           = nonlinear::LineSearchTrialDisposition::non_finite_residual;
                                rejectionSource            = "normalized residual";
                                minimumJacobianDeterminant = preparationResult->minimumJacobianDeterminant;
                                ++diagnostics.nonFiniteLineSearchTrials;
                            } else {
                                residualAvailable                   = true;
                                minimumJacobianDeterminant          = preparationResult->minimumJacobianDeterminant;
                                const Clock::time_point metricStart = Clock::now();
                                trialMetric                         = EvaluateMetric(state.trialNormalizedResidual);
                                metricSeconds = std::chrono::duration<double>(Clock::now() - metricStart).count();
                                timings.metricEvaluationSeconds += metricSeconds;
                                diagnostics.totalMetricEvaluationSeconds += metricSeconds;
                                observedTrialMetric = trialMetric;
                                if (!detail::AllStellarRanksSatisfy(ValidMetric(trialMetric), communicator)) {
                                    rejectionReason = StellarEquilibriumFailureReason::non_finite_residual;
                                    rejectionMessage =
                                        "The backtracking line search produced a non-finite residual metric on at "
                                        "least one rank.";
                                    trialDisposition = nonlinear::LineSearchTrialDisposition::non_finite_residual;
                                    rejectionSource  = "residual metric";
                                    ++diagnostics.nonFiniteLineSearchTrials;
                                } else if (
                                    detail::AllStellarRanksSatisfy(
                                        AcceptStep(acceptedMetric, trialMetric, stepLength, options.backtracking),
                                        communicator
                                    )
                                ) {
                                    const Clock::time_point refreshStart = Clock::now();
                                    int localCommitStatus                = 0;
                                    std::exception_ptr commitFailure;
                                    try {
                                        state.PrepareTrialCommit();
                                    } catch (...) {
                                        localCommitStatus = 1;
                                        commitFailure     = std::current_exception();
                                    }
                                    if (detail::MaximumStellarRankValue(localCommitStatus, communicator) != 0) {
                                        if (commitFailure != nullptr) {
                                            std::rethrow_exception(commitFailure);
                                        }
                                        throw std::runtime_error(
                                            "Preparing to commit a Newton step failed on another rank."
                                        );
                                    }
                                    const double refreshSeconds =
                                        std::chrono::duration<double>(Clock::now() - refreshStart).count();
                                    timings.preconditionerRefreshSeconds += refreshSeconds;
                                    diagnostics.totalPreconditionerRefreshSeconds += refreshSeconds;

                                    accepted           = true;
                                    acceptedStepLength = stepLength;
                                    trialDisposition   = nonlinear::LineSearchTrialDisposition::accepted;
                                } else {
                                    rejectionReason = StellarEquilibriumFailureReason::globalization_failure;
                                    rejectionMessage =
                                        "The candidate Newton step did not provide sufficient metric decrease on "
                                        "every rank.";
                                    trialDisposition = nonlinear::LineSearchTrialDisposition::insufficient_decrease;
                                    rejectionSource  = "residual metric";
                                    ++diagnostics.insufficientDecreaseTrials;
                                }
                            }
                        }

                        timings.observerSeconds += NotifyLineSearchTrial(
                            iteration, trial, stepLength, trialDisposition, rejectionSource, observedTrialMetric,
                            minimumJacobianDeterminant, preparationSeconds, metricSeconds, residualAvailable, state
                        );

                        if (accepted) {
                            if (!minimumJacobianDeterminant.has_value() ||
                                !std::isfinite(*minimumJacobianDeterminant) || *minimumJacobianDeterminant <= 0.0) {
                                throw std::logic_error(
                                    "An accepted Newton trial did not retain a valid geometry certificate."
                                );
                            }
                            state.CommitPreparedTrial(*minimumJacobianDeterminant);
                            state.AdvanceCorrectionWarmStart(stepLength);
                            break;
                        }

                        stepLength = nonlinear::detail::NextBacktrackingStepLength(
                            stepLength, state.acceptedMinimumJacobianDeterminant, minimumJacobianDeterminant,
                            invertedGeometryRejection, options.backtracking
                        );
                        if (stepLength < options.backtracking.minimumStepLength) {
                            break;
                        }
                    }
                } catch (...) {
                    timings.lineSearchSeconds = DurationExcludingObserver(lineSearchStart, timings.observerSeconds);
                    state.RestoreAccepted();
                    throw;
                }
                timings.lineSearchSeconds = DurationExcludingObserver(lineSearchStart, timings.observerSeconds);
                diagnostics.totalLineSearchSeconds += timings.lineSearchSeconds;

                if (!accepted) {
                    const Clock::time_point rollbackStart = Clock::now();
                    state.RestoreAccepted();
                    timings.rollbackSeconds = std::chrono::duration<double>(Clock::now() - rollbackStart).count();
                    diagnostics.totalRollbackSeconds += timings.rollbackSeconds;
                    diagnostics.finalResidualNorm = acceptedMetric.residualNorm;
                    NotifyAfter(
                        iteration, DispositionFor(rejectionReason), false, 0.0, lineSearchTrials,
                        diagnostics.initialResidualNorm, previousMetric, acceptedMetric, linearReport, timings, state
                    );
                    return Failure(state, std::move(diagnostics), rejectionReason, std::move(rejectionMessage));
                }

                acceptedMetric           = trialMetric;
                nextLineSearchStepLength = ExpandedLineSearchStepLength(acceptedStepLength, options.backtracking);
                ++diagnostics.acceptedNonlinearIterations;
                diagnostics.finalResidualNorm      = acceptedMetric.residualNorm;
                diagnostics.lastAcceptedStepLength = acceptedStepLength;
                const bool converged =
                    detail::AllStellarRanksSatisfy(acceptedMetric.residualNorm <= convergenceThreshold, communicator);
                const bool reachedIterationLimit = !converged && iteration + 1 == options.maximumIterations;
                NotifyAfter(
                    iteration,
                    converged               ? nonlinear::IterationDisposition::converged
                    : reachedIterationLimit ? nonlinear::IterationDisposition::iteration_limit
                                            : nonlinear::IterationDisposition::accepted,
                    true, acceptedStepLength, lineSearchTrials, diagnostics.initialResidualNorm, previousMetric,
                    acceptedMetric, linearReport, timings, state
                );

                if (converged) {
                    return Success(state, std::move(diagnostics));
                }
                if (reachedIterationLimit) {
                    return Failure(
                        state, std::move(diagnostics), StellarEquilibriumFailureReason::iteration_limit,
                        "The damped Newton method reached its nonlinear iteration limit."
                    );
                }
            }

            return Failure(
                state, std::move(diagnostics), StellarEquilibriumFailureReason::iteration_limit,
                "The damped Newton method reached its nonlinear iteration limit."
            );
        }

    private:
        using State = typename ContextType::State;

        class EvaluationGuard final {
        public:
            explicit EvaluationGuard(bool &active) : m_active(active) {
                if (m_active) {
                    throw std::logic_error("A stellar-equilibrium solver cannot evaluate reentrantly.");
                }
                m_active = true;
            }

            ~EvaluationGuard() {
                m_active = false;
            }

            EvaluationGuard(const EvaluationGuard &)            = delete;
            EvaluationGuard &operator=(const EvaluationGuard &) = delete;

        private:
            bool &m_active;
        };

        [[nodiscard]] nonlinear::MetricEvaluation EvaluateMetric(const mfem::Vector &residual) const {
            using nonlinear::getMetric;
            return getMetric(m_newton.metric(), residual, m_context->communicator());
        }

        [[nodiscard]] static bool ValidMetric(const nonlinear::MetricEvaluation &metric) noexcept {
            return std::isfinite(metric.residualNorm) && metric.residualNorm >= 0.0 && std::isfinite(metric.merit) &&
                   metric.merit >= 0.0;
        }

        [[nodiscard]] static bool AllFinite(const mfem::Vector &values) noexcept {
            for (int index = 0; index < values.Size(); ++index) {
                if (!std::isfinite(values(index))) {
                    return false;
                }
            }
            return true;
        }

        [[nodiscard]] static double RelativeResidual(
            const double residualNorm,
            const double initialResidualNorm
        ) noexcept {
            if (initialResidualNorm > 0.0) {
                return residualNorm / initialResidualNorm;
            }
            return residualNorm == 0.0 ? 0.0 : std::numeric_limits<double>::infinity();
        }

        [[nodiscard]] static double DurationExcludingObserver(
            const Clock::time_point start,
            const double observerSeconds
        ) noexcept {
            const double elapsed = std::chrono::duration<double>(Clock::now() - start).count();
            return elapsed > observerSeconds ? elapsed - observerSeconds : 0.0;
        }

        [[nodiscard]] static bool AcceptStep(
            const nonlinear::MetricEvaluation &accepted,
            const nonlinear::MetricEvaluation &trial,
            const double stepLength,
            const nonlinear::BacktrackingOptions &options
        ) noexcept {
            return trial.merit <= (1.0 - options.sufficientDecrease * stepLength) * accepted.merit;
        }

        [[nodiscard]] static double ExpandedLineSearchStepLength(
            const double acceptedStepLength,
            const nonlinear::BacktrackingOptions &options
        ) noexcept {
            const double expandedStepLength = acceptedStepLength / options.contractionFactor;
            return expandedStepLength < options.initialStepLength ? expandedStepLength : options.initialStepLength;
        }

        [[nodiscard]] static nonlinear::IterationDisposition
        DispositionFor(const StellarEquilibriumFailureReason reason) noexcept {
            switch (reason) {
            case StellarEquilibriumFailureReason::inadmissible_state:
                return nonlinear::IterationDisposition::inadmissible_state;
            case StellarEquilibriumFailureReason::non_finite_state:
                return nonlinear::IterationDisposition::non_finite_state;
            case StellarEquilibriumFailureReason::non_finite_residual:
                return nonlinear::IterationDisposition::non_finite_residual;
            case StellarEquilibriumFailureReason::linear_solve_failure:
                return nonlinear::IterationDisposition::linear_solve_failure;
            case StellarEquilibriumFailureReason::stagnation:
                return nonlinear::IterationDisposition::stagnation;
            case StellarEquilibriumFailureReason::iteration_limit:
                return nonlinear::IterationDisposition::iteration_limit;
            default:
                return nonlinear::IterationDisposition::globalization_failure;
            }
        }

        void NotifyBefore(
            const int iteration,
            const double initialResidualNorm,
            const nonlinear::MetricEvaluation &metric,
            const State &state
        ) {
            if constexpr (nonlinear::detail::ObservesBeforeIteration<ObserverType>) {
                const nonlinear::BeforeIteration event{
                    .iteration            = iteration,
                    .initialResidualNorm  = initialResidualNorm,
                    .residualNorm         = metric.residualNorm,
                    .relativeResidualNorm = RelativeResidual(metric.residualNorm, initialResidualNorm),
                    .merit                = metric.merit,
                    .communicator         = state.Problem().GetCommunicator(),
                    .physicalState        = detail::ReadOnlySpan(state.AcceptedPhysicalState()),
                    .normalizedState      = detail::ReadOnlySpan(state.acceptedNormalizedState),
                    .normalizedResidual   = detail::ReadOnlySpan(state.acceptedNormalizedResidual)
                };
                nonlinear::detail::InvokeBeforeIteration(m_observer, event);
            }
        }

        [[nodiscard]] double NotifyLineSearchTrial(
            const int iteration,
            const int trial,
            const double stepLength,
            const nonlinear::LineSearchTrialDisposition disposition,
            const std::string_view rejectionSource,
            const std::optional<nonlinear::MetricEvaluation> &metric,
            const std::optional<double> minimumJacobianDeterminant,
            const double preparationSeconds,
            const double metricSeconds,
            const bool residualAvailable,
            const State &state
        ) {
            if constexpr (nonlinear::detail::ObservesLineSearchTrial<ObserverType>) {
                const nonlinear::AfterLineSearchTrial event{
                    .iteration                   = iteration,
                    .trial                       = trial,
                    .stepLength                  = stepLength,
                    .disposition                 = disposition,
                    .rejectionSource             = rejectionSource,
                    .metric                      = metric,
                    .minimumJacobianDeterminant  = minimumJacobianDeterminant,
                    .preparationSeconds          = preparationSeconds,
                    .metricSeconds               = metricSeconds,
                    .communicator                = state.Problem().GetCommunicator(),
                    .candidatePhysicalState      = detail::ReadOnlySpan(state.candidatePhysicalState),
                    .candidateNormalizedState    = detail::ReadOnlySpan(state.trialNormalizedState),
                    .candidateNormalizedResidual = residualAvailable
                                                       ? detail::ReadOnlySpan(state.trialNormalizedResidual)
                                                       : std::span<const mfem::real_t>{}
                };
                const Clock::time_point observerStart = Clock::now();
                nonlinear::detail::InvokeAfterLineSearchTrial(m_observer, event);
                return std::chrono::duration<double>(Clock::now() - observerStart).count();
            }
            return 0.0;
        }

        void NotifyAfter(
            const int iteration,
            const nonlinear::IterationDisposition disposition,
            const bool stepAccepted,
            const double acceptedStepLength,
            const int lineSearchTrials,
            const double initialResidualNorm,
            const nonlinear::MetricEvaluation &previousMetric,
            const nonlinear::MetricEvaluation &metric,
            const LinearSolveReport &linearReport,
            const IterationTimings &timings,
            const State &state
        ) {
            if constexpr (nonlinear::detail::ObservesAfterIteration<ObserverType>) {
                const nonlinear::AfterIteration event{
                    .iteration                    = iteration,
                    .disposition                  = disposition,
                    .stepAccepted                 = stepAccepted,
                    .acceptedStepLength           = acceptedStepLength,
                    .lineSearchTrials             = lineSearchTrials,
                    .initialResidualNorm          = initialResidualNorm,
                    .previousResidualNorm         = previousMetric.residualNorm,
                    .residualNorm                 = metric.residualNorm,
                    .relativeResidualNorm         = RelativeResidual(metric.residualNorm, initialResidualNorm),
                    .merit                        = metric.merit,
                    .iterationSeconds             = DurationExcludingObserver(timings.start, timings.observerSeconds),
                    .geometryPreflightSeconds     = timings.geometryPreflightSeconds,
                    .lineSearchSeconds            = timings.lineSearchSeconds,
                    .trialPreparationSeconds      = timings.trialPreparationSeconds,
                    .metricEvaluationSeconds      = timings.metricEvaluationSeconds,
                    .preconditionerRefreshSeconds = timings.preconditionerRefreshSeconds,
                    .rollbackSeconds              = timings.rollbackSeconds,
                    .geometryPreflight            = timings.geometryPreflight,
                    .linearSolve                  = linearReport,
                    .communicator                 = state.Problem().GetCommunicator(),
                    .physicalState                = detail::ReadOnlySpan(state.AcceptedPhysicalState()),
                    .normalizedState              = detail::ReadOnlySpan(state.acceptedNormalizedState),
                    .normalizedResidual           = detail::ReadOnlySpan(state.acceptedNormalizedResidual)
                };
                nonlinear::detail::InvokeAfterIteration(m_observer, event);
            }
        }

        [[nodiscard]] EvaluationReport Success(
            State &state,
            StellarEquilibriumEvaluationDiagnostics diagnostics
        ) {
            return detail::StellarEvaluationReportAccess<ProblemType>::Success(state.storage, std::move(diagnostics));
        }

        [[nodiscard]] EvaluationReport Failure(
            State &state,
            StellarEquilibriumEvaluationDiagnostics diagnostics,
            const StellarEquilibriumFailureReason reason,
            std::string message
        ) {
            return detail::StellarEvaluationReportAccess<ProblemType>::Failure(
                state.storage, std::move(diagnostics), reason, std::move(message)
            );
        }

        ContextType *m_context;
        NewtonType m_newton;
        ObserverType m_observer;
        bool m_evaluating{false};
    };

    template <
        StellarEquilibriumContextType Context,
        typename NewtonConfiguration,
        typename Observer>
        requires nonlinear::NewtonConfiguration<std::remove_cvref_t<NewtonConfiguration>> &&
                 nonlinear::NewtonObserver<std::remove_cvref_t<Observer>> &&
                 std::move_constructible<std::remove_cvref_t<Observer>>
    [[nodiscard]] auto make(
        Context &context,
        NewtonConfiguration newtonConfiguration,
        Observer observer
    ) {
        using Solver = StellarEquilibriumSolver<
            std::remove_cvref_t<Context>, std::remove_cvref_t<NewtonConfiguration>, std::remove_cvref_t<Observer>>;
        return Solver{context, std::move(newtonConfiguration), std::move(observer)};
    }

    template <
        StellarEquilibriumContextType Context,
        typename NewtonConfiguration>
        requires nonlinear::NewtonConfiguration<std::remove_cvref_t<NewtonConfiguration>>
    [[nodiscard]] auto make(
        Context &context,
        NewtonConfiguration newtonConfiguration
    ) {
        using Solver = StellarEquilibriumSolver<
            std::remove_cvref_t<Context>, std::remove_cvref_t<NewtonConfiguration>, nonlinear::NoObserver>;
        return Solver{context, std::move(newtonConfiguration), nonlinear::NoObserver{}};
    }

    template <
        StellarEquilibriumContextType Context,
        typename NewtonConfiguration,
        typename Observer>
    auto make(
        Context &&,
        NewtonConfiguration,
        Observer
    ) = delete;

    template <
        StellarEquilibriumContextType Context,
        typename NewtonConfiguration>
    auto make(
        Context &&,
        NewtonConfiguration
    ) = delete;
} // namespace mean_field::solver