initial lumped mass matrix implementation

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,

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_<MassMatrixType>(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_<OmegaHConservativeProjection>(m, "OmegaHConservativeProjection")
     .def(py::init([](const FunctionSpace& source_space,
-                     const FunctionSpace& target_space) {
-           return std::make_unique<OmegaHConservativeProjection>(source_space,
-                                                                 target_space);
+                     const FunctionSpace& target_space,
+                     MassMatrixType mass_matrix_type) {
+           return std::make_unique<OmegaHConservativeProjection>(
+             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<OmegaHControlVariateProjection>(
-             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.")

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
 )
 
 

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;
 };
 

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_);

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

src/pcms/transfer/omega_h_conservative_projection.cpp

@@ -41,7 +41,8 @@ void CheckApplyCompatible(const Field<Real>& source, const Field<Real>& 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<const OmegaHLagrangeLayout>(
       source_space.GetLayout())),
     target_layout_(std::dynamic_pointer_cast<const OmegaHLagrangeLayout>(
@@ -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<GalerkinProjectionSolver>(mass_integrator,
                                                        *rhs_integrator_);
   target_values_ = Kokkos::View<Real**, DeviceMemorySpace>(

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 <Kokkos_Core.hpp>
 #include <memory>
@@ -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<Real>
 {
 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;

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<const OmegaHLagrangeLayout>(
       target_space.GetLayout())),
     rhs_integrator_(std::make_unique<OmegaHMonteCarloRHSIntegrator>(
@@ -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<GalerkinProjectionSolver>(mass_integrator,
                                                        *rhs_integrator_);
 }

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 <cstdint>
@@ -37,11 +38,11 @@ class OmegaHMonteCarloRHSIntegrator;
 class OmegaHControlVariateProjection : public TransferOperator<Real>
 {
 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;

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<const LO*, DeviceMemorySpace>& global_to_local,
-  const Kokkos::View<PetscInt*, DeviceMemorySpace>& coo_rows,
+  bool lumped, const Kokkos::View<PetscInt*, DeviceMemorySpace>& coo_rows,
   const Kokkos::View<PetscInt*, DeviceMemorySpace>& coo_cols)
 {
   Kokkos::parallel_for(
@@ -57,25 +60,28 @@ 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<PetscInt>(global_to_local(verts[i]));
-          coo_cols(idx) = static_cast<PetscInt>(global_to_local(verts[j]));
+          const auto row = static_cast<PetscInt>(global_to_local(verts[i]));
+          coo_rows(idx) = row;
+          coo_cols(idx) =
+            lumped ? row : static_cast<PetscInt>(global_to_local(verts[j]));
         }
       }
     });
 }
 
 } // namespace
 
-OmegaHMassIntegrator::OmegaHMassIntegrator(const FunctionSpace& target_space)
+OmegaHMassIntegrator::OmegaHMassIntegrator(const FunctionSpace& target_space,
+                                           MassMatrixType mass_type)
   : OmegaHMassIntegrator(std::dynamic_pointer_cast<const OmegaHLagrangeLayout>(
                            target_space.GetLayout()),
-                         target_space.GetCoordinateSystem())
+                         target_space.GetCoordinateSystem(), mass_type)
 {
 }
 
 OmegaHMassIntegrator::OmegaHMassIntegrator(
   std::shared_ptr<const OmegaHLagrangeLayout> 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<PetscInt>(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<PetscInt>(nelems) * 9;
   Kokkos::View<PetscInt*, DeviceMemorySpace> coo_rows("mass_coo_rows", nnz);
   Kokkos::View<PetscInt*, DeviceMemorySpace> 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<BilinearFormIntegrator> BuildOmegaHMassIntegrator(
-  const FunctionSpace& target_space)
+  const FunctionSpace& target_space, MassMatrixType mass_type)
 {
-  return std::make_unique<OmegaHMassIntegrator>(target_space);
+  return std::make_unique<OmegaHMassIntegrator>(target_space, mass_type);
 }
 
 } // namespace pcms

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 <memory>
 #include <vector>
 
@@ -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<const OmegaHLagrangeLayout> 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<BilinearFormIntegrator> BuildOmegaHMassIntegrator(
-  const FunctionSpace& target_space);
+  const FunctionSpace& target_space,
+  MassMatrixType mass_type = MassMatrixType::Consistent);
 
 } // namespace pcms