tboudreaux/MeanField
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Add fixed-work mass-consistent gravity Schur refinement

Browse Source
This commit is contained in:
Emily Boudreaux
2026-09-17 16:19:22 -04:00
parent 018fc1c992
commit b18f41f88f
3 changed files with 531 additions and 7 deletions
Download Patch File Download Diff File Expand all files Collapse all files

298
tests/preconditioning/gravity_field.cpp
View File

@@ -1,6 +1,8 @@
#include <algorithm>
#include <array>
#include <bit>
#include <cmath>
#include <cstdint>
#include <limits>
#include <stdexcept>
#include <type_traits>
@@ -26,6 +28,10 @@ namespace {
GravityFieldBlock<backend::MatrixFreeChebyshev, FixedAMG, preconditioning::GravityApproximateLDU>;
using AdaptiveGravityLDU =
preconditioning::GravityFieldBlock<backend::Diagonal, AdaptiveAMG, preconditioning::GravityApproximateLDU>;
using ConsistentGravityLDU = preconditioning::
GravityFieldBlock<backend::MatrixFreeChebyshev, FixedAMG, preconditioning::GravityConsistentSchurLDU>;
using AdaptiveConsistentGravityLDU = preconditioning::
GravityFieldBlock<backend::MatrixFreeChebyshev, AdaptiveAMG, preconditioning::GravityConsistentSchurLDU>;
using DensityIdentity =
preconditioning::IdentityBlock<blocks::density::mass::value, blocks::density::mass::residual>;
@@ -141,6 +147,72 @@ namespace {
}
}
struct DenseSchurRefinementFixture final {
static mfem::DenseMatrix makeMass() {
mfem::DenseMatrix result(3);
result = 0.0;
result(0, 0) = 4.0;
result(1, 1) = 3.0;
result(2, 2) = 2.0;
result(0, 1) = result(1, 0) = 1.0;
result(1, 2) = result(2, 1) = 1.0;
return result;
}
static mfem::DenseMatrix makeDivergence() {
mfem::DenseMatrix result(2, 3);
result(0, 0) = 1.0;
result(0, 1) = 2.0;
result(0, 2) = -1.0;
result(1, 0) = 2.0;
result(1, 1) = -1.0;
result(1, 2) = 1.0;
return result;
}
static mfem::DenseMatrix makeSurrogate() {
// B diag(M)^-1 B^T differs from B M^-1 B^T. Its inverse
// does not commute with the latter, so a reversed correction
// order cannot pass the rational-reference checks below.
mfem::DenseMatrix result(2);
result(0, 0) = 25.0 / 12.0;
result(0, 1) = result(1, 0) = -2.0 / 3.0;
result(1, 1) = 11.0 / 6.0;
return result;
}
mfem::DenseMatrix mass{makeMass()};
mfem::DenseMatrix divergence{makeDivergence()};
mfem::DenseMatrix surrogate{makeSurrogate()};
backend::PreparedDenseDirect massInverse{backend::DenseDirect{}, mass};
backend::PreparedDenseDirect schurInverse{backend::DenseDirect{}, surrogate};
static mfem::Vector rightHandSide() {
mfem::Vector result(5);
result(0) = 3.0;
result(1) = -2.0;
result(2) = 1.0;
result(3) = 4.0;
result(4) = -3.0;
return result;
}
mfem::Vector residual(const mfem::Vector &action, const mfem::Vector &rhs) const {
mfem::Vector result(5);
mfem::Vector flux(const_cast<mfem::real_t *>(action.GetData()), 3);
mfem::Vector potential(const_cast<mfem::real_t *>(action.GetData()) + 3, 2);
mfem::Vector gradientResidual(result.GetData(), 3);
mfem::Vector potentialResidual(result.GetData() + 3, 2);
mfem::Vector transposeAction(3);
mass.Mult(flux, gradientResidual);
divergence.MultTranspose(potential, transposeAction);
gradientResidual += transposeAction;
divergence.Mult(flux, potentialResidual);
result -= rhs;
return result;
}
};
struct PreparedGeometry final {
mean_field::fem::FEM finiteElements;
gravity_context::GravityFieldGeometryContext context;
@@ -175,9 +247,15 @@ TEST_CASE(
STATIC_CHECK(preconditioning::PreconditionerComponent<FixedGravityLDU>);
STATIC_CHECK(preconditioning::PreconditionerComponent<ChebyshevGravityLDU>);
STATIC_CHECK(preconditioning::PreconditionerComponent<AdaptiveGravityLDU>);
STATIC_CHECK(preconditioning::PreconditionerComponent<ConsistentGravityLDU>);
STATIC_CHECK(preconditioning::GravityFactorizationPolicy<preconditioning::GravityConsistentSchurLDU>);
STATIC_CHECK(preconditioning::backend::ArnoldiAdmissible<typename FixedGravityLDU::BackendType>);
STATIC_CHECK(preconditioning::backend::ArnoldiAdmissible<typename ChebyshevGravityLDU::BackendType>);
STATIC_CHECK_FALSE(preconditioning::backend::ArnoldiAdmissible<typename AdaptiveGravityLDU::BackendType>);
STATIC_CHECK(preconditioning::backend::ArnoldiAdmissible<typename ConsistentGravityLDU::BackendType>);
STATIC_CHECK_FALSE(
preconditioning::backend::ArnoldiAdmissible<typename AdaptiveConsistentGravityLDU::BackendType>
);
STATIC_CHECK(preconditioning::CompletePreconditionerFor<FixedGravityPlan, Form>);
STATIC_CHECK(preconditioning::CompatiblePreconditionerFor<FixedGravityPlan, Form, JacobianForm>);
STATIC_CHECK(preconditioning::StationaryLinearPreconditionerPlan<FixedGravityPlan>);
@@ -216,6 +294,123 @@ TEST_CASE(
tags::preconditioning_gravity_unit
) {
checkExactDenseRecovery(preconditioning::GravityApproximateLDU{});
checkExactDenseRecovery(preconditioning::GravityConsistentSchurLDU{.correctionSweeps = 2, .relaxation = 0.5});
}
TEST_CASE(
"Mass-Consistent Gravity Schur Refinement Preserves The Baseline And Signed Algebra",
tags::preconditioning_gravity_unit
) {
DenseSchurRefinementFixture fixture;
const mfem::Vector rhs = fixture.rightHandSide();
preconditioning::GravityFactorizationOperator baseline(
preconditioning::GravityApproximateLDU{}, fixture.massInverse, fixture.schurInverse, fixture.divergence
);
mfem::Vector baselineAction(5);
baseline.Mult(rhs, baselineAction);
const double baselineResidual = fixture.residual(baselineAction, rhs).Norml2();
// These rational answers were obtained from the two dense saddle-point
// blocks, with S_R=[3,-2;-2,59/18] and omega=1/2. They exercise a genuinely
// inexact Schur inverse rather than an already-exact scalar example.
const std::array<std::array<double, 5>, 2> expected{{
{-458734.0 / 4782969.0, 6011089.0 / 4782969.0, -7518059.0 / 4782969.0,
-646046.0 / 531441.0, 888176.0 / 531441.0},
{-631839955.0 / 6973568802.0, 4376327111.0 / 3486784401.0, -10946970773.0 / 6973568802.0,
-472857412.0 / 387420489.0, 644634595.0 / 387420489.0}
}};
double previousResidual = baselineResidual;
for (const int sweeps : {0, 1, 2}) {
CAPTURE(sweeps);
preconditioning::GravityFactorizationOperator refined(
preconditioning::GravityConsistentSchurLDU{.correctionSweeps = sweeps, .relaxation = 0.5},
fixture.massInverse, fixture.schurInverse, fixture.divergence
);
mfem::Vector action(5);
action = std::numeric_limits<double>::quiet_NaN();
const auto *const actionStorage = action.GetData();
refined.Mult(rhs, action);
CHECK(action.GetData() == actionStorage);
for (int index = 0; index < action.Size(); ++index) {
if (sweeps == 0) {
CHECK(std::bit_cast<std::uint64_t>(action(index)) ==
std::bit_cast<std::uint64_t>(baselineAction(index)));
} else {
CHECK(action(index) == Catch::Approx(expected[sweeps - 1][index]).margin(2.0e-13));
}
}
const auto &statistics = refined.GetStatistics();
CHECK(statistics.applications == 1);
CHECK(statistics.massInverseApplications == static_cast<std::uint64_t>(2 + sweeps));
CHECK(statistics.potentialSchurApplications == static_cast<std::uint64_t>(1 + sweeps));
CHECK(statistics.divergenceApplications == static_cast<std::uint64_t>(1 + sweeps));
CHECK(statistics.transposeDivergenceApplications == static_cast<std::uint64_t>(1 + sweeps));
CHECK(statistics.schurRefinementSweeps == static_cast<std::uint64_t>(sweeps));
const mfem::Vector residual = fixture.residual(action, rhs);
for (int index = 0; index < 3; ++index) CHECK(std::abs(residual(index)) < 2.0e-13);
if (sweeps > 0) {
CHECK(residual.Norml2() < previousResidual);
previousResidual = residual.Norml2();
}
}
}
TEST_CASE(
"Fixed Gravity Schur Refinement Is Linear And Has No Application History",
tags::preconditioning_gravity_unit
) {
DenseSchurRefinementFixture fixture;
preconditioning::GravityFactorizationOperator refined(
preconditioning::GravityConsistentSchurLDU{.correctionSweeps = 2, .relaxation = 0.5},
fixture.massInverse, fixture.schurInverse, fixture.divergence
);
const mfem::Vector x = fixture.rightHandSide();
mfem::Vector y(5), combination(x), px(5), py(5), actual(5), expected(5), repeat(5);
for (int index = 0; index < y.Size(); ++index) y(index) = std::sin(0.7 * (index + 1));
combination *= -2.3;
combination.Add(0.4, y);
refined.Mult(x, px);
refined.Mult(y, py);
refined.Mult(combination, actual);
expected = px;
expected *= -2.3;
expected.Add(0.4, py);
actual -= expected;
CHECK(actual.Norml2() < 2.0e-13 * expected.Norml2());
repeat = std::numeric_limits<double>::quiet_NaN();
refined.Mult(x, repeat);
for (int index = 0; index < repeat.Size(); ++index) CHECK(repeat(index) == px(index));
combination = 0.0;
refined.Mult(combination, actual);
CHECK(actual.Norml2() == 0.0);
CHECK(refined.GetStatistics().applications == 5);
CHECK(refined.GetStatistics().schurRefinementSweeps == 10);
}
TEST_CASE(
"Gravity Schur Refinement Rejects Invalid Fixed-Work Options",
tags::preconditioning_gravity_unit
) {
DenseSchurRefinementFixture fixture;
const std::array<preconditioning::GravityConsistentSchurLDU, 5> invalid{{
{.correctionSweeps = -1, .relaxation = 1.0},
{.correctionSweeps = 0, .relaxation = 0.0},
{.correctionSweeps = 1, .relaxation = -0.5},
{.correctionSweeps = 1, .relaxation = std::numeric_limits<double>::infinity()},
{.correctionSweeps = 1, .relaxation = std::numeric_limits<double>::quiet_NaN()}
}};
for (const auto policy : invalid) {
CAPTURE(policy.correctionSweeps, policy.relaxation);
CHECK_THROWS_AS(
(preconditioning::GravityFactorizationOperator(
policy, fixture.massInverse, fixture.schurInverse, fixture.divergence
)),
std::invalid_argument
);
}
}
TEST_CASE(
@@ -398,6 +593,109 @@ TEST_CASE(
}
}
TEST_CASE(
"Prepared Gravity Schur Refinement Reuses Sparse Couplings Across Geometry Revisions",
tags::preconditioning_gravity_integration
) {
const auto arguments = test_utils::setup_args();
PreparedGeometry geometry(arguments);
const auto makeBlock = [](const int sweeps) {
return preconditioning::GravityFieldBlock(
backend::MatrixFreeChebyshev{.order = 2, .powerIterations = 10},
FixedAMG{backend::FixedCycles{.cycles = 1}},
preconditioning::GravityConsistentSchurLDU{.correctionSweeps = sweeps, .relaxation = 0.5}
);
};
auto baseline = preconditioning::prepare(geometry.finiteElements, geometry.context, makeBlock(0));
auto refined = preconditioning::prepare(geometry.finiteElements, geometry.context, makeBlock(2));
const MPI_Comm communicator = geometry.finiteElements.mesh->GetComm();
const auto gradientMap = geometry.context.GetMassOperator().GetFluxMap();
const auto potentialMap = geometry.context.GetSourceOperator().GetPotentialMap();
const mfem::Vector rhs = gravity_prepared_test_utils::make_deterministic_vector(refined.Width(), 0.29);
mfem::Vector actual(refined.Height()), expected(refined.Height());
const auto compareApplications = [&] {
const auto sparse = preconditioning::assembleGravityDivergenceSurrogate(geometry.finiteElements);
preconditioning::ReducedGravityDivergenceOperator sparseReduced(*sparse, gradientMap, potentialMap);
preconditioning::ReducedGravityDivergenceOperator matrixFreeReduced(
geometry.context.GetDivergenceOperator(), gradientMap, potentialMap
);
for (const double phase : {0.19, 0.53, 0.87}) {
const auto flux = gravity_prepared_test_utils::make_deterministic_vector(sparseReduced.Width(), phase);
const auto potential =
gravity_prepared_test_utils::make_deterministic_vector(sparseReduced.Height(), phase + 0.17);
mfem::Vector sparseAction(sparseReduced.Height()), matrixFreeAction(matrixFreeReduced.Height());
mfem::Vector sparseTranspose(sparseReduced.Width()), matrixFreeTranspose(matrixFreeReduced.Width());
sparseReduced.Mult(flux, sparseAction);
matrixFreeReduced.Mult(flux, matrixFreeAction);
sparseReduced.MultTranspose(potential, sparseTranspose);
matrixFreeReduced.MultTranspose(potential, matrixFreeTranspose);
CHECK(gravity_prepared_test_utils::relative_error(sparseAction, matrixFreeAction, communicator) < 2.0e-12);
CHECK(gravity_prepared_test_utils::relative_error(sparseTranspose, matrixFreeTranspose, communicator) <
2.0e-12);
}
// Reuse precisely the same child solvers, so zero-sweep equivalence
// checks the path itself rather than independent AMG setup roundoff.
preconditioning::ReducedFieldSolverAdapter baselineSchur(baseline.GetPotentialSchurInverse(), potentialMap);
preconditioning::GravityFactorizationOperator baselineReference(
preconditioning::GravityApproximateLDU{}, baseline.GetMassInverse(), baselineSchur, matrixFreeReduced
);
baseline.Mult(rhs, actual);
baselineReference.Mult(rhs, expected);
for (int index = 0; index < actual.Size(); ++index) {
CHECK(std::bit_cast<std::uint64_t>(actual(index)) == std::bit_cast<std::uint64_t>(expected(index)));
}
preconditioning::ReducedFieldSolverAdapter refinedSchur(refined.GetPotentialSchurInverse(), potentialMap);
preconditioning::GravityFactorizationOperator refinedReference(
makeBlock(2).factorizationPolicy(), refined.GetMassInverse(), refinedSchur, matrixFreeReduced
);
refined.Mult(rhs, actual);
refinedReference.Mult(rhs, expected);
CHECK(gravity_prepared_test_utils::relative_error(actual, expected, communicator) < 2.0e-10);
};
compareApplications();
CHECK(refined.GetDivergenceAssemblyCount() == 1);
CHECK_FALSE(refined.Refresh(geometry.finiteElements, geometry.context).DidAnyWork());
const mfem::Vector displacementTrue = gravity_prepared_test_utils::make_displacement(geometry.finiteElements, 0.4);
const mfem::Vector displacement = geometry.context.GetDisplacementMap().gather(displacementTrue);
geometry.context.PreparePrimal(displacement, {.value = 1}, {.value = 2});
CHECK_THROWS_AS(refined.Mult(rhs, actual), std::logic_error);
const auto geometryRefresh = refined.Refresh(geometry.finiteElements, geometry.context);
(void) baseline.Refresh(geometry.finiteElements, geometry.context);
CHECK(geometryRefresh.geometryChanged);
CHECK_FALSE(geometryRefresh.rebuiltDivergenceBinding);
CHECK(refined.GetDivergenceAssemblyCount() == 1);
compareApplications();
geometry.context.PreparePrimal(displacement, {.value = 2}, {.value = 2});
CHECK_THROWS_AS(refined.Mult(rhs, actual), std::logic_error);
const auto discretizationRefresh = refined.Refresh(geometry.finiteElements, geometry.context);
(void) baseline.Refresh(geometry.finiteElements, geometry.context);
CHECK(discretizationRefresh.rebuiltDivergenceBinding);
CHECK(refined.GetDivergenceAssemblyCount() == 2);
compareApplications();
int rank = 0, rankCount = 0;
MPI_Comm_rank(communicator, &rank);
MPI_Comm_size(communicator, &rankCount);
if (rankCount > 1) {
const auto checkCollectiveRejection = [&](const preconditioning::GravityConsistentSchurLDU policy) {
const auto block = preconditioning::GravityFieldBlock(
backend::Diagonal{}, FixedAMG{backend::FixedCycles{.cycles = 1}}, policy
);
CHECK_THROWS_AS(
preconditioning::prepare(geometry.finiteElements, geometry.context, block), std::invalid_argument
);
};
checkCollectiveRejection({.correctionSweeps = rank == 0 ? -1 : 1, .relaxation = 0.5});
checkCollectiveRejection({.correctionSweeps = rank == 0 ? 0 : 1, .relaxation = 0.5});
checkCollectiveRejection({.correctionSweeps = 1, .relaxation = rank == 0 ? 0.5 : 1.0});
}
}
TEST_CASE(
"Fixed Gravity Divergence Cache Preserves Changing Mass Weights And Invalidates",
tags::preconditioning_gravity_integration