From 995d49b81347a1c68609e248a4d21c485bf2d9d5 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Thu, 17 Sep 2026 16:35:10 -0400 Subject: [PATCH] Prototype deterministic coupled shape coarse inverse for frozen benchmarks --- experiments/plateau_coupled_coarse.hpp | 697 ++++++++++++++++++++++++ tests/solver/plateau_coupled_coarse.cpp | 392 +++++++++++++ 2 files changed, 1089 insertions(+) create mode 100644 experiments/plateau_coupled_coarse.hpp create mode 100644 tests/solver/plateau_coupled_coarse.cpp diff --git a/experiments/plateau_coupled_coarse.hpp b/experiments/plateau_coupled_coarse.hpp new file mode 100644 index 0000000..e854c3d --- /dev/null +++ b/experiments/plateau_coupled_coarse.hpp @@ -0,0 +1,697 @@ +#pragma once + +// Include after importing mean_field. This is an experiment-local inverse; +// it does not replace or prepare the authoritative normalized Jacobian. +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +namespace experiment::plateau_diagnostics { + struct CoupledShapeCoarseOptions final { + double relativeRankTolerance{1.e-10}; + double minimumSubspaceOverlap{1.e-8}; + }; + + enum class CoupledShapeCoarseFallback { none, no_independent_modes, insufficient_subspace_overlap }; + enum class CoupledShapeCoarseModeStatus { retained, dependent_source, dependent_filtered_image }; + + struct CoupledShapeCoarseModeReport final { + // QR admission is retained even if the complete space later falls back. + CoupledShapeCoarseModeStatus status{CoupledShapeCoarseModeStatus::dependent_source}; + double normalizedInputNorm{0.0}; + double sourceImageNorm{0.0}; + double sourceIndependentFraction{0.0}; + double filteredImageNorm{0.0}; + double filteredIndependentFraction{0.0}; + }; + + struct CoupledShapeCoarseStatistics final { + int inputModes{0}; + int independentModes{0}; + int activeModes{0}; + int sourceDependentModes{0}; + int imageDependentModes{0}; + CoupledShapeCoarseFallback fallback{CoupledShapeCoarseFallback::no_independent_modes}; + double reciprocalOverlapCondition{std::numeric_limits::quiet_NaN()}; + double minimumOverlapLowerBound{std::numeric_limits::quiet_NaN()}; + std::uint64_t setups{0}; + std::uint64_t setupOperatorApplications{0}; + std::uint64_t setupInverseApplications{0}; + double lastSetupSeconds{0.0}; + double lastSetupOperatorSeconds{0.0}; + double lastSetupInverseSeconds{0.0}; + double totalSetupSeconds{0.0}; + // Application clocks are rank-local; aggregate them outside the timed + // solve. The only application reductions validate and project vectors. + std::uint64_t applications{0}; + std::uint64_t inverseApplications{0}; + std::uint64_t failedApplications{0}; + double applicationSecondsLocal{0.0}; + double inverseSecondsLocal{0.0}; + std::size_t persistentVectorBytes{0}; + }; + + template + concept CoupledCoarseNormalizedOperator = + std::derived_from && + requires(const Candidate &operation, const mfem::Vector &physical, mfem::Vector &normalized) { + { operation.IsPrepared() } -> std::same_as; + { operation.GetProblem().GetPreparationGeneration() } -> std::convertible_to; + { operation.NormalizeState(physical, normalized) } -> std::same_as; + }; + + template + concept CoupledCoarseStationaryInverse = + std::derived_from && + mean_field::solver::LinearPreconditionerApplicationContractAvailable && + (mean_field::solver::linearPreconditionerApplicationContract == + mean_field::preconditioning::ApplicationContract::stationary_linear) && + requires(const Candidate &inverse) { + { inverse.IsCurrent() } -> std::same_as; + }; + + namespace coupled_shape_detail { + template + using RootBlock = std::remove_cvref_t().GetBlock())>; + template + using MaterialBlock = typename RootBlock::StructureComponent::MaterialSurfaceComponent; + + template + concept SupportsShapeModes = requires(const Problem &problem, const Inverse &inverse) { + typename MaterialBlock::Factorization; + requires std::same_as< + typename MaterialBlock::Factorization, + mean_field::preconditioning::SurfaceThenMaterialTriangular>; + { inverse.GetProblem() } -> std::same_as; + { inverse.IsCurrent() } -> std::same_as; + inverse.GetCoordinateMap().GetCorrectionRanges(); + inverse.GetGroupedPreconditioner().GetStructurePreconditioner().GetCrossCouplings(); + inverse.GetGroupedPreconditioner() + .GetStructurePreconditioner() + .GetMaterialSurfacePreconditioner() + .GetCoupledOperator(); + inverse.GetGroupedPreconditioner().GetStructurePreconditioner().GetGravityPreconditioner(); + { problem.GetPreparationGeneration() } -> std::convertible_to; + { + problem.GetPhysicalOperator().GetDomainDeformation().surfaceDeformationPrescription().radialDirection( + 0, 0 + ) + } -> std::same_as; + { + problem.GetPhysicalOperator().GetDomainDeformation().surfaceDeformationPrescription().referenceRadius(0) + } -> std::same_as; + }; + } // namespace coupled_shape_detail + + // These are geometry modes, not Newton-RHS samples. Complete their material + // and gravity response with existing fixed inverses; neither completion is + // an exact elimination. Scalar unknowns stay zero, while the later full-J + // images retain every scalar equation, including an optional central border. + // Construction performs nine explicit surface-to-material actions, material + // inverse applications, material-to-gravity actions, and gravity applications. + // The material inverse also performs its usual internal coupling actions. + template < + class Problem, + class PhysicalInverse> + requires coupled_shape_detail::SupportsShapeModes< + Problem, + PhysicalInverse> + [[nodiscard]] std::vector MakeCoupledShapeModes( + const Problem &problem, + const PhysicalInverse &inverse + ) { + using namespace mean_field; + using namespace utils::blocks; + using Material = coupled_shape_detail::MaterialBlock; + using Root = coupled_shape_detail::RootBlock; + using Corrections = typename Material::CorrectionBlocks; + using Residuals = typename Material::ResidualBlocks; + using Gravity = typename Root::StructureComponent::GravityComponent; + // The borrowed coupling operator exposes this exact three-field packing. + // Constraining the descriptor prevents a superficially compatible custom + // material recipe from silently assigning a physical mode to wrong rows. + static_assert( + std::same_as< + Corrections, + type_list> + ); + static_assert(std::same_as< + Residuals, type_list< + density::mass::residual, surface_deformation::shape_equilibrium::residual, + enthalpy::specific::residual>>); + static_assert( + std::same_as< + typename Gravity::CorrectionBlocks, type_list> + ); + static_assert( + std::same_as< + typename Gravity::ResidualBlocks, type_list> + ); + constexpr auto densityRow = type_index_v; + constexpr auto enthalpyRow = type_index_v; + constexpr auto shapeColumn = type_index_v; + constexpr auto materialFirst = type_index_v; + constexpr auto gravityFirst = type_index_v; + const auto communicator = problem.GetCommunicator(); + if (communicator == MPI_COMM_NULL) + throw std::invalid_argument("Physical shape modes require a valid communicator."); + const auto all = [communicator](bool condition) { + int local = condition, global = 0; + if (MPI_Allreduce(&local, &global, 1, MPI_INT, MPI_MIN, communicator) != MPI_SUCCESS) + throw std::runtime_error("Could not validate collective physical shape modes."); + return global != 0; + }; + if (!all( + &inverse.GetProblem() == &problem && inverse.IsCurrent() && + problem.GetPhysicalOperator().IsPrepared() && !inverse.iterative_mode + )) + throw std::logic_error("Physical shape modes require the current inverse of this problem."); + const auto generation = problem.GetPreparationGeneration(); + const auto &structure = inverse.GetGroupedPreconditioner().GetStructurePreconditioner(); + const auto &material = structure.GetMaterialSurfacePreconditioner(); + const auto &gravityInverse = structure.GetGravityPreconditioner(); + const auto &couplings = material.GetCoupledOperator(); + const auto &cross = structure.GetCrossCouplings(); + const auto &offsets = couplings.GetOffsets(); + const auto &coordinates = inverse.GetCoordinateMap(); + const auto &ranges = coordinates.GetCorrectionRanges(); + const auto &radial = problem.GetPhysicalOperator().GetDomainDeformation().surfaceDeformationPrescription(); + const auto materialOffset = ranges[materialFirst].preconditionerOffset; + const auto gravityOffset = ranges[gravityFirst].preconditionerOffset; + const auto shapeSize = offsets[shapeColumn + 1] - offsets[shapeColumn]; + if (!all( + radial.spatialDimension() == 3 && radial.parameterCount() == shapeSize && + material.Height() == material.Width() && gravityInverse.Height() == gravityInverse.Width() && + material.Height() == cross.MaterialSize() && gravityInverse.Height() == cross.GravitySize() && + materialOffset == 0 && gravityOffset == material.Height() && + gravityOffset + gravityInverse.Height() <= coordinates.PreconditionerCorrectionSize() && + !material.iterative_mode && !gravityInverse.iterative_mode + )) + throw std::invalid_argument("Physical shape modes require compatible radial three-dimensional blocks."); + mfem::Vector materialRhs(material.Width()), materialAction(material.Height()); + mfem::Vector gravityRhs(gravityInverse.Width()), gravityAction(gravityInverse.Height()); + mfem::Vector grouped(coordinates.PreconditionerCorrectionSize()); + mfem::Vector shape(shapeSize); + mfem::Vector densityRhs(materialRhs, offsets[densityRow], offsets[densityRow + 1] - offsets[densityRow]); + mfem::Vector enthalpyRhs(materialRhs, offsets[enthalpyRow], offsets[enthalpyRow + 1] - offsets[enthalpyRow]); + mfem::Vector shapeAction(materialAction, offsets[shapeColumn], shapeSize); + mfem::Vector groupedMaterial(grouped, materialOffset, materialAction.Size()); + mfem::Vector groupedGravity(grouped, gravityOffset, gravityAction.Size()); + std::vector result; + result.reserve(9); + for (int mode = 0; mode < 9; ++mode) { + bool finite = true; + for (int i = 0; i < shapeSize; ++i) { + const double x = radial.radialDirection(i, 0), y = radial.radialDirection(i, 1), + z = radial.radialDirection(i, 2), radius = radial.referenceRadius(i); + // All real harmonics through l=2, up to irrelevant constants. + // Their span is rotation invariant; individual columns are not + // asserted to be exact translation or homology null directions. + const std::array angular{1.0, x, y, z, x * y, x * z, y * z, x * x - y * y, 3 * z * z - 1}; + shape(i) = radius * angular[mode]; + finite = finite && std::isfinite(shape(i)) && radius > 0.0; + } + if (!all(finite)) + throw std::runtime_error("A physical radial shape mode is nonfinite or has invalid radius."); + materialRhs = 0.0; + couplings.ApplySurfaceToMaterial(shape, densityRhs, enthalpyRhs); + densityRhs *= -1.0; + enthalpyRhs *= -1.0; + densityRhs.SyncAliasMemory(materialRhs); + enthalpyRhs.SyncAliasMemory(materialRhs); + // Zero shape RHS in this policy holds its returned shape at zero. + // We then restore the prescribed shape rather than solving A_qq. + material.Mult(materialRhs, materialAction); + shapeAction = shape; + shapeAction.SyncAliasMemory(materialAction); + cross.ApplyMaterialToGravity(materialAction, gravityRhs); + gravityRhs *= -1.0; + gravityInverse.Mult(gravityRhs, gravityAction); + grouped = 0.0; + groupedMaterial = materialAction; + groupedGravity = gravityAction; + groupedMaterial.SyncAliasMemory(grouped); + groupedGravity.SyncAliasMemory(grouped); + result.emplace_back(problem.StateSize()); + coordinates.UnpackCorrection(grouped, result.back()); + finite = true; + for (int i = 0; i < result.back().Size(); ++i) + finite = finite && std::isfinite(result.back()(i)); + if (!all(finite && inverse.IsCurrent() && generation == problem.GetPreparationGeneration())) + throw std::runtime_error("Physical shape completion failed or its preparation changed."); + } + return result; + } + + /* + * The caller supplies deterministic, physically completed shape directions + * before exposing a Newton RHS. With A the full normalized Jacobian, form + * V = orth(A T), U = P0 V, and A U = Q R. Then apply + * + * P r = P0 [r + (V R^-1 - Q) Q^T r]. + * + * This equals U R^-1 Q^T r + P0 (I-Q Q^T)r and is exact on range(A U) + * in exact arithmetic. The P0 A filtering can change the supplied physical + * correction space; this is not an exact completed-shape Schur inverse. + * + * For stationary P0, the extra right factor is nonsingular iff Q^T V is, + * provided R is nonsingular. Orthonormal V and Q let the small overlap + * matrix expose near-orthogonal subspaces. A failed overlap check selects + * P0 explicitly; it never regularizes A or adds a physical constraint. + * Invertibility/quality of P0 itself and full-system true-defect acceptance + * remain the caller's responsibility. No solver improvement is assumed. + */ + template + class PreparedCoupledShapeCoarse final : public mfem::Solver { + public: + static constexpr auto applicationContract = mean_field::preconditioning::ApplicationContract::stationary_linear; + static constexpr int maximumModes = 16; + + PreparedCoupledShapeCoarse( + const Operation &operation, + const Inverse &inverse, + std::span physicalModes, + MPI_Comm communicator, + CoupledShapeCoarseOptions options = {} + ) + : mfem::Solver( + operation.Width(), + operation.Height() + ), + m_operation(&operation), + m_inverse(&inverse), + m_communicator(communicator), + m_options(options), + m_residualWork(operation.Height()), + m_actionWork(operation.Width()) { + Refresh(physicalModes); + } + + PreparedCoupledShapeCoarse(const PreparedCoupledShapeCoarse &) = delete; + PreparedCoupledShapeCoarse &operator=(const PreparedCoupledShapeCoarse &) = delete; + PreparedCoupledShapeCoarse(PreparedCoupledShapeCoarse &&) = delete; + PreparedCoupledShapeCoarse &operator=(PreparedCoupledShapeCoarse &&) = delete; + + void SetOperator(const mfem::Operator &operation) override { + if (&operation != m_operation) + throw std::invalid_argument("A coupled coarse inverse cannot bind a different Jacobian."); + RequireCurrent(); + } + + [[nodiscard]] bool IsCurrent() const { + return m_ready && m_operation->IsPrepared() && m_inverse->IsCurrent() && !m_inverse->iterative_mode && + m_generation == m_operation->GetProblem().GetPreparationGeneration(); + } + + [[nodiscard]] bool UsesCoarseCorrection() const noexcept { + return !m_sources.empty(); + } + [[nodiscard]] MPI_Comm GetCommunicator() const noexcept { + return m_communicator; + } + [[nodiscard]] const CoupledShapeCoarseStatistics &GetStatistics() const noexcept { + return m_statistics; + } + [[nodiscard]] std::span GetModeReports() const noexcept { + return m_reports; + } + + // These normalized residual-space columns permit independent checks: + // for each v, compare P(A(P0 v)) with P0 v using the authoritative A. + [[nodiscard]] std::span GetSourceBasis() const noexcept { + return m_sources; + } + + void Refresh(std::span physicalModes) { + m_ready = false; + const double start = MPI_Wtime(); + ValidatePreparation(physicalModes); + m_generation = m_operation->GetProblem().GetPreparationGeneration(); + m_sources.clear(); + m_images.clear(); + m_reports.clear(); + m_sources.reserve(physicalModes.size()); + m_images.reserve(physicalModes.size()); + m_reports.reserve(physicalModes.size()); + m_upper.fill(0.0); + m_statistics.inputModes = static_cast(physicalModes.size()); + m_statistics.independentModes = m_statistics.activeModes = 0; + m_statistics.sourceDependentModes = m_statistics.imageDependentModes = 0; + m_statistics.fallback = CoupledShapeCoarseFallback::none; + m_statistics.reciprocalOverlapCondition = std::numeric_limits::quiet_NaN(); + m_statistics.minimumOverlapLowerBound = std::numeric_limits::quiet_NaN(); + double operatorSeconds = 0.0, inverseSeconds = 0.0; + mfem::Vector value(Height()), source(Width()), image(Width()); + + for (const auto &physical : physicalModes) { + CoupledShapeCoarseModeReport report; + m_operation->NormalizeState(physical, value); + RequireFinite(value); + report.normalizedInputNorm = Norm(value); + if (report.normalizedInputNorm == 0.0) { + ++m_statistics.sourceDependentModes; + m_reports.push_back(report); + continue; + } + value /= report.normalizedInputNorm; + double actionStart = MPI_Wtime(); + m_operation->Mult(value, source); + operatorSeconds += MPI_Wtime() - actionStart; + ++m_statistics.setupOperatorApplications; + RequireFinite(source); + report.sourceImageNorm = Norm(source); + std::array ignored{}; + Reorthogonalize(source, m_sources, ignored); + const double sourceNorm = Norm(source); + report.sourceIndependentFraction = Fraction(sourceNorm, report.sourceImageNorm); + if (sourceNorm == 0.0 || report.sourceIndependentFraction <= m_options.relativeRankTolerance) { + ++m_statistics.sourceDependentModes; + m_reports.push_back(report); + continue; + } + source /= sourceNorm; + actionStart = MPI_Wtime(); + m_inverse->Mult(source, value); + inverseSeconds += MPI_Wtime() - actionStart; + ++m_statistics.setupInverseApplications; + RequireFinite(value); + actionStart = MPI_Wtime(); + m_operation->Mult(value, image); + operatorSeconds += MPI_Wtime() - actionStart; + ++m_statistics.setupOperatorApplications; + RequireFinite(image); + report.filteredImageNorm = Norm(image); + std::array upperColumn{}; + Reorthogonalize(image, m_images, upperColumn); + const double imageNorm = Norm(image); + report.filteredIndependentFraction = Fraction(imageNorm, report.filteredImageNorm); + if (imageNorm == 0.0 || report.filteredIndependentFraction <= m_options.relativeRankTolerance) { + report.status = CoupledShapeCoarseModeStatus::dependent_filtered_image; + ++m_statistics.imageDependentModes; + m_reports.push_back(report); + continue; + } + const int column = static_cast(m_sources.size()); + for (int row = 0; row < column; ++row) + m_upper[row * maximumModes + column] = upperColumn[row]; + m_upper[column * maximumModes + column] = imageNorm; + image /= imageNorm; + m_sources.emplace_back(source); + m_images.emplace_back(image); + report.status = CoupledShapeCoarseModeStatus::retained; + m_reports.push_back(report); + } + m_statistics.independentModes = static_cast(m_sources.size()); + if (m_sources.empty()) { + m_statistics.fallback = CoupledShapeCoarseFallback::no_independent_modes; + } else { + const auto [condition, overlap] = CheckOverlap(); + m_statistics.reciprocalOverlapCondition = condition; + m_statistics.minimumOverlapLowerBound = overlap; + if (!std::isfinite(overlap) || overlap < m_options.minimumSubspaceOverlap) { + m_statistics.fallback = CoupledShapeCoarseFallback::insufficient_subspace_overlap; + m_sources.clear(); + m_images.clear(); + } + } + m_statistics.activeModes = static_cast(m_sources.size()); + // Payload bytes exclude allocator metadata and borrowed objects. + m_statistics.persistentVectorBytes = + sizeof(mfem::real_t) * (static_cast(m_residualWork.Size()) + m_actionWork.Size() + + 2 * m_sources.size() * static_cast(Width())); + if (!All( + m_operation->IsPrepared() && m_inverse->IsCurrent() && + m_generation == m_operation->GetProblem().GetPreparationGeneration() + )) + throw std::logic_error("The coupled coarse preparation changed during setup."); + m_statistics.lastSetupOperatorSeconds = Maximum(operatorSeconds); + m_statistics.lastSetupInverseSeconds = Maximum(inverseSeconds); + m_statistics.lastSetupSeconds = Maximum(MPI_Wtime() - start); + m_statistics.totalSetupSeconds += m_statistics.lastSetupSeconds; + ++m_statistics.setups; + m_ready = true; + } + + void Mult( + const mfem::Vector &rightHandSide, + mfem::Vector &action + ) const override { + const double start = MPI_Wtime(); + ++m_statistics.applications; + try { + RequireCurrent(); + if (!All(rightHandSide.Size() == Width() && action.Size() == Height())) + throw std::invalid_argument("A coupled coarse application requires preallocated root vectors."); + RequireFinite(rightHandSide); + m_residualWork = rightHandSide; + const int count = static_cast(m_sources.size()); + if (count != 0) { + std::array projection{}, coefficients{}; + Project(rightHandSide, m_images, projection); + coefficients = projection; + for (int row = count - 1; row >= 0; --row) { + for (int column = row + 1; column < count; ++column) + coefficients[row] -= m_upper[row * maximumModes + column] * coefficients[column]; + coefficients[row] /= m_upper[row * maximumModes + row]; + } + for (int column = 0; column < count; ++column) { + m_residualWork.Add(coefficients[column], m_sources[column]); + m_residualWork.Add(-projection[column], m_images[column]); + } + RequireFinite(m_residualWork); + } + const double inverseStart = MPI_Wtime(); + ++m_statistics.inverseApplications; + m_inverse->Mult(m_residualWork, m_actionWork); + m_statistics.inverseSecondsLocal += MPI_Wtime() - inverseStart; + RequireFinite(m_actionWork); + RequireCurrent(); + action = m_actionWork; + m_statistics.applicationSecondsLocal += MPI_Wtime() - start; + } catch (...) { + ++m_statistics.failedApplications; + m_statistics.applicationSecondsLocal += MPI_Wtime() - start; + throw; + } + } + + private: + [[nodiscard]] static double Fraction( + double independent, + double original + ) { + return original == 0.0 ? 0.0 : independent / original; + } + + [[nodiscard]] bool All(bool value) const { + int local = value, global = 0; + if (MPI_Allreduce(&local, &global, 1, MPI_INT, MPI_MIN, m_communicator) != MPI_SUCCESS) + throw std::runtime_error("Could not validate a collective coupled coarse operation."); + return global != 0; + } + + [[nodiscard]] double Maximum(double value) const { + double result = 0.0; + if (MPI_Allreduce(&value, &result, 1, MPI_DOUBLE, MPI_MAX, m_communicator) != MPI_SUCCESS) + throw std::runtime_error("Could not measure coupled coarse setup time."); + return result; + } + + void RequireCurrent() const { + if (!All(IsCurrent())) + throw std::logic_error("The coupled coarse inverse is stale; refresh its dependencies and basis."); + } + + void RequireFinite(const mfem::Vector &vector) const { + bool finite = true; + for (int i = 0; i < vector.Size(); ++i) + finite = finite && std::isfinite(vector(i)); + if (!All(finite)) + throw std::runtime_error("A coupled coarse vector is nonfinite."); + } + + [[nodiscard]] double Norm(const mfem::Vector &vector) const { + const double local = vector * vector; + double global = 0.0; + if (MPI_Allreduce(&local, &global, 1, MPI_DOUBLE, MPI_SUM, m_communicator) != MPI_SUCCESS) + throw std::runtime_error("Could not compute a coupled coarse norm."); + if (!std::isfinite(global) || global < 0.0) + throw std::runtime_error("A coupled coarse norm is nonfinite."); + return std::sqrt(global); + } + + void Project( + const mfem::Vector &vector, + const std::vector &basis, + std::array< + double, + maximumModes> &coefficients + ) const { + std::array local{}; + for (std::size_t i = 0; i < basis.size(); ++i) + local[i] = basis[i] * vector; + if (MPI_Allreduce( + local.data(), coefficients.data(), static_cast(basis.size()), MPI_DOUBLE, MPI_SUM, + m_communicator + ) != MPI_SUCCESS) + throw std::runtime_error("Could not project a coupled coarse vector."); + for (std::size_t i = 0; i < basis.size(); ++i) + if (!std::isfinite(coefficients[i])) + throw std::runtime_error("A coupled coarse projection is nonfinite."); + } + + void Reorthogonalize( + mfem::Vector &vector, + const std::vector &basis, + std::array< + double, + maximumModes> &coefficients + ) const { + for (int pass = 0; pass < 2; ++pass) { + std::array projection{}; + Project(vector, basis, projection); + for (std::size_t i = 0; i < basis.size(); ++i) { + coefficients[i] += projection[i]; + vector.Add(-projection[i], basis[i]); + } + } + } + + void ValidatePreparation(std::span modes) const { + if (m_communicator == MPI_COMM_NULL) + throw std::invalid_argument("A coupled coarse inverse requires a valid communicator."); + const int count = modes.size() <= maximumModes ? static_cast(modes.size()) : maximumModes + 1; + int minimumCount = 0, maximumCount = 0; + if (MPI_Allreduce(&count, &minimumCount, 1, MPI_INT, MPI_MIN, m_communicator) != MPI_SUCCESS || + MPI_Allreduce(&count, &maximumCount, 1, MPI_INT, MPI_MAX, m_communicator) != MPI_SUCCESS) + throw std::runtime_error("Could not compare coupled coarse mode counts."); + const std::array options{m_options.relativeRankTolerance, m_options.minimumSubspaceOverlap}; + std::array minimumOptions{}, maximumOptions{}; + if (MPI_Allreduce(options.data(), minimumOptions.data(), 2, MPI_DOUBLE, MPI_MIN, m_communicator) != + MPI_SUCCESS || + MPI_Allreduce(options.data(), maximumOptions.data(), 2, MPI_DOUBLE, MPI_MAX, m_communicator) != + MPI_SUCCESS) + throw std::runtime_error("Could not compare coupled coarse preparation options."); + bool valid = count > 0 && count <= maximumModes && minimumCount == maximumCount && + minimumOptions == maximumOptions && std::isfinite(options[0]) && options[0] > 0.0 && + options[0] < 1.0 && std::isfinite(options[1]) && options[1] > 0.0 && options[1] < 1.0 && + m_operation->IsPrepared() && m_inverse->IsCurrent() && !m_inverse->iterative_mode && + m_operation->Height() == m_operation->Width() && m_inverse->Height() == Height() && + m_inverse->Width() == Width(); + if constexpr (requires { m_inverse->GetNormalizedJacobian(); }) + valid = valid && &m_inverse->GetNormalizedJacobian() == m_operation; + if (!All(valid)) + throw std::invalid_argument( + "Coupled coarse setup needs current compatible objects and common valid options." + ); + const long long localSize = Width(); + long long globalSize = 0; + if (MPI_Allreduce(&localSize, &globalSize, 1, MPI_LONG_LONG, MPI_SUM, m_communicator) != MPI_SUCCESS || + globalSize <= 0) + throw std::invalid_argument("A coupled coarse system must be globally nonempty."); + for (const auto &mode : modes) { + if (!All(mode.Size() == Height())) + throw std::invalid_argument("Coupled coarse modes must use complete physical root coordinates."); + RequireFinite(mode); + } + } + + [[nodiscard]] std::pair< + double, + double> + CheckOverlap() const { + const int count = static_cast(m_sources.size()); + std::array local{}, matrix{}, inverse{}; + for (int row = 0; row < count; ++row) + for (int column = 0; column < count; ++column) + local[row * count + column] = m_images[row] * m_sources[column]; + if (MPI_Allreduce(local.data(), matrix.data(), count * count, MPI_DOUBLE, MPI_SUM, m_communicator) != + MPI_SUCCESS) + throw std::runtime_error("Could not assemble the coupled coarse overlap matrix."); + for (int index = 0; index < count * count; ++index) + if (!std::isfinite(matrix[index])) + return {0.0, 0.0}; + double norm = 0.0; + for (int row = 0; row < count; ++row) { + inverse[row * count + row] = 1.0; + double sum = 0.0; + for (int column = 0; column < count; ++column) + sum += std::abs(matrix[row * count + column]); + norm = std::max(norm, sum); + } + if (!std::isfinite(norm) || norm == 0.0) + return {0.0, 0.0}; + for (int column = 0; column < count; ++column) { + int pivot = column; + for (int row = column + 1; row < count; ++row) + if (std::abs(matrix[row * count + column]) > std::abs(matrix[pivot * count + column])) + pivot = row; + if (!std::isfinite(matrix[pivot * count + column]) || + std::abs(matrix[pivot * count + column]) <= std::numeric_limits::epsilon() * count * norm) + return {0.0, 0.0}; + for (int j = 0; j < count; ++j) { + std::swap(matrix[column * count + j], matrix[pivot * count + j]); + std::swap(inverse[column * count + j], inverse[pivot * count + j]); + } + const double diagonal = matrix[column * count + column]; + for (int j = 0; j < count; ++j) { + matrix[column * count + j] /= diagonal; + inverse[column * count + j] /= diagonal; + } + for (int row = 0; row < count; ++row) + if (row != column) { + const double factor = matrix[row * count + column]; + for (int j = 0; j < count; ++j) { + matrix[row * count + j] -= factor * matrix[column * count + j]; + inverse[row * count + j] -= factor * inverse[column * count + j]; + } + } + } + double inverseOne = 0.0, inverseInfinity = 0.0; + for (int index = 0; index < count * count; ++index) + if (!std::isfinite(inverse[index])) + return {0.0, 0.0}; + for (int i = 0; i < count; ++i) { + double rowSum = 0.0, columnSum = 0.0; + for (int j = 0; j < count; ++j) { + rowSum += std::abs(inverse[i * count + j]); + columnSum += std::abs(inverse[j * count + i]); + } + inverseInfinity = std::max(inverseInfinity, rowSum); + inverseOne = std::max(inverseOne, columnSum); + } + if (!std::isfinite(inverseOne) || !std::isfinite(inverseInfinity)) + return {0.0, 0.0}; + // Condition number alone misses epsilon*I. Since ||S^-1||_2 is + // bounded by sqrt(||S^-1||_1 ||S^-1||_inf), this also checks the + // absolute overlap of the two orthonormal subspaces. The computed + // inverse supplies a floating-point estimate, not an interval proof. + return {1.0 / (norm * inverseInfinity), (1.0 / std::sqrt(inverseOne)) / std::sqrt(inverseInfinity)}; + } + + const Operation *m_operation; + const Inverse *m_inverse; + MPI_Comm m_communicator; + CoupledShapeCoarseOptions m_options; + std::uint64_t m_generation{0}; + bool m_ready{false}; + std::vector m_sources, m_images; + std::vector m_reports; + std::array m_upper{}; + mutable mfem::Vector m_residualWork, m_actionWork; + mutable CoupledShapeCoarseStatistics m_statistics; + }; +} // namespace experiment::plateau_diagnostics diff --git a/tests/solver/plateau_coupled_coarse.cpp b/tests/solver/plateau_coupled_coarse.cpp new file mode 100644 index 0000000..f8eb02a --- /dev/null +++ b/tests/solver/plateau_coupled_coarse.cpp @@ -0,0 +1,392 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +import mean_field; + +#include "../../experiments/plateau_coupled_coarse.hpp" + +namespace coarse_test { + namespace probe = ::experiment::plateau_diagnostics; + using namespace mean_field; + + [[nodiscard]] mfem::Vector Vector(std::initializer_list values) { + mfem::Vector result(static_cast(values.size())); + int index = 0; + for (double value : values) + result(index++) = value; + return result; + } + + [[nodiscard]] mfem::DenseMatrix Identity(int size) { + mfem::DenseMatrix result(size); + result = 0.0; + for (int i = 0; i < size; ++i) + result(i, i) = 1.0; + return result; + } + + [[nodiscard]] double Norm(const mfem::Vector &value) { + double local = value * value, global = 0.0; + MPI_Allreduce(&local, &global, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + return std::sqrt(global); + } + + [[nodiscard]] double Difference( + const mfem::Vector &left, + const mfem::Vector &right + ) { + mfem::Vector difference(left); + difference -= right; + return Norm(difference); + } + + class Operation final : public mfem::Operator { + public: + explicit Operation(mfem::DenseMatrix matrix) + : mfem::Operator( + matrix.Height(), + matrix.Width() + ), + matrix(std::move(matrix)), + factors(Width()) { + factors = 1.0; + } + void Mult( + const mfem::Vector &value, + mfem::Vector &action + ) const override { + ++applications; + matrix.Mult(value, action); + } + void NormalizeState( + const mfem::Vector &physical, + mfem::Vector &normalized + ) const { + ++normalizations; + for (int i = 0; i < Width(); ++i) + normalized(i) = physical(i) * factors(i); + } + bool IsPrepared() const { + return true; + } + const Operation &GetProblem() const { + return *this; + } + std::uint64_t GetPreparationGeneration() const { + return generation; + } + void Prepare() { + ++generation; + } + mfem::DenseMatrix matrix; + mfem::Vector factors; + mutable std::uint64_t applications{0}, normalizations{0}; + std::uint64_t generation{1}; + }; + + class Inverse : public mfem::Solver { + public: + static constexpr auto applicationContract = preconditioning::ApplicationContract::stationary_linear; + explicit Inverse(const Operation &operation) + : mfem::Solver(operation.Width()), + matrix(Identity(operation.Width())), + operation(operation) { + Refresh(); + } + void SetOperator(const mfem::Operator &) override { + } + void Mult( + const mfem::Vector &rhs, + mfem::Vector &action + ) const override { + ++applications; + matrix.Mult(rhs, action); + if (produceNonfinite && action.Size()) + action(0) = std::numeric_limits::infinity(); + } + void Refresh() { + generation = operation.generation; + } + bool IsCurrent() const { + return generation == operation.generation; + } + mfem::DenseMatrix matrix; + const Operation &operation; + bool produceNonfinite{false}; + mutable std::uint64_t applications{0}; + std::uint64_t generation{0}; + }; + + class FlexibleInverse final : public Inverse { + public: + using Inverse::Inverse; + static constexpr auto applicationContract = preconditioning::ApplicationContract::flexible; + }; + + static_assert(probe::CoupledCoarseNormalizedOperator); + static_assert(probe::CoupledCoarseStationaryInverse); + static_assert(!probe::CoupledCoarseStationaryInverse); + static_assert(!probe::CoupledCoarseNormalizedOperator); + + using OrdinaryModel = model::StellarModel>; + using CentralModel = model::StellarModel>; + template using Problem = equilibrium::StellarEquilibriumProblem; + template + using NormalizedOperation = normalization::NormalizedStellarEquilibriumOperator>; + template + using PhysicalInverse = + preconditioning::PreparedStellarInverseType>; + template + using NormalizedInverse = normalization::NormalizedStellarPreconditioner, PhysicalInverse>; + template + using Coarse = probe::PreparedCoupledShapeCoarse, NormalizedInverse>; + + // The experiment borrows the existing radial geometry rather than inventing + // a second nodal-coordinate map. Instantiation below checks actual solver + // types with optional generated borders, without a costly stellar fixture. + template + concept BorrowsRadialSurface = requires(const Problem &problem) { + { + problem.GetPhysicalOperator().GetDomainDeformation().surfaceDeformationPrescription().radialDirection(0, 0) + } -> std::same_as; + { + problem.GetPhysicalOperator().GetDomainDeformation().surfaceDeformationPrescription().referenceRadius(0) + } -> std::same_as; + requires std::is_lvalue_reference_v< + decltype(problem.GetPhysicalOperator().GetDomainDeformation().surfaceDeformationPrescription())>; + }; + static_assert(BorrowsRadialSurface); + static_assert(BorrowsRadialSurface); + static_assert(solver::LinearPreconditionerApplicationContractAvailable>); + static_assert(solver::LinearPreconditionerApplicationContractAvailable>); +} // namespace coarse_test + +template class experiment::plateau_diagnostics::PreparedCoupledShapeCoarse< + coarse_test::NormalizedOperation, + coarse_test::NormalizedInverse>; +template class experiment::plateau_diagnostics::PreparedCoupledShapeCoarse< + coarse_test::NormalizedOperation, + coarse_test::NormalizedInverse>; + +template std::vector experiment::plateau_diagnostics::MakeCoupledShapeModes( + const coarse_test::Problem &, + const coarse_test::PhysicalInverse & +); +template std::vector experiment::plateau_diagnostics::MakeCoupledShapeModes( + const coarse_test::Problem &, + const coarse_test::PhysicalInverse & +); + +TEST_CASE( + "Filtered shape coarse inverse is exact on its source space and linear", + "[solver][coupled_shape_coarse]" +) { + using namespace coarse_test; + mfem::DenseMatrix matrix(3); + matrix = 0.0; + matrix(0, 0) = 3.0; + matrix(0, 1) = 0.7; + matrix(0, 2) = 0.2; + matrix(1, 0) = -0.2; + matrix(1, 1) = 2.0; + matrix(1, 2) = 0.4; + matrix(2, 0) = 0.5; + matrix(2, 2) = 1.3; + Operation operation(matrix); + operation.factors = Vector({0.2, 3.0, 10.0}); + Inverse inverse(operation); + inverse.matrix(0, 0) = 0.4; + inverse.matrix(1, 1) = 0.5; + inverse.matrix(2, 2) = 0.8; + std::vector modes{Vector({1.0, 0.2, -0.1}), Vector({0.3, 1.0, 0.4}), Vector({1.0, 0.2, -0.1})}; + probe::PreparedCoupledShapeCoarse coarse(operation, inverse, std::span(modes), MPI_COMM_WORLD); + REQUIRE(coarse.UsesCoarseCorrection()); + REQUIRE(coarse.GetSourceBasis().size() == 2); + CHECK(coarse.GetStatistics().activeModes == 2); + CHECK(coarse.GetStatistics().sourceDependentModes == 1); + CHECK(coarse.GetStatistics().setupOperatorApplications == 5); + CHECK(coarse.GetStatistics().setupInverseApplications == 2); + CHECK(operation.normalizations == 3); + CHECK(coarse.GetStatistics().persistentVectorBytes == 18 * sizeof(mfem::real_t)); + + mfem::Vector value(3), rhs(3), actual(3), expected(3); + operation.NormalizeState(modes[0], value); + operation.Mult(value, rhs); + rhs /= Norm(rhs); + CHECK(Difference(coarse.GetSourceBasis()[0], rhs) < 2.e-14); + for (const auto &source : coarse.GetSourceBasis()) { + inverse.Mult(source, expected); + operation.Mult(expected, rhs); + const auto before = operation.applications; + coarse.Mult(rhs, actual); + CHECK(operation.applications == before); + CHECK(Difference(actual, expected) < 2.e-12 * std::max(1.0, Norm(expected))); + } + const auto before = operation.applications; + const auto inverseBefore = inverse.applications; + rhs = 0.0; + actual = 17.0; + coarse.Mult(rhs, actual); + CHECK(Norm(actual) == 0.0); + CHECK(operation.applications == before); + CHECK(inverse.applications == inverseBefore + 1); + + const auto first = Vector({0.2, -0.8, 1.3}), second = Vector({-0.4, 0.1, 0.7}); + mfem::Vector firstAction(3), secondAction(3), combination(first); + combination *= 0.7; + combination.Add(-1.2, second); + coarse.Mult(first, firstAction); + coarse.Mult(second, secondAction); + coarse.Mult(combination, actual); + expected = firstAction; + expected *= 0.7; + expected.Add(-1.2, secondAction); + CHECK(Difference(actual, expected) < 2.e-12); + coarse.Mult(combination, combination); + CHECK(Difference(combination, actual) < 2.e-14); + CHECK(operation.applications == before); +} + +TEST_CASE( + "Filtered coarse overlap rejects a near-right-angle even at rank one", + "[solver][coupled_shape_coarse]" +) { + using namespace coarse_test; + for (const double cosine : {0.25, 1.e-12, 0.0}) { + const double sine = std::sqrt(1.0 - cosine * cosine); + mfem::DenseMatrix rotation(2); + rotation(0, 0) = rotation(1, 1) = cosine; + rotation(0, 1) = -sine; + rotation(1, 0) = sine; + Operation operation(rotation); + Inverse inverse(operation); + std::array modes{Vector({1.0, 0.0})}; + probe::PreparedCoupledShapeCoarse coarse( + operation, inverse, std::span(modes), MPI_COMM_WORLD + ); + CHECK(coarse.GetStatistics().independentModes == 1); + if (cosine == 0.25) { + CHECK(coarse.UsesCoarseCorrection()); + CHECK(coarse.GetStatistics().minimumOverlapLowerBound == Catch::Approx(cosine).epsilon(1.e-12)); + } else { + CHECK_FALSE(coarse.UsesCoarseCorrection()); + CHECK(coarse.GetStatistics().fallback == probe::CoupledShapeCoarseFallback::insufficient_subspace_overlap); + if (cosine != 0.0) { + CHECK(coarse.GetStatistics().reciprocalOverlapCondition == Catch::Approx(1.0)); + CHECK(coarse.GetStatistics().minimumOverlapLowerBound < 1.e-8); + } + const auto rhs = Vector({0.3, 0.7}); + mfem::Vector result(2); + const auto before = operation.applications; + coarse.Mult(rhs, result); + CHECK(Difference(result, rhs) == 0.0); + CHECK(operation.applications == before); + CHECK(coarse.GetStatistics().inverseApplications == 1); + } + } +} + +TEST_CASE( + "Coarse setup drops a null filtered image without altering equations", + "[solver][coupled_shape_coarse]" +) { + using namespace coarse_test; + mfem::DenseMatrix matrix(2); + matrix = 0.0; + matrix(0, 1) = 1.0; + Operation operation(matrix); + Inverse inverse(operation); + std::array modes{Vector({0.0, 1.0})}; + probe::PreparedCoupledShapeCoarse coarse(operation, inverse, std::span(modes), MPI_COMM_WORLD); + CHECK_FALSE(coarse.UsesCoarseCorrection()); + CHECK(coarse.GetStatistics().imageDependentModes == 1); + CHECK(coarse.GetStatistics().fallback == probe::CoupledShapeCoarseFallback::no_independent_modes); + CHECK(operation.matrix(0, 0) == 0.0); + CHECK(operation.matrix(0, 1) == 1.0); + CHECK(operation.matrix(1, 0) == 0.0); + CHECK(operation.matrix(1, 1) == 0.0); +} + +TEST_CASE( + "Coarse inverse rejects stale preparation and preserves failed outputs", + "[solver][coupled_shape_coarse]" +) { + using namespace coarse_test; + Operation operation(Identity(2)); + Inverse inverse(operation); + std::array modes{Vector({1.0, 0.5})}; + probe::PreparedCoupledShapeCoarse coarse(operation, inverse, std::span(modes), MPI_COMM_WORLD); + const auto rhs = Vector({0.3, 0.2}); + mfem::Vector result(2); + result = 9.0; + operation.Prepare(); + CHECK_FALSE(coarse.IsCurrent()); + CHECK_THROWS_AS(coarse.Mult(rhs, result), std::logic_error); + CHECK(result(0) == 9.0); + CHECK(result(1) == 9.0); + inverse.Refresh(); + CHECK_FALSE(coarse.IsCurrent()); + coarse.Refresh(modes); + CHECK(coarse.IsCurrent()); + int rank = 0; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + inverse.produceNonfinite = rank == 0; + CHECK_THROWS_AS(coarse.Mult(rhs, result), std::runtime_error); + CHECK(result(0) == 9.0); + CHECK(result(1) == 9.0); + inverse.produceNonfinite = false; + coarse.Mult(rhs, result); + CHECK(Difference(result, rhs) < 1.e-14); + CHECK(coarse.GetStatistics().failedApplications == 2); + CHECK(coarse.GetStatistics().setups == 2); +} + +TEST_CASE( + "Coarse preparation supports locally empty ranks and rejects a globally empty system", + "[solver][coupled_shape_coarse]" +) { + using namespace coarse_test; + int rank = 0; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + const int localSize = rank == 0 ? 2 : 0; + Operation operation(Identity(localSize)); + Inverse inverse(operation); + std::array modes{mfem::Vector(localSize)}; + modes[0] = 1.0; + probe::PreparedCoupledShapeCoarse coarse(operation, inverse, std::span(modes), MPI_COMM_WORLD); + REQUIRE(coarse.UsesCoarseCorrection()); + mfem::Vector result(localSize); + coarse.Mult(modes[0], result); + CHECK(Difference(result, modes[0]) < 1.e-14); + Operation empty(Identity(0)); + Inverse emptyInverse(empty); + std::array emptyMode{mfem::Vector(0)}; + CHECK_THROWS_AS( + (probe::PreparedCoupledShapeCoarse( + empty, emptyInverse, std::span(emptyMode), MPI_COMM_WORLD + )), + std::invalid_argument + ); +}