Fix flat array to rank-2 array conversion

src/pcms/transfer/CMakeLists.txt

@@ -21,7 +21,8 @@ set(PCMS_FIELD_TRANSFER_HEADERS
 
 set(PCMS_FIELD_TRANSFER_SOURCES
     adj_search.cpp
-    load_vector_integrator.cpp
+    load_vector_integrator_mi.cpp
+    load_vector_integrator_mc.cpp
     mesh_intersection.cpp
     mls_interpolation.cpp
     interpolation_base.cpp)

src/pcms/transfer/calculate_load_vector.cpp

@@ -5,7 +5,7 @@
 
 namespace pcms
 {
-PetscErrorCode calculateLoadVector(Omega_h::Mesh& target_mesh,
+PetscErrorCode calculateLoadVectorMI(Omega_h::Mesh& target_mesh,
                                    Omega_h::Mesh& source_mesh,
                                    const IntersectionResults& intersection,
                                    const Omega_h::Reals& source_values,
@@ -21,7 +21,52 @@ PetscErrorCode calculateLoadVector(Omega_h::Mesh& target_mesh,
 
   // Fill COO values
   auto elmLoadVector =
-    buildLoadVector(target_mesh, source_mesh, intersection, source_values);
+    buildLoadVectorMI(target_mesh, source_mesh, intersection, source_values);
+
+  auto hostElmLoadVector = Kokkos::create_mirror_view(elmLoadVector);
+  Kokkos::deep_copy(hostElmLoadVector, elmLoadVector);
+  PetscCheck(static_cast<PetscInt>(hostElmLoadVector.extent(0)) == nnz,
+             PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
+             "Element load vector size (%d) does not match COO nnz (%d)",
+             static_cast<PetscInt>(hostElmLoadVector.extent(0)), nnz);
+
+  for (PetscInt e = 0; e < nnz; ++e) {
+    coo_vals[e] = hostElmLoadVector(e);
+  }
+
+  // create vector with preallocated COO structure
+  Vec vec;
+  PetscCall(createSeqVec(PETSC_COMM_WORLD, target_mesh.nverts(), &vec));
+  PetscCall(VecSetPreallocationCOO(vec, nnz, coo_i));
+  PetscCall(VecSetValuesCOO(vec, coo_vals, ADD_VALUES));
+  PetscCall(PetscFree(coo_i));
+  PetscCall(PetscFree(coo_vals));
+
+  if (target_mesh.nelems() < 10) {
+    PetscCall(VecView(vec, PETSC_VIEWER_STDOUT_WORLD));
+  }
+
+  *loadVec_out = vec;
+  PetscFunctionReturn(PETSC_SUCCESS);
+}
+
+PetscErrorCode calculateLoadVectorMC(
+  Omega_h::Mesh& target_mesh, const Omega_h::Reals& field_values_at_points,
+  const int npoints_each_tri, SamplingMethod method, Vec* loadVec_out,
+  const std::string sobol_filename)
+{
+
+  PetscFunctionBeginUser;
+  PetscInt nnz = 0;
+  PetscInt* coo_i = nullptr;
+  PetscCall(build_linear_triangle_coo_rows(target_mesh, &coo_i, &nnz));
+  PetscScalar* coo_vals = nullptr;
+  PetscCall(PetscMalloc1(nnz, &coo_vals));
+
+  // Fill COO values
+  auto elmLoadVector =
+    buildLoadVectorMC(target_mesh, field_values_at_points, npoints_each_tri,
+                      method, sobol_filename);
 
   auto hostElmLoadVector = Kokkos::create_mirror_view(elmLoadVector);
   Kokkos::deep_copy(hostElmLoadVector, elmLoadVector);

src/pcms/transfer/calculate_load_vector.hpp

@@ -13,14 +13,18 @@
 #define PCMS_TRANSFER_CALCULATE_LOAD_VECTOR_HPP
 #include <Omega_h_mesh.hpp>
 #include <pcms/transfer/mesh_intersection.hpp>
+#include <pcms/transfer/load_vector_integrator.hpp>
 #include <petscvec.h>
 
+namespace pcms
+{
+
 /**
- * @brief Assembles the global load vector.
+ * @brief Assembles the global load vector using mesh intersection method.
  *
  * This function computes the unassembled local load vector contributions for
- * each triangular element in the target mesh using `buildLoadVector()` and then
- * assembles them into a global PETSc vector in COO format.
+ * each triangular element in the target mesh using `buildLoadVectorMI()` and
+ * then assembles them into a global PETSc vector in COO format.
  *
  *
  * @param target_mesh The target Omega_h mesh to which the scalar field is being
@@ -40,20 +44,66 @@
  * @note
  * - Works for 2D linear triangular elements.
  * - Uses COO-style preallocation and insertion into the PETSc vector.
- * - Internally calls `buildLoadVector()` to compute per-element contributions.
+ * - Internally calls `buildLoadVectorMI()` to compute per-element
+ * contributions.
  * - The resulting vector is used as the right-hand side (RHS) in a projection
  * solve.
  *
- * @see buildLoadVector,IntersectionResults
+ * @see buildLoadVectorMI,IntersectionResults
  */
-
-namespace pcms
-{
 // FIXME use PCMS error handling rather than returning a PETSC error code
-PetscErrorCode calculateLoadVector(Omega_h::Mesh& target_mesh,
-                                   Omega_h::Mesh& source_mesh,
-                                   const IntersectionResults& intersection,
-                                   const Omega_h::Reals& source_values,
-                                   Vec* loadVec_out);
+PetscErrorCode calculateLoadVectorMI(Omega_h::Mesh& target_mesh,
+                                     Omega_h::Mesh& source_mesh,
+                                     const IntersectionResults& intersection,
+                                     const Omega_h::Reals& source_values,
+                                     Vec* loadVec_out);
+
+/**
+ * @brief Assembles the global load vector using Monte Carlo integration.
+ *
+ * This function computes the unassembled local load vector contributions for
+ * each element in the target mesh using Monte Carlo integration and
+ * then assembles them into a global PETSc vector in COO format.
+ *
+ * The element-local contributions are computed from source-field values already
+ * evaluated at the sampled physical points in each target element. The sampling
+ * pattern is defined on the reference triangle and may be generated either by
+ * uniform random sampling or by precomputed Sobol barycentric samples.
+ *
+ * @param target_mesh The target Omega_h mesh to which the scalar field is being
+ * projected.
+ * @param field_values_at_points Source-field values evaluated at the sampled
+ * physical points in the target mesh. The data is expected to be stored
+ * element-by-element, with total size
+ * `target_mesh.nelems() * npoints_each_tri`.
+ * @param npoints_each_tri Number of sample points used in each target triangle
+ * for the Monte Carlo integration.
+ * @param method Sampling method used to define the reference barycentric sample
+ * coordinates. Supported options include uniform random sampling and Sobol
+ * sampling.
+ * @param sobol_filename Path to the file containing precomputed Sobol
+ * barycentric samples when `method` is `SamplingMethod::SOBOL`.
+ * @param[out] loadVec_out Pointer to a PETSc Vec where the assembled load
+ * vector will be stored.
+ *
+ * @return PetscErrorCode Returns PETSC_SUCCESS if successful, or an appropriate
+ * PETSc error code otherwise.
+ *
+ * @note
+ * - Works for 2D linear triangular elements.
+ * - Uses COO-style preallocation and insertion into the PETSc vector.
+ * - Internally computes per-element Monte Carlo load vector contributions
+ *   before assembling the global vector.
+ * - For linear triangular elements, the barycentric coordinates are equal to
+ *   the local shape-function values used in the Monte Carlo estimator.
+ * - The resulting vector is used as the right-hand side (RHS) in a projection
+ *   solve.
+ *
+ * @see buildLoadVectorMC, SamplingMethod
+ */
+PetscErrorCode calculateLoadVectorMC(
+  Omega_h::Mesh& target_mesh, const Omega_h::Reals& field_values_at_points,
+  const int npoints_each_tri, SamplingMethod method, Vec* loadVec_out,
+  const std::string sobol_filename = "");
 } // namespace pcms
 #endif // PCMS_TRANSFER_CALCULATE_LOAD_VECTOR_HPP

src/pcms/transfer/calculate_mass_matrix.cpp

@@ -36,7 +36,7 @@ PetscErrorCode calculateMassMatrix(Omega_h::Mesh& mesh, Mat* mass_out)
 {
   PetscFunctionBeginUser;
 
-  MeshField::OmegahMeshField<DefaultExecutionSpace, 2,
+  MeshField::OmegahMeshField<Kokkos::DefaultExecutionSpace, 2,
                              MeshField::KokkosController>
     omf(mesh);
 

src/pcms/transfer/conservative_projection_solver.cpp

@@ -48,16 +48,22 @@ static Vec solveLinearSystem(Mat A, Vec b)
   ierr = KSPSolve(ksp, b, x);
   CHKERRABORT(PETSC_COMM_WORLD, ierr);
 
-  /// compute and print condition number estimate
-  PetscReal smax = 0.0, smin = 0.0;
-  ierr = KSPComputeExtremeSingularValues(ksp, &smax, &smin);
-  if (!ierr && smin > 0.0) {
-    PetscPrintf(PETSC_COMM_WORLD,
-                "Estimated condition number of matrix A: %.6e\n", smax / smin);
-  } else {
-    PetscPrintf(PETSC_COMM_WORLD,
-                "Condition number estimate unavailable (smin <= 0 or error)\n");
-  }
+  KSPConvergedReason reason;
+  PetscInt its = 0;
+  PetscReal rnorm = 0.0;
+
+  ierr = KSPGetConvergedReason(ksp, &reason);
+  CHKERRABORT(PETSC_COMM_WORLD, ierr);
+
+  ierr = KSPGetIterationNumber(ksp, &its);
+  CHKERRABORT(PETSC_COMM_WORLD, ierr);
+
+  ierr = KSPGetResidualNorm(ksp, &rnorm);
+  CHKERRABORT(PETSC_COMM_WORLD, ierr);
+
+  PetscPrintf(PETSC_COMM_WORLD,
+              "KSP reason=%d iterations=%d residual=%e\n",
+              (int)reason, (int)its, (double)rnorm);
 
   ierr = KSPDestroy(&ksp);
   CHKERRABORT(PETSC_COMM_WORLD, ierr);
@@ -86,10 +92,9 @@ static Omega_h::Reals vecToOmegaHReals(Vec vec)
   return Omega_h::Reals(values_host);
 }
 
-Omega_h::Reals solveGalerkinProjection(Omega_h::Mesh& target_mesh,
-                                       Omega_h::Mesh& source_mesh,
-                                       const IntersectionResults& intersection,
-                                       const Omega_h::Reals& source_values)
+Omega_h::Reals solveGalerkinProjectionMI(
+  Omega_h::Mesh& target_mesh, Omega_h::Mesh& source_mesh,
+  const IntersectionResults& intersection, const Omega_h::Reals& source_values)
 {
   if ((PetscInt)source_values.size() !=
       source_mesh.coords().size() / source_mesh.dim()) {
@@ -105,8 +110,8 @@ Omega_h::Reals solveGalerkinProjection(Omega_h::Mesh& target_mesh,
   CHKERRABORT(PETSC_COMM_WORLD, ierr);
 
   Vec vec;
-  ierr = calculateLoadVector(target_mesh, source_mesh, intersection,
-                             source_values, &vec);
+  ierr = calculateLoadVectorMI(target_mesh, source_mesh, intersection,
+                               source_values, &vec);
   CHKERRABORT(PETSC_COMM_WORLD, ierr);
 
   Vec x = solveLinearSystem(mass, vec);
@@ -123,15 +128,46 @@ Omega_h::Reals solveGalerkinProjection(Omega_h::Mesh& target_mesh,
 
   return solution_vector;
 }
-Omega_h::Reals rhsVectorMI(Omega_h::Mesh& target_mesh,
-                           Omega_h::Mesh& source_mesh,
-                           const IntersectionResults& intersection,
-                           const Omega_h::Reals& source_values)
+
+Omega_h::Reals solveGalerkinProjectionMC(
+  Omega_h::Mesh& target_mesh, const Omega_h::Reals& field_values_at_points,
+  const int npoints_each_tri, SamplingMethod method,
+  const std::string sobol_filename)
+{
+
+  Mat mass;
+  PetscErrorCode ierr = calculateMassMatrix(target_mesh, &mass);
+  CHKERRABORT(PETSC_COMM_WORLD, ierr);
+
+  Vec vec;
+  ierr = calculateLoadVectorMC(target_mesh, field_values_at_points,
+                               npoints_each_tri, method, &vec, sobol_filename);
+  CHKERRABORT(PETSC_COMM_WORLD, ierr);
+
+  Vec x = solveLinearSystem(mass, vec);
+  auto solution_vector = vecToOmegaHReals(x);
+
+  ierr = VecDestroy(&x);
+  CHKERRABORT(PETSC_COMM_WORLD, ierr);
+
+  ierr = MatDestroy(&mass);
+  CHKERRABORT(PETSC_COMM_WORLD, ierr);
+
+  ierr = VecDestroy(&vec);
+  CHKERRABORT(PETSC_COMM_WORLD, ierr);
+
+  return solution_vector;
+}
+
+Omega_h::Reals computeRhsVectorMI(Omega_h::Mesh& target_mesh,
+                                  Omega_h::Mesh& source_mesh,
+                                  const IntersectionResults& intersection,
+                                  const Omega_h::Reals& source_values)
 {
   Vec vec;
   PetscErrorCode ierr;
-  ierr = calculateLoadVector(target_mesh, source_mesh, intersection,
-                             source_values, &vec);
+  ierr = calculateLoadVectorMI(target_mesh, source_mesh, intersection,
+                               source_values, &vec);
   CHKERRABORT(PETSC_COMM_WORLD, ierr);
 
   auto rhsvector = vecToOmegaHReals(vec);
@@ -141,4 +177,26 @@ Omega_h::Reals rhsVectorMI(Omega_h::Mesh& target_mesh,
 
   return rhsvector;
 }
+
+Omega_h::Reals computeRhsVectorMC(Omega_h::Mesh& target_mesh,
+                                  const Omega_h::Reals& field_values_at_points,
+                                  const int npoints_each_tri,
+                                  SamplingMethod method,
+                                  const std::string& sobol_filename)
+{
+
+    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;
+}	
+
 } // namespace pcms