diff --git a/examples/python-api-example/example_conservative_transfer.py b/examples/python-api-example/example_conservative_transfer.py index d6e6fd54..89f6e212 100644 --- a/examples/python-api-example/example_conservative_transfer.py +++ b/examples/python-api-example/example_conservative_transfer.py @@ -216,7 +216,20 @@ def main(plot_path=None): report("mesh intersection (exact quadrature)", mi_target, target_mesh, target_coords, reference_values, source_integral) - # 3b. Monte Carlo (control-variate) conservative projection. + # 3b. Mesh-intersection projection with a row-sum lumped mass matrix. + # The linear solve degenerates to an exact diagonal scaling; the + # integral is still conserved exactly, but the max nodal error + # grows relative to the consistent mass matrix. + lumped_target = target_space.create_field() + lumped_projection = pcms.OmegaHConservativeProjection( + source_space, target_space, + mass_matrix_type=pcms.MassMatrixType.Lumped, + ) + lumped_projection.apply(source_field, lumped_target) + report("mesh intersection (lumped mass)", lumped_target, target_mesh, + target_coords, reference_values, source_integral) + + # 3c. Monte Carlo (control-variate) conservative projection. mc_target = target_space.create_field() monte_carlo = pcms.OmegaHControlVariateProjection( source_space, diff --git a/src/pcms/pythonapi/bind_transfer_field.cpp b/src/pcms/pythonapi/bind_transfer_field.cpp index cc599bb0..1947a8bb 100644 --- a/src/pcms/pythonapi/bind_transfer_field.cpp +++ b/src/pcms/pythonapi/bind_transfer_field.cpp @@ -11,6 +11,7 @@ #include "numpy_array_transform.h" #include "pcms/utility/types.h" #if defined(PCMS_ENABLE_PETSC) && defined(PCMS_ENABLE_MESHFIELDS) +#include "pcms/transfer/mass_matrix_type.hpp" #include "pcms/transfer/omega_h_conservative_projection.hpp" #include "pcms/transfer/omega_h_control_variate_projection.hpp" #include "pcms/transfer/omega_h_mc_rhs_integrator.hpp" @@ -79,16 +80,28 @@ void bind_transfer_field_module(py::module& m) "Copy source field data to target field (same layout required)."); #if defined(PCMS_ENABLE_PETSC) && defined(PCMS_ENABLE_MESHFIELDS) + // Mass-matrix formulation for the Galerkin projections. Bound before the + // projection classes so it can be used as a default argument value. + py::enum_(m, "MassMatrixType") + .value("Consistent", MassMatrixType::Consistent, + "Full Galerkin mass matrix") + .value("Lumped", MassMatrixType::Lumped, + "Row-sum lumped diagonal mass matrix (exact diagonal solve; " + "conserves the integral but adds lumping error)") + .export_values(); + // Conservative L2 (Galerkin) projection between order-1 Lagrange spaces on // Omega_h 2D simplex meshes. The RHS is integrated exactly over the // intersection of the source and target meshes. py::class_(m, "OmegaHConservativeProjection") .def(py::init([](const FunctionSpace& source_space, - const FunctionSpace& target_space) { - return std::make_unique(source_space, - target_space); + const FunctionSpace& target_space, + MassMatrixType mass_matrix_type) { + return std::make_unique( + source_space, target_space, mass_matrix_type); }), py::arg("source_space"), py::arg("target_space"), + py::arg("mass_matrix_type") = MassMatrixType::Consistent, "Construct a mesh-intersection conservative projection. Mesh " "intersection, quadrature setup, and mass-matrix factorization happen " "here and are cached; call apply() repeatedly.") @@ -114,14 +127,17 @@ void bind_transfer_field_module(py::module& m) "OmegaHControlVariateProjection") .def(py::init([](const FunctionSpace& source_space, const FunctionSpace& target_space, int samples_per_element, - MonteCarloSampling sampling, std::uint64_t seed) { + MonteCarloSampling sampling, std::uint64_t seed, + MassMatrixType mass_matrix_type) { return std::make_unique( - source_space, target_space, samples_per_element, sampling, seed); + source_space, target_space, samples_per_element, sampling, seed, + mass_matrix_type); }), py::arg("source_space"), py::arg("target_space"), py::arg("samples_per_element"), py::arg("sampling") = MonteCarloSampling::UniformRandom, py::arg("seed") = std::uint64_t(8675309), + py::arg("mass_matrix_type") = MassMatrixType::Consistent, "Construct a Monte Carlo (control-variate) conservative projection. " "Sample-point generation, the control-variate interpolator, and the " "mass-matrix factorization are cached; call apply() repeatedly.") diff --git a/src/pcms/transfer/CMakeLists.txt b/src/pcms/transfer/CMakeLists.txt index c3983e0b..ad62383b 100644 --- a/src/pcms/transfer/CMakeLists.txt +++ b/src/pcms/transfer/CMakeLists.txt @@ -7,6 +7,8 @@ set(PCMS_FIELD_TRANSFER_HEADERS interpolator.h copy.h transfer_operator.hpp + mass_matrix_type.hpp + monte_carlo_sampling.hpp ) diff --git a/src/pcms/transfer/bilinear_form_integrator.hpp b/src/pcms/transfer/bilinear_form_integrator.hpp index 060667bd..df8ace69 100644 --- a/src/pcms/transfer/bilinear_form_integrator.hpp +++ b/src/pcms/transfer/bilinear_form_integrator.hpp @@ -14,6 +14,10 @@ class BilinearFormIntegrator // counting keeps the matrix alive until the KSP is destroyed). virtual Mat GetMatrix() const noexcept = 0; + // True when the assembled matrix is diagonal; lets solvers pick an exact + // one-application diagonal solve (e.g. preonly + Jacobi). + virtual bool IsDiagonal() const noexcept { return false; } + virtual ~BilinearFormIntegrator() noexcept = default; }; diff --git a/src/pcms/transfer/conservative_projection_solver.cpp b/src/pcms/transfer/conservative_projection_solver.cpp index e6955580..fc0ad8ec 100644 --- a/src/pcms/transfer/conservative_projection_solver.cpp +++ b/src/pcms/transfer/conservative_projection_solver.cpp @@ -54,6 +54,20 @@ GalerkinProjectionSolver::GalerkinProjectionSolver( CHKERRABORT(PETSC_COMM_WORLD, ierr); ierr = KSPSetOperators(ksp_, A, A); CHKERRABORT(PETSC_COMM_WORLD, ierr); + if (mass_integrator.IsDiagonal()) { + // Jacobi preconditioning applies inv(diag(A)), which for a diagonal + // matrix is exactly inv(A). Pairing it with KSPPREONLY therefore gives + // the exact solution in a single application — no Krylov iterations, no + // factorization, no convergence tolerance in play. Command line options + // can still override via KSPSetFromOptions below. + ierr = KSPSetType(ksp_, KSPPREONLY); + CHKERRABORT(PETSC_COMM_WORLD, ierr); + PC pc = nullptr; + ierr = KSPGetPC(ksp_, &pc); + CHKERRABORT(PETSC_COMM_WORLD, ierr); + ierr = PCSetType(pc, PCJACOBI); + CHKERRABORT(PETSC_COMM_WORLD, ierr); + } ierr = KSPSetFromOptions(ksp_); CHKERRABORT(PETSC_COMM_WORLD, ierr); ierr = KSPSetUp(ksp_); diff --git a/src/pcms/transfer/mass_matrix_type.hpp b/src/pcms/transfer/mass_matrix_type.hpp new file mode 100644 index 00000000..15898813 --- /dev/null +++ b/src/pcms/transfer/mass_matrix_type.hpp @@ -0,0 +1,16 @@ +#ifndef PCMS_TRANSFER_MASS_MATRIX_TYPE_HPP +#define PCMS_TRANSFER_MASS_MATRIX_TYPE_HPP + +namespace pcms +{ + +// Mass-matrix formulation used by Galerkin projections. +enum class MassMatrixType +{ + Consistent, // full Galerkin mass matrix + Lumped, // row-sum lumped diagonal mass matrix +}; + +} // namespace pcms + +#endif // PCMS_TRANSFER_MASS_MATRIX_TYPE_HPP diff --git a/src/pcms/transfer/omega_h_conservative_projection.cpp b/src/pcms/transfer/omega_h_conservative_projection.cpp index 8230da22..31b1b9e7 100644 --- a/src/pcms/transfer/omega_h_conservative_projection.cpp +++ b/src/pcms/transfer/omega_h_conservative_projection.cpp @@ -41,7 +41,8 @@ void CheckApplyCompatible(const Field& source, const Field& target, } // namespace OmegaHConservativeProjection::OmegaHConservativeProjection( - const FunctionSpace& source_space, const FunctionSpace& target_space) + const FunctionSpace& source_space, const FunctionSpace& target_space, + MassMatrixType mass_matrix_type) : source_layout_(std::dynamic_pointer_cast( source_space.GetLayout())), target_layout_(std::dynamic_pointer_cast( @@ -57,8 +58,8 @@ OmegaHConservativeProjection::OmegaHConservativeProjection( // Mass integrator is only needed to build the solver; PETSc reference-counts // the matrix so it remains alive inside the KSP after this scope ends. - OmegaHMassIntegrator mass_integrator(target_layout_, - target_space.GetCoordinateSystem()); + OmegaHMassIntegrator mass_integrator( + target_layout_, target_space.GetCoordinateSystem(), mass_matrix_type); solver_ = std::make_unique(mass_integrator, *rhs_integrator_); target_values_ = Kokkos::View( diff --git a/src/pcms/transfer/omega_h_conservative_projection.hpp b/src/pcms/transfer/omega_h_conservative_projection.hpp index 18b7a2cb..2b2abf3d 100644 --- a/src/pcms/transfer/omega_h_conservative_projection.hpp +++ b/src/pcms/transfer/omega_h_conservative_projection.hpp @@ -4,6 +4,7 @@ #include "pcms/field/function_space.h" #include "pcms/field/layout/omega_h_lagrange.h" #include "pcms/field/point_evaluator.h" +#include "pcms/transfer/mass_matrix_type.hpp" #include "pcms/transfer/transfer_operator.hpp" #include #include @@ -21,7 +22,10 @@ class OmegaHIntersectionRHSIntegrator; // (OmegaHIntersectionRHSIntegrator). // - Localizes integration points in the source mesh (PointEvaluator). // - Assembles the target mass matrix and factors it via KSP -// (GalerkinProjectionSolver). +// (GalerkinProjectionSolver). With MassMatrixType::Lumped the mass +// matrix is row-sum lumped and the KSP degenerates to an exact diagonal +// scaling; conservation is unaffected, but linear fields are no longer +// reproduced exactly on order-1 targets. // // Apply() (per-call cost): // - Evaluates the source field at the fixed integration points. @@ -31,8 +35,9 @@ class OmegaHIntersectionRHSIntegrator; class OmegaHConservativeProjection : public TransferOperator { public: - OmegaHConservativeProjection(const FunctionSpace& source_space, - const FunctionSpace& target_space); + OmegaHConservativeProjection( + const FunctionSpace& source_space, const FunctionSpace& target_space, + MassMatrixType mass_matrix_type = MassMatrixType::Consistent); // Defined in the .cpp where GalerkinProjectionSolver is a complete type. ~OmegaHConservativeProjection() override; diff --git a/src/pcms/transfer/omega_h_control_variate_projection.cpp b/src/pcms/transfer/omega_h_control_variate_projection.cpp index 256f92a6..f253f74e 100644 --- a/src/pcms/transfer/omega_h_control_variate_projection.cpp +++ b/src/pcms/transfer/omega_h_control_variate_projection.cpp @@ -11,7 +11,8 @@ namespace pcms OmegaHControlVariateProjection::OmegaHControlVariateProjection( const FunctionSpace& source_space, const FunctionSpace& target_space, - int samples_per_element, MonteCarloSampling sampling, uint64_t seed) + int samples_per_element, MonteCarloSampling sampling, uint64_t seed, + MassMatrixType mass_matrix_type) : target_layout_(std::dynamic_pointer_cast( target_space.GetLayout())), rhs_integrator_(std::make_unique( @@ -28,8 +29,8 @@ OmegaHControlVariateProjection::OmegaHControlVariateProjection( // Mass integrator is only needed to build the solver; PETSc reference-counts // the matrix so it remains alive inside the KSP after this scope ends. - OmegaHMassIntegrator mass_integrator(target_layout_, - target_space.GetCoordinateSystem()); + OmegaHMassIntegrator mass_integrator( + target_layout_, target_space.GetCoordinateSystem(), mass_matrix_type); solver_ = std::make_unique(mass_integrator, *rhs_integrator_); } diff --git a/src/pcms/transfer/omega_h_control_variate_projection.hpp b/src/pcms/transfer/omega_h_control_variate_projection.hpp index 4196ea5d..0ec79dea 100644 --- a/src/pcms/transfer/omega_h_control_variate_projection.hpp +++ b/src/pcms/transfer/omega_h_control_variate_projection.hpp @@ -6,6 +6,7 @@ #include "pcms/field/layout/omega_h_lagrange.h" #include "pcms/field/point_evaluator.h" #include "pcms/transfer/interpolator.h" +#include "pcms/transfer/mass_matrix_type.hpp" #include "pcms/transfer/monte_carlo_sampling.hpp" #include "pcms/transfer/transfer_operator.hpp" #include @@ -37,11 +38,11 @@ class OmegaHMonteCarloRHSIntegrator; class OmegaHControlVariateProjection : public TransferOperator { public: - OmegaHControlVariateProjection(const FunctionSpace& source_space, - const FunctionSpace& target_space, - int samples_per_element, - MonteCarloSampling sampling, - uint64_t seed = 8675309); + OmegaHControlVariateProjection( + const FunctionSpace& source_space, const FunctionSpace& target_space, + int samples_per_element, MonteCarloSampling sampling, + uint64_t seed = 8675309, + MassMatrixType mass_matrix_type = MassMatrixType::Consistent); // Defined in the .cpp where GalerkinProjectionSolver is a complete type. ~OmegaHControlVariateProjection() override; diff --git a/src/pcms/transfer/omega_h_mass_integrator.cpp b/src/pcms/transfer/omega_h_mass_integrator.cpp index 29cd90e2..1a6fc5f6 100644 --- a/src/pcms/transfer/omega_h_mass_integrator.cpp +++ b/src/pcms/transfer/omega_h_mass_integrator.cpp @@ -44,11 +44,14 @@ void FillP0MassCoo( } // Fills the 3x3-block COO sparsity pattern of a P1 (linear) mass matrix: each -// element contributes a dense block coupling its three vertex DOFs. +// element contributes a dense block coupling its three vertex DOFs. When +// lumped, every block entry is mapped onto the row's diagonal instead; PETSc's +// COO assembly sums repeated indices, so the unmodified element values +// accumulate into the row-sum lumped diagonal. void FillP1MassCooPattern( int nelems, const Omega_h::LOs& faces2nodes, const Kokkos::View& global_to_local, - const Kokkos::View& coo_rows, + bool lumped, const Kokkos::View& coo_rows, const Kokkos::View& coo_cols) { Kokkos::parallel_for( @@ -57,8 +60,10 @@ void FillP1MassCooPattern( for (int i = 0; i < 3; ++i) { for (int j = 0; j < 3; ++j) { const int idx = e * 9 + i * 3 + j; - coo_rows(idx) = static_cast(global_to_local(verts[i])); - coo_cols(idx) = static_cast(global_to_local(verts[j])); + const auto row = static_cast(global_to_local(verts[i])); + coo_rows(idx) = row; + coo_cols(idx) = + lumped ? row : static_cast(global_to_local(verts[j])); } } }); @@ -66,16 +71,17 @@ void FillP1MassCooPattern( } // namespace -OmegaHMassIntegrator::OmegaHMassIntegrator(const FunctionSpace& target_space) +OmegaHMassIntegrator::OmegaHMassIntegrator(const FunctionSpace& target_space, + MassMatrixType mass_type) : OmegaHMassIntegrator(std::dynamic_pointer_cast( target_space.GetLayout()), - target_space.GetCoordinateSystem()) + target_space.GetCoordinateSystem(), mass_type) { } OmegaHMassIntegrator::OmegaHMassIntegrator( std::shared_ptr target_layout, - CoordinateSystem coordinate_system) + CoordinateSystem coordinate_system, MassMatrixType mass_type) { detail::CheckOmegaHScalarLagrangeLayout(coordinate_system, target_layout, "OmegaHMassIntegrator", "target"); @@ -85,6 +91,8 @@ OmegaHMassIntegrator::OmegaHMassIntegrator( const PetscInt num_dofs = static_cast(target_layout->GetNumOwnedDofHolder()); const int nelems = mesh.nelems(); + diagonal_ = + target_layout->GetOrder() == 0 || mass_type == MassMatrixType::Lumped; if (target_layout->GetOrder() == 0) { // P0 target: piecewise-constant basis functions have disjoint support, so @@ -132,8 +140,8 @@ OmegaHMassIntegrator::OmegaHMassIntegrator( const PetscInt nnz = static_cast(nelems) * 9; Kokkos::View coo_rows("mass_coo_rows", nnz); Kokkos::View coo_cols("mass_coo_cols", nnz); - FillP1MassCooPattern(nelems, faces2nodes, global_to_local, coo_rows, - coo_cols); + FillP1MassCooPattern(nelems, faces2nodes, global_to_local, + mass_type == MassMatrixType::Lumped, coo_rows, coo_cols); // Create sparse matrix, preallocate with COO pattern, then bulk-set values. // elm_mass_dev is in the same element-major order as coo_rows/coo_cols, so @@ -164,10 +172,15 @@ Mat OmegaHMassIntegrator::GetMatrix() const noexcept return mat_; } +bool OmegaHMassIntegrator::IsDiagonal() const noexcept +{ + return diagonal_; +} + std::unique_ptr BuildOmegaHMassIntegrator( - const FunctionSpace& target_space) + const FunctionSpace& target_space, MassMatrixType mass_type) { - return std::make_unique(target_space); + return std::make_unique(target_space, mass_type); } } // namespace pcms diff --git a/src/pcms/transfer/omega_h_mass_integrator.hpp b/src/pcms/transfer/omega_h_mass_integrator.hpp index 54cc7a66..b130d088 100644 --- a/src/pcms/transfer/omega_h_mass_integrator.hpp +++ b/src/pcms/transfer/omega_h_mass_integrator.hpp @@ -4,6 +4,7 @@ #include "pcms/field/function_space.h" #include "pcms/field/layout/omega_h_lagrange.h" #include "pcms/transfer/bilinear_form_integrator.hpp" +#include "pcms/transfer/mass_matrix_type.hpp" #include #include @@ -12,7 +13,9 @@ namespace pcms // Mass-matrix bilinear form integrator for Lagrange spaces on Omega_h 2D // simplex meshes. Supports order-0 (diagonal, M_ee = area(e)) and order-1 -// (consistent 3x3 element blocks) target spaces. +// (consistent 3x3 element blocks) target spaces. With +// MassMatrixType::Lumped the order-1 matrix is row-sum lumped into a +// diagonal matrix (order-0 is already diagonal, so lumping is a no-op). // // The matrix is assembled once at construction into a PETSc sparse AIJ matrix. // Sparsity is derived directly from the element node GIDs — no separate COO @@ -22,23 +25,29 @@ namespace pcms class OmegaHMassIntegrator : public BilinearFormIntegrator { public: - explicit OmegaHMassIntegrator(const FunctionSpace& target_space); + explicit OmegaHMassIntegrator( + const FunctionSpace& target_space, + MassMatrixType mass_type = MassMatrixType::Consistent); OmegaHMassIntegrator( std::shared_ptr target_layout, - CoordinateSystem coordinate_system); + CoordinateSystem coordinate_system, + MassMatrixType mass_type = MassMatrixType::Consistent); ~OmegaHMassIntegrator(); Mat GetMatrix() const noexcept override; + bool IsDiagonal() const noexcept override; private: Mat mat_ = nullptr; + bool diagonal_ = false; }; // Builds an OmegaHMassIntegrator for the given target space. // target_space must use OmegaHLagrangeLayout, scalar, Cartesian, order-0 or // order-1. std::unique_ptr BuildOmegaHMassIntegrator( - const FunctionSpace& target_space); + const FunctionSpace& target_space, + MassMatrixType mass_type = MassMatrixType::Consistent); } // namespace pcms diff --git a/test/test_mesh_intersection_field_transfer.cpp b/test/test_mesh_intersection_field_transfer.cpp index fbb34ddb..6b4cd8ae 100644 --- a/test/test_mesh_intersection_field_transfer.cpp +++ b/test/test_mesh_intersection_field_transfer.cpp @@ -1,5 +1,6 @@ #include #include +#include #include #include #include @@ -9,7 +10,9 @@ #include #include "field_test_utils.h" +#include #include +#include #include namespace @@ -192,9 +195,11 @@ TEST_CASE("OmegaHConservativeProjection reproduces constant and linear fields", // For a field that does not live in the target space the projection is not an // interpolation, but conservation of the integral is the defining property of -// the conservative (Galerkin) projection and must still hold exactly. This -// also exercises a second Apply with a different source field to confirm the -// cached integrators/evaluator/factorization are reused without stale state. +// the conservative (Galerkin) projection and must still hold exactly — for +// the lumped mass matrix too, since its row sums (and hence the discrete +// integral identity) match the consistent matrix. This also exercises a +// second Apply with a different source field to confirm the cached +// integrators/evaluator/factorization are reused without stale state. TEST_CASE("OmegaHConservativeProjection conserves the integral", "[transfer][mesh_intersection]") { @@ -211,7 +216,12 @@ TEST_CASE("OmegaHConservativeProjection conserves the integral", auto source = source_space->CreateFunction(); auto target = target_space->CreateFunction(); - pcms::OmegaHConservativeProjection projection(*source_space, *target_space); + const auto mass_type = + GENERATE(pcms::MassMatrixType::Consistent, pcms::MassMatrixType::Lumped); + CAPTURE(static_cast(mass_type)); + + pcms::OmegaHConservativeProjection projection(*source_space, *target_space, + mass_type); pcms::test::SetField( source, OMEGA_H_LAMBDA(pcms::Real x, pcms::Real y) { @@ -231,6 +241,79 @@ TEST_CASE("OmegaHConservativeProjection conserves the integral", Catch::Approx(IntegrateP1Field(source_mesh, source)).margin(1e-12)); } +// Row-sum lumping keeps constants exact (the lumped diagonal is exactly the +// row sum, so M_L c = M_C c for constant c) and keeps the integral conserved, +// but it is no longer a true L2 projection: linear fields are smeared toward +// patch averages instead of being reproduced nodally. The error bound below +// guards against the mass_matrix_type option being silently dropped anywhere +// between the operator ctor and the solver. +TEST_CASE("OmegaHConservativeProjection with lumped mass reproduces constants " + "but not linears", + "[transfer][mesh_intersection]") +{ + Omega_h::Library lib; + + Omega_h::Mesh source_mesh = + BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3})); + Omega_h::Mesh target_mesh = + BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3})); + + auto source_space = MakeP1Space(source_mesh); + auto target_space = MakeP1Space(target_mesh); + + auto source = source_space->CreateFunction(); + auto target = target_space->CreateFunction(); + + pcms::OmegaHConservativeProjection projection(*source_space, *target_space, + pcms::MassMatrixType::Lumped); + + SECTION("constant field is preserved and conserved") + { + const double c = 2.0; + pcms::test::SetField( + source, OMEGA_H_LAMBDA(pcms::Real, pcms::Real) { return c; }); + + projection.Apply(source, target); + + const auto target_values = + pcms::FlattenToRank1View(target.GetDOFHolderDataHost()); + REQUIRE(static_cast(target_values.size()) == + target_mesh.nverts()); + for (Omega_h::LO i = 0; i < target_mesh.nverts(); ++i) { + REQUIRE(target_values[i] == Catch::Approx(c).margin(1e-10)); + } + REQUIRE(IntegrateP1Field(target_mesh, target) == + Catch::Approx(IntegrateP1Field(source_mesh, source)).margin(1e-10)); + } + + SECTION("linear field conserves the integral but is not nodally exact") + { + pcms::test::SetField( + source, OMEGA_H_LAMBDA(pcms::Real x, pcms::Real y) { return x + y; }); + + projection.Apply(source, target); + + REQUIRE(IntegrateP1Field(target_mesh, target) == + Catch::Approx(IntegrateP1Field(source_mesh, source)).margin(1e-9)); + + // On this mesh the lumped solution at the corner vertices is the + // phi-weighted patch average of x + y, which is off by 0.5 — far above + // solver tolerance. A consistent mass matrix reproduces the field to + // 1e-9 (see the exact-reproduction test above), so a large error here + // proves the lumped path was actually taken. + const auto target_values = + pcms::FlattenToRank1View(target.GetDOFHolderDataHost()); + const auto tgt_coords_h = + Omega_h::HostRead(target_mesh.coords()); + double max_err = 0.0; + for (Omega_h::LO i = 0; i < target_mesh.nverts(); ++i) { + const double expected = tgt_coords_h[2 * i + 0] + tgt_coords_h[2 * i + 1]; + max_err = std::max(max_err, std::abs(target_values[i] - expected)); + } + CHECK(max_err > 1e-3); + } +} + TEST_CASE("OmegaHConservativeProjection writes reordered target GIDs in local " "DOF order", "[transfer][mesh_intersection]") diff --git a/test/test_omega_h_mass_integrator.cpp b/test/test_omega_h_mass_integrator.cpp index c9013f9d..a49de590 100644 --- a/test/test_omega_h_mass_integrator.cpp +++ b/test/test_omega_h_mass_integrator.cpp @@ -4,6 +4,7 @@ #include #include #include +#include #include #include #include @@ -166,6 +167,108 @@ TEST_CASE("OmegaHMassIntegrator: row sums match lumped mass", // mat is owned by integrator; no MatDestroy here. } +TEST_CASE("OmegaHMassIntegrator: lumped diagonal equals consistent row sums", + "[mass_integrator]") +{ + // Row-sum lumping: M_L(i,i) = sum_j M_C(i,j) = integral(N_i dx), all + // off-diagonal entries zero. The trace equals the domain area. + Omega_h::Library lib; + auto mesh = BuildUnitSquare(lib); + + auto space = pcms::LagrangeFunctionSpace::FromMesh( + mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", + pcms::LagrangeFunctionSpace::Backend::OmegaH); + + const auto layout = + std::dynamic_pointer_cast( + space->GetLayout()); + const auto& gids = layout->GetGidsHost(); + const int nverts = mesh.nverts(); + + auto consistent = pcms::BuildOmegaHMassIntegrator(*space); + auto lumped = + pcms::BuildOmegaHMassIntegrator(*space, pcms::MassMatrixType::Lumped); + CHECK_FALSE(consistent->IsDiagonal()); + CHECK(lumped->IsDiagonal()); + Mat consistent_mat = consistent->GetMatrix(); + Mat lumped_mat = lumped->GetMatrix(); + + pcms::Real trace = 0.0; + for (int v = 0; v < nverts; ++v) { + const PetscInt row = static_cast(gids[v]); + pcms::Real row_sum = 0.0; + for (int u = 0; u < nverts; ++u) { + PetscInt r = row; + PetscInt c = static_cast(gids[u]); + PetscScalar val = 0.0; + MatGetValues(consistent_mat, 1, &r, 1, &c, &val); + row_sum += static_cast(val); + + PetscScalar lumped_val = 0.0; + MatGetValues(lumped_mat, 1, &r, 1, &c, &lumped_val); + CAPTURE(v, u, row_sum, lumped_val); + if (r == c) { + trace += static_cast(lumped_val); + } else { + CHECK(static_cast(lumped_val) == + Catch::Approx(0.0).margin(1e-14)); + } + } + PetscInt r = row; + PetscScalar diag = 0.0; + MatGetValues(lumped_mat, 1, &r, 1, &r, &diag); + CAPTURE(v, row_sum, diag); + CHECK(static_cast(diag) == + Catch::Approx(row_sum).epsilon(1e-12)); + } + + // Trace == area of domain == 1. + CHECK(trace == Catch::Approx(1.0).epsilon(1e-10)); + // matrices are owned by the integrators; no MatDestroy here. +} + +TEST_CASE("OmegaHMassIntegrator: P0 lumped equals consistent", + "[mass_integrator]") +{ + // P0 basis functions have disjoint support, so the consistent mass matrix + // is already diagonal (M_ee = area(e)) and lumping is a no-op. + Omega_h::Library lib; + auto mesh = BuildUnitSquare(lib); + + auto space = pcms::LagrangeFunctionSpace::FromMesh( + mesh, 0, 1, pcms::CoordinateSystem::Cartesian, "global", + pcms::LagrangeFunctionSpace::Backend::OmegaH); + + auto consistent = pcms::BuildOmegaHMassIntegrator(*space); + auto lumped = + pcms::BuildOmegaHMassIntegrator(*space, pcms::MassMatrixType::Lumped); + CHECK(consistent->IsDiagonal()); + CHECK(lumped->IsDiagonal()); + Mat consistent_mat = consistent->GetMatrix(); + Mat lumped_mat = lumped->GetMatrix(); + + const auto elem_areas = Omega_h::measure_elements_real(&mesh); + const auto elem_areas_h = Omega_h::HostRead(elem_areas); + const auto layout = + std::dynamic_pointer_cast( + space->GetLayout()); + const auto& gids = layout->GetGidsHost(); + + for (int e = 0; e < mesh.nelems(); ++e) { + PetscInt r = static_cast(gids[e]); + PetscScalar consistent_val = 0.0; + PetscScalar lumped_val = 0.0; + MatGetValues(consistent_mat, 1, &r, 1, &r, &consistent_val); + MatGetValues(lumped_mat, 1, &r, 1, &r, &lumped_val); + CAPTURE(e, elem_areas_h[e]); + CHECK(static_cast(consistent_val) == + Catch::Approx(elem_areas_h[e]).epsilon(1e-12)); + CHECK(static_cast(lumped_val) == + Catch::Approx(elem_areas_h[e]).epsilon(1e-12)); + } + // matrices are owned by the integrators; no MatDestroy here. +} + TEST_CASE("OmegaHMassIntegrator: rejects invalid layouts", "[mass_integrator]") { Omega_h::Library lib; diff --git a/test/test_omega_h_mc_rhs_integrator.cpp b/test/test_omega_h_mc_rhs_integrator.cpp index 43172fce..0cc33008 100644 --- a/test/test_omega_h_mc_rhs_integrator.cpp +++ b/test/test_omega_h_mc_rhs_integrator.cpp @@ -1,5 +1,6 @@ #include #include +#include #include #include #include @@ -151,6 +152,8 @@ TEST_CASE("OmegaHControlVariateProjection: exact for fields in the target " // For a globally linear field the control variate (target-space interpolant // of the source field) equals the source field, so the sampled residual is // identically zero and the projection is exact regardless of sample count. + // That holds for any invertible mass matrix (delta = M^-1 0 = 0), so the + // lumped variant must be exact here too. Omega_h::Library lib; auto source_mesh = BuildUnitSquare(lib, 1); auto target_mesh = BuildUnitSquare(lib, 0); @@ -163,9 +166,13 @@ TEST_CASE("OmegaHControlVariateProjection: exact for fields in the target " source_field, OMEGA_H_LAMBDA(pcms::Real x, pcms::Real y) { return 3.0 * x - y + 0.25; }); + const auto mass_type = + GENERATE(pcms::MassMatrixType::Consistent, pcms::MassMatrixType::Lumped); + CAPTURE(static_cast(mass_type)); + pcms::OmegaHControlVariateProjection projection( *source_space, *target_space, /*samples_per_element=*/4, - pcms::MonteCarloSampling::UniformRandom); + pcms::MonteCarloSampling::UniformRandom, /*seed=*/8675309, mass_type); projection.Apply(source_field, target_field); const auto values =