Conservative field transfer (integrator + Monte-Carlo projection) by jacobmerson · Pull Request #335 · SCOREC/pcms

jacobmerson · 2026-07-12T16:31:22Z

Adds conservative field-transfer support built on the integrator-based (Galerkin) and Monte-Carlo / control-variate projection operators, replacing the old calculate_mass_matrix / load_vector implementation.

Highlights:

Integrator-based conservative (Galerkin) projection transfer
Monte-Carlo / control-variate conservative projection
DOF holders indexed by a compact GID permutation; coupler separates GID-layout and field messages
Python bindings expose conservative transfer

Base is jacobmerson/various-fixes (#340 ).

Copilot

Pull request overview

This PR refactors conservative field transfer to use reusable integrator abstractions (mass + RHS) and a cached Galerkin projection solver, adds Monte Carlo (control-variate) conservative projection, and updates coupler/GID handling and Python bindings accordingly.

Changes:

Replace legacy calculate_mass_matrix / calculate_load_vector and load-vector integrator paths with BilinearFormIntegrator/LinearFormIntegrator + GalerkinProjectionSolver.
Add Omega_h-specific conservative RHS integrators (mesh-intersection quadrature) and Monte Carlo RHS integrator, plus a control-variate projection operator.
Introduce compact GID permutations and adjust field exchange planning/serialization; expose new transfer operators via Python and add new unit tests + example.

Reviewed changes

Copilot reviewed 53 out of 53 changed files in this pull request and generated 7 comments.

Show a summary per file

File	Description
test/test_omega_h_mc_rhs_integrator.cpp	Adds tests for Monte Carlo RHS integrator and control-variate projection behavior.
test/test_omega_h_mass_integrator.cpp	Adds tests validating Omega_h mass-matrix assembly vs MeshField reference.
test/test_omega_h_intersection_rhs_integrator.cpp	Adds tests for intersection-based conservative RHS integrator correctness.
test/test_omega_h_form_integrator_utils.cpp	Adds tests for polygon vertex dedup/relink helper used in intersection integration.
test/test_mesh_intersection_field_transfer.cpp	Reworks/extends conservative projection tests (P0/P1, reordered/sparse GIDs, conservation).
test/test_load_vector.cpp	Removes legacy load vector test tied to deleted implementation.
test/test_field_exchange_planner.cpp	Adds tests ensuring exchange plans exclude GID headers and headers are inserted correctly.
test/CMakeLists.txt	Registers new tests and updates conditional test lists.
src/pcms/transfer/omega_h_mc_rhs_integrator.hpp	Defines Monte Carlo RHS integrator interface and sampling enum.
src/pcms/transfer/omega_h_mc_rhs_integrator.cpp	Implements Monte Carlo sample generation and COO RHS assembly.
src/pcms/transfer/omega_h_mass_integrator.hpp	Defines Omega_h mass-matrix integrator (P0/P1) and builder.
src/pcms/transfer/omega_h_mass_integrator.cpp	Implements PETSc COO mass-matrix assembly for Omega_h layouts.
src/pcms/transfer/omega_h_intersection_rhs_integrator.hpp	Defines intersection-based conservative RHS integrator and builder.
src/pcms/transfer/omega_h_intersection_rhs_integrator.cpp	Implements mesh-intersection quadrature point generation and RHS assembly.
src/pcms/transfer/omega_h_form_integrator_utils.hpp	Adds shared validation, intersection traversal, quadrature staging, and target basis helpers.
src/pcms/transfer/omega_h_control_variate_projection.hpp	Declares control-variate Monte Carlo conservative projection operator.
src/pcms/transfer/omega_h_control_variate_projection.cpp	Implements control-variate projection workflow (interpolate + residual projection).
src/pcms/transfer/omega_h_conservative_projection.hpp	Updates conservative projection to cache integrators/evaluator/solver (no per-call refactorization).
src/pcms/transfer/omega_h_conservative_projection.cpp	Implements cached projection Apply() using solver + permutation scatter.
src/pcms/transfer/mass_matrix_integrator.hpp	Renames helper to `buildElementMassMatrix` and fixes override signature.
src/pcms/transfer/linear_form_integrator.hpp	Introduces LinearFormIntegrator interface for RHS assembly.
src/pcms/transfer/bilinear_form_integrator.hpp	Introduces BilinearFormIntegrator interface for matrix assembly.
src/pcms/transfer/integration_point_set.hpp	Adds IntegrationPointSet wrapper for ordered integration coordinates.
src/pcms/transfer/conservative_projection_solver.hpp	Refactors solver into cached `GalerkinProjectionSolver` using integrator interfaces.
src/pcms/transfer/conservative_projection_solver.cpp	Implements cached KSP setup and Solve overloads (evaluator-driven and pre-sampled).
src/pcms/transfer/load_vector_integrator.hpp	Removes legacy supermesh load-vector integrator API.
src/pcms/transfer/load_vector_integrator.cpp	Removes legacy load-vector integrator implementation.
src/pcms/transfer/calculate_mass_matrix.hpp	Removes legacy PETSc mass matrix assembly API.
src/pcms/transfer/calculate_mass_matrix.cpp	Removes legacy mass matrix assembly implementation.
src/pcms/transfer/calculate_load_vector.hpp	Removes legacy PETSc load vector assembly API.
src/pcms/transfer/calculate_load_vector.cpp	Removes legacy load vector assembly implementation.
src/pcms/transfer/interpolator.h	Fixes include path for transfer operator header.
src/pcms/transfer/CMakeLists.txt	Wires in new transfer headers/sources and removes deprecated ones.
src/pcms/pythonapi/pythonapi.cpp	Initializes PETSc at import time when PETSc+MeshFields are enabled.
src/pcms/pythonapi/CMakeLists.txt	Links PETSc/MeshFields into Python module when needed for new transfer bindings.
src/pcms/pythonapi/bind_transfer_field.cpp	Exposes conservative and Monte Carlo/control-variate projections + sampling enum to Python.
src/pcms/pythonapi/bind_omega_h.cpp	Adds VTU read helper and fixes mesh construction in parallel VTK read.
src/pcms/pythonapi/bind_field_base.cpp	Improves numpy contiguity handling; exposes FunctionSpace::mesh and Lagrange backend selection.
src/pcms/localization/adj_search.cpp	Adjusts handling of exterior points and caps adaptive radius loop.
src/pcms/field/layout/omega_h_lagrange.h	Adds device GID accessor.
src/pcms/field/layout/omega_h_lagrange.cpp	Builds global-to-local permutation eagerly in constructors; implements GetGids().
src/pcms/field/gid_permutation.hpp	Adds utilities for unmapped/sorted GID permutations and permutation construction.
src/pcms/field/gid_permutation.cpp	Implements permutation utilities and uniqueness checks.
src/pcms/field/field_layout.h	Adds global-to-local permutation API and eager-build helper/caches.
src/pcms/field/field_layout.cpp	Implements permutation cache building and accessors.
src/pcms/field/element_dispatch.h	Adds runtime-to-compile-time FE order dispatch helper.
src/pcms/field/CMakeLists.txt	Adds new field sources/headers to build.
src/pcms/coupler/serializer/xgc.h	Skips unmapped permutation entries during (de)serialization.
src/pcms/coupler/field_serializer.h	Skips unmapped permutation entries during (de)serialization and sorted extraction.
src/pcms/coupler/field_layout_communicator.cpp	Adjusts outgoing GID message layout/size to account for per-message headers.
src/pcms/coupler/field_exchange_planner.h	Clarifies exchange plan offsets/permutation exclude GID headers.
src/pcms/coupler/field_exchange_planner.cpp	Refactors GID message handling to insert/parse headers and support unmapped holders.
examples/python-api-example/example_conservative_transfer.py	Adds end-to-end Python example for conservative + Monte Carlo transfer and mixed-order cases.

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

+          static_cast<PetscScalar>(coeffs(i * ndof + j)) * f;
+      }
+    });
+


jacobmerson · 2026-07-13T05:05:39Z

There are a few GPU issues I'm still sorting out.

…trix / load_vector) Squashed original commits: 0a27b53, 0c0e697

…nsfer Squashed original commits: bad8434, eca15d1, b4af186, 26f35c6, b8baf3b, 6b1d7f2, dba7966, 5baf365, 4b00ac1 (mesh_fields_rhs_integrator churn from 26f35c6/b8baf3b dropped)

Squashed original commits: 844dc1c, b4a4545, 629419f (sobol_sequence churn dropped)

…y buffers Squashed original commits: e7488c6, 4e68c80

…fter #308) Squashed original commit: 75a3045

Currently the petsc-based solver is only needed for conservative projection which requires MeshFields. So petsc is only linked when PCMS_ENABLE_PETSC and PCMS_ENABLE_MESHFIELDS

Add a shared, format-agnostic GID-to-local permutation and use it so the transfer operators handle sparse or reordered global IDs correctly. This allows us to reuse the global-to-local mapping for things outside of the serialization such as setting up the petsc arrays.

…olders Split the GID (layout) message from the field (data) message and reuse the field GID-permutation primitives so exchanges handle overlap masks and reordered global IDs without corrupting payloads.

The python pcms module links libpetsc.so dynamically when python_api is ON, so PETSc lib dir must be on LD_LIBRARY_PATH for import to succeed.

clang-tidy bugprone-implicit-widening-of-multiplication-result (treated as error) flagged int*int products used to size 64-bit Views in the intersection/MC RHS integrators and the field-exchange-planner test. Perform the multiplication in size_t via an explicit cast.

Forward-declare the RHS integrators

We had kokkos lambdas in private functions which does not compile.

First pass is to just make a host mirror. Long term, we need to ask PETSc team if there is a better way to handle this.

jacobmerson · 2026-07-13T07:28:02Z

GPU builds are now working/tests passing (at least on othello)

Sichao25

Overall looks good to me.

Sichao25 · 2026-07-21T19:22:16Z

+// concrete integrator implementation. It owns the coordinate buffer; callers
+// receive a non-owning CoordinateView from GetCoordinates().
+template <typename MemorySpace>
+class IntegrationPointSet


This looks like a thin wrapper that just owns the coordinates and provides a coordinate view. Do we need this instead of just a plain Kokkos view?

Yeah, we can drop this. I originally was thinking of using that as an abstract base, but the use-case I had in my mind was more theoretical than executed in the codebase.

Sichao25 · 2026-07-21T19:27:37Z

-        if (owned[i]) {
+        // A negative permutation entry marks a holder outside the exchange
+        // (non-owned, or owned but outside the overlap region); it has no slot.
+        if (owned[i] && permutation[i] >= 0) {


The permutation is constructed when an holder is owned and in the overlap domain. If so, does the permutation itself already include the 'owned' info? Is there a use case where we'd need to check both?

I believe that's correct. with the current version of the permutation, we explicitly have negative values for any non-owned entities. We probably need to formalize this somehow. At least through comments, but we may want to add an assert in the loop and add test cases that check that the serialized array is negative for non-owning entries.

Sichao25 · 2026-07-21T19:49:05Z

+      Omega_h::Real area =
+        Kokkos::fabs(Omega_h::triangle_area_from_basis(basis));
+
+      const double eps_area = abs_tol + rel_tol * poly_area;


This might be irrelevant to the implementation here, but I'm wondering how abs_tol and rel_tol are determined? Is there a rule for this? I want to double-check because I remember a mismatch between the fuzz value （already be removed from point search and epsilon from meshfields causing a bug in some cases.

At the moment, I believe these are hard-coded, but this code will eventually need to be made more robust. @abhiyan123 may have more details on why those triangles are skipped, but I'm pretty sure it was to avoid degeneracies that can happen with slivers, etc.

Long term, we should look into simulation of simplicity https://dl.acm.org/doi/10.1145/77635.77639 , or other methods that can help to avoid these cases without tolerances.

Sichao25 · 2026-07-21T20:31:56Z

The OmegaHMassIntegrator test failed. I will take a look.

Sichao25 · 2026-07-22T21:14:27Z

I've added a few commits addressing the test failures and the comments. If it looks good, it's ready to merge.

jacobmerson mentioned this pull request Jul 12, 2026

stochastic galerkin projection routines #300

Closed

jacobmerson requested a review from Copilot July 12, 2026 16:42

Copilot started reviewing on behalf of jacobmerson July 12, 2026 16:43 View session

jacobmerson mentioned this pull request Jul 12, 2026

Field/Function split: FieldFactory + two-level naming, remove AddLayout #336

Open

Copilot AI reviewed Jul 12, 2026

View reviewed changes

jacobmerson force-pushed the jacobmerson/conservative-transfer branch from 7dfbf1e to 5bda307 Compare July 12, 2026 18:29

jacobmerson force-pushed the jacobmerson/various-fixes branch from b4d81b3 to 6c13b7c Compare July 12, 2026 18:29

jacobmerson added 15 commits July 13, 2026 03:07

transfer: remove old-style conservative projection (calculate_mass_ma…

b53cc15

…trix / load_vector) Squashed original commits: 0a27b53, 0c0e697

transfer: add integrator-based conservative (Galerkin) projection tra…

0bb00c9

…nsfer Squashed original commits: bad8434, eca15d1, b4af186, 26f35c6, b8baf3b, 6b1d7f2, dba7966, 5baf365, 4b00ac1 (mesh_fields_rhs_integrator churn from 26f35c6/b8baf3b dropped)

transfer: add Monte-Carlo / control-variate conservative projection

08d6143

Squashed original commits: 844dc1c, b4a4545, 629419f (sobol_sequence churn dropped)

pythonapi: expose conservative transfer; flatten DOF fill for ND nump…

2e58ed9

…y buffers Squashed original commits: e7488c6, 4e68c80

transfer: ad-hoc adjacency-search snapping fixes (TEMPORARY; remove a…

28aac51

…fter #308) Squashed original commit: 75a3045

fix petsc include guards in python bindings

328e674

Currently the petsc-based solver is only needed for conservative projection which requires MeshFields. So petsc is only linked when PCMS_ENABLE_PETSC and PCMS_ENABLE_MESHFIELDS

coupler: separate GID-layout and field messages; skip non-exchanged h…

af9c175

…olders Split the GID (layout) message from the field (data) message and reuse the field GID-permutation primitives so exchanges handle overlap masks and reordered global IDs without corrupting payloads.

clang-format field_exchange_planner.cpp

d1f768d

ci: add PETSc lib dir to LD_LIBRARY_PATH for python API tests

ab51122

The python pcms module links libpetsc.so dynamically when python_api is ON, so PETSc lib dir must be on LD_LIBRARY_PATH for import to succeed.

transfer: keep PETSc/MeshField backends out of public headers

9c1ff5c

Forward-declare the RHS integrators

Fix GPU build for integrators

a44f8c6

We had kokkos lambdas in private functions which does not compile.

extract out more kokkos lambdas to addressable functions

b473bb3

Petsc functions to set the sparsity patterns take CPU pointers

c72ccad

First pass is to just make a host mirror. Long term, we need to ask PETSc team if there is a better way to handle this.

jacobmerson force-pushed the jacobmerson/various-fixes branch from 6c13b7c to 9559d9d Compare July 13, 2026 07:11

jacobmerson force-pushed the jacobmerson/conservative-transfer branch from 2730ec6 to c72ccad Compare July 13, 2026 07:11

jacobmerson changed the base branch from jacobmerson/various-fixes to develop July 21, 2026 09:14

Merge branch 'develop' into jacobmerson/conservative-transfer

95fdba3

jacobmerson requested a deployment to asics July 21, 2026 09:20 — with GitHub Actions Waiting

Reformat with clang-format 18.1.3 to match CI

30d99b9

jacobmerson had a problem deploying to asics July 21, 2026 14:17 — with GitHub Actions Failure

Sichao25 reviewed Jul 21, 2026

View reviewed changes

debug device memory access in mass integrator test

c640350

Sichao25 requested a deployment to asics July 21, 2026 22:41 — with GitHub Actions Waiting

remove integration_point_set

d0831b7

Sichao25 temporarily deployed to asics July 22, 2026 01:39 — with GitHub Actions Inactive

remove owned check since it is included by negative permutation

3a0c843

Sichao25 deployed to asics July 22, 2026 17:19 — with GitHub Actions Active

Uh oh!

Conversation

jacobmerson commented Jul 12, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Reviewed changes

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

jacobmerson commented Jul 13, 2026

Uh oh!

jacobmerson commented Jul 13, 2026

Uh oh!

Sichao25 left a comment

Choose a reason for hiding this comment

Uh oh!

Sichao25 Jul 21, 2026

Choose a reason for hiding this comment

Uh oh!

jacobmerson Jul 21, 2026

Choose a reason for hiding this comment

Uh oh!

Sichao25 Jul 21, 2026

Choose a reason for hiding this comment

Uh oh!

jacobmerson Jul 21, 2026

Choose a reason for hiding this comment

Uh oh!

Sichao25 Jul 21, 2026

Choose a reason for hiding this comment

Uh oh!

jacobmerson Jul 21, 2026

Choose a reason for hiding this comment

Uh oh!

Sichao25 commented Jul 21, 2026

Uh oh!

Sichao25 commented Jul 22, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

jacobmerson commented Jul 12, 2026 •

edited

Loading