From 237485de5102d70ea72e82720456c7decfabb553 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Thu, 17 Sep 2026 09:39:39 -0400 Subject: [PATCH] Expose cached gravity target for prepared trial verification --- .../interface/solver/gravity_subsystem.cppm | 55 ++++++++++++------- tests/solver/gravity_subsystem.cpp | 9 ++- 2 files changed, 42 insertions(+), 22 deletions(-) diff --git a/libmeanfield/interface/solver/gravity_subsystem.cppm b/libmeanfield/interface/solver/gravity_subsystem.cppm index da2f77a..857979c 100644 --- a/libmeanfield/interface/solver/gravity_subsystem.cppm +++ b/libmeanfield/interface/solver/gravity_subsystem.cppm @@ -561,6 +561,23 @@ export namespace mean_field::solver { return m_baseGeneration; } + void BuildLinearPrediction( + const double alpha, + mfem::Vector &gravityResidual + ) const { + RequireCollectively( + m_hasBase && gravityResidual.Size() == m_baseResidual.Size(), + "A gravity prediction requires a captured base and a preallocated gravity residual." + ); + RequireStepLength(alpha); + bool finite = true; + for (int index = 0; index < gravityResidual.Size(); ++index) { + gravityResidual(index) = m_baseResidual(index) + alpha * m_baseAction(index); + finite = finite && std::isfinite(gravityResidual(index)); + } + RequireCollectively(finite, "The captured gravity linear prediction is nonfinite."); + } + // The caller owns candidate preparation and commits only after a new // complete residual evaluation. A failed inner solve never modifies // an accepted state; its output remains a diagnostic correction only. @@ -572,29 +589,12 @@ export namespace mean_field::solver { ) { const auto &partition = m_subsystem->GetPartition(); RequireCollectively( - m_hasBase && m_subsystem->IsCurrent() && std::isfinite(alpha) && alpha >= 0.0 && - trialResidual.Size() == partition.RootResidualSize() && + m_hasBase && m_subsystem->IsCurrent() && trialResidual.Size() == partition.RootResidualSize() && gravityCorrection.Size() == partition.GravityValueSize(), "Gravity curvature correction requires a captured base, a current " "candidate, and valid vectors." ); - double minimumAlpha = 0.0; - double maximumAlpha = 0.0; - if (MPI_Allreduce(&alpha, &minimumAlpha, 1, MPI_DOUBLE, MPI_MIN, m_subsystem->GetCommunicator()) != - MPI_SUCCESS || - MPI_Allreduce(&alpha, &maximumAlpha, 1, MPI_DOUBLE, MPI_MAX, m_subsystem->GetCommunicator()) != - MPI_SUCCESS) { - throw std::runtime_error( - "Gravity curvature correction could not verify " - "its distributed step length." - ); - } - if (minimumAlpha != maximumAlpha) { - throw std::invalid_argument( - "Gravity curvature correction requires the " - "same step length on every rank." - ); - } + RequireStepLength(alpha); partition.GatherGravityResiduals(trialResidual, m_rightHandSide); for (int index = 0; index < m_rightHandSide.Size(); ++index) { m_rightHandSide(index) = m_baseResidual(index) + alpha * m_baseAction(index) - m_rightHandSide(index); @@ -604,6 +604,23 @@ export namespace mean_field::solver { } private: + void RequireStepLength(const double alpha) const { + RequireCollectively( + std::isfinite(alpha) && alpha >= 0.0, "Gravity curvature requires a finite, non-negative step length." + ); + double minimumAlpha = 0.0; + double maximumAlpha = 0.0; + if (MPI_Allreduce(&alpha, &minimumAlpha, 1, MPI_DOUBLE, MPI_MIN, m_subsystem->GetCommunicator()) != + MPI_SUCCESS || + MPI_Allreduce(&alpha, &maximumAlpha, 1, MPI_DOUBLE, MPI_MAX, m_subsystem->GetCommunicator()) != + MPI_SUCCESS) { + throw std::runtime_error("Gravity curvature could not verify its distributed step length."); + } + if (minimumAlpha != maximumAlpha) { + throw std::invalid_argument("Gravity curvature requires the same step length on every rank."); + } + } + void RequireCollectively( const bool valid, const char *message diff --git a/tests/solver/gravity_subsystem.cpp b/tests/solver/gravity_subsystem.cpp index 2fc3fbe..d56b99b 100644 --- a/tests/solver/gravity_subsystem.cpp +++ b/tests/solver/gravity_subsystem.cpp @@ -341,7 +341,7 @@ TEST_CASE( operation, inverse, solver::GravityBlockPartition(2, 2, {1}, {0}), MPI_COMM_WORLD ); solver::PreparedGravityCurvatureCorrection correction(subsystem); - mfem::Vector baseResidual(2), direction(2), trialResidual(2), gravityStep(1); + mfem::Vector baseResidual(2), direction(2), trialResidual(2), gravityStep(1), prediction(1); baseResidual(0) = 0.75; baseResidual(1) = 0.75; direction(0) = 0.125; @@ -352,8 +352,10 @@ TEST_CASE( operation.Prepare(jacobian); inverse.Refresh(); subsystem.Refresh(); - trialResidual(0) = 0.75 + alpha * 1.75; - trialResidual(1) = 0.75 + alpha * 0.625; + trialResidual(0) = 0.75 + alpha * 1.75; + trialResidual(1) = 0.75 + alpha * 0.625; + correction.BuildLinearPrediction(alpha, prediction); + CHECK(prediction(0) == trialResidual(0)); const auto report = correction.SolveCorrection(alpha, trialResidual, gravityStep, {}); REQUIRE(report.Converged()); CHECK(report.rightHandSideNorm == 0.0); @@ -362,6 +364,7 @@ TEST_CASE( } correction.Invalidate(); CHECK_FALSE(correction.HasBase()); + CHECK_THROWS_AS(correction.BuildLinearPrediction(1.0, prediction), std::logic_error); CHECK_THROWS_AS(correction.SolveCorrection(1.0, trialResidual, gravityStep, {}), std::logic_error); }