stochastic galerkin projection routines by abhiyanpaudel · Pull Request #300 · SCOREC/pcms

abhiyanpaudel · 2026-03-20T04:06:05Z

This PR adds the stochastic Galerkin projection workflow. Currently, this works for 2D linear elements. A follow-up cleanup would be to remove the current file-based Sobol input path and add a function for Sobol sampling.

Copilot

Pull request overview

This PR introduces a stochastic (Monte Carlo / Sobol-sampled) workflow for assembling the Galerkin projection load vector on 2D linear triangular meshes, alongside renaming the existing mesh-intersection (exact) projection APIs for clarity.

Changes:

Add Monte Carlo sampling + point-localization utilities and a Monte Carlo load-vector integrator.
Split/rename mesh-intersection routines to *MI variants and add *MC solver/RHS variants.
Add unit tests covering sampling/mapping/localization/field-evaluation and end-to-end Monte Carlo projection.

Reviewed changes

Copilot reviewed 13 out of 13 changed files in this pull request and generated 11 comments.

Show a summary per file

File	Description
test/test_monte_carlo_projection_functions.cpp	New tests for Sobol file read, barycentric sampling, mapping, localization, and field evaluation pipeline.
test/test_monte_carlo_field_transfer.cpp	New end-to-end PETSc-enabled test for Monte Carlo projection (constant/linear fields).
test/test_mesh_intersection_field_transfer.cpp	Update to renamed mesh-intersection projection entrypoint (`solveGalerkinProjectionMI`).
test/test_load_vector.cpp	Update to renamed mesh-intersection load-vector builder (`buildLoadVectorMI`).
test/CMakeLists.txt	Register new test sources (Omega_h + PETSc gates).
src/pcms/transfer/load_vector_integrator_mi.cpp	Rename mesh-intersection integrator entrypoint to `buildLoadVectorMI`.
src/pcms/transfer/load_vector_integrator_mc.cpp	New Monte Carlo sampling, mapping, localization, field-eval, and `buildLoadVectorMC`.
src/pcms/transfer/load_vector_integrator.hpp	Add MC API declarations + `SamplingMethod`; rename MI entrypoint declaration.
src/pcms/transfer/conservative_projection_solver.hpp	Add `solveGalerkinProjectionMI/MC` and `computeRhsVectorMI/MC` APIs + updated documentation.
src/pcms/transfer/conservative_projection_solver.cpp	Implement MI/MC projection solves and RHS-vector helpers using PETSc.
src/pcms/transfer/calculate_load_vector.hpp	Rename MI assembler to `calculateLoadVectorMI`; add MC assembler API.
src/pcms/transfer/calculate_load_vector.cpp	Implement `calculateLoadVectorMI` and new `calculateLoadVectorMC`.
src/pcms/transfer/CMakeLists.txt	Build both MI + MC integrator sources; keep PETSc pieces behind `PCMS_ENABLE_PETSC`.

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

Copilot · 2026-03-20T04:13:14Z

test/test_monte_carlo_projection_functions.cpp

+  const std::string filename = "tmp_sobol_barycentric_samples.txt";
+  std::ofstream out(filename);
+  REQUIRE(out.is_open());
+  out << "l0 l1 l2\n";
+  out << "1.0 0.0 0.0\n";
+  out << "0.0 1.0 0.0\n";
+  out << "0.2 0.3 0.5\n";
+  out.close();
+


This test uses a fixed temporary filename (tmp_sobol_barycentric_samples.txt). If tests are executed in parallel (or rerun after a crash), this can cause collisions or stale data. Consider using a unique temp path and/or an RAII cleanup guard so the file is removed even when assertions fail.

Copilot · 2026-03-20T04:13:14Z

src/pcms/transfer/load_vector_integrator_mc.cpp

+
+  int nx = edge_length[0] * Kokkos::sqrt(mesh.nelems());
+  int ny = edge_length[1] * Kokkos::sqrt(mesh.nelems());
+


localize_points_in_mesh computes nx/ny by casting a floating expression to int. For small bounding boxes (edge length < 1) this can truncate to 0, which later causes division-by-zero in UniformGrid (divisions used as a divisor). Clamp nx/ny to at least 1 (and consider using Omega_h::LO/std::max).

Suggested change

if (nx < 1) nx = 1;

if (ny < 1) ny = 1;

Copilot · 2026-03-20T04:13:15Z

src/pcms/transfer/load_vector_integrator.hpp

+#include <Omega_h_mesh.hpp>
+#include <fstream>
+#include <iostream>
+#include <sstream>


load_vector_integrator.hpp declares APIs using std::string, but it does not include <string>. Please add the missing include to make the header self-contained, and consider removing heavy/unused includes like <iostream> (include it in the .cpp that needs it instead).

Suggested change

#include <sstream>

#include <sstream>

#include <string>

Copilot · 2026-03-20T04:13:15Z

src/pcms/transfer/load_vector_integrator_mc.cpp

+  const auto npoints = points.extent(0);
+


const auto npoints = points.extent(0); is unused, which will trigger warnings in this TU. Please remove it or use it for validation (e.g., early return on zero points).

Suggested change

const auto npoints = points.extent(0);

Copilot · 2026-03-20T04:13:15Z

src/pcms/transfer/load_vector_integrator_mc.cpp

+  const auto& npoints = results.extent(0);
+
+  const auto& faces2nodes = mesh.ask_down(Omega_h::FACE, Omega_h::VERT).ab2b;
+  const auto& coordinates = mesh.coords();


const auto& coordinates = mesh.coords(); is unused in evaluate_field_from_point_localization. Please remove it to avoid warnings and keep the function focused.

Suggested change

const auto& coordinates = mesh.coords();

Copilot · 2026-03-20T04:13:16Z

src/pcms/transfer/conservative_projection_solver.cpp

+
+    Vec vec;
+    PetscErrorCode ierr;
+    ierr =
+      calculateLoadVectorMC(target_mesh, field_values_at_points,
+                            npoints_each_tri, method, &vec, sobol_filename);
+    CHKERRABORT(PETSC_COMM_WORLD, ierr);
+
+    auto rhsvector = vecToOmegaHReals(vec);
+    ierr = VecDestroy(&vec);
+    CHKERRABORT(PETSC_COMM_WORLD, ierr);
+
+    return rhsvector;
+}	
+


computeRhsVectorMC has inconsistent indentation and a trailing tab after the closing brace, which will fight auto-formatting and make diffs noisy. Please reformat this block to match the rest of the file.

Suggested change

Vec vec;

PetscErrorCode ierr;

ierr =

calculateLoadVectorMC(target_mesh, field_values_at_points,

npoints_each_tri, method, &vec, sobol_filename);

CHKERRABORT(PETSC_COMM_WORLD, ierr);

auto rhsvector = vecToOmegaHReals(vec);

ierr = VecDestroy(&vec);

CHKERRABORT(PETSC_COMM_WORLD, ierr);

return rhsvector;

}

Vec vec;

PetscErrorCode ierr;

ierr =

calculateLoadVectorMC(target_mesh, field_values_at_points,

npoints_each_tri, method, &vec, sobol_filename);

CHKERRABORT(PETSC_COMM_WORLD, ierr);

auto rhsvector = vecToOmegaHReals(vec);

ierr = VecDestroy(&vec);

CHKERRABORT(PETSC_COMM_WORLD, ierr);

return rhsvector;

}

Copilot · 2026-03-20T04:13:16Z

test/test_monte_carlo_field_transfer.cpp

+#include <pcms/transfer/conservative_projection_solver.hpp>
+#include <pcms/transfer/load_vector_integrator.hpp>
+
+#include <fstream>


This test uses std::function (see evaluate_field_at_barycentric_samples) but does not include <functional>. Please add the missing include to avoid relying on transitive headers.

Suggested change

#include <fstream>

#include <fstream>

#include <functional>

Copilot · 2026-03-20T04:13:16Z

test/test_monte_carlo_field_transfer.cpp

+  auto bary = make_sobol_barycentric_samples();
+  const int npoints_each_tri = bary.extent(0);
+  const std::string sobol_filename = "tmp_mc_barycentric_samples.txt";
+  write_barycentric_samples_file(sobol_filename, bary);
+


The test writes tmp_mc_barycentric_samples.txt but never deletes it. This can pollute the working directory and can also cause collisions when tests run concurrently. Add cleanup (ideally via an RAII guard) and consider using a unique filename per test run.

Copilot · 2026-03-20T04:13:16Z

test/CMakeLists.txt

      test_interpolation_class.cpp
      test_mesh_geometry.cpp
+	  test_monte_carlo_projection_functions.cpp
      test_omega_h_field2_outofbounds.cpp)


This added entry is indented with a tab, which is inconsistent with the surrounding CMake formatting and may fail format checks. Replace the tab with spaces to match the rest of the file.

Copilot · 2026-03-20T04:13:17Z

test/CMakeLists.txt

    list(APPEND PCMS_UNIT_TEST_SOURCES
-      test_mesh_intersection_field_transfer.cpp)
+      test_mesh_intersection_field_transfer.cpp
+  	  test_monte_carlo_field_transfer.cpp)
  endif()


This list append adds an entry with tab indentation. Use spaces to match the existing CMake formatting so automated formatters don't keep rewriting this line.

Sichao25

Looks good overall. Just commented a few minor issues.

Sichao25 · 2026-03-23T17:27:13Z

src/pcms/transfer/calculate_load_vector.cpp

 namespace pcms
 {
-PetscErrorCode calculateLoadVector(Omega_h::Mesh& target_mesh,
+PetscErrorCode calculateLoadVectorMI(Omega_h::Mesh& target_mesh,


If I haven't missed anything, the calculateLoadVectorMI and calculateLoadVectorMC share roughly 90% of their code, there may be room to refactor and eliminate the duplication.

Sichao25 · 2026-03-23T17:32:51Z

src/pcms/transfer/calculate_load_vector.cpp

+    buildLoadVectorMI(target_mesh, source_mesh, intersection, source_values);
+
+  auto hostElmLoadVector = Kokkos::create_mirror_view(elmLoadVector);
+  Kokkos::deep_copy(hostElmLoadVector, elmLoadVector);


There is a create_mirror_view_and_copy function.

Sichao25 · 2026-03-23T17:36:43Z

src/pcms/transfer/load_vector_integrator_mc.cpp

+{
+  std::ifstream infile(file_path);
+  if (!infile) {
+    throw std::runtime_error("Could not open sobol sample file");


We already have pcms_error.

abhiyan123 added 30 commits March 17, 2026 13:43

add mc conservative projection file

4bf033a

add mc helper functions declaration

d2eb0af

change function name from buildLoadVector to buildLoadVectorMI

8f5148a

fix the load vector integrator name

9094916

add mc load vector build declaration

a95460b

add header

e0b9fcd

delete pcms_mc_proj_utils.hpp

26f1a4f

rename pcms_mc_proj_utils.cpp -> load_vector_integrator_mc.cpp

12c50c7

rename load_vector_integrator.cpp -> load_vector_integrator_mi.cpp

21fdaba

add mc load vector build

ab31a90

add build load vector with mc approach

54213bb

change the parameter sequence

ecd8234

add mc solver function

56c0160

change parameter sequence

e489bea

update variants of galerkin projection solve

17644e3

add namespace

d542524

update filenames in CMakeLists.txt

aea8cb6

update solveGalerkinProjection to solveGalerkinProjectionMI

15897c1

update buildLoadVector to buildLoadVectorMI

5643ec2

add mc test functions

6356c81

update mc test file in CMakeLists.txt

bb20610

add semicolon

6d6e7e1

correct the class name

65a09c9

update variable

4f664bc

correct class name

af6e728

change tri_id to element_id

a6207e7

change tri_id to element_id

d51b3cd

minor changes

73d2143

fix the function name

4e626c9

update header

5dca135

abhiyan123 added 14 commits March 19, 2026 13:42

update header

e0915d1

update function

0310306

minor changes

4556bdc

add test case for monte carlo field transfer

24eccd5

add filename in CMakeLists.txt

083144c

update header

30785f1

update solveGalerkinProjectionMC parameters

82f2fa8

update args in a solveGalerkinProjectionMC function

7a9785d

create a host view instead of device

5a3d44a

update the API docs for MC projection

80f8aba

add sobol samples

e0c310d

added more sobol barycentric points

86285ff

change triangular element to just element

a92a3d1

minor changes

228d58a

abhiyanpaudel had a problem deploying to asics March 20, 2026 04:06 — with GitHub Actions Failure

abhiyanpaudel requested review from Copilot and jacobmerson March 20, 2026 04:06

Copilot started reviewing on behalf of abhiyanpaudel March 20, 2026 04:08 View session

Copilot AI reviewed Mar 20, 2026

View reviewed changes

abhiyanpaudel had a problem deploying to asics March 23, 2026 17:51 — with GitHub Actions Failure

Sichao25 reviewed Mar 23, 2026

View reviewed changes


		int nx = edge_length[0] * Kokkos::sqrt(mesh.nelems());
		int ny = edge_length[1] * Kokkos::sqrt(mesh.nelems());

Conversation

abhiyanpaudel commented Mar 20, 2026

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Reviewed changes

Uh oh!

Copilot AI Mar 20, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Mar 20, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Mar 20, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Mar 20, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Mar 20, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Mar 20, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Mar 20, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Mar 20, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Mar 20, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Mar 20, 2026

Choose a reason for hiding this comment

Uh oh!

Sichao25 left a comment

Choose a reason for hiding this comment

Uh oh!

Sichao25 Mar 23, 2026

Choose a reason for hiding this comment

Uh oh!

Sichao25 Mar 23, 2026

Choose a reason for hiding this comment

Uh oh!

Sichao25 Mar 23, 2026

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants