diff --git a/experiments/GRAVITY_SUBSYSTEM_IMPLEMENTATION_2026-09-17.md b/experiments/GRAVITY_SUBSYSTEM_IMPLEMENTATION_2026-09-17.md new file mode 100644 index 0000000..7bdc769 --- /dev/null +++ b/experiments/GRAVITY_SUBSYSTEM_IMPLEMENTATION_2026-09-17.md @@ -0,0 +1,85 @@ +# Prepared gravity subsystem and curvature correction + +This branch introduces solver components; the complete equilibrium equations, +normalization, and model specification remain the source of truth. + +## Block and lifetime contracts + +`MakeGravityBlockPartition
(layout)` selects the gravity-gradient and +gravity-potential blocks by their compiled value and residual types. Its +constraints reject forms lacking those blocks. The complement is retained in +root order, including any generated angular-momentum or central-density border. +There are no string identifiers or fixed root offsets in this selection. + +`GravityBlockPartition` also has a checked explicit-index constructor for +algebraic tests and other compiled adapters. Local partitions may be empty; +gravity and retained local systems must be square. A globally empty gravity +system is rejected by the prepared linear backend. Gather/injection require +preallocated, nonoverlapping storage. + +`PreparedGravitySubsystem` owns restricted operator/inverse workspaces and the +existing prepared linear backend. It borrows the full normalized operator, +inverse, and communicator at stable addresses. It is neither movable nor +copyable. Preparation records the root preparation generation. After candidate +preparation the caller must refresh its physical inverse and explicitly refresh +the subsystem; stale applications and solves fail. The subsystem does not +prepare a candidate, change an accepted state, or commit a step. + +The exact gravity action is the principal restriction of the root Jacobian. +The preconditioner is the corresponding restriction of the supplied root +inverse. This reuses existing normalization, gravity boundary conditions, and +preconditioning recipes, but does not assert that the restricted inverse is an +exact gravity inverse. Every solve retains the backend's independently checked +true-defect gate. Both wrappers preallocate full-vector workspaces; applications +perform gather/scatter and the existing root actions without allocating their +own vectors. Full root actions may do unnecessary non-gravity work; a specialized +gravity action/inverse adapter remains a future optimization requiring equality +tests against this reference restriction. + +Attempted solves, completed inner iterations, failed calls, completed solve time, +and the last completed solve report are available without history allocation. +Contract exceptions count as failed calls but do not trigger new collectives to +measure elapsed time during unwinding. + +The existing FGMRES backend now admits ranks with zero local equations while +requiring a nonempty global system. `GetConfiguration()` returns its immutable +prescription so another prepared operator can reuse the same backend options. + +## Curvature correction + +`PreparedGravityCurvatureCorrection` captures `G(x)` and `J_G(x) p` before any +candidate preparation. For a prepared candidate `x + alpha p` it solves + +``` +G_z(x + alpha p) delta_z = G(x) + alpha J_G(x) p - G(x + alpha p). +``` + +Only the compact gravity correction is returned. The caller may inject it only +after a converged inner report, then must prepare and evaluate the complete +candidate again before applying its usual acceptance rules. In particular, +improving gravity does not certify the hydrostatic, material, or shape residual. + +For smooth equations with uniformly bounded inverse gravity block, the right +side and correction are second order in `alpha`. This preserves the base and +Newton tangent even when the base does not lie on the gravity solution manifold. +It differs from projecting the candidate to `G=0`, which can introduce a finite +correction at zero step when the base gravity residual is nonzero. Inexact inner +tolerances control cancellation of the remainder; they cannot replace the outer +equilibrium and virial checks. + +The owner captures a new base for each new Newton direction and invalidates the +cache when that direction is discarded. The stored base generation is exposed +for diagnostics; preparing trials intentionally does not invalidate the captured +base arrays. Reusing those arrays after replacing the base or direction violates +the caller's explicit lifecycle contract. + +## Test coverage + +`tests/solver/gravity_subsystem.cpp` covers typed forms with/without a central +border, noncontiguous and differently ordered row/value maps, invalid maps and +overlapping storage, exact restricted solves, true-defect failure, nonfinite +inputs, stale operator/inverse preparation, locally empty MPI partitions, +globally empty rejection, preserved backend options, zero affine curvature at +an off-manifold base, and the second-order correction at curved trial geometry. +The root integration campaign compiles these components and tests real stellar +states; this branch does not claim a new validated stellar equilibrium. diff --git a/libmeanfield/interface/solver/gravity_subsystem.cppm b/libmeanfield/interface/solver/gravity_subsystem.cppm new file mode 100644 index 0000000..da2f77a --- /dev/null +++ b/libmeanfield/interface/solver/gravity_subsystem.cppm @@ -0,0 +1,632 @@ +module; + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +export module mean_field:solver.gravity_subsystem; + +export import :solver.linear_backend; +export import :utils.blocks; + +export namespace mean_field::solver { + /* The maps are local true-dof maps. Gravity and retained coordinates need + * not be contiguous or share the residual ordering. Only the discretized + * sizes are runtime data; the production factory selects block identities + * from the compiled equation form. */ + class GravityBlockPartition final { + public: + GravityBlockPartition( + const int rootValueSize, + const int rootResidualSize, + std::vector gravityValues, + std::vector gravityResiduals + ) + : m_rootValueSize(rootValueSize), + m_rootResidualSize(rootResidualSize), + m_gravityValues(std::move(gravityValues)), + m_gravityResiduals(std::move(gravityResiduals)), + m_retainedValues(Complement( + rootValueSize, + m_gravityValues + )), + m_retainedResiduals(Complement( + rootResidualSize, + m_gravityResiduals + )) { + if (m_gravityValues.size() != m_gravityResiduals.size() || + m_retainedValues.size() != m_retainedResiduals.size()) { + throw std::invalid_argument("Gravity elimination requires square gravity and retained blocks."); + } + } + + [[nodiscard]] int RootValueSize() const noexcept { + return m_rootValueSize; + } + [[nodiscard]] int RootResidualSize() const noexcept { + return m_rootResidualSize; + } + [[nodiscard]] int GravityValueSize() const noexcept { + return static_cast(m_gravityValues.size()); + } + [[nodiscard]] int GravityResidualSize() const noexcept { + return static_cast(m_gravityResiduals.size()); + } + [[nodiscard]] int RetainedValueSize() const noexcept { + return static_cast(m_retainedValues.size()); + } + [[nodiscard]] int RetainedResidualSize() const noexcept { + return static_cast(m_retainedResiduals.size()); + } + + void GatherGravityValues( + const mfem::Vector &full, + mfem::Vector &compact + ) const { + Gather(m_gravityValues, m_rootValueSize, full, compact); + } + void GatherGravityResiduals( + const mfem::Vector &full, + mfem::Vector &compact + ) const { + Gather(m_gravityResiduals, m_rootResidualSize, full, compact); + } + void GatherRetainedValues( + const mfem::Vector &full, + mfem::Vector &compact + ) const { + Gather(m_retainedValues, m_rootValueSize, full, compact); + } + void GatherRetainedResiduals( + const mfem::Vector &full, + mfem::Vector &compact + ) const { + Gather(m_retainedResiduals, m_rootResidualSize, full, compact); + } + void InjectGravityValues( + const mfem::Vector &compact, + mfem::Vector &full + ) const { + Inject(m_gravityValues, m_rootValueSize, compact, full); + } + void InjectGravityResiduals( + const mfem::Vector &compact, + mfem::Vector &full + ) const { + Inject(m_gravityResiduals, m_rootResidualSize, compact, full); + } + void InjectRetainedValues( + const mfem::Vector &compact, + mfem::Vector &full + ) const { + Inject(m_retainedValues, m_rootValueSize, compact, full); + } + void InjectRetainedResiduals( + const mfem::Vector &compact, + mfem::Vector &full + ) const { + Inject(m_retainedResiduals, m_rootResidualSize, compact, full); + } + void AddGravityValues( + const mfem::Vector &compact, + mfem::Vector &full + ) const { + VerifySizes(m_gravityValues, m_rootValueSize, full, compact); + for (std::size_t index = 0; index < m_gravityValues.size(); ++index) { + full(m_gravityValues[index]) += compact(static_cast(index)); + } + } + + private: + [[nodiscard]] static std::vector Complement( + const int size, + const std::vector &selected + ) { + if (size < 0) { + throw std::invalid_argument("A gravity partition requires non-negative root dimensions."); + } + std::vector used(static_cast(size), false); + for (const int index : selected) { + if (index < 0 || index >= size || used[static_cast(index)]) { + throw std::invalid_argument("Gravity partition indices must be unique and in range."); + } + used[static_cast(index)] = true; + } + std::vector complement; + complement.reserve(static_cast(size) - selected.size()); + for (int index = 0; index < size; ++index) { + if (!used[static_cast(index)]) { + complement.push_back(index); + } + } + return complement; + } + + static void VerifySizes( + const std::vector &indices, + const int fullSize, + const mfem::Vector &full, + const mfem::Vector &compact + ) { + if (full.Size() != fullSize || compact.Size() != static_cast(indices.size())) { + throw std::invalid_argument("Gravity partition vectors must have their preallocated sizes."); + } + const std::less before; + if (compact.Size() != 0 && full.Size() != 0 && before(full.GetData(), compact.GetData() + compact.Size()) && + before(compact.GetData(), full.GetData() + full.Size())) { + throw std::invalid_argument("Gravity gather and injection require distinct vector storage."); + } + } + + static void Gather( + const std::vector &indices, + const int fullSize, + const mfem::Vector &full, + mfem::Vector &compact + ) { + VerifySizes(indices, fullSize, full, compact); + for (std::size_t index = 0; index < indices.size(); ++index) { + compact(static_cast(index)) = full(indices[index]); + } + } + + static void Inject( + const std::vector &indices, + const int fullSize, + const mfem::Vector &compact, + mfem::Vector &full + ) { + VerifySizes(indices, fullSize, full, compact); + full = 0.0; + for (std::size_t index = 0; index < indices.size(); ++index) { + full(indices[index]) = compact(static_cast(index)); + } + } + + int m_rootValueSize; + int m_rootResidualSize; + std::vector m_gravityValues; + std::vector m_gravityResiduals; + std::vector m_retainedValues; + std::vector m_retainedResiduals; + }; + + template + concept GravityPartitionForm = requires { + typename Form::value_blocks; + typename Form::residual_blocks; + requires utils::blocks::block_form_is_valid_v; + requires utils::blocks::contains_type_v; + requires utils::blocks::contains_type_v; + requires utils::blocks::contains_type_v< + utils::blocks::gravity::gradient::residual, typename Form::residual_blocks>; + requires utils::blocks::contains_type_v< + utils::blocks::gravity::poisson::residual, typename Form::residual_blocks>; + }; + + template + [[nodiscard]] GravityBlockPartition MakeGravityBlockPartition(const utils::blocks::form_layout &layout) { + using namespace utils::blocks; + std::vector values; + std::vector residuals; + const auto append = [](std::vector &indices, const int offset, const int size) { + for (int index = 0; index < size; ++index) { + indices.push_back(offset + index); + } + }; + const auto value = [&]() { + constexpr auto block = value_block>{}; + append(values, layout.offset(block), layout.size(block)); + }; + const auto residual = [&]() { + constexpr auto block = residual_block>{}; + append(residuals, layout.offset(block), layout.size(block)); + }; + value.template operator()(); + value.template operator()(); + residual.template operator()(); + residual.template operator()(); + return { + layout.value_offsets().Last(), layout.residual_offsets().Last(), std::move(values), std::move(residuals) + }; + } + + template + concept GravitySubsystemOperation = + std::derived_from && requires(const Operation &operation) { + { operation.IsPrepared() } -> std::same_as; + { operation.GetProblem().GetPreparationGeneration() } -> std::convertible_to; + }; + + template + concept GravitySubsystemPreconditioner = + std::derived_from && + LinearPreconditionerApplicationContractAvailable && requires(const Preconditioner &inverse) { + { inverse.IsCurrent() } -> std::same_as; + }; +} // namespace mean_field::solver + +namespace mean_field::solver::detail { + template class GravityRestrictedOperator final : public mfem::Operator { + public: + GravityRestrictedOperator( + const Operation &operation, + const GravityBlockPartition &partition + ) + : mfem::Operator( + partition.GravityResidualSize(), + partition.GravityValueSize() + ), + m_operation(std::addressof(operation)), + m_partition(std::addressof(partition)), + m_fullDirection(partition.RootValueSize()), + m_fullAction(partition.RootResidualSize()) { + if (operation.Width() != partition.RootValueSize() || operation.Height() != partition.RootResidualSize()) { + throw std::invalid_argument("The gravity partition does not match its root operator."); + } + Refresh(); + } + + void Refresh() { + if (!m_operation->IsPrepared()) { + throw std::logic_error("Gravity requires a prepared root linearization before refresh."); + } + m_generation = m_operation->GetProblem().GetPreparationGeneration(); + } + + [[nodiscard]] bool IsCurrent() const { + return m_operation->IsPrepared() && m_generation == m_operation->GetProblem().GetPreparationGeneration(); + } + + void Mult( + const mfem::Vector &direction, + mfem::Vector &action + ) const override { + if (!IsCurrent()) { + throw std::logic_error("The gravity linearization must be refreshed after preparation."); + } + m_partition->InjectGravityValues(direction, m_fullDirection); + m_operation->Mult(m_fullDirection, m_fullAction); + m_partition->GatherGravityResiduals(m_fullAction, action); + } + + private: + const Operation *m_operation; + const GravityBlockPartition *m_partition; + std::uint64_t m_generation{0}; + mutable mfem::Vector m_fullDirection; + mutable mfem::Vector m_fullAction; + }; + + template + class GravityRestrictedInverse final : public mfem::Solver { + public: + static constexpr auto applicationContract = linearPreconditionerApplicationContract; + + GravityRestrictedInverse( + const Preconditioner &inverse, + const GravityRestrictedOperator &operation, + const GravityBlockPartition &partition + ) + : mfem::Solver(partition.GravityValueSize()), + m_inverse(std::addressof(inverse)), + m_operation(std::addressof(operation)), + m_partition(std::addressof(partition)), + m_fullResidual(partition.RootResidualSize()), + m_fullCorrection(partition.RootValueSize()) { + if (inverse.Width() != partition.RootResidualSize() || inverse.Height() != partition.RootValueSize()) { + throw std::invalid_argument("The gravity partition does not match its root inverse."); + } + } + + [[nodiscard]] bool IsCurrent() const { + return m_operation->IsCurrent() && m_inverse->IsCurrent(); + } + + void SetOperator(const mfem::Operator &operation) override { + if (std::addressof(operation) != m_operation) { + throw std::invalid_argument( + "A prepared gravity inverse cannot be " + "rebound to a different operator." + ); + } + } + + void Mult( + const mfem::Vector &residual, + mfem::Vector &correction + ) const override { + if (!IsCurrent()) { + throw std::logic_error("The gravity inverse requires a current root inverse."); + } + m_partition->InjectGravityResiduals(residual, m_fullResidual); + m_inverse->Mult(m_fullResidual, m_fullCorrection); + m_partition->GatherGravityValues(m_fullCorrection, correction); + } + + private: + const Preconditioner *m_inverse; + const GravityRestrictedOperator *m_operation; + const GravityBlockPartition *m_partition; + mutable mfem::Vector m_fullResidual; + mutable mfem::Vector m_fullCorrection; + }; +} // namespace mean_field::solver::detail + +export namespace mean_field::solver { + /* Restricting the supplied inverse is a preconditioner, not an assertion + * that its principal block equals G_z^{-1}. The backend still verifies the + * true gravity defect. Keeping the exact root operator as the source of + * every action also retains its boundary rows and normalization. */ + template < + GravitySubsystemOperation Operation, + GravitySubsystemPreconditioner Preconditioner, + LinearBackendConfiguration Backend = linear::FGMRES> + requires LinearBackendRuntimeAvailableFor< + Backend, + detail::GravityRestrictedOperator, + detail::GravityRestrictedInverse> + class PreparedGravitySubsystem final { + public: + using OperatorType = Operation; + using PreconditionerType = Preconditioner; + using RestrictedOperator = detail::GravityRestrictedOperator; + using RestrictedInverse = detail::GravityRestrictedInverse; + using BackendType = PreparedLinearBackendType; + + PreparedGravitySubsystem( + const Operation &operation, + const Preconditioner &preconditioner, + GravityBlockPartition partition, + const MPI_Comm communicator, + Backend backend = {} + ) + : m_fullOperator(std::addressof(operation)), + m_fullPreconditioner(std::addressof(preconditioner)), + m_partition(std::move(partition)), + m_operator( + operation, + m_partition + ), + m_inverse( + preconditioner, + m_operator, + m_partition + ), + m_backend(prepareLinearBackend( + std::move(backend), + m_operator, + m_inverse, + communicator + )) { + if (!IsCurrent()) { + throw std::logic_error("Gravity subsystem preparation requires a current inverse."); + } + } + + PreparedGravitySubsystem(const PreparedGravitySubsystem &) = delete; + PreparedGravitySubsystem &operator=(const PreparedGravitySubsystem &) = delete; + PreparedGravitySubsystem(PreparedGravitySubsystem &&) = delete; + PreparedGravitySubsystem &operator=(PreparedGravitySubsystem &&) = delete; + + void Refresh() { + if (!m_fullPreconditioner->IsCurrent()) { + throw std::logic_error("Refresh the root inverse before refreshing its gravity subsystem."); + } + m_operator.Refresh(); + } + + [[nodiscard]] bool IsCurrent() const { + return m_operator.IsCurrent() && m_inverse.IsCurrent() && m_backend.IsReady(); + } + [[nodiscard]] const GravityBlockPartition &GetPartition() const noexcept { + return m_partition; + } + [[nodiscard]] const RestrictedOperator &GetOperator() const noexcept { + return m_operator; + } + [[nodiscard]] const Operation &GetFullOperator() const noexcept { + return *m_fullOperator; + } + [[nodiscard]] const Preconditioner &GetFullPreconditioner() const noexcept { + return *m_fullPreconditioner; + } + [[nodiscard]] MPI_Comm GetCommunicator() const noexcept { + return m_backend.GetCommunicator(); + } + [[nodiscard]] std::uint64_t GetSolveCount() const noexcept { + return m_solveCount; + } + [[nodiscard]] std::uint64_t GetIterationCount() const noexcept { + return m_iterationCount; + } + [[nodiscard]] std::uint64_t GetFailureCount() const noexcept { + return m_failureCount; + } + [[nodiscard]] double GetSolveSeconds() const noexcept { + return m_solveSeconds; + } + [[nodiscard]] const std::optional &GetLastSolveReport() const noexcept { + return m_lastSolve; + } + + [[nodiscard]] LinearSolveReport Solve( + const mfem::Vector &rightHandSide, + mfem::Vector &correction, + const LinearSolveControl &control + ) { + ++m_solveCount; + try { + const auto report = m_backend.Solve(rightHandSide, correction, control); + m_iterationCount += static_cast(report.iterations); + m_solveSeconds += report.solveSeconds; + if (!report.Converged()) { + ++m_failureCount; + } + m_lastSolve = report; + return report; + } catch (...) { + ++m_failureCount; + throw; + } + } + + private: + const Operation *m_fullOperator; + const Preconditioner *m_fullPreconditioner; + GravityBlockPartition m_partition; + RestrictedOperator m_operator; + RestrictedInverse m_inverse; + BackendType m_backend; + // Times and iteration counts cover completed calls; failed contracts + // are counted without attempting another collective during unwinding. + std::uint64_t m_solveCount{0}; + std::uint64_t m_iterationCount{0}; + std::uint64_t m_failureCount{0}; + double m_solveSeconds{0.0}; + std::optional m_lastSolve; + }; + + template + concept GravitySubsystemRuntime = requires( + Subsystem &subsystem, + const Subsystem &constant, + const mfem::Vector &rhs, + mfem::Vector &correction, + const LinearSolveControl &control + ) { + { constant.IsCurrent() } -> std::same_as; + { constant.GetPartition() } -> std::same_as; + { constant.GetCommunicator() } -> std::same_as; + { subsystem.Solve(rhs, correction, control) } -> std::same_as; + { constant.GetFullOperator().Mult(rhs, correction) } -> std::same_as; + { constant.GetFullOperator().GetProblem().GetPreparationGeneration() } -> std::convertible_to; + }; + + template class PreparedGravityCurvatureCorrection final { + public: + explicit PreparedGravityCurvatureCorrection(Subsystem &subsystem) + : m_subsystem(std::addressof(subsystem)), + m_baseResidual(subsystem.GetPartition().GravityResidualSize()), + m_baseAction(subsystem.GetPartition().GravityResidualSize()), + m_fullAction(subsystem.GetPartition().RootResidualSize()), + m_rightHandSide(subsystem.GetPartition().GravityResidualSize()) { + } + + PreparedGravityCurvatureCorrection(const PreparedGravityCurvatureCorrection &) = delete; + PreparedGravityCurvatureCorrection &operator=(const PreparedGravityCurvatureCorrection &) = delete; + PreparedGravityCurvatureCorrection(PreparedGravityCurvatureCorrection &&) = delete; + PreparedGravityCurvatureCorrection &operator=(PreparedGravityCurvatureCorrection &&) = delete; + + // Capture must precede trial preparation: using a candidate Jacobian + // here would introduce a first-order change into the Newton tangent. + void CaptureBase( + const mfem::Vector &baseResidual, + const mfem::Vector &baseDirection + ) { + m_hasBase = false; + const auto &partition = m_subsystem->GetPartition(); + RequireCollectively( + m_subsystem->IsCurrent() && baseResidual.Size() == partition.RootResidualSize() && + baseDirection.Size() == partition.RootValueSize(), + "Gravity curvature capture requires a current base and " + "complete vectors." + ); + partition.GatherGravityResiduals(baseResidual, m_baseResidual); + m_subsystem->GetFullOperator().Mult(baseDirection, m_fullAction); + partition.GatherGravityResiduals(m_fullAction, m_baseAction); + m_baseGeneration = m_subsystem->GetFullOperator().GetProblem().GetPreparationGeneration(); + m_hasBase = true; + } + + void Invalidate() noexcept { + m_hasBase = false; + } + [[nodiscard]] bool HasBase() const noexcept { + return m_hasBase; + } + [[nodiscard]] std::uint64_t GetBaseGeneration() const noexcept { + return m_baseGeneration; + } + + // 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. + [[nodiscard]] LinearSolveReport SolveCorrection( + const double alpha, + const mfem::Vector &trialResidual, + mfem::Vector &gravityCorrection, + const LinearSolveControl &control + ) { + const auto &partition = m_subsystem->GetPartition(); + RequireCollectively( + m_hasBase && m_subsystem->IsCurrent() && std::isfinite(alpha) && alpha >= 0.0 && + 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." + ); + } + 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); + } + gravityCorrection = 0.0; + return m_subsystem->Solve(m_rightHandSide, gravityCorrection, control); + } + + private: + void RequireCollectively( + const bool valid, + const char *message + ) const { + int local = valid ? 1 : 0; + int global = 0; + if (MPI_Allreduce(&local, &global, 1, MPI_INT, MPI_MIN, m_subsystem->GetCommunicator()) != MPI_SUCCESS) { + throw std::runtime_error( + "Gravity curvature correction could not verify " + "distributed validity." + ); + } + if (global == 0) { + throw std::logic_error(message); + } + } + + Subsystem *m_subsystem; + mfem::Vector m_baseResidual; + mfem::Vector m_baseAction; + mfem::Vector m_fullAction; + mfem::Vector m_rightHandSide; + std::uint64_t m_baseGeneration{0}; + bool m_hasBase{false}; + }; +} // namespace mean_field::solver diff --git a/libmeanfield/interface/solver/linear_backend.cppm b/libmeanfield/interface/solver/linear_backend.cppm index 88f2db2..4d1f198 100644 --- a/libmeanfield/interface/solver/linear_backend.cppm +++ b/libmeanfield/interface/solver/linear_backend.cppm @@ -664,12 +664,25 @@ export namespace mean_field::solver::linear { if (m_communicator == MPI_COMM_NULL) { throw std::invalid_argument("MFEM FGMRES requires a non-null MPI communicator."); } - if (operation.Height() <= 0 || operation.Width() <= 0 || operation.Height() != operation.Width() || - preconditioner.Height() != operation.Width() || preconditioner.Width() != operation.Height()) { + const bool localDimensionsAreValid = + operation.Height() >= 0 && operation.Width() >= 0 && operation.Height() == operation.Width() && + preconditioner.Height() == operation.Width() && preconditioner.Width() == operation.Height(); + if (!AllRanksAgree(localDimensionsAreValid, m_communicator)) { throw std::invalid_argument( - "MFEM FGMRES requires compatible square operator and preconditioner dimensions." + "MFEM FGMRES requires compatible square " + "operator and preconditioner dimensions." ); } + // A partition may own no local unknowns while still taking + // part in distributed operator applications and reductions. + const long long localSize = operation.Height(); + long long globalSize = 0; + if (MPI_Allreduce(&localSize, &globalSize, 1, MPI_LONG_LONG, MPI_SUM, m_communicator) != MPI_SUCCESS) { + throw std::runtime_error("MFEM FGMRES could not determine its global system size."); + } + if (globalSize == 0) { + throw std::invalid_argument("MFEM FGMRES requires a globally nonempty system."); + } m_rightHandSide = 0.0; m_operationAction = 0.0; @@ -699,6 +712,10 @@ export namespace mean_field::solver::linear { return m_communicator; } + [[nodiscard]] FGMRES GetConfiguration() const { + return FGMRES{m_options}; + } + [[nodiscard]] bool IsReady() const { return m_operation != nullptr && m_preconditioner != nullptr && m_communicator != MPI_COMM_NULL && m_operation->Height() == m_operation->Width() && diff --git a/tests/solver/gravity_subsystem.cpp b/tests/solver/gravity_subsystem.cpp new file mode 100644 index 0000000..2fc3fbe --- /dev/null +++ b/tests/solver/gravity_subsystem.cpp @@ -0,0 +1,411 @@ +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +import mean_field; + +namespace { + using namespace mean_field; + + class PreparedDenseOperator final : public mfem::Operator { + public: + explicit PreparedDenseOperator(mfem::DenseMatrix matrix) + : mfem::Operator( + matrix.Height(), + matrix.Width() + ), + m_matrix(std::move(matrix)) { + } + + void Mult( + const mfem::Vector &direction, + mfem::Vector &action + ) const override { + if (!m_prepared) { + throw std::logic_error("Test operator is unprepared."); + } + m_matrix.Mult(direction, action); + } + void Prepare(mfem::DenseMatrix matrix) { + m_matrix = std::move(matrix); + m_prepared = true; + ++m_generation; + } + void Invalidate() { + m_prepared = false; + } + [[nodiscard]] bool IsPrepared() const { + return m_prepared; + } + [[nodiscard]] const PreparedDenseOperator &GetProblem() const { + return *this; + } + [[nodiscard]] std::uint64_t GetPreparationGeneration() const { + return m_generation; + } + + private: + mfem::DenseMatrix m_matrix; + std::uint64_t m_generation{1}; + bool m_prepared{true}; + }; + + class PreparedPermutationInverse final : public mfem::Solver { + public: + static constexpr auto applicationContract = preconditioning::ApplicationContract::stationary_linear; + + explicit PreparedPermutationInverse( + const PreparedDenseOperator &operation, + std::vector residualForValue = {} + ) + : mfem::Solver(operation.Width()), + m_operation(operation), + m_residualForValue(std::move(residualForValue)) { + if (m_residualForValue.empty()) { + for (int index = 0; index < Width(); ++index) { + m_residualForValue.push_back(index); + } + } + Refresh(); + } + void SetOperator(const mfem::Operator &) override { + } + void Mult( + const mfem::Vector &rhs, + mfem::Vector &action + ) const override { + if (!IsCurrent()) { + throw std::logic_error("Test inverse is stale."); + } + for (int index = 0; index < Width(); ++index) { + action(index) = rhs(m_residualForValue[static_cast(index)]); + } + } + void Refresh() { + m_generation = m_operation.GetPreparationGeneration(); + } + [[nodiscard]] bool IsCurrent() const { + return m_operation.IsPrepared() && m_generation == m_operation.GetPreparationGeneration(); + } + + private: + const PreparedDenseOperator &m_operation; + std::vector m_residualForValue; + std::uint64_t m_generation{0}; + }; + + [[nodiscard]] mfem::DenseMatrix Diagonal(const std::initializer_list entries) { + mfem::DenseMatrix result(static_cast(entries.size())); + result = 0.0; + int index = 0; + for (const double entry : entries) { + result(index, index) = entry; + ++index; + } + return result; + } + + [[nodiscard]] mfem::DenseMatrix CurvedJacobian( + const double q, + const double z + ) { + mfem::DenseMatrix result(2); + result(0, 0) = 2.0 * q * z - 3.0 * q * q; + result(0, 1) = 1.0 + q * q; + result(1, 0) = 1.0; + result(1, 1) = 2.0; + return result; + } + + [[nodiscard]] mfem::Vector CurvedResidual( + const double q, + const double z + ) { + mfem::Vector result(2); + result(0) = (1.0 + q * q) * z - q * q * q; + result(1) = q + 2.0 * z; + return result; + } + + using PlainForm = utils::blocks::surface_deformed_stellar_equilibrium_form; + using CentralForm = utils::blocks::central_density_bordered_stellar_equilibrium_form; + static_assert(solver::GravityPartitionForm); + static_assert(solver::GravityPartitionForm); + static_assert(!solver::GravityPartitionForm); + static_assert(solver::GravitySubsystemOperation); + static_assert(solver::GravitySubsystemPreconditioner); + static_assert(!solver::GravitySubsystemOperation); + static_assert(!solver::GravitySubsystemPreconditioner); +} // namespace + +TEST_CASE( + "Gravity partition preserves noncontiguous coordinate and equation order", + "[solver][gravity_subsystem]" +) { + using Catch::Approx; + using namespace mean_field; + solver::GravityBlockPartition partition(5, 5, {3, 1}, {0, 4}); + mfem::Vector full(5), compact(2), retained(3), assembled(5); + for (int index = 0; index < 5; ++index) { + full(index) = 10.0 + index; + } + partition.GatherGravityValues(full, compact); + CHECK(compact(0) == 13.0); + CHECK(compact(1) == 11.0); + partition.InjectGravityResiduals(compact, assembled); + CHECK(assembled(0) == 13.0); + CHECK(assembled(4) == 11.0); + CHECK(assembled(1) == 0.0); + partition.GatherRetainedValues(full, retained); + CHECK(retained(0) == 10.0); + CHECK(retained(1) == 12.0); + CHECK(retained(2) == 14.0); + partition.InjectRetainedValues(retained, assembled); + partition.AddGravityValues(compact, assembled); + assembled -= full; + CHECK(assembled.Norml2() == 0.0); + CHECK_THROWS_AS(solver::GravityBlockPartition(3, 3, {0, 0}, {0, 1}), std::invalid_argument); + CHECK_THROWS_AS(solver::GravityBlockPartition(3, 3, {3}, {0}), std::invalid_argument); + CHECK_THROWS_AS(solver::GravityBlockPartition(3, 3, {0}, {0, 1}), std::invalid_argument); + mfem::Vector aliased(full.GetData() + 1, 2); + CHECK_THROWS_AS(partition.InjectGravityValues(aliased, full), std::invalid_argument); +} + +TEST_CASE( + "Gravity partition compiles from forms with and without a central border", + "[solver][gravity_subsystem]" +) { + using namespace mean_field; + const auto check = []() { + std::array valueSizes{}; + std::array residualSizes{}; + valueSizes.fill(1); + residualSizes.fill(1); + utils::blocks::form_layout layout(valueSizes, residualSizes); + const auto partition = solver::MakeGravityBlockPartition(layout); + CHECK(partition.GravityValueSize() == 2); + CHECK(partition.GravityResidualSize() == 2); + CHECK(partition.RetainedValueSize() == Form::value_block_count - 2); + CHECK(partition.RetainedResidualSize() == Form::residual_block_count - 2); + }; + check.template operator()(); + check.template operator()(); +} + +TEST_CASE( + "Gravity subsystem solves the actual restricted operator and rejects " + "stale preparation", + "[solver][gravity_subsystem]" +) { + using Catch::Approx; + using namespace mean_field; + PreparedDenseOperator operation(Diagonal({2.0, 3.0, 7.0})); + PreparedPermutationInverse inverse(operation); + solver::PreparedGravitySubsystem subsystem( + operation, inverse, solver::GravityBlockPartition(3, 3, {1, 0}, {1, 0}), MPI_COMM_WORLD + ); + static_assert(solver::GravitySubsystemRuntime); + mfem::Vector rhs(2), correction(2), action(2); + rhs(0) = 9.0; + rhs(1) = -8.0; + correction = 0.0; + const solver::LinearSolveControl control{.relativeTolerance = 1.0e-12, .maximumIterations = 8}; + const auto report = subsystem.Solve(rhs, correction, control); + REQUIRE(report.Converged()); + CHECK(correction(0) == Approx(3.0).margin(1.0e-12)); + CHECK(correction(1) == Approx(-4.0).margin(1.0e-12)); + subsystem.GetOperator().Mult(correction, action); + action -= rhs; + double localSquared = action * action; + double globalSquared = 0.0; + MPI_Allreduce(&localSquared, &globalSquared, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + CHECK(report.trueResidualNorm == Approx(std::sqrt(globalSquared)).margin(1.0e-14)); + + operation.Prepare(Diagonal({4.0, 6.0, 7.0})); + CHECK_FALSE(subsystem.IsCurrent()); + CHECK_THROWS_AS(subsystem.Solve(rhs, correction, control), std::logic_error); + CHECK_THROWS_AS(subsystem.Refresh(), std::logic_error); + inverse.Refresh(); + CHECK_FALSE(subsystem.IsCurrent()); + subsystem.Refresh(); + REQUIRE(subsystem.Solve(rhs, correction, control).Converged()); + CHECK(correction(0) == Approx(1.5).margin(1.0e-12)); + CHECK(correction(1) == Approx(-2.0).margin(1.0e-12)); + operation.Invalidate(); + CHECK_FALSE(subsystem.IsCurrent()); + CHECK_THROWS_AS(subsystem.Refresh(), std::logic_error); +} + +TEST_CASE( + "Gravity inner solve reports an unresolved true defect instead of " + "certifying a correction", + "[solver][gravity_subsystem]" +) { + using namespace mean_field; + PreparedDenseOperator operation(Diagonal({2.0, 3.0, 7.0})); + PreparedPermutationInverse inverse(operation); + solver::PreparedGravitySubsystem subsystem( + operation, inverse, solver::GravityBlockPartition(3, 3, {0, 1}, {0, 1}), MPI_COMM_WORLD + ); + mfem::Vector rhs(2), correction(2); + rhs = 1.0; + correction = 0.0; + const solver::LinearSolveControl control{.relativeTolerance = 1.0e-12, .maximumIterations = 1}; + const auto report = subsystem.Solve(rhs, correction, control); + CHECK_FALSE(report.Converged()); + CHECK(report.status == solver::LinearSolveStatus::maximum_iterations); + CHECK(report.trueResidualNorm > control.ConvergenceThreshold(report.rightHandSideNorm)); + CHECK(report.relativeTrueResidualNorm > 0.1); + CHECK(subsystem.GetSolveCount() == 1); + CHECK(subsystem.GetIterationCount() == static_cast(report.iterations)); + CHECK(subsystem.GetFailureCount() == 1); + CHECK(subsystem.GetSolveSeconds() == report.solveSeconds); + REQUIRE(subsystem.GetLastSolveReport().has_value()); + CHECK(subsystem.GetLastSolveReport()->trueResidualNorm == report.trueResidualNorm); + rhs(0) = std::numeric_limits::quiet_NaN(); + CHECK_THROWS_AS(subsystem.Solve(rhs, correction, control), std::invalid_argument); + CHECK(subsystem.GetSolveCount() == 2); + CHECK(subsystem.GetFailureCount() == 2); +} + +TEST_CASE( + "Distributed gravity solves admit empty local partitions but reject " + "a globally empty system", + "[solver][gravity_subsystem]" +) { + using Catch::Approx; + using namespace mean_field; + int rank = 0; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + const int localSize = rank == 0 ? 1 : 0; + mfem::DenseMatrix matrix(localSize); + if (localSize != 0) { + matrix(0, 0) = 4.0; + } + PreparedDenseOperator operation(matrix); + PreparedPermutationInverse inverse(operation); + std::vector indices = localSize != 0 ? std::vector{0} : std::vector{}; + solver::PreparedGravitySubsystem subsystem( + operation, inverse, solver::GravityBlockPartition(localSize, localSize, indices, indices), MPI_COMM_WORLD + ); + mfem::Vector rhs(localSize), correction(localSize); + rhs = 12.0; + correction = 0.0; + const auto report = subsystem.Solve(rhs, correction, {}); + REQUIRE(report.Converged()); + CHECK(report.rightHandSideNorm == 12.0); + if (localSize != 0) { + CHECK(correction(0) == Approx(3.0)); + } + + auto backend = solver::linear::prepareLinearBackend( + solver::linear::FGMRES{{.restartLength = 13, .printLevel = -1}}, operation, inverse, MPI_COMM_WORLD + ); + CHECK(backend.GetConfiguration().GetOptions().restartLength == 13); + CHECK(backend.GetConfiguration().GetOptions().printLevel == -1); + + PreparedDenseOperator emptyOperation(mfem::DenseMatrix(0)); + PreparedPermutationInverse emptyInverse(emptyOperation); + CHECK_THROWS_AS( + (solver::PreparedGravitySubsystem( + emptyOperation, emptyInverse, solver::GravityBlockPartition(0, 0, {}, {}), MPI_COMM_WORLD + )), + std::invalid_argument + ); +} + +TEST_CASE( + "Affine gravity receives no curvature correction even away from " + "gravity equilibrium", + "[solver][gravity_subsystem]" +) { + using namespace mean_field; + mfem::DenseMatrix jacobian(2); + jacobian(0, 0) = 2.0; + jacobian(0, 1) = 3.0; + jacobian(1, 0) = 1.0; + jacobian(1, 1) = 1.0; + PreparedDenseOperator operation(jacobian); + PreparedPermutationInverse inverse(operation, {1, 0}); + solver::PreparedGravitySubsystem subsystem( + 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); + baseResidual(0) = 0.75; + baseResidual(1) = 0.75; + direction(0) = 0.125; + direction(1) = 0.5; + correction.CaptureBase(baseResidual, direction); + const auto baseGeneration = correction.GetBaseGeneration(); + for (const double alpha : {1.0, 0.5, 0.25, 0.125, 0.0}) { + operation.Prepare(jacobian); + inverse.Refresh(); + subsystem.Refresh(); + trialResidual(0) = 0.75 + alpha * 1.75; + trialResidual(1) = 0.75 + alpha * 0.625; + const auto report = correction.SolveCorrection(alpha, trialResidual, gravityStep, {}); + REQUIRE(report.Converged()); + CHECK(report.rightHandSideNorm == 0.0); + CHECK(gravityStep.Norml2() == 0.0); + CHECK(correction.GetBaseGeneration() == baseGeneration); + } + correction.Invalidate(); + CHECK_FALSE(correction.HasBase()); + CHECK_THROWS_AS(correction.SolveCorrection(1.0, trialResidual, gravityStep, {}), std::logic_error); +} + +TEST_CASE( + "Gravity curvature correction preserves the frozen tangent and is " + "second order", + "[solver][gravity_subsystem]" +) { + using Catch::Approx; + using namespace mean_field; + constexpr double baseQ = 0.4, baseZ = 1.2, deltaQ = 0.3, deltaZ = -0.2; + PreparedDenseOperator operation(CurvedJacobian(baseQ, baseZ)); + PreparedPermutationInverse inverse(operation, {1, 0}); + solver::PreparedGravitySubsystem subsystem( + operation, inverse, solver::GravityBlockPartition(2, 2, {1}, {0}), MPI_COMM_WORLD + ); + solver::PreparedGravityCurvatureCorrection correction(subsystem); + const auto baseResidual = CurvedResidual(baseQ, baseZ); + mfem::Vector direction(2), action(2), gravityStep(1), candidate(2); + direction(0) = deltaQ; + direction(1) = deltaZ; + operation.Mult(direction, action); + correction.CaptureBase(baseResidual, direction); + double previousStep = 0.0; + for (const double alpha : {0.02, 0.01, 0.005, 0.0025}) { + candidate(0) = baseQ + alpha * deltaQ; + candidate(1) = baseZ + alpha * deltaZ; + operation.Prepare(CurvedJacobian(candidate(0), candidate(1))); + inverse.Refresh(); + subsystem.Refresh(); + const auto trialResidual = CurvedResidual(candidate(0), candidate(1)); + const auto report = correction.SolveCorrection( + alpha, trialResidual, gravityStep, {.relativeTolerance = 1.0e-12, .maximumIterations = 4} + ); + REQUIRE(report.Converged()); + subsystem.GetPartition().AddGravityValues(gravityStep, candidate); + const auto correctedResidual = CurvedResidual(candidate(0), candidate(1)); + CHECK(correctedResidual(0) == Approx(baseResidual(0) + alpha * action(0)).margin(2.0e-15)); + CHECK(candidate(0) == baseQ + alpha * deltaQ); + if (previousStep != 0.0) { + CHECK(std::abs(gravityStep(0) / previousStep - 0.25) < 0.004); + } + previousStep = gravityStep(0); + } + CHECK(baseResidual(0) == CurvedResidual(baseQ, baseZ)(0)); + CHECK_THROWS_AS(correction.SolveCorrection(-1.0, baseResidual, gravityStep, {}), std::logic_error); +}