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

Reuse prepared gravity force for stellar cross preconditioning

Browse Source
This commit is contained in:
Emily Boudreaux
2026-09-17 16:28:56 -04:00
parent 7f7d7b1cfa
commit fac147774a
3 changed files with 131 additions and 2 deletions
Download Patch File Download Diff File Expand all files Collapse all files

79
tests/preconditioning/stellar_structure.cpp
View File

@@ -413,6 +413,70 @@ TEST_CASE(
CHECK(relativeError(crossGravityPotential, expectedGravityPotential) <= 2.0e-12);
CHECK(relativeError(crossAction, expected) <= 2.0e-12);
const auto checkPreparedGravityColumn = [&](const int revision) {
CAPTURE(revision);
const auto &mechanics = physical.GetDisplacementOperator();
const auto preparationCount = mechanics.GetGravityOperator().GetResidualPreparationCount();
mfem::Vector gravityInput(cross.GravitySize()), actual(cross.MaterialSize()), legacy(cross.MaterialSize());
mfem::Vector repeated(cross.MaterialSize()), rootInput(physical.Width()), rootAction(physical.Height());
mfem::Vector mechanicalReference(physical.GetDomainDeformation().volumeDisplacementSize());
mfem::Vector zeroEnthalpy(physical.GetBarotropicClosureOperator().GetEnthalpySize());
zeroEnthalpy = 0.0;
for (const double phase : {0.19, 0.53, 0.87}) {
CAPTURE(phase);
for (int index = 0; index < gravityInput.Size(); ++index) {
gravityInput(index) = 0.01 * std::sin(phase * (index + 1));
}
const mfem::Vector gradientInput(gravityInput.GetData(), gravityOffsets[1]);
const mfem::Vector potentialInput(
gravityInput.GetData() + gravityOffsets[1], gravityOffsets[2] - gravityOffsets[1]
);
legacy = 0.0;
mfem::Vector legacySurface(legacy, materialOffsets[1], materialOffsets[2] - materialOffsets[1]);
mfem::Vector legacyEnthalpy(legacy, materialOffsets[2], materialOffsets[3] - materialOffsets[2]);
// Keep the independent legacy column as the numerical reference;
// the production cross action now uses the cached complete kernel.
mechanics.ApplyGravityGradientJacobianAction(gradientInput, mechanicalReference);
physical.GetDomainDeformation().applyJacobianTranspose(
physical.GetSurfaceDeformationParameters(), mechanicalReference, legacySurface
);
physical.GetHydrostaticOperator().ApplyGravityPotentialJacobianAction(potentialInput, legacyEnthalpy);
physical.GetSurfaceConstraintOperator().ApplyJacobianRows(zeroEnthalpy, legacyEnthalpy);
legacySurface.SyncAliasMemory(legacy);
legacyEnthalpy.SyncAliasMemory(legacy);
cross.ApplyGravityToMaterial(gravityInput, actual);
const mfem::Vector actualSurface(
actual.GetData() + materialOffsets[1], materialOffsets[2] - materialOffsets[1]
);
const mfem::Vector actualEnthalpy(
actual.GetData() + materialOffsets[2], materialOffsets[3] - materialOffsets[2]
);
INFO("legacy/cached surface error = " << relativeError(actualSurface, legacySurface)
<< ", cached norm = " << actualSurface.Norml2() << ", legacy norm = " << legacySurface.Norml2());
CHECK(relativeError(actualSurface, legacySurface) <= 2.0e-12);
CHECK(relativeError(actualEnthalpy, legacyEnthalpy) <= 2.0e-12);
CHECK(relativeError(actual, legacy) <= 2.0e-12);
rootInput = 0.0;
const auto inputView = physical.GetRootManifest().stateView(rootInput);
mfem::Vector rootGradient = inputView.block(blocks::gravity_field.gradient_term);
mfem::Vector rootPotential = inputView.block(blocks::gravity_field.poisson_term);
rootGradient = gradientInput;
rootPotential = potentialInput;
rootGradient.SyncAliasMemory(rootInput);
rootPotential.SyncAliasMemory(rootInput);
physical.Mult(rootInput, rootAction);
const auto actionView = physical.GetRootManifest().residualView(rootAction);
CHECK(relativeError(actualSurface,
actionView.block(blocks::surface_deformation_field.shape_equilibrium_term)) <= 2.0e-12);
CHECK(relativeError(actualEnthalpy, actionView.block(blocks::enthalpy_field.specific_term)) <= 2.0e-12);
cross.ApplyGravityToMaterial(gravityInput, repeated);
CHECK(relativeError(actual, repeated) <= 2.0e-14);
}
CHECK(mechanics.GetGravityOperator().GetResidualPreparationCount() == preparationCount);
};
checkPreparedGravityColumn(1);
auto materialBlock = preconditioning::materialSurfaceBlock(problem);
auto gravityBlock = preconditioning::GravityFieldBlock(
backend::Diagonal{}, FixedAMG{backend::FixedCycles{.cycles = 1}}, preconditioning::GravityApproximateLDU{}
@@ -440,4 +504,19 @@ TEST_CASE(
const auto unchanged = prepared.Refresh();
CHECK_FALSE(unchanged.DidAnyWork());
CHECK(prepared.IsCurrent());
mfem::Vector changedState(projected.values);
const auto changedView = problem.GetManifest().stateView(changedState);
mfem::Vector changedDensity = changedView.block(blocks::density_field.mass_term);
mfem::Vector changedSurface = changedView.block(blocks::surface_deformation_field.parameters_term);
changedDensity *= 1.001;
for (int index = 0; index < changedSurface.Size(); ++index) {
changedSurface(index) += 1.0e-7 * std::sin(0.31 * (index + 1));
}
changedDensity.SyncAliasMemory(changedState);
changedSurface.SyncAliasMemory(changedState);
problem.Prepare(changedState, makeDependencies(2), zeroRotation());
// The cross operator borrows the owning physical operator, whose prepared
// geometry and density are replaced in place. It must not retain old data.
checkPreparedGravityColumn(2);
}