initial lumped mass matrix implementation - #353
Conversation
|
fixes #354 |
Fuad-HH
left a comment
Fuad-HH
left a comment
There was a problem hiding this comment.
Everything looks good to me. Here are some comments which are NOT very important.
| copy.h | ||
| transfer_operator.hpp | ||
| mass_matrix_type.hpp | ||
| monte_carlo_sampling.hpp |
There was a problem hiding this comment.
Is this related to this PR or just an older bug?
There was a problem hiding this comment.
This is an old bug that the header was missing.
There was a problem hiding this comment.
Should we add a new header file just for this enum class? I am not sure which existing header will be the best to move it to though.
There was a problem hiding this comment.
Probably not needed.
| max_err = std::max(max_err, std::abs(target_values[i] - expected)); | ||
| } | ||
| CHECK(max_err > 1e-3); | ||
| } |
There was a problem hiding this comment.
We should also assert the boundedness of the interpolated field for this lumped case.
|
I guess, the iterative method described with equation 36 in the mesh intersection paper for minimal diffusive bounded interpolant is also important and can be built on this. However, as we discussed before, diffusiveness is not an urgent issue for XGC-Degas2 coupling. |
There was a problem hiding this comment.
Pull request overview
Adds a selectable lumped (row-sum diagonal) mass-matrix formulation to the Omega_h Galerkin projection operators, enabling diagonal mass solves (and associated bounded/conservative workflows) while preserving the existing consistent-mass behavior by default. This threads a new MassMatrixType option through C++ and Python APIs, updates solver setup for diagonal systems, and expands tests/examples to cover the new mode.
Changes:
- Introduces
pcms::MassMatrixTypeand plumbs it throughOmegaHMassIntegrator,OmegaHConservativeProjection, andOmegaHControlVariateProjection. - Adds a diagonal-matrix fast path in
GalerkinProjectionSolver(KSPPREONLY + Jacobi) based on a newBilinearFormIntegrator::IsDiagonal()hook. - Extends unit tests and the Python example/bindings to exercise lumped vs consistent mass behavior.
Reviewed changes
Copilot reviewed 15 out of 15 changed files in this pull request and generated 1 comment.
Show a summary per file
| File | Description |
|---|---|
| test/test_omega_h_mc_rhs_integrator.cpp | Parameterizes an exactness test over consistent vs lumped mass in the control-variate projection. |
| test/test_omega_h_mass_integrator.cpp | Adds coverage verifying lumped diagonal equals consistent row sums and that P0 lumping is a no-op. |
| test/test_mesh_intersection_field_transfer.cpp | Adds lumped-vs-consistent conservation coverage and a targeted test demonstrating lumped non-reproduction of linears. |
| src/pcms/transfer/omega_h_mass_integrator.hpp | Adds MassMatrixType option and IsDiagonal() reporting for mass integrators. |
| src/pcms/transfer/omega_h_mass_integrator.cpp | Implements row-sum lumping via COO diagonal mapping and tracks diagonal-ness. |
| src/pcms/transfer/omega_h_control_variate_projection.hpp | Extends constructor to accept MassMatrixType (default consistent). |
| src/pcms/transfer/omega_h_control_variate_projection.cpp | Passes MassMatrixType through to the mass integrator used by the solver. |
| src/pcms/transfer/omega_h_conservative_projection.hpp | Extends constructor to accept MassMatrixType (default consistent) and documents behavior impact. |
| src/pcms/transfer/omega_h_conservative_projection.cpp | Passes MassMatrixType through to the mass integrator used by the solver. |
| src/pcms/transfer/mass_matrix_type.hpp | New enum defining consistent vs row-sum lumped mass-matrix formulations. |
| src/pcms/transfer/conservative_projection_solver.cpp | Selects an exact diagonal solve configuration when the mass matrix is diagonal. |
| src/pcms/transfer/CMakeLists.txt | Installs the new header(s) in the transfer file set. |
| src/pcms/transfer/bilinear_form_integrator.hpp | Adds IsDiagonal() virtual hook (default false) for solver optimization. |
| src/pcms/pythonapi/bind_transfer_field.cpp | Exposes MassMatrixType to Python and adds it to projection constructors with defaults. |
| examples/python-api-example/example_conservative_transfer.py | Demonstrates using the lumped mass option from Python and contrasts expected error behavior. |
💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.
| 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); |
| 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); | ||
| } |
There was a problem hiding this comment.
Please replace the PETSC_COMM_WORLDs with PETSC_COMM_SELF s here.
Agree. I did a bit of work on this too. Once we get the lumped mass matrix version in place, I can open a PR with those additions. It's directly additive on this since it requires use of the lumped mass. |
Bounded interpolant in conservative field interpolation.