From c393847c6f6d7730d6964632b9cd506384ed0bf6 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Thu, 17 Sep 2026 09:45:05 -0400 Subject: [PATCH] Integrate typed Newton strategies with checked trial gravity and feasibility --- CMakeLists.txt | 10 + experiments/plateau_diagnostic.cpp | 77 +++- libmeanfield/interface/mean_field.cppm | 7 + libmeanfield/interface/solver/newton.cppm | 21 +- .../interface/solver/newton_numerics.cppm | 83 ++++ .../interface/solver/stellar_context.cppm | 68 ++- .../solver/stellar_equilibrium_types.cppm | 2 + .../solver/stellar_newton_strategy.cppm | 396 ++++++++++++++++++ 8 files changed, 643 insertions(+), 21 deletions(-) create mode 100644 libmeanfield/interface/solver/newton_numerics.cppm create mode 100644 libmeanfield/interface/solver/stellar_newton_strategy.cppm diff --git a/CMakeLists.txt b/CMakeLists.txt index 7f1d9ca..50969f1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -152,6 +152,12 @@ target_sources(mean_field libmeanfield/interface/solver/fields.cppm libmeanfield/interface/solver/linear_backend.cppm libmeanfield/interface/solver/newton.cppm + libmeanfield/interface/solver/enthalpy_admissibility.cppm + libmeanfield/interface/solver/gravity_subsystem.cppm + libmeanfield/interface/solver/gravity_reduction.cppm + libmeanfield/interface/solver/stellar_enthalpy_strategy.cppm + libmeanfield/interface/solver/newton_numerics.cppm + libmeanfield/interface/solver/stellar_newton_strategy.cppm libmeanfield/interface/solver/preconditioning_diagnostics.cppm libmeanfield/interface/solver/stellar_equilibrium_types.cppm libmeanfield/interface/solver/stellar_structure.cppm @@ -365,6 +371,10 @@ add_executable(tests tests/solver/stellar_equilibrium_architecture.cpp tests/solver/stellar_equilibrium_architecture_internal.cpp tests/solver/stellar_equilibrium_runtime.cpp + tests/solver/enthalpy_admissibility.cpp + tests/solver/gravity_subsystem.cpp + tests/solver/gravity_reduction.cpp + tests/solver/stellar_numerical_strategies.cpp ) target_link_libraries(tests PRIVATE mean_field test_mod Catch2::Catch2 Boost::boost) diff --git a/experiments/plateau_diagnostic.cpp b/experiments/plateau_diagnostic.cpp index bb6bf58..66b40f6 100644 --- a/experiments/plateau_diagnostic.cpp +++ b/experiments/plateau_diagnostic.cpp @@ -20,6 +20,8 @@ import mean_field; #include "plateau_preconditioner_diagnostics.hpp" #include "plateau_pressure_obstruction_diagnostics.hpp" #include "plateau_shape_curvature_diagnostics.hpp" +#include "plateau_numerical_strategy_diagnostics.hpp" +#include "plateau_positivity_strategy_diagnostics.hpp" using namespace mean_field; namespace { @@ -28,6 +30,10 @@ struct Options { double n = 1., angularMomentum = 0., atol = 0., rtol = 1e-8, linearTolerance = .03; int maxNewton = 12, maxLinear = 200, restart = 40, probeRestart = 40; bool probes = false, noCentral = false; + bool strategyProbe = false, positivityProbe = false; + std::string strategy = "baseline"; + double gravityTolerance = 1e-10; + int gravityIterations = 300, gravityRestart = 80; bool pressureProbe = false, shapeCurvature = false, surfaceTrace = false; std::vector captureIterations; std::string mesh = "sandbox.smesh"; @@ -70,6 +76,18 @@ Options Parse(int argc, char **argv) { o.probeDirection = v(); else if (a == "--replay-candidate") o.replayCandidate = v(); + else if (a == "--strategy") + o.strategy = v(); + else if (a == "--gravity-tol") + o.gravityTolerance = std::stod(v()); + else if (a == "--gravity-iterations") + o.gravityIterations = std::stoi(v()); + else if (a == "--gravity-restart") + o.gravityRestart = std::stoi(v()); + else if (a == "--strategy-probe") + o.strategyProbe = true; + else if (a == "--positivity-probe") + o.positivityProbe = true; else if (a == "--pressure-probe") o.pressureProbe = true; else if (a == "--shape-curvature") @@ -100,10 +118,15 @@ Options Parse(int argc, char **argv) { throw std::invalid_argument("Invalid options"); if (o.probeRestart != 40 && o.probeRestart != 160) throw std::invalid_argument("--probe-restart must be 40 or 160"); - if ((o.pressureProbe || o.shapeCurvature || o.surfaceTrace) && (o.probeState.empty() || o.probeDirection.empty())) + if ((o.pressureProbe || o.shapeCurvature || o.surfaceTrace || o.strategyProbe || o.positivityProbe) && (o.probeState.empty() || o.probeDirection.empty())) throw std::invalid_argument("Focused probes require --probe-state and --probe-direction"); if ((!o.probeDirection.empty() || !o.replayCandidate.empty()) && o.probeState.empty()) throw std::invalid_argument("A supplied direction or candidate requires --probe-state"); + if (o.strategy != "baseline" && o.strategy != "positivity" && o.strategy != "curvature" && + o.strategy != "reduced" && o.strategy != "positivity-curvature" && o.strategy != "positivity-reduced") + throw std::invalid_argument("Unknown numerical strategy"); + if (!std::isfinite(o.gravityTolerance) || o.gravityTolerance <= 0.0 || o.gravityIterations < 1 || o.gravityRestart < 1) + throw std::invalid_argument("Invalid gravity solve controls"); if (std::filesystem::exists(o.output)) throw std::runtime_error("Output directory already exists"); std::filesystem::create_directories(o.output); @@ -227,6 +250,10 @@ void FocusedProbes(State &s, const fem::FEM &f, const Options &o, const mfem::Ve for (const auto &[key, value] : experiment::plateau_diagnostics::MeasureEnthalpyStep(s, f, o.n, &physicalBase)) bounds << key << ',' << value << '\n'; bounds.flush(); + if (o.strategyProbe) + experiment::plateau_diagnostics::ProbeNumericalStrategies(s, f, base, direction, o.output); + if (o.positivityProbe) + experiment::plateau_diagnostics::ProbePolynomialEnthalpy(s, f, base, direction, o.output); if (o.surfaceTrace) experiment::plateau_diagnostics::ProbeSurfaceTrace(s, f, o.n, base, direction, o.output); if (o.pressureProbe) @@ -357,6 +384,8 @@ template int RunContext(Context &context, const Options &o, cons << "\nfrozen_probe_backend=" << (o.probeDirection.empty() ? "direct_mfem_fgmres_cold" : "supplied_direction_no_linear_solve") << "\nlinear_tol=" << o.linearTolerance << "\nmax_linear=" << o.maxLinear << "\nmax_newton=" << o.maxNewton << "\natol=" << o.atol << "\nrtol=" << o.rtol + << "\nstrategy=" << o.strategy << "\ngravity_tolerance=" << o.gravityTolerance + << "\ngravity_iterations=" << o.gravityIterations << "\ngravity_restart=" << o.gravityRestart << "\ncompiler=" << __VERSION__ << "\ncompiled=" << __DATE__ << ' ' << __TIME__ << '\n'; metadata.flush(); if (!o.probeState.empty()) { @@ -537,18 +566,44 @@ template int RunContext(Context &context, const Options &o, cons }); bool converged = false; if (o.maxNewton > 0) { - auto report = [&]() { - auto solve = solver::make(context, - solver::nonlinear::Newton(solver::nonlinear::NewtonOptions{ - .relativeTolerance = o.rtol, - .absoluteTolerance = o.atol, - .maximumIterations = o.maxNewton, - .linearSolve = {.relativeTolerance = o.linearTolerance, - .absoluteTolerance = 0., - .maximumIterations = o.maxLinear}}), - observer); + namespace nl = solver::nonlinear; + const nl::NewtonOptions options{ + .relativeTolerance = o.rtol, .absoluteTolerance = o.atol, .maximumIterations = o.maxNewton, + .linearSolve = {.relativeTolerance = o.linearTolerance, .absoluteTolerance = 0., .maximumIterations = o.maxLinear}}; + const nl::GravitySolveOptions gravity{ + .control = {.relativeTolerance = o.gravityTolerance, .absoluteTolerance = 1e-14, + .maximumIterations = o.gravityIterations}, + .backend = {.restartLength = o.gravityRestart}}; + auto run = [&](auto numerics) { + auto solve = solver::make(context, nl::Newton(options, nl::NormalizedResidualMetric{}, numerics), observer); return solve.evaluate(); + }; + auto report = [&] { + if (o.strategy == "positivity") return run(nl::NewtonNumerics{nl::PolynomialEnthalpy{}, nl::MonolithicGravity{}}); + if (o.strategy == "curvature") return run(nl::NewtonNumerics{nl::SampledEnthalpy{}, nl::GravityCurvatureCorrection{gravity}}); + if (o.strategy == "reduced") return run(nl::NewtonNumerics{nl::SampledEnthalpy{}, nl::EliminatedGravity{gravity}}); + if (o.strategy == "positivity-curvature") return run(nl::NewtonNumerics{nl::PolynomialEnthalpy{}, nl::GravityCurvatureCorrection{gravity}}); + if (o.strategy == "positivity-reduced") return run(nl::NewtonNumerics{nl::PolynomialEnthalpy{}, nl::EliminatedGravity{gravity}}); + return run(nl::NewtonNumerics{}); }(); + const auto &strategyDiagnostics = report.diagnostics().numericalStrategies; + auto strategyMetrics = File(o.output / "numerical_strategies.csv"); + strategyMetrics << "metric,value\n"; + const auto metric = [&](const char *name, auto value) { strategyMetrics << name << ',' << value << '\n'; }; + metric("polynomial_checks", strategyDiagnostics.polynomialChecks); + metric("direction_repairs", strategyDiagnostics.directionRepairs); + metric("last_repair_fraction", strategyDiagnostics.lastRepairFraction); + metric("last_merit_directional_derivative", strategyDiagnostics.lastMeritDirectionalDerivative); + metric("gravity_solves", strategyDiagnostics.gravitySolves); + metric("gravity_iterations", strategyDiagnostics.gravityIterations); + metric("gravity_failures", strategyDiagnostics.gravityFailures); + metric("gravity_solve_seconds", strategyDiagnostics.gravitySolveSeconds); + metric("last_gravity_true_residual", strategyDiagnostics.lastGravityTrueResidual); + metric("last_gravity_correction_norm", strategyDiagnostics.lastGravityCorrectionNorm); + metric("last_gravity_prepared_defect", strategyDiagnostics.lastGravityPreparedDefect); + metric("initial_gravity_residual", strategyDiagnostics.initialGravityResidual); + metric("completed_initial_gravity_residual", strategyDiagnostics.completedInitialGravityResidual); + metric("last_lifted_relative_true_residual", strategyDiagnostics.lastLiftedRelativeTrueResidual); converged = report.converged(); metadata << "converged=" << converged << "\naccepted_steps=" << report.completedNonlinearIterations() << '\n'; if (!converged) diff --git a/libmeanfield/interface/mean_field.cppm b/libmeanfield/interface/mean_field.cppm index da2daf7..2b02681 100644 --- a/libmeanfield/interface/mean_field.cppm +++ b/libmeanfield/interface/mean_field.cppm @@ -103,3 +103,10 @@ export import :equilibrium.stellar_discretization; export import :operators.stellar_equilibrium_problem; export import :seed.stellar_equilibrium_projection; export import :operators.stellar_equilibrium_system; + +export import :solver.enthalpy_admissibility; +export import :solver.gravity_subsystem; +export import :solver.gravity_reduction; + +export import :solver.stellar_enthalpy_strategy; +export import :solver.newton_numerics; diff --git a/libmeanfield/interface/solver/newton.cppm b/libmeanfield/interface/solver/newton.cppm index 0a5c199..97eaf7a 100644 --- a/libmeanfield/interface/solver/newton.cppm +++ b/libmeanfield/interface/solver/newton.cppm @@ -19,6 +19,7 @@ export module mean_field:solver.newton; export import :deformation.safe_newton_step; export import :solver.linear_backend; +export import :solver.newton_numerics; export namespace mean_field::solver::nonlinear { /* @@ -143,11 +144,12 @@ export namespace mean_field::solver::nonlinear { return {.residualNorm = residualNorm, .merit = 0.5 * residualNorm * residualNorm}; } - template + template > requires std::move_constructible> class Newton final { public: using MetricType = std::remove_cvref_t; + using NumericsType = std::remove_cvref_t; Newton() requires std::default_initializable @@ -167,10 +169,12 @@ export namespace mean_field::solver::nonlinear { Newton( NewtonOptions options, - MetricType metric + MetricType metric, + NumericsType numerics = {} ) : m_options(std::move(options)), - m_metric(std::move(metric)) { + m_metric(std::move(metric)), + m_numerics(std::move(numerics)) { m_options.Validate(); } @@ -182,9 +186,12 @@ export namespace mean_field::solver::nonlinear { return m_metric; } + [[nodiscard]] const NumericsType &numerics() const noexcept { return m_numerics; } + private: NewtonOptions m_options; [[no_unique_address]] MetricType m_metric; + [[no_unique_address]] NumericsType m_numerics; }; Newton() -> Newton; @@ -196,9 +203,12 @@ export namespace mean_field::solver::nonlinear { Metric ) -> Newton>; + template + Newton(NewtonOptions, Metric, Numerics) -> Newton, std::remove_cvref_t>; + template struct IsNewtonConfiguration : std::false_type { }; - template struct IsNewtonConfiguration> : std::true_type { }; + template struct IsNewtonConfiguration> : std::true_type { }; template concept NewtonConfiguration = IsNewtonConfiguration>::value; @@ -237,7 +247,8 @@ export namespace mean_field::solver::nonlinear { inadmissible_state, non_finite_state, non_finite_residual, - insufficient_decrease + insufficient_decrease, + auxiliary_solve_failure }; struct AfterLineSearchTrial final { diff --git a/libmeanfield/interface/solver/newton_numerics.cppm b/libmeanfield/interface/solver/newton_numerics.cppm new file mode 100644 index 0000000..fd724fd --- /dev/null +++ b/libmeanfield/interface/solver/newton_numerics.cppm @@ -0,0 +1,83 @@ +module; + +#include +#include +#include + +export module mean_field:solver.newton_numerics; + +export import :solver.enthalpy_admissibility; +export import :solver.linear_backend; + +export namespace mean_field::solver::nonlinear { + struct SampledEnthalpy final { }; + struct PolynomialEnthalpy final { + EnthalpyAdmissibilityOptions options{}; + }; + + struct MonolithicGravity final { }; + + // These tolerances control an auxiliary equation solve, not the physical + // model's acceptance threshold. The complete residual remains authoritative. + struct GravitySolveOptions final { + LinearSolveControl control{.relativeTolerance = 1.0e-10, .absoluteTolerance = 1.0e-14, + .maximumIterations = 300}; + linear::FGMRESOptions backend{.restartLength = 80}; + + void Validate() const { control.Validate(); backend.Validate(); } + }; + + struct GravityCurvatureCorrection final { + GravitySolveOptions gravity{}; + }; + + struct EliminatedGravity final { + GravitySolveOptions gravity{}; + }; + + template + concept EnthalpyNumerics = std::same_as || + std::same_as; + + template + concept GravityNumerics = std::same_as || + std::same_as || + std::same_as; + + // Types select algorithms; tolerances and the mesh-dependent work sizes + // remain runtime data. Adding a strategy cannot silently change a model's + // equation set or impose an additional specification such as rho_c. + template + struct NewtonNumerics final { + using EnthalpyType = Enthalpy; + using GravityType = Gravity; + [[no_unique_address]] Enthalpy enthalpy{}; + [[no_unique_address]] Gravity gravity{}; + }; + + template + NewtonNumerics(Enthalpy, Gravity) -> NewtonNumerics; + + template struct IsNewtonNumerics : std::false_type { }; + template + struct IsNewtonNumerics> : std::true_type { }; + template + concept NewtonNumericalConfiguration = IsNewtonNumerics>::value; + + struct NumericalStrategyDiagnostics final { + int polynomialChecks{0}; + int directionRepairs{0}; + double lastRepairFraction{0.0}; + double lastMeritDirectionalDerivative{0.0}; + int gravitySolves{0}; + int gravityIterations{0}; + int gravityFailures{0}; + double gravitySolveSeconds{0.0}; + double lastGravityTrueResidual{0.0}; + double lastGravityCorrectionNorm{0.0}; + double lastGravityPreparedDefect{0.0}; + double initialGravityResidual{0.0}; + double completedInitialGravityResidual{0.0}; + double lastLiftedRelativeTrueResidual{0.0}; + }; +} diff --git a/libmeanfield/interface/solver/stellar_context.cppm b/libmeanfield/interface/solver/stellar_context.cppm index 4e20623..9947b8b 100644 --- a/libmeanfield/interface/solver/stellar_context.cppm +++ b/libmeanfield/interface/solver/stellar_context.cppm @@ -32,6 +32,7 @@ export import :seed.stellar_equilibrium_projection; export import :solver.linear_backend; export import :solver.newton; export import :solver.stellar_structure; +export import :solver.stellar_newton_strategy; export namespace mean_field::solver { template @@ -523,6 +524,10 @@ export namespace mean_field::solver { std::unique_ptr physicalInverse; std::unique_ptr normalizedInverse; std::unique_ptr linearBackend; + // Only destruction is erased: strategy construction and all + // numerical calls retain their statically checked concrete types. + // Declaring this last keeps borrowed dependencies alive on teardown. + std::unique_ptr numericalStrategyStorage; private: [[nodiscard]] static ProblemType *ProblemPointer(const std::shared_ptr &candidate) noexcept { @@ -1209,9 +1214,37 @@ export namespace mean_field::solver { const auto &options = m_newton.options(); options.Validate(); detail::RequireCollectivelyIdenticalNewtonOptions(options, communicator); + detail::RequireCollectivelyIdenticalNumerics(m_newton.numerics(), communicator); state.BeginEvaluation(); + using Strategy = detail::PreparedNewtonStrategy, + typename NewtonType::NumericsType>; + static_assert(!Strategy::altersDirection || + std::same_as, + "Polynomial direction repair currently requires the normalized residual merit."); + state.numericalStrategyStorage.reset(); + auto preparedStrategy = std::make_unique(state, m_newton.numerics()); + auto &strategy = *preparedStrategy; + state.numericalStrategyStorage = std::move(preparedStrategy); StellarEquilibriumEvaluationDiagnostics diagnostics{}; + if (!strategy.CheckState(state.acceptedNormalizedState, diagnostics.numericalStrategies)) { + const auto initial = EvaluateMetric(state.acceptedNormalizedResidual); + diagnostics.initialResidualNorm = diagnostics.finalResidualNorm = initial.residualNorm; + return Failure(state, std::move(diagnostics), StellarEquilibriumFailureReason::inadmissible_state, + "The initial enthalpy polynomial could not be certified nonnegative."); + } + try { + if (!strategy.Initialize(state, diagnostics.numericalStrategies)) { + state.RestoreAccepted(); + const auto initial = EvaluateMetric(state.acceptedNormalizedResidual); + diagnostics.initialResidualNorm = diagnostics.finalResidualNorm = initial.residualNorm; + return Failure(state, std::move(diagnostics), StellarEquilibriumFailureReason::linear_solve_failure, + "Initial gravity completion failed its checked inner solve or prepared residual."); + } + } catch (...) { + state.RestoreAccepted(); + throw; + } nonlinear::MetricEvaluation acceptedMetric = EvaluateMetric(state.acceptedNormalizedResidual); diagnostics.initialResidualNorm = acceptedMetric.residualNorm; diagnostics.finalResidualNorm = acceptedMetric.residualNorm; @@ -1240,9 +1273,7 @@ export namespace mean_field::solver { state.linearRightHandSide = state.acceptedNormalizedResidual; state.linearRightHandSide *= -1.0; - LinearSolveReport linearReport = state.linearBackend->Solve( - state.linearRightHandSide, state.normalizedCorrection, options.linearSolve - ); + LinearSolveReport linearReport = strategy.Solve(state, options.linearSolve, diagnostics.numericalStrategies); diagnostics.lastLinearSolve = linearReport; diagnostics.totalLinearSolveSeconds += linearReport.solveSeconds; @@ -1275,6 +1306,13 @@ export namespace mean_field::solver { ); } + if (!strategy.CaptureDirection(state, diagnostics.numericalStrategies)) { + NotifyAfter(iteration, nonlinear::IterationDisposition::globalization_failure, false, 0.0, 0, + diagnostics.initialResidualNorm, acceptedMetric, acceptedMetric, linearReport, timings, state); + return Failure(state, std::move(diagnostics), StellarEquilibriumFailureReason::globalization_failure, + "No certified enthalpy repair with a decreasing linear merit was available."); + } + const nonlinear::MetricEvaluation previousMetric = acceptedMetric; if (options.backtracking.geometryPreflight) { const Clock::time_point geometryPreflightStart = Clock::now(); @@ -1345,9 +1383,15 @@ export namespace mean_field::solver { trialDisposition = nonlinear::LineSearchTrialDisposition::non_finite_state; rejectionSource = "candidate state"; ++diagnostics.nonFiniteLineSearchTrials; + } else if (!strategy.CheckState(state.trialNormalizedState, diagnostics.numericalStrategies)) { + rejectionReason = StellarEquilibriumFailureReason::inadmissible_state; + rejectionMessage = "The candidate enthalpy polynomial could not be certified nonnegative."; + trialDisposition = nonlinear::LineSearchTrialDisposition::inadmissible_state; + rejectionSource = "enthalpy polynomial"; + ++diagnostics.inadmissibleLineSearchTrials; } else { const Clock::time_point preparationStart = Clock::now(); - auto preparationResult = state.PrepareTrial(); + auto preparationResult = strategy.PrepareTrial(state, stepLength, diagnostics.numericalStrategies); preparationSeconds = std::chrono::duration(Clock::now() - preparationStart).count(); timings.trialPreparationSeconds += preparationSeconds; diagnostics.totalTrialPreparationSeconds += preparationSeconds; @@ -1422,6 +1466,12 @@ export namespace mean_field::solver { minimumJacobianDeterminant = preparationRejection.minimumJacobianDeterminant; } } + } else if (!strategy.TrialSolveSucceeded()) { + rejectionReason = StellarEquilibriumFailureReason::linear_solve_failure; + rejectionMessage = "The candidate gravity correction failed its checked solve or prepared residual."; + trialDisposition = nonlinear::LineSearchTrialDisposition::auxiliary_solve_failure; + rejectionSource = "gravity correction"; + minimumJacobianDeterminant = preparationResult->minimumJacobianDeterminant; } else if (!detail::AllStellarRanksSatisfy( AllFinite(state.trialNormalizedResidual), communicator )) { @@ -1453,7 +1503,15 @@ export namespace mean_field::solver { ++diagnostics.nonFiniteLineSearchTrials; } else if ( detail::AllStellarRanksSatisfy( - AcceptStep(acceptedMetric, trialMetric, stepLength, options.backtracking), + [&] { + if constexpr (Strategy::altersDirection) { + return trialMetric.merit <= acceptedMetric.merit + + options.backtracking.sufficientDecrease * stepLength * + diagnostics.numericalStrategies.lastMeritDirectionalDerivative; + } else { + return AcceptStep(acceptedMetric, trialMetric, stepLength, options.backtracking); + } + }(), communicator ) ) { diff --git a/libmeanfield/interface/solver/stellar_equilibrium_types.cppm b/libmeanfield/interface/solver/stellar_equilibrium_types.cppm index abb6482..2484dda 100644 --- a/libmeanfield/interface/solver/stellar_equilibrium_types.cppm +++ b/libmeanfield/interface/solver/stellar_equilibrium_types.cppm @@ -8,6 +8,7 @@ export module mean_field:solver.stellar_equilibrium_types; export import :deformation.safe_newton_step; export import :solver.linear_backend; +export import :solver.newton_numerics; export namespace mean_field::solver { enum class StellarEquilibriumFailureReason : std::uint8_t { @@ -60,5 +61,6 @@ export namespace mean_field::solver { double totalRollbackSeconds{0.0}; std::optional lastGeometryPreflight; std::optional lastLinearSolve; + nonlinear::NumericalStrategyDiagnostics numericalStrategies{}; }; } // namespace mean_field::solver diff --git a/libmeanfield/interface/solver/stellar_newton_strategy.cppm b/libmeanfield/interface/solver/stellar_newton_strategy.cppm new file mode 100644 index 0000000..448527b --- /dev/null +++ b/libmeanfield/interface/solver/stellar_newton_strategy.cppm @@ -0,0 +1,396 @@ +module; + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +export module mean_field:solver.stellar_newton_strategy; + +export import :solver.newton_numerics; +export import :solver.gravity_reduction; +export import :solver.stellar_enthalpy_strategy; + +namespace mean_field::solver::detail { + struct NewtonStrategyStorage { + virtual ~NewtonStrategyStorage() = default; + }; + + inline double StrategyDot(const mfem::Vector &a, const mfem::Vector &b, const MPI_Comm communicator) { + if (a.Size() != b.Size()) { throw std::invalid_argument("Strategy dot product has incompatible vectors."); } + double local = a * b; + double global = 0.0; + if (MPI_Allreduce(&local, &global, 1, MPI_DOUBLE, MPI_SUM, communicator) != MPI_SUCCESS) { + throw std::runtime_error("A numerical strategy could not reduce its vector product."); + } + return global; + } + + inline double StrategyNorm(const mfem::Vector &vector, const MPI_Comm communicator) { + return std::sqrt(StrategyDot(vector, vector, communicator)); + } + + template + void RequireCollectivelyIdenticalNumerics(const Recipe &recipe, const MPI_Comm communicator) { + using Enthalpy = typename Recipe::EnthalpyType; + using Gravity = typename Recipe::GravityType; + std::array integers{}; + std::array reals{}; + integers[0] = std::same_as; + integers[1] = std::same_as ? 1 : + std::same_as ? 2 : 0; + if constexpr (std::same_as) { + integers[2] = recipe.enthalpy.options.maximumSubdivisionDepth; + integers[3] = recipe.enthalpy.options.maximumSubcells; + } + if constexpr (!std::same_as) { + const auto &options = recipe.gravity.gravity; + integers[4] = options.control.maximumIterations; + integers[5] = options.backend.restartLength; + integers[6] = options.backend.printLevel; + reals = {options.control.relativeTolerance, options.control.absoluteTolerance}; + } + integers[7] = std::isfinite(reals[0]) && std::isfinite(reals[1]); + std::array minimumIntegers{}, maximumIntegers{}; + std::array minimumReals{}, maximumReals{}; + if (MPI_Allreduce(integers.data(), minimumIntegers.data(), 8, MPI_INT, MPI_MIN, communicator) != MPI_SUCCESS || + MPI_Allreduce(integers.data(), maximumIntegers.data(), 8, MPI_INT, MPI_MAX, communicator) != MPI_SUCCESS || + MPI_Allreduce(reals.data(), minimumReals.data(), 2, MPI_DOUBLE, MPI_MIN, communicator) != MPI_SUCCESS || + MPI_Allreduce(reals.data(), maximumReals.data(), 2, MPI_DOUBLE, MPI_MAX, communicator) != MPI_SUCCESS) { + throw std::runtime_error("Newton could not compare numerical strategies across ranks."); + } + if (minimumIntegers != maximumIntegers || minimumReals != maximumReals || + minimumIntegers[7] == 0) { + throw std::invalid_argument("Newton numerical strategies and controls must agree on every rank."); + } + if constexpr (!std::same_as) { recipe.gravity.gravity.Validate(); } + } + + template class PreparedGravityStrategy; + + template + class PreparedGravityStrategy final { + public: + PreparedGravityStrategy(State &, const nonlinear::MonolithicGravity &) { } + bool Initialize(State &, nonlinear::NumericalStrategyDiagnostics &) { return true; } + LinearSolveReport Solve(State &state, const LinearSolveControl &control) { + return state.linearBackend->Solve(state.linearRightHandSide, state.normalizedCorrection, control); + } + bool CaptureDirection(State &) { return true; } + auto PrepareTrial(State &state, double) { return state.PrepareTrial(); } + bool TrialSolveSucceeded() const noexcept { return true; } + void CopyDiagnostics(nonlinear::NumericalStrategyDiagnostics &) const { } + }; + + template + class StellarGravityWorkspace { + public: + using Operation = std::remove_cvref_t().normalizedOperator)>; + using Inverse = std::remove_cvref_t().normalizedInverse)>; + using Problem = std::remove_cvref_t().Problem())>; + using Subsystem = PreparedGravitySubsystem; + + StellarGravityWorkspace(State &state, const nonlinear::GravitySolveOptions &options) + : m_options(options), + subsystem(*state.normalizedOperator, *state.normalizedInverse, + MakeGravityBlockPartition(state.Problem().GetManifest().layout()), + state.Problem().GetCommunicator(), linear::FGMRES(options.backend)), + gravityResidual(subsystem.GetPartition().GravityResidualSize()), + gravityCorrection(subsystem.GetPartition().GravityValueSize()), + expectedResidual(subsystem.GetPartition().GravityResidualSize()), + gravityAction(subsystem.GetPartition().GravityResidualSize()), + desiredResidual(subsystem.GetPartition().GravityResidualSize()), + fullCorrection(state.Problem().StateSize()), fullAction(state.Problem().EquationSize()) { + m_options.Validate(); + } + + bool TrialSolveSucceeded() const noexcept { return trialSolveSucceeded; } + + void CopyDiagnostics(nonlinear::NumericalStrategyDiagnostics &report) const { + report.gravitySolves = subsystem.GetSolveCount(); + report.gravityIterations = subsystem.GetIterationCount(); + report.gravityFailures = subsystem.GetFailureCount(); + report.gravitySolveSeconds = subsystem.GetSolveSeconds(); + if (const auto &last = subsystem.GetLastSolveReport(); last.has_value()) { + report.lastGravityTrueResidual = last->trueResidualNorm; + } + report.lastGravityCorrectionNorm = lastCorrectionNorm; + report.lastGravityPreparedDefect = lastPreparedDefect; + report.lastLiftedRelativeTrueResidual = liftedRelativeDefect; + } + + protected: + void Refresh(State &state) { + // Preparing a trial invalidates both normalized and physical + // inverses. Refreshing only the restricted wrapper would hide this. + state.physicalInverse->Refresh(); + subsystem.Refresh(); + } + + bool SolveCompletion(State &state) { + subsystem.GetPartition().GatherGravityResiduals(state.trialNormalizedResidual, gravityResidual); + gravityResidual *= -1.0; + desiredResidual = 0.0; + gravityCorrection = 0.0; + lastSolve = subsystem.Solve(gravityResidual, gravityCorrection, m_options.control); + return lastSolve.Converged(); + } + + auto ApplyCorrection(State &state) { + const auto &partition = subsystem.GetPartition(); + partition.InjectGravityValues(gravityCorrection, fullCorrection); + state.normalizedOperator->Mult(fullCorrection, fullAction); + partition.GatherGravityResiduals(state.trialNormalizedResidual, expectedResidual); + partition.GatherGravityResiduals(fullAction, gravityAction); + expectedResidual += gravityAction; + lastCorrectionNorm = StrategyNorm(gravityCorrection, subsystem.GetCommunicator()); + partition.AddGravityValues(gravityCorrection, state.trialNormalizedState); + state.normalizedOperator->DenormalizeState(state.trialNormalizedState, state.candidatePhysicalState); + auto preparation = state.PrepareTrial(); + if (preparation) { + partition.GatherGravityResiduals(state.trialNormalizedResidual, gravityResidual); + const double scale = StrategyNorm(gravityResidual, subsystem.GetCommunicator()) + + StrategyNorm(expectedResidual, subsystem.GetCommunicator()) + + StrategyNorm(desiredResidual, subsystem.GetCommunicator()) + + StrategyNorm(gravityAction, subsystem.GetCommunicator()); + expectedResidual -= gravityResidual; + const double disagreement = StrategyNorm(expectedResidual, subsystem.GetCommunicator()); + gravityResidual -= desiredResidual; + const double targetDefect = StrategyNorm(gravityResidual, subsystem.GetCommunicator()); + const double tolerance = m_options.control.ConvergenceThreshold(lastSolve.rightHandSideNorm) + + 256.0 * std::numeric_limits::epsilon() * scale; + // G is affine in gravity for fixed retained variables. Verify + // that property against a fresh preparation, including its + // boundary rows, before trusting a completed candidate. + trialSolveSucceeded = std::isfinite(disagreement) && disagreement <= tolerance && + std::isfinite(targetDefect) && targetDefect <= tolerance; + lastPreparedDefect = targetDefect; + } + return preparation; + } + + nonlinear::GravitySolveOptions m_options; + Subsystem subsystem; + mfem::Vector gravityResidual, gravityCorrection, expectedResidual, gravityAction, desiredResidual; + mfem::Vector fullCorrection, fullAction; + LinearSolveReport lastSolve{}; + bool trialSolveSucceeded{true}; + double lastCorrectionNorm{0.0}; + double lastPreparedDefect{0.0}; + double liftedRelativeDefect{0.0}; + }; + + template + class PreparedGravityStrategy final + : public StellarGravityWorkspace { + using Base = StellarGravityWorkspace; + public: + PreparedGravityStrategy(State &state, const nonlinear::GravityCurvatureCorrection &policy) + : Base(state, policy.gravity), m_curvature(this->subsystem) { } + bool Initialize(State &, nonlinear::NumericalStrategyDiagnostics &) { return true; } + LinearSolveReport Solve(State &state, const LinearSolveControl &control) { + return state.linearBackend->Solve(state.linearRightHandSide, state.normalizedCorrection, control); + } + bool CaptureDirection(State &state) { + this->subsystem.Refresh(); + m_curvature.CaptureBase(state.acceptedNormalizedResidual, state.normalizedCorrection); + return true; + } + auto PrepareTrial(State &state, const double alpha) { + this->trialSolveSucceeded = true; + auto preparation = state.PrepareTrial(); + if (!preparation) { return preparation; } + this->Refresh(state); + m_curvature.BuildLinearPrediction(alpha, this->desiredResidual); + this->lastSolve = m_curvature.SolveCorrection(alpha, state.trialNormalizedResidual, + this->gravityCorrection, this->m_options.control); + this->trialSolveSucceeded = this->lastSolve.Converged(); + if (!this->trialSolveSucceeded) { return preparation; } + return this->ApplyCorrection(state); + } + private: + PreparedGravityCurvatureCorrection m_curvature; + }; + + template + class PreparedGravityStrategy final + : public StellarGravityWorkspace { + using Base = StellarGravityWorkspace; + using ReducedOperator = GravityReducedOperator; + using ReducedInverse = GravityReducedPreconditioner; + using Configuration = std::remove_cvref_t().linearBackend->GetConfiguration())>; + using Backend = PreparedLinearBackendType; + public: + PreparedGravityStrategy(State &state, const nonlinear::EliminatedGravity &policy) + : Base(state, policy.gravity), m_operator(this->subsystem, policy.gravity.control), + m_inverse(m_operator, *state.normalizedInverse), + m_backend(prepareLinearBackend(state.linearBackend->GetConfiguration(), m_operator, m_inverse, + state.Problem().GetCommunicator())), + m_rhs(m_operator.Height()), m_correction(m_operator.Width()) { } + + bool Initialize(State &state, nonlinear::NumericalStrategyDiagnostics &diagnostics) { + const auto &partition = this->subsystem.GetPartition(); + partition.GatherGravityResiduals(state.acceptedNormalizedResidual, this->gravityResidual); + diagnostics.initialGravityResidual = StrategyNorm(this->gravityResidual, this->subsystem.GetCommunicator()); + state.trialNormalizedState = state.acceptedNormalizedState; + state.trialNormalizedResidual = state.acceptedNormalizedResidual; + this->trialSolveSucceeded = this->SolveCompletion(state); + if (!this->trialSolveSucceeded) { return false; } + const auto preparation = this->ApplyCorrection(state); + if (!preparation || !this->trialSolveSucceeded || + !std::isfinite(StrategyNorm(state.trialNormalizedResidual, this->subsystem.GetCommunicator()))) { return false; } + partition.GatherGravityResiduals(state.trialNormalizedResidual, this->gravityResidual); + diagnostics.completedInitialGravityResidual = + StrategyNorm(this->gravityResidual, this->subsystem.GetCommunicator()); + state.PrepareTrialCommit(); + state.CommitPreparedTrial(preparation->minimumJacobianDeterminant); + this->subsystem.Refresh(); + return true; + } + + LinearSolveReport Solve(State &state, const LinearSolveControl &control) { + this->subsystem.Refresh(); + try { + auto inner = m_operator.TryCondenseRightHandSide(state.linearRightHandSide, m_rhs); + if (!inner.Converged()) { return inner; } + this->subsystem.GetPartition().GatherRetainedValues(state.normalizedCorrection, m_correction); + auto report = m_backend.Solve(m_rhs, m_correction, control); + if (!report.Converged()) { return report; } + inner = m_operator.TryLiftCorrection(m_correction, state.linearRightHandSide, state.normalizedCorrection); + if (!inner.Converged()) { return inner; } + state.normalizedOperator->Mult(state.normalizedCorrection, this->fullAction); + this->fullAction -= state.linearRightHandSide; + const auto communicator = this->subsystem.GetCommunicator(); + const double rhsNorm = StrategyNorm(state.linearRightHandSide, communicator); + const double defect = StrategyNorm(this->fullAction, communicator); + this->liftedRelativeDefect = rhsNorm > 0.0 ? defect / rhsNorm : (defect == 0.0 ? 0.0 : INFINITY); + // Inner Krylov truncation makes Schur actions approximate. + // The original equations decide whether the lifted result is + // actually an admissible inexact Newton correction. + // Keep backend recurrence/initial/true norms in their reduced + // coordinates. The separate lifted diagnostic reports the + // original-system check; these are different linear systems. + if (!std::isfinite(defect) || defect > control.ConvergenceThreshold(rhsNorm)) { + report.status = std::isfinite(defect) ? LinearSolveStatus::backend_failure : LinearSolveStatus::non_finite; + } + return report; + } catch (const GravityReductionFailure &failure) { + return failure.GetReport(); + } + } + + bool CaptureDirection(State &state) { + // A positivity repair changes retained coordinates. Lift that + // repaired tangent again so gravity follows the reduced manifold. + this->subsystem.GetPartition().GatherRetainedValues(state.normalizedCorrection, m_correction); + this->lastSolve = m_operator.TryLiftCorrection(m_correction, state.linearRightHandSide, + state.normalizedCorrection); + return this->lastSolve.Converged(); + } + + auto PrepareTrial(State &state, double) { + this->trialSolveSucceeded = true; + auto preparation = state.PrepareTrial(); + if (!preparation) { return preparation; } + this->Refresh(state); + this->trialSolveSucceeded = this->SolveCompletion(state); + if (!this->trialSolveSucceeded) { return preparation; } + return this->ApplyCorrection(state); + } + private: + ReducedOperator m_operator; + ReducedInverse m_inverse; + Backend m_backend; + mfem::Vector m_rhs, m_correction; + }; + + template class PreparedEnthalpyStrategy; + + template + class PreparedEnthalpyStrategy final { + public: + PreparedEnthalpyStrategy(State &, const nonlinear::SampledEnthalpy &) { } + bool Check(const mfem::Vector &, nonlinear::NumericalStrategyDiagnostics &) { return true; } + bool Repair(State &, nonlinear::NumericalStrategyDiagnostics &) { return true; } + }; + + template + class PreparedEnthalpyStrategy final { + public: + PreparedEnthalpyStrategy(State &state, const nonlinear::PolynomialEnthalpy &policy) + : m_prepared(state, policy.options) { } + bool Check(const mfem::Vector &state, nonlinear::NumericalStrategyDiagnostics &diagnostics) { + ++diagnostics.polynomialChecks; + return m_prepared.CheckState(state).status == PolynomialSignStatus::certified_nonnegative; + } + bool Repair(State &state, nonlinear::NumericalStrategyDiagnostics &diagnostics) { + const auto report = m_prepared.RepairDirection(state); + diagnostics.lastRepairFraction = report.chosenFraction; + diagnostics.lastMeritDirectionalDerivative = report.meritDirectionalDerivative; + if (report.chosenFraction > 0.0) { ++diagnostics.directionRepairs; } + return report.status == EnthalpyDirectionRepairStatus::feasible && + std::isfinite(report.meritDirectionalDerivative) && report.meritDirectionalDerivative < 0.0; + } + private: + PreparedStellarEnthalpyStrategy m_prepared; + }; + + template + class PreparedNewtonStrategy final : public NewtonStrategyStorage { + public: + using Enthalpy = typename Recipe::EnthalpyType; + using Gravity = typename Recipe::GravityType; + static constexpr bool altersDirection = std::same_as; + + PreparedNewtonStrategy(State &state, const Recipe &recipe) + : m_enthalpy(state, recipe.enthalpy), m_gravity(state, recipe.gravity) { } + + bool CheckState(const mfem::Vector &state, nonlinear::NumericalStrategyDiagnostics &diagnostics) { + return m_enthalpy.Check(state, diagnostics); + } + bool Initialize(State &state, nonlinear::NumericalStrategyDiagnostics &diagnostics) { + const bool ready = m_gravity.Initialize(state, diagnostics); + m_gravity.CopyDiagnostics(diagnostics); + return ready; + } + LinearSolveReport Solve(State &state, const LinearSolveControl &control, + nonlinear::NumericalStrategyDiagnostics &diagnostics) { + auto report = m_gravity.Solve(state, control); + m_gravity.CopyDiagnostics(diagnostics); + return report; + } + bool CaptureDirection(State &state, nonlinear::NumericalStrategyDiagnostics &diagnostics) { + if (!m_enthalpy.Repair(state, diagnostics)) { return false; } + if constexpr (!std::same_as || altersDirection) { + if (!m_gravity.CaptureDirection(state)) { m_gravity.CopyDiagnostics(diagnostics); return false; } + } + if constexpr (altersDirection && std::same_as) { + mfem::Vector action(state.Problem().EquationSize()); + state.normalizedOperator->Mult(state.normalizedCorrection, action); + diagnostics.lastMeritDirectionalDerivative = + StrategyDot(state.acceptedNormalizedResidual, action, state.Problem().GetCommunicator()); + if (!std::isfinite(diagnostics.lastMeritDirectionalDerivative) || + diagnostics.lastMeritDirectionalDerivative >= 0.0) { return false; } + } + m_gravity.CopyDiagnostics(diagnostics); + return true; + } + auto PrepareTrial(State &state, const double alpha, nonlinear::NumericalStrategyDiagnostics &diagnostics) { + auto preparation = m_gravity.PrepareTrial(state, alpha); + m_gravity.CopyDiagnostics(diagnostics); + return preparation; + } + bool TrialSolveSucceeded() const noexcept { return m_gravity.TrialSolveSucceeded(); } + private: + PreparedEnthalpyStrategy m_enthalpy; + PreparedGravityStrategy m_gravity; + }; +}