From 522a2cfc021c8ba0b803b4d732c8e56461f10f70 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Sat, 19 Sep 2026 13:45:57 -0400 Subject: [PATCH] Test normalized unbordered adapter and report its storage scope --- .../plateau_preconditioner_benchmark.hpp | 3 + tests/solver/plateau_surface_border_schur.cpp | 93 ++++++++++++++++++- 2 files changed, 95 insertions(+), 1 deletion(-) diff --git a/experiments/plateau_preconditioner_benchmark.hpp b/experiments/plateau_preconditioner_benchmark.hpp index 4d3ff88..83bfd0c 100644 --- a/experiments/plateau_preconditioner_benchmark.hpp +++ b/experiments/plateau_preconditioner_benchmark.hpp @@ -354,6 +354,8 @@ void BenchmarkPreconditioners(State &state, const std::filesystem::path &output, metric("setup_jacobian_seconds", stats.setupOperatorSeconds); metric("setup_inverse_seconds", stats.setupInverseSeconds); metric("persistent_vector_bytes", stats.persistentVectorBytes); + metric("unbordered_adapter_vector_bytes", sizeof(mfem::real_t) * + 4 * static_cast(state.normalizedOperator->Height())); metric("applications", stats.applications); metric("internal_jacobian_applications", stats.operatorApplications); metric("internal_inverse_applications", stats.inverseApplications); @@ -643,6 +645,7 @@ void BenchmarkPreconditioners(State &state, const std::filesystem::path &output, << "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" + << "The unbordered adapter is frozen at construction; rebuild it and the Schur wrapper for a new Jacobian.\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/tests/solver/plateau_surface_border_schur.cpp b/tests/solver/plateau_surface_border_schur.cpp index 868838b..6afd6d4 100644 --- a/tests/solver/plateau_surface_border_schur.cpp +++ b/tests/solver/plateau_surface_border_schur.cpp @@ -45,7 +45,11 @@ namespace surface_border_schur_test { class Operation final : public mfem::Operator { public: explicit Operation(mfem::DenseMatrix matrix) - : mfem::Operator(matrix.Height()), matrix(std::move(matrix)) {} + : mfem::Operator(matrix.Height()), matrix(std::move(matrix)), + residualScales(Height()), stateScales(Width()) { + residualScales = 1.0; + stateScales = 1.0; + } void Mult(const mfem::Vector &input, mfem::Vector &output) const override { ++applications; matrix.Mult(input, output); @@ -53,7 +57,16 @@ namespace surface_border_schur_test { bool IsPrepared() const { return prepared; } const Operation &GetProblem() const { return *this; } std::uint64_t GetPreparationGeneration() const { return generation; } + void DenormalizeResidual(const mfem::Vector &input, mfem::Vector &output) const { + for (int i = 0; i < input.Size(); ++i) + output(i) = input(i) * residualScales(i); + } + void NormalizeState(const mfem::Vector &input, mfem::Vector &output) const { + for (int i = 0; i < input.Size(); ++i) + output(i) = input(i) / stateScales(i); + } mfem::DenseMatrix matrix; + mfem::Vector residualScales, stateScales; bool prepared{true}; std::uint64_t generation{1}; mutable std::uint64_t applications{0}; @@ -124,6 +137,7 @@ namespace surface_border_schur_test { namespace blocks = utils::blocks; struct StructureDescriptor { + using BackendType = preconditioning::backend::Diagonal; using CorrectionBlocks = blocks::type_list; using ResidualBlocks = blocks::type_list 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; } + void PackResidual(const mfem::Vector &input, mfem::Vector &output) const { + for (const auto &range : rows) + for (int i = 0; i < range.size; ++i) + output(range.preconditionerOffset + i) = input(range.equilibriumOffset + i); + } + void UnpackCorrection(const mfem::Vector &input, mfem::Vector &output) const { + for (const auto &range : columns) + for (int i = 0; i < range.size; ++i) + output(range.equilibriumOffset + i) = input(range.preconditionerOffset + i); + } }; struct CoordinateBinding { RootDescriptor root; @@ -154,6 +178,73 @@ namespace surface_border_schur_test { const auto &GetBlock() const { return root; } const auto &GetCoordinateMap() const { return coordinates; } }; + + class PackedStructure final : public mfem::Solver { + public: + PackedStructure() : mfem::Solver(6), matrix(6) { + matrix = 0.0; + for (int i = 0; i < 6; ++i) + matrix(i, i) = i + 2.0; + matrix(0, 5) = 0.4; + matrix(4, 1) = -0.6; + } + void SetOperator(const mfem::Operator &) override {} + void Mult(const mfem::Vector &input, mfem::Vector &output) const override { + ++applications; + matrix.Mult(input, output); + } + mfem::DenseMatrix matrix; + mutable unsigned applications{0}; + }; + struct GroupedStructure { + PackedStructure structure; + const auto &GetStructurePreconditioner() const { return structure; } + }; + class PhysicalStructureBinding final : public mfem::Solver, public CoordinateBinding { + public: + static constexpr auto applicationContract = preconditioning::ApplicationContract::stationary_linear; + explicit PhysicalStructureBinding(const Operation &operation) + : mfem::Solver(operation.Width()), operation(operation), generation(operation.generation) {} + void SetOperator(const mfem::Operator &) override {} + void Mult(const mfem::Vector &, mfem::Vector &) const override { + throw std::logic_error("The adapter must not apply the scalar-bordered inverse."); + } + const auto &GetProblem() const { return operation; } + const auto &GetGroupedPreconditioner() const { return grouped; } + bool IsCurrent() const { return generation == operation.generation; } + const Operation &operation; + std::uint64_t generation; + GroupedStructure grouped; + }; +} + +TEST_CASE("Surface Border Schur Structure Adapter Preserves Nonunit Normalization And Typed Packing", + "[solver][surface_border_schur][serial]") { + using namespace surface_border_schur_test; + mfem::DenseMatrix matrix(8); + matrix = 0.0; + Operation operation(std::move(matrix)); + operation.residualScales = Vector({2, 3, 5, 7, 11, 13, 17, 19}); + operation.stateScales = Vector({23, 29, 31, 37, 41, 43, 47, 53}); + PhysicalStructureBinding physical(operation); + probe::NormalizedUnborderedStructureInverse adapter(operation, physical); + STATIC_CHECK(probe::SurfaceBorderSchurInverse); + auto rhs = Vector({1, 2, 3, 4, 5, 6, 7, 8}); + const auto expected = Vector({42.0 / 23, 112.0 / 29, 275.0 / 31, 3.0 / 37, + 0, 240.4 / 43, 45.0 / 47, 0}); + mfem::Vector observed(8); + adapter.Mult(rhs, observed); + CHECK(Difference(observed, expected) < 1.e-13); + CHECK(physical.grouped.structure.applications == 1); + adapter.Mult(rhs, rhs); + CHECK(Difference(rhs, expected) < 1.e-13); + physical.grouped.structure.iterative_mode = true; + CHECK_FALSE(adapter.IsCurrent()); + physical.grouped.structure.iterative_mode = false; + ++operation.generation; + physical.generation = operation.generation; + CHECK_FALSE(adapter.IsCurrent()); + CHECK_THROWS_AS(adapter.Mult(rhs, observed), std::logic_error); } TEST_CASE("Surface Border Schur Coordinates Include All Native Modes And Generated Central Border",