diff --git a/CMakeLists.txt b/CMakeLists.txt index 9bc0c81..4254996 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -389,6 +389,7 @@ add_executable(tests tests/solver/stellar_equilibrium_runtime.cpp tests/solver/enthalpy_admissibility.cpp tests/solver/plateau_coupled_coarse.cpp + tests/solver/plateau_surface_border_schur.cpp tests/solver/gravity_subsystem.cpp tests/solver/gravity_compact_applications.cpp tests/solver/gravity_reduction.cpp @@ -417,6 +418,10 @@ add_executable(continuation_contract_tests EXCLUDE_FROM_ALL tests/solver/stellar_numerical_strategies.cpp tests/solver/continuation_snapshot.cpp) target_link_libraries(continuation_contract_tests PRIVATE mean_field test_mod Catch2::Catch2 Boost::boost) +add_executable(surface_border_schur_contract_tests EXCLUDE_FROM_ALL + tests/test_main.cpp + tests/solver/plateau_surface_border_schur.cpp) +target_link_libraries(surface_border_schur_contract_tests PRIVATE mean_field test_mod Catch2::Catch2 Boost::boost) # A deliberately opt-in API workbench. It is compiled explicitly during API # verification but temporary user experiments do not break the default build. diff --git a/experiments/plateau_preconditioner_benchmark.hpp b/experiments/plateau_preconditioner_benchmark.hpp index 33b1ceb..4d3ff88 100644 --- a/experiments/plateau_preconditioner_benchmark.hpp +++ b/experiments/plateau_preconditioner_benchmark.hpp @@ -14,6 +14,7 @@ #include #include "plateau_coupled_coarse.hpp" +#include "plateau_surface_border_schur.hpp" #include #include #include @@ -31,7 +32,9 @@ enum class PreconditionerBenchmarkVariant { gravity_schur, stellar_ldu_gravity_schur, coupled_shape_coarse, - gravity_schur_coarse + gravity_schur_coarse, + surface_border_schur, + gravity_schur_surface_border_schur }; inline PreconditionerBenchmarkVariant ParsePreconditionerBenchmarkVariant(const std::string &name) { @@ -44,6 +47,9 @@ inline PreconditionerBenchmarkVariant ParsePreconditionerBenchmarkVariant(const if (name == "stellar-ldu-gravity-schur") return PreconditionerBenchmarkVariant::stellar_ldu_gravity_schur; if (name == "coupled-shape-coarse") return PreconditionerBenchmarkVariant::coupled_shape_coarse; if (name == "gravity-schur-coarse") return PreconditionerBenchmarkVariant::gravity_schur_coarse; + if (name == "surface-border-schur") return PreconditionerBenchmarkVariant::surface_border_schur; + if (name == "gravity-schur-surface-border-schur") + return PreconditionerBenchmarkVariant::gravity_schur_surface_border_schur; throw std::invalid_argument("Unknown preconditioner benchmark variant: " + name); } @@ -271,6 +277,12 @@ void BenchmarkPreconditioners(State &state, const std::filesystem::path &output, "true_check_jacobian_applications,krylov_other_seconds\n"; auto coarseReport = detail::Open(output / "preconditioner_coarse.csv"); coarseReport << "role,variant,stage,metric,value\n"; + auto schurReport = detail::Open(output / "preconditioner_surface_border_schur.csv"); + schurReport << "role,variant,stage,metric,value\n"; + auto schurCoordinates = detail::Open(output / "preconditioner_surface_border_coordinates.csv"); + schurCoordinates << "role,variant,coarse_index,residual_index,correction_index\n"; + auto schurMatrices = detail::Open(output / "preconditioner_surface_border_matrices.csv"); + schurMatrices << "role,variant,row,column,uneliminated,schur\n"; auto work = detail::Open(output / "preconditioner_work.csv"); work << "role,variant,stage,metric,count\n"; auto actions = detail::Open(output / "preconditioner_actions.csv"); @@ -285,7 +297,7 @@ void BenchmarkPreconditioners(State &state, const std::filesystem::path &output, history << "role,variant,iteration,reported_norm,final\n"; const auto run = [&]( + bool RefinedGravity = false, bool UseCoarse = false, bool UseSurfaceBorderSchur = false>( const std::string &role, const std::string &variant) { const auto candidateStart = detail::Clock::now(); auto physical = detail::Prepare(state.Problem(), options); @@ -319,8 +331,38 @@ void BenchmarkPreconditioners(State &state, const std::filesystem::path &output, const double bindingSeconds = detail::Seconds(sectionStart); using NormalizedOperation = std::remove_reference_t; using Coarse = PreparedCoupledShapeCoarse; + using UnborderedStructure = NormalizedUnborderedStructureInverse; + using SurfaceBorderSchur = PreparedSurfaceBorderSchur; std::unique_ptr coarse; + std::unique_ptr unborderedStructure; + std::unique_ptr surfaceBorderSchur; + static_assert(!(UseCoarse && UseSurfaceBorderSchur)); mfem::Solver *applicationInverse = &normalized; + const auto recordSurfaceBorderSchur = [&](const char *stage) { + if (!surfaceBorderSchur) return; + const auto &stats = surfaceBorderSchur->GetStatistics(); + const auto metric = [&](const char *name, auto value) { + schurReport << role << ',' << variant << ',' << stage << ',' << name << ',' << value << '\n'; + }; + metric("status", static_cast(stats.status)); + metric("coarse_size", stats.coarseSize); + metric("reciprocal_condition", stats.reciprocalCondition); + metric("minimum_scaled_pivot", stats.minimumScaledPivot); + metric("setup_jacobian_applications", stats.setupOperatorApplications); + metric("setup_inverse_applications", stats.setupInverseApplications); + metric("setup_seconds", stats.lastSetupSeconds); + metric("setup_jacobian_seconds", stats.setupOperatorSeconds); + metric("setup_inverse_seconds", stats.setupInverseSeconds); + metric("persistent_vector_bytes", stats.persistentVectorBytes); + metric("applications", stats.applications); + metric("internal_jacobian_applications", stats.operatorApplications); + metric("internal_inverse_applications", stats.inverseApplications); + metric("internal_jacobian_seconds", stats.operatorSeconds); + metric("internal_inverse_seconds", stats.inverseSeconds); + metric("application_seconds", stats.applicationSeconds); + metric("failed_applications", stats.failedApplications); + schurReport.flush(); + }; const auto recordCoarse = [&](const char *stage) { if (!coarse) return; const auto &stats = coarse->GetStatistics(); @@ -398,6 +440,34 @@ void BenchmarkPreconditioners(State &state, const std::filesystem::path &output, throw std::runtime_error("The actual FEM coarse inverse failed its exact-subspace check."); } + if constexpr (UseSurfaceBorderSchur) { + sectionStart = detail::Clock::now(); + unborderedStructure = std::make_unique(*state.normalizedOperator, physical); + surfaceBorderSchur = std::make_unique( + *state.normalizedOperator, *unborderedStructure, MakeSurfaceBorderSchurCoordinates(physical), + state.Problem().GetCommunicator()); + const double extraConstructionSeconds = detail::Seconds(sectionStart); + constructionSeconds += extraConstructionSeconds; + const auto &coordinates = surfaceBorderSchur->GetCoordinates(); + const auto &schur = surfaceBorderSchur->GetSchurMatrix(); + const auto &uneliminated = surfaceBorderSchur->GetUneliminatedCoarseMatrix(); + for (int row = 0; row < schur.Height(); ++row) { + schurCoordinates << role << ',' << variant << ',' << row << ',' + << coordinates.residualIndices[row] << ',' + << coordinates.correctionIndices[row] << '\n'; + for (int column = 0; column < schur.Width(); ++column) + schurMatrices << role << ',' << variant << ',' << row << ',' << column << ',' + << uneliminated(row, column) << ',' << schur(row, column) << '\n'; + } + schurCoordinates.flush(); + schurMatrices.flush(); + recordSurfaceBorderSchur("constructed"); + recordWork("after_surface_border_setup"); + if (!surfaceBorderSchur->IsCurrent()) + throw std::runtime_error("The measured surface/scalar Schur matrix is singular or ill-conditioned; no fallback was applied."); + applicationInverse = surfaceBorderSchur.get(); + } + mfem::Vector inverseAction(rhs.Size()), inverseJacobianAction(rhs.Size()); sectionStart = detail::Clock::now(); applicationInverse->Mult(rhs, inverseAction); @@ -412,11 +482,23 @@ void BenchmarkPreconditioners(State &state, const std::filesystem::path &output, if (!detail::Finite(inverseJacobianAction) || !std::isfinite(actionNorm)) throw std::runtime_error("The diagnostic inverse Jacobian image or its norm is nonfinite."); inverseJacobianAction -= rhs; + if (surfaceBorderSchur) { + double squared = 0.0; + for (int row : surfaceBorderSchur->GetCoordinates().residualIndices) + squared += inverseJacobianAction(row) * inverseJacobianAction(row); + const double coarseRelative = std::sqrt(squared) / rhsNorm; + schurReport << role << ',' << variant + << ",validation,coarse_equation_defect_relative_to_full_rhs," << coarseRelative << '\n'; + schurReport.flush(); + if (!std::isfinite(coarseRelative) || coarseRelative > 1.e-8) + throw std::runtime_error("The Schur inverse failed the authoritative-J coarse-equation check."); + } mfem::Vector referenceDifference(inverseAction); referenceDifference -= referenceCorrection; const double referenceRelative = detail::Relative(referenceDifference.Norml2(), std::max(inverseAction.Norml2(), referenceCorrection.Norml2())); const bool requireEquivalent = std::same_as && !MaterialLDU && !UseCoarse && + !UseSurfaceBorderSchur && (!RefinedGravity || options.schurCorrectionSweeps == 0); actions << role << ',' << variant << ',' << inverseAction.Norml2() << ',' << actionNorm << ',' << inverseJacobianAction.Norml2() << ',' << inverseJacobianAction.Norml2() / rhsNorm << ',' @@ -447,6 +529,7 @@ void BenchmarkPreconditioners(State &state, const std::filesystem::path &output, const double krylovSetupSeconds = detail::Seconds(sectionStart); recordWork("before_krylov"); recordCoarse("before_krylov"); + recordSurfaceBorderSchur("before_krylov"); countedJacobian.ResetStatistics(); countedInverse.ResetStatistics(); sectionStart = detail::Clock::now(); @@ -454,6 +537,7 @@ void BenchmarkPreconditioners(State &state, const std::filesystem::path &output, const double krylovSeconds = detail::Seconds(sectionStart); recordWork("after_krylov"); recordCoarse("after_krylov"); + recordSurfaceBorderSchur("after_krylov"); const auto jacobianStats = countedJacobian.GetStatistics(); const auto inverseStats = countedInverse.GetStatistics(); const auto lifecycle = countedInverse.GetLifecycleStatistics(); @@ -537,6 +621,10 @@ void BenchmarkPreconditioners(State &state, const std::filesystem::path &output, run.template operator()("candidate", options.variant); break; case PreconditionerBenchmarkVariant::gravity_schur_coarse: run.template operator()("candidate", options.variant); break; + case PreconditionerBenchmarkVariant::surface_border_schur: + run.template operator()("candidate", options.variant); break; + case PreconditionerBenchmarkVariant::gravity_schur_surface_border_schur: + run.template operator()("candidate", options.variant); break; } auto scope = detail::Open(output / "preconditioner_scope.txt"); scope << "Both solves use fresh physical inverses, separately completed scalar borders, the same physical J,\n" @@ -549,6 +637,12 @@ void BenchmarkPreconditioners(State &state, const std::filesystem::path &output, << "unchanged_state_refresh_seconds measures an immediate Refresh; it is not a changed-state refresh benchmark.\n" << "Each candidate receives one diagnostic inverse/Jv application before the timed cold Krylov solve.\n" << "Application timers/counts in the summary include only Krylov work; the independent true check is separate.\n" + << "Surface/scalar Schur setup is included in construction; its internal J applications are separately counted.\n" + << "For that inverse, outer Krylov J counts exclude internal Schur J work; inverse seconds include it.\n" + << "Schur status codes are 0=unprepared, 1=ready, 2=singular, 3=ill-conditioned; rejected setup has no fallback.\n" + << "The Schur partition uses every native surface parameter and generated scalar block, including a central border if present.\n" + << "Its field inverse restricts the unbordered structure inverse, not an exact field elimination.\n" + << "Existing scalar-border construction is retained in setup cost but its inverse is not applied by the new wrapper.\n" << "setup_plus_solve_seconds includes construction, normalization binding, Krylov setup, and Krylov solve.\n" << "setup_solve_and_true_check_seconds adds the separately measured original-J true-residual gate.\n" << "measured_candidate_seconds also includes diagnostic refresh, coarse exactness checks, and inverse/Jv comparison.\n" diff --git a/experiments/plateau_surface_border_schur.hpp b/experiments/plateau_surface_border_schur.hpp new file mode 100644 index 0000000..e83fb73 --- /dev/null +++ b/experiments/plateau_surface_border_schur.hpp @@ -0,0 +1,432 @@ +#pragma once + +// Include after import mean_field. The experiment changes only the inverse +// supplied to Krylov; every column below comes from the authoritative Jacobian. +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +namespace experiment::plateau_diagnostics { + struct SurfaceBorderSchurCoordinates final { + // Equation and unknown ordering need not coincide in a generated form. + std::vector residualIndices; + std::vector correctionIndices; + }; + + template + [[nodiscard]] SurfaceBorderSchurCoordinates MakeSurfaceBorderSchurCoordinates(const PhysicalInverse &inverse) { + using namespace mean_field; + using namespace utils::blocks; + using Root = std::remove_cvref_t; + using Structure = typename Root::StructureComponent; + constexpr auto shapeColumn = type_index_v; + constexpr auto shapeRow = type_index_v; + static_assert(Root::borderValueArity == Root::borderResidualArity); + static_assert(shapeColumn < Structure::CorrectionBlocks::size); + static_assert(shapeRow < Structure::ResidualBlocks::size); + const auto &map = inverse.GetCoordinateMap(); + const auto &columns = map.GetCorrectionRanges(); + const auto &rows = map.GetResidualRanges(); + SurfaceBorderSchurCoordinates result; + const auto append = [](const auto &range, std::vector &indices) { + for (int i = 0; i < range.size; ++i) + indices.push_back(range.equilibriumOffset + i); + }; + append(columns[shapeColumn], result.correctionIndices); + append(rows[shapeRow], result.residualIndices); + // Generated scalar blocks follow the complete structure descriptor. + // This includes the optional central phase condition without naming a + // particular model or interpreting a runtime stable-ID string. + for (std::size_t i = Structure::CorrectionBlocks::size; i < columns.size(); ++i) + append(columns[i], result.correctionIndices); + for (std::size_t i = Structure::ResidualBlocks::size; i < rows.size(); ++i) + append(rows[i], result.residualIndices); + return result; + } + + template + concept SurfaceBorderSchurOperation = std::derived_from && + requires(const Candidate &operation) { + { operation.IsPrepared() } -> std::same_as; + { operation.GetProblem().GetPreparationGeneration() } -> std::convertible_to; + }; + + template + concept SurfaceBorderSchurInverse = 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; + }; + + // Reusing the complete scalar-bordered inverse here would impose zero + // scalar-constraint RHS before elimination. Its field restriction can be + // singular even when A_FF is invertible (a zero scalar diagonal suffices). + // Borrow only the original structure inverse, preserving its typed packing + // and the exact physical/normalized coordinate maps. + template + class NormalizedUnborderedStructureInverse final : public mfem::Solver { + using Root = std::remove_cvref_t().GetBlock())>; + using Backend = typename Root::StructureComponent::BackendType; + public: + static constexpr auto applicationContract = mean_field::preconditioning::backend::applicationContract; + static_assert(applicationContract == mean_field::preconditioning::ApplicationContract::stationary_linear); + + NormalizedUnborderedStructureInverse(const Operation &operation, const PhysicalInverse &inverse) + : mfem::Solver(operation.Width(), operation.Height()), m_operation(&operation), m_inverse(&inverse), + m_generation(operation.GetProblem().GetPreparationGeneration()), m_physicalRhs(Width()), + m_physicalAction(Height()), m_groupedRhs(Width()), m_groupedAction(Height()) { + if (&operation.GetProblem() != &inverse.GetProblem() || inverse.Height() != Height() || + inverse.Width() != Width() || !IsCurrent()) + throw std::invalid_argument("The structure adapter requires the current inverse of the same problem."); + } + void SetOperator(const mfem::Operator &operation) override { + if (&operation != m_operation || !IsCurrent()) + throw std::invalid_argument("The structure adapter cannot bind another or stale Jacobian."); + } + [[nodiscard]] bool IsCurrent() const { + return m_operation->IsPrepared() && m_inverse->IsCurrent() && !m_inverse->iterative_mode && + !m_inverse->GetGroupedPreconditioner().GetStructurePreconditioner().iterative_mode && + m_generation == m_operation->GetProblem().GetPreparationGeneration(); + } + [[nodiscard]] const Operation &GetNormalizedJacobian() const noexcept { return *m_operation; } + void Mult(const mfem::Vector &rhs, mfem::Vector &action) const override { + if (!IsCurrent() || rhs.Size() != Width() || action.Size() != Height()) + throw std::logic_error("The structure adapter requires current objects and compatible vectors."); + const auto &structure = m_inverse->GetGroupedPreconditioner().GetStructurePreconditioner(); + const auto &map = m_inverse->GetCoordinateMap(); + m_operation->DenormalizeResidual(rhs, m_physicalRhs); + map.PackResidual(m_physicalRhs, m_groupedRhs); + m_groupedAction = 0.0; + const mfem::Vector fieldRhs(m_groupedRhs, 0, structure.Width()); + mfem::Vector fieldAction(m_groupedAction, 0, structure.Height()); + structure.Mult(fieldRhs, fieldAction); + fieldAction.SyncAliasMemory(m_groupedAction); + map.UnpackCorrection(m_groupedAction, m_physicalAction); + m_operation->NormalizeState(m_physicalAction, action); + } + private: + const Operation *m_operation; + const PhysicalInverse *m_inverse; + std::uint64_t m_generation; + mutable mfem::Vector m_physicalRhs, m_physicalAction, m_groupedRhs, m_groupedAction; + }; + + enum class SurfaceBorderSchurStatus { unprepared, ready, singular, ill_conditioned }; + + struct SurfaceBorderSchurOptions final { + double minimumReciprocalCondition{1.e-14}; + }; + + struct SurfaceBorderSchurStatistics final { + SurfaceBorderSchurStatus status{SurfaceBorderSchurStatus::unprepared}; + int coarseSize{0}; + std::uint64_t setups{0}, setupOperatorApplications{0}, setupInverseApplications{0}; + double reciprocalCondition{0.0}; + double minimumScaledPivot{0.0}; + double lastSetupSeconds{0.0}, setupOperatorSeconds{0.0}, setupInverseSeconds{0.0}; + std::uint64_t applications{0}, operatorApplications{0}, inverseApplications{0}, failedApplications{0}; + double applicationSeconds{0.0}, operatorSeconds{0.0}, inverseSeconds{0.0}; + std::size_t persistentVectorBytes{0}; + }; + + /* With C the complete surface/scalar block and F every remaining field, + * B = R_F P0 E_F approximates A_FF^{-1}. Store + * + * Z = E_C - E_F B A_FC, S = R_C A Z. + * + * Application gives y = E_F B r_F and x = y + Z S^{-1}(r_C-R_C A y). + * The coarse equations are satisfied exactly (up to arithmetic), even + * when B is approximate. The complete inverse is exact only when B is. + * Restricting the existing inverse defines B; its internal border response + * is not claimed to be an exact elimination of the frozen field block. + * No symmetry, positive definiteness, or reduced physical equation is + * assumed. A singular/unsafe S is reported and cannot be applied. + */ + template + class PreparedSurfaceBorderSchur final : public mfem::Solver { + public: + static constexpr auto applicationContract = mean_field::preconditioning::ApplicationContract::stationary_linear; + static constexpr std::size_t maximumCoarseSize = 128; + + PreparedSurfaceBorderSchur(const Operation &operation, const Inverse &inverse, + SurfaceBorderSchurCoordinates coordinates, MPI_Comm communicator, + SurfaceBorderSchurOptions options = {}) + : mfem::Solver(operation.Width(), operation.Height()), m_operation(&operation), m_inverse(&inverse), + m_coordinates(std::move(coordinates)), m_communicator(communicator), m_options(options), + m_rhs(operation.Height()), m_fieldAction(operation.Width()), m_image(operation.Height()) { + Validate(); + Refresh(); + } + + PreparedSurfaceBorderSchur(const PreparedSurfaceBorderSchur &) = delete; + PreparedSurfaceBorderSchur &operator=(const PreparedSurfaceBorderSchur &) = delete; + PreparedSurfaceBorderSchur(PreparedSurfaceBorderSchur &&) = delete; + PreparedSurfaceBorderSchur &operator=(PreparedSurfaceBorderSchur &&) = delete; + + void SetOperator(const mfem::Operator &operation) override { + if (&operation != m_operation) + throw std::invalid_argument("The surface-border Schur inverse cannot bind a different Jacobian."); + RequireCurrent(); + } + + [[nodiscard]] bool IsCurrent() const { + return m_statistics.status == SurfaceBorderSchurStatus::ready && DependenciesCurrent(); + } + [[nodiscard]] const SurfaceBorderSchurStatistics &GetStatistics() const noexcept { return m_statistics; } + [[nodiscard]] const SurfaceBorderSchurCoordinates &GetCoordinates() const noexcept { return m_coordinates; } + [[nodiscard]] const mfem::DenseMatrix &GetSchurMatrix() const noexcept { return m_schur; } + [[nodiscard]] const mfem::DenseMatrix &GetUneliminatedCoarseMatrix() const noexcept { return m_coarse; } + [[nodiscard]] std::span GetCompletedColumns() const noexcept { return m_columns; } + + void Refresh() { + const double start = MPI_Wtime(); + Validate(); + m_statistics.status = SurfaceBorderSchurStatus::unprepared; + m_generation = m_operation->GetProblem().GetPreparationGeneration(); + const int count = static_cast(m_coordinates.correctionIndices.size()); + m_statistics.coarseSize = count; + m_columns.clear(); + m_columns.reserve(count); + m_schur.SetSize(count); + m_coarse.SetSize(count); + m_solveWork.SetSize(count); + mfem::Vector basis(Height()), image(Width()), response(Height()); + for (int column = 0; column < count; ++column) { + basis = 0.0; + basis(m_coordinates.correctionIndices[column]) = 1.0; + ApplyOperation(basis, image, true); + for (int row = 0; row < count; ++row) + m_coarse(row, column) = image(m_coordinates.residualIndices[row]); + ApplyFields(image, response, true); + response *= -1.0; + response += basis; + m_columns.emplace_back(response); + ApplyOperation(response, image, true); + for (int row = 0; row < count; ++row) + m_schur(row, column) = image(m_coordinates.residualIndices[row]); + } + Factor(); + if (!DependenciesCurrent()) + throw std::logic_error("The surface-border Schur dependencies changed during preparation."); + m_statistics.lastSetupSeconds = MPI_Wtime() - start; + m_statistics.persistentVectorBytes = sizeof(mfem::real_t) * + (static_cast(count) * Height() + 3 * static_cast(Height()) + count); + ++m_statistics.setups; + } + + void Mult(const mfem::Vector &rightHandSide, mfem::Vector &action) const override { + const double start = MPI_Wtime(); + ++m_statistics.applications; + try { + RequireCurrent(); + if (rightHandSide.Size() != Width() || action.Size() != Height()) + throw std::invalid_argument("The surface-border Schur inverse requires preallocated root vectors."); + RequireFinite(rightHandSide); + ApplyFields(rightHandSide, m_fieldAction, false); + ApplyOperation(m_fieldAction, m_image, false); + for (int row = 0; row < m_solveWork.Size(); ++row) + m_solveWork(row) = rightHandSide(m_coordinates.residualIndices[row]) - + m_image(m_coordinates.residualIndices[row]); + SolveFactored(m_solveWork); + // No write to the caller's output occurs before all RHS uses, + // so in-place application has the same algebra as separate buffers. + action = m_fieldAction; + for (int column = 0; column < m_solveWork.Size(); ++column) + action.Add(m_solveWork(column), m_columns[column]); + RequireFinite(action); + RequireCurrent(); + m_statistics.applicationSeconds += MPI_Wtime() - start; + } catch (...) { + ++m_statistics.failedApplications; + m_statistics.applicationSeconds += MPI_Wtime() - start; + throw; + } + } + + private: + void Validate() const { + int ranks = 0; + if (m_communicator == MPI_COMM_NULL || MPI_Comm_size(m_communicator, &ranks) != MPI_SUCCESS || ranks != 1) + throw std::invalid_argument("The surface-border Schur experiment currently requires one MPI rank."); + if (!m_operation->IsPrepared() || !m_inverse->IsCurrent() || m_inverse->iterative_mode || + m_operation->Height() != m_operation->Width() || m_inverse->Height() != Height() || + m_inverse->Width() != Width()) + throw std::invalid_argument("The surface-border Schur experiment requires current compatible objects."); + if constexpr (requires { m_inverse->GetNormalizedJacobian(); }) { + if (&m_inverse->GetNormalizedJacobian() != m_operation) + throw std::invalid_argument("The inverse belongs to a different normalized Jacobian."); + } + if (!std::isfinite(m_options.minimumReciprocalCondition) || + m_options.minimumReciprocalCondition <= 0.0 || m_options.minimumReciprocalCondition >= 1.0) + throw std::invalid_argument("The Schur reciprocal-condition threshold must lie strictly between zero and one."); + if (m_coordinates.residualIndices.empty() || + m_coordinates.residualIndices.size() != m_coordinates.correctionIndices.size() || + m_coordinates.residualIndices.size() >= static_cast(Height())) + throw std::invalid_argument("The Schur partition requires nonempty fields and square coarse coordinates."); + if (m_coordinates.residualIndices.size() > maximumCoarseSize) + throw std::invalid_argument("The bounded surface-border experiment supports at most 128 coarse coordinates."); + const auto verify = [this](const auto &indices) { + std::vector used(Height(), false); + for (int index : indices) { + if (index < 0 || index >= Height() || used[index]) + throw std::invalid_argument("Schur coordinate indices must be distinct and in range."); + used[index] = true; + } + }; + verify(m_coordinates.residualIndices); + verify(m_coordinates.correctionIndices); + } + + [[nodiscard]] bool DependenciesCurrent() const { + return m_operation->IsPrepared() && m_inverse->IsCurrent() && !m_inverse->iterative_mode && + m_generation == m_operation->GetProblem().GetPreparationGeneration(); + } + void RequireCurrent() const { + if (!IsCurrent()) + throw std::logic_error("The surface-border Schur inverse is stale, singular, or ill-conditioned."); + } + static void RequireFinite(const mfem::Vector &vector) { + for (int i = 0; i < vector.Size(); ++i) + if (!std::isfinite(vector(i))) + throw std::runtime_error("The surface-border Schur experiment encountered a nonfinite vector."); + } + void ApplyOperation(const mfem::Vector &input, mfem::Vector &output, bool setup) const { + const double start = MPI_Wtime(); + m_operation->Mult(input, output); + const double seconds = MPI_Wtime() - start; + if (setup) { + ++m_statistics.setupOperatorApplications; + m_statistics.setupOperatorSeconds += seconds; + } else { + ++m_statistics.operatorApplications; + m_statistics.operatorSeconds += seconds; + } + RequireFinite(output); + } + void ApplyFields(const mfem::Vector &input, mfem::Vector &output, bool setup) const { + m_rhs = input; + for (int index : m_coordinates.residualIndices) + m_rhs(index) = 0.0; + const double start = MPI_Wtime(); + m_inverse->Mult(m_rhs, output); + const double seconds = MPI_Wtime() - start; + if (setup) { + ++m_statistics.setupInverseApplications; + m_statistics.setupInverseSeconds += seconds; + } else { + ++m_statistics.inverseApplications; + m_statistics.inverseSeconds += seconds; + } + RequireFinite(output); + for (int index : m_coordinates.correctionIndices) + output(index) = 0.0; + } + + void Factor() { + const int count = m_schur.Height(); + m_lu = m_schur; + m_pivots.resize(count); + m_scale = 0.0; + for (int i = 0; i < count * count; ++i) { + if (!std::isfinite(m_lu.Data()[i])) + throw std::runtime_error("The coarse Schur matrix is nonfinite."); + m_scale = std::max(m_scale, std::abs(m_lu.Data()[i])); + } + m_statistics.reciprocalCondition = 0.0; + m_statistics.minimumScaledPivot = 0.0; + m_statistics.status = SurfaceBorderSchurStatus::singular; + if (m_scale == 0.0) + return; + for (int i = 0; i < count * count; ++i) + m_lu.Data()[i] /= m_scale; + double norm = 0.0; + for (int row = 0; row < count; ++row) { + double sum = 0.0; + for (int column = 0; column < count; ++column) + sum += std::abs(m_lu(row, column)); + norm = std::max(norm, sum); + } + m_statistics.minimumScaledPivot = std::numeric_limits::infinity(); + for (int column = 0; column < count; ++column) { + int pivot = column; + for (int row = column + 1; row < count; ++row) + if (std::abs(m_lu(row, column)) > std::abs(m_lu(pivot, column))) + pivot = row; + m_pivots[column] = pivot; + const double magnitude = std::abs(m_lu(pivot, column)); + m_statistics.minimumScaledPivot = std::min(m_statistics.minimumScaledPivot, magnitude); + if (!std::isfinite(magnitude) || magnitude == 0.0) + return; + if (pivot != column) + for (int j = 0; j < count; ++j) + std::swap(m_lu(column, j), m_lu(pivot, j)); + for (int row = column + 1; row < count; ++row) { + m_lu(row, column) /= m_lu(column, column); + for (int j = column + 1; j < count; ++j) + m_lu(row, j) -= m_lu(row, column) * m_lu(column, j); + } + } + // A small explicit inverse supplies an infinity-norm condition + // estimate without regularizing a physically important coarse mode. + mfem::Vector column(count), sums(count); + sums = 0.0; + for (int j = 0; j < count; ++j) { + column = 0.0; + column(j) = m_scale; + SolveFactored(column); + for (int i = 0; i < count; ++i) + sums(i) += std::abs(column(i)); + } + const double inverseNorm = sums.Max(); + m_statistics.reciprocalCondition = 1.0 / (norm * inverseNorm); + m_statistics.status = std::isfinite(m_statistics.reciprocalCondition) && + m_statistics.reciprocalCondition >= m_options.minimumReciprocalCondition + ? SurfaceBorderSchurStatus::ready : SurfaceBorderSchurStatus::ill_conditioned; + } + + void SolveFactored(mfem::Vector &values) const { + const int count = values.Size(); + values /= m_scale; + for (int i = 0; i < count; ++i) { + if (m_pivots[i] != i) + std::swap(values(i), values(m_pivots[i])); + } + for (int i = 0; i < count; ++i) + for (int j = 0; j < i; ++j) + values(i) -= m_lu(i, j) * values(j); + for (int i = count - 1; i >= 0; --i) { + for (int j = i + 1; j < count; ++j) + values(i) -= m_lu(i, j) * values(j); + values(i) /= m_lu(i, i); + } + } + + const Operation *m_operation; + const Inverse *m_inverse; + SurfaceBorderSchurCoordinates m_coordinates; + MPI_Comm m_communicator; + SurfaceBorderSchurOptions m_options; + std::uint64_t m_generation{0}; + mfem::DenseMatrix m_schur, m_coarse, m_lu; + std::vector m_columns; + std::vector m_pivots; + double m_scale{0.0}; + mutable mfem::Vector m_rhs, m_fieldAction, m_image, m_solveWork; + mutable SurfaceBorderSchurStatistics m_statistics; + }; +} // namespace experiment::plateau_diagnostics diff --git a/tests/solver/plateau_surface_border_schur.cpp b/tests/solver/plateau_surface_border_schur.cpp new file mode 100644 index 0000000..868838b --- /dev/null +++ b/tests/solver/plateau_surface_border_schur.cpp @@ -0,0 +1,336 @@ +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +import mean_field; + +#include "../../experiments/plateau_surface_border_schur.hpp" + +namespace surface_border_schur_test { + namespace probe = experiment::plateau_diagnostics; + using namespace mean_field; + + mfem::DenseMatrix Matrix(int size, std::initializer_list entries) { + mfem::DenseMatrix result(size); + auto entry = entries.begin(); + for (int row = 0; row < size; ++row) + for (int column = 0; column < size; ++column) + result(row, column) = *entry++; + return result; + } + + mfem::Vector Vector(std::initializer_list values) { + mfem::Vector result(values.size()); + int i = 0; + for (double value : values) + result(i++) = value; + return result; + } + + double Difference(const mfem::Vector &left, const mfem::Vector &right) { + mfem::Vector difference(left); + difference -= right; + return difference.Norml2(); + } + + class Operation final : public mfem::Operator { + public: + explicit Operation(mfem::DenseMatrix matrix) + : mfem::Operator(matrix.Height()), matrix(std::move(matrix)) {} + void Mult(const mfem::Vector &input, mfem::Vector &output) const override { + ++applications; + matrix.Mult(input, output); + } + bool IsPrepared() const { return prepared; } + const Operation &GetProblem() const { return *this; } + std::uint64_t GetPreparationGeneration() const { return generation; } + mfem::DenseMatrix matrix; + bool prepared{true}; + std::uint64_t generation{1}; + mutable std::uint64_t applications{0}; + }; + + class Inverse : public mfem::Solver { + public: + static constexpr auto applicationContract = preconditioning::ApplicationContract::stationary_linear; + Inverse(const Operation &operation, mfem::DenseMatrix matrix) + : mfem::Solver(matrix.Height()), operation(operation), matrix(std::move(matrix)), + generation(operation.generation) {} + void SetOperator(const mfem::Operator &) override {} + void Mult(const mfem::Vector &input, mfem::Vector &output) const override { + ++applications; + matrix.Mult(input, output); + if (nonfinite) + output(0) = std::numeric_limits::infinity(); + } + bool IsCurrent() const { return generation == operation.generation; } + const Operation &GetNormalizedJacobian() const { return operation; } + const Operation &operation; + mfem::DenseMatrix matrix; + std::uint64_t generation; + bool nonfinite{false}; + mutable std::uint64_t applications{0}; + }; + + class FlexibleInverse : public Inverse { + public: + using Inverse::Inverse; + static constexpr auto applicationContract = preconditioning::ApplicationContract::flexible; + }; + static_assert(probe::SurfaceBorderSchurOperation); + static_assert(probe::SurfaceBorderSchurInverse); + static_assert(!probe::SurfaceBorderSchurInverse); + static_assert(!probe::SurfaceBorderSchurOperation); + + using Schur = probe::PreparedSurfaceBorderSchur; + static_assert(probe::SurfaceBorderSchurInverse); + static_assert(!std::is_copy_constructible_v); + + // Fields occupy different equation/unknown positions. A contiguous or + // same-offset assumption would silently factor the wrong block here. + probe::SurfaceBorderSchurCoordinates Coordinates() { return {{1, 3}, {0, 2}}; } + + mfem::DenseMatrix CoupledMatrix() { + return Matrix(4, {0.3, 4, -1, 1, + 0, 2, 2, 0.4, + 0.5, 2, 0.7, 3, + 3, -0.2, 1, 0.9}); + } + + mfem::DenseMatrix FieldInverse(bool exact) { + mfem::DenseMatrix result(4); + result = 0.0; + result(1, 0) = exact ? 0.3 : 0.2; + result(1, 2) = exact ? -0.1 : 0.0; + result(3, 0) = exact ? -0.2 : 0.0; + result(3, 2) = exact ? 0.4 : 0.2; + // These coarse contributions must be excluded, not accidentally used + // to complete the field response twice. + result(0, 0) = 11; + result(2, 2) = -17; + result(1, 1) = 23; + result(3, 3) = -31; + return result; + } + + namespace blocks = utils::blocks; + struct StructureDescriptor { + using CorrectionBlocks = blocks::type_list; + using ResidualBlocks = blocks::type_list; + }; + struct RootDescriptor { + using StructureComponent = StructureDescriptor; + using CorrectionBlocks = blocks::type_list; + using ResidualBlocks = blocks::type_list; + static constexpr std::size_t borderValueArity = 2, borderResidualArity = 2; + }; + struct CoordinateMap { + using Range = preconditioning::EquilibriumCoordinateRange; + std::array columns{{{5, 0, 2}, {1, 2, 3}, {0, 5, 1}, {4, 6, 1}, {7, 7, 1}}}; + std::array rows{{{6, 0, 1}, {2, 1, 3}, {0, 4, 2}, {5, 6, 1}, {7, 7, 1}}}; + const auto &GetCorrectionRanges() const { return columns; } + const auto &GetResidualRanges() const { return rows; } + }; + struct CoordinateBinding { + RootDescriptor root; + CoordinateMap coordinates; + const auto &GetBlock() const { return root; } + const auto &GetCoordinateMap() const { return coordinates; } + }; +} + +TEST_CASE("Surface Border Schur Coordinates Include All Native Modes And Generated Central Border", + "[solver][surface_border_schur][serial]") { + using namespace surface_border_schur_test; + const auto coordinates = probe::MakeSurfaceBorderSchurCoordinates(CoordinateBinding{}); + CHECK(coordinates.correctionIndices == std::vector{1, 2, 3, 4, 7}); + CHECK(coordinates.residualIndices == std::vector{2, 3, 4, 5, 7}); +} + +TEST_CASE("Surface Border Schur With Exact Field Inverse Solves Full Nonsymmetric System", + "[solver][surface_border_schur][serial]") { + using namespace surface_border_schur_test; + Operation operation(CoupledMatrix()); + Inverse inverse(operation, FieldInverse(true)); + Schur schur(operation, inverse, Coordinates(), MPI_COMM_WORLD); + REQUIRE(schur.IsCurrent()); + CHECK(schur.GetStatistics().setupOperatorApplications == 4); + CHECK(schur.GetStatistics().setupInverseApplications == 2); + const auto rhs = Vector({1, -2, 3, 4}); + mfem::Vector action(4), expected(4), defect(4); + mfem::DenseMatrixInverse direct(operation.matrix); + direct.Mult(rhs, expected); + schur.Mult(rhs, action); + CHECK(Difference(action, expected) < 1.e-13); + operation.Mult(action, defect); + CHECK(Difference(defect, rhs) < 1.e-13); + CHECK(schur.GetStatistics().operatorApplications == 1); + CHECK(schur.GetStatistics().inverseApplications == 1); + for (int column = 0; column < 2; ++column) { + const auto &completed = schur.GetCompletedColumns()[column]; + operation.Mult(completed, defect); + CHECK(std::abs(defect(0)) < 1.e-13); + CHECK(std::abs(defect(2)) < 1.e-13); + CHECK(completed(Coordinates().correctionIndices[column]) == 1.0); + } +} + +TEST_CASE("Surface Border Schur With Approximate Fields Enforces Every Coarse Equation", + "[solver][surface_border_schur][serial]") { + using namespace surface_border_schur_test; + Operation operation(CoupledMatrix()); + Inverse inverse(operation, FieldInverse(false)); + Schur schur(operation, inverse, Coordinates(), MPI_COMM_WORLD); + REQUIRE(schur.IsCurrent()); + const auto rhs = Vector({1, -2, 3, 4}); + mfem::Vector action(4), image(4); + schur.Mult(rhs, action); + operation.Mult(action, image); + for (int index : Coordinates().residualIndices) + CHECK(std::abs(image(index) - rhs(index)) < 1.e-13); + CHECK(Difference(image, rhs) > 1.e-3); + // Verify S = D - A_CF B A_FC independently in this explicit packing. + const auto &s = schur.GetSchurMatrix(); + CHECK(s(0, 0) == Catch::Approx(-0.16)); + CHECK(s(0, 1) == Catch::Approx(2.344)); + CHECK(s(1, 0) == Catch::Approx(2.922)); + CHECK(s(1, 1) == Catch::Approx(0.834)); +} + +TEST_CASE("Surface Border Schur Is Stationary Linear And Supports Aliasing", + "[solver][surface_border_schur][serial]") { + using namespace surface_border_schur_test; + Operation operation(CoupledMatrix()); + Inverse inverse(operation, FieldInverse(false)); + Schur schur(operation, inverse, Coordinates(), MPI_COMM_WORLD); + const auto first = Vector({1, 2, -1, 3}), second = Vector({-2, 4, 0.5, 1}); + mfem::Vector a(4), b(4), combined(first), expected(4), observed(4); + schur.Mult(first, a); + schur.Mult(second, b); + combined *= 2.3; + combined.Add(-0.7, second); + expected = a; + expected *= 2.3; + expected.Add(-0.7, b); + schur.Mult(combined, observed); + CHECK(Difference(observed, expected) < 1.e-13); + schur.Mult(combined, combined); + CHECK(Difference(combined, expected) < 1.e-13); + schur.Mult(first, observed); + CHECK(Difference(observed, a) == 0.0); +} + +TEST_CASE("Surface Border Schur Reports Singular And Ill Conditioned Coarse Blocks Without Regularization", + "[solver][surface_border_schur][serial]") { + using namespace surface_border_schur_test; + Operation operation(Matrix(3, {1, 0, 0, 0, 0, 0, 0, 0, 1})); + Inverse inverse(operation, Matrix(3, {1, 0, 0, 0, 0, 0, 0, 0, 0})); + Schur singular(operation, inverse, {{1, 2}, {1, 2}}, MPI_COMM_WORLD); + CHECK(singular.GetStatistics().status == probe::SurfaceBorderSchurStatus::singular); + CHECK_FALSE(singular.IsCurrent()); + auto rhs = Vector({1, 2, 3}); + mfem::Vector action(3); + CHECK_THROWS_AS(singular.Mult(rhs, action), std::logic_error); + operation.matrix(1, 1) = 1.e-16; + ++operation.generation; + inverse.generation = operation.generation; + singular.Refresh(); + CHECK(singular.GetStatistics().status == probe::SurfaceBorderSchurStatus::ill_conditioned); + CHECK(singular.GetStatistics().reciprocalCondition == Catch::Approx(1.e-16)); + CHECK_FALSE(singular.IsCurrent()); + CHECK_THROWS_AS(singular.SetOperator(operation), std::logic_error); + CHECK(singular.GetSchurMatrix()(0, 0) == 1.e-16); +} + +TEST_CASE("Surface Border Schur Handles Pivoting And Uniform Matrix Scale", + "[solver][surface_border_schur][serial]") { + using namespace surface_border_schur_test; + for (double scale : {1.e-100, 1.0, 1.e100}) { + auto matrix = Matrix(3, {2, 0, 0, 0, 0, 2, 0, 3, 1}); + matrix *= scale; + Operation operation(std::move(matrix)); + auto base = Matrix(3, {0.5, 0, 0, 0, 0, 0, 0, 0, 0}); + base *= 1.0 / scale; + Inverse inverse(operation, std::move(base)); + Schur schur(operation, inverse, {{1, 2}, {1, 2}}, MPI_COMM_WORLD); + REQUIRE(schur.IsCurrent()); + auto rhs = Vector({1, 2, 4}); + rhs *= scale; + mfem::Vector action(3); + schur.Mult(rhs, action); + CHECK(Difference(action, Vector({0.5, 1, 1})) < 1.e-13); + CHECK(schur.GetStatistics().reciprocalCondition == Catch::Approx(0.5)); + } +} + +TEST_CASE("Surface Border Schur Requires Unbordered Fields For A Zero Scalar Diagonal", + "[solver][surface_border_schur][serial]") { + using namespace surface_border_schur_test; + Operation operation(Matrix(2, {1, 1, 1, 0})); + Inverse alreadyBordered(operation, Matrix(2, {0, 1, 1, -1})); + Schur wrongRestriction(operation, alreadyBordered, {{1}, {1}}, MPI_COMM_WORLD); + CHECK(wrongRestriction.GetStatistics().status == probe::SurfaceBorderSchurStatus::singular); + Inverse unbordered(operation, Matrix(2, {1, 0, 0, 0})); + Schur schur(operation, unbordered, {{1}, {1}}, MPI_COMM_WORLD); + REQUIRE(schur.IsCurrent()); + auto rhs = Vector({3, 2}); + mfem::Vector action(2); + schur.Mult(rhs, action); + CHECK(Difference(action, Vector({2, 1})) < 1.e-13); +} + +TEST_CASE("Surface Border Schur Rejects Invalid Partitions And Stale Or Nonfinite Dependencies", + "[solver][surface_border_schur][serial]") { + using namespace surface_border_schur_test; + Operation operation(CoupledMatrix()); + Inverse inverse(operation, FieldInverse(true)); + const auto reject = [&](probe::SurfaceBorderSchurCoordinates coordinates) { + CHECK_THROWS_AS(Schur(operation, inverse, std::move(coordinates), MPI_COMM_WORLD), std::invalid_argument); + }; + reject({{}, {}}); + reject({{1}, {0, 2}}); + reject({{1, 1}, {0, 2}}); + reject({{1, 3}, {0, 0}}); + reject({{1, 4}, {0, 2}}); + reject({{0, 1, 2, 3}, {0, 1, 2, 3}}); + CHECK_THROWS_AS(Schur(operation, inverse, Coordinates(), MPI_COMM_NULL), std::invalid_argument); + Schur schur(operation, inverse, Coordinates(), MPI_COMM_WORLD); + auto rhs = Vector({1, 2, 3, 4}); + mfem::Vector action(4); + Operation other(CoupledMatrix()); + CHECK_THROWS_AS(schur.SetOperator(other), std::invalid_argument); + CHECK_THROWS_AS(Schur(other, inverse, Coordinates(), MPI_COMM_WORLD), std::invalid_argument); + ++operation.generation; + CHECK_FALSE(schur.IsCurrent()); + CHECK_THROWS_AS(schur.Mult(rhs, action), std::logic_error); + inverse.generation = operation.generation; + CHECK_FALSE(schur.IsCurrent()); + schur.Refresh(); + REQUIRE(schur.IsCurrent()); + inverse.iterative_mode = true; + CHECK_FALSE(schur.IsCurrent()); + CHECK_THROWS_AS(schur.Mult(rhs, action), std::logic_error); + inverse.iterative_mode = false; + inverse.nonfinite = true; + CHECK_THROWS_AS(schur.Mult(rhs, action), std::runtime_error); + inverse.nonfinite = false; + rhs(0) = std::numeric_limits::quiet_NaN(); + CHECK_THROWS_AS(schur.Mult(rhs, action), std::runtime_error); +}