diff --git a/CMakeLists.txt b/CMakeLists.txt index 42f85da..a0f6f18 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -155,6 +155,7 @@ target_sources(mean_field libmeanfield/interface/solver/enthalpy_admissibility.cppm libmeanfield/interface/solver/gravity_subsystem.cppm libmeanfield/interface/solver/gravity_reduction.cppm + libmeanfield/interface/solver/coupled_trial_correction.cppm libmeanfield/interface/solver/stellar_enthalpy_strategy.cppm libmeanfield/interface/solver/newton_numerics.cppm libmeanfield/interface/solver/stellar_newton_strategy.cppm @@ -376,6 +377,7 @@ add_executable(tests tests/solver/gravity_subsystem.cpp tests/solver/gravity_compact_applications.cpp tests/solver/gravity_reduction.cpp + tests/solver/coupled_trial_correction.cpp tests/solver/stellar_numerical_strategies.cpp ) diff --git a/experiments/COUPLED_TRIAL_CORRECTION_2026-09-18.md b/experiments/COUPLED_TRIAL_CORRECTION_2026-09-18.md new file mode 100644 index 0000000..891e56f --- /dev/null +++ b/experiments/COUPLED_TRIAL_CORRECTION_2026-09-18.md @@ -0,0 +1,80 @@ +# Bounded coupled second-order trial correction + +Branch `research/coupled-trial-correction-20260918`, based on the completed +nonlinear obstruction investigation at `2944dd5`. The configured worktree +retains its existing name, `worktrees/preconditioner-improvement`, so that +build reuse does not disturb the preceding branch or immutable executables. + +## Hypothesis + +At the recorded n = 1.5, J = 0 endpoint, tenfold tighter linear accuracy did +not improve the useful trial range. The full directional derivative was +consistent with the prepared residual. Gravity-only correction reduced the +remainder norm but increased its adverse projection through hydrostatic, +EOS and surface balance. The next hypothesis is that correcting the entire +trial remainder can preserve the intended Newton response over larger steps. + +For the captured base residual F0 and verified Newton direction p, preserve +L(alpha) = F0 + alpha J0 p. At a prepared candidate y, solve + + J(y) c = -(F(y) - L(alpha)), + +then prepare and evaluate y+c. All equations, normalization factors and +ordinary full-residual acceptance thresholds remain unchanged. A separate +auxiliary true-defect check is required; the base Newton report is never +replaced by a correction report. The locally second-order argument requires +a bounded inverse, which is not assumed to hold uniformly in this model. + +## Implementation contracts + +- A third, defaulted compile-time trial-path policy composes with the + established enthalpy and gravity policies. Existing two-policy recipes + retain their types and behavior. The new policy is incompatible with the + eliminated-gravity path until a corresponding reduced-manifold contract + is designed; unsupported combinations fail at compile time. +- The auxiliary solve has its own backend, a finite iteration budget and an + independent application of the complete candidate Jacobian. Its initial + options are relative tolerance 0.03, absolute tolerance 1e-14, at most160 + iterations and restart160. +- Limit the correction norm to alpha||p|| by default. The initial work limit + is one actual auxiliary solve per Newton iteration, and one solve per + selected candidate. Runtime expert controls change these bounds explicitly. +- Select correction only after a prepared candidate fails the existing + acceptance predicate. Reuse the exact configured metric for selecting + between original and corrected candidates. An already acceptable step + incurs no coupled solve; an unprepared candidate cannot invoke one. +- Recheck the same selected admissibility policy and actual geometry after + correction. Retain the original candidate if the auxiliary solve fails, + the norm cap rejects it, preparation/admissibility fails, or its metric + worsens. Restoring a fallback includes fresh operator preparation, not + merely copying state vectors. +- The exported numerical core and trial transaction are shared by production + and frozen diagnostics. The physics-facing model API is unchanged. + +## Initial screen + +The committed `coupled_trial_screen_2026-09-18.json` records complete controls +and immutable input paths. All heavy builds and numerical cases run serially. +Evidence is written to `diagnostics/coupled_trial_2026-09-18/` in the task +workspace; individual campaign directories are immutable. + +Use the previously saved late eta0.03 direction, independently reapplying the +full Jacobian at the restored base before using it. This avoids conflating a +new linear solve with the new trial map. Test both raw and gravity-corrected +starting candidates at alpha1/16,1/64,1/256,1/1024. The middle two previously +prepared but failed the merit test. The largest failed pressure preparation; +the smallest already decreased the residual and must not trigger correction. + +The seed control uses its saved eta0.03 direction and alpha1 with the +gravity-corrected path. Its candidate already passes; it should preserve that +result without auxiliary coupled work. Each frozen alpha is an independent +trial with a fresh per-iteration budget, not a simulated live line search. + +Record inputs/factors/layout, F0,p,J0p, the actual auxiliary RHS and J(y)c, +candidate and returned state/residual vectors, true defects, norm-cap gates, +selection/fallback decisions, restoration invariants and complete work. +Useful frozen progress is a prerequisite for a bounded nonlinear campaign. +Any accepted model must also be checked for virial balance and internal +admissibility; a lower residual alone is insufficient. + +Implementation and test results are pending. diff --git a/experiments/coupled_trial_screen_2026-09-18.json b/experiments/coupled_trial_screen_2026-09-18.json new file mode 100644 index 0000000..03b915e --- /dev/null +++ b/experiments/coupled_trial_screen_2026-09-18.json @@ -0,0 +1,138 @@ +[ + { + "id": "late_curvature_coupled", + "arguments": [ + "--index", + "1.5", + "--angular-momentum", + "0", + "--no-central", + "--max-newton", + "0", + "--coupled-trial-probe", + "--strategy", + "curvature-coupled", + "--probe-state", + "/Users/tboudreaux/.codex/.chatgpt-projects/g-p-6a6b3f2194ec8191bebcd59075897f1d/diagnostics/nonlinear_step_2026-09-18/late_v1/n1p5_late_eta003/nonlinear_step_physical_base.txt", + "--probe-direction", + "/Users/tboudreaux/.codex/.chatgpt-projects/g-p-6a6b3f2194ec8191bebcd59075897f1d/diagnostics/nonlinear_step_2026-09-18/late_v1/n1p5_late_eta003/nonlinear_step_direction_full.txt", + "--linear-tol", + "0.03", + "--restart", + "160", + "--preconditioner-variant", + "gravity-schur", + "--schur-sweeps", + "1", + "--schur-relaxation", + "0.5", + "--gravity-tol", + "1e-10", + "--gravity-iterations", + "300", + "--gravity-restart", + "80", + "--coupled-tol", + "0.03", + "--coupled-iterations", + "160", + "--coupled-restart", + "160", + "--coupled-max-ratio", + "1", + "--coupled-attempts", + "1" + ] + }, + { + "id": "late_raw_coupled", + "arguments": [ + "--index", + "1.5", + "--angular-momentum", + "0", + "--no-central", + "--max-newton", + "0", + "--coupled-trial-probe", + "--strategy", + "coupled", + "--probe-state", + "/Users/tboudreaux/.codex/.chatgpt-projects/g-p-6a6b3f2194ec8191bebcd59075897f1d/diagnostics/nonlinear_step_2026-09-18/late_v1/n1p5_late_eta003/nonlinear_step_physical_base.txt", + "--probe-direction", + "/Users/tboudreaux/.codex/.chatgpt-projects/g-p-6a6b3f2194ec8191bebcd59075897f1d/diagnostics/nonlinear_step_2026-09-18/late_v1/n1p5_late_eta003/nonlinear_step_direction_full.txt", + "--linear-tol", + "0.03", + "--restart", + "160", + "--preconditioner-variant", + "gravity-schur", + "--schur-sweeps", + "1", + "--schur-relaxation", + "0.5", + "--gravity-tol", + "1e-10", + "--gravity-iterations", + "300", + "--gravity-restart", + "80", + "--coupled-tol", + "0.03", + "--coupled-iterations", + "160", + "--coupled-restart", + "160", + "--coupled-max-ratio", + "1", + "--coupled-attempts", + "1" + ] + }, + { + "id": "seed_curvature_coupled", + "arguments": [ + "--index", + "1.5", + "--angular-momentum", + "0", + "--no-central", + "--max-newton", + "0", + "--coupled-trial-probe", + "--strategy", + "curvature-coupled", + "--probe-state", + "/Users/tboudreaux/.codex/.chatgpt-projects/g-p-6a6b3f2194ec8191bebcd59075897f1d/diagnostics/nonlinear_step_2026-09-18/controls_v1/n1p5_seed_eta003/nonlinear_step_physical_base.txt", + "--probe-direction", + "/Users/tboudreaux/.codex/.chatgpt-projects/g-p-6a6b3f2194ec8191bebcd59075897f1d/diagnostics/nonlinear_step_2026-09-18/controls_v1/n1p5_seed_eta003/nonlinear_step_direction_full.txt", + "--linear-tol", + "0.03", + "--restart", + "160", + "--preconditioner-variant", + "gravity-schur", + "--schur-sweeps", + "1", + "--schur-relaxation", + "0.5", + "--gravity-tol", + "1e-10", + "--gravity-iterations", + "300", + "--gravity-restart", + "80", + "--coupled-tol", + "0.03", + "--coupled-iterations", + "160", + "--coupled-restart", + "160", + "--coupled-max-ratio", + "1", + "--coupled-attempts", + "1", + "--coupled-seed-control" + ] + } +] diff --git a/experiments/plateau_coupled_trial_diagnostics.hpp b/experiments/plateau_coupled_trial_diagnostics.hpp new file mode 100644 index 0000000..96687a8 --- /dev/null +++ b/experiments/plateau_coupled_trial_diagnostics.hpp @@ -0,0 +1,471 @@ +#pragma once + +// Included after import mean_field. The experiment owns observations only; +// the prepared production transaction owns correction and fallback decisions. +#include "plateau_nonlinear_step_diagnostics.hpp" + +#include +#include + +namespace experiment::plateau_diagnostics { +namespace coupled_trial_detail { +// Expose only the production workspace operations needed by this frozen probe. +// Keeping the affine gravity gate here would let the comparison silently drift +// away from the trial path used by the nonlinear solver. +template +struct GravityWorkspace final : mean_field::solver::detail::StellarGravityWorkspace { + using Base = mean_field::solver::detail::StellarGravityWorkspace; + using Base::Base; + using Base::ApplyCorrection; + using Base::BeginGravityAttempt; + using Base::Refresh; + using Base::desiredResidual; + using Base::fullCorrection; + using Base::gravityCorrection; + using Base::lastAffineDisagreement; + using Base::lastPreparedDefect; + using Base::lastPreparedTargetThreshold; + using Base::lastSolve; + using Base::subsystem; + using Base::trialSolveSucceeded; +}; +} + +struct CoupledTrialProbeOptions final { + double outerLinearTolerance{.03}; + mean_field::solver::nonlinear::GravitySolveOptions gravity{}; + mean_field::solver::CoupledTrialCorrectionOptions coupled{}; + bool gravityPreconditionPath{true}; + bool seedControl{false}; + double sufficientDecrease{1.e-4}; +}; + +template +void ProbeCoupledTrial(State &state, const mfem::Vector &normalizedBase, + const mfem::Vector &suppliedDirection, const std::filesystem::path &output, + CoupledTrialProbeOptions options = {}) { + using namespace mean_field; + namespace d = nonlinear_step_detail; + using GravityWorkspace = coupled_trial_detail::GravityWorkspace; + using GravityCorrection = solver::PreparedGravityCurvatureCorrection; + using Transaction = solver::detail::PreparedCoupledTrialTransaction; + const auto communicator = state.Problem().GetCommunicator(); + int ranks = 0; + MPI_Comm_size(communicator, &ranks); + const int size = state.Problem().StateSize(); + if (ranks != 1 || size != state.Problem().EquationSize() || normalizedBase.Size() != size || + suppliedDirection.Size() != size || !d::Finite(normalizedBase) || !d::Finite(suppliedDirection) || + !std::isfinite(options.outerLinearTolerance) || options.outerLinearTolerance <= 0. || + options.outerLinearTolerance >= 1. || !std::isfinite(options.sufficientDecrease) || + options.sufficientDecrease <= 0. || options.sufficientDecrease >= 1.) + throw std::invalid_argument("Coupled-trial probes require one rank, compatible finite vectors and valid controls."); + options.gravity.Validate(); + options.coupled.Validate(); + const auto start = d::Clock::now(); + const double nan = std::numeric_limits::quiet_NaN(); + const mfem::Vector base(normalizedBase), direction(suppliedDirection); + const auto &normalization = state.normalizedOperator->GetNormalization(); + const mfem::Vector factors(normalization.StateFactors()), residualFactors(normalization.ResidualFactors()); + const mfem::Vector savedAccepted(state.acceptedNormalizedState), + savedAcceptedResidual(state.acceptedNormalizedResidual), savedPhysical(state.AcceptedPhysicalState()), + savedCorrection(state.normalizedCorrection), savedPhysicalCorrection(state.physicalCorrection), + savedVolumeDirection(state.volumeDisplacementDirection); + const double savedMinimumJacobian = state.acceptedMinimumJacobianDeterminant; + std::filesystem::create_directories(output); + const auto file = [&](const std::string &name) { return output / ("coupled_trial_" + name); }; + auto summary = d::Open(file("summary.csv")); + auto preparations = d::Open(file("preparations.csv")); + auto trials = d::Open(file("trials.csv")); + auto blocks = d::Open(file("blocks.csv")); + auto manifest = d::Open(file("layout.csv")); + summary << "metric,value\n"; + preparations << "sequence,kind,sample,alpha,prepared,rejection_reason,rejection_stage,seconds,residual_norm\n"; + trials << "sample,alpha,gravity_path,raw_prepared,raw_rejection_reason,raw_rejection_stage,raw_norm,raw_merit," + "gravity_attempted,gravity_reported_gate,gravity_prepared_gate,gravity_target_defect," + "gravity_affine_disagreement,gravity_target_threshold,before_prepared,before_usable,before_norm," + "before_merit,before_accepted,soc_called,soc_status,soc_attempted,aux_status,aux_reported_gate," + "aux_iterations,aux_rhs_norm,aux_reported_true_norm,aux_verified_true_norm,aux_verified_gate," + "aux_threshold,aux_verification_applied,correction_norm,correction_cap,correction_ratio," + "corrected_metric_observed,corrected_norm,corrected_merit,corrected_accepted,selected_prepared," + "selected_norm,selected_merit,selected_accepted,selected_is_original,fallback_residual_difference," + "raw_seconds,gravity_seconds,soc_seconds,soc_refresh_seconds,soc_solve_seconds," + "soc_verification_seconds,soc_preparation_seconds,soc_preparations,soc_restorations," + "aux_operator_applications,aux_inverse_applications,verified_action_defect_norm," + "capture_seconds,total_sample_seconds,aux_core_status,rhs_saved,correction_saved,verified_vectors_saved\n"; + blocks << "sample,alpha,stage,block,offset,size,residual_norm\n"; + manifest << "kind,block,offset,size\n"; + for (const auto &value : state.Problem().GetManifest().valueBlocks()) + manifest << "value," << value.stableId << ',' << value.offset << ',' << value.size << '\n'; + for (const auto &row : state.Problem().GetManifest().residualBlocks()) + manifest << "residual," << row.stableId << ',' << row.offset << ',' << row.size << '\n'; + manifest.flush(); + const auto metric = [&](const char *name, auto value) { summary << name << ',' << value << '\n'; }; + metric("schema_version", 1); + metric("independent_frozen_alphas", 1); + metric("accepted_commits", 0); + metric("outer_linear_tolerance", options.outerLinearTolerance); + metric("sufficient_decrease", options.sufficientDecrease); + metric("gravity_path", options.gravityPreconditionPath); + metric("seed_control", options.seedControl); + metric("coupled_maximum_correction_ratio", options.coupled.maximumCorrectionRatio); + metric("coupled_relative_tolerance", options.coupled.control.relativeTolerance); + metric("coupled_absolute_tolerance", options.coupled.control.absoluteTolerance); + metric("coupled_maximum_iterations", options.coupled.control.maximumIterations); + metric("coupled_restart", options.coupled.backend.restartLength); + metric("gravity_relative_tolerance", options.gravity.control.relativeTolerance); + metric("gravity_absolute_tolerance", options.gravity.control.absoluteTolerance); + metric("gravity_maximum_iterations", options.gravity.control.maximumIterations); + metric("gravity_restart", options.gravity.backend.restartLength); + metric("polynomial_direction_repairs", 0); + metric("sampled_enthalpy_policy", 1); + metric("coupled_attempt_budget_per_independent_alpha", 1); + metric("soc_status_not_attempted", int(solver::CoupledTrialAttemptStatus::not_attempted)); + metric("soc_status_original_accepted", int(solver::CoupledTrialAttemptStatus::original_accepted)); + metric("soc_status_budget_exhausted", int(solver::CoupledTrialAttemptStatus::budget_exhausted)); + metric("soc_status_kept", int(solver::CoupledTrialAttemptStatus::kept)); + metric("soc_status_linear_failure", int(solver::CoupledTrialAttemptStatus::linear_failure)); + metric("soc_status_norm_limit", int(solver::CoupledTrialAttemptStatus::norm_limit)); + metric("soc_status_non_finite", int(solver::CoupledTrialAttemptStatus::non_finite)); + metric("soc_status_inadmissible", int(solver::CoupledTrialAttemptStatus::inadmissible)); + metric("soc_status_not_improved", int(solver::CoupledTrialAttemptStatus::not_improved)); + std::vector alphas = options.seedControl ? std::vector{1.} + : std::vector{1. / 16., 1. / 64., 1. / 256., 1. / 1024.}; + for (std::size_t sample = 0; sample < alphas.size(); ++sample) + summary << "alpha_" << sample << ',' << alphas[sample] << '\n'; + mfem::Vector residual(size), baseAction(size); + bool captured = false, haveBaseAction = false, finished = false; + int preparationsCount = 0, refreshes = 0, directActions = 0; + double preparationSeconds = 0., refreshSeconds = 0., actionSeconds = 0., restorationError = 0.; + const auto prepare = [&](const mfem::Vector &candidate, const char *kind, int sample, double alpha) { + state.trialNormalizedState = candidate; + state.normalizedOperator->DenormalizeState(candidate, state.candidatePhysicalState); + const auto started = d::Clock::now(); + auto result = state.PrepareTrial(); + const double seconds = d::Seconds(started); + preparationSeconds += seconds; + preparations << ++preparationsCount << ',' << kind << ',' << sample << ',' << alpha << ',' << bool(result) + << ',' << (result ? -1 : int(result.error().reason)) << ',' + << (result ? -1 : int(result.error().stage)) << ',' << seconds << ',' + << (result ? d::Norm(state.trialNormalizedResidual) : nan) << '\n'; + preparations.flush(); + if (result && !d::Finite(state.trialNormalizedResidual)) + throw std::runtime_error("Coupled-trial preparation produced a nonfinite residual."); + return result; + }; + const auto refresh = [&]() { + const auto started = d::Clock::now(); + state.physicalInverse->Refresh(); + refreshSeconds += d::Seconds(started); + ++refreshes; + }; + const auto restore = [&](const char *kind) { + if (!prepare(base, kind, -1, 0.)) throw std::runtime_error("The coupled-trial base cannot be restored."); + if (!captured) { residual = state.trialNormalizedResidual; captured = true; } + restorationError = std::max(restorationError, d::Difference(residual, state.trialNormalizedResidual)); + refresh(); + }; + const auto action = [&](const mfem::Vector &input, mfem::Vector &result) { + const auto started = d::Clock::now(); + state.normalizedOperator->Mult(input, result); + actionSeconds += d::Seconds(started); + ++directActions; + if (!d::Finite(result)) throw std::runtime_error("The coupled-trial full-J action is nonfinite."); + }; + const auto evaluate = [&](const mfem::Vector &value) { + return solver::nonlinear::getMetric(solver::nonlinear::NormalizedResidualMetric{}, value, communicator); + }; + const auto blockRows = [&](std::size_t sample, double alpha, const char *stage, const mfem::Vector &value) { + blocks << sample << ',' << alpha << ',' << stage << ",all,0," << size << ',' << d::Norm(value) << '\n'; + for (const auto &row : state.Problem().GetManifest().residualBlocks()) + blocks << sample << ',' << alpha << ',' << stage << ',' << row.stableId << ',' << row.offset << ',' + << row.size << ',' << d::Norm(value, row.offset, row.size) << '\n'; + blocks.flush(); + }; + const auto finish = [&]() { + // Output is append-only evidence. A failed cleanup must not retry a + // previously written filename or replace the original failure record. + finished = true; + const bool directionsNeverChanged = d::Identical(savedCorrection, state.normalizedCorrection) && + d::Identical(savedPhysicalCorrection, state.physicalCorrection) && + d::Identical(savedVolumeDirection, state.volumeDisplacementDirection); + state.normalizedCorrection = savedCorrection; + state.physicalCorrection = savedPhysicalCorrection; + state.volumeDisplacementDirection = savedVolumeDirection; + restore("final_base_restore"); + d::Save(file("restored_residual.txt"), state.trialNormalizedResidual); + mfem::Vector restoredAction(size); + action(direction, restoredAction); + d::Save(file("restored_action.txt"), restoredAction); + const bool acceptedIdentical = d::Identical(savedAccepted, state.acceptedNormalizedState) && + d::Identical(savedAcceptedResidual, state.acceptedNormalizedResidual) && + d::Identical(savedPhysical, state.AcceptedPhysicalState()) && + savedMinimumJacobian == state.acceptedMinimumJacobianDeterminant; + const bool factorsIdentical = d::Identical(factors, normalization.StateFactors()) && + d::Identical(residualFactors, normalization.ResidualFactors()); + const bool directionsIdentical = d::Identical(savedCorrection, state.normalizedCorrection) && + d::Identical(savedPhysicalCorrection, state.physicalCorrection) && + d::Identical(savedVolumeDirection, state.volumeDisplacementDirection); + metric("accepted_buffers_bitwise_identical", acceptedIdentical); + metric("normalization_factors_bitwise_identical", factorsIdentical); + metric("direction_buffers_bitwise_identical", directionsIdentical); + metric("direction_buffers_unchanged_before_restoration", directionsNeverChanged); + metric("maximum_base_restoration_error", restorationError); + metric("restored_base_action_error", haveBaseAction ? d::Difference(baseAction, restoredAction) : nan); + metric("fresh_operator_and_inverse_current", state.normalizedOperator->IsPrepared() && + state.physicalInverse->IsCurrent() && state.normalizedInverse->IsCurrent()); + metric("explicit_probe_preparations", preparationsCount); + metric("explicit_probe_preparation_seconds", preparationSeconds); + metric("explicit_probe_inverse_refreshes", refreshes); + metric("explicit_probe_inverse_refresh_seconds", refreshSeconds); + metric("explicit_probe_operator_applications", directActions); + metric("explicit_probe_operator_seconds", actionSeconds); + metric("total_seconds", d::Seconds(start)); + summary.flush(); + if (!acceptedIdentical || !factorsIdentical || !directionsIdentical || !directionsNeverChanged) + throw std::runtime_error("The coupled-trial probe changed accepted buffers, directions or normalization."); + }; + + try { + restore("initial_base_restore"); + mfem::Vector physicalBase(size); + state.normalizedOperator->DenormalizeState(base, physicalBase); + action(direction, baseAction); + haveBaseAction = true; + mfem::Vector defect(baseAction); + defect += residual; + const auto baseMetric = evaluate(residual); + const double outerThreshold = options.outerLinearTolerance * baseMetric.residualNorm; + const bool outerGate = d::Norm(defect) <= outerThreshold; + d::Save(file("base_normalized.txt"), base); + d::Save(file("base_physical.txt"), physicalBase); + d::Save(file("base_residual.txt"), residual); + d::Save(file("state_factors.txt"), factors); + d::Save(file("residual_factors.txt"), residualFactors); + d::Save(file("direction.txt"), direction); + d::Save(file("base_action.txt"), baseAction); + d::Save(file("base_defect.txt"), defect); + metric("base_residual_norm", baseMetric.residualNorm); + metric("base_merit", baseMetric.merit); + metric("direction_norm", d::Norm(direction)); + metric("base_merit_slope", d::Dot(residual, baseAction, 0, size)); + metric("outer_verified_defect_norm", d::Norm(defect)); + metric("outer_verified_relative_defect", d::Ratio(d::Norm(defect), baseMetric.residualNorm)); + metric("outer_verified_threshold", outerThreshold); + metric("outer_verified_gate", outerGate); + metric("outer_gate_failure_skips_trials", !outerGate); + summary.flush(); + if (!outerGate) { metric("probe_exception", 0); finish(); return; } + + Transaction transaction(state, options.coupled, 1); + std::optional gravity; + std::optional gravityCorrection; + if (options.gravityPreconditionPath) { + gravity.emplace(state, options.gravity); + gravityCorrection.emplace(gravity->subsystem); + metric("compact_gravity_operator", gravity->subsystem.UsesCompactOperator()); + metric("compact_gravity_inverse", gravity->subsystem.UsesCompactPreconditioner()); + } + for (std::size_t sample = 0; sample < alphas.size(); ++sample) { + const auto sampleStart = d::Clock::now(); + const double alpha = alphas[sample]; + const std::string suffix = "_a" + std::to_string(sample) + ".txt"; + restore("capture_base"); + const auto captureStart = d::Clock::now(); + // These are independent frozen trials. Resetting here does not + // represent several attempts inside one production line search. + transaction.CaptureBase(residual, direction); + transaction.BeginTrial(); + if (gravity) { + gravity->subsystem.Refresh(); + gravityCorrection->CaptureBase(residual, direction); + gravity->BeginGravityAttempt(); + } + const double captureSeconds = d::Seconds(captureStart); + mfem::Vector prediction(residual), rawState(base); + prediction.Add(alpha, baseAction); + rawState.Add(alpha, direction); + d::Save(file("prediction" + suffix), prediction); + d::Save(file("raw_state" + suffix), rawState); + const auto rawStart = d::Clock::now(); + auto preparation = prepare(rawState, "raw", static_cast(sample), alpha); + const double rawSeconds = d::Seconds(rawStart); + const bool rawPrepared = bool(preparation); + const int rawReason = rawPrepared ? -1 : int(preparation.error().reason); + const int rawStage = rawPrepared ? -1 : int(preparation.error().stage); + solver::nonlinear::MetricEvaluation rawMetric{nan, nan}; + if (rawPrepared) { + rawMetric = evaluate(state.trialNormalizedResidual); + d::Save(file("raw_residual" + suffix), state.trialNormalizedResidual); + blockRows(sample, alpha, "raw", state.trialNormalizedResidual); + } + bool gravityAttempted = false, gravityReportedGate = !gravity, gravityPreparedGate = !gravity; + double gravitySeconds = 0.; + if (rawPrepared && gravity) { + const auto gravityStart = d::Clock::now(); + gravity->Refresh(state); + gravityCorrection->BuildLinearPrediction(alpha, gravity->desiredResidual); + gravityAttempted = true; + gravity->lastSolve = gravityCorrection->SolveCorrection( + alpha, state.trialNormalizedResidual, gravity->gravityCorrection, options.gravity.control); + gravity->trialSolveSucceeded = gravity->lastSolve.Converged(); + gravityReportedGate = gravity->trialSolveSucceeded; + if (gravityReportedGate) { + // This is the same production affine-gravity assessment, + // including its contribution-aware arithmetic allowance. + preparation = gravity->ApplyCorrection(state); + gravityPreparedGate = bool(preparation) && gravity->TrialSolveSucceeded(); + d::Save(file("gravity_correction" + suffix), gravity->fullCorrection); + } + gravitySeconds = d::Seconds(gravityStart); + } + const bool beforePrepared = bool(preparation); + const bool usable = beforePrepared && gravityReportedGate && gravityPreparedGate; + mfem::Vector beforeState(state.trialNormalizedState), beforeResidual; + solver::nonlinear::MetricEvaluation beforeMetric{nan, nan}; + if (beforePrepared) { + beforeResidual = state.trialNormalizedResidual; + beforeMetric = evaluate(beforeResidual); + d::Save(file("before_state" + suffix), beforeState); + d::Save(file("before_residual" + suffix), beforeResidual); + blockRows(sample, alpha, "before", beforeResidual); + } + const auto accepts = [&](const solver::nonlinear::MetricEvaluation &candidate) { + return std::isfinite(candidate.merit) && candidate.merit >= 0. && + candidate.merit <= (1. - options.sufficientDecrease * alpha) * baseMetric.merit; + }; + const auto better = [](const auto &candidate, const auto &original) { + return std::isfinite(candidate.merit) && candidate.merit >= 0. && + std::isfinite(candidate.residualNorm) && candidate.residualNorm >= 0. && + candidate.merit < original.merit; + }; + const bool beforeAccepted = usable && accepts(beforeMetric); + auto selectedMetric = beforeMetric; + bool correctedObserved = false; + mfem::Vector observedCorrectedState, observedCorrectedResidual; + solver::nonlinear::MetricEvaluation observedCorrectedMetric{nan, nan}; + const auto observeMetric = [&](const mfem::Vector &value) { + const auto result = evaluate(value); + // This callback observes the actual prepared correction before + // a possible merit-based rollback. It neither reconstructs nor + // reapplies the remedy and does not change its selection rule. + if (!d::Identical(state.trialNormalizedState, beforeState)) { + observedCorrectedState = state.trialNormalizedState; + observedCorrectedResidual = value; + observedCorrectedMetric = result; + correctedObserved = true; + } + return result; + }; + const auto previousDiagnostics = transaction.GetDiagnostics(); + double socSeconds = 0.; + if (usable) { + const auto socStart = d::Clock::now(); + transaction.TryImprove(state, alpha, preparation, selectedMetric, observeMetric, + accepts, better, [](const mfem::Vector &) { return true; }); + socSeconds = d::Seconds(socStart); + } + const auto &diagnostics = transaction.GetDiagnostics(); + const auto &report = diagnostics.currentSolve; + const bool attempted = report.has_value(); + const bool rhsSaved = attempted && d::Finite(transaction.GetCore().GetRightHandSide()); + const bool correctionSaved = attempted && d::Finite(transaction.GetCorrection()); + const bool verifiedVectorsSaved = attempted && report->verificationApplied && rhsSaved && + d::Finite(transaction.GetCore().GetVerifiedAction()) && d::Finite(transaction.GetCore().GetVerifiedDefect()); + double savedDefectNorm = nan; + if (attempted) { + if (rhsSaved) d::Save(file("soc_rhs" + suffix), transaction.GetCore().GetRightHandSide()); + if (correctionSaved) d::Save(file("soc_correction" + suffix), transaction.GetCorrection()); + if (verifiedVectorsSaved) { + d::Save(file("soc_action" + suffix), transaction.GetCore().GetVerifiedAction()); + mfem::Vector savedDefect(transaction.GetCore().GetVerifiedAction()); + savedDefect -= transaction.GetCore().GetRightHandSide(); + savedDefectNorm = d::Norm(savedDefect); + d::Save(file("soc_defect" + suffix), savedDefect); + } + } + if (correctedObserved) { + d::Save(file("corrected_state" + suffix), observedCorrectedState); + d::Save(file("corrected_residual" + suffix), observedCorrectedResidual); + mfem::Vector physical(size); + state.normalizedOperator->DenormalizeState(observedCorrectedState, physical); + d::Save(file("corrected_physical" + suffix), physical); + blockRows(sample, alpha, "corrected", observedCorrectedResidual); + } + const bool selectedPrepared = bool(preparation); + const bool selectedIsOriginal = beforePrepared && d::Identical(beforeState, state.trialNormalizedState); + double fallbackDifference = nan; + if (selectedPrepared) { + selectedMetric = evaluate(state.trialNormalizedResidual); + d::Save(file("selected_state" + suffix), state.trialNormalizedState); + d::Save(file("selected_residual" + suffix), state.trialNormalizedResidual); + d::Save(file("selected_physical" + suffix), state.candidatePhysicalState); + blockRows(sample, alpha, "selected", state.trialNormalizedResidual); + if (selectedIsOriginal) fallbackDifference = d::Difference(beforeResidual, state.trialNormalizedResidual); + } + const bool selectedAccepted = usable && selectedPrepared && accepts(selectedMetric); + trials << sample << ',' << alpha << ',' << options.gravityPreconditionPath << ',' << rawPrepared << ',' + << rawReason << ',' << rawStage << ',' << rawMetric.residualNorm << ',' << rawMetric.merit << ',' + << gravityAttempted << ',' << gravityReportedGate << ',' << gravityPreparedGate << ',' + << (gravityAttempted ? gravity->lastPreparedDefect : nan) << ',' + << (gravityAttempted ? gravity->lastAffineDisagreement : nan) << ',' + << (gravityAttempted ? gravity->lastPreparedTargetThreshold : nan) << ',' + << beforePrepared << ',' << usable << ',' << beforeMetric.residualNorm << ',' << beforeMetric.merit + << ',' << beforeAccepted << ',' << usable << ',' << int(diagnostics.status) << ',' << attempted + << ',' << (attempted ? int(report->linearSolve.status) : -1) << ',' + << (attempted && report->linearSolve.Converged()) << ',' + << (attempted ? report->linearSolve.iterations : 0) << ',' + << (attempted ? report->linearSolve.rightHandSideNorm : nan) << ',' + << (attempted ? report->linearSolve.trueResidualNorm : nan) << ',' + << (attempted ? report->verifiedTrueResidualNorm : nan) << ',' + << (attempted && report->verifiedGate) << ',' + << (attempted ? report->convergenceThreshold : nan) << ',' + << (attempted && report->verificationApplied) << ',' + << (attempted ? report->correctionNorm : nan) << ',' + << (attempted ? report->maximumCorrectionNorm : nan) << ',' + << (attempted ? d::Ratio(report->correctionNorm, alpha * d::Norm(direction)) : nan) << ',' + << correctedObserved << ',' << observedCorrectedMetric.residualNorm << ',' + << observedCorrectedMetric.merit << ',' << (correctedObserved && accepts(observedCorrectedMetric)) << ',' + << selectedPrepared << ',' << selectedMetric.residualNorm << ',' << selectedMetric.merit << ',' + << selectedAccepted << ',' << selectedIsOriginal << ',' << fallbackDifference << ',' + << rawSeconds << ',' << gravitySeconds << ',' << socSeconds << ',' + << diagnostics.refreshSeconds - previousDiagnostics.refreshSeconds << ',' + << diagnostics.solveSeconds - previousDiagnostics.solveSeconds << ',' + << diagnostics.verificationSeconds - previousDiagnostics.verificationSeconds << ',' + << diagnostics.preparationSeconds - previousDiagnostics.preparationSeconds << ',' + << diagnostics.preparations - previousDiagnostics.preparations << ',' + << diagnostics.restorations - previousDiagnostics.restorations << ',' + << (attempted ? report->linearSolve.operatorApplications : 0) << ',' + << (attempted ? report->linearSolve.inversePreconditionerApplications : 0) << ',' + << savedDefectNorm << ',' << captureSeconds << ',' << d::Seconds(sampleStart) << ',' + << (attempted ? int(report->status) : -1) << ',' << rhsSaved << ',' << correctionSaved << ',' + << verifiedVectorsSaved << '\n'; + trials.flush(); + summary.flush(); + } + const auto &totals = transaction.GetDiagnostics(); + metric("soc_attempts", totals.attempts); + metric("soc_kept", totals.kept); + metric("soc_fallbacks", totals.fallbacks); + metric("soc_iterations", totals.iterations); + metric("soc_solve_seconds", totals.solveSeconds); + metric("soc_verification_seconds", totals.verificationSeconds); + metric("soc_preparation_seconds", totals.preparationSeconds); + metric("soc_refresh_seconds", totals.refreshSeconds); + metric("soc_preparations", totals.preparations); + metric("soc_restorations", totals.restorations); + if (gravity) { + metric("gravity_solves", gravity->subsystem.GetSolveCount()); + metric("gravity_iterations", gravity->subsystem.GetIterationCount()); + metric("gravity_failures", gravity->subsystem.GetFailureCount()); + metric("gravity_solve_seconds", gravity->subsystem.GetSolveSeconds()); + } + } catch (...) { + const auto failure = std::current_exception(); + metric("probe_exception", 1); + if (!finished) finish(); + std::rethrow_exception(failure); + } + metric("probe_exception", 0); + finish(); +} +} // namespace experiment::plateau_diagnostics diff --git a/experiments/plateau_coupled_trial_physical.hpp b/experiments/plateau_coupled_trial_physical.hpp new file mode 100644 index 0000000..2fed596 --- /dev/null +++ b/experiments/plateau_coupled_trial_physical.hpp @@ -0,0 +1,185 @@ +#pragma once + +// Include after import mean_field. Physical replay is deliberately separate +// from timed trial selection: these quadratures do not influence acceptance. +#include "plateau_physical_diagnostics.hpp" + +#include +#include +#include +#include +#include +#include +#include + +namespace experiment::plateau_diagnostics { +namespace coupled_physical_detail { +inline bool Finite(const mfem::Vector &values) { + for (int i = 0; i < values.Size(); ++i) + if (!std::isfinite(values(i))) return false; + return true; +} + +inline bool Identical(const mfem::Vector &a, const mfem::Vector &b) { + return a.Size() == b.Size() && + (a.Size() == 0 || std::memcmp(a.GetData(), b.GetData(), sizeof(mfem::real_t) * a.Size()) == 0); +} + +inline mfem::Vector Load(const std::filesystem::path &path, int size) { + std::ifstream input(path); + long long storedSize = -1; + if (!(input >> storedSize) || storedSize != size) + throw std::invalid_argument("Coupled physical snapshot has an incompatible size: " + path.string()); + mfem::Vector values(size); + for (int i = 0; i < size; ++i) + if (!(input >> values(i)) || !std::isfinite(values(i))) + throw std::invalid_argument("Coupled physical snapshot is incomplete or nonfinite: " + path.string()); + std::string trailing; + if (input >> trailing || input.bad()) + throw std::invalid_argument("Coupled physical snapshot has trailing data or a read error: " + path.string()); + return values; +} + +// PhysicalState's structural constructor uses these historical accepted-state +// member names. This adapter supplies an explicitly prepared, uncommitted +// snapshot instead; it never aliases a writable accepted vector. Reusing that +// constructor preserves the same reduced-field scatter and Piola transforms. +template +struct PreparedSnapshotView final { + State &state; + decltype(std::declval().normalizedOperator) &normalizedOperator; + const mfem::Vector &acceptedNormalizedResidual; + + explicit PreparedSnapshotView(State &source) + : state(source), normalizedOperator(source.normalizedOperator), + acceptedNormalizedResidual(source.trialNormalizedResidual) { + if (!normalizedOperator->IsPrepared() || !Finite(acceptedNormalizedResidual)) + throw std::logic_error("Physical snapshot reconstruction requires a fresh prepared trial."); + } + decltype(auto) Problem() const { return state.Problem(); } + const mfem::Vector &AcceptedPhysicalState() const { return normalizedOperator->GetPhysicalState(); } +}; +} // namespace coupled_physical_detail + +template +void MeasureCoupledTrialSnapshots(State &state, const mean_field::fem::FEM &fem, + const std::filesystem::path &output, double n, double K) { + namespace d = coupled_physical_detail; + int ranks = 0; + MPI_Comm_size(state.Problem().GetCommunicator(), &ranks); + if (ranks != 1 || !std::isfinite(n) || n < 1. || !std::isfinite(K) || K <= 0.) + throw std::invalid_argument("Coupled physical replay requires one rank, finite n>=1 and finite K>0."); + + constexpr std::string_view prefix = "coupled_trial_selected_physical_a"; + constexpr std::string_view suffix = ".txt"; + std::vector> snapshots; + for (const auto &entry : std::filesystem::directory_iterator(output)) { + const std::string name = entry.path().filename().string(); + if (!entry.is_regular_file() || !name.starts_with(prefix) || !name.ends_with(suffix)) continue; + const std::string digits = name.substr(prefix.size(), name.size() - prefix.size() - suffix.size()); + if (digits.empty() || !std::all_of(digits.begin(), digits.end(), [](char ch) { return ch >= '0' && ch <= '9'; })) + throw std::invalid_argument("Coupled physical snapshot has an invalid sample index: " + name); + snapshots.emplace_back(std::stoi(digits), entry.path()); + } + std::sort(snapshots.begin(), snapshots.end(), [](const auto &a, const auto &b) { return a.first < b.first; }); + for (std::size_t i = 1; i < snapshots.size(); ++i) + if (snapshots[i - 1].first == snapshots[i].first) + throw std::invalid_argument("Coupled physical replay found duplicate sample indices."); + + const auto path = output / "coupled_trial_physical.csv"; + if (std::filesystem::exists(path)) + throw std::invalid_argument("Coupled physical replay output already exists: " + path.string()); + std::ofstream stream(path); + stream.exceptions(std::ios::failbit | std::ios::badbit); + stream << std::setprecision(17) << "sample,source,quadrature,metric,value\n"; + const auto write = [&](int sample, const std::string &source, int quadrature, const std::string &metric, double value) { + stream << sample << ',' << source << ',' << quadrature << ',' << metric << ',' << value << '\n'; + }; + const mfem::Vector accepted(state.acceptedNormalizedState), acceptedResidual(state.acceptedNormalizedResidual), + acceptedPhysical(state.AcceptedPhysicalState()), correction(state.normalizedCorrection), + physicalCorrection(state.physicalCorrection), volumeDirection(state.volumeDisplacementDirection); + const auto &normalization = state.normalizedOperator->GetNormalization(); + const mfem::Vector stateFactors(normalization.StateFactors()), residualFactors(normalization.ResidualFactors()); + const double minimumJacobian = state.acceptedMinimumJacobianDeterminant; + const auto acceptedUnchanged = [&] { + return d::Identical(accepted, state.acceptedNormalizedState) && + d::Identical(acceptedResidual, state.acceptedNormalizedResidual) && + d::Identical(acceptedPhysical, state.AcceptedPhysicalState()) && + minimumJacobian == state.acceptedMinimumJacobianDeterminant; + }; + bool restorationAttempted = false; + const auto restore = [&] { + restorationAttempted = true; + const bool before = acceptedUnchanged(); + state.RestoreAccepted(); + const bool after = acceptedUnchanged(); + const bool factorsUnchanged = d::Identical(stateFactors, normalization.StateFactors()) && + d::Identical(residualFactors, normalization.ResidualFactors()); + const bool directionsUnchanged = d::Identical(correction, state.normalizedCorrection) && + d::Identical(physicalCorrection, state.physicalCorrection) && + d::Identical(volumeDirection, state.volumeDisplacementDirection); + const bool ready = state.normalizedOperator->IsPrepared() && state.physicalInverse->IsCurrent() && + state.normalizedInverse->IsCurrent(); + write(-1, "restoration", 0, "accepted_buffers_unchanged_before_restoration", before); + write(-1, "restoration", 0, "accepted_buffers_unchanged_after_restoration", after); + write(-1, "restoration", 0, "normalization_factors_unchanged", factorsUnchanged); + write(-1, "restoration", 0, "direction_buffers_unchanged", directionsUnchanged); + write(-1, "restoration", 0, "accepted_operator_and_inverse_current", ready); + stream.flush(); + if (!before || !after || !factorsUnchanged || !directionsUnchanged || !ready) + throw std::runtime_error("Coupled physical replay failed its accepted-state restoration invariant."); + }; + + try { + const double targetMass = state.Problem().GetStellarModel() + .template specification().targetMass().value(); + for (const auto &[sample, snapshotPath] : snapshots) { + const std::string source = snapshotPath.filename().string(); + const auto saved = d::Load(snapshotPath, state.Problem().StateSize()); + state.normalizedOperator->NormalizeState(saved, state.trialNormalizedState); + state.normalizedOperator->DenormalizeState(state.trialNormalizedState, state.candidatePhysicalState); + const auto preparation = state.PrepareTrial(); + write(sample, source, 0, "fresh_preparation_passed", bool(preparation)); + if (!preparation) { + write(sample, source, 0, "rejection_reason", int(preparation.error().reason)); + write(sample, source, 0, "rejection_stage", int(preparation.error().stage)); + throw std::runtime_error("Selected coupled trial is inadmissible during physical replay: " + source); + } + const auto &prepared = state.normalizedOperator->GetPhysicalState(); + if (prepared.Size() != saved.Size() || !d::Finite(prepared) || !d::Finite(state.trialNormalizedResidual)) + throw std::runtime_error("Selected coupled trial produced invalid prepared physical data: " + source); + long double squaredDifference = 0.; + double maximumDifference = 0., maximumScaledDifference = 0.; + for (int i = 0; i < saved.Size(); ++i) { + const double difference = std::abs(prepared(i) - saved(i)); + squaredDifference += static_cast(difference) * difference; + maximumDifference = std::max(maximumDifference, difference); + maximumScaledDifference = std::max(maximumScaledDifference, difference / std::max(1., std::abs(saved(i)))); + } + d::PreparedSnapshotView view(state); + experiment::polytrope_validation::PhysicalState physical(view, fem); + for (const int quadrature : {14, 18}) { + auto measurements = MeasurePhysical(physical, n, K, quadrature); + measurements["saved_prepared_l2_difference"] = static_cast(std::sqrt(squaredDifference)); + measurements["saved_prepared_maximum_difference"] = maximumDifference; + measurements["saved_prepared_maximum_scaled_difference"] = maximumScaledDifference; + measurements["minimum_preparation_jacobian"] = preparation->minimumJacobianDeterminant; + measurements["normalized_bordered_residual"] = physical.normalizedBorderedResidualNorm; + measurements["normalized_unbordered_residual"] = physical.normalizedUnborderedResidualNorm; + measurements["normalized_border_action"] = physical.normalizedCentralBorderActionNorm; + measurements["central_border"] = physical.centralBorder; + measurements["target_mass"] = targetMass; + measurements["mass_absolute_error"] = std::abs(measurements.at("mass") - targetMass); + for (const auto &[metric, value] : measurements) write(sample, source, quadrature, metric, value); + stream.flush(); + } + } + write(-1, "replay", 0, "selected_snapshots_measured", static_cast(snapshots.size())); + restore(); + } catch (...) { + const auto failure = std::current_exception(); + if (!restorationAttempted) restore(); + std::rethrow_exception(failure); + } +} +} // namespace experiment::plateau_diagnostics diff --git a/experiments/plateau_diagnostic.cpp b/experiments/plateau_diagnostic.cpp index c0723bf..941558c 100644 --- a/experiments/plateau_diagnostic.cpp +++ b/experiments/plateau_diagnostic.cpp @@ -24,6 +24,8 @@ import mean_field; #include "plateau_shape_curvature_diagnostics.hpp" #include "plateau_numerical_strategy_diagnostics.hpp" #include "plateau_nonlinear_step_diagnostics.hpp" +#include "plateau_coupled_trial_diagnostics.hpp" +#include "plateau_coupled_trial_physical.hpp" #include "plateau_gravity_application_diagnostics.hpp" #include "plateau_positivity_strategy_diagnostics.hpp" #include "plateau_polynomial_certificate.hpp" @@ -37,6 +39,9 @@ struct Options { bool probes = false, noCentral = false; bool strategyProbe = false, positivityProbe = false, positivityCertificate = false; bool preconditionerProbe = false, nonlinearStepProbe = false; + bool coupledTrialProbe = false, coupledSeedControl = false; + double coupledTolerance = .03, coupledMaximumRatio = 1.; + int coupledIterations = 160, coupledRestart = 160, coupledAttempts = 1; std::string preconditionerVariant = "independent"; int schurCorrectionSweeps = 1; double schurRelaxation = 0.5; @@ -103,6 +108,20 @@ Options Parse(int argc, char **argv) { o.preconditionerProbe = true; else if (a == "--nonlinear-step-probe") o.nonlinearStepProbe = true; + else if (a == "--coupled-trial-probe") + o.coupledTrialProbe = true; + else if (a == "--coupled-seed-control") + o.coupledSeedControl = true; + else if (a == "--coupled-tol") + o.coupledTolerance = std::stod(v()); + else if (a == "--coupled-iterations") + o.coupledIterations = std::stoi(v()); + else if (a == "--coupled-restart") + o.coupledRestart = std::stoi(v()); + else if (a == "--coupled-max-ratio") + o.coupledMaximumRatio = std::stod(v()); + else if (a == "--coupled-attempts") + o.coupledAttempts = std::stoi(v()); else if (a == "--preconditioner-variant") o.preconditionerVariant = v(); else if (a == "--schur-sweeps") @@ -139,17 +158,17 @@ 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.strategyProbe || o.positivityProbe) && (o.probeState.empty() || o.probeDirection.empty())) + if ((o.pressureProbe || o.shapeCurvature || o.surfaceTrace || o.strategyProbe || o.positivityProbe || o.coupledTrialProbe) && (o.probeState.empty() || o.probeDirection.empty())) throw std::invalid_argument("Focused probes require --probe-state and --probe-direction"); if (o.positivityCertificate && o.probeState.empty()) throw std::invalid_argument("--positivity-certificate requires --probe-state"); if (o.positivityCertificate && (o.probes || o.pressureProbe || o.shapeCurvature || o.surfaceTrace || - o.strategyProbe || o.positivityProbe || o.preconditionerProbe || o.nonlinearStepProbe || !o.probeDirection.empty() || !o.replayCandidate.empty())) + o.strategyProbe || o.positivityProbe || o.preconditionerProbe || o.nonlinearStepProbe || o.coupledTrialProbe || !o.probeDirection.empty() || !o.replayCandidate.empty())) throw std::invalid_argument("--positivity-certificate is a certificate-only mode; omit other probe options"); if (o.preconditionerProbe) { experiment::plateau_diagnostics::ParsePreconditionerBenchmarkVariant(o.preconditionerVariant); if (o.maxNewton != 0 || o.probes || o.pressureProbe || o.shapeCurvature || o.surfaceTrace || - o.strategyProbe || o.positivityProbe || o.positivityCertificate || o.nonlinearStepProbe || !o.probeDirection.empty() || + o.strategyProbe || o.positivityProbe || o.positivityCertificate || o.nonlinearStepProbe || o.coupledTrialProbe || !o.probeDirection.empty() || !o.replayCandidate.empty() || !o.captureIterations.empty() || o.strategy != "baseline") throw std::invalid_argument("--preconditioner-probe requires --max-newton 0 and no nonlinear/probe strategies"); if (!std::isfinite(o.linearTolerance) || o.linearTolerance <= 0.0) @@ -157,25 +176,37 @@ Options Parse(int argc, char **argv) { } if (o.nonlinearStepProbe) { if (o.maxNewton != 0 || o.probes || o.pressureProbe || o.shapeCurvature || o.surfaceTrace || - o.strategyProbe || o.positivityProbe || o.positivityCertificate || o.preconditionerProbe || + o.strategyProbe || o.positivityProbe || o.positivityCertificate || o.preconditionerProbe || o.coupledTrialProbe || !o.probeDirection.empty() || !o.replayCandidate.empty() || !o.captureIterations.empty() || o.strategy != "baseline") throw std::invalid_argument("--nonlinear-step-probe requires --max-newton 0 and no other probes or strategy; it compares raw and curvature paths itself"); if (!std::isfinite(o.linearTolerance) || o.linearTolerance <= 0.0) throw std::invalid_argument("A nonlinear-step probe requires a finite positive --linear-tol"); } + if (o.coupledTrialProbe) { + if (o.maxNewton != 0 || o.probes || o.pressureProbe || o.shapeCurvature || o.surfaceTrace || + o.strategyProbe || o.positivityProbe || o.positivityCertificate || o.preconditionerProbe || + o.nonlinearStepProbe || !o.replayCandidate.empty() || !o.captureIterations.empty() || + o.coupledAttempts != 1 || (o.strategy != "coupled" && o.strategy != "curvature-coupled")) + throw std::invalid_argument("--coupled-trial-probe requires --max-newton 0 and strategy coupled or curvature-coupled, without other probes"); + if (!std::isfinite(o.linearTolerance) || o.linearTolerance <= 0.) + throw std::invalid_argument("A coupled trial probe requires a finite positive outer --linear-tol"); + } else if (o.coupledSeedControl) { + throw std::invalid_argument("--coupled-seed-control requires --coupled-trial-probe"); + } if (!o.preconditionerProbe) { if (o.preconditionerVariant != "independent" && o.preconditionerVariant != "gravity-schur" && o.preconditionerVariant != "stellar-ldu" && o.preconditionerVariant != "stellar-ldu-gravity-schur") throw std::invalid_argument("Newton preconditioner variants are independent, gravity-schur, stellar-ldu, or stellar-ldu-gravity-schur"); if (o.preconditionerVariant != "independent") { if (o.probes || o.pressureProbe || o.shapeCurvature || o.surfaceTrace || o.strategyProbe || - o.positivityProbe || o.positivityCertificate || (!o.probeState.empty() && !o.nonlinearStepProbe) || !o.probeDirection.empty() || + o.positivityProbe || o.positivityCertificate || (!o.probeState.empty() && !o.nonlinearStepProbe && !o.coupledTrialProbe) || + (!o.probeDirection.empty() && !o.coupledTrialProbe) || !o.replayCandidate.empty()) throw std::invalid_argument("A custom Newton preconditioner cannot be combined with frozen or certificate probes"); if (o.strategy != "baseline" && o.strategy != "positivity" && o.strategy != "curvature" && - o.strategy != "positivity-curvature") - throw std::invalid_argument("Custom Newton preconditioners support baseline, positivity, curvature, or positivity-curvature"); + o.strategy != "positivity-curvature" && o.strategy != "coupled" && o.strategy != "curvature-coupled") + throw std::invalid_argument("Unsupported custom Newton preconditioner strategy"); if (!std::isfinite(o.linearTolerance) || o.linearTolerance <= 0.0 || o.schurCorrectionSweeps < 0 || !std::isfinite(o.schurRelaxation) || o.schurRelaxation <= 0.0) throw std::invalid_argument("Invalid custom Newton preconditioner controls"); @@ -184,10 +215,14 @@ Options Parse(int argc, char **argv) { 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") + o.strategy != "reduced" && o.strategy != "positivity-curvature" && o.strategy != "positivity-reduced" && + o.strategy != "coupled" && o.strategy != "curvature-coupled") 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::isfinite(o.coupledTolerance) || o.coupledTolerance <= 0. || o.coupledIterations < 1 || + o.coupledRestart < 1 || !std::isfinite(o.coupledMaximumRatio) || o.coupledMaximumRatio <= 0. || o.coupledAttempts < 1) + throw std::invalid_argument("Invalid coupled correction controls"); if (std::filesystem::exists(o.output)) throw std::runtime_error("Output directory already exists"); std::filesystem::create_directories(o.output); @@ -453,6 +488,7 @@ int RunContext(Context &context, const Options &o, const Scales &scales) { << "\nM=1\nR_reference=1\nxi1=" << scales.xi << "\nmu1=" << scales.mu << "\ncentral_phase=" << !o.noCentral << "\nrestart=" << o.restart << "\nprobe_restart=" << o.probeRestart << "\nfrozen_probe_backend=" << (o.preconditionerProbe ? "direct_mfem_fgmres_preconditioner_benchmark" : + o.coupledTrialProbe ? "supplied_outer_direction_verified_with_full_J_and_auxiliary_fgmres" : o.nonlinearStepProbe ? "prepared_fgmres_cold_with_independent_true_gate" : o.positivityCertificate ? "polynomial_certificate_no_linear_solve" : o.probeDirection.empty() ? "direct_mfem_fgmres_cold" : "supplied_direction_no_linear_solve") @@ -462,6 +498,13 @@ int RunContext(Context &context, const Options &o, const Scales &scales) { << "\ngravity_iterations=" << o.gravityIterations << "\ngravity_restart=" << o.gravityRestart << "\npreconditioner_probe=" << o.preconditionerProbe << "\nnonlinear_step_probe=" << o.nonlinearStepProbe + << "\ncoupled_trial_probe=" << o.coupledTrialProbe + << "\ncoupled_seed_control=" << o.coupledSeedControl + << "\ncoupled_tolerance=" << o.coupledTolerance + << "\ncoupled_iterations=" << o.coupledIterations + << "\ncoupled_restart=" << o.coupledRestart + << "\ncoupled_maximum_ratio=" << o.coupledMaximumRatio + << "\ncoupled_attempts_per_iteration=" << o.coupledAttempts << "\npreconditioner_variant=" << o.preconditionerVariant << "\ncontext_preconditioner_variant=" << (o.preconditionerProbe ? "independent" : o.preconditionerVariant) << "\ncontext_gravity_schur_sweeps=" @@ -470,6 +513,32 @@ int RunContext(Context &context, const Options &o, const Scales &scales) { << "\nschur_correction_sweeps=" << o.schurCorrectionSweeps << "\nschur_relaxation=" << o.schurRelaxation << "\ncompiler=" << __VERSION__ << "\ncompiled=" << __DATE__ << ' ' << __TIME__ << '\n'; metadata.flush(); + if (o.coupledTrialProbe) { + metadata << "coupled_trial_input=physical_snapshot_and_supplied_normalized_direction\n" + << "coupled_trial_admissibility=sampled_enthalpy_without_repair\n" + << "coupled_trial_alphas_are_independent=1\n" + << "probe_state=" << o.probeState << "\nprobe_direction=" << o.probeDirection << '\n'; + metadata.flush(); + solver::detail::StellarEquilibriumContextDiagnostics::WithState(context, [&](auto &s, const fem::FEM &f) { + const auto physical = LoadVector(o.probeState, s.AcceptedPhysicalState().Size()); + mfem::Vector base; + s.normalizedOperator->NormalizeState(physical, base); + const auto direction = LoadVector(o.probeDirection, base.Size()); + experiment::plateau_diagnostics::ProbeCoupledTrial(s, base, direction, o.output, { + .outerLinearTolerance = o.linearTolerance, + .gravity = {.control = {.relativeTolerance = o.gravityTolerance, .absoluteTolerance = 1e-14, + .maximumIterations = o.gravityIterations}, + .backend = {.restartLength = o.gravityRestart, .printLevel = -1}}, + .coupled = {.control = {.relativeTolerance = o.coupledTolerance, .absoluteTolerance = 1e-14, + .maximumIterations = o.coupledIterations}, + .backend = {.restartLength = o.coupledRestart, .printLevel = -1}, + .maximumCorrectionRatio = o.coupledMaximumRatio}, + .gravityPreconditionPath = o.strategy == "curvature-coupled", + .seedControl = o.coupledSeedControl}); + experiment::plateau_diagnostics::MeasureCoupledTrialSnapshots(s, f, o.output, o.n, scales.K); + }); + return 0; + } if (o.nonlinearStepProbe) { metadata << "nonlinear_step_input=" << (o.probeState.empty() ? "lane_emden_seed" : "physical_snapshot") << "\nnonlinear_step_paths=raw,curvature\nnonlinear_step_positivity_repair=0\n"; @@ -600,6 +669,12 @@ int RunContext(Context &context, const Options &o, const Scales &scales) { blocklog << "iteration,stage,block,normalized_norm,physical_norm,max_abs_normalized\n"; auto triallog = File(o.output / "trials.csv"); triallog << "iteration,trial,alpha,disposition,source,F,min_jacobian\n"; + auto coupledLog = File(o.output / "coupled_trials.csv"); + coupledLog << "iteration,trial,alpha,diagnostics_available,status,attempts_total,kept_total,fallbacks_total," + "iterations_total,preparations_total,restorations_total,solve_seconds_total,verification_seconds_total," + "preparation_seconds_total,refresh_seconds_total,original_merit,corrected_merit,has_solve," + "linear_status,linear_iterations,rhs_norm,reported_true_norm,verified_true_norm,threshold,verified_gate," + "correction_norm,maximum_correction_norm\n"; auto stateConsistency = File(o.output / "trial_state_consistency.csv"); stateConsistency << "iteration,trial,normalized_physical_roundtrip_error,relative_error\n"; mfem::Vector captureBase; @@ -643,6 +718,23 @@ int RunContext(Context &context, const Options &o, const Scales &scales) { std::cout << "Newton " << e.iteration << " F=" << e.residualNorm << std::endl; }, [&](const solver::nonlinear::AfterLineSearchTrial &e) { + coupledLog << e.iteration << ',' << e.trial << ',' << e.stepLength << ',' << (e.numericalStrategies != nullptr); + if (e.numericalStrategies) { + const auto &c = e.numericalStrategies->coupledTrial; + coupledLog << ',' << static_cast(c.status) << ',' << c.attempts << ',' << c.kept << ',' + << c.fallbacks << ',' << c.iterations << ',' << c.preparations << ',' << c.restorations << ',' + << c.solveSeconds << ',' << c.verificationSeconds << ',' << c.preparationSeconds << ',' + << c.refreshSeconds << ',' << c.originalMerit << ',' << c.correctedMerit << ',' << c.currentSolve.has_value(); + if (c.currentSolve) { + const auto &r = *c.currentSolve; + coupledLog << ',' << static_cast(r.linearSolve.status) << ',' << r.linearSolve.iterations << ',' + << r.linearSolve.rightHandSideNorm << ',' << r.linearSolve.trueResidualNorm << ',' + << r.verifiedTrueResidualNorm << ',' << r.convergenceThreshold << ',' << r.verifiedGate << ',' + << r.correctionNorm << ',' << r.maximumCorrectionNorm; + } else coupledLog << ",-1,0,nan,nan,nan,nan,0,nan,nan"; + } else coupledLog << ",-1,0,0,0,0,0,0,0,0,0,0,nan,nan,0,-1,0,nan,nan,nan,nan,0,nan,nan"; + coupledLog << '\n'; + coupledLog.flush(); if (e.candidatePhysicalState.size() != e.candidateNormalizedState.size() || e.candidateNormalizedState.size() != static_cast(stateFactors.Size())) throw std::runtime_error("Trial observer exposed incompatible state representations"); @@ -731,6 +823,12 @@ int RunContext(Context &context, const Options &o, const Scales &scales) { .control = {.relativeTolerance = o.gravityTolerance, .absoluteTolerance = 1e-14, .maximumIterations = o.gravityIterations}, .backend = {.restartLength = o.gravityRestart, .printLevel = -1}}; + const nl::CoupledCurvatureCorrection coupled{ + .correction = {.control = {.relativeTolerance = o.coupledTolerance, .absoluteTolerance = 1e-14, + .maximumIterations = o.coupledIterations}, + .backend = {.restartLength = o.coupledRestart, .printLevel = -1}, + .maximumCorrectionRatio = o.coupledMaximumRatio}, + .maximumAttemptsPerIteration = o.coupledAttempts}; auto run = [&](auto numerics) { auto solve = solver::make(context, nl::Newton(options, nl::NormalizedResidualMetric{}, numerics), observer); return solve.evaluate(); @@ -738,6 +836,8 @@ int RunContext(Context &context, const Options &o, const Scales &scales) { 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 == "coupled") return run(nl::NewtonNumerics{nl::SampledEnthalpy{}, nl::MonolithicGravity{}, coupled}); + if (o.strategy == "curvature-coupled") return run(nl::NewtonNumerics{nl::SampledEnthalpy{}, nl::GravityCurvatureCorrection{gravity}, coupled}); if constexpr (FullDefaultStack) 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}}); @@ -770,6 +870,17 @@ int RunContext(Context &context, const Options &o, const Scales &scales) { metric("initial_gravity_residual", strategyDiagnostics.initialGravityResidual); metric("completed_initial_gravity_residual", strategyDiagnostics.completedInitialGravityResidual); metric("last_lifted_relative_true_residual", strategyDiagnostics.lastLiftedRelativeTrueResidual); + const auto &coupledDiagnostics = strategyDiagnostics.coupledTrial; + metric("coupled_attempts", coupledDiagnostics.attempts); + metric("coupled_kept", coupledDiagnostics.kept); + metric("coupled_fallbacks", coupledDiagnostics.fallbacks); + metric("coupled_iterations", coupledDiagnostics.iterations); + metric("coupled_preparations", coupledDiagnostics.preparations); + metric("coupled_restorations", coupledDiagnostics.restorations); + metric("coupled_solve_seconds", coupledDiagnostics.solveSeconds); + metric("coupled_verification_seconds", coupledDiagnostics.verificationSeconds); + metric("coupled_preparation_seconds", coupledDiagnostics.preparationSeconds); + metric("coupled_refresh_seconds", coupledDiagnostics.refreshSeconds); 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 2b02681..f6e3096 100644 --- a/libmeanfield/interface/mean_field.cppm +++ b/libmeanfield/interface/mean_field.cppm @@ -107,6 +107,7 @@ export import :operators.stellar_equilibrium_system; export import :solver.enthalpy_admissibility; export import :solver.gravity_subsystem; export import :solver.gravity_reduction; +export import :solver.coupled_trial_correction; export import :solver.stellar_enthalpy_strategy; export import :solver.newton_numerics; diff --git a/libmeanfield/interface/solver/coupled_trial_correction.cppm b/libmeanfield/interface/solver/coupled_trial_correction.cppm new file mode 100644 index 0000000..cd85809 --- /dev/null +++ b/libmeanfield/interface/solver/coupled_trial_correction.cppm @@ -0,0 +1,457 @@ +module; + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +export module mean_field:solver.coupled_trial_correction; + +export import :solver.linear_backend; + +export namespace mean_field::solver { + struct CoupledTrialCorrectionOptions final { + LinearSolveControl control{.relativeTolerance = .03, .absoluteTolerance = 1.e-14, .maximumIterations = 160}; + linear::FGMRESOptions backend{.restartLength = 160, .printLevel = -1}; + double maximumCorrectionRatio{1.0}; + + void Validate() const { + control.Validate(); + backend.Validate(); + if (control.relativeTolerance <= 0.0 || !std::isfinite(maximumCorrectionRatio) || + maximumCorrectionRatio <= 0.0) + throw std::invalid_argument("Coupled trial correction requires positive relative tolerance and norm cap."); + } + }; + + enum class CoupledTrialCorrectionStatus : std::uint8_t { admissible, linear_failure, norm_limit, non_finite }; + + struct CoupledTrialCorrectionReport final { + CoupledTrialCorrectionStatus status{CoupledTrialCorrectionStatus::linear_failure}; + LinearSolveReport linearSolve{}; + double verifiedTrueResidualNorm{std::numeric_limits::quiet_NaN()}; + double convergenceThreshold{std::numeric_limits::quiet_NaN()}; + double correctionNorm{std::numeric_limits::quiet_NaN()}; + double maximumCorrectionNorm{std::numeric_limits::quiet_NaN()}; + double verificationSeconds{0.0}; + bool verifiedGate{false}; + bool verificationApplied{false}; + [[nodiscard]] bool Admissible() const noexcept { return status == CoupledTrialCorrectionStatus::admissible; } + }; + + enum class CoupledTrialAttemptStatus : std::uint8_t { + not_attempted, original_accepted, budget_exhausted, kept, linear_failure, norm_limit, + non_finite, inadmissible, not_improved + }; + + struct CoupledTrialDiagnostics final { + std::uint64_t attempts{0}, kept{0}, fallbacks{0}, iterations{0}, preparations{0}, restorations{0}; + double solveSeconds{0.0}, verificationSeconds{0.0}, preparationSeconds{0.0}, refreshSeconds{0.0}; + // These fields describe the current trial; BeginTrial clears them even + // when an earlier preparation rejection prevents another auxiliary solve. + CoupledTrialAttemptStatus status{CoupledTrialAttemptStatus::not_attempted}; + std::optional currentSolve{}; + double originalMerit{std::numeric_limits::quiet_NaN()}; + double correctedMerit{std::numeric_limits::quiet_NaN()}; + }; +} + +namespace mean_field::solver::detail { + template + concept CoupledGenerationAvailable = requires(const Operation &operation) { + { operation.GetPreparationGeneration() } -> std::convertible_to; + } || requires(const Operation &operation) { + { operation.GetProblem().GetPreparationGeneration() } -> std::convertible_to; + }; + + template + std::uint64_t CoupledGeneration(const Operation &operation) { + if constexpr (requires { operation.GetPreparationGeneration(); }) return operation.GetPreparationGeneration(); + else return operation.GetProblem().GetPreparationGeneration(); + } + + inline bool CoupledFinite(const mfem::Vector &values) { + for (int index = 0; index < values.Size(); ++index) + if (!std::isfinite(values(index))) return false; + return true; + } + inline bool CoupledAll(bool condition, MPI_Comm communicator) { + int local = condition, global = 0; + if (MPI_Allreduce(&local, &global, 1, MPI_INT, MPI_MIN, communicator) != MPI_SUCCESS) + throw std::runtime_error("A coupled trial correction could not reduce its validity check."); + return global != 0; + } + inline double CoupledNorm(const mfem::Vector &values, MPI_Comm communicator) { + const double local = values * values; + double global = 0.0; + if (MPI_Allreduce(&local, &global, 1, MPI_DOUBLE, MPI_SUM, communicator) != MPI_SUCCESS) + throw std::runtime_error("A coupled trial correction could not reduce its norm."); + return std::sqrt(global); + } + inline void CoupledRequire(bool condition, MPI_Comm communicator, const char *message) { + if (!CoupledAll(condition, communicator)) throw std::invalid_argument(message); + } + inline CoupledTrialCorrectionOptions PrepareCoupledOptions(CoupledTrialCorrectionOptions options, MPI_Comm communicator) { + int initialized = 0, finalized = 0; + if (MPI_Initialized(&initialized) != MPI_SUCCESS || !initialized || + MPI_Finalized(&finalized) != MPI_SUCCESS || finalized) + throw std::logic_error("A coupled trial correction requires initialized MPI that has not been finalized."); + if (communicator == MPI_COMM_NULL) throw std::invalid_argument("A coupled trial correction requires a communicator."); + bool valid = true; + try { options.Validate(); } catch (const std::invalid_argument &) { valid = false; } + CoupledRequire(valid, communicator, "Coupled trial correction controls are invalid on at least one rank."); + const std::array integers{options.control.maximumIterations, options.backend.restartLength, options.backend.printLevel}; + const std::array reals{options.control.relativeTolerance, options.control.absoluteTolerance, options.maximumCorrectionRatio}; + std::array imin{}, imax{}; + std::array rmin{}, rmax{}; + if (MPI_Allreduce(integers.data(), imin.data(), 3, MPI_INT, MPI_MIN, communicator) != MPI_SUCCESS || + MPI_Allreduce(integers.data(), imax.data(), 3, MPI_INT, MPI_MAX, communicator) != MPI_SUCCESS || + MPI_Allreduce(reals.data(), rmin.data(), 3, MPI_DOUBLE, MPI_MIN, communicator) != MPI_SUCCESS || + MPI_Allreduce(reals.data(), rmax.data(), 3, MPI_DOUBLE, MPI_MAX, communicator) != MPI_SUCCESS) + throw std::runtime_error("Coupled trial correction controls could not be compared across ranks."); + if (imin != imax || rmin != rmax) + throw std::invalid_argument("Coupled trial correction controls must agree across ranks."); + return options; + } + template + CoupledTrialCorrectionOptions PrepareCoupledSystemOptions(CoupledTrialCorrectionOptions options, + const Operation &operation, const Inverse &inverse, MPI_Comm communicator) { + options = PrepareCoupledOptions(options, communicator); + CoupledRequire(operation.Height() >= 0 && operation.Width() == operation.Height() && + inverse.Height() == operation.Width() && inverse.Width() == operation.Height(), communicator, + "A coupled correction requires compatible square operator and inverse dimensions."); + const long long localSize = operation.Height(); + long long globalSize = 0; + if (MPI_Allreduce(&localSize, &globalSize, 1, MPI_LONG_LONG, MPI_SUM, communicator) != MPI_SUCCESS) + throw std::runtime_error("A coupled correction could not determine its global size."); + if (globalSize == 0) throw std::invalid_argument("A coupled correction requires a globally nonempty system."); + return options; + } + inline linear::FGMRES PrepareCoupledBackend(const CoupledTrialCorrectionOptions &options, MPI_Comm communicator) { + return linear::FGMRES(PrepareCoupledOptions(options, communicator).backend); + } +} + +export namespace mean_field::solver { + template + concept CoupledTrialOperation = std::derived_from && + detail::CoupledGenerationAvailable && requires(const Operation &operation) { + { operation.IsPrepared() } -> std::same_as; + }; + + template + concept CoupledTrialInverse = std::derived_from && requires(const Inverse &inverse) { + { inverse.IsCurrent() } -> std::same_as; + }; + + // Capture fixes the affine target at the accepted linearization. Later + // SolveCorrection borrows the freshly prepared candidate Jacobian; neither + // operation owns or mutates a model state or the outer Newton direction. + template + requires LinearBackendRuntimeAvailableFor + class PreparedCoupledTrialCorrection final { + public: + PreparedCoupledTrialCorrection(const Operation &operation, Inverse &inverse, MPI_Comm communicator, + CoupledTrialCorrectionOptions options = {}) + requires std::same_as + : PreparedCoupledTrialCorrection(operation, inverse, communicator, options, + detail::PrepareCoupledBackend(options, communicator)) { } + + PreparedCoupledTrialCorrection(const Operation &operation, Inverse &inverse, MPI_Comm communicator, + CoupledTrialCorrectionOptions options, Backend backend) + : m_operation(std::addressof(operation)), m_inverse(std::addressof(inverse)), m_communicator(communicator), + m_options(detail::PrepareCoupledSystemOptions(options, operation, inverse, communicator)), + m_backend(prepareLinearBackend(std::move(backend), operation, inverse, communicator)), + m_baseResidual(operation.Height()), m_baseAction(operation.Height()), m_rightHandSide(operation.Height()), + m_verifiedAction(operation.Height()), m_defect(operation.Height()) { + detail::CoupledRequire(IsCurrent(), communicator, "Coupled correction preparation requires a current operator and inverse."); + } + PreparedCoupledTrialCorrection(const PreparedCoupledTrialCorrection &) = delete; + PreparedCoupledTrialCorrection &operator=(const PreparedCoupledTrialCorrection &) = delete; + PreparedCoupledTrialCorrection(PreparedCoupledTrialCorrection &&) = delete; + PreparedCoupledTrialCorrection &operator=(PreparedCoupledTrialCorrection &&) = delete; + + [[nodiscard]] bool IsCurrent() const { + return m_operation->IsPrepared() && m_inverse->IsCurrent() && m_backend.IsReady() && + m_operation->Height() == m_operation->Width() && + m_inverse->Height() == m_operation->Width() && m_inverse->Width() == m_operation->Height() && + m_baseResidual.Size() == m_operation->Height() && m_baseAction.Size() == m_operation->Height(); + } + [[nodiscard]] bool HasBase() const noexcept { return m_hasBase; } + [[nodiscard]] std::uint64_t GetBaseGeneration() const noexcept { return m_baseGeneration; } + [[nodiscard]] double GetBaseDirectionNorm() const noexcept { return m_directionNorm; } + [[nodiscard]] MPI_Comm GetCommunicator() const noexcept { return m_communicator; } + [[nodiscard]] const mfem::Vector &GetRightHandSide() const noexcept { return m_rightHandSide; } + [[nodiscard]] const mfem::Vector &GetVerifiedAction() const noexcept { return m_verifiedAction; } + [[nodiscard]] const mfem::Vector &GetVerifiedDefect() const noexcept { return m_defect; } + [[nodiscard]] const auto &GetLastReport() const noexcept { return m_lastReport; } + [[nodiscard]] const CoupledTrialCorrectionOptions &GetOptions() const noexcept { return m_options; } + void Invalidate() noexcept { m_hasBase = false; } + + void CaptureBase(const mfem::Vector &residual, const mfem::Vector &direction) { + m_hasBase = false; + detail::CoupledRequire(IsCurrent() && residual.Size() == m_operation->Height() && + direction.Size() == m_operation->Width() && detail::CoupledFinite(residual) && detail::CoupledFinite(direction), + m_communicator, "Coupled correction capture requires a current base and finite compatible vectors."); + m_baseResidual = residual; + m_operation->Mult(direction, m_baseAction); + m_directionNorm = detail::CoupledNorm(direction, m_communicator); + detail::CoupledRequire(detail::CoupledFinite(m_baseAction) && std::isfinite(m_directionNorm), m_communicator, + "The captured coupled correction tangent is nonfinite."); + m_baseGeneration = detail::CoupledGeneration(*m_operation); + m_hasBase = true; + } + + void BuildLinearPrediction(double alpha, mfem::Vector &prediction) const { + detail::CoupledRequire(m_hasBase && prediction.Size() == m_baseResidual.Size(), m_communicator, + "A coupled prediction requires a captured base and compatible output."); + RequireAlpha(alpha); + prediction = m_baseResidual; + prediction.Add(alpha, m_baseAction); + detail::CoupledRequire(detail::CoupledFinite(prediction), m_communicator, "The coupled affine target is nonfinite."); + } + + [[nodiscard]] CoupledTrialCorrectionReport SolveCorrection(double alpha, const mfem::Vector &trialResidual, + mfem::Vector &correction) { + detail::CoupledRequire(m_hasBase && IsCurrent() && detail::CoupledGeneration(*m_operation) >= m_baseGeneration && + trialResidual.Size() == m_operation->Height() && correction.Size() == m_operation->Width(), + m_communicator, "Coupled correction requires a captured base, current candidate and compatible vectors."); + RequireAlpha(alpha); + CoupledTrialCorrectionReport report; + const auto finish = [&]() { m_lastReport = report; return report; }; + correction = 0.0; + m_verifiedAction = std::numeric_limits::quiet_NaN(); + m_defect = std::numeric_limits::quiet_NaN(); + report.maximumCorrectionNorm = m_options.maximumCorrectionRatio * alpha * m_directionNorm; + BuildLinearPrediction(alpha, m_rightHandSide); + m_rightHandSide -= trialResidual; + if (!detail::CoupledAll(detail::CoupledFinite(m_rightHandSide) && std::isfinite(report.maximumCorrectionNorm), m_communicator)) { + report.status = CoupledTrialCorrectionStatus::non_finite; + return finish(); + } + const double rhsNorm = detail::CoupledNorm(m_rightHandSide, m_communicator); + const double relativeThreshold = m_options.control.relativeTolerance * rhsNorm; + if (!std::isfinite(rhsNorm) || !std::isfinite(relativeThreshold)) { + report.status = CoupledTrialCorrectionStatus::non_finite; + return finish(); + } + report.convergenceThreshold = m_options.control.ConvergenceThreshold(rhsNorm); + report.linearSolve = m_backend.Solve(m_rightHandSide, correction, m_options.control); + if (!detail::CoupledAll(correction.Size() == m_operation->Width() && detail::CoupledFinite(correction), m_communicator)) { + report.status = CoupledTrialCorrectionStatus::non_finite; + return finish(); + } + const auto start = std::chrono::steady_clock::now(); + m_operation->Mult(correction, m_verifiedAction); + report.verificationApplied = true; + m_defect = m_verifiedAction; + m_defect -= m_rightHandSide; + const bool finite = detail::CoupledAll(detail::CoupledFinite(m_defect), m_communicator); + report.verifiedTrueResidualNorm = finite ? detail::CoupledNorm(m_defect, m_communicator) : + std::numeric_limits::quiet_NaN(); + report.correctionNorm = detail::CoupledNorm(correction, m_communicator); + report.verificationSeconds = std::chrono::duration(std::chrono::steady_clock::now() - start).count(); + report.verifiedGate = std::isfinite(report.convergenceThreshold) && report.convergenceThreshold >= 0.0 && + std::isfinite(report.verifiedTrueResidualNorm) && + report.verifiedTrueResidualNorm <= report.convergenceThreshold; + const bool solved = detail::CoupledAll(report.linearSolve.Converged(), m_communicator); + if (!finite || !std::isfinite(report.verifiedTrueResidualNorm) || !std::isfinite(report.correctionNorm)) + report.status = CoupledTrialCorrectionStatus::non_finite; + else if (!solved || !report.verifiedGate) report.status = CoupledTrialCorrectionStatus::linear_failure; + else if (report.correctionNorm > report.maximumCorrectionNorm) report.status = CoupledTrialCorrectionStatus::norm_limit; + else report.status = CoupledTrialCorrectionStatus::admissible; + return finish(); + } + + private: + void RequireAlpha(double alpha) const { + detail::CoupledRequire(std::isfinite(alpha) && alpha >= 0.0, m_communicator, + "Coupled correction requires a finite nonnegative step length."); + double minimum = 0.0, maximum = 0.0; + if (MPI_Allreduce(&alpha, &minimum, 1, MPI_DOUBLE, MPI_MIN, m_communicator) != MPI_SUCCESS || + MPI_Allreduce(&alpha, &maximum, 1, MPI_DOUBLE, MPI_MAX, m_communicator) != MPI_SUCCESS) + throw std::runtime_error("A coupled correction could not compare step lengths."); + if (minimum != maximum) throw std::invalid_argument("Coupled correction step lengths must agree across ranks."); + } + const Operation *m_operation; + Inverse *m_inverse; + MPI_Comm m_communicator; + CoupledTrialCorrectionOptions m_options; + PreparedLinearBackendType m_backend; + mfem::Vector m_baseResidual, m_baseAction, m_rightHandSide, m_verifiedAction, m_defect; + double m_directionNorm{0.0}; + std::uint64_t m_baseGeneration{0}; + bool m_hasBase{false}; + std::optional m_lastReport{}; + }; +} + +// The state protocol is deliberately structural: production and focused probes +// execute the same transaction while their typed callbacks choose the metric +// and enthalpy policy. This helper never reads or writes accepted-state buffers. +export namespace mean_field::solver::detail { + template class PreparedCoupledTrialTransaction final { + using Operation = std::remove_cvref_t().normalizedOperator)>; + using Inverse = std::remove_cvref_t().normalizedInverse)>; + using Clock = std::chrono::steady_clock; + public: + PreparedCoupledTrialTransaction(State &state, const CoupledTrialCorrectionOptions &options, + int maximumAttemptsPerIteration = 1) + : m_core(*state.normalizedOperator, *state.normalizedInverse, state.Problem().GetCommunicator(), options), + m_maximumAttempts(maximumAttemptsPerIteration), m_originalState(state.normalizedOperator->Width()), + m_originalResidual(state.normalizedOperator->Height()), m_originalPhysical(state.normalizedOperator->Width()), + m_correction(state.normalizedOperator->Width()) { + CoupledRequire(maximumAttemptsPerIteration > 0, m_core.GetCommunicator(), "A coupled trial attempt budget must be positive."); + int minimum = 0, maximum = 0; + if (MPI_Allreduce(&maximumAttemptsPerIteration, &minimum, 1, MPI_INT, MPI_MIN, m_core.GetCommunicator()) != MPI_SUCCESS || + MPI_Allreduce(&maximumAttemptsPerIteration, &maximum, 1, MPI_INT, MPI_MAX, m_core.GetCommunicator()) != MPI_SUCCESS) + throw std::runtime_error("Coupled trial attempt budgets could not be compared."); + if (minimum != maximum) throw std::invalid_argument("Coupled trial attempt budgets must agree across ranks."); + } + void CaptureBase(const mfem::Vector &residual, const mfem::Vector &direction) { + m_core.CaptureBase(residual, direction); + m_attemptsThisIteration = 0; + } + void BeginTrial() noexcept { + m_diagnostics.status = CoupledTrialAttemptStatus::not_attempted; + m_diagnostics.currentSolve.reset(); + m_diagnostics.originalMerit = m_diagnostics.correctedMerit = std::numeric_limits::quiet_NaN(); + } + [[nodiscard]] const CoupledTrialDiagnostics &GetDiagnostics() const noexcept { return m_diagnostics; } + [[nodiscard]] const auto &GetCore() const noexcept { return m_core; } + [[nodiscard]] const mfem::Vector &GetCorrection() const noexcept { return m_correction; } + + template + void TryImprove(State &state, double alpha, Preparation &preparation, Metric &metric, + Evaluate &&evaluate, Accepts &&accepts, Better &&better, CheckState &&checkState) { + const auto communicator = m_core.GetCommunicator(); + CoupledRequire(preparation.has_value() && state.normalizedOperator->IsPrepared() && + CoupledFinite(state.trialNormalizedState) && CoupledFinite(state.trialNormalizedResidual) && + std::isfinite(metric.merit), communicator, + "Coupled correction requires a successfully prepared original trial."); + m_diagnostics.originalMerit = metric.merit; + if (CoupledAll(accepts(metric), communicator)) { + m_diagnostics.status = CoupledTrialAttemptStatus::original_accepted; + return; + } + if (m_attemptsThisIteration >= m_maximumAttempts) { + m_diagnostics.status = CoupledTrialAttemptStatus::budget_exhausted; + return; + } + m_originalState = state.trialNormalizedState; + m_originalResidual = state.trialNormalizedResidual; + m_originalPhysical = state.candidatePhysicalState; + Refresh(state); + ++m_attemptsThisIteration; + ++m_diagnostics.attempts; + m_diagnostics.currentSolve = m_core.SolveCorrection(alpha, m_originalResidual, m_correction); + const auto &report = *m_diagnostics.currentSolve; + m_diagnostics.iterations += static_cast(report.linearSolve.iterations); + m_diagnostics.solveSeconds += report.linearSolve.solveSeconds; + m_diagnostics.verificationSeconds += report.verificationSeconds; + if (!report.Admissible()) { + ++m_diagnostics.fallbacks; + switch (report.status) { + case CoupledTrialCorrectionStatus::linear_failure: m_diagnostics.status = CoupledTrialAttemptStatus::linear_failure; break; + case CoupledTrialCorrectionStatus::norm_limit: m_diagnostics.status = CoupledTrialAttemptStatus::norm_limit; break; + default: m_diagnostics.status = CoupledTrialAttemptStatus::non_finite; break; + } + // No candidate mutation has occurred; its prepared operator and + // freshly refreshed inverse still describe the original trial. + return; + } + const Metric originalMetric = metric; + bool restorationStarted = false; + auto restore = [&]() { + restorationStarted = true; + state.trialNormalizedState = m_originalState; + state.candidatePhysicalState = m_originalPhysical; + auto restored = Prepare(state); + ++m_diagnostics.restorations; + if (!CoupledAll(restored.has_value() && CoupledFinite(state.trialNormalizedResidual), communicator)) + throw std::runtime_error("A coupled trial fallback could not reprepare its original candidate."); + preparation = std::move(restored); + Refresh(state); + metric = evaluate(state.trialNormalizedResidual); + ++m_diagnostics.fallbacks; + }; + try { + state.trialNormalizedState += m_correction; + state.normalizedOperator->DenormalizeState(state.trialNormalizedState, state.candidatePhysicalState); + if (!CoupledAll(CoupledFinite(state.trialNormalizedState) && CoupledFinite(state.candidatePhysicalState), communicator)) { + m_diagnostics.status = CoupledTrialAttemptStatus::non_finite; + restore(); + return; + } + if (!CoupledAll(checkState(state.trialNormalizedState), communicator)) { + m_diagnostics.status = CoupledTrialAttemptStatus::inadmissible; + restore(); + return; + } + auto corrected = Prepare(state); + if (!CoupledAll(corrected.has_value(), communicator)) { + m_diagnostics.status = CoupledTrialAttemptStatus::inadmissible; + restore(); + return; + } + if (!CoupledAll(CoupledFinite(state.trialNormalizedResidual), communicator)) { + m_diagnostics.status = CoupledTrialAttemptStatus::non_finite; + restore(); + return; + } + const auto correctedMetric = evaluate(state.trialNormalizedResidual); + m_diagnostics.correctedMerit = correctedMetric.merit; + if (!CoupledAll(better(correctedMetric, originalMetric), communicator)) { + m_diagnostics.status = CoupledTrialAttemptStatus::not_improved; + restore(); + return; + } + preparation = std::move(corrected); + metric = correctedMetric; + m_diagnostics.status = CoupledTrialAttemptStatus::kept; + ++m_diagnostics.kept; + // This is a nonlinear candidate, not an affine gravity completion. + // Its fresh full residual is judged by the caller's normal metric. + // As for an ordinary prepared trial, the caller refreshes the + // inverse only if it commits or needs another auxiliary solve. + } catch (...) { + // Numerical rejection is handled above. For collectively + // propagated exceptions, restore the original prepared trial + // before rethrowing. A failed restoration itself propagates to + // the owning Newton context, which restores the accepted base. + if (!restorationStarted) restore(); + throw; + } + } + + private: + void Refresh(State &state) { + const auto start = Clock::now(); + state.physicalInverse->Refresh(); + m_diagnostics.refreshSeconds += std::chrono::duration(Clock::now() - start).count(); + } + auto Prepare(State &state) { + const auto start = Clock::now(); + auto result = state.PrepareTrial(); + m_diagnostics.preparationSeconds += std::chrono::duration(Clock::now() - start).count(); + ++m_diagnostics.preparations; + return result; + } + PreparedCoupledTrialCorrection m_core; + int m_maximumAttempts, m_attemptsThisIteration{0}; + mfem::Vector m_originalState, m_originalResidual, m_originalPhysical, m_correction; + CoupledTrialDiagnostics m_diagnostics{}; + }; +} diff --git a/libmeanfield/interface/solver/newton.cppm b/libmeanfield/interface/solver/newton.cppm index 2a13485..4be7dcf 100644 --- a/libmeanfield/interface/solver/newton.cppm +++ b/libmeanfield/interface/solver/newton.cppm @@ -268,6 +268,9 @@ export namespace mean_field::solver::nonlinear { // A curved trial path cannot recover its Newton tangent by dividing // the candidate displacement by alpha. Observers borrow it directly. std::span normalizedDirection{}; + // Auxiliary trial solves must remain distinct from the Newton report. + // Like the state spans, this diagnostic view is borrowed for the callback. + const NumericalStrategyDiagnostics *numericalStrategies{nullptr}; }; struct AfterIteration final { diff --git a/libmeanfield/interface/solver/newton_numerics.cppm b/libmeanfield/interface/solver/newton_numerics.cppm index 027d8c9..618ccc3 100644 --- a/libmeanfield/interface/solver/newton_numerics.cppm +++ b/libmeanfield/interface/solver/newton_numerics.cppm @@ -5,6 +5,7 @@ module; #include #include #include +#include #include #include @@ -12,6 +13,7 @@ export module mean_field:solver.newton_numerics; export import :solver.enthalpy_admissibility; export import :solver.linear_backend; +export import :solver.coupled_trial_correction; export namespace mean_field::solver::nonlinear { struct SampledEnthalpy final { }; @@ -39,6 +41,18 @@ export namespace mean_field::solver::nonlinear { GravitySolveOptions gravity{}; }; + struct UncorrectedTrial final { }; + struct CoupledCurvatureCorrection final { + CoupledTrialCorrectionOptions correction{}; + int maximumAttemptsPerIteration{1}; + + void Validate() const { + correction.Validate(); + if (maximumAttemptsPerIteration <= 0) + throw std::invalid_argument("A coupled trial policy requires a positive attempt budget."); + } + }; + template concept EnthalpyNumerics = std::same_as || std::same_as; @@ -48,23 +62,33 @@ export namespace mean_field::solver::nonlinear { std::same_as || std::same_as; + template + concept TrialPathNumerics = 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 + template + requires (!std::same_as || std::same_as) struct NewtonNumerics final { using EnthalpyType = Enthalpy; using GravityType = Gravity; + using TrialType = Trial; [[no_unique_address]] Enthalpy enthalpy{}; [[no_unique_address]] Gravity gravity{}; + [[no_unique_address]] Trial trial{}; }; template NewtonNumerics(Enthalpy, Gravity) -> NewtonNumerics; + template + NewtonNumerics(Enthalpy, Gravity, Trial) -> NewtonNumerics; template struct IsNewtonNumerics : std::false_type { }; - template - struct IsNewtonNumerics> : std::true_type { }; + template + struct IsNewtonNumerics> : std::true_type { }; template concept NewtonNumericalConfiguration = IsNewtonNumerics>::value; @@ -102,6 +126,9 @@ export namespace mean_field::solver::nonlinear { double initialGravityResidual{0.0}; double completedInitialGravityResidual{0.0}; double lastLiftedRelativeTrueResidual{0.0}; + // Gravity diagnostics refer to the preceding gravity-only preparation + // when this transaction subsequently keeps a fully coupled candidate. + CoupledTrialDiagnostics coupledTrial{}; }; } diff --git a/libmeanfield/interface/solver/stellar_context.cppm b/libmeanfield/interface/solver/stellar_context.cppm index 92b3689..f2a629a 100644 --- a/libmeanfield/interface/solver/stellar_context.cppm +++ b/libmeanfield/interface/solver/stellar_context.cppm @@ -1359,6 +1359,7 @@ export namespace mean_field::solver { const Clock::time_point lineSearchStart = Clock::now(); try { for (int trial = 0; trial < options.backtracking.maximumTrials; ++trial) { + strategy.BeginTrial(diagnostics.numericalStrategies); ++lineSearchTrials; ++diagnostics.totalLineSearchTrials; nonlinear::LineSearchTrialDisposition trialDisposition = @@ -1497,6 +1498,41 @@ export namespace mean_field::solver { metricSeconds = std::chrono::duration(Clock::now() - metricStart).count(); timings.metricEvaluationSeconds += metricSeconds; diagnostics.totalMetricEvaluationSeconds += metricSeconds; + const auto acceptsTrialMetric = [&](const nonlinear::MetricEvaluation &candidateMetric) { + if constexpr (Strategy::altersDirection) { + return candidateMetric.merit <= acceptedMetric.merit + + options.backtracking.sufficientDecrease * stepLength * + diagnostics.numericalStrategies.lastMeritDirectionalDerivative; + } else { + return AcceptStep(acceptedMetric, candidateMetric, stepLength, options.backtracking); + } + }; + if constexpr (Strategy::correctsTrial) { + if (detail::AllStellarRanksSatisfy(ValidMetric(trialMetric), communicator)) { + const Clock::time_point correctionStart = Clock::now(); + const double previousMetricSeconds = metricSeconds; + strategy.TryImproveTrial(state, stepLength, preparationResult, trialMetric, + diagnostics.numericalStrategies, + [&](const mfem::Vector &residual) { + const Clock::time_point start = Clock::now(); + auto metric = EvaluateMetric(residual); + const double seconds = std::chrono::duration(Clock::now() - start).count(); + metricSeconds += seconds; + timings.metricEvaluationSeconds += seconds; + diagnostics.totalMetricEvaluationSeconds += seconds; + return metric; + }, acceptsTrialMetric, + [&](const nonlinear::MetricEvaluation &candidate, const nonlinear::MetricEvaluation &original) { + return ValidMetric(candidate) && candidate.merit < original.merit; + }); + const double seconds = std::chrono::duration(Clock::now() - correctionStart).count() - + (metricSeconds - previousMetricSeconds); + preparationSeconds += seconds; + timings.trialPreparationSeconds += seconds; + diagnostics.totalTrialPreparationSeconds += seconds; + minimumJacobianDeterminant = preparationResult->minimumJacobianDeterminant; + } + } observedTrialMetric = trialMetric; if (!detail::AllStellarRanksSatisfy(ValidMetric(trialMetric), communicator)) { rejectionReason = StellarEquilibriumFailureReason::non_finite_residual; @@ -1508,15 +1544,7 @@ export namespace mean_field::solver { ++diagnostics.nonFiniteLineSearchTrials; } else if ( detail::AllStellarRanksSatisfy( - [&] { - if constexpr (Strategy::altersDirection) { - return trialMetric.merit <= acceptedMetric.merit + - options.backtracking.sufficientDecrease * stepLength * - diagnostics.numericalStrategies.lastMeritDirectionalDerivative; - } else { - return AcceptStep(acceptedMetric, trialMetric, stepLength, options.backtracking); - } - }(), + acceptsTrialMetric(trialMetric), communicator ) ) { @@ -1559,7 +1587,8 @@ export namespace mean_field::solver { timings.observerSeconds += NotifyLineSearchTrial( iteration, trial, stepLength, trialDisposition, rejectionSource, observedTrialMetric, - minimumJacobianDeterminant, preparationSeconds, metricSeconds, residualAvailable, state + minimumJacobianDeterminant, preparationSeconds, metricSeconds, residualAvailable, state, + diagnostics.numericalStrategies ); if (accepted) { @@ -1767,7 +1796,8 @@ export namespace mean_field::solver { const double preparationSeconds, const double metricSeconds, const bool residualAvailable, - const State &state + const State &state, + const nonlinear::NumericalStrategyDiagnostics &numericalStrategies ) { if constexpr (nonlinear::detail::ObservesLineSearchTrial) { const nonlinear::AfterLineSearchTrial event{ @@ -1786,7 +1816,8 @@ export namespace mean_field::solver { .candidateNormalizedResidual = residualAvailable ? detail::ReadOnlySpan(state.trialNormalizedResidual) : std::span{}, - .normalizedDirection = detail::ReadOnlySpan(state.normalizedCorrection) + .normalizedDirection = detail::ReadOnlySpan(state.normalizedCorrection), + .numericalStrategies = std::addressof(numericalStrategies) }; const Clock::time_point observerStart = Clock::now(); nonlinear::detail::InvokeAfterLineSearchTrial(m_observer, event); diff --git a/libmeanfield/interface/solver/stellar_newton_strategy.cppm b/libmeanfield/interface/solver/stellar_newton_strategy.cppm index f1d569d..07170b5 100644 --- a/libmeanfield/interface/solver/stellar_newton_strategy.cppm +++ b/libmeanfield/interface/solver/stellar_newton_strategy.cppm @@ -43,8 +43,9 @@ namespace mean_field::solver::detail { 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{}; + using Trial = typename Recipe::TrialType; + std::array integers{}; + std::array reals{}; integers[0] = std::same_as; integers[1] = std::same_as ? 1 : std::same_as ? 2 : 0; @@ -57,15 +58,27 @@ namespace mean_field::solver::detail { integers[4] = options.control.maximumIterations; integers[5] = options.backend.restartLength; integers[6] = options.backend.printLevel; - reals = {options.control.relativeTolerance, options.control.absoluteTolerance}; + reals[0] = options.control.relativeTolerance; + reals[1] = 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) { + integers[8] = std::same_as; + if constexpr (std::same_as) { + const auto &options = recipe.trial.correction; + integers[9] = options.control.maximumIterations; + integers[10] = options.backend.restartLength; + integers[11] = options.backend.printLevel; + integers[12] = recipe.trial.maximumAttemptsPerIteration; + reals[2] = options.control.relativeTolerance; + reals[3] = options.control.absoluteTolerance; + reals[4] = options.maximumCorrectionRatio; + } + integers[7] = std::all_of(reals.begin(), reals.end(), [](double value) { return std::isfinite(value); }); + std::array minimumIntegers{}, maximumIntegers{}; + std::array minimumReals{}, maximumReals{}; + if (MPI_Allreduce(integers.data(), minimumIntegers.data(), 13, MPI_INT, MPI_MIN, communicator) != MPI_SUCCESS || + MPI_Allreduce(integers.data(), maximumIntegers.data(), 13, MPI_INT, MPI_MAX, communicator) != MPI_SUCCESS || + MPI_Allreduce(reals.data(), minimumReals.data(), 5, MPI_DOUBLE, MPI_MIN, communicator) != MPI_SUCCESS || + MPI_Allreduce(reals.data(), maximumReals.data(), 5, MPI_DOUBLE, MPI_MAX, communicator) != MPI_SUCCESS) { throw std::runtime_error("Newton could not compare numerical strategies across ranks."); } if (minimumIntegers != maximumIntegers || minimumReals != maximumReals || @@ -73,6 +86,7 @@ namespace mean_field::solver::detail { throw std::invalid_argument("Newton numerical strategies and controls must agree on every rank."); } if constexpr (!std::same_as) { recipe.gravity.gravity.Validate(); } + if constexpr (std::same_as) { recipe.trial.Validate(); } } template class PreparedGravityStrategy; @@ -91,7 +105,9 @@ namespace mean_field::solver::detail { void CopyDiagnostics(nonlinear::NumericalStrategyDiagnostics &) const { } }; - template + // Export the prepared workspace so bounded diagnostics can exercise this + // same gravity completion and arithmetic gate without reproducing them. + export template class StellarGravityWorkspace { public: using Operation = std::remove_cvref_t().normalizedOperator)>; @@ -406,15 +422,34 @@ namespace mean_field::solver::detail { PreparedStellarEnthalpyStrategy m_prepared; }; + template struct PreparedTrialPath; + template struct PreparedTrialPath { + PreparedTrialPath(State &, const nonlinear::UncorrectedTrial &) { } + }; + template struct PreparedTrialPath { + PreparedTrialPath(State &state, const nonlinear::CoupledCurvatureCorrection &policy) + : transaction(state, policy.correction, policy.maximumAttemptsPerIteration) { } + PreparedCoupledTrialTransaction transaction; + }; + template class PreparedNewtonStrategy final : public NewtonStrategyStorage { public: using Enthalpy = typename Recipe::EnthalpyType; using Gravity = typename Recipe::GravityType; + using Trial = typename Recipe::TrialType; static constexpr bool altersDirection = std::same_as; + static constexpr bool correctsTrial = std::same_as; PreparedNewtonStrategy(State &state, const Recipe &recipe) - : m_enthalpy(state, recipe.enthalpy), m_gravity(state, recipe.gravity) { } + : m_enthalpy(state, recipe.enthalpy), m_gravity(state, recipe.gravity), m_trial(state, recipe.trial) { } + + void BeginTrial(nonlinear::NumericalStrategyDiagnostics &diagnostics) { + if constexpr (correctsTrial) { + m_trial.transaction.BeginTrial(); + diagnostics.coupledTrial = m_trial.transaction.GetDiagnostics(); + } + } bool CheckState(const mfem::Vector &state, nonlinear::NumericalStrategyDiagnostics &diagnostics) { return m_enthalpy.Check(state, diagnostics); @@ -451,6 +486,12 @@ namespace mean_field::solver::detail { if (!std::isfinite(diagnostics.lastMeritDirectionalDerivative) || diagnostics.lastMeritDirectionalDerivative >= 0.0) { return false; } } + if constexpr (correctsTrial) { + // Repair has finished and the accepted Jacobian is still current. + // Every later trial, including a gravity-corrected one, uses this + // same full affine target and the original-direction norm cap. + m_trial.transaction.CaptureBase(state.acceptedNormalizedResidual, state.normalizedCorrection); + } return true; } auto PrepareTrial(State &state, const double alpha, nonlinear::NumericalStrategyDiagnostics &diagnostics) { @@ -460,9 +501,19 @@ namespace mean_field::solver::detail { } bool TrialSolveSucceeded() const noexcept { return m_gravity.TrialSolveSucceeded(); } bool DirectionSolveFailed() const noexcept { return m_directionSolveFailed; } + template + void TryImproveTrial(State &state, double alpha, Preparation &preparation, Metric &metric, + nonlinear::NumericalStrategyDiagnostics &diagnostics, + Evaluate &&evaluate, Accepts &&accepts, Better &&better) requires(correctsTrial) { + m_trial.transaction.TryImprove(state, alpha, preparation, metric, + std::forward(evaluate), std::forward(accepts), std::forward(better), + [&](const mfem::Vector &candidate) { return m_enthalpy.Check(candidate, diagnostics); }); + diagnostics.coupledTrial = m_trial.transaction.GetDiagnostics(); + } private: bool m_directionSolveFailed{false}; PreparedEnthalpyStrategy m_enthalpy; PreparedGravityStrategy m_gravity; + [[no_unique_address]] PreparedTrialPath m_trial; }; } diff --git a/tests/solver/coupled_trial_correction.cpp b/tests/solver/coupled_trial_correction.cpp new file mode 100644 index 0000000..4e0f411 --- /dev/null +++ b/tests/solver/coupled_trial_correction.cpp @@ -0,0 +1,791 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +import mean_field; + +namespace coupled_trial_correction_test { + using namespace mean_field; + using namespace mean_field::solver; + + class DenseOperation final : public mfem::Operator { + public: + explicit DenseOperation(mfem::DenseMatrix matrix) + : mfem::Operator(matrix.Height(), matrix.Width()), m_matrix(std::move(matrix)) { } + + void Mult(const mfem::Vector &input, mfem::Vector &output) const override { + if (!prepared) throw std::logic_error("The test operator is not prepared."); + output.SetSize(Height()); + m_matrix.Mult(input, output); + ++applications; + } + void DenormalizeState(const mfem::Vector &normalized, mfem::Vector &physical) const { physical = normalized; } + void Prepare(mfem::DenseMatrix matrix) { + m_matrix = std::move(matrix); + prepared = true; + ++generation; + } + [[nodiscard]] bool IsPrepared() const noexcept { return prepared; } + [[nodiscard]] std::uint64_t GetPreparationGeneration() const noexcept { return generation; } + [[nodiscard]] const DenseOperation &GetProblem() const noexcept { return *this; } + [[nodiscard]] const mfem::DenseMatrix &Matrix() const noexcept { return m_matrix; } + + bool prepared{true}; + mutable int applications{0}; + std::uint64_t generation{1}; + private: + mfem::DenseMatrix m_matrix; + }; + + class ExactInverse final : public mfem::Solver { + public: + static constexpr auto applicationContract = preconditioning::ApplicationContract::stationary_linear; + explicit ExactInverse(const DenseOperation &operation) + : mfem::Solver(operation.Width()), m_operation(operation) { Refresh(); } + void SetOperator(const mfem::Operator &) override { } + void Mult(const mfem::Vector &input, mfem::Vector &output) const override { + if (!IsCurrent()) throw std::logic_error("The test inverse is stale."); + output.SetSize(Height()); + if (Height() > 0) { + mfem::DenseMatrixInverse inverse(m_operation.Matrix()); + inverse.Mult(input, output); + if (Height() > 1) output(1) += extraSecondComponent * input(0); + } + ++applications; + } + void Refresh() { m_generation = m_operation.GetPreparationGeneration(); } + [[nodiscard]] bool IsCurrent() const noexcept { + return m_operation.IsPrepared() && m_generation == m_operation.GetPreparationGeneration(); + } + mutable int applications{0}; + double extraSecondComponent{0.0}; + private: + const DenseOperation &m_operation; + std::uint64_t m_generation{}; + }; + + [[nodiscard]] mfem::DenseMatrix identity(const int size) { + mfem::DenseMatrix matrix(size); + matrix = 0.0; + for (int i = 0; i < size; ++i) matrix(i, i) = 1.0; + return matrix; + } + + [[nodiscard]] mfem::Vector coupledResidual(const mfem::Vector &state) { + mfem::Vector result(3); + result(0) = 1.0 + state(0) + 10.0 * state(1) * state(2); + result(1) = 1.0 + state(1); + result(2) = 1.0 + state(2); + return result; + } + + [[nodiscard]] mfem::DenseMatrix coupledJacobian(const mfem::Vector &state) { + auto matrix = identity(3); + matrix(0, 1) = 10.0 * state(2); + matrix(0, 2) = 10.0 * state(1); + return matrix; + } + + [[nodiscard]] CoupledTrialCorrectionOptions exactOptions() { + return {.control = {.relativeTolerance = 1.e-12, .absoluteTolerance = 1.e-14, .maximumIterations = 6}, + .backend = {.restartLength = 4, .printLevel = -1}, .maximumCorrectionRatio = 4.0}; + } + + enum class BackendBehavior { exact, inaccurate, failed_status, non_finite }; + struct BackendProbe { + BackendBehavior behavior{BackendBehavior::exact}; + int calls{0}; + bool allInitialGuessesZero{true}; + }; + struct TestBackend final : LinearBackendConfigurationTag { + static constexpr auto supportedPreconditionerContract = preconditioning::ApplicationContract::flexible; + std::shared_ptr probe; + explicit TestBackend(std::shared_ptr value) : probe(std::move(value)) { } + }; + + template + class PreparedTestBackend final { + public: + PreparedTestBackend(const Operation &operation, Inverse &inverse, const MPI_Comm communicator, + std::shared_ptr probe) + : m_operation(operation), m_inverse(inverse), m_communicator(communicator), m_probe(std::move(probe)) { } + [[nodiscard]] const Operation &GetOperator() const noexcept { return m_operation; } + [[nodiscard]] const Inverse &GetPreconditioner() const noexcept { return m_inverse; } + [[nodiscard]] MPI_Comm GetCommunicator() const noexcept { return m_communicator; } + [[nodiscard]] bool IsReady() const noexcept { return m_communicator != MPI_COMM_NULL; } + [[nodiscard]] int RightHandSideSize() const noexcept { return m_operation.Height(); } + [[nodiscard]] int CorrectionSize() const noexcept { return m_operation.Width(); } + [[nodiscard]] LinearSolveReport Solve(const mfem::Vector &rhs, mfem::Vector &correction, + const LinearSolveControl &control) { + ++m_probe->calls; + for (int i = 0; i < correction.Size(); ++i) + m_probe->allInitialGuessesZero &= correction(i) == 0.0; + m_inverse.Mult(rhs, correction); + if (m_probe->behavior == BackendBehavior::inaccurate) correction *= 0.5; + if (m_probe->behavior == BackendBehavior::non_finite && correction.Size() > 0) + correction(0) = std::numeric_limits::quiet_NaN(); + double local = rhs * rhs, global = 0.0; + if (MPI_Allreduce(&local, &global, 1, MPI_DOUBLE, MPI_SUM, m_communicator) != MPI_SUCCESS) + throw std::runtime_error("The test backend could not reduce its norm."); + // A backend claiming exact convergence is deliberately insufficient + // evidence: the correction owner must verify the actual J action. + return {.status = m_probe->behavior == BackendBehavior::failed_status ? + LinearSolveStatus::maximum_iterations : LinearSolveStatus::converged, + .control = control, .iterations = 1, .rightHandSideNorm = std::sqrt(global), + .trueResidualNorm = 0.0, .relativeTrueResidualNorm = 0.0}; + } + private: + const Operation &m_operation; + Inverse &m_inverse; + MPI_Comm m_communicator; + std::shared_ptr m_probe; + }; + + template + [[nodiscard]] auto prepareLinearBackend(TestBackend recipe, const Operation &operation, Inverse &inverse, + const MPI_Comm communicator) { + return PreparedTestBackend(operation, inverse, communicator, std::move(recipe.probe)); + } + + static_assert(CoupledTrialOperation); + static_assert(CoupledTrialInverse); + static_assert(!CoupledTrialOperation); + static_assert(!CoupledTrialInverse); + static_assert(LinearBackendRuntimeAvailableFor); + using ScriptedCorrection = PreparedCoupledTrialCorrection; +} // namespace coupled_trial_correction_test + +TEST_CASE("A coupled trial correction preserves the tangent and removes mixed nonlinear curvature", + "[solver][coupled_trial_correction]") { + using namespace coupled_trial_correction_test; + using Catch::Approx; + mfem::Vector base(3), direction(3), correction(3), prediction(3); + base = 0.0; + direction = -1.0; + const auto baseResidual = coupledResidual(base); + DenseOperation operation(coupledJacobian(base)); + ExactInverse inverse(operation); + PreparedCoupledTrialCorrection prepared(operation, inverse, MPI_COMM_SELF, exactOptions()); + prepared.CaptureBase(baseResidual, direction); + const auto capturedGeneration = prepared.GetBaseGeneration(); + CHECK(prepared.GetBaseDirectionNorm() == Approx(std::sqrt(3.0))); + double previousScaledDisplacement = 0.0; + for (const double alpha : {.5, .25, .125, .0625}) { + CAPTURE(alpha); + mfem::Vector trial(direction); + trial *= alpha; + const auto rawResidual = coupledResidual(trial); + operation.Prepare(coupledJacobian(trial)); + inverse.Refresh(); + correction = 17.0; + const auto report = prepared.SolveCorrection(alpha, rawResidual, correction); + REQUIRE(report.Admissible()); + REQUIRE(report.verifiedGate); + CHECK(report.linearSolve.Converged()); + CHECK(report.linearSolve.iterations <= exactOptions().control.maximumIterations); + CHECK(report.verifiedTrueResidualNorm <= report.convergenceThreshold); + CHECK(correction(0) == Approx(-10.0 * alpha * alpha).margin(1.e-13)); + CHECK(std::abs(correction(1)) < 1.e-13); + CHECK(std::abs(correction(2)) < 1.e-13); + trial += correction; + const auto correctedResidual = coupledResidual(trial); + prepared.BuildLinearPrediction(alpha, prediction); + mfem::Vector remainder(correctedResidual); + remainder -= prediction; + CHECK(remainder.Norml2() < 1.e-12); + CHECK(correctedResidual.Norml2() < rawResidual.Norml2()); + if (alpha == .5) { + CHECK(rawResidual.Norml2() > baseResidual.Norml2()); + CHECK(correctedResidual.Norml2() == Approx(.5 * baseResidual.Norml2())); + } + trial /= alpha; + trial -= direction; + const double scaledDisplacement = trial.Norml2(); + CHECK(scaledDisplacement == Approx(10.0 * alpha)); + if (previousScaledDisplacement > 0.0) + CHECK(scaledDisplacement == Approx(previousScaledDisplacement / 2.0)); + previousScaledDisplacement = scaledDisplacement; + CHECK(prepared.GetBaseGeneration() == capturedGeneration); + for (int i = 0; i < 3; ++i) CHECK(direction(i) == -1.0); + } +} + +TEST_CASE("Coupled correction captures the base action while solving with the current trial Jacobian", + "[solver][coupled_trial_correction][lifetime]") { + using namespace coupled_trial_correction_test; + using Catch::Approx; + DenseOperation operation(identity(2)); + ExactInverse inverse(operation); + auto probe = std::make_shared(); + ScriptedCorrection prepared(operation, inverse, MPI_COMM_SELF, exactOptions(), TestBackend(probe)); + mfem::Vector f(2), p(2), prediction(2), correction(2), trialResidual(2); + f = 1.0; + p = -1.0; + prepared.CaptureBase(f, p); + const auto generation = prepared.GetBaseGeneration(); + auto matrix = identity(2); + matrix(0, 0) = 2.0; + matrix(1, 1) = 4.0; + operation.Prepare(matrix); + const int applications = operation.applications; + prepared.BuildLinearPrediction(.5, prediction); + CHECK(operation.applications == applications); + CHECK(prediction(0) == .5); + CHECK(prediction(1) == .5); + trialResidual = prediction; + trialResidual(0) += .2; + trialResidual(1) += .4; + REQUIRE_THROWS_AS(prepared.SolveCorrection(.5, trialResidual, correction), std::invalid_argument); + inverse.Refresh(); + const auto report = prepared.SolveCorrection(.5, trialResidual, correction); + REQUIRE(report.Admissible()); + CHECK(correction(0) == Approx(-.1)); + CHECK(correction(1) == Approx(-.1)); + CHECK(prepared.GetBaseGeneration() == generation); + CHECK(probe->allInitialGuessesZero); + prepared.Invalidate(); + CHECK_FALSE(prepared.HasBase()); + CHECK_THROWS_AS(prepared.BuildLinearPrediction(.5, prediction), std::invalid_argument); +} + +TEST_CASE("Coupled correction independently rejects a false auxiliary convergence report", + "[solver][coupled_trial_correction][failure]") { + using namespace coupled_trial_correction_test; + DenseOperation operation(identity(2)); + ExactInverse inverse(operation); + auto probe = std::make_shared(); + ScriptedCorrection prepared(operation, inverse, MPI_COMM_SELF, exactOptions(), TestBackend(probe)); + mfem::Vector f(2), p(2), trial(2), correction(2); + f = 1.0; + p = -1.0; + trial = .75; + prepared.CaptureBase(f, p); + SECTION("The returned correction fails the independently computed true gate") { + probe->behavior = BackendBehavior::inaccurate; + const auto report = prepared.SolveCorrection(.5, trial, correction); + CHECK_FALSE(report.Admissible()); + CHECK(report.status == CoupledTrialCorrectionStatus::linear_failure); + CHECK(report.linearSolve.Converged()); + CHECK_FALSE(report.verifiedGate); + CHECK(report.verifiedTrueResidualNorm > report.convergenceThreshold); + } + SECTION("An accurate correction does not erase backend nonconvergence") { + probe->behavior = BackendBehavior::failed_status; + const auto report = prepared.SolveCorrection(.5, trial, correction); + CHECK_FALSE(report.Admissible()); + CHECK(report.status == CoupledTrialCorrectionStatus::linear_failure); + CHECK_FALSE(report.linearSolve.Converged()); + } + SECTION("A nonfinite returned correction is never admissible") { + probe->behavior = BackendBehavior::non_finite; + const auto report = prepared.SolveCorrection(.5, trial, correction); + CHECK_FALSE(report.Admissible()); + CHECK(report.status == CoupledTrialCorrectionStatus::non_finite); + CHECK_FALSE(report.verifiedGate); + } + CHECK(probe->calls == 1); + CHECK(probe->allInitialGuessesZero); + for (int i = 0; i < 2; ++i) { + CHECK(f(i) == 1.0); + CHECK(p(i) == -1.0); + CHECK(trial(i) == .75); + } +} + +TEST_CASE("Coupled correction bounds work and correction size without changing the Newton direction", + "[solver][coupled_trial_correction][safeguard]") { + using namespace coupled_trial_correction_test; + using Catch::Approx; + DenseOperation operation(identity(2)); + ExactInverse inverse(operation); + auto options = exactOptions(); + options.maximumCorrectionRatio = .1; + auto probe = std::make_shared(); + ScriptedCorrection prepared(operation, inverse, MPI_COMM_SELF, options, TestBackend(probe)); + mfem::Vector f(2), p(2), trial(2), correction(2); + f = 1.0; + p = -1.0; + trial = .9; + prepared.CaptureBase(f, p); + const auto report = prepared.SolveCorrection(.5, trial, correction); + CHECK(report.status == CoupledTrialCorrectionStatus::norm_limit); + CHECK_FALSE(report.Admissible()); + CHECK(report.verifiedGate); + CHECK(report.maximumCorrectionNorm == Approx(.1 * .5 * std::sqrt(2.0))); + CHECK(report.correctionNorm > report.maximumCorrectionNorm); + CHECK(report.linearSolve.iterations <= options.control.maximumIterations); + CHECK(probe->calls == 1); + CHECK(probe->allInitialGuessesZero); + CHECK(p(0) == -1.0); + CHECK(p(1) == -1.0); +} + +TEST_CASE("Zero coupled remainder does not manufacture an auxiliary displacement", + "[solver][coupled_trial_correction][safeguard]") { + using namespace coupled_trial_correction_test; + DenseOperation operation(identity(2)); + ExactInverse inverse(operation); + auto probe = std::make_shared(); + ScriptedCorrection prepared(operation, inverse, MPI_COMM_SELF, exactOptions(), TestBackend(probe)); + mfem::Vector f(2), p(2), trial(2), correction(2); + f = 1.0; + p = -1.0; + trial = .5; + correction = 900.0; + prepared.CaptureBase(f, p); + const auto report = prepared.SolveCorrection(.5, trial, correction); + CHECK(report.Admissible()); + CHECK(report.verifiedGate); + CHECK(correction.Norml2() == 0.0); + CHECK(report.correctionNorm == 0.0); + CHECK(report.verifiedTrueResidualNorm == 0.0); + CHECK(probe->calls <= 1); + CHECK(probe->allInitialGuessesZero); +} + +TEST_CASE("Coupled correction rejects invalid controls before numerical work", + "[solver][coupled_trial_correction][configuration]") { + using namespace coupled_trial_correction_test; + auto options = exactOptions(); + SECTION("A nonfinite cap cannot define a bounded correction") { + options.maximumCorrectionRatio = std::numeric_limits::infinity(); + } + SECTION("A negative cap has no feasible correction size") { options.maximumCorrectionRatio = -1.0; } + SECTION("A zero cap cannot enable the correction strategy") { options.maximumCorrectionRatio = 0.0; } + SECTION("The auxiliary iteration budget must be finite and positive") { options.control.maximumIterations = 0; } + SECTION("The auxiliary tolerance cannot be NaN") { + options.control.relativeTolerance = std::numeric_limits::quiet_NaN(); + } + CHECK_THROWS_AS(options.Validate(), std::invalid_argument); +} + +TEST_CASE("Coupled correction decisions remain collective with an empty peer rank", + "[solver][coupled_trial_correction][mpi]") { + using namespace coupled_trial_correction_test; + using Catch::Approx; + int rank = 0; + REQUIRE(MPI_Comm_rank(MPI_COMM_WORLD, &rank) == MPI_SUCCESS); + // Only one rank owns this tiny synthetic problem. A distributed stop must + // not rely on every rank having a local entry to inspect. + const int localSize = rank == 0 ? 2 : 0; + DenseOperation operation(identity(localSize)); + ExactInverse inverse(operation); + auto probe = std::make_shared(); + ScriptedCorrection prepared(operation, inverse, MPI_COMM_WORLD, exactOptions(), TestBackend(probe)); + mfem::Vector f(localSize), p(localSize), trial(localSize), correction(localSize); + f = 1.0; + p = -1.0; + trial = .75; + prepared.CaptureBase(f, p); + SECTION("Finite global norms include the owning rank exactly once") { + const auto report = prepared.SolveCorrection(.5, trial, correction); + CHECK(report.Admissible()); + CHECK(report.correctionNorm == Approx(std::sqrt(2.0) * .25)); + CHECK(probe->calls == 1); + } + SECTION("One nonfinite trial prevents the backend call on every rank") { + if (rank == 0) trial(0) = std::numeric_limits::quiet_NaN(); + const auto report = prepared.SolveCorrection(.5, trial, correction); + CHECK(report.status == CoupledTrialCorrectionStatus::non_finite); + CHECK_FALSE(report.Admissible()); + CHECK(probe->calls == 0); + } + SECTION("One nonfinite correction fails the independently checked global gate") { + probe->behavior = BackendBehavior::non_finite; + const auto report = prepared.SolveCorrection(.5, trial, correction); + CHECK(report.status == CoupledTrialCorrectionStatus::non_finite); + CHECK_FALSE(report.Admissible()); + CHECK_FALSE(report.verifiedGate); + CHECK(probe->calls == 1); + } +} + +namespace coupled_trial_correction_test { + struct ToyProblem final { + MPI_Comm communicator{MPI_COMM_SELF}; + [[nodiscard]] MPI_Comm GetCommunicator() const noexcept { return communicator; } + }; + + struct ToyState final { + DenseOperation operation{identity(3)}; + ExactInverse inverse{operation}; + DenseOperation *normalizedOperator{&operation}; + ExactInverse *normalizedInverse{&inverse}; + ExactInverse *physicalInverse{&inverse}; + mfem::Vector trialNormalizedState = mfem::Vector(3); + mfem::Vector trialNormalizedResidual = mfem::Vector(3); + mfem::Vector candidatePhysicalState = mfem::Vector(3); + mfem::Vector acceptedState = mfem::Vector(3); + mfem::Vector acceptedResidual = mfem::Vector(3); + mfem::Vector normalizedCorrection = mfem::Vector(3); + LinearSolveReport outerReport{.status = LinearSolveStatus::converged, .iterations = 11, + .trueResidualNorm = .0123}; + ToyProblem problem{}; + int preparations{0}; + bool quarticBarrier{false}, rejectCorrectedPreparation{false}, nonfiniteCorrectedResidual{false}; + bool throwCorrectedPreparation{false}; + + ToyState() { + trialNormalizedState = 0.0; + candidatePhysicalState = 0.0; + trialNormalizedResidual = 1.0; + acceptedState = 0.0; + acceptedResidual = 1.0; + normalizedCorrection = -1.0; + } + [[nodiscard]] ToyProblem &Problem() noexcept { return problem; } + [[nodiscard]] const ToyProblem &Problem() const noexcept { return problem; } + [[nodiscard]] std::expected PrepareTrial() { + ++preparations; + candidatePhysicalState = trialNormalizedState; + if (throwCorrectedPreparation && trialNormalizedState(0) < -2.0) { + operation.prepared = false; + throw std::runtime_error("Synthetic corrected preparation failure."); + } + if (rejectCorrectedPreparation && trialNormalizedState(0) < -2.0) { + operation.prepared = false; + return std::unexpected(7); + } + trialNormalizedResidual = coupledResidual(trialNormalizedState); + auto matrix = coupledJacobian(trialNormalizedState); + if (quarticBarrier) { + // The extra term and its derivative vanish at the base and at + // alpha=.5. It makes the solved correction harmful without + // supplying the transaction an inconsistent trial Jacobian. + const double x = trialNormalizedState(0); + trialNormalizedResidual(0) += x*x*(x+.5)*(x+.5); + matrix(0, 0) += 2.0*x*(x+.5)*(2.0*x+.5); + } + operation.Prepare(matrix); + if (nonfiniteCorrectedResidual && trialNormalizedState(0) < -2.0) + trialNormalizedResidual(0) = std::numeric_limits::quiet_NaN(); + return {}; + } + [[nodiscard]] std::expected PrepareAlpha(const double alpha) { + trialNormalizedState = normalizedCorrection; + trialNormalizedState *= alpha; + return PrepareTrial(); + } + }; + + [[nodiscard]] nonlinear::MetricEvaluation toyMetric(const mfem::Vector &residual) { + const double norm = residual.Norml2(); + return {.residualNorm = norm, .merit = .5*norm*norm}; + } + [[nodiscard]] bool toyAccepts(const nonlinear::MetricEvaluation &metric) { + return std::isfinite(metric.merit) && metric.merit < 1.5; + } + [[nodiscard]] bool toyImproves(const nonlinear::MetricEvaluation &candidate, + const nonlinear::MetricEvaluation &original) { + return std::isfinite(candidate.merit) && candidate.merit < original.merit; + } + [[nodiscard]] bool admissibleToy(const mfem::Vector &) { return true; } + + void checkOuterData(const ToyState &state) { + for (int i = 0; i < 3; ++i) { + CHECK(state.acceptedState(i) == 0.0); + CHECK(state.acceptedResidual(i) == 1.0); + CHECK(state.normalizedCorrection(i) == -1.0); + } + CHECK(state.outerReport.status == LinearSolveStatus::converged); + CHECK(state.outerReport.iterations == 11); + CHECK(state.outerReport.trueResidualNorm == .0123); + } + void checkRestoredTrial(const ToyState &state, const mfem::Vector &originalState, + const mfem::Vector &originalResidual, const mfem::Vector &originalAction) { + for (int i = 0; i < 3; ++i) { + CHECK(state.trialNormalizedState(i) == originalState(i)); + CHECK(state.candidatePhysicalState(i) == originalState(i)); + CHECK(state.trialNormalizedResidual(i) == originalResidual(i)); + } + CHECK(state.operation.IsPrepared()); + CHECK(state.inverse.IsCurrent()); + mfem::Vector action(3); + state.operation.Mult(state.normalizedCorrection, action); + for (int i = 0; i < 3; ++i) CHECK(action(i) == originalAction(i)); + checkOuterData(state); + } +} + +TEST_CASE("The production coupled transaction keeps an improved full residual and leaves outer data alone", + "[solver][coupled_trial_correction][transaction]") { + using namespace coupled_trial_correction_test; + using Catch::Approx; + ToyState state; + solver::detail::PreparedCoupledTrialTransaction transaction(state, exactOptions()); + transaction.CaptureBase(state.acceptedResidual, state.normalizedCorrection); + auto preparation = state.PrepareAlpha(.5); + REQUIRE(preparation.has_value()); + auto metric = toyMetric(state.trialNormalizedResidual); + transaction.BeginTrial(); + transaction.TryImprove(state, .5, preparation, metric, toyMetric, toyAccepts, toyImproves, admissibleToy); + const auto &diagnostics = transaction.GetDiagnostics(); + CHECK(diagnostics.status == CoupledTrialAttemptStatus::kept); + REQUIRE(diagnostics.currentSolve.has_value()); + CHECK(diagnostics.currentSolve->Admissible()); + CHECK(diagnostics.currentSolve->verifiedGate); + CHECK(diagnostics.attempts == 1); + CHECK(diagnostics.kept == 1); + CHECK(diagnostics.fallbacks == 0); + CHECK(diagnostics.preparations == 1); + CHECK(diagnostics.restorations == 0); + CHECK(metric.merit == Approx(.375)); + CHECK(metric.merit < diagnostics.originalMerit); + REQUIRE(preparation.has_value()); + for (int i = 0; i < 3; ++i) CHECK(state.trialNormalizedResidual(i) == Approx(.5)); + checkOuterData(state); + transaction.BeginTrial(); + CHECK(transaction.GetDiagnostics().status == CoupledTrialAttemptStatus::not_attempted); + CHECK_FALSE(transaction.GetDiagnostics().currentSolve.has_value()); + CHECK(std::isnan(transaction.GetDiagnostics().originalMerit)); + CHECK(transaction.GetDiagnostics().attempts == 1); +} + +TEST_CASE("The production coupled transaction restores rejected corrections and their prepared operators", + "[solver][coupled_trial_correction][transaction][rollback]") { + using namespace coupled_trial_correction_test; + ToyState state; + bool rejectAdmissibility = false; + CoupledTrialAttemptStatus expected = CoupledTrialAttemptStatus::inadmissible; + SECTION("Corrected physical preparation rejects the candidate") { state.rejectCorrectedPreparation = true; } + SECTION("Corrected polynomial admissibility rejects before preparation") { rejectAdmissibility = true; } + SECTION("A consistent nonlinear map makes the solved correction worse") { + state.quarticBarrier = true; + expected = CoupledTrialAttemptStatus::not_improved; + } + SECTION("A successfully prepared nonfinite residual is unusable") { + state.nonfiniteCorrectedResidual = true; + expected = CoupledTrialAttemptStatus::non_finite; + } + solver::detail::PreparedCoupledTrialTransaction transaction(state, exactOptions()); + transaction.CaptureBase(state.acceptedResidual, state.normalizedCorrection); + auto preparation = state.PrepareAlpha(.5); + REQUIRE(preparation.has_value()); + const mfem::Vector originalState(state.trialNormalizedState), originalResidual(state.trialNormalizedResidual); + mfem::Vector originalAction(3); + state.operation.Mult(state.normalizedCorrection, originalAction); + auto metric = toyMetric(originalResidual); + const double originalMerit = metric.merit; + transaction.BeginTrial(); + transaction.TryImprove(state, .5, preparation, metric, toyMetric, toyAccepts, toyImproves, + [rejectAdmissibility](const mfem::Vector &candidate) { return !rejectAdmissibility || candidate(0) >= -2.0; }); + const auto &diagnostics = transaction.GetDiagnostics(); + CHECK(diagnostics.status == expected); + REQUIRE(diagnostics.currentSolve.has_value()); + CHECK(diagnostics.currentSolve->Admissible()); + CHECK(diagnostics.attempts == 1); + CHECK(diagnostics.kept == 0); + CHECK(diagnostics.fallbacks == 1); + CHECK(diagnostics.restorations == 1); + CHECK(diagnostics.preparations == (rejectAdmissibility ? 1 : 2)); + CHECK(metric.merit == originalMerit); + if (expected == CoupledTrialAttemptStatus::not_improved) + CHECK(diagnostics.correctedMerit > diagnostics.originalMerit); + REQUIRE(preparation.has_value()); + checkRestoredTrial(state, originalState, originalResidual, originalAction); +} + +TEST_CASE("The coupled transaction skips acceptable trials and bounds attempts per captured direction", + "[solver][coupled_trial_correction][transaction][budget]") { + using namespace coupled_trial_correction_test; + ToyState state; + state.quarticBarrier = true; + solver::detail::PreparedCoupledTrialTransaction transaction(state, exactOptions(), 1); + transaction.CaptureBase(state.acceptedResidual, state.normalizedCorrection); + auto preparation = state.PrepareAlpha(.05); + auto metric = toyMetric(state.trialNormalizedResidual); + REQUIRE(toyAccepts(metric)); + const int before = state.preparations; + transaction.BeginTrial(); + transaction.TryImprove(state, .05, preparation, metric, toyMetric, toyAccepts, toyImproves, admissibleToy); + CHECK(transaction.GetDiagnostics().status == CoupledTrialAttemptStatus::original_accepted); + CHECK(transaction.GetDiagnostics().attempts == 0); + CHECK(state.preparations == before); + CHECK_FALSE(transaction.GetDiagnostics().currentSolve.has_value()); + + preparation = state.PrepareAlpha(.5); + metric = toyMetric(state.trialNormalizedResidual); + transaction.BeginTrial(); + transaction.TryImprove(state, .5, preparation, metric, toyMetric, toyAccepts, toyImproves, admissibleToy); + CHECK(transaction.GetDiagnostics().status == CoupledTrialAttemptStatus::not_improved); + CHECK(transaction.GetDiagnostics().attempts == 1); + const auto iterations = transaction.GetDiagnostics().iterations; + transaction.BeginTrial(); + transaction.TryImprove(state, .5, preparation, metric, toyMetric, toyAccepts, toyImproves, admissibleToy); + CHECK(transaction.GetDiagnostics().status == CoupledTrialAttemptStatus::budget_exhausted); + CHECK(transaction.GetDiagnostics().attempts == 1); + CHECK(transaction.GetDiagnostics().iterations == iterations); + CHECK_FALSE(transaction.GetDiagnostics().currentSolve.has_value()); + checkOuterData(state); +} + +TEST_CASE("A rank-local corrected admissibility rejection restores the original trial everywhere", + "[solver][coupled_trial_correction][transaction][rollback][mpi]") { + using namespace coupled_trial_correction_test; + int rank = 0; + REQUIRE(MPI_Comm_rank(MPI_COMM_WORLD, &rank) == MPI_SUCCESS); + ToyState state; + state.problem.communicator = MPI_COMM_WORLD; + solver::detail::PreparedCoupledTrialTransaction transaction(state, exactOptions()); + transaction.CaptureBase(state.acceptedResidual, state.normalizedCorrection); + auto preparation = state.PrepareAlpha(.5); + const mfem::Vector originalState(state.trialNormalizedState), originalResidual(state.trialNormalizedResidual); + mfem::Vector originalAction(3); + state.operation.Mult(state.normalizedCorrection, originalAction); + auto metric = toyMetric(originalResidual); + transaction.BeginTrial(); + transaction.TryImprove(state, .5, preparation, metric, toyMetric, toyAccepts, toyImproves, + [rank](const mfem::Vector &) { return rank != 0; }); + CHECK(transaction.GetDiagnostics().status == CoupledTrialAttemptStatus::inadmissible); + CHECK(transaction.GetDiagnostics().restorations == 1); + checkRestoredTrial(state, originalState, originalResidual, originalAction); +} + +TEST_CASE("A coupled auxiliary gate or cap failure preserves the already prepared original candidate", + "[solver][coupled_trial_correction][transaction][failure]") { + using namespace coupled_trial_correction_test; + ToyState state; + auto options = exactOptions(); + auto expected = CoupledTrialAttemptStatus::norm_limit; + SECTION("A correction outside the cap is rejected without mutation") { + options.maximumCorrectionRatio = .1; + } + SECTION("A bounded unsuccessful Krylov solve is rejected without mutation") { + // This stationary perturbation prevents one Krylov vector from spanning + // the exact correction. The failure comes from a real verified solve. + state.inverse.extraSecondComponent = .1; + options.control.maximumIterations = 1; + options.backend.restartLength = 1; + expected = CoupledTrialAttemptStatus::linear_failure; + } + solver::detail::PreparedCoupledTrialTransaction transaction(state, options); + transaction.CaptureBase(state.acceptedResidual, state.normalizedCorrection); + auto preparation = state.PrepareAlpha(.5); + const mfem::Vector originalState(state.trialNormalizedState), originalResidual(state.trialNormalizedResidual); + mfem::Vector originalAction(3); + state.operation.Mult(state.normalizedCorrection, originalAction); + auto metric = toyMetric(originalResidual); + transaction.BeginTrial(); + transaction.TryImprove(state, .5, preparation, metric, toyMetric, toyAccepts, toyImproves, admissibleToy); + const auto &diagnostics = transaction.GetDiagnostics(); + CHECK(diagnostics.status == expected); + REQUIRE(diagnostics.currentSolve.has_value()); + CHECK_FALSE(diagnostics.currentSolve->Admissible()); + CHECK(diagnostics.attempts == 1); + CHECK(diagnostics.fallbacks == 1); + CHECK(diagnostics.preparations == 0); + CHECK(diagnostics.restorations == 0); + CHECK(diagnostics.iterations <= static_cast(options.control.maximumIterations)); + if (expected == CoupledTrialAttemptStatus::linear_failure) + CHECK_FALSE(diagnostics.currentSolve->verifiedGate); + checkRestoredTrial(state, originalState, originalResidual, originalAction); +} + +TEST_CASE("Coupled correction options and step lengths must agree across ranks", + "[solver][coupled_trial_correction][configuration][mpi]") { + using namespace coupled_trial_correction_test; + int rank = 0, ranks = 0; + REQUIRE(MPI_Comm_rank(MPI_COMM_WORLD, &rank) == MPI_SUCCESS); + REQUIRE(MPI_Comm_size(MPI_COMM_WORLD, &ranks) == MPI_SUCCESS); + DenseOperation operation(identity(2)); + ExactInverse inverse(operation); + auto options = exactOptions(); + SECTION("A rank-local invalid cap fails collectively") { + if (rank == 0) options.maximumCorrectionRatio = std::numeric_limits::quiet_NaN(); + const auto create = [&] { return PreparedCoupledTrialCorrection(operation, inverse, MPI_COMM_WORLD, options); }; + CHECK_THROWS_AS(create(), std::invalid_argument); + } + SECTION("Different valid caps fail collectively") { + if (ranks < 2) return; + options.maximumCorrectionRatio = rank == 0 ? 1.0 : 2.0; + const auto create = [&] { return PreparedCoupledTrialCorrection(operation, inverse, MPI_COMM_WORLD, options); }; + CHECK_THROWS_AS(create(), std::invalid_argument); + } + SECTION("Different step lengths cannot form one global trial") { + if (ranks < 2) return; + PreparedCoupledTrialCorrection prepared(operation, inverse, MPI_COMM_WORLD, options); + mfem::Vector f(2), p(2), prediction(2); + f = 1.0; + p = -1.0; + prepared.CaptureBase(f, p); + CHECK_THROWS_AS(prepared.BuildLinearPrediction(rank == 0 ? .5 : .25, prediction), std::invalid_argument); + } +} + +TEST_CASE("Finite entries with an overflowing auxiliary norm never authorize a correction", + "[solver][coupled_trial_correction][failure]") { + using namespace coupled_trial_correction_test; + DenseOperation operation(identity(2)); + ExactInverse inverse(operation); + auto probe = std::make_shared(); + ScriptedCorrection prepared(operation, inverse, MPI_COMM_SELF, exactOptions(), TestBackend(probe)); + mfem::Vector f(2), p(2), trial(2), correction(2); + f = 1.0; + p = -1.0; + trial = 1.e200; + prepared.CaptureBase(f, p); + const auto report = prepared.SolveCorrection(.5, trial, correction); + CHECK(report.status == CoupledTrialCorrectionStatus::non_finite); + CHECK_FALSE(report.Admissible()); + CHECK_FALSE(report.verifiedGate); + CHECK(probe->calls == 0); +} + +TEST_CASE("Exceptions after coupled mutation restore the original prepared candidate before escaping", + "[solver][coupled_trial_correction][transaction][rollback][exception]") { + using namespace coupled_trial_correction_test; + ToyState state; + enum class FailureSite { admissibility, preparation, metric, comparison }; + FailureSite site = FailureSite::admissibility; + SECTION("The admissibility callback throws") { site = FailureSite::admissibility; } + SECTION("Physical preparation invalidates its operator and throws") { + site = FailureSite::preparation; + state.throwCorrectedPreparation = true; + } + SECTION("Evaluation of the corrected metric throws") { site = FailureSite::metric; } + SECTION("Comparison of the corrected full merit throws") { site = FailureSite::comparison; } + solver::detail::PreparedCoupledTrialTransaction transaction(state, exactOptions()); + transaction.CaptureBase(state.acceptedResidual, state.normalizedCorrection); + auto preparation = state.PrepareAlpha(.5); + const mfem::Vector originalState(state.trialNormalizedState), originalResidual(state.trialNormalizedResidual); + mfem::Vector originalAction(3); + state.operation.Mult(state.normalizedCorrection, originalAction); + auto metric = toyMetric(originalResidual); + const double originalMerit = metric.merit; + int evaluations = 0; + transaction.BeginTrial(); + const auto attempt = [&] { + transaction.TryImprove(state, .5, preparation, metric, + [&](const mfem::Vector &residual) { + if (++evaluations == 1 && site == FailureSite::metric) + throw std::runtime_error("Synthetic corrected metric failure."); + return toyMetric(residual); + }, toyAccepts, + [&](const auto &candidate, const auto &original) { + if (site == FailureSite::comparison) + throw std::runtime_error("Synthetic corrected comparison failure."); + return toyImproves(candidate, original); + }, [&](const mfem::Vector &) { + if (site == FailureSite::admissibility) + throw std::runtime_error("Synthetic corrected admissibility failure."); + return true; + }); + }; + CHECK_THROWS_AS(attempt(), std::runtime_error); + CHECK(transaction.GetDiagnostics().kept == 0); + CHECK(transaction.GetDiagnostics().restorations == 1); + CHECK(metric.merit == originalMerit); + REQUIRE(preparation.has_value()); + checkRestoredTrial(state, originalState, originalResidual, originalAction); +} diff --git a/tests/solver/stellar_numerical_strategies.cpp b/tests/solver/stellar_numerical_strategies.cpp index ac3a37e..ca187ec 100644 --- a/tests/solver/stellar_numerical_strategies.cpp +++ b/tests/solver/stellar_numerical_strategies.cpp @@ -21,6 +21,7 @@ namespace stellar_numerical_strategy_test { struct UnregisteredEnthalpy { }; struct UnregisteredGravity { }; + struct UnregisteredTrial { }; template concept MakesNumerics = requires { typename nonlinear::NewtonNumerics; }; static_assert(!MakesNumerics); @@ -32,6 +33,21 @@ namespace stellar_numerical_strategy_test { nonlinear::NormalizedResidualMetric{})), nonlinear::Newton<>>); static_assert(MakesNumerics); static_assert(MakesNumerics); + template + concept MakesTrialNumerics = requires { typename nonlinear::NewtonNumerics; }; + static_assert(std::same_as::TrialType, nonlinear::UncorrectedTrial>); + static_assert(std::same_as>); + static_assert(MakesTrialNumerics); + static_assert(MakesTrialNumerics); + static_assert(!MakesTrialNumerics); + static_assert(!MakesTrialNumerics); + static_assert(!MakesTrialNumerics); + [[nodiscard]] double matchingCentralDensity() { return std::numbers::pi_v * utils::MASS / (4.0 * utils::RADIUS * utils::RADIUS * utils::RADIUS); @@ -284,4 +300,47 @@ TEST_CASE("Polynomial admissibility is invoked by each production gravity strate check(nonlinear::EliminatedGravity{permissiveGravity()}); REQUIRE(context.isReady()); checkPreserved(original, snapshot(context)); + { + const auto report = evaluate(context, nonlinear::NewtonNumerics{polynomial, + nonlinear::MonolithicGravity{}, nonlinear::CoupledCurvatureCorrection{}}); + CHECK_FALSE(report.converged()); + CHECK(report.diagnostics().numericalStrategies.polynomialChecks > 0); + CHECK(report.diagnostics().acceptedNonlinearIterations == 0); + CHECK(report.lastAcceptedCheckpointView().valid()); + } + REQUIRE(context.isReady()); + checkPreserved(original, snapshot(context)); +} + +TEST_CASE("Coupled trial selection keeps the outer linear report and rejected stellar checkpoint", + "[solver][numerical_strategies][coupled_trial_correction][rollback]") { + using namespace stellar_numerical_strategy_test; + auto context = makeContext(); + const auto original = snapshot(context); + const auto check = [&](auto gravity) { + const nonlinear::CoupledCurvatureCorrection correction{}; + const auto report = evaluate(context, nonlinear::NewtonNumerics{ + nonlinear::SampledEnthalpy{}, gravity, correction}); + CHECK_FALSE(report.converged()); + CHECK(report.diagnostics().acceptedNonlinearIterations == 0); + REQUIRE(report.diagnostics().lastLinearSolve.has_value()); + const auto &outer = *report.diagnostics().lastLinearSolve; + CHECK(outer.iterations == 0); + CHECK(outer.control.absoluteTolerance == 1.e100); + CHECK(report.diagnostics().numericalStrategies.lastLinearReportSpace == nonlinear::LinearReportSpace::full_system); + const auto &auxiliary = report.diagnostics().numericalStrategies.coupledTrial; + REQUIRE(auxiliary.currentSolve.has_value()); + CHECK(auxiliary.attempts == 1); + CHECK(auxiliary.kept == 0); + CHECK(auxiliary.fallbacks == 1); + CHECK(auxiliary.currentSolve->linearSolve.control.relativeTolerance == correction.correction.control.relativeTolerance); + CHECK(auxiliary.currentSolve->linearSolve.control.absoluteTolerance == correction.correction.control.absoluteTolerance); + CHECK(report.lastAcceptedCheckpointView().valid()); + }; + check(nonlinear::MonolithicGravity{}); + REQUIRE(context.isReady()); + checkPreserved(original, snapshot(context)); + check(nonlinear::GravityCurvatureCorrection{permissiveGravity()}); + REQUIRE(context.isReady()); + checkPreserved(original, snapshot(context)); }