diff --git a/experiments/PRECONDITIONER_CACHED_CROSS_ACTION_2026-09-18.md b/experiments/PRECONDITIONER_CACHED_CROSS_ACTION_2026-09-18.md new file mode 100644 index 0000000..bc4850e --- /dev/null +++ b/experiments/PRECONDITIONER_CACHED_CROSS_ACTION_2026-09-18.md @@ -0,0 +1,40 @@ +# Prepared gravity-to-material preconditioner action + +The stellar-structure gravity-to-material coupling used +`PreparedDisplacementResidualOperator::ApplyGravityGradientJacobianAction`, which +delegates to the standalone gravity-gradient force column. That column evaluates +the legacy geometry kernel during each application. The full stellar Jacobian +instead uses the already prepared complete gravity-force action, including its +cached tensor path. + +The cross preconditioner now calls that same prepared gravity-force action with +zero density and displacement variations. These zero vectors are allocated and +initialized once with the cross operator. For a fixed prepared state, + +\[ + D F_{\mathrm{gravity\ force}}[0,0,\delta q] + = (D_q F_{\mathrm{gravity\ force}})\delta q. +\] + +The existing domain-deformation adjoint then maps this mechanical action into the +shape row. Gravity-potential coupling into hydrostatic equilibrium and pressure +boundary-row replacement remain unchanged. Pressure and rotational complete +actions are not invoked: they have no gravity-gradient column. + +This is a preconditioner execution-path change. It changes no equation, +quadrature selection, cross-block approximation, factorization policy, default, +or full Newton-Jacobian implementation. Prepared and legacy force evaluations +may differ by floating-point evaluation order; numerical equality is tested with +the existing `2e-12` scaled norm criterion. + +The material-to-gravity cross action already uses the same prepared gravity +Jacobian as the full stellar operator and is unchanged. This patch does not +optimize other material/surface policies or standalone partial-column APIs. + +The existing stellar-structure FEM integration test now compares three dense +gravity directions against both the legacy partial-column calculation and the +restricted full Jacobian, checking surface and enthalpy rows separately. It +repeats after density and shape changes, checks repeated application after full-J +scratch use, and verifies that these applications do not reprepare the force. +The tests are supplied for the coordinated build; no timing improvement is +claimed before measurement on the frozen benchmark state. diff --git a/libmeanfield/interface/preconditioning/stellar_structure.cppm b/libmeanfield/interface/preconditioning/stellar_structure.cppm index 3c39294..f631f81 100644 --- a/libmeanfield/interface/preconditioning/stellar_structure.cppm +++ b/libmeanfield/interface/preconditioning/stellar_structure.cppm @@ -413,6 +413,8 @@ export namespace mean_field::preconditioning { m_gravityDirection(operation.GetGravityJacobianOperator().Width()), m_volumeDisplacement(operation.GetDomainDeformation().volumeDisplacementSize()), m_mechanicalAction(operation.GetDomainDeformation().volumeDisplacementSize()), + m_zeroDensity(operation.GetGravityContext().GetDensityMap().reduced_size()), + m_zeroDisplacement(operation.GetGravityContext().GetDisplacementMap().reduced_size()), m_zeroEnthalpy(operation.GetBarotropicClosureOperator().GetEnthalpySize()) { const auto &context = operation.GetGravityContext(); m_materialOffsets[0] = 0; @@ -425,6 +427,8 @@ export namespace mean_field::preconditioning { m_combinedOffsets[0] = 0; m_combinedOffsets[1] = MaterialSize(); m_combinedOffsets[2] = Height(); + m_zeroDensity = 0.0; + m_zeroDisplacement = 0.0; m_zeroEnthalpy = 0.0; } @@ -488,8 +492,12 @@ export namespace mean_field::preconditioning { auto enthalpyAction = MaterialBlock(materialAction, 2); densityAction = 0.0; - m_operation->GetDisplacementOperator().ApplyGravityGradientJacobianAction( - gravityGradientDirection, m_mechanicalAction + // Restrict the same prepared tensor action used by the full Jacobian. + // The standalone gravity-gradient column still uses a legacy kernel + // that rebuilds geometry during every application. Zero variations + // select this column without changing its quadrature or derivative. + m_operation->GetDisplacementOperator().GetGravityOperator().ApplyCompleteJacobianAction( + m_zeroDensity, m_zeroDisplacement, gravityGradientDirection, m_mechanicalAction ); m_operation->GetDomainDeformation().applyJacobianTranspose( m_operation->GetSurfaceDeformationParameters(), m_mechanicalAction, surfaceAction @@ -615,6 +623,8 @@ export namespace mean_field::preconditioning { mutable mfem::Vector m_gravityDirection; mutable mfem::Vector m_volumeDisplacement; mutable mfem::Vector m_mechanicalAction; + mfem::Vector m_zeroDensity; + mfem::Vector m_zeroDisplacement; mfem::Vector m_zeroEnthalpy; }; diff --git a/tests/preconditioning/stellar_structure.cpp b/tests/preconditioning/stellar_structure.cpp index 117e873..b9506ff 100644 --- a/tests/preconditioning/stellar_structure.cpp +++ b/tests/preconditioning/stellar_structure.cpp @@ -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); }