From b53cc15b1ce00df4c48e027d357883fab2cc6275 Mon Sep 17 00:00:00 2001 From: Jacob Merson Date: Sat, 11 Jul 2026 11:11:56 -0400 Subject: [PATCH 01/15] transfer: remove old-style conservative projection (calculate_mass_matrix / load_vector) Squashed original commits: 0a27b53, 0c0e697 --- src/pcms/transfer/calculate_load_vector.cpp | 52 --- src/pcms/transfer/calculate_load_vector.hpp | 59 --- src/pcms/transfer/calculate_mass_matrix.cpp | 74 ---- src/pcms/transfer/calculate_mass_matrix.hpp | 34 -- src/pcms/transfer/load_vector_integrator.cpp | 430 ------------------- src/pcms/transfer/load_vector_integrator.hpp | 220 ---------- test/test_load_vector.cpp | 143 ------ 7 files changed, 1012 deletions(-) delete mode 100644 src/pcms/transfer/calculate_load_vector.cpp delete mode 100644 src/pcms/transfer/calculate_load_vector.hpp delete mode 100644 src/pcms/transfer/calculate_mass_matrix.cpp delete mode 100644 src/pcms/transfer/calculate_mass_matrix.hpp delete mode 100644 src/pcms/transfer/load_vector_integrator.cpp delete mode 100644 src/pcms/transfer/load_vector_integrator.hpp delete mode 100644 test/test_load_vector.cpp diff --git a/src/pcms/transfer/calculate_load_vector.cpp b/src/pcms/transfer/calculate_load_vector.cpp deleted file mode 100644 index f351a4ef..00000000 --- a/src/pcms/transfer/calculate_load_vector.cpp +++ /dev/null @@ -1,52 +0,0 @@ -#include "pcms/transfer/calculate_load_vector.hpp" -#include "pcms/transfer/coo_assembly_utils.hpp" -#include "pcms/transfer/load_vector_integrator.hpp" -#include "pcms/transfer/petsc_utils.hpp" - -namespace pcms -{ -PetscErrorCode calculateLoadVector(Omega_h::Mesh& target_mesh, - Omega_h::Mesh& source_mesh, - const IntersectionResults& intersection, - const Omega_h::Reals& source_values, - Vec* loadVec_out) -{ - - 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 = - buildLoadVector(target_mesh, source_mesh, intersection, source_values); - - auto hostElmLoadVector = Kokkos::create_mirror_view(elmLoadVector); - Kokkos::deep_copy(hostElmLoadVector, elmLoadVector); - PetscCheck(static_cast(hostElmLoadVector.extent(0)) == nnz, - PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, - "Element load vector size (%d) does not match COO nnz (%d)", - static_cast(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); -} -} // namespace pcms diff --git a/src/pcms/transfer/calculate_load_vector.hpp b/src/pcms/transfer/calculate_load_vector.hpp deleted file mode 100644 index b1968338..00000000 --- a/src/pcms/transfer/calculate_load_vector.hpp +++ /dev/null @@ -1,59 +0,0 @@ -/** - * @file calculate_load_vector.hpp - * @brief Routines for assembling global load vector in conservative field - * projection. - * - * Provides functionality to compute and assemble the global load vector - * used in Galerkin-based conservative field transfer between non-conforming - * meshes. - * - */ - -#ifndef PCMS_TRANSFER_CALCULATE_LOAD_VECTOR_HPP -#define PCMS_TRANSFER_CALCULATE_LOAD_VECTOR_HPP -#include -#include -#include - -/** - * @brief Assembles the global load vector. - * - * 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. - * - * - * @param target_mesh The target Omega_h mesh to which the scalar field is being - * projected. - * @param source_mesh The source Omega_h mesh containing the original scalar - * field values. - * @param intersection Precomputed intersection data for each target element. - * Includes the number and indices of intersecting source - * elements. - * @param source_values Nodal scalar field values defined on the source mesh. - * @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 calls `buildLoadVector()` to compute per-element contributions. - * - The resulting vector is used as the right-hand side (RHS) in a projection - * solve. - * - * @see buildLoadVector,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); -} // namespace pcms -#endif // PCMS_TRANSFER_CALCULATE_LOAD_VECTOR_HPP diff --git a/src/pcms/transfer/calculate_mass_matrix.cpp b/src/pcms/transfer/calculate_mass_matrix.cpp deleted file mode 100644 index 2b54e4c6..00000000 --- a/src/pcms/transfer/calculate_mass_matrix.cpp +++ /dev/null @@ -1,74 +0,0 @@ -#include "pcms/transfer/calculate_mass_matrix.hpp" -#include "pcms/transfer/coo_assembly_utils.hpp" -#include "pcms/transfer/mass_matrix_integrator.hpp" -#include "pcms/transfer/petsc_utils.hpp" -#include "pcms/utility/memory_spaces.h" - -namespace pcms -{ -/** - * @brief Creates a PETSc matrix based on mesh connectivity - * - * This function creates a sparse matrix with the proper sparsity pattern - * according to the mesh connectivity. The matrix size corresponds to the - * number of vertices in the mesh. - * - * @param mesh The Omega_h mesh to create the matrix from - * @param[out] A Pointer to the PETSc matrix to be created - * @return PetscErrorCode PETSc error code (PETSC_SUCCESS if successful) - */ -static PetscErrorCode create_linear_triangle_coo_matrix(Omega_h::Mesh& mesh, - Mat* A) -{ - PetscInt* coo_rows = nullptr; - PetscInt* coo_cols = nullptr; - PetscInt matSize = 0; - PetscFunctionBeginUser; - PetscCall(createSeqAIJMat(PETSC_COMM_WORLD, mesh.nverts(), mesh.nverts(), 0, - nullptr, A)); - PetscCall( - build_linear_triangle_coo_rows_cols(mesh, &coo_rows, &coo_cols, &matSize)); - PetscCall(MatSetPreallocationCOO(*A, matSize, coo_rows, coo_cols)); - PetscCall(PetscFree2(coo_rows, coo_cols)); - PetscFunctionReturn(PETSC_SUCCESS); -} -PetscErrorCode calculateMassMatrix(Omega_h::Mesh& mesh, Mat* mass_out) -{ - PetscFunctionBeginUser; - - MeshField::OmegahMeshField - omf(mesh); - - const auto ShapeOrder = 1; - auto coordField = omf.getCoordField(); - const auto [shp, map] = - MeshField::Omegah::getTriangleElement(mesh); - MeshField::FieldElement coordFe(mesh.nelems(), coordField, shp, map); - - auto elmMassMatrix = buildMassMatrix(mesh, coordFe); - const PetscInt expected_nnz = 9 * mesh.nelems(); - PetscCheck(static_cast(elmMassMatrix.extent(0)) == expected_nnz, - PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, - "Element mass data size (%d) does not match expected linear COO " - "entries (%d)", - static_cast(elmMassMatrix.extent(0)), expected_nnz); - - Mat mass; - PetscCall(create_linear_triangle_coo_matrix(mesh, &mass)); - PetscCall(MatZeroEntries(mass)); - PetscCall( - MatSetValuesCOO(mass, elmMassMatrix.data(), - INSERT_VALUES)); // FIXME fails here on gpu, calls into host - // implementation... AFAIK, petsc checks - // the type of the input array of values to - // decide which backend to use... - // - if (mesh.nelems() < 10) { - PetscCall(MatView(mass, PETSC_VIEWER_STDOUT_WORLD)); - } - - *mass_out = mass; - PetscFunctionReturn(PETSC_SUCCESS); -} -} // namespace pcms diff --git a/src/pcms/transfer/calculate_mass_matrix.hpp b/src/pcms/transfer/calculate_mass_matrix.hpp deleted file mode 100644 index 0889fdc8..00000000 --- a/src/pcms/transfer/calculate_mass_matrix.hpp +++ /dev/null @@ -1,34 +0,0 @@ -/** - * @file calculateMassMatrix.hpp - * @brief Functions for calculating mass matrices on finite element meshes - * @author [Cameron Smith] - * @date April, 2025 - * - * This file contains functions for creating and computing mass matrices - * for finite element calculations using Omega_h mesh structures and PETSc. - */ - -#ifndef PCMS_TRANSFER_CALCULATE_MASS_MATRIX_HPP -#define PCMS_TRANSFER_CALCULATE_MASS_MATRIX_HPP - -#include -#include - -namespace pcms -{ -/** - * @brief Calculates the mass matrix for a given mesh - * - * This function constructs a mass matrix based on the provided mesh using - * a finite element approach. It creates coordinate field elements, builds - * the mass matrix using the massMatrixIntegrator, and sets up the PETSc matrix - * with appropriate values. - * - * @param mesh The Omega_h mesh to calculate the mass matrix for - * @param[out] mass_out Pointer to the resulting mass matrix - * @return PetscErrorCode PETSc error code (PETSC_SUCCESS if successful) - */ - -PetscErrorCode calculateMassMatrix(Omega_h::Mesh& mesh, Mat* mass_out); -} // namespace pcms -#endif // PCMS_TRANSFER_CALCULATE_MASS_MATRIX_HPP diff --git a/src/pcms/transfer/load_vector_integrator.cpp b/src/pcms/transfer/load_vector_integrator.cpp deleted file mode 100644 index 221e3905..00000000 --- a/src/pcms/transfer/load_vector_integrator.cpp +++ /dev/null @@ -1,430 +0,0 @@ -#include "pcms/transfer/load_vector_integrator.hpp" - -namespace pcms -{ - -/** - * @brief Converts barycentric coordinates to global (physical) coordinates. - * - * Given barycentric coordinates within a 2D triangle and the coordinates of - * the triangle's vertices, this function computes the corresponding global - * position. - * - * @param barycentric_coord The barycentric coordinates \f$(\lambda_1, - * \lambda_2, \lambda_3)\f$ of the point. - * @param verts_coord The coordinates of the triangle's three vertices in global - * space. - * @return The 2D global coordinates corresponding to the given barycentric - * position. - */ - -[[nodiscard]] OMEGA_H_INLINE Omega_h::Vector<2> global_from_barycentric( - const MeshField::Vector3& barycentric_coord, - const Omega_h::Few, 3>& verts_coord) -{ - Omega_h::Vector<2> real_coords = {0.0, 0.0}; - - for (int i = 0; i < 3; ++i) { - real_coords[0] += barycentric_coord[i] * verts_coord[i][0]; - real_coords[1] += barycentric_coord[i] * verts_coord[i][1]; - } - return real_coords; -} - -/** - * @brief Computes the barycentric coordinates of a 2D point with respect to a - * triangle. - * - * Given a point in global (x, y) coordinates and the coordinates of the three - * vertices of a triangle, this function evaluates the barycentric coordinates - * \f$(\lambda_1, \lambda_2, \lambda_3)\f$ of the point with respect to that - * triangle. - * - * @param point The 2D global coordinates of the point to evaluate (in - * Omega_h::Vector<2> format). - * @param verts_coord The vertex coordinates of the triangle (in r3d::Vector<2> - * format). - * @return A vector of three barycentric coordinates corresponding to the input - * point. - * - */ - -[[nodiscard]] OMEGA_H_INLINE Omega_h::Vector<3> evaluate_barycentric( - const Omega_h::Vector<2>& point, - const r3d::Few, 3>& verts_coord) -{ - Omega_h::Few, 3> omegah_vector; - for (int i = 0; i < 3; ++i) { - omegah_vector[i][0] = verts_coord[i][0]; - omegah_vector[i][1] = verts_coord[i][1]; - } - - auto barycentric_coordinate = - Omega_h::barycentric_from_global<2, 2>(point, omegah_vector); - - return barycentric_coordinate; -} -/** - * @brief Evaluates the value of a linear function at a given point using - * barycentric coordinates. - * - * This function computes the interpolated value of a nodal scalar field over a - * triangle, using barycentric coordinates within the specified element. - * - * @param nodal_values The global array of nodal field values. - * @param faces2nodes The element-to-node connectivity array. - * @param bary_coords The barycentric coordinates of the evaluation point within - * the triangle. - * @param elm_id The ID of the triangle element being evaluated. - * @return The interpolated function value at the given point. - * - * @note This function assumes linear (3-node) triangular element. - */ - -[[nodiscard]] OMEGA_H_INLINE double evaluate_function_value( - const Omega_h::Reals& nodal_values, const Omega_h::LOs& faces2nodes, - const Omega_h::Vector<3>& bary_coords, const int elm_id) -{ - - double value = 0; - const auto elm_verts = Omega_h::gather_verts<3>(faces2nodes, elm_id); - for (int i = 0; i < 3; ++i) { - int nid = elm_verts[i]; - value += nodal_values[nid] * bary_coords[i]; - } - - return value; -} - -/** - * @brief Deduplicate, reorder, and orient a 2D polygon produced by r3d. - * - * This function performs all necessary cleanup and reconstruction of a polygon - * returned by `r3d::intersect_simplices()`, which may contain: - * - duplicated vertices, - * - unordered vertex lists, - * - invalid or inconsistent neighbor links (`pnbrs`), - * - negative orientation (CW instead of CCW). - * - * The cleanup proceeds with the following stages: - * - * **(1) Geometric deduplication:** - * Vertices whose coordinates are equal within a tolerance `tol` are collapsed - * into a single unique vertex. A compacted vertex list is built. - * - * **(2) CCW vertex reordering:** - * The remaining unique vertices are sorted by their polar angle around the - * polygon centroid. This yields a globally consistent counter-clockwise (CCW) - * - * - * **(3) Rebuilding neighbor links:** - * After sorting, each vertex's two neighbors (`pnbrs[0]` and `pnbrs[1]`) are - * reassigned to form a closed CCW cycle: - * - * pnbrs[0] = previous vertex in CCW order - * pnbrs[1] = next vertex in CCW order - * - * - * - * @param poly The polygon to clean, and reorder. - * - * @param tol Tolerance for geometric duplicate detection (default: 1e-12). - * - * @return The number of unique, CCW-ordered vertices remaining in the polygon. - * - * @see r3d::Polytope - * @see r3d::measure - * @see r3d::intersect_simplices - */ - -[[nodiscard]] OMEGA_H_INLINE int remove_duplicate_vertices_and_fix_links( - r3d::Polytope<2>& poly, const double tol = 1e-12) -{ - - const int old_n = poly.nverts; - int new_n = 0; - - // geometric deupliactes filtering - for (int i = 0; i < old_n; ++i) { - const auto& pi = poly.verts[i].pos; - bool dup = false; - for (int j = 0; j < new_n; ++j) { - const auto& pj = poly.verts[j].pos; - if (Kokkos::fabs(pi[0] - pj[0]) < tol && - Kokkos::fabs(pi[1] - pj[1]) < tol) { - dup = true; - break; - } - } - if (!dup) { - poly.verts[new_n] = poly.verts[i]; - ++new_n; - } - } - poly.nverts = new_n; - - // CCW reorder verts by angle about centroid - if (new_n >= 3) { - // centroid - double cx = 0.0; - double cy = 0.0; - for (int i = 0; i < new_n; ++i) { - cx += poly.verts[i].pos[0]; - cy += poly.verts[i].pos[1]; - } - cx /= new_n; - cy /= new_n; - - // angle sort - int order[r3d::MaxVerts<2>::value]; - for (int i = 0; i < new_n; ++i) { - order[i] = i; - } - - for (int i = 0; i < new_n - 1; ++i) { - for (int j = i + 1; j < new_n; ++j) { - const double a1 = Kokkos::atan2(poly.verts[order[i]].pos[1] - cy, - poly.verts[order[i]].pos[0] - cx); - const double a2 = Kokkos::atan2(poly.verts[order[j]].pos[1] - cy, - poly.verts[order[j]].pos[0] - cx); - if (a1 > a2) { - int t = order[i]; - order[i] = order[j]; - order[j] = t; - } - } - } - - r3d::Vertex<2> tmp[r3d::MaxVerts<2>::value]; - for (int i = 0; i < new_n; ++i) { - tmp[i] = poly.verts[order[i]]; - } - - for (int i = 0; i < new_n; ++i) { - poly.verts[i] = tmp[i]; - } - - // set circular neighbors consistent with verts[] order - for (int i = 0; i < new_n; ++i) { - const int prev = (i - 1 + new_n) % new_n; - const int next = (i + 1) % new_n; - poly.verts[i].pnbrs[0] = prev; // CCW prev - poly.verts[i].pnbrs[1] = next; // CCW next - } - - // ensure positive orientation (if measure is still negative, reverse) - double area = r3d::measure(poly); - if (area < 0.0) { - // reverse verts and relink - for (int i = 0; i < new_n / 2; ++i) { - r3d::Vertex<2> t = poly.verts[i]; - poly.verts[i] = poly.verts[new_n - 1 - i]; - poly.verts[new_n - 1 - i] = t; - } - for (int i = 0; i < new_n; ++i) { - const int prev = (i - 1 + new_n) % new_n; - const int next = (i + 1) % new_n; - poly.verts[i].pnbrs[0] = prev; - poly.verts[i].pnbrs[1] = next; - } - } - } else { - // clear links - for (int i = 0; i < new_n; ++i) { - poly.verts[i].pnbrs[0] = -1; - poly.verts[i].pnbrs[1] = -1; - } - } - - return new_n; -} - -template -OMEGA_H_INLINE void for_each_intersection_subtriangle( - const int elm, const IntersectionResults& intersection, - const Omega_h::Reals& tgt_coords, const Omega_h::Reals& src_coords, - const Omega_h::LOs& tgt_faces2nodes, const Omega_h::LOs& src_faces2nodes, - TriangleOp&& op) -{ - auto tgt_elm_vert_coords = - get_vert_coords_of_elem(tgt_coords, tgt_faces2nodes, elm); - const int start = intersection.tgt2src_offsets[elm]; - const int end = intersection.tgt2src_offsets[elm + 1]; - - for (int i = start; i < end; ++i) { - const int current_src_elm = intersection.tgt2src_indices[i]; - auto src_elm_vert_coords = - get_vert_coords_of_elem(src_coords, src_faces2nodes, current_src_elm); - r3d::Polytope<2> poly; - r3d::intersect_simplices(poly, tgt_elm_vert_coords, src_elm_vert_coords); - auto nverts = remove_duplicate_vertices_and_fix_links(poly, 1e-12); - ; - auto poly_area = r3d::measure(poly); - - for (int j = 1; j < nverts - 1; ++j) { - // build triangle from poly.verts[0], poly.verts[j], - // poly.verts[j+1] - auto& p0 = poly.verts[0].pos; - auto& p1 = poly.verts[j].pos; - auto& p2 = poly.verts[j + 1].pos; - - Omega_h::Few, 3> tri_coords; - tri_coords[0] = {p0[0], p0[1]}; - tri_coords[1] = {p1[0], p1[1]}; - tri_coords[2] = {p2[0], p2[1]}; - - Omega_h::Few, 2> basis; - basis[0] = tri_coords[1] - tri_coords[0]; - basis[1] = tri_coords[2] - tri_coords[0]; - - Omega_h::Real area = - Kokkos::fabs(Omega_h::triangle_area_from_basis(basis)); - - const double EPS_AREA = abs_tol + rel_tol * poly_area; - if (area <= EPS_AREA) - continue; // drops duplicates and colinear/degenerates - - op(tri_coords, tgt_elm_vert_coords, src_elm_vert_coords, current_src_elm, - area); - } - } -} - -Kokkos::View buildLoadVector( - Omega_h::Mesh& target_mesh, Omega_h::Mesh& source_mesh, - const IntersectionResults& intersection, const Omega_h::Reals& source_values) -{ - - const auto& tgt_coords = target_mesh.coords(); - const auto& src_coords = source_mesh.coords(); - const auto& tgt_faces2nodes = - target_mesh.ask_down(Omega_h::FACE, Omega_h::VERT).ab2b; - const auto& src_faces2nodes = - source_mesh.ask_down(Omega_h::FACE, Omega_h::VERT).ab2b; - - IntegrationData<2> integrationPoints; - int npts = integrationPoints.size(); - - // TODO: Make it generalised; hardcoded for liner 2D - Kokkos::View elmLoadVector( - "elmLoadVector", static_cast(target_mesh.nelems()) * 3); - Kokkos::parallel_for( - "calculate load vector", target_mesh.nelems(), - KOKKOS_LAMBDA(const int& elm) { - Omega_h::Vector<3> part_integration = {0.0, 0.0, 0.0}; - for_each_intersection_subtriangle( - elm, intersection, tgt_coords, src_coords, tgt_faces2nodes, - src_faces2nodes, - [&](const Omega_h::Few, 3>& tri_coords, - const r3d::Few, 3>& tgt_elm_vert_coords, - const r3d::Few, 3>& src_elm_vert_coords, - const int current_src_elm, const Omega_h::Real area) { - for (int ip = 0; ip < npts; ++ip) { - auto bary = integrationPoints.bary_coords(ip); - auto weight = integrationPoints.weights(ip); - - // convert barycentric to real coords in triangle - auto real_coords = global_from_barycentric(bary, tri_coords); - - // evaluate shape function (barycentric wrt target for linear) - auto shape_fn = - evaluate_barycentric(real_coords, tgt_elm_vert_coords); - - // evaluate function at point (barycentric wrt source for linear) - auto src_bary = - evaluate_barycentric(real_coords, src_elm_vert_coords); - auto fval = evaluate_function_value(source_values, src_faces2nodes, - src_bary, current_src_elm); - - // integration - for (int k = 0; k < 3; ++k) { - part_integration[k] += shape_fn[k] * fval * weight * 2 * area; - } - } - }); - - for (int j = 0; j < 3; ++j) { - elmLoadVector(elm * 3 + j) = part_integration[j]; - } - }); - - return elmLoadVector; -} -Errors evaluate_proj_and_cons_errors(Omega_h::Mesh& target_mesh, - Omega_h::Mesh& source_mesh, - const IntersectionResults& intersection, - const Omega_h::Reals& target_values, - const Omega_h::Reals& source_values) -{ - - const auto& tgt_coords = target_mesh.coords(); - const auto& src_coords = source_mesh.coords(); - const auto& tgt_faces2nodes = - target_mesh.ask_down(Omega_h::FACE, Omega_h::VERT).ab2b; - const auto& src_faces2nodes = - source_mesh.ask_down(Omega_h::FACE, Omega_h::VERT).ab2b; - - IntegrationData<2> integrationPoints; - int npts = integrationPoints.size(); - - constexpr double EPS_DEN = 1e-30; - - Kokkos::View accum("accum", 4); - Kokkos::deep_copy(accum, 0.0); - - Kokkos::parallel_for( - "evaluate relative errors", target_mesh.nelems(), - KOKKOS_LAMBDA(const int& elm) { - double N2 = 0.0, D2 = 0.0, C = 0.0, QD = 0.0; - for_each_intersection_subtriangle( - elm, intersection, tgt_coords, src_coords, tgt_faces2nodes, - src_faces2nodes, - [&](const Omega_h::Few, 3>& tri_coords, - const r3d::Few, 3>& tgt_elm_vert_coords, - const r3d::Few, 3>& src_elm_vert_coords, - const int current_src_elm, const Omega_h::Real area) { - for (int ip = 0; ip < npts; ++ip) { - auto bary = integrationPoints.bary_coords(ip); - auto weight = integrationPoints.weights(ip); - - // convert barycentric to real coords in triangle - auto real_coords = global_from_barycentric(bary, tri_coords); - auto tgt_bary = - evaluate_barycentric(real_coords, tgt_elm_vert_coords); - - // evaluate shape function (barycentric wrt target for linear) - auto tgtVal = evaluate_function_value( - target_values, tgt_faces2nodes, tgt_bary, elm); - - // evaluate function at point (barycentric wrt source for linear) - auto src_bary = - evaluate_barycentric(real_coords, src_elm_vert_coords); - auto srcVal = evaluate_function_value( - source_values, src_faces2nodes, src_bary, current_src_elm); - - // integration - auto diff = srcVal - tgtVal; - auto w = 2 * weight * area; - N2 += diff * diff * w; - D2 += srcVal * srcVal * w; - C += diff * w; - QD += srcVal * w; - } - }); - - Kokkos::atomic_add(&accum(0), N2); - Kokkos::atomic_add(&accum(1), D2); - Kokkos::atomic_add(&accum(2), C); - Kokkos::atomic_add(&accum(3), QD); - }); - - auto h_accum = Kokkos::create_mirror(accum); - Kokkos::deep_copy(h_accum, accum); - const double proj_err = - Kokkos::sqrt(h_accum(0)) / Kokkos::max(Kokkos::sqrt(h_accum(1)), EPS_DEN); - const double cons_err = - Kokkos::fabs(h_accum(2)) / Kokkos::max(Kokkos::fabs(h_accum(3)), EPS_DEN); - - return Errors{.proj_err = proj_err, .cons_err = cons_err}; -} -} // namespace pcms diff --git a/src/pcms/transfer/load_vector_integrator.hpp b/src/pcms/transfer/load_vector_integrator.hpp deleted file mode 100644 index 4e1d2fad..00000000 --- a/src/pcms/transfer/load_vector_integrator.hpp +++ /dev/null @@ -1,220 +0,0 @@ -/** - * @file load_vector_integrator.hpp - * @brief Functions for computing load vectors in conservative field projection. - * - * This file implements routines to compute the element-wise load vector - * (right-hand side) contributions used in Galerkin projection of scalar - * fields from a source mesh to a target mesh. - * - * The integration is performed over polygonal intersections (supermesh) between - * source and target elements using barycentric quadrature. The resulting values - * represent unassembled local contributions that can later be combined into a - * global load vector. - * - * @note - * - Assumes 2D linear triangular meshes. - * - Intersection data is provided via the `IntersectionResults` structure. - */ -#ifndef PCMS_TRANSFER_LOAD_VECTOR_INTEGRATOR_HPP -#define PCMS_TRANSFER_LOAD_VECTOR_INTEGRATOR_HPP - -#include -#include -#include -#include -#include -#include - -namespace pcms -{ -/** - * @brief Computes the load vector for each target element in the conservative - * field transfer. - * - * This routine is used for constructing the right-hand side (RHS) of the - * conservative field transfer formulation, projecting field quantities from the - * source mesh to the target mesh. - * - * The underlying algorithm computes contributions to the load vector - * using geometric intersection data between source and target elements. - * - * @note Currently this method works for a two-dimensional linear triangles. - */ - -/** - * @brief Provides barycentric integration points and weights for a triangle - * element. - * - * This templated struct stores the barycentric coordinates - * and quadrature weights for performing numerical integration over a reference - * triangle. It is used for integrating functions over elements in the - * conservative field transfer. - * - * @tparam order The quadrature order (number of integration points and - * polynomial accuracy). - */ -template -struct IntegrationData -{ - // Barycentric coordinates of integration points - Kokkos::View bary_coords; - - // Quadrature weights associated with each integration point - Kokkos::View weights; - - /** - * @brief Constructs the integration data for a given quadrature order - * - * Initializes barycentric coordinates and weights using - * MeshField's predefined triangle quadrature rules. - */ - IntegrationData() - { - auto ip_vec = MeshField::getIntegrationPoints(order); - std::size_t num_ip = ip_vec.size(); - - bary_coords = Kokkos::View("bary_coords", num_ip); - weights = Kokkos::View("weights", num_ip); - - auto bary_coords_host = Kokkos::create_mirror_view(bary_coords); - auto weights_host = Kokkos::create_mirror_view(weights); - - for (std::size_t i = 0; i < num_ip; ++i) { - bary_coords_host(i) = ip_vec[i].param; - weights_host(i) = ip_vec[i].weight; - } - - Kokkos::deep_copy(bary_coords, bary_coords_host); - Kokkos::deep_copy(weights, weights_host); - } - - /** - * @brief Returns the number of integration points - * - * @return Number of integration points for the selected order. - */ - - int size() const { return bary_coords.extent(0); } -}; - -/** - * @brief Computes the per-element RHS load vectors for conservative field - * projection from source to target mesh. - * - * This function computes local (element-wise) right-hand side (RHS) - * contributions for the Galerkin projection of a scalar field from the source - * mesh to the target mesh. It integrates over the polygonal intersection - * regions between each target element and its intersecting source elements - * using barycentric quadrature. - * - * The output is a flat array containing unassembled load vector contributions - * at the nodes of each target triangle. - * - * @param target_mesh The target mesh object receiving the projected scalar - * field. - * @param source_mesh The source mesh object containing the original scalar - * field values. - * @param intersection Precomputed intersection data for each target element. - * Includes the number and indices of intersecting source - * elements. - * @param source_values Scalar field values defined at the nodes of the source - * mesh. - * - * @return A Kokkos view containing per-element load vectors. - * Each triangle contributes 3 values (one per node), so the view has - * size 3 × (number of target elements). - * - * @note - * - This function assumes 2D linear triangular elements. - * - Degenerate or near-zero-area intersection polygons are skipped. - * - Each polygon is triangulated using a fan structure and integrated using - * barycentric quadrature rules. - * - The returned vector must be assembled into a global RHS vector in a later - * step. - * - * @see evaluate_barycentric, evaluate_function_value, global_from_barycentric - * @see IntersectionResults - */ - -Kokkos::View buildLoadVector( - Omega_h::Mesh& target_mesh, Omega_h::Mesh& source_mesh, - const IntersectionResults& intersection, const Omega_h::Reals& source_values); -/// Holds projection and conservation error metrics returned by -/// evaluate_pro_and_cons_errors(). -struct Errors -{ - double proj_err; ///< L2 projection error computed on the supermesh. - double cons_err; ///< Relative conservation error over the supermesh. -}; - -/** - * @brief Computes projection and conservation errors over the supermesh for - * scalar field transfer. - * - * This function quantifies the accuracy and conservation properties of - * conservative field transfer between nonconforming meshes using - * supermesh-based integration. It returns a struct containing two error - * metrics: - * - * - **Projection Error (`proj_err`)** — Measures the L2 norm of the difference - * between the projected source field and the target field over the - * supermesh.This reflects how accurately the field has been projected. - * - * - **Conservation Error (`cons_err`)** — Relative difference in the integrated - * field values between source and target representations. This captures - * conservation loss across the transfer. - * - * ### Mathematical Definitions: - * \f[ - * \text{proj\_err} = \frac{ \| q_D - q_T \|_{L_2(\Omega_S)} } - * { \| q_D \|_{L_2(\Omega_S)} }, \quad - * \text{cons\_err} = \frac{ \left| \int_{\Omega_S} q_D - \int_{\Omega_S} q_T - * \right| } { \left| \int_{\Omega_S} q_D \right| } - * \f] - * - * where: - * - \f$q_D\f$ is the scalar field defined on the source mesh (mesh from where - * the field is defined), - * - \f$q_T\f$ is the projected field on the target mesh (mesh to where the - * field is projected), - * - \f$\Omega_S\f$ is the supermesh formed by polygonal intersections of source - * and target elements. - * - * Integration is performed by triangulating each intersection region and - * applying barycentric quadrature. Degenerate or near-zero-area triangles are - * skipped based on area tolerance. - * - * - * @param target_mesh The target mesh object receiving the projected scalar - * field. - * @param source_mesh The source mesh object containing the original scalar - * field values. - * @param intersection Precomputed intersection data for each target element. - * Includes the number and indices of intersecting source - * elements. - * @param target_values Nodal scalar field values evaluated on the target mesh - * using galerkin projection. - * @param source_values Scalar field values defined at the nodes of the source - * mesh. - * - * - * @return A struct containing: - * - `proj_err`: Exact L2 projection error over the supermesh. - * - `cons_err`: Relative conservation error over the supermesh. - * - * @note - * - Assumes 2D linear (P1) triangular elements. - * - Ideal for validating conservative transfer schemes or testing projection - * fidelity. - * - * @see IntersectionResults, buildLoadVector - */ - -Errors evaluate_proj_and_cons_errors(Omega_h::Mesh& target_mesh, - Omega_h::Mesh& source_mesh, - const IntersectionResults& intersection, - const Omega_h::Reals& target_values, - const Omega_h::Reals& source_values); -} // namespace pcms - -#endif // PCMS_TRANSFER_LOAD_VECTOR_INTEGRATOR_HPP diff --git a/test/test_load_vector.cpp b/test/test_load_vector.cpp deleted file mode 100644 index 1a38df50..00000000 --- a/test/test_load_vector.cpp +++ /dev/null @@ -1,143 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -TEST_CASE("Load vector computation on intersected regions", "[load_vector]") -{ - - Omega_h::Library lib; - - // 2D coordinates (x,y) : 4 vertices - Omega_h::Reals coords({ - 0.0, 0.0, // v0 - 1.0, 0.0, // v1 - 1.0, 1.0, // v2 - 0.0, 1.0 // v3 - }); - - // Target Mesh with two triangles - // Two triangles, CCW - // T0: (v0,v1,v3) = (0,1,3) - // T1: (v1,v2,v3) = (1,2,3) - Omega_h::LOs ev2v_target({0, 1, 3, 1, 2, 3}); - - Omega_h::Mesh target_mesh(&lib); - Omega_h::build_from_elems_and_coords(&target_mesh, OMEGA_H_SIMPLEX, 2, - ev2v_target, coords); - - // Source Mesh with two triangles - // Two triangles, CCW - // T0: (v0,v1,v3) = (0,1,2) - // T1: (v1,v2,v3) = (0,2,3) - Omega_h::LOs ev2v_source({0, 1, 2, 0, 2, 3}); - - Omega_h::Mesh source_mesh(&lib); - Omega_h::build_from_elems_and_coords(&source_mesh, OMEGA_H_SIMPLEX, 2, - ev2v_source, coords); - - REQUIRE(source_mesh.dim() == 2); - REQUIRE(target_mesh.dim() == 2); - - int num_tgt_elems = target_mesh.nelems(); - - auto intersection = pcms::intersectTargets(source_mesh, target_mesh); - SECTION("check localization routine for coincident cases") - { - Kokkos::View points("test_points", 3); - auto points_h = Kokkos::create_mirror_view(points); - points_h(0, 0) = 0.0; - points_h(0, 1) = 0.0; - points_h(1, 0) = 1.0; - points_h(1, 1) = 0.0; - points_h(2, 0) = 1.0; - points_h(2, 1) = 1.0; - - Kokkos::deep_copy(points, points_h); - pcms::GridPointSearch2D search_cell(target_mesh, 20, 20); - - auto d_results = search_cell(points); - - auto h_results = Kokkos::create_mirror_view(d_results); - - Kokkos::deep_copy(h_results, d_results); - - REQUIRE(h_results.extent(0) == 3); - - for (int i = 0; i < h_results.extent(0); ++i) { - auto r = h_results(i); - std::cout << "tri_id = " << r.element_id << "\n"; - std::cout << "dim = " << static_cast(r.dimensionality) << "\n"; - // std::cout << "parametric = " << r.parametric_coords << "\n"; - INFO("tri_id : " << r.element_id); - REQUIRE(r.element_id >= 0); - } - } - - SECTION("Basic shape and consistency of result") - { - - // Fill dummy values at each vertex of source - Omega_h::Write values(source_mesh.nverts(), 1.0); - - auto load_vector = - pcms::buildLoadVector(target_mesh, source_mesh, intersection, values); - - auto load_vector_host = Kokkos::create_mirror(load_vector); - Kokkos::deep_copy(load_vector_host, load_vector); - - REQUIRE(static_cast(load_vector_host.extent(0)) == num_tgt_elems * 3); - - for (int i = 0; i < load_vector_host.extent(0); ++i) - REQUIRE(load_vector_host(i) >= - 0.0); // Since function is 1.0 and everything is positive - } - - SECTION("Integration is zero if source values are zero") - { - Omega_h::Write zero_field(source_mesh.nverts(), 0.0); - - auto load_vector = - pcms::buildLoadVector(target_mesh, source_mesh, intersection, zero_field); - - auto load_vector_host = Kokkos::create_mirror_view(load_vector); - Kokkos::deep_copy(load_vector_host, load_vector); - - for (int i = 0; i < load_vector_host.extent(0); ++i) { - REQUIRE(load_vector_host(i) == Catch::Approx(0.0)); - } - } - - SECTION("load vector computed after the intersection of simple target and " - "source elements") - { - // the source elements are triangle1 (0,0), (1,0) & (1,1) and triangle2 - // (0,0), (1,1) & (0,1) the target elements are triangle1 (0,0), (1,0) & - // (0,1) and triangle2 (1,0), (1,1) & (0,1) - - Omega_h::Write constant_field(source_mesh.nverts(), 2.0); - - auto load_vector = pcms::buildLoadVector(target_mesh, source_mesh, - intersection, constant_field); - - auto load_vector_host = Kokkos::create_mirror_view(load_vector); - Kokkos::deep_copy(load_vector_host, load_vector); - - double expected_load_vector[6] = {0.333333, 0.333333, 0.333333, - 0.333333, 0.333333, 0.333333}; - double tolerance = 1e-6; - for (int i = 0; i < load_vector_host.extent(0); ++i) { - - CAPTURE(i, expected_load_vector[i], load_vector_host[i], tolerance); - CHECK_THAT(expected_load_vector[i], - Catch::Matchers::WithinAbs(load_vector_host[i], tolerance)); - } - } -} From 0bb00c945d1572a7777fed141f256201fa7e2700 Mon Sep 17 00:00:00 2001 From: Jacob Merson Date: Sat, 11 Jul 2026 11:11:58 -0400 Subject: [PATCH 02/15] transfer: add integrator-based conservative (Galerkin) projection transfer Squashed original commits: bad8434, eca15d1, b4af186, 26f35c6, b8baf3b, 6b1d7f2, dba7966, 5baf365, 4b00ac1 (mesh_fields_rhs_integrator churn from 26f35c6/b8baf3b dropped) --- src/pcms/field/CMakeLists.txt | 1 + src/pcms/field/element_dispatch.h | 37 ++ src/pcms/field/layout/omega_h_lagrange.cpp | 5 + src/pcms/field/layout/omega_h_lagrange.h | 1 + src/pcms/transfer/CMakeLists.txt | 21 +- .../transfer/bilinear_form_integrator.hpp | 22 + .../conservative_projection_solver.cpp | 149 +++---- .../conservative_projection_solver.hpp | 116 ++--- src/pcms/transfer/integration_point_set.hpp | 55 +++ src/pcms/transfer/linear_form_integrator.hpp | 34 ++ src/pcms/transfer/mass_matrix_integrator.hpp | 6 +- .../omega_h_conservative_projection.cpp | 61 +-- .../omega_h_conservative_projection.hpp | 28 +- .../omega_h_form_integrator_utils.hpp | 339 +++++++++++++++ .../omega_h_intersection_rhs_integrator.cpp | 273 ++++++++++++ .../omega_h_intersection_rhs_integrator.hpp | 95 ++++ src/pcms/transfer/omega_h_mass_integrator.cpp | 138 ++++++ src/pcms/transfer/omega_h_mass_integrator.hpp | 45 ++ test/CMakeLists.txt | 5 +- .../test_mesh_intersection_field_transfer.cpp | 408 +++++++++++------- test/test_omega_h_form_integrator_utils.cpp | 80 ++++ ...st_omega_h_intersection_rhs_integrator.cpp | 337 +++++++++++++++ test/test_omega_h_mass_integrator.cpp | 187 ++++++++ 23 files changed, 2079 insertions(+), 364 deletions(-) create mode 100644 src/pcms/field/element_dispatch.h create mode 100644 src/pcms/transfer/bilinear_form_integrator.hpp create mode 100644 src/pcms/transfer/integration_point_set.hpp create mode 100644 src/pcms/transfer/linear_form_integrator.hpp create mode 100644 src/pcms/transfer/omega_h_form_integrator_utils.hpp create mode 100644 src/pcms/transfer/omega_h_intersection_rhs_integrator.cpp create mode 100644 src/pcms/transfer/omega_h_intersection_rhs_integrator.hpp create mode 100644 src/pcms/transfer/omega_h_mass_integrator.cpp create mode 100644 src/pcms/transfer/omega_h_mass_integrator.hpp create mode 100644 test/test_omega_h_form_integrator_utils.cpp create mode 100644 test/test_omega_h_intersection_rhs_integrator.cpp create mode 100644 test/test_omega_h_mass_integrator.cpp diff --git a/src/pcms/field/CMakeLists.txt b/src/pcms/field/CMakeLists.txt index 5d37a858..b37865be 100644 --- a/src/pcms/field/CMakeLists.txt +++ b/src/pcms/field/CMakeLists.txt @@ -3,6 +3,7 @@ find_package(KokkosKernels REQUIRED) set(PCMS_FIELD_HEADERS field.h coordinate_system.h + element_dispatch.h evaluation_request.h field_layout.h field_metadata.h diff --git a/src/pcms/field/element_dispatch.h b/src/pcms/field/element_dispatch.h new file mode 100644 index 00000000..67ac1fad --- /dev/null +++ b/src/pcms/field/element_dispatch.h @@ -0,0 +1,37 @@ +#ifndef PCMS_FIELD_ELEMENT_DISPATCH_H +#define PCMS_FIELD_ELEMENT_DISPATCH_H + +#include "pcms/utility/assert.h" +#include +#include + +namespace pcms::detail +{ + +// Runtime -> compile-time bridge for finite-element order. +// +// Turns a runtime Lagrange order into a compile-time std::integral_constant and +// invokes the generic callable `f` with it, so the callee can instantiate +// order-templated kernels (shape functions, DOF-per-element counts) without +// scattering `if (order == ...)` branches through device code. This is the same +// dispatch shape as MakeMeshFieldBackend's switch, factored out so both the +// transfer integrators and the MeshField evaluator backend can share it. +// +// Supported orders: 0 (one DOF per element) and 1 (one DOF per vertex). Add a +// case here when a higher-order element (P2, ...) lands; every consumer picks +// it up at once. `f` must return the same type for every supported order. +template +auto DispatchByOrder(int order, F&& f) +{ + switch (order) { + case 0: return f(std::integral_constant{}); + case 1: return f(std::integral_constant{}); + default: + throw pcms_error("DispatchByOrder: unsupported element order " + + std::to_string(order)); + } +} + +} // namespace pcms::detail + +#endif // PCMS_FIELD_ELEMENT_DISPATCH_H diff --git a/src/pcms/field/layout/omega_h_lagrange.cpp b/src/pcms/field/layout/omega_h_lagrange.cpp index a2749a6e..1b2f14a4 100644 --- a/src/pcms/field/layout/omega_h_lagrange.cpp +++ b/src/pcms/field/layout/omega_h_lagrange.cpp @@ -223,6 +223,11 @@ GlobalIDView OmegaHLagrangeLayout::GetGidsHost() const return GlobalIDView(gids_host_.data(), gids_host_.size()); } +GlobalIDView OmegaHLagrangeLayout::GetGids() const +{ + return GlobalIDView(gids_.data(), gids_.size()); +} + CoordinateView OmegaHLagrangeLayout::GetDOFHolderCoordinates() const { diff --git a/src/pcms/field/layout/omega_h_lagrange.h b/src/pcms/field/layout/omega_h_lagrange.h index f750b47a..68154979 100644 --- a/src/pcms/field/layout/omega_h_lagrange.h +++ b/src/pcms/field/layout/omega_h_lagrange.h @@ -35,6 +35,7 @@ class OmegaHLagrangeLayout : public FieldLayout Rank1View GetOwnedHost() const override; GlobalIDView GetGidsHost() const override; + GlobalIDView GetGids() const; CoordinateView GetDOFHolderCoordinates() const override; [[nodiscard]] bool IsDistributed() const override; diff --git a/src/pcms/transfer/CMakeLists.txt b/src/pcms/transfer/CMakeLists.txt index 216b5759..2de13931 100644 --- a/src/pcms/transfer/CMakeLists.txt +++ b/src/pcms/transfer/CMakeLists.txt @@ -7,6 +7,7 @@ set(PCMS_FIELD_TRANSFER_HEADERS interpolator.h copy.h transfer_operator.hpp + integration_point_set.hpp ) @@ -17,25 +18,29 @@ set(PCMS_FIELD_TRANSFER_SOURCES if(PCMS_ENABLE_MESHFIELDS) list(APPEND PCMS_FIELD_TRANSFER_HEADERS - load_vector_integrator.hpp mass_matrix_integrator.hpp) - list(APPEND PCMS_FIELD_TRANSFER_SOURCES - load_vector_integrator.cpp) endif() if(PCMS_ENABLE_PETSC AND PCMS_ENABLE_MESHFIELDS) list(APPEND PCMS_FIELD_TRANSFER_HEADERS - calculate_load_vector.hpp - calculate_mass_matrix.hpp + linear_form_integrator.hpp + bilinear_form_integrator.hpp conservative_projection_solver.hpp omega_h_conservative_projection.hpp + omega_h_intersection_rhs_integrator.hpp + omega_h_mass_integrator.hpp + omega_h_mc_rhs_integrator.hpp + omega_h_control_variate_projection.hpp coo_assembly_utils.hpp - petsc_utils.hpp) + petsc_utils.hpp + omega_h_form_integrator_utils.hpp) list(APPEND PCMS_FIELD_TRANSFER_SOURCES - calculate_load_vector.cpp - calculate_mass_matrix.cpp conservative_projection_solver.cpp omega_h_conservative_projection.cpp + omega_h_intersection_rhs_integrator.cpp + omega_h_mass_integrator.cpp + omega_h_mc_rhs_integrator.cpp + omega_h_control_variate_projection.cpp petsc_utils.cpp coo_assembly_utils.cpp) endif() diff --git a/src/pcms/transfer/bilinear_form_integrator.hpp b/src/pcms/transfer/bilinear_form_integrator.hpp new file mode 100644 index 00000000..060667bd --- /dev/null +++ b/src/pcms/transfer/bilinear_form_integrator.hpp @@ -0,0 +1,22 @@ +#ifndef PCMS_TRANSFER_BILINEAR_FORM_INTEGRATOR_H +#define PCMS_TRANSFER_BILINEAR_FORM_INTEGRATOR_H + +#include + +namespace pcms +{ + +class BilinearFormIntegrator +{ +public: + // Returns the internally-assembled matrix. The integrator owns the matrix + // and its lifetime must exceed any KSP that references it (PETSc's reference + // counting keeps the matrix alive until the KSP is destroyed). + virtual Mat GetMatrix() const noexcept = 0; + + virtual ~BilinearFormIntegrator() noexcept = default; +}; + +} // namespace pcms + +#endif // PCMS_TRANSFER_BILINEAR_FORM_INTEGRATOR_H diff --git a/src/pcms/transfer/conservative_projection_solver.cpp b/src/pcms/transfer/conservative_projection_solver.cpp index 148454f8..e15555e3 100644 --- a/src/pcms/transfer/conservative_projection_solver.cpp +++ b/src/pcms/transfer/conservative_projection_solver.cpp @@ -1,69 +1,12 @@ #include #include "pcms/transfer/conservative_projection_solver.hpp" +#include "pcms/utility/arrays.h" #include "pcms/transfer/petsc_utils.hpp" -#include "pcms/transfer/calculate_load_vector.hpp" -#include "pcms/transfer/calculate_mass_matrix.hpp" +#include namespace pcms { -/** - * @brief Solves a linear system Ax = b using PETSc's KSP solvers - * - * Uses PETSc's Krylov Subspace solvers to find x in Ax = b. - * The solver can be configured through PETSc runtime options. - * - * @param A The system matrix - * @param b The right-hand side vector - * @return Vec Solution vector x - */ -static Vec solveLinearSystem(Mat A, Vec b) -{ - PetscInt m, n; - PetscErrorCode ierr; - - ierr = MatGetSize(A, &m, &n); - CHKERRABORT(PETSC_COMM_WORLD, ierr); - - Vec x; - ierr = createSeqVec(PETSC_COMM_WORLD, n, &x); - CHKERRABORT(PETSC_COMM_WORLD, ierr); - - KSP ksp; - ierr = KSPCreate(PETSC_COMM_WORLD, &ksp); - CHKERRABORT(PETSC_COMM_WORLD, ierr); - - ierr = KSPSetOperators(ksp, A, A); - CHKERRABORT(PETSC_COMM_WORLD, ierr); - - ierr = KSPSetComputeSingularValues(ksp, PETSC_TRUE); - CHKERRABORT(PETSC_COMM_WORLD, ierr); - - ierr = KSPSetFromOptions(ksp); - CHKERRABORT(PETSC_COMM_WORLD, ierr); - - ierr = KSPSetUp(ksp); - CHKERRABORT(PETSC_COMM_WORLD, ierr); - - 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"); - } - - ierr = KSPDestroy(&ksp); - CHKERRABORT(PETSC_COMM_WORLD, ierr); - - return x; -} static Omega_h::Reals vecToOmegaHReals(Vec vec) { @@ -86,59 +29,71 @@ 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) -{ - if ((PetscInt)source_values.size() != - source_mesh.coords().size() / source_mesh.dim()) { - std::cerr << "ERROR: source_values size (" << source_values.size() - << ") doesn't match expected size (" - << source_mesh.coords().size() / source_mesh.dim() << ")" - << std::endl; - throw std::runtime_error("source_values length mismatch"); - } +// --------------------------------------------------------------------------- +// GalerkinProjectionSolver +// --------------------------------------------------------------------------- - Mat mass; - PetscErrorCode ierr = calculateMassMatrix(target_mesh, &mass); - CHKERRABORT(PETSC_COMM_WORLD, ierr); +GalerkinProjectionSolver::GalerkinProjectionSolver( + const BilinearFormIntegrator& mass_integrator, + LinearFormIntegrator& rhs_integrator) + : rhs_integrator_(&rhs_integrator) +{ + Mat A = mass_integrator.GetMatrix(); - Vec vec; - ierr = calculateLoadVector(target_mesh, source_mesh, intersection, - source_values, &vec); + PetscInt m = 0, n = 0; + PetscErrorCode ierr = MatGetSize(A, &m, &n); CHKERRABORT(PETSC_COMM_WORLD, ierr); + nverts_ = m; - Vec x = solveLinearSystem(mass, vec); - auto solution_vector = vecToOmegaHReals(x); + const std::size_t num_pts = + rhs_integrator_->GetIntegrationPoints().NumPoints(); + sampled_values_ = + Kokkos::View("rhs_sampled", num_pts, 1); - ierr = VecDestroy(&x); + ierr = KSPCreate(PETSC_COMM_WORLD, &ksp_); CHKERRABORT(PETSC_COMM_WORLD, ierr); - - ierr = MatDestroy(&mass); + ierr = KSPSetOperators(ksp_, A, A); CHKERRABORT(PETSC_COMM_WORLD, ierr); - - ierr = VecDestroy(&vec); + ierr = KSPSetFromOptions(ksp_); CHKERRABORT(PETSC_COMM_WORLD, ierr); + ierr = KSPSetUp(ksp_); + CHKERRABORT(PETSC_COMM_WORLD, ierr); +} - return solution_vector; +GalerkinProjectionSolver::~GalerkinProjectionSolver() +{ + if (ksp_) { + KSPDestroy(&ksp_); + } + // mat_ is owned by the BilinearFormIntegrator; KSP holds its own reference. +} + +Omega_h::Reals GalerkinProjectionSolver::Solve( + const PointEvaluator& evaluator, const Field& source_field) const +{ + evaluator.Evaluate(source_field, MakeRank2View(sampled_values_)); + return Solve(MakeConstRank2View(sampled_values_)); } -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 GalerkinProjectionSolver::Solve( + Rank2View sampled_values) const { - Vec vec; - PetscErrorCode ierr; - ierr = calculateLoadVector(target_mesh, source_mesh, intersection, - source_values, &vec); + rhs_integrator_->Assemble(sampled_values); + Vec rhs_vector = rhs_integrator_->GetVector(); + + Vec solution = nullptr; + PetscErrorCode ierr = createSeqVec(PETSC_COMM_WORLD, nverts_, &solution); CHKERRABORT(PETSC_COMM_WORLD, ierr); - auto rhsvector = vecToOmegaHReals(vec); + ierr = KSPSolve(ksp_, rhs_vector, solution); + CHKERRABORT(PETSC_COMM_WORLD, ierr); - ierr = VecDestroy(&vec); + auto result = vecToOmegaHReals(solution); + + ierr = VecDestroy(&solution); CHKERRABORT(PETSC_COMM_WORLD, ierr); - return rhsvector; + return result; } + } // namespace pcms diff --git a/src/pcms/transfer/conservative_projection_solver.hpp b/src/pcms/transfer/conservative_projection_solver.hpp index d3e46a37..f0f2df0e 100644 --- a/src/pcms/transfer/conservative_projection_solver.hpp +++ b/src/pcms/transfer/conservative_projection_solver.hpp @@ -1,72 +1,72 @@ -/** - * @file conservative_projection_solver.hpp - * @brief Solves the conservative projection of scalar fields between - * non-matching meshes. - * - * Provides the main interface to perform Galerkin projection of scalar fields - * from a source mesh to a target mesh using conservative transfer using a - * supermesh generated from mesh intersections. - * - * The solver computes the right-hand side (load vector), assembles the mass - * matrix, and solves the resulting linear system to obtain projected nodal - * values. - * - */ - #ifndef PCMS_TRANSFER_CONSERVATIVE_PROJECTION_SOLVER_HPP #define PCMS_TRANSFER_CONSERVATIVE_PROJECTION_SOLVER_HPP +#include #include -#include +#include -#include +#include +#include +#include +#include +#include namespace pcms { -/** - * @brief Solves a conservative galerkin projection problem to transfer scalar - * field values onto a target mesh. - * - * This function assembles and solves a linear system of the form: - * \f[ - * M \cdot x = f - * \f] - * where: - * - \f$M\f$ is the mass matrix on the target mesh (based on P1 finite - * elements), - * - \f$f\f$ is the load vector computed on the supermesh, - * - \f$x\f$ is the unknown nodal field on the target mesh (solution). - * - * The method computes the conservative field transfer between two non-matching - * meshes using mesh intersections (supermesh). - * - * ### Algorithm Steps: - * 1. Compute and assemble mass matrix and load vector - * 2. Solve the linear system using PETSc. - * 3. Return the solution as a nodal field on the target mesh. - * - * @param target_mesh The Omega_h mesh where the field is projected. - * @param source_mesh The Omega_h mesh containing the original field data. - * @param intersection Precomputed intersection information between source and - * target meshes. - * @param source_values Nodal scalar field values on the source mesh. - * - * @return A vector of nodal values on the target mesh after projection - * (Omega_h::Reals). - * - * - */ +// Cached Galerkin projection solver for M*x = f. +// +// The mass matrix is obtained from mass_integrator (which owns and assembled +// it) and the KSP is set up (preconditioned/factored) at construction. The RHS +// integrator is also bound at construction: its integration points fix the +// size of the cached sample buffer, which is allocated once and reused. Each +// Solve() call only fills that buffer, assembles a fresh RHS vector, and calls +// KSPSolve, reusing the cached factorization. +// +// Binding one RHS integrator per solver matches every current consumer (each +// builds M and the integrator together and reuses them in lockstep) and lets +// the sample buffer be sized exactly once. Reusing a factored M with a +// different integrator is intentionally not supported; if that need arises, +// add a setter that rebinds the integrator and resizes the buffer. +// +// Non-copyable and non-movable because it owns PETSc handles. +class GalerkinProjectionSolver +{ +public: + // Retrieves the assembled matrix from mass_integrator, sets up the KSP, and + // sizes the sample buffer from rhs_integrator's integration points. PETSc + // reference-counts the matrix, so mass_integrator need not outlive this + // solver, but rhs_integrator must (it is assembled into on every Solve). + GalerkinProjectionSolver(const BilinearFormIntegrator& mass_integrator, + LinearFormIntegrator& rhs_integrator); + ~GalerkinProjectionSolver(); + + GalerkinProjectionSolver(const GalerkinProjectionSolver&) = delete; + GalerkinProjectionSolver& operator=(const GalerkinProjectionSolver&) = delete; + + // Evaluates source_field at the bound integrator's integration points via + // evaluator into the cached buffer, assembles the RHS, solves M*x = f using + // the cached KSP, and returns nodal values on the target mesh. + // + // May be called repeatedly with different source fields, provided the target + // mesh (and therefore M) is unchanged. + Omega_h::Reals Solve(const PointEvaluator& evaluator, + const Field& source_field) const; + + // Assembles the RHS from values already sampled at the bound integrator's + // integration points (shape [num_points][num_components]) and solves + // M*x = f using the cached KSP. Used by workflows that pre-process the + // sampled values (e.g. control-variate residuals) before assembly. + Omega_h::Reals Solve( + Rank2View sampled_values) const; -Omega_h::Reals solveGalerkinProjection(Omega_h::Mesh& target_mesh, - Omega_h::Mesh& source_mesh, - const IntersectionResults& intersection, - const Omega_h::Reals& source_values); +private: + PetscInt nverts_; + KSP ksp_ = nullptr; + LinearFormIntegrator* rhs_integrator_; + mutable Kokkos::View sampled_values_; +}; -Omega_h::Reals rhsVectorMI(Omega_h::Mesh& target_mesh, - Omega_h::Mesh& source_mesh, - const IntersectionResults& intersection, - const Omega_h::Reals& source_values); } // namespace pcms #endif // PCMS_TRANSFER_CONSERVATIVE_PROJECTION_SOLVER_HPP diff --git a/src/pcms/transfer/integration_point_set.hpp b/src/pcms/transfer/integration_point_set.hpp new file mode 100644 index 00000000..3a91c9e1 --- /dev/null +++ b/src/pcms/transfer/integration_point_set.hpp @@ -0,0 +1,55 @@ +#ifndef PCMS_TRANSFER_INTEGRATION_POINT_SET_H +#define PCMS_TRANSFER_INTEGRATION_POINT_SET_H + +#include "pcms/field/coordinate_system.h" +#include "pcms/utility/arrays.h" +#include "pcms/utility/types.h" + +#include +#include +#include + +namespace pcms +{ + +// An ordered collection of global coordinates at which a source field must +// be evaluated for integration. The ordering is significant: sampled values +// passed back to a LinearFormIntegrator must use the same point order as +// stored here. +// +// This type intentionally does not expose quadrature weights, basis values, +// entity ids, or any other assembly metadata. Those remain private to the +// concrete integrator implementation. It owns the coordinate buffer; callers +// receive a non-owning CoordinateView from GetCoordinates(). +template +class IntegrationPointSet +{ +public: + IntegrationPointSet(CoordinateSystem coordinate_system, + Kokkos::View coords) + : coordinate_system_(coordinate_system), coords_(std::move(coords)) + { + } + + // Global coordinates of all integration points, shape [num_points][dim]. + CoordinateView GetCoordinates() const + { + return CoordinateView(coordinate_system_, + MakeConstRank2View(coords_)); + } + + // Number of integration points. + [[nodiscard]] + std::size_t NumPoints() const + { + return coords_.extent(0); + } + +private: + CoordinateSystem coordinate_system_; + Kokkos::View coords_; +}; + +} // namespace pcms + +#endif // PCMS_TRANSFER_INTEGRATION_POINT_SET_H diff --git a/src/pcms/transfer/linear_form_integrator.hpp b/src/pcms/transfer/linear_form_integrator.hpp new file mode 100644 index 00000000..38835a8a --- /dev/null +++ b/src/pcms/transfer/linear_form_integrator.hpp @@ -0,0 +1,34 @@ +#ifndef PCMS_TRANSFER_LINEAR_FORM_INTEGRATOR_H +#define PCMS_TRANSFER_LINEAR_FORM_INTEGRATOR_H + +#include "pcms/transfer/integration_point_set.hpp" +#include "pcms/utility/arrays.h" +#include + +namespace pcms +{ + +class LinearFormIntegrator +{ +public: + // Returns the ordered integration points. The index ordering must match the + // first dimension of sampled_values passed to Assemble. + virtual const IntegrationPointSet& GetIntegrationPoints() + const noexcept = 0; + + // Returns the internally-owned vector. The integrator owns the vector; + // callers receive a non-owning handle. Valid after construction. + virtual Vec GetVector() const noexcept = 0; + + // Zeros the owned vector, scatter-adds weighted contributions from the + // sampled source-field values, and finalizes assembly. + // sampled_values shape: [num_integration_points][num_components] + virtual void Assemble( + Rank2View sampled_values) = 0; + + virtual ~LinearFormIntegrator() noexcept = default; +}; + +} // namespace pcms + +#endif // PCMS_TRANSFER_LINEAR_FORM_INTEGRATOR_H diff --git a/src/pcms/transfer/mass_matrix_integrator.hpp b/src/pcms/transfer/mass_matrix_integrator.hpp index 1a2b4b59..7d36c975 100644 --- a/src/pcms/transfer/mass_matrix_integrator.hpp +++ b/src/pcms/transfer/mass_matrix_integrator.hpp @@ -31,7 +31,7 @@ class MassMatrixIntegrator : public MeshField::Integrator } void atPoints(Kokkos::View p, Kokkos::View w, - Kokkos::View dV) + Kokkos::View dV) override { // std::cerr << "MassMatrixIntegrator::atPoints(...)\n"; const size_t numPtsPerElem = p.extent(0) / mesh.nelems(); @@ -77,8 +77,8 @@ class MassMatrixIntegrator : public MeshField::Integrator }; template -Kokkos::View buildMassMatrix(Omega_h::Mesh& mesh, - FieldElement& coordFe) +Kokkos::View buildElementMassMatrix(Omega_h::Mesh& mesh, + FieldElement& coordFe) { MassMatrixIntegrator mmi(mesh, coordFe); mmi.process(coordFe); diff --git a/src/pcms/transfer/omega_h_conservative_projection.cpp b/src/pcms/transfer/omega_h_conservative_projection.cpp index c7ae7e65..f5d911ec 100644 --- a/src/pcms/transfer/omega_h_conservative_projection.cpp +++ b/src/pcms/transfer/omega_h_conservative_projection.cpp @@ -1,6 +1,7 @@ #include "pcms/transfer/omega_h_conservative_projection.hpp" +#include "pcms/transfer/conservative_projection_solver.hpp" +#include "pcms/transfer/omega_h_mass_integrator.hpp" #include "pcms/utility/arrays.h" -#include "pcms/utility/assert.h" #include #include @@ -10,28 +11,6 @@ namespace pcms namespace { -void CheckSupportedLayout( - const FunctionSpace& space, - const std::shared_ptr& layout, const char* role) -{ - if (layout == nullptr) { - throw pcms_error(std::string("OmegaHConservativeProjection: ") + role + - " space must use OmegaHLagrangeLayout"); - } - if (layout->GetOrder() != 1) { - throw pcms_error(std::string("OmegaHConservativeProjection: ") + role + - " space must be order-1"); - } - if (layout->GetNumComponents() != 1) { - throw pcms_error(std::string("OmegaHConservativeProjection: ") + role + - " space must have exactly one component"); - } - if (space.GetCoordinateSystem() != CoordinateSystem::Cartesian) { - throw pcms_error(std::string("OmegaHConservativeProjection: ") + role + - " space must use Cartesian coordinates"); - } -} - void CheckApplyCompatible(const Field& source, const Field& target, const OmegaHLagrangeLayout& source_layout, const OmegaHLagrangeLayout& target_layout) @@ -58,15 +37,6 @@ void CheckApplyCompatible(const Field& source, const Field& target, } } -Omega_h::Reals MakeOmegaHReals(Rank1View values) -{ - Omega_h::HostWrite values_host(values.size()); - for (size_t i = 0; i < values.size(); ++i) { - values_host[i] = values[i]; - } - return Omega_h::Reals(values_host); -} - } // namespace OmegaHConservativeProjection::OmegaHConservativeProjection( @@ -76,23 +46,32 @@ OmegaHConservativeProjection::OmegaHConservativeProjection( target_layout_(std::dynamic_pointer_cast( target_space.GetLayout())) { - CheckSupportedLayout(source_space, source_layout_, "source"); - CheckSupportedLayout(target_space, target_layout_, "target"); + rhs_integrator_ = std::make_unique( + source_layout_, source_space.GetCoordinateSystem(), target_layout_, + target_space.GetCoordinateSystem()); + + evaluator_ = + source_space.CreatePointEvaluator(EvaluationRequest::FromCoordinates( + rhs_integrator_->GetIntegrationPoints().GetCoordinates())); - intersections_ = - intersectTargets(source_layout_->GetMesh(), target_layout_->GetMesh()); + // 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()); + solver_ = std::make_unique(mass_integrator, + *rhs_integrator_); } +// Defined here so that GalerkinProjectionSolver (forward-declared in the +// header) is a complete type when unique_ptr's destructor is instantiated. +OmegaHConservativeProjection::~OmegaHConservativeProjection() = default; + void OmegaHConservativeProjection::Apply(const Field& source, Field& target) const { CheckApplyCompatible(source, target, *source_layout_, *target_layout_); - const auto source_values = - MakeOmegaHReals(FlattenToRank1View(source.GetDOFHolderDataHost())); - const auto target_values = solveGalerkinProjection( - target_layout_->GetMesh(), source_layout_->GetMesh(), intersections_, - source_values); + const auto target_values = solver_->Solve(*evaluator_, source); auto target_values_h = Omega_h::HostRead(target_values); target.SetDOFHolderDataHost(Rank2View( target_values_h.data(), target_layout_->GetNumOwnedDofHolder(), diff --git a/src/pcms/transfer/omega_h_conservative_projection.hpp b/src/pcms/transfer/omega_h_conservative_projection.hpp index 386cb260..114fc2a0 100644 --- a/src/pcms/transfer/omega_h_conservative_projection.hpp +++ b/src/pcms/transfer/omega_h_conservative_projection.hpp @@ -3,30 +3,48 @@ #include "pcms/field/function_space.h" #include "pcms/field/layout/omega_h_lagrange.h" -#include "pcms/transfer/conservative_projection_solver.hpp" +#include "pcms/field/point_evaluator.h" +#include "pcms/transfer/omega_h_intersection_rhs_integrator.hpp" #include "pcms/transfer/transfer_operator.hpp" #include namespace pcms { +// Forward declaration keeps PETSc headers out of this header. +class GalerkinProjectionSolver; + // Conservative Galerkin projection between Omega_h order-1 Lagrange spaces. // -// Construction caches the source-target mesh intersections. Apply() assembles -// the projection system for the current source coefficients, solves it, and -// writes the projected scalar field into the target Field. +// Construction (one-time cost): +// - Computes mesh intersections and quadrature data +// (OmegaHIntersectionRHSIntegrator). +// - Localizes integration points in the source mesh (PointEvaluator). +// - Assembles the target mass matrix and factors it via KSP +// (GalerkinProjectionSolver). +// +// Apply() (per-call cost): +// - Evaluates the source field at the fixed integration points. +// - Assembles a fresh RHS vector. +// - Calls KSPSolve reusing the cached factorization. +// No mesh intersection, no matrix assembly, no refactorization per call. class OmegaHConservativeProjection : public TransferOperator { public: OmegaHConservativeProjection(const FunctionSpace& source_space, const FunctionSpace& target_space); + // Defined in the .cpp where GalerkinProjectionSolver is a complete type. + ~OmegaHConservativeProjection() override; + void Apply(const Field& source, Field& target) const override; private: std::shared_ptr source_layout_; std::shared_ptr target_layout_; - IntersectionResults intersections_; + std::unique_ptr rhs_integrator_; + std::unique_ptr> evaluator_; + std::unique_ptr solver_; }; } // namespace pcms diff --git a/src/pcms/transfer/omega_h_form_integrator_utils.hpp b/src/pcms/transfer/omega_h_form_integrator_utils.hpp new file mode 100644 index 00000000..3294e45e --- /dev/null +++ b/src/pcms/transfer/omega_h_form_integrator_utils.hpp @@ -0,0 +1,339 @@ +#ifndef PCMS_TRANSFER_OMEGA_H_FORM_INTEGRATOR_UTILS_HPP +#define PCMS_TRANSFER_OMEGA_H_FORM_INTEGRATOR_UTILS_HPP + +#include "pcms/field/function_space.h" +#include "pcms/field/layout/omega_h_lagrange.h" +#include "pcms/transfer/mesh_intersection.hpp" +#include "pcms/utility/assert.h" +#include +#include +#include +#include + +namespace pcms::detail +{ + +// Shared checks for a scalar Cartesian Lagrange space on a 2D simplex mesh, +// independent of order. Order is validated separately by the callers below. +inline void CheckOmegaHScalarSimplex2DLayout( + CoordinateSystem coordinate_system, + const std::shared_ptr& layout, + const char* context, const char* role) +{ + if (layout == nullptr) { + throw pcms_error(std::string(context) + ": " + role + + " space must use OmegaHLagrangeLayout"); + } + if (layout->GetNumComponents() != 1) { + throw pcms_error(std::string(context) + ": " + role + + " space must have exactly one component"); + } + if (coordinate_system != CoordinateSystem::Cartesian) { + throw pcms_error(std::string(context) + ": " + role + + " space must use Cartesian coordinates"); + } + const Omega_h::Mesh& mesh = layout->GetMesh(); + if (mesh.dim() != 2) { + throw pcms_error(std::string(context) + ": " + role + " mesh must be 2D"); + } + if (mesh.family() != OMEGA_H_SIMPLEX) { + throw pcms_error(std::string(context) + ": " + role + + " mesh must be a simplex (triangle) mesh"); + } +} + +// Strict order-1 check (used where only P1 is supported, e.g. Monte-Carlo RHS). +inline void CheckOmegaHScalarP1Layout( + CoordinateSystem coordinate_system, + const std::shared_ptr& layout, + const char* context, const char* role) +{ + CheckOmegaHScalarSimplex2DLayout(coordinate_system, layout, context, role); + if (layout->GetOrder() != 1) { + throw pcms_error(std::string(context) + ": " + role + + " space must be order-1"); + } +} + +// Conservative-projection check: any supported Lagrange order (P0 or P1). The +// intersection integrator handles source and target orders independently, so +// this replaces the strict P1 requirement on those paths. +inline void CheckOmegaHScalarLagrangeLayout( + CoordinateSystem coordinate_system, + const std::shared_ptr& layout, + const char* context, const char* role) +{ + CheckOmegaHScalarSimplex2DLayout(coordinate_system, layout, context, role); + const int order = layout->GetOrder(); + if (order != 0 && order != 1) { + throw pcms_error(std::string(context) + ": " + role + + " space must be order-0 or order-1"); + } +} + +[[nodiscard]] OMEGA_H_INLINE Omega_h::Vector<2> GlobalFromBarycentric( + const MeshField::Vector3& barycentric_coord, + const Omega_h::Few, 3>& verts_coord) +{ + Omega_h::Vector<2> real_coords = {0.0, 0.0}; + for (int i = 0; i < 3; ++i) { + real_coords[0] += barycentric_coord[i] * verts_coord[i][0]; + real_coords[1] += barycentric_coord[i] * verts_coord[i][1]; + } + return real_coords; +} + +[[nodiscard]] OMEGA_H_INLINE Omega_h::Vector<3> EvaluateBarycentric( + const Omega_h::Vector<2>& point, + const r3d::Few, 3>& verts_coord) +{ + Omega_h::Few, 3> omegah_vector; + for (int i = 0; i < 3; ++i) { + omegah_vector[i][0] = verts_coord[i][0]; + omegah_vector[i][1] = verts_coord[i][1]; + } + return Omega_h::barycentric_from_global<2, 2>(point, omegah_vector); +} + +[[nodiscard]] OMEGA_H_INLINE int RemoveDuplicateVerticesAndFixLinks( + r3d::Polytope<2>& poly, const double tol = 1e-12) +{ + const int old_n = poly.nverts; + int new_n = 0; + + for (int i = 0; i < old_n; ++i) { + const auto& pi = poly.verts[i].pos; + bool dup = false; + for (int j = 0; j < new_n; ++j) { + const auto& pj = poly.verts[j].pos; + if (Kokkos::fabs(pi[0] - pj[0]) < tol && + Kokkos::fabs(pi[1] - pj[1]) < tol) { + dup = true; + break; + } + } + if (!dup) { + poly.verts[new_n] = poly.verts[i]; + ++new_n; + } + } + poly.nverts = new_n; + + if (new_n >= 3) { + double cx = 0.0; + double cy = 0.0; + for (int i = 0; i < new_n; ++i) { + cx += poly.verts[i].pos[0]; + cy += poly.verts[i].pos[1]; + } + cx /= new_n; + cy /= new_n; + + int order[r3d::MaxVerts<2>::value]; + for (int i = 0; i < new_n; ++i) { + order[i] = i; + } + + for (int i = 0; i < new_n - 1; ++i) { + for (int j = i + 1; j < new_n; ++j) { + const double a1 = Kokkos::atan2(poly.verts[order[i]].pos[1] - cy, + poly.verts[order[i]].pos[0] - cx); + const double a2 = Kokkos::atan2(poly.verts[order[j]].pos[1] - cy, + poly.verts[order[j]].pos[0] - cx); + if (a1 > a2) { + int t = order[i]; + order[i] = order[j]; + order[j] = t; + } + } + } + + r3d::Vertex<2> tmp[r3d::MaxVerts<2>::value]; + for (int i = 0; i < new_n; ++i) { + tmp[i] = poly.verts[order[i]]; + } + + for (int i = 0; i < new_n; ++i) { + poly.verts[i] = tmp[i]; + } + + for (int i = 0; i < new_n; ++i) { + const int prev = (i - 1 + new_n) % new_n; + const int next = (i + 1) % new_n; + poly.verts[i].pnbrs[0] = prev; + poly.verts[i].pnbrs[1] = next; + } + + double area = r3d::measure(poly); + if (area < 0.0) { + for (int i = 0; i < new_n / 2; ++i) { + r3d::Vertex<2> t = poly.verts[i]; + poly.verts[i] = poly.verts[new_n - 1 - i]; + poly.verts[new_n - 1 - i] = t; + } + for (int i = 0; i < new_n; ++i) { + const int prev = (i - 1 + new_n) % new_n; + const int next = (i + 1) % new_n; + poly.verts[i].pnbrs[0] = prev; + poly.verts[i].pnbrs[1] = next; + } + } + } else { + for (int i = 0; i < new_n; ++i) { + poly.verts[i].pnbrs[0] = -1; + poly.verts[i].pnbrs[1] = -1; + } + } + + return new_n; +} + +template +OMEGA_H_INLINE void ForEachIntersectionSubtriangle( + const int elm, const IntersectionResults& intersection, + const Omega_h::Reals& tgt_coords, const Omega_h::Reals& src_coords, + const Omega_h::LOs& tgt_faces2nodes, const Omega_h::LOs& src_faces2nodes, + TriangleOp&& op) +{ + auto tgt_elm_vert_coords = + get_vert_coords_of_elem(tgt_coords, tgt_faces2nodes, elm); + const int start = intersection.tgt2src_offsets[elm]; + const int end = intersection.tgt2src_offsets[elm + 1]; + + for (int i = start; i < end; ++i) { + const int current_src_elm = intersection.tgt2src_indices[i]; + auto src_elm_vert_coords = + get_vert_coords_of_elem(src_coords, src_faces2nodes, current_src_elm); + r3d::Polytope<2> poly; + r3d::intersect_simplices(poly, tgt_elm_vert_coords, src_elm_vert_coords); + auto nverts = RemoveDuplicateVerticesAndFixLinks(poly, 1e-12); + auto poly_area = r3d::measure(poly); + + for (int j = 1; j < nverts - 1; ++j) { + auto& p0 = poly.verts[0].pos; + auto& p1 = poly.verts[j].pos; + auto& p2 = poly.verts[j + 1].pos; + + Omega_h::Few, 3> tri_coords; + tri_coords[0] = {p0[0], p0[1]}; + tri_coords[1] = {p1[0], p1[1]}; + tri_coords[2] = {p2[0], p2[1]}; + + Omega_h::Few, 2> basis; + basis[0] = tri_coords[1] - tri_coords[0]; + basis[1] = tri_coords[2] - tri_coords[0]; + + Omega_h::Real area = + Kokkos::fabs(Omega_h::triangle_area_from_basis(basis)); + + const double eps_area = abs_tol + rel_tol * poly_area; + if (area <= eps_area) { + continue; + } + + op(tri_coords, tgt_elm_vert_coords, src_elm_vert_coords, current_src_elm, + area); + } + } +} + +// Barycentric integration points and weights for a reference triangle, taken +// from MeshField's predefined triangle quadrature rules and staged on device +// for use in element integration kernels. +// +// MeshField::getIntegrationPoints returns a host std::vector, which cannot be +// dereferenced inside a device kernel, so the (tiny) rule is copied into device +// Kokkos views once at construction. +// +// The quadrature order is a runtime argument because the required polynomial +// accuracy depends on the source and target element orders (degree = +// source_order + target_order), which are only known at construction. +struct IntegrationData +{ + Kokkos::View bary_coords; // barycentric coordinates + Kokkos::View weights; // quadrature weights + + explicit IntegrationData(int order) + { + auto ip_vec = MeshField::getIntegrationPoints(order); + const std::size_t num_ip = ip_vec.size(); + + bary_coords = Kokkos::View("bary_coords", num_ip); + weights = Kokkos::View("weights", num_ip); + + auto bary_coords_host = Kokkos::create_mirror_view(bary_coords); + auto weights_host = Kokkos::create_mirror_view(weights); + for (std::size_t i = 0; i < num_ip; ++i) { + bary_coords_host(i) = ip_vec[i].param; + weights_host(i) = ip_vec[i].weight; + } + Kokkos::deep_copy(bary_coords, bary_coords_host); + Kokkos::deep_copy(weights, weights_host); + } + + int size() const { return bary_coords.extent(0); } +}; + +// Target Lagrange basis on a triangle, parameterized by element order, for the +// conservative-projection RHS assembly. Order 0 is a single element-constant +// DOF; order 1 is the three vertex (barycentric) DOFs. Higher orders slot in as +// additional specializations, kept in lock-step with element_dispatch.h. +// +// Each specialization provides, for a target element `elm` with local vertex +// ids `verts` and vertex coordinates `tgt_verts`: +// ndof number of local target DOFs +// Gid(...) global DOF id for local dof k (element id for P0, vertex +// id for P1) — matching OmegaHLagrangeLayout::GetGids() +// Values(pt, ...) basis values at the (global) integration point pt +template +struct TargetTriBasis; + +template <> +struct TargetTriBasis<0> +{ + static constexpr int ndof = 1; + + KOKKOS_INLINE_FUNCTION static GO Gid(const GO* gids, int elm, + const Omega_h::Few&, + int /*k*/) + { + return gids[elm]; + } + + KOKKOS_INLINE_FUNCTION static void Values( + const Omega_h::Vector<2>&, const Omega_h::Few, 3>&, + Omega_h::Real out[ndof]) + { + out[0] = 1.0; + } +}; + +template <> +struct TargetTriBasis<1> +{ + static constexpr int ndof = 3; + + KOKKOS_INLINE_FUNCTION static GO Gid( + const GO* gids, int /*elm*/, const Omega_h::Few& verts, + int k) + { + return gids[verts[k]]; + } + + // P1 basis functions are the barycentric coordinates of the target element + // evaluated at the (global) integration point. + KOKKOS_INLINE_FUNCTION static void Values( + const Omega_h::Vector<2>& pt, + const Omega_h::Few, 3>& tgt_verts, + Omega_h::Real out[ndof]) + { + const auto bary = Omega_h::barycentric_from_global<2, 2>(pt, tgt_verts); + out[0] = bary[0]; + out[1] = bary[1]; + out[2] = bary[2]; + } +}; + +} // namespace pcms::detail + +#endif // PCMS_TRANSFER_OMEGA_H_FORM_INTEGRATOR_UTILS_HPP diff --git a/src/pcms/transfer/omega_h_intersection_rhs_integrator.cpp b/src/pcms/transfer/omega_h_intersection_rhs_integrator.cpp new file mode 100644 index 00000000..10103981 --- /dev/null +++ b/src/pcms/transfer/omega_h_intersection_rhs_integrator.cpp @@ -0,0 +1,273 @@ +#include "pcms/transfer/omega_h_intersection_rhs_integrator.hpp" +#include "pcms/field/element_dispatch.h" +#include "pcms/transfer/omega_h_form_integrator_utils.hpp" +#include "pcms/transfer/petsc_utils.hpp" +#include "pcms/utility/assert.h" +#include +#include +#include +#include + +namespace pcms +{ + +// --------------------------------------------------------------------------- +// OmegaHIntersectionRHSIntegrator::BuildData +// --------------------------------------------------------------------------- + +OmegaHIntersectionRHSIntegrator::Data +OmegaHIntersectionRHSIntegrator::BuildData(const FunctionSpace& source_space, + const FunctionSpace& target_space) +{ + return BuildData(std::dynamic_pointer_cast( + source_space.GetLayout()), + source_space.GetCoordinateSystem(), + std::dynamic_pointer_cast( + target_space.GetLayout()), + target_space.GetCoordinateSystem()); +} + +OmegaHIntersectionRHSIntegrator::Data +OmegaHIntersectionRHSIntegrator::BuildData( + std::shared_ptr source_layout, + CoordinateSystem source_coordinate_system, + std::shared_ptr target_layout, + CoordinateSystem target_coordinate_system) +{ + detail::CheckOmegaHScalarLagrangeLayout( + source_coordinate_system, source_layout, "OmegaHIntersectionRHSIntegrator", + "source"); + detail::CheckOmegaHScalarLagrangeLayout( + target_coordinate_system, target_layout, "OmegaHIntersectionRHSIntegrator", + "target"); + + // The integrand f_src * phi_target has polynomial degree source_order + + // target_order on each intersection subtriangle; integrate it exactly (with a + // 1-point floor so a P0->P0 pair still gets a valid rule). + const int quad_order = + std::max(1, source_layout->GetOrder() + target_layout->GetOrder()); + + return detail::DispatchByOrder(target_layout->GetOrder(), [&](auto order_c) { + constexpr int TgtOrder = decltype(order_c)::value; + return BuildDataImpl(*source_layout, *target_layout, quad_order); + }); +} + +template +OmegaHIntersectionRHSIntegrator::Data +OmegaHIntersectionRHSIntegrator::BuildDataImpl( + const OmegaHLagrangeLayout& source_layout, + const OmegaHLagrangeLayout& target_layout, int quad_order) +{ + using Basis = detail::TargetTriBasis; + constexpr int ndof = Basis::ndof; + + Omega_h::Mesh& source_mesh = source_layout.GetMesh(); + Omega_h::Mesh& target_mesh = target_layout.GetMesh(); + + const auto intersections = intersectTargets(source_mesh, target_mesh); + + const auto& tgt_coords = target_mesh.coords(); + const auto& tgt_faces2nodes = + target_mesh.ask_down(Omega_h::FACE, Omega_h::VERT).ab2b; + const auto& src_coords = source_mesh.coords(); + const auto& src_faces2nodes = + source_mesh.ask_down(Omega_h::FACE, Omega_h::VERT).ab2b; + + detail::IntegrationData ip_data(quad_order); + const int npts = ip_data.size(); + auto bary_coords = ip_data.bary_coords; // device view + auto weights = ip_data.weights; // device view + + const auto device_gids = target_layout.GetGids(); + const GO* gids_ptr = device_gids.data_handle(); + + const Omega_h::LOs tgt2src_offsets = intersections.tgt2src_offsets; + const Omega_h::LOs tgt2src_indices = intersections.tgt2src_indices; + const int nelems = target_mesh.nelems(); + + // Pass 1: count integration points per target element. + Omega_h::Write ip_counts(nelems, 0, "rhs_ip_counts"); + Kokkos::parallel_for( + "rhs_count", nelems, KOKKOS_LAMBDA(int elm) { + int count = 0; + detail::ForEachIntersectionSubtriangle( + elm, {tgt2src_offsets, tgt2src_indices}, tgt_coords, src_coords, + tgt_faces2nodes, src_faces2nodes, + [&](const Omega_h::Few, 3>&, + const r3d::Few, 3>&, + const r3d::Few, 3>&, int, + Omega_h::Real) { count += npts; }); + ip_counts[elm] = count; + }); + Kokkos::fence(); + + const auto ip_offsets = Omega_h::offset_scan( + Omega_h::Read(ip_counts), "rhs_ip_offsets"); + const int num_pts = static_cast(ip_offsets.last()); + + // Pass 2: fill coords, node_gids, and coeffs on device. node_gids/coeffs hold + // ndof (target DOFs per element) entries per integration point. + Kokkos::View coords("rhs_coords", num_pts, 2); + Kokkos::View node_gids("rhs_node_gids", + num_pts * ndof); + Kokkos::View coeffs("rhs_coeffs", num_pts * ndof); + + Kokkos::parallel_for( + "rhs_fill", nelems, KOKKOS_LAMBDA(int elm) { + const auto tgt_verts = Omega_h::gather_verts<3>(tgt_faces2nodes, elm); + const Omega_h::Matrix<2, 3> tgt_vert_mat = + Omega_h::gather_vectors<3, 2>(tgt_coords, tgt_verts); + Omega_h::Few, 3> tgt_omh; + for (int i = 0; i < 3; ++i) + tgt_omh[i] = {tgt_vert_mat[i][0], tgt_vert_mat[i][1]}; + + int ip_local = 0; + const int offset = ip_offsets[elm]; + + detail::ForEachIntersectionSubtriangle( + elm, {tgt2src_offsets, tgt2src_indices}, tgt_coords, src_coords, + tgt_faces2nodes, src_faces2nodes, + [&](const Omega_h::Few, 3>& tri, + const r3d::Few, 3>&, + const r3d::Few, 3>&, int, Omega_h::Real area) { + for (int ip_idx = 0; ip_idx < npts; ++ip_idx) { + const auto bary = bary_coords(ip_idx); + const double w = weights(ip_idx); + const auto pt = detail::GlobalFromBarycentric(bary, tri); + + Omega_h::Real basis[ndof]; + Basis::Values(pt, tgt_omh, basis); + + const int global_ip = offset + ip_local; + coords(global_ip, 0) = pt[0]; + coords(global_ip, 1) = pt[1]; + for (int k = 0; k < ndof; ++k) { + node_gids(global_ip * ndof + k) = + static_cast(Basis::Gid(gids_ptr, elm, tgt_verts, k)); + coeffs(global_ip * ndof + k) = basis[k] * w * 2.0 * area; + } + ++ip_local; + } + }); + }); + Kokkos::fence(); + + Data d; + d.coords = std::move(coords); + d.node_gids = std::move(node_gids); + d.coeffs = std::move(coeffs); + d.ndof_per_elem = ndof; + d.num_target_dofs = + static_cast(target_layout.GetNumGlobalDofHolder()); + return d; +} + +// --------------------------------------------------------------------------- +// OmegaHIntersectionRHSIntegrator constructors +// --------------------------------------------------------------------------- + +OmegaHIntersectionRHSIntegrator::OmegaHIntersectionRHSIntegrator( + const FunctionSpace& source_space, const FunctionSpace& target_space) + : OmegaHIntersectionRHSIntegrator( + std::dynamic_pointer_cast( + source_space.GetLayout()), + source_space.GetCoordinateSystem(), + std::dynamic_pointer_cast( + target_space.GetLayout()), + target_space.GetCoordinateSystem()) +{ +} + +OmegaHIntersectionRHSIntegrator::OmegaHIntersectionRHSIntegrator( + std::shared_ptr source_layout, + CoordinateSystem source_coordinate_system, + std::shared_ptr target_layout, + CoordinateSystem target_coordinate_system) + : OmegaHIntersectionRHSIntegrator( + BuildData(std::move(source_layout), source_coordinate_system, + std::move(target_layout), target_coordinate_system)) +{ +} + +OmegaHIntersectionRHSIntegrator::OmegaHIntersectionRHSIntegrator(Data data) + : point_set_(CoordinateSystem::Cartesian, std::move(data.coords)), + node_gids_(std::move(data.node_gids)), + coeffs_(std::move(data.coeffs)), + ndof_per_elem_(data.ndof_per_elem) +{ + const PetscInt nnz = static_cast(node_gids_.extent(0)); + PetscErrorCode ierr = + createSeqVec(PETSC_COMM_WORLD, data.num_target_dofs, &vec_); + CHKERRABORT(PETSC_COMM_WORLD, ierr); + ierr = VecSetPreallocationCOO(vec_, nnz, node_gids_.data()); + CHKERRABORT(PETSC_COMM_WORLD, ierr); +} + +OmegaHIntersectionRHSIntegrator::~OmegaHIntersectionRHSIntegrator() +{ + if (vec_) { + VecDestroy(&vec_); + } +} + +// --------------------------------------------------------------------------- +// OmegaHIntersectionRHSIntegrator public interface +// --------------------------------------------------------------------------- + +const IntegrationPointSet& +OmegaHIntersectionRHSIntegrator::GetIntegrationPoints() const noexcept +{ + return point_set_; +} + +Vec OmegaHIntersectionRHSIntegrator::GetVector() const noexcept +{ + return vec_; +} + +void OmegaHIntersectionRHSIntegrator::Assemble( + Rank2View sampled_values) +{ + const int ndof = ndof_per_elem_; + const std::size_t num_pts = + static_cast(node_gids_.extent(0)) / ndof; + PCMS_ALWAYS_ASSERT(static_cast(sampled_values.extent(0)) == + num_pts); + PCMS_ALWAYS_ASSERT(sampled_values.extent(1) >= 1); + + PetscErrorCode ierr = VecZeroEntries(vec_); + CHKERRABORT(PETSC_COMM_WORLD, ierr); + + auto sv = Kokkos::View>( + sampled_values.data_handle(), sampled_values.extent(0), + sampled_values.extent(1)); + Kokkos::View coo_vals("rhs_coo_vals", + num_pts * ndof); + auto coeffs = coeffs_; + Kokkos::parallel_for( + "rhs_coo_vals", static_cast(num_pts), KOKKOS_LAMBDA(int i) { + const PetscScalar f = static_cast(sv(i, 0)); + for (int j = 0; j < ndof; ++j) { + coo_vals(i * ndof + j) = + static_cast(coeffs(i * ndof + j)) * f; + } + }); + + ierr = VecSetValuesCOO(vec_, coo_vals.data(), ADD_VALUES); + CHKERRABORT(PETSC_COMM_WORLD, ierr); +} + +// --------------------------------------------------------------------------- +// Builder +// --------------------------------------------------------------------------- + +std::unique_ptr BuildOmegaHConservativeRHSIntegrator( + const FunctionSpace& source_space, const FunctionSpace& target_space) +{ + return std::make_unique(source_space, + target_space); +} + +} // namespace pcms diff --git a/src/pcms/transfer/omega_h_intersection_rhs_integrator.hpp b/src/pcms/transfer/omega_h_intersection_rhs_integrator.hpp new file mode 100644 index 00000000..45e36205 --- /dev/null +++ b/src/pcms/transfer/omega_h_intersection_rhs_integrator.hpp @@ -0,0 +1,95 @@ +#ifndef PCMS_TRANSFER_OMEGA_H_INTERSECTION_RHS_INTEGRATOR_HPP +#define PCMS_TRANSFER_OMEGA_H_INTERSECTION_RHS_INTEGRATOR_HPP + +#include "pcms/field/function_space.h" +#include "pcms/field/layout/omega_h_lagrange.h" +#include "pcms/transfer/integration_point_set.hpp" +#include "pcms/transfer/linear_form_integrator.hpp" +#include "pcms/transfer/mesh_intersection.hpp" +#include +#include + +namespace pcms +{ + +// Conservative-projection RHS integrator for Lagrange spaces on Omega_h 2D +// simplex meshes. Source and target orders are handled independently: the +// source enters only through pointwise samples (any supported order), while the +// target order selects the basis and DOF layout (P0: one element-constant DOF; +// P1: three vertex DOFs). The quadrature order is source_order + target_order. +// +// Construction validates both spaces, computes mesh intersections, and runs a +// two-pass device kernel (count then fill) to build the quadrature data and +// allocate the owned RHS vector sized to the target DOF count. +// +// Assemble(sampled_values): +// Zeros the owned vector, then for each integration point i and local target +// DOF j: +// contribution = stored_coeff[i][j] * sampled_values(i, 0) +// Scatter-adds via VecSetValuesCOO with ADD_VALUES. +// GetVector() returns the assembled vector (non-owning handle). +class OmegaHIntersectionRHSIntegrator : public LinearFormIntegrator +{ +public: + OmegaHIntersectionRHSIntegrator(const FunctionSpace& source_space, + const FunctionSpace& target_space); + OmegaHIntersectionRHSIntegrator( + std::shared_ptr source_layout, + CoordinateSystem source_coordinate_system, + std::shared_ptr target_layout, + CoordinateSystem target_coordinate_system); + ~OmegaHIntersectionRHSIntegrator(); + + Vec GetVector() const noexcept override; + + const IntegrationPointSet& GetIntegrationPoints() + const noexcept override; + + void Assemble( + Rank2View sampled_values) override; + +private: + // Internal data bundle produced by the two-pass device kernel. node_gids and + // coeffs are laid out with ndof_per_elem entries per integration point. + struct Data + { + Kokkos::View coords; // [num_pts, 2] + Kokkos::View node_gids; // [num_pts*ndof] + Kokkos::View coeffs; // [num_pts*ndof] + int ndof_per_elem = 0; // target DOFs/element + PetscInt num_target_dofs = 0; // target DOF count + }; + + static Data BuildData(const FunctionSpace& source_space, + const FunctionSpace& target_space); + static Data BuildData( + std::shared_ptr source_layout, + CoordinateSystem source_coordinate_system, + std::shared_ptr target_layout, + CoordinateSystem target_coordinate_system); + + // Order-templated two-pass builder, selected at runtime on the target order. + template + static Data BuildDataImpl(const OmegaHLagrangeLayout& source_layout, + const OmegaHLagrangeLayout& target_layout, + int quad_order); + + OmegaHIntersectionRHSIntegrator(Data data); + + Vec vec_ = nullptr; + IntegrationPointSet point_set_; + Kokkos::View + node_gids_; // COO idx, num_pts*ndof + Kokkos::View coeffs_; // weights, num_pts*ndof + int ndof_per_elem_ = 0; // target DOFs/element +}; + +// Builds an OmegaHIntersectionRHSIntegrator for conservative L2 projection. +// Both spaces must use OmegaHLagrangeLayout, scalar, Cartesian, order-0 or +// order-1 (source and target orders may differ). +std::unique_ptr BuildOmegaHConservativeRHSIntegrator( + const FunctionSpace& source_space, const FunctionSpace& target_space); + +} // namespace pcms + +#endif // PCMS_TRANSFER_OMEGA_H_INTERSECTION_RHS_INTEGRATOR_HPP diff --git a/src/pcms/transfer/omega_h_mass_integrator.cpp b/src/pcms/transfer/omega_h_mass_integrator.cpp new file mode 100644 index 00000000..e5fe259d --- /dev/null +++ b/src/pcms/transfer/omega_h_mass_integrator.cpp @@ -0,0 +1,138 @@ +#include "pcms/transfer/omega_h_mass_integrator.hpp" +#include "pcms/transfer/mass_matrix_integrator.hpp" +#include "pcms/transfer/omega_h_form_integrator_utils.hpp" +#include "pcms/transfer/petsc_utils.hpp" +#include "pcms/utility/assert.h" +#include "pcms/utility/memory_spaces.h" +#include +#include +#include +#include +#include +#include + +namespace pcms +{ + +OmegaHMassIntegrator::OmegaHMassIntegrator(const FunctionSpace& target_space) + : OmegaHMassIntegrator(std::dynamic_pointer_cast( + target_space.GetLayout()), + target_space.GetCoordinateSystem()) +{ +} + +OmegaHMassIntegrator::OmegaHMassIntegrator( + std::shared_ptr target_layout, + CoordinateSystem coordinate_system) +{ + detail::CheckOmegaHScalarLagrangeLayout(coordinate_system, target_layout, + "OmegaHMassIntegrator", "target"); + + Omega_h::Mesh& mesh = target_layout->GetMesh(); + const auto device_gids = target_layout->GetGids(); + const GO* gids_ptr = device_gids.data_handle(); + const PetscInt num_dofs = + static_cast(target_layout->GetNumGlobalDofHolder()); + const int nelems = mesh.nelems(); + + if (target_layout->GetOrder() == 0) { + // P0 target: piecewise-constant basis functions have disjoint support, so + // the mass matrix is diagonal with M_ee = area(e). One COO entry per + // element on its own (element-id) diagonal. + const auto& coords = mesh.coords(); + const auto& faces2nodes = mesh.ask_down(Omega_h::FACE, Omega_h::VERT).ab2b; + + const PetscInt nnz = static_cast(nelems); + Kokkos::View coo_rows("mass_coo_rows", nnz); + Kokkos::View coo_cols("mass_coo_cols", nnz); + Kokkos::View vals("mass_vals", nnz); + Kokkos::parallel_for( + "mass_p0_diag", nelems, KOKKOS_LAMBDA(int e) { + const auto verts = Omega_h::gather_verts<3>(faces2nodes, e); + const Omega_h::Matrix<2, 3> vm = + Omega_h::gather_vectors<3, 2>(coords, verts); + Omega_h::Few, 2> basis; + basis[0] = vm[1] - vm[0]; + basis[1] = vm[2] - vm[0]; + const Omega_h::Real area = + Kokkos::fabs(Omega_h::triangle_area_from_basis(basis)); + const PetscInt g = static_cast(gids_ptr[e]); + coo_rows(e) = g; + coo_cols(e) = g; + vals(e) = static_cast(area); + }); + + PetscErrorCode ierr = + createSeqAIJMat(PETSC_COMM_WORLD, num_dofs, num_dofs, 0, nullptr, &mat_); + CHKERRABORT(PETSC_COMM_WORLD, ierr); + ierr = MatSetPreallocationCOO(mat_, nnz, coo_rows.data(), coo_cols.data()); + CHKERRABORT(PETSC_COMM_WORLD, ierr); + ierr = MatSetValuesCOO(mat_, vals.data(), INSERT_VALUES); + CHKERRABORT(PETSC_COMM_WORLD, ierr); + return; + } + + // P1 target: consistent mass matrix assembled from MeshField per-element 3x3 + // blocks. (Higher MeshField orders extend this branch via + // getTriangleElement.) + MeshField::OmegahMeshField + omf(mesh); + auto coordField = omf.getCoordField(); + const auto [shp, map] = MeshField::Omegah::getTriangleElement<1>(mesh); + MeshField::FieldElement coordFe(mesh.nelems(), coordField, shp, map); + auto elm_mass_dev = buildElementMassMatrix(mesh, coordFe); + + // Build COO sparsity pattern on device: each element contributes a 3x3 block. + const auto& faces2nodes = mesh.ask_down(Omega_h::FACE, Omega_h::VERT).ab2b; + + const PetscInt nnz = static_cast(nelems) * 9; + Kokkos::View coo_rows("mass_coo_rows", nnz); + Kokkos::View coo_cols("mass_coo_cols", nnz); + Kokkos::parallel_for( + "mass_coo_pattern", nelems, KOKKOS_LAMBDA(int e) { + const auto verts = Omega_h::gather_verts<3>(faces2nodes, e); + 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(gids_ptr[verts[i]]); + coo_cols(idx) = static_cast(gids_ptr[verts[j]]); + } + } + }); + + // 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 + // it can be passed directly to MatSetValuesCOO — no host copy needed. + PetscErrorCode ierr = + createSeqAIJMat(PETSC_COMM_WORLD, num_dofs, num_dofs, 0, nullptr, &mat_); + CHKERRABORT(PETSC_COMM_WORLD, ierr); + ierr = MatSetPreallocationCOO(mat_, nnz, coo_rows.data(), coo_cols.data()); + CHKERRABORT(PETSC_COMM_WORLD, ierr); + ierr = MatSetValuesCOO(mat_, elm_mass_dev.data(), INSERT_VALUES); + CHKERRABORT(PETSC_COMM_WORLD, ierr); +} + +OmegaHMassIntegrator::~OmegaHMassIntegrator() +{ + if (mat_) { + MatDestroy(&mat_); + } +} + +Mat OmegaHMassIntegrator::GetMatrix() const noexcept +{ + return mat_; +} + +// --------------------------------------------------------------------------- +// Builder +// --------------------------------------------------------------------------- + +std::unique_ptr BuildOmegaHMassIntegrator( + const FunctionSpace& target_space) +{ + return std::make_unique(target_space); +} + +} // namespace pcms diff --git a/src/pcms/transfer/omega_h_mass_integrator.hpp b/src/pcms/transfer/omega_h_mass_integrator.hpp new file mode 100644 index 00000000..54cc7a66 --- /dev/null +++ b/src/pcms/transfer/omega_h_mass_integrator.hpp @@ -0,0 +1,45 @@ +#ifndef PCMS_TRANSFER_OMEGA_H_MASS_INTEGRATOR_HPP +#define PCMS_TRANSFER_OMEGA_H_MASS_INTEGRATOR_HPP + +#include "pcms/field/function_space.h" +#include "pcms/field/layout/omega_h_lagrange.h" +#include "pcms/transfer/bilinear_form_integrator.hpp" +#include +#include + +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. +// +// 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 +// utility is needed. The integrator owns the matrix; callers receive a +// non-owning handle via GetMatrix(). PETSc reference-counts the matrix, so a +// KSP that calls KSPSetOperators(GetMatrix()) safely outlives this object. +class OmegaHMassIntegrator : public BilinearFormIntegrator +{ +public: + explicit OmegaHMassIntegrator(const FunctionSpace& target_space); + OmegaHMassIntegrator( + std::shared_ptr target_layout, + CoordinateSystem coordinate_system); + ~OmegaHMassIntegrator(); + + Mat GetMatrix() const noexcept override; + +private: + Mat mat_ = nullptr; +}; + +// Builds an OmegaHMassIntegrator for the given target space. +// target_space must use OmegaHLagrangeLayout, scalar, Cartesian, order-0 or +// order-1. +std::unique_ptr BuildOmegaHMassIntegrator( + const FunctionSpace& target_space); + +} // namespace pcms + +#endif // PCMS_TRANSFER_OMEGA_H_MASS_INTEGRATOR_HPP diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 1c603c45..914ab322 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -404,11 +404,14 @@ if(Catch2_FOUND) if(PCMS_ENABLE_MESHFIELDS) list(APPEND PCMS_UNIT_TEST_SOURCES - test_load_vector.cpp) + test_omega_h_form_integrator_utils.cpp) endif() if(PCMS_ENABLE_PETSC AND PCMS_ENABLE_MESHFIELDS) list(APPEND PCMS_UNIT_TEST_SOURCES + test_omega_h_intersection_rhs_integrator.cpp + test_omega_h_mass_integrator.cpp + test_omega_h_mc_rhs_integrator.cpp test_mesh_intersection_field_transfer.cpp) endif() diff --git a/test/test_mesh_intersection_field_transfer.cpp b/test/test_mesh_intersection_field_transfer.cpp index bd34af54..a66579e7 100644 --- a/test/test_mesh_intersection_field_transfer.cpp +++ b/test/test_mesh_intersection_field_transfer.cpp @@ -5,223 +5,329 @@ #include #include -#include -#include #include #include #include "field_test_utils.h" +#include +#include + namespace { -double integrate_linear_field(Omega_h::Mesh& mesh, const Omega_h::Reals& u) +// Builds a unit-square 2D simplex mesh from the given element connectivity and +// adds the geometric classification tags required to build an Omega_h-backed +// Lagrange function space. +Omega_h::Mesh BuildUnitSquare(Omega_h::Library& lib, const Omega_h::LOs& ev2v) +{ + const Omega_h::Reals coords({ + 0.0, 0.0, // v0 + 1.0, 0.0, // v1 + 1.0, 1.0, // v2 + 0.0, 1.0 // v3 + }); + Omega_h::Mesh mesh(&lib); + Omega_h::build_from_elems_and_coords(&mesh, OMEGA_H_SIMPLEX, 2, ev2v, coords); + for (Omega_h::Int dim = 0; dim <= 2; ++dim) { + mesh.add_tag( + dim, "class_dim", 1, + Omega_h::Read(mesh.nents(dim), Omega_h::I8(dim))); + mesh.add_tag( + dim, "class_id", 1, + Omega_h::Read(mesh.nents(dim), Omega_h::ClassId(0))); + } + return mesh; +} + +pcms::LagrangeFunctionSpace MakeP1Space(Omega_h::Mesh& mesh) +{ + return pcms::LagrangeFunctionSpace::FromMesh( + mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", + pcms::LagrangeFunctionSpace::Backend::OmegaH); +} + +pcms::LagrangeFunctionSpace MakeP0Space(Omega_h::Mesh& mesh) +{ + return pcms::LagrangeFunctionSpace::FromMesh( + mesh, 0, 1, pcms::CoordinateSystem::Cartesian, "global", + pcms::LagrangeFunctionSpace::Backend::OmegaH); +} + +// Integral over the mesh of an order-0 (piecewise-constant) field whose values +// are stored in element order. Exact by construction. +double IntegrateP0Field(Omega_h::Mesh& mesh, + const pcms::Field& field) { + const auto values = pcms::FlattenToRank1View(field.GetDOFHolderDataHost()); const auto elem_areas = Omega_h::measure_elements_real(&mesh); - const auto elem_verts = mesh.ask_elem_verts(); - const auto nverts = mesh.nverts(); + const auto elem_areas_h = Omega_h::HostRead(elem_areas); + + double integral = 0.0; + for (Omega_h::LO e = 0; e < mesh.nelems(); ++e) { + integral += elem_areas_h[e] * values[e]; + } + return integral; +} - REQUIRE(static_cast(u.size()) == nverts); +// Per-element centroid coordinates in element order. +std::vector> ElementCentroids(Omega_h::Mesh& mesh) +{ + const auto coords_h = Omega_h::HostRead(mesh.coords()); + const auto elem_verts_h = + Omega_h::HostRead(mesh.ask_elem_verts()); + std::vector> centroids(mesh.nelems()); + for (Omega_h::LO e = 0; e < mesh.nelems(); ++e) { + double cx = 0.0, cy = 0.0; + for (int k = 0; k < 3; ++k) { + const Omega_h::LO v = elem_verts_h[3 * e + k]; + cx += coords_h[2 * v + 0]; + cy += coords_h[2 * v + 1]; + } + centroids[e] = {cx / 3.0, cy / 3.0}; + } + return centroids; +} +// Integral over the mesh of an order-1 nodal field whose values are stored in +// vertex order (as produced by Field::GetDOFHolderDataHost for the +// Omega_h backend), computed exactly via per-element vertex averaging. +double IntegrateP1Field(Omega_h::Mesh& mesh, + const pcms::Field& field) +{ + const auto values = pcms::FlattenToRank1View(field.GetDOFHolderDataHost()); + const auto elem_areas = Omega_h::measure_elements_real(&mesh); + const auto elem_verts = mesh.ask_elem_verts(); const auto elem_areas_h = Omega_h::HostRead(elem_areas); const auto elem_verts_h = Omega_h::HostRead(elem_verts); - const auto u_h = Omega_h::HostRead(u); double integral = 0.0; for (Omega_h::LO e = 0; e < mesh.nelems(); ++e) { const Omega_h::LO v0 = elem_verts_h[3 * e + 0]; const Omega_h::LO v1 = elem_verts_h[3 * e + 1]; const Omega_h::LO v2 = elem_verts_h[3 * e + 2]; - - const double area = elem_areas_h[e]; - const double avg = (u_h[v0] + u_h[v1] + u_h[v2]) / 3.0; - integral += area * avg; + const double avg = (values[v0] + values[v1] + values[v2]) / 3.0; + integral += elem_areas_h[e] * avg; } return integral; } -Omega_h::Reals make_omega_h_reals( - pcms::Rank1View values) -{ - Omega_h::HostWrite values_h(values.size()); - for (size_t i = 0; i < values.size(); ++i) { - values_h[i] = values[i]; - } - return Omega_h::Reals(values_h); -} - } // namespace -TEST_CASE("mesh intersection linear/constant conservation", +// Fields that already live in the target order-1 space (constants and affine +// functions) must be reproduced exactly by the L2 projection, and the integral +// must be conserved. These properties hold without any reference solver, so +// they pin down correctness on their own. +TEST_CASE("OmegaHConservativeProjection reproduces constant and linear fields", "[transfer][mesh_intersection]") { Omega_h::Library lib; - Omega_h::Reals coords({ - 0.0, 0.0, // v0 - 1.0, 0.0, // v1 - 1.0, 1.0, // v2 - 0.0, 1.0 // v3 - }); - - // Source mesh: diagonal (v0-v2) - Omega_h::LOs ev2v_source({0, 1, 2, 0, 2, 3}); - Omega_h::Mesh source_mesh(&lib); - Omega_h::build_from_elems_and_coords(&source_mesh, OMEGA_H_SIMPLEX, 2, - ev2v_source, coords); + // Source and target triangulate the same square along opposite diagonals. + Omega_h::Mesh source_mesh = + BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3})); + Omega_h::Mesh target_mesh = + BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3})); - // Target mesh: opposite diagonal (v1-v3) - Omega_h::LOs ev2v_target({0, 1, 3, 1, 2, 3}); - Omega_h::Mesh target_mesh(&lib); - Omega_h::build_from_elems_and_coords(&target_mesh, OMEGA_H_SIMPLEX, 2, - ev2v_target, coords); + auto source_space = MakeP1Space(source_mesh); + auto target_space = MakeP1Space(target_mesh); - const auto src_coords = source_mesh.coords(); + auto source = source_space.CreateField(); + auto target = target_space.CreateField(); - auto intersections = pcms::intersectTargets(source_mesh, target_mesh); + pcms::OmegaHConservativeProjection projection(source_space, target_space); SECTION("constant field is preserved and conserved") { const double c = 2.0; - Omega_h::Write source_const(source_mesh.nverts()); - Omega_h::parallel_for( - source_const.size(), OMEGA_H_LAMBDA(int i) { source_const[i] = c; }); + pcms::test::SetField( + source, OMEGA_H_LAMBDA(pcms::Real, pcms::Real) { return c; }); - auto projected = pcms::solveGalerkinProjection(target_mesh, source_mesh, - intersections, source_const); - auto projected_h = Omega_h::HostRead(projected); + projection.Apply(source, target); - REQUIRE(static_cast(projected.size()) == target_mesh.nverts()); + const auto target_values = + pcms::FlattenToRank1View(target.GetDOFHolderDataHost()); + REQUIRE(static_cast(target_values.size()) == + target_mesh.nverts()); for (Omega_h::LO i = 0; i < target_mesh.nverts(); ++i) { - REQUIRE(projected_h[i] == Catch::Approx(c).margin(1e-10)); + REQUIRE(target_values[i] == Catch::Approx(c).margin(1e-10)); } - - const double source_integral = - integrate_linear_field(source_mesh, source_const); - const double target_integral = - integrate_linear_field(target_mesh, projected); - REQUIRE(target_integral == Catch::Approx(source_integral).margin(1e-10)); + REQUIRE(IntegrateP1Field(target_mesh, target) == + Catch::Approx(IntegrateP1Field(source_mesh, source)).margin(1e-10)); } - SECTION("linear field is reproduced on target vertices") + SECTION("linear field is reproduced on target vertices and conserved") { - Omega_h::Write source_linear(source_mesh.nverts()); - Omega_h::parallel_for( - source_mesh.nverts(), OMEGA_H_LAMBDA(int i) { - const double x = src_coords[2 * i + 0]; - const double y = src_coords[2 * i + 1]; - source_linear[i] = x + y; - }); - - auto projected = pcms::solveGalerkinProjection( - target_mesh, source_mesh, intersections, source_linear); - auto projected_h = Omega_h::HostRead(projected); + pcms::test::SetField( + source, OMEGA_H_LAMBDA(pcms::Real x, pcms::Real y) { return x + y; }); + + projection.Apply(source, target); + + const auto target_values = + pcms::FlattenToRank1View(target.GetDOFHolderDataHost()); const auto tgt_coords_h = Omega_h::HostRead(target_mesh.coords()); - for (Omega_h::LO i = 0; i < target_mesh.nverts(); ++i) { const double expected = tgt_coords_h[2 * i + 0] + tgt_coords_h[2 * i + 1]; - REQUIRE(projected_h[i] == Catch::Approx(expected).margin(1e-9)); + REQUIRE(target_values[i] == Catch::Approx(expected).margin(1e-9)); } + REQUIRE(IntegrateP1Field(target_mesh, target) == + Catch::Approx(IntegrateP1Field(source_mesh, source)).margin(1e-9)); } } -TEST_CASE("OmegaHConservativeProjection matches conservative projection solver", +// For a field that does not live in the target space the projection is not an +// interpolation, but conservation of the integral is the defining property of +// the conservative (Galerkin) projection and must still hold exactly. This +// also exercises a second Apply with a different source field to confirm the +// cached integrators/evaluator/factorization are reused without stale state. +TEST_CASE("OmegaHConservativeProjection conserves the integral", "[transfer][mesh_intersection]") { Omega_h::Library lib; - Omega_h::Reals coords({0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 1.0}); - - Omega_h::LOs ev2v_source({0, 1, 2, 0, 2, 3}); - Omega_h::Mesh source_mesh(&lib); - Omega_h::build_from_elems_and_coords(&source_mesh, OMEGA_H_SIMPLEX, 2, - ev2v_source, coords); - - Omega_h::LOs ev2v_target({0, 1, 3, 1, 2, 3}); - Omega_h::Mesh target_mesh(&lib); - Omega_h::build_from_elems_and_coords(&target_mesh, OMEGA_H_SIMPLEX, 2, - ev2v_target, coords); - - // Add classification for all dimensions - // Vertices (dim 0) - source_mesh.add_tag( - 0, "class_dim", 1, - Omega_h::Read(source_mesh.nverts(), Omega_h::I8(0))); - source_mesh.add_tag( - 0, "class_id", 1, - Omega_h::Read(source_mesh.nverts(), Omega_h::ClassId(0))); - - // Edges (dim 1) - source_mesh.add_tag( - 1, "class_dim", 1, - Omega_h::Read(source_mesh.nedges(), Omega_h::I8(1))); - source_mesh.add_tag( - 1, "class_id", 1, - Omega_h::Read(source_mesh.nedges(), Omega_h::ClassId(0))); - - // Faces (dim 2) - source_mesh.add_tag( - 2, "class_dim", 1, - Omega_h::Read(source_mesh.nelems(), Omega_h::I8(2))); - source_mesh.add_tag( - 2, "class_id", 1, - Omega_h::Read(source_mesh.nelems(), Omega_h::ClassId(0))); - - // Add classification for all dimensions - // Vertices (dim 0) - target_mesh.add_tag( - 0, "class_dim", 1, - Omega_h::Read(target_mesh.nverts(), Omega_h::I8(0))); - target_mesh.add_tag( - 0, "class_id", 1, - Omega_h::Read(target_mesh.nverts(), Omega_h::ClassId(0))); - - // Edges (dim 1) - target_mesh.add_tag( - 1, "class_dim", 1, - Omega_h::Read(target_mesh.nedges(), Omega_h::I8(1))); - target_mesh.add_tag( - 1, "class_id", 1, - Omega_h::Read(target_mesh.nedges(), Omega_h::ClassId(0))); - - // Faces (dim 2) - target_mesh.add_tag( - 2, "class_dim", 1, - Omega_h::Read(target_mesh.nelems(), Omega_h::I8(2))); - target_mesh.add_tag( - 2, "class_id", 1, - Omega_h::Read(target_mesh.nelems(), Omega_h::ClassId(0))); - - auto source_space = pcms::LagrangeFunctionSpace::FromMesh( - source_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", - pcms::LagrangeFunctionSpace::Backend::OmegaH); - auto target_space = pcms::LagrangeFunctionSpace::FromMesh( - target_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", - pcms::LagrangeFunctionSpace::Backend::OmegaH); + Omega_h::Mesh source_mesh = + BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3})); + Omega_h::Mesh target_mesh = + BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3})); + + auto source_space = MakeP1Space(source_mesh); + auto target_space = MakeP1Space(target_mesh); auto source = source_space.CreateField(); auto target = target_space.CreateField(); + pcms::OmegaHConservativeProjection projection(source_space, target_space); + pcms::test::SetField( source, OMEGA_H_LAMBDA(pcms::Real x, pcms::Real y) { return x * x + x * y + 0.5 * y * y; }); + projection.Apply(source, target); + REQUIRE(IntegrateP1Field(target_mesh, target) == + Catch::Approx(IntegrateP1Field(source_mesh, source)).margin(1e-12)); + + // Second Apply with a different source field — verifies cached state is + // reused correctly and is not stale. + pcms::test::SetField( + source, + OMEGA_H_LAMBDA(pcms::Real x, pcms::Real y) { return 2.0 * x - y + 0.5; }); + projection.Apply(source, target); + REQUIRE(IntegrateP1Field(target_mesh, target) == + Catch::Approx(IntegrateP1Field(source_mesh, source)).margin(1e-12)); +} + +// P0 (piecewise-constant) source -> P1 target. A constant lives in P1, so it is +// reproduced exactly; a non-constant P0 source is not reproduced but its +// integral must still be conserved. +TEST_CASE("OmegaHConservativeProjection P0 source to P1 target", + "[transfer][mesh_intersection]") +{ + Omega_h::Library lib; + + Omega_h::Mesh source_mesh = + BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3})); + Omega_h::Mesh target_mesh = + BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3})); - const auto intersections = pcms::intersectTargets(source_mesh, target_mesh); - const auto expected = - pcms::solveGalerkinProjection(target_mesh, source_mesh, intersections, - make_omega_h_reals(pcms::FlattenToRank1View( - source.GetDOFHolderDataHost()))); - const auto expected_h = Omega_h::HostRead(expected); + auto source_space = MakeP0Space(source_mesh); + auto target_space = MakeP1Space(target_mesh); + + auto source = source_space.CreateField(); + auto target = target_space.CreateField(); pcms::OmegaHConservativeProjection projection(source_space, target_space); - projection.Apply(source, target); - const auto target_values = - pcms::FlattenToRank1View(target.GetDOFHolderDataHost()); - REQUIRE(static_cast(target_values.size()) == - target_mesh.nverts()); - for (Omega_h::LO i = 0; i < target_mesh.nverts(); ++i) { - REQUIRE(target_values[i] == Catch::Approx(expected_h[i]).margin(1e-12)); + SECTION("constant field is preserved and conserved") + { + const double c = 3.5; + pcms::test::SetField( + source, OMEGA_H_LAMBDA(pcms::Real, pcms::Real) { return c; }); + + projection.Apply(source, target); + + const auto target_values = + pcms::FlattenToRank1View(target.GetDOFHolderDataHost()); + REQUIRE(static_cast(target_values.size()) == + target_mesh.nverts()); + for (Omega_h::LO i = 0; i < target_mesh.nverts(); ++i) { + REQUIRE(target_values[i] == Catch::Approx(c).margin(1e-10)); + } + REQUIRE(IntegrateP1Field(target_mesh, target) == + Catch::Approx(IntegrateP0Field(source_mesh, source)).margin(1e-10)); + } + + SECTION("non-constant P0 source conserves the integral") + { + pcms::test::SetField( + source, + OMEGA_H_LAMBDA(pcms::Real x, pcms::Real y) { return x + 2.0 * y; }); + + projection.Apply(source, target); + + REQUIRE(IntegrateP1Field(target_mesh, target) == + Catch::Approx(IntegrateP0Field(source_mesh, source)).margin(1e-10)); + } +} + +// P1 source -> P0 (piecewise-constant) target. A constant is reproduced +// exactly; a linear source projects to its cell average, which for a linear +// function equals the value at each element centroid. The integral is conserved +// in every case. +TEST_CASE("OmegaHConservativeProjection P1 source to P0 target", + "[transfer][mesh_intersection]") +{ + Omega_h::Library lib; + + Omega_h::Mesh source_mesh = + BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3})); + Omega_h::Mesh target_mesh = + BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3})); + + auto source_space = MakeP1Space(source_mesh); + auto target_space = MakeP0Space(target_mesh); + + auto source = source_space.CreateField(); + auto target = target_space.CreateField(); + + pcms::OmegaHConservativeProjection projection(source_space, target_space); + + SECTION("constant field is preserved and conserved") + { + const double c = -1.25; + pcms::test::SetField( + source, OMEGA_H_LAMBDA(pcms::Real, pcms::Real) { return c; }); + + projection.Apply(source, target); + + const auto target_values = + pcms::FlattenToRank1View(target.GetDOFHolderDataHost()); + REQUIRE(static_cast(target_values.size()) == + target_mesh.nelems()); + for (Omega_h::LO e = 0; e < target_mesh.nelems(); ++e) { + REQUIRE(target_values[e] == Catch::Approx(c).margin(1e-10)); + } + REQUIRE(IntegrateP0Field(target_mesh, target) == + Catch::Approx(IntegrateP1Field(source_mesh, source)).margin(1e-10)); + } + + SECTION("linear field projects to cell average and conserves the integral") + { + const auto f = [](double x, double y) { return x + 2.0 * y; }; + pcms::test::SetField( + source, + OMEGA_H_LAMBDA(pcms::Real x, pcms::Real y) { return x + 2.0 * y; }); + + projection.Apply(source, target); + + const auto target_values = + pcms::FlattenToRank1View(target.GetDOFHolderDataHost()); + const auto centroids = ElementCentroids(target_mesh); + REQUIRE(static_cast(target_values.size()) == centroids.size()); + for (std::size_t e = 0; e < centroids.size(); ++e) { + const double expected = f(centroids[e][0], centroids[e][1]); + REQUIRE(target_values[e] == Catch::Approx(expected).margin(1e-9)); + } + REQUIRE(IntegrateP0Field(target_mesh, target) == + Catch::Approx(IntegrateP1Field(source_mesh, source)).margin(1e-9)); } } diff --git a/test/test_omega_h_form_integrator_utils.cpp b/test/test_omega_h_form_integrator_utils.cpp new file mode 100644 index 00000000..91cedac7 --- /dev/null +++ b/test/test_omega_h_form_integrator_utils.cpp @@ -0,0 +1,80 @@ +#include +#include +#include + +namespace +{ + +r3d::Polytope<2> MakePoly(int nverts, + std::initializer_list> pts) +{ + r3d::Polytope<2> poly; + poly.nverts = nverts; + int i = 0; + for (auto [x, y] : pts) { + poly.verts[i].pos[0] = x; + poly.verts[i].pos[1] = y; + poly.verts[i].pnbrs[0] = -1; + poly.verts[i].pnbrs[1] = -1; + ++i; + } + return poly; +} + +} // namespace + +TEST_CASE( + "RemoveDuplicateVerticesAndFixLinks: no duplicates leaves triangle intact", + "[form_integrator_utils]") +{ + // A clean CCW triangle — nothing to remove, area should equal 0.5. + auto poly = MakePoly(3, {{0.0, 0.0}, {1.0, 0.0}, {0.0, 1.0}}); + const int n = pcms::detail::RemoveDuplicateVerticesAndFixLinks(poly); + REQUIRE(n == 3); + REQUIRE(poly.nverts == 3); + // The function sets correct pnbrs, so r3d::measure is valid after the call. + const double area = r3d::measure(poly); + REQUIRE(area == Catch::Approx(0.5).epsilon(1e-12)); +} + +TEST_CASE("RemoveDuplicateVerticesAndFixLinks: one duplicate vertex is removed", + "[form_integrator_utils]") +{ + // 4-vertex polygon; last vertex duplicates the first within the default tol. + // After dedup: 3 unique vertices, area of the underlying triangle preserved. + auto poly = MakePoly(4, {{0.0, 0.0}, {1.0, 0.0}, {0.0, 1.0}, {5e-13, 0.0}}); + const int n = pcms::detail::RemoveDuplicateVerticesAndFixLinks(poly); + REQUIRE(n == 3); + REQUIRE(poly.nverts == 3); + const double area = r3d::measure(poly); + REQUIRE(area == Catch::Approx(0.5).epsilon(1e-10)); +} + +TEST_CASE( + "RemoveDuplicateVerticesAndFixLinks: polygon becomes degenerate after dedup", + "[form_integrator_utils]") +{ + // 3 vertices where two positions are identical — only 2 unique remain. + // The function must mark the polygon degenerate: all pnbrs set to -1. + auto poly = MakePoly(3, {{0.0, 0.0}, {1.0, 0.0}, {0.0, 0.0}}); + const int n = pcms::detail::RemoveDuplicateVerticesAndFixLinks(poly); + REQUIRE(n == 2); + REQUIRE(poly.nverts == 2); + for (int i = 0; i < 2; ++i) { + CAPTURE(i); + CHECK(poly.verts[i].pnbrs[0] == -1); + CHECK(poly.verts[i].pnbrs[1] == -1); + } +} + +TEST_CASE( + "RemoveDuplicateVerticesAndFixLinks: near-tolerance vertex is retained", + "[form_integrator_utils]") +{ + // Vertex at (0, 2e-12) differs from (0, 0) by 2e-12 in y, which exceeds the + // default tol=1e-12, so it must NOT be treated as a duplicate. + auto poly = MakePoly(3, {{0.0, 0.0}, {1.0, 0.0}, {0.0, 2e-12}}); + const int n = pcms::detail::RemoveDuplicateVerticesAndFixLinks(poly); + REQUIRE(n == 3); + REQUIRE(poly.nverts == 3); +} diff --git a/test/test_omega_h_intersection_rhs_integrator.cpp b/test/test_omega_h_intersection_rhs_integrator.cpp new file mode 100644 index 00000000..b778f2ae --- /dev/null +++ b/test/test_omega_h_intersection_rhs_integrator.cpp @@ -0,0 +1,337 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "field_test_utils.h" +#include +#include +#include +#include + +// --------------------------------------------------------------------------- +// Helpers +// --------------------------------------------------------------------------- + +namespace +{ + +// Build a unit-square 2D simplex mesh. +// diagonal=0: T0=(0,1,3), T1=(1,2,3) +// diagonal=1: T0=(0,1,2), T1=(0,2,3) +Omega_h::Mesh BuildUnitSquare(Omega_h::Library& lib, int diagonal) +{ + const Omega_h::Reals coords({0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 1.0}); + Omega_h::LOs ev2v = (diagonal == 0) ? Omega_h::LOs({0, 1, 3, 1, 2, 3}) + : Omega_h::LOs({0, 1, 2, 0, 2, 3}); + Omega_h::Mesh mesh(&lib); + Omega_h::build_from_elems_and_coords(&mesh, OMEGA_H_SIMPLEX, 2, ev2v, coords); + mesh.add_tag( + 0, "class_dim", 1, + Omega_h::Read(mesh.nverts(), Omega_h::I8(0))); + mesh.add_tag( + 0, "class_id", 1, + Omega_h::Read(mesh.nverts(), Omega_h::ClassId(0))); + mesh.add_tag( + 1, "class_dim", 1, + Omega_h::Read(mesh.nedges(), Omega_h::I8(1))); + mesh.add_tag( + 1, "class_id", 1, + Omega_h::Read(mesh.nedges(), Omega_h::ClassId(0))); + mesh.add_tag( + 2, "class_dim", 1, + Omega_h::Read(mesh.nelems(), Omega_h::I8(2))); + mesh.add_tag( + 2, "class_id", 1, + Omega_h::Read(mesh.nelems(), Omega_h::ClassId(0))); + return mesh; +} + +// Independent reference for the assembled conservative load vector +// b_j = \int phi_j^target f dx +// when f is exactly representable in the target order-1 space (constant or +// affine): the integral then equals (M g)_j, where M is the target consistent +// mass matrix and g_k = f(target node k). Assembled here from per-element mass +// matrices on the target mesh alone — no intersection or source-mesh machinery +// — so it cross-checks the integrator instead of restating its computation. +// +// The result is keyed by global DOF id to match the PETSc row indexing of the +// integrator's assembled vector. +std::unordered_map ExpectedLoadByGid( + Omega_h::Mesh& target_mesh, const pcms::OmegaHLagrangeLayout& target_layout, + const std::vector& g_by_vert) +{ + const auto coords_h = Omega_h::HostRead(target_mesh.coords()); + const auto faces2nodes_h = Omega_h::HostRead( + target_mesh.ask_down(Omega_h::FACE, Omega_h::VERT).ab2b); + const auto& gids = target_layout.GetGidsHost(); + + std::unordered_map b; + for (int e = 0; e < target_mesh.nelems(); ++e) { + const int v[3] = {faces2nodes_h[3 * e + 0], faces2nodes_h[3 * e + 1], + faces2nodes_h[3 * e + 2]}; + const double x0 = coords_h[2 * v[0]], y0 = coords_h[2 * v[0] + 1]; + const double x1 = coords_h[2 * v[1]], y1 = coords_h[2 * v[1] + 1]; + const double x2 = coords_h[2 * v[2]], y2 = coords_h[2 * v[2] + 1]; + const double area = + 0.5 * std::abs((x1 - x0) * (y2 - y0) - (x2 - x0) * (y1 - y0)); + + // Consistent P1 element mass matrix M^e = (area/12) [[2,1,1],[1,2,1], + // [1,1,2]], so (M^e g)_j = (area/12) (g_j + sum_k g_k). + const double g_sum = g_by_vert[v[0]] + g_by_vert[v[1]] + g_by_vert[v[2]]; + for (int j = 0; j < 3; ++j) { + const double bj = (area / 12.0) * (g_by_vert[v[j]] + g_sum); + b[static_cast(gids[v[j]])] += static_cast(bj); + } + } + return b; +} + +// Assembles the integrator's load vector for source_field and compares it, +// per global DOF id, against the supplied expected values. +void CheckAssembledLoadMatches( + pcms::FunctionSpace& source_space, pcms::LinearFormIntegrator& integrator, + const pcms::Field& source_field, + const std::unordered_map& expected) +{ + const auto& pts = integrator.GetIntegrationPoints(); + const std::size_t npts = pts.GetCoordinates().GetValues().extent(0); + + auto evaluator = source_space.CreatePointEvaluator( + pcms::EvaluationRequest::FromCoordinates(pts.GetCoordinates())); + Kokkos::View sampled("sampled", npts, + 1); + evaluator->Evaluate(source_field, pcms::MakeRank2View(sampled)); + + integrator.Assemble(pcms::MakeConstRank2View(sampled)); + Vec vec = integrator.GetVector(); + + const PetscScalar* vec_array = nullptr; + VecGetArrayRead(vec, &vec_array); + for (const auto& [gid, ref_val] : expected) { + const pcms::Real got = static_cast(vec_array[gid]); + CAPTURE(gid, ref_val, got); + CHECK(got == Catch::Approx(ref_val).epsilon(1e-10)); + } + VecRestoreArrayRead(vec, &vec_array); + // vec is owned by integrator; no VecDestroy here. +} + +} // namespace + +// --------------------------------------------------------------------------- +// Tests +// --------------------------------------------------------------------------- + +TEST_CASE("OmegaHIntersectionRHSIntegrator: integration points lie inside " + "target domain", + "[rhs_integrator]") +{ + Omega_h::Library lib; + auto source_mesh = BuildUnitSquare(lib, 1); + auto target_mesh = BuildUnitSquare(lib, 0); + + auto source_space = pcms::LagrangeFunctionSpace::FromMesh( + source_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", + pcms::LagrangeFunctionSpace::Backend::OmegaH); + auto target_space = pcms::LagrangeFunctionSpace::FromMesh( + target_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", + pcms::LagrangeFunctionSpace::Backend::OmegaH); + + auto integrator = + pcms::BuildOmegaHConservativeRHSIntegrator(source_space, target_space); + + const auto raw_coords = + integrator->GetIntegrationPoints().GetCoordinates().GetValues(); + auto raw_coords_view = Kokkos::View>( + raw_coords.data_handle(), raw_coords.extent(0), raw_coords.extent(1)); + auto raw_coords_host = + Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, raw_coords_view); + const std::size_t n = raw_coords_host.extent(0); + + REQUIRE(n > 0); + for (std::size_t i = 0; i < n; ++i) { + CAPTURE(i, raw_coords_host(i, 0), raw_coords_host(i, 1)); + CHECK(raw_coords_host(i, 0) >= -1e-12); + CHECK(raw_coords_host(i, 0) <= 1.0 + 1e-12); + CHECK(raw_coords_host(i, 1) >= -1e-12); + CHECK(raw_coords_host(i, 1) <= 1.0 + 1e-12); + } +} + +TEST_CASE("OmegaHIntersectionRHSIntegrator: zero source field gives zero " + "assembled vector", + "[rhs_integrator]") +{ + Omega_h::Library lib; + auto source_mesh = BuildUnitSquare(lib, 1); + auto target_mesh = BuildUnitSquare(lib, 0); + + auto source_space = pcms::LagrangeFunctionSpace::FromMesh( + source_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", + pcms::LagrangeFunctionSpace::Backend::OmegaH); + auto target_space = pcms::LagrangeFunctionSpace::FromMesh( + target_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", + pcms::LagrangeFunctionSpace::Backend::OmegaH); + + auto source_field = source_space.CreateField(); + // Default-constructed field has zero data. + + auto integrator = + pcms::BuildOmegaHConservativeRHSIntegrator(source_space, target_space); + const auto& pts = integrator->GetIntegrationPoints(); + const std::size_t npts = pts.GetCoordinates().GetValues().extent(0); + + auto evaluator = source_space.CreatePointEvaluator( + pcms::EvaluationRequest::FromCoordinates(pts.GetCoordinates())); + + Kokkos::View sampled("sampled", npts, + 1); + evaluator->Evaluate(source_field, pcms::MakeRank2View(sampled)); + + integrator->Assemble(pcms::MakeConstRank2View(sampled)); + Vec vec = integrator->GetVector(); + + PetscReal norm = 0.0; + VecNorm(vec, NORM_INFINITY, &norm); + INFO("vector inf-norm=" << norm); + CHECK(static_cast(norm) == Catch::Approx(0.0).margin(1e-14)); + // vec is owned by integrator; no VecDestroy here. +} + +TEST_CASE("OmegaHIntersectionRHSIntegrator: constant field matches " + "consistent-mass load", + "[rhs_integrator]") +{ + // For a constant source field the assembled load vector must equal the + // target consistent mass matrix applied to the constant nodal vector. + Omega_h::Library lib; + auto source_mesh = BuildUnitSquare(lib, 1); + auto target_mesh = BuildUnitSquare(lib, 0); + + auto source_space = pcms::LagrangeFunctionSpace::FromMesh( + source_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", + pcms::LagrangeFunctionSpace::Backend::OmegaH); + auto target_space = pcms::LagrangeFunctionSpace::FromMesh( + target_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", + pcms::LagrangeFunctionSpace::Backend::OmegaH); + + const double c = 2.0; + auto source_field = source_space.CreateField(); + pcms::test::SetField( + source_field, OMEGA_H_LAMBDA(pcms::Real, pcms::Real) { return c; }); + + const auto target_layout = + std::dynamic_pointer_cast( + target_space.GetLayout()); + const std::vector g(target_mesh.nverts(), c); + const auto expected = ExpectedLoadByGid(target_mesh, *target_layout, g); + + auto integrator = + pcms::BuildOmegaHConservativeRHSIntegrator(source_space, target_space); + CheckAssembledLoadMatches(source_space, *integrator, source_field, expected); +} + +TEST_CASE( + "OmegaHIntersectionRHSIntegrator: linear field matches consistent-mass load", + "[rhs_integrator]") +{ + // f(x,y) = x + y exercises non-trivial quadrature paths that a constant + // field cannot catch (a constant integrates exactly under any rule, while a + // linear field requires at least 1st-order accuracy). The field is affine, + // hence exactly representable in the target P1 space, so the assembled load + // vector must equal the target consistent mass matrix applied to the nodal + // values of x + y. + Omega_h::Library lib; + auto source_mesh = BuildUnitSquare(lib, 1); + auto target_mesh = BuildUnitSquare(lib, 0); + + auto source_space = pcms::LagrangeFunctionSpace::FromMesh( + source_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", + pcms::LagrangeFunctionSpace::Backend::OmegaH); + auto target_space = pcms::LagrangeFunctionSpace::FromMesh( + target_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", + pcms::LagrangeFunctionSpace::Backend::OmegaH); + + auto source_field = source_space.CreateField(); + pcms::test::SetField( + source_field, OMEGA_H_LAMBDA(pcms::Real x, pcms::Real y) { return x + y; }); + + const auto target_layout = + std::dynamic_pointer_cast( + target_space.GetLayout()); + const auto tgt_coords_h = + Omega_h::HostRead(target_mesh.coords()); + std::vector g(target_mesh.nverts()); + for (int i = 0; i < target_mesh.nverts(); ++i) { + g[i] = tgt_coords_h[2 * i + 0] + tgt_coords_h[2 * i + 1]; + } + const auto expected = ExpectedLoadByGid(target_mesh, *target_layout, g); + + auto integrator = + pcms::BuildOmegaHConservativeRHSIntegrator(source_space, target_space); + CheckAssembledLoadMatches(source_space, *integrator, source_field, expected); +} + +TEST_CASE("OmegaHIntersectionRHSIntegrator: rejects invalid layouts", + "[rhs_integrator]") +{ + Omega_h::Library lib; + auto source_mesh = BuildUnitSquare(lib, 1); + auto target_mesh = BuildUnitSquare(lib, 0); + + SECTION("multi-component source space throws") + { + auto source_space = pcms::LagrangeFunctionSpace::FromMesh( + source_mesh, 1, 2, pcms::CoordinateSystem::Cartesian, "global", + pcms::LagrangeFunctionSpace::Backend::OmegaH); + auto target_space = pcms::LagrangeFunctionSpace::FromMesh( + target_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", + pcms::LagrangeFunctionSpace::Backend::OmegaH); + REQUIRE_THROWS( + pcms::BuildOmegaHConservativeRHSIntegrator(source_space, target_space)); + } + + SECTION("multi-component target space throws") + { + auto source_space = pcms::LagrangeFunctionSpace::FromMesh( + source_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", + pcms::LagrangeFunctionSpace::Backend::OmegaH); + auto target_space = pcms::LagrangeFunctionSpace::FromMesh( + target_mesh, 1, 2, pcms::CoordinateSystem::Cartesian, "global", + pcms::LagrangeFunctionSpace::Backend::OmegaH); + REQUIRE_THROWS( + pcms::BuildOmegaHConservativeRHSIntegrator(source_space, target_space)); + } + + SECTION("non-Cartesian source coordinate system throws") + { + auto source_space = pcms::LagrangeFunctionSpace::FromMesh( + source_mesh, 1, 1, pcms::CoordinateSystem::Cylindrical, "global", + pcms::LagrangeFunctionSpace::Backend::OmegaH); + auto target_space = pcms::LagrangeFunctionSpace::FromMesh( + target_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", + pcms::LagrangeFunctionSpace::Backend::OmegaH); + REQUIRE_THROWS( + pcms::BuildOmegaHConservativeRHSIntegrator(source_space, target_space)); + } + + SECTION("non-Cartesian target coordinate system throws") + { + auto source_space = pcms::LagrangeFunctionSpace::FromMesh( + source_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", + pcms::LagrangeFunctionSpace::Backend::OmegaH); + auto target_space = pcms::LagrangeFunctionSpace::FromMesh( + target_mesh, 1, 1, pcms::CoordinateSystem::Cylindrical, "global", + pcms::LagrangeFunctionSpace::Backend::OmegaH); + REQUIRE_THROWS( + pcms::BuildOmegaHConservativeRHSIntegrator(source_space, target_space)); + } +} diff --git a/test/test_omega_h_mass_integrator.cpp b/test/test_omega_h_mass_integrator.cpp new file mode 100644 index 00000000..b916651c --- /dev/null +++ b/test/test_omega_h_mass_integrator.cpp @@ -0,0 +1,187 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// --------------------------------------------------------------------------- +// Helpers +// --------------------------------------------------------------------------- + +namespace +{ + +Omega_h::Mesh BuildUnitSquare(Omega_h::Library& lib) +{ + const Omega_h::Reals coords({0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 1.0}); + Omega_h::LOs ev2v({0, 1, 3, 1, 2, 3}); + Omega_h::Mesh mesh(&lib); + Omega_h::build_from_elems_and_coords(&mesh, OMEGA_H_SIMPLEX, 2, ev2v, coords); + mesh.add_tag( + 0, "class_dim", 1, + Omega_h::Read(mesh.nverts(), Omega_h::I8(0))); + mesh.add_tag( + 0, "class_id", 1, + Omega_h::Read(mesh.nverts(), Omega_h::ClassId(0))); + mesh.add_tag( + 1, "class_dim", 1, + Omega_h::Read(mesh.nedges(), Omega_h::I8(1))); + mesh.add_tag( + 1, "class_id", 1, + Omega_h::Read(mesh.nedges(), Omega_h::ClassId(0))); + mesh.add_tag( + 2, "class_dim", 1, + Omega_h::Read(mesh.nelems(), Omega_h::I8(2))); + mesh.add_tag( + 2, "class_id", 1, + Omega_h::Read(mesh.nelems(), Omega_h::ClassId(0))); + return mesh; +} + +std::map, pcms::Real> BuildReferenceMassMap( + Omega_h::Mesh& mesh, const pcms::FunctionSpace& space) +{ + using namespace pcms; + + MeshField::OmegahMeshField + omf(mesh); + auto coordField = omf.getCoordField(); + const auto [shp, map] = MeshField::Omegah::getTriangleElement<1>(mesh); + MeshField::FieldElement coordFe(mesh.nelems(), coordField, shp, map); + auto elm_mass_dev = buildElementMassMatrix(mesh, coordFe); + auto elm_mass_host = + Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, elm_mass_dev); + + const auto layout = + std::dynamic_pointer_cast(space.GetLayout()); + REQUIRE(layout != nullptr); + const auto& faces2nodes = mesh.ask_down(Omega_h::FACE, Omega_h::VERT).ab2b; + const auto& gids = layout->GetGidsHost(); + + std::map, Real> ref; + for (int e = 0; e < mesh.nelems(); ++e) { + const auto verts = Omega_h::gather_verts<3>(faces2nodes, e); + for (int i = 0; i < 3; ++i) { + const GO row = static_cast(gids[verts[i]]); + for (int j = 0; j < 3; ++j) { + const GO col = static_cast(gids[verts[j]]); + ref[{row, col}] += elm_mass_host(e * 9 + i * 3 + j); + } + } + } + return ref; +} + +} // namespace + +// --------------------------------------------------------------------------- +// Tests +// --------------------------------------------------------------------------- + +TEST_CASE( + "OmegaHMassIntegrator: assembled matrix matches buildElementMassMatrix", + "[mass_integrator]") +{ + Omega_h::Library lib; + auto mesh = BuildUnitSquare(lib); + + auto space = pcms::LagrangeFunctionSpace::FromMesh( + mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", + pcms::LagrangeFunctionSpace::Backend::OmegaH); + + const auto ref = BuildReferenceMassMap(mesh, space); + + auto integrator = pcms::BuildOmegaHMassIntegrator(space); + Mat mat = integrator->GetMatrix(); + + for (const auto& [key, ref_val] : ref) { + PetscInt r = static_cast(key.first); + PetscInt c = static_cast(key.second); + PetscScalar got = 0.0; + MatGetValues(mat, 1, &r, 1, &c, &got); + CAPTURE(key.first, key.second, ref_val, got); + CHECK(static_cast(got) == + Catch::Approx(ref_val).epsilon(1e-12)); + } + // mat is owned by integrator; no MatDestroy here. +} + +TEST_CASE("OmegaHMassIntegrator: row sums match lumped mass", + "[mass_integrator]") +{ + // For a consistent mass matrix M, the row sums equal the lumped mass at + // each node: row_sum(i) = integral(N_i dx). For a regular mesh of unit + // area with uniform nodal distribution the sum over all nodes equals 1. + Omega_h::Library lib; + auto mesh = BuildUnitSquare(lib); + + auto space = pcms::LagrangeFunctionSpace::FromMesh( + mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", + pcms::LagrangeFunctionSpace::Backend::OmegaH); + + const auto layout = + std::dynamic_pointer_cast( + space.GetLayout()); + const auto& gids = layout->GetGidsHost(); + const int nverts = mesh.nverts(); + + auto integrator = pcms::BuildOmegaHMassIntegrator(space); + Mat mat = integrator->GetMatrix(); + + // Sum all row sums. + pcms::Real grand_total = 0.0; + for (int v = 0; v < nverts; ++v) { + const pcms::GO row = static_cast(gids[v]); + pcms::Real row_sum = 0.0; + for (int u = 0; u < nverts; ++u) { + const pcms::GO col = static_cast(gids[u]); + PetscInt r = static_cast(row); + PetscInt c = static_cast(col); + PetscScalar val = 0.0; + MatGetValues(mat, 1, &r, 1, &c, &val); + row_sum += static_cast(val); + } + CAPTURE(v, row, row_sum); + CHECK(row_sum > 0.0); // each row sum must be positive + grand_total += row_sum; + } + + // Total == area of domain == 1. + CHECK(grand_total == Catch::Approx(1.0).epsilon(1e-10)); + // mat is owned by integrator; no MatDestroy here. +} + +TEST_CASE("OmegaHMassIntegrator: rejects invalid layouts", "[mass_integrator]") +{ + Omega_h::Library lib; + auto mesh = BuildUnitSquare(lib); + + SECTION("multi-component space throws") + { + auto space = pcms::LagrangeFunctionSpace::FromMesh( + mesh, 1, 2, pcms::CoordinateSystem::Cartesian, "global", + pcms::LagrangeFunctionSpace::Backend::OmegaH); + REQUIRE_THROWS(pcms::BuildOmegaHMassIntegrator(space)); + } + + SECTION("non-Cartesian coordinate system throws") + { + auto space = pcms::LagrangeFunctionSpace::FromMesh( + mesh, 1, 1, pcms::CoordinateSystem::Cylindrical, "global", + pcms::LagrangeFunctionSpace::Backend::OmegaH); + REQUIRE_THROWS(pcms::BuildOmegaHMassIntegrator(space)); + } +} From 08d6143e2b35262508d000b67e4ef962668a6365 Mon Sep 17 00:00:00 2001 From: Jacob Merson Date: Sat, 11 Jul 2026 11:11:58 -0400 Subject: [PATCH 03/15] transfer: add Monte-Carlo / control-variate conservative projection Squashed original commits: 844dc1c, b4a4545, 629419f (sobol_sequence churn dropped) --- src/pcms/transfer/interpolator.h | 2 +- .../omega_h_control_variate_projection.cpp | 89 +++++++ .../omega_h_control_variate_projection.hpp | 65 +++++ .../transfer/omega_h_mc_rhs_integrator.cpp | 211 +++++++++++++++ .../transfer/omega_h_mc_rhs_integrator.hpp | 84 ++++++ test/test_omega_h_mc_rhs_integrator.cpp | 242 ++++++++++++++++++ 6 files changed, 692 insertions(+), 1 deletion(-) create mode 100644 src/pcms/transfer/omega_h_control_variate_projection.cpp create mode 100644 src/pcms/transfer/omega_h_control_variate_projection.hpp create mode 100644 src/pcms/transfer/omega_h_mc_rhs_integrator.cpp create mode 100644 src/pcms/transfer/omega_h_mc_rhs_integrator.hpp create mode 100644 test/test_omega_h_mc_rhs_integrator.cpp diff --git a/src/pcms/transfer/interpolator.h b/src/pcms/transfer/interpolator.h index a8b5b07e..25ef3e1e 100644 --- a/src/pcms/transfer/interpolator.h +++ b/src/pcms/transfer/interpolator.h @@ -12,7 +12,7 @@ #include "pcms/utility/types.h" #include #include -#include +#include "pcms/transfer/transfer_operator.hpp" namespace pcms { diff --git a/src/pcms/transfer/omega_h_control_variate_projection.cpp b/src/pcms/transfer/omega_h_control_variate_projection.cpp new file mode 100644 index 00000000..c22b3f71 --- /dev/null +++ b/src/pcms/transfer/omega_h_control_variate_projection.cpp @@ -0,0 +1,89 @@ +#include "pcms/transfer/omega_h_control_variate_projection.hpp" +#include "pcms/transfer/conservative_projection_solver.hpp" +#include "pcms/transfer/omega_h_mass_integrator.hpp" +#include "pcms/utility/arrays.h" +#include +#include + +namespace pcms +{ + +OmegaHControlVariateProjection::OmegaHControlVariateProjection( + const FunctionSpace& source_space, const FunctionSpace& target_space, + int samples_per_element, MonteCarloSampling sampling, uint64_t seed) + : target_layout_(std::dynamic_pointer_cast( + target_space.GetLayout())), + rhs_integrator_(std::make_unique( + target_layout_, target_space.GetCoordinateSystem(), samples_per_element, + sampling, seed)), + interpolator_(source_space, target_space), + control_variate_(target_space.CreateField()) +{ + const auto sample_coords = + rhs_integrator_->GetIntegrationPoints().GetCoordinates(); + source_at_samples_ = source_space.CreatePointEvaluator( + EvaluationRequest::FromCoordinates(sample_coords)); + control_variate_at_samples_ = target_space.CreatePointEvaluator( + EvaluationRequest::FromCoordinates(sample_coords)); + + // 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()); + solver_ = std::make_unique(mass_integrator, + *rhs_integrator_); +} + +// Defined here so that GalerkinProjectionSolver (forward-declared in the +// header) is a complete type when unique_ptr's destructor is instantiated. +OmegaHControlVariateProjection::~OmegaHControlVariateProjection() = default; + +void OmegaHControlVariateProjection::Apply(const Field& source, + Field& target) const +{ + // 1. Control variate: interpolate the source field onto the target space. + interpolator_.Apply(source, control_variate_); + + // 2. Sample the source field and the control variate at the fixed Monte + // Carlo sample points; the stochastic RHS integrates the residual. + const std::size_t num_samples = + rhs_integrator_->GetIntegrationPoints().NumPoints(); + Kokkos::View f_samples("cv_f_samples", num_samples, + 1); + Kokkos::View residual("cv_residual", num_samples, + 1); + source_at_samples_->Evaluate(source, MakeRank2View(f_samples)); + control_variate_at_samples_->Evaluate(control_variate_, + MakeRank2View(residual)); + Kokkos::parallel_for( + "cv_residual", Kokkos::RangePolicy(0, num_samples), + KOKKOS_LAMBDA(int i) { + residual(i, 0) = f_samples(i, 0) - residual(i, 0); + }); + + // 3. Project the residual: M * delta = r. + const auto delta = solver_->Solve(MakeConstRank2View(residual)); + + // 4. x = g + delta. delta is indexed by global id (PETSc row), the DOF + // holder by local vertex id. + const auto g_nodal = control_variate_.GetDOFHolderData(); + auto g_view = Kokkos::View>( + g_nodal.data_handle(), g_nodal.extent(0)); + const auto device_gids = target_layout_->GetGids(); + const GO* gids_ptr = device_gids.data_handle(); + const int nverts = static_cast(g_view.extent(0)); + + Kokkos::View result("cv_result", nverts); + Kokkos::parallel_for( + "cv_add_correction", Kokkos::RangePolicy(0, nverts), + KOKKOS_LAMBDA(int i) { + result(i) = g_view(i) + delta[static_cast(gids_ptr[i])]; + }); + Kokkos::fence(); + + target.SetDOFHolderData( + Rank2View(result.data(), nverts, 1)); +} + +} // namespace pcms diff --git a/src/pcms/transfer/omega_h_control_variate_projection.hpp b/src/pcms/transfer/omega_h_control_variate_projection.hpp new file mode 100644 index 00000000..dd3c93ff --- /dev/null +++ b/src/pcms/transfer/omega_h_control_variate_projection.hpp @@ -0,0 +1,65 @@ +#ifndef PCMS_TRANSFER_OMEGA_H_CONTROL_VARIATE_PROJECTION_HPP +#define PCMS_TRANSFER_OMEGA_H_CONTROL_VARIATE_PROJECTION_HPP + +#include "pcms/field/field.h" +#include "pcms/field/function_space.h" +#include "pcms/field/layout/omega_h_lagrange.h" +#include "pcms/field/point_evaluator.h" +#include "pcms/transfer/interpolator.h" +#include "pcms/transfer/omega_h_mc_rhs_integrator.hpp" +#include "pcms/transfer/transfer_operator.hpp" +#include +#include + +namespace pcms +{ + +// Forward declaration keeps PETSc headers out of this header. +class GalerkinProjectionSolver; + +// Variance-reduced Monte Carlo Galerkin projection using a control variate. +// +// The control variate g is the interpolation of the source field onto the +// target space (Interpolator: evaluate the source field at the target nodal +// locations, interpolate with the target basis). Since g lives in the target +// space, its L2 projection is g itself, so only the residual f - g needs the +// stochastic RHS: +// M x = M g + \int phi (f - g) => x = g + M^{-1} r, +// r_j ~= (|T| / N) * sum_i phi_j(x_i) (f(x_i) - g(x_i)) +// Because g approximates f, the sampled residual has far smaller variance +// than sampling f directly (and is exactly zero when f is in the target +// space). +// +// Apply() per-call cost: one source-field evaluation at the target nodes +// (interpolation), one source-field and one target-field evaluation at the +// fixed sample points, RHS assembly, and a KSPSolve with the cached +// factorization. +class OmegaHControlVariateProjection : public TransferOperator +{ +public: + OmegaHControlVariateProjection(const FunctionSpace& source_space, + const FunctionSpace& target_space, + int samples_per_element, + MonteCarloSampling sampling, + uint64_t seed = 8675309); + + // Defined in the .cpp where GalerkinProjectionSolver is a complete type. + ~OmegaHControlVariateProjection() override; + + void Apply(const Field& source, Field& target) const override; + +private: + std::shared_ptr target_layout_; + std::unique_ptr rhs_integrator_; + Interpolator interpolator_; + // Scratch field on the target space holding the interpolated control + // variate; overwritten on every Apply. + mutable Field control_variate_; + std::unique_ptr> source_at_samples_; + std::unique_ptr> control_variate_at_samples_; + std::unique_ptr solver_; +}; + +} // namespace pcms + +#endif // PCMS_TRANSFER_OMEGA_H_CONTROL_VARIATE_PROJECTION_HPP diff --git a/src/pcms/transfer/omega_h_mc_rhs_integrator.cpp b/src/pcms/transfer/omega_h_mc_rhs_integrator.cpp new file mode 100644 index 00000000..212539fc --- /dev/null +++ b/src/pcms/transfer/omega_h_mc_rhs_integrator.cpp @@ -0,0 +1,211 @@ +#include "pcms/transfer/omega_h_mc_rhs_integrator.hpp" +#include "pcms/transfer/omega_h_form_integrator_utils.hpp" +#include "pcms/transfer/petsc_utils.hpp" +#include "pcms/utility/assert.h" +#include +#include +#include + +namespace pcms +{ + +namespace +{ + +// Maps a uniform point on the unit square to uniform barycentric coordinates +// Shape Distributions (ACM Transactions on Graphics, Vol. 21, No. 4, October +// 2002.) page 814 Eq 1 +KOKKOS_INLINE_FUNCTION Omega_h::Vector<3> UniformTriangleBarycentric(Real u, + Real v) +{ + const Real s = Kokkos::sqrt(u); + return {1.0 - s, s * (1.0 - v), s * v}; +} + +// Fills coords, node_gids, and coeffs for all samples of one target element. +// unit_sample_at(s, u, v) provides the s-th unit-square sample. +template +KOKKOS_INLINE_FUNCTION void FillElementSamples( + const int elm, const int samples_per_element, + const Omega_h::Reals& mesh_coords, const Omega_h::LOs& faces2nodes, + const GO* gids, const Kokkos::View& coords, + const Kokkos::View& node_gids, + const Kokkos::View& coeffs, + const UnitSampleAt& unit_sample_at) +{ + const auto verts = Omega_h::gather_verts<3>(faces2nodes, elm); + const auto vert_coords = Omega_h::gather_vectors<3, 2>(mesh_coords, verts); + Omega_h::Few, 2> basis; + basis[0] = vert_coords[1] - vert_coords[0]; + basis[1] = vert_coords[2] - vert_coords[0]; + const Real area = Kokkos::fabs(Omega_h::triangle_area_from_basis(basis)); + const Real weight = area / samples_per_element; + + for (int s = 0; s < samples_per_element; ++s) { + Real u = 0.0; + Real v = 0.0; + unit_sample_at(s, u, v); + const auto bary = UniformTriangleBarycentric(u, v); + + const int i = elm * samples_per_element + s; + Real x = 0.0; + Real y = 0.0; + for (int k = 0; k < 3; ++k) { + x += bary[k] * vert_coords[k][0]; + y += bary[k] * vert_coords[k][1]; + node_gids(i * 3 + k) = static_cast(gids[verts[k]]); + coeffs(i * 3 + k) = bary[k] * weight; + } + coords(i, 0) = x; + coords(i, 1) = y; + } +} + +} // namespace + +// --------------------------------------------------------------------------- +// OmegaHMonteCarloRHSIntegrator::BuildData +// --------------------------------------------------------------------------- + +OmegaHMonteCarloRHSIntegrator::Data OmegaHMonteCarloRHSIntegrator::BuildData( + const std::shared_ptr& target_layout, + CoordinateSystem target_coordinate_system, int samples_per_element, + MonteCarloSampling sampling, uint64_t seed) +{ + detail::CheckOmegaHScalarP1Layout(target_coordinate_system, target_layout, + "OmegaHMonteCarloRHSIntegrator", "target"); + if (samples_per_element <= 0) { + throw pcms_error( + "OmegaHMonteCarloRHSIntegrator: samples_per_element must be positive"); + } + + Omega_h::Mesh& mesh = target_layout->GetMesh(); + + const int nelems = mesh.nelems(); + const int num_samples = nelems * samples_per_element; + const auto mesh_coords = mesh.coords(); + const auto faces2nodes = mesh.ask_down(Omega_h::FACE, Omega_h::VERT).ab2b; + const auto device_gids = target_layout->GetGids(); + const GO* gids_ptr = device_gids.data_handle(); + + Data d; + d.coords = + Kokkos::View("mc_rhs_coords", num_samples, 2); + d.node_gids = Kokkos::View("mc_rhs_node_gids", + num_samples * 3); + d.coeffs = + Kokkos::View("mc_rhs_coeffs", num_samples * 3); + d.nverts = static_cast(mesh.nverts()); + + const int npe = samples_per_element; + auto coords = d.coords; + auto node_gids = d.node_gids; + auto coeffs = d.coeffs; + + Kokkos::Random_XorShift64_Pool pool(seed); + Kokkos::parallel_for( + "mc_rhs_fill_random", Kokkos::RangePolicy(0, nelems), + KOKKOS_LAMBDA(int elm) { + auto gen = pool.get_state(); + FillElementSamples(elm, npe, mesh_coords, faces2nodes, gids_ptr, coords, + node_gids, coeffs, [&](int /*s*/, Real& u, Real& v) { + u = gen.drand(); + v = gen.drand(); + }); + pool.free_state(gen); + }); + Kokkos::fence(); + + return d; +} + +// --------------------------------------------------------------------------- +// OmegaHMonteCarloRHSIntegrator constructors +// --------------------------------------------------------------------------- + +OmegaHMonteCarloRHSIntegrator::OmegaHMonteCarloRHSIntegrator( + const FunctionSpace& target_space, int samples_per_element, + MonteCarloSampling sampling, uint64_t seed) + : OmegaHMonteCarloRHSIntegrator( + std::dynamic_pointer_cast( + target_space.GetLayout()), + target_space.GetCoordinateSystem(), samples_per_element, sampling, seed) +{ +} + +OmegaHMonteCarloRHSIntegrator::OmegaHMonteCarloRHSIntegrator( + std::shared_ptr target_layout, + CoordinateSystem target_coordinate_system, int samples_per_element, + MonteCarloSampling sampling, uint64_t seed) + : OmegaHMonteCarloRHSIntegrator( + BuildData(target_layout, target_coordinate_system, samples_per_element, + sampling, seed)) +{ +} + +OmegaHMonteCarloRHSIntegrator::OmegaHMonteCarloRHSIntegrator(Data data) + : point_set_(CoordinateSystem::Cartesian, std::move(data.coords)), + node_gids_(std::move(data.node_gids)), + coeffs_(std::move(data.coeffs)) +{ + const PetscInt nnz = static_cast(node_gids_.extent(0)); + PetscErrorCode ierr = createSeqVec(PETSC_COMM_WORLD, data.nverts, &vec_); + CHKERRABORT(PETSC_COMM_WORLD, ierr); + ierr = VecSetPreallocationCOO(vec_, nnz, node_gids_.data()); + CHKERRABORT(PETSC_COMM_WORLD, ierr); +} + +OmegaHMonteCarloRHSIntegrator::~OmegaHMonteCarloRHSIntegrator() +{ + if (vec_) { + VecDestroy(&vec_); + } +} + +// --------------------------------------------------------------------------- +// OmegaHMonteCarloRHSIntegrator public interface +// --------------------------------------------------------------------------- + +const IntegrationPointSet& +OmegaHMonteCarloRHSIntegrator::GetIntegrationPoints() const noexcept +{ + return point_set_; +} + +Vec OmegaHMonteCarloRHSIntegrator::GetVector() const noexcept +{ + return vec_; +} + +void OmegaHMonteCarloRHSIntegrator::Assemble( + Rank2View sampled_values) +{ + const std::size_t num_samples = + static_cast(node_gids_.extent(0) / 3); + PCMS_ALWAYS_ASSERT(static_cast(sampled_values.extent(0)) == + num_samples); + PCMS_ALWAYS_ASSERT(sampled_values.extent(1) >= 1); + + PetscErrorCode ierr = VecZeroEntries(vec_); + CHKERRABORT(PETSC_COMM_WORLD, ierr); + + auto sv = Kokkos::View>( + sampled_values.data_handle(), sampled_values.extent(0), + sampled_values.extent(1)); + Kokkos::View coo_vals("mc_rhs_coo_vals", + num_samples * 3); + auto coeffs = coeffs_; + Kokkos::parallel_for( + "mc_rhs_coo_vals", static_cast(num_samples), KOKKOS_LAMBDA(int i) { + const PetscScalar f = static_cast(sv(i, 0)); + for (int j = 0; j < 3; ++j) { + coo_vals(i * 3 + j) = static_cast(coeffs(i * 3 + j)) * f; + } + }); + + ierr = VecSetValuesCOO(vec_, coo_vals.data(), ADD_VALUES); + CHKERRABORT(PETSC_COMM_WORLD, ierr); +} + +} // namespace pcms diff --git a/src/pcms/transfer/omega_h_mc_rhs_integrator.hpp b/src/pcms/transfer/omega_h_mc_rhs_integrator.hpp new file mode 100644 index 00000000..1fed3ef4 --- /dev/null +++ b/src/pcms/transfer/omega_h_mc_rhs_integrator.hpp @@ -0,0 +1,84 @@ +#ifndef PCMS_TRANSFER_OMEGA_H_MC_RHS_INTEGRATOR_HPP +#define PCMS_TRANSFER_OMEGA_H_MC_RHS_INTEGRATOR_HPP + +#include "pcms/field/function_space.h" +#include "pcms/field/layout/omega_h_lagrange.h" +#include "pcms/transfer/integration_point_set.hpp" +#include "pcms/transfer/linear_form_integrator.hpp" +#include +#include +#include + +namespace pcms +{ + +// Strategy for generating sample points on each target element. +enum class MonteCarloSampling +{ + UniformRandom, // independent uniform draws per element (XorShift64 pool) +}; + +// Monte Carlo RHS integrator for conservative L2 projection onto order-1 +// Lagrange spaces on Omega_h 2D simplex meshes. +// +// Instead of intersection-based quadrature, each load-vector entry +// b_j = \int phi_j f dx +// is estimated per target element T from a uniform sample x_1..x_N over T: +// b_j ~= (|T| / N) * sum_i phi_j(x_i) f(x_i) +// Unit-square samples (pseudo-random) are mapped to each triangle +// with the area-preserving sqrt parametrization. Sample locations are fixed +// at construction and exposed through GetIntegrationPoints(); Assemble() +// consumes source-field values sampled there, exactly like +// OmegaHIntersectionRHSIntegrator, so this drops into +// GalerkinProjectionSolver::Solve. +class OmegaHMonteCarloRHSIntegrator : public LinearFormIntegrator +{ +public: + OmegaHMonteCarloRHSIntegrator(const FunctionSpace& target_space, + int samples_per_element, + MonteCarloSampling sampling, + uint64_t seed = 8675309); + OmegaHMonteCarloRHSIntegrator( + std::shared_ptr target_layout, + CoordinateSystem target_coordinate_system, int samples_per_element, + MonteCarloSampling sampling, uint64_t seed = 8675309); + ~OmegaHMonteCarloRHSIntegrator(); + + OmegaHMonteCarloRHSIntegrator(const OmegaHMonteCarloRHSIntegrator&) = delete; + OmegaHMonteCarloRHSIntegrator& operator=( + const OmegaHMonteCarloRHSIntegrator&) = delete; + + Vec GetVector() const noexcept override; + + const IntegrationPointSet& GetIntegrationPoints() + const noexcept override; + + void Assemble( + Rank2View sampled_values) override; + +private: + // Internal data bundle produced by the sample-generation kernel. + struct Data + { + Kokkos::View coords; // [num_samples, 2] + Kokkos::View node_gids; // [num_samples * 3] + Kokkos::View coeffs; // [num_samples * 3] + PetscInt nverts = 0; // target DOF count + }; + + static Data BuildData( + const std::shared_ptr& target_layout, + CoordinateSystem target_coordinate_system, int samples_per_element, + MonteCarloSampling sampling, uint64_t seed); + + explicit OmegaHMonteCarloRHSIntegrator(Data data); + + Vec vec_ = nullptr; + IntegrationPointSet point_set_; + Kokkos::View node_gids_; // COO indices + Kokkos::View coeffs_; // basis * |T| / N +}; + +} // namespace pcms + +#endif // PCMS_TRANSFER_OMEGA_H_MC_RHS_INTEGRATOR_HPP diff --git a/test/test_omega_h_mc_rhs_integrator.cpp b/test/test_omega_h_mc_rhs_integrator.cpp new file mode 100644 index 00000000..a26095ce --- /dev/null +++ b/test/test_omega_h_mc_rhs_integrator.cpp @@ -0,0 +1,242 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "field_test_utils.h" +#include +#include + +// --------------------------------------------------------------------------- +// Helpers +// --------------------------------------------------------------------------- + +namespace +{ + +// Build a unit-square 2D simplex mesh. +// diagonal=0: T0=(0,1,3), T1=(1,2,3) +// diagonal=1: T0=(0,1,2), T1=(0,2,3) +Omega_h::Mesh BuildUnitSquare(Omega_h::Library& lib, int diagonal) +{ + const Omega_h::Reals coords({0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 1.0}); + Omega_h::LOs ev2v = (diagonal == 0) ? Omega_h::LOs({0, 1, 3, 1, 2, 3}) + : Omega_h::LOs({0, 1, 2, 0, 2, 3}); + Omega_h::Mesh mesh(&lib); + Omega_h::build_from_elems_and_coords(&mesh, OMEGA_H_SIMPLEX, 2, ev2v, coords); + mesh.add_tag( + 0, "class_dim", 1, + Omega_h::Read(mesh.nverts(), Omega_h::I8(0))); + mesh.add_tag( + 0, "class_id", 1, + Omega_h::Read(mesh.nverts(), Omega_h::ClassId(0))); + mesh.add_tag( + 1, "class_dim", 1, + Omega_h::Read(mesh.nedges(), Omega_h::I8(1))); + mesh.add_tag( + 1, "class_id", 1, + Omega_h::Read(mesh.nedges(), Omega_h::ClassId(0))); + mesh.add_tag( + 2, "class_dim", 1, + Omega_h::Read(mesh.nelems(), Omega_h::I8(2))); + mesh.add_tag( + 2, "class_id", 1, + Omega_h::Read(mesh.nelems(), Omega_h::ClassId(0))); + return mesh; +} + +pcms::LagrangeFunctionSpace MakeP1Space(Omega_h::Mesh& mesh) +{ + return pcms::LagrangeFunctionSpace::FromMesh( + mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", + pcms::LagrangeFunctionSpace::Backend::OmegaH); +} + +// Evaluates source_field at the integrator's sample points and assembles. +void EvaluateAndAssemble(pcms::LinearFormIntegrator& integrator, + const pcms::FunctionSpace& source_space, + const pcms::Field& source_field) +{ + const auto& pts = integrator.GetIntegrationPoints(); + const std::size_t npts = pts.GetCoordinates().GetValues().extent(0); + auto evaluator = source_space.CreatePointEvaluator( + pcms::EvaluationRequest::FromCoordinates(pts.GetCoordinates())); + Kokkos::View sampled("sampled", npts, + 1); + evaluator->Evaluate(source_field, pcms::MakeRank2View(sampled)); + integrator.Assemble(pcms::MakeConstRank2View(sampled)); +} + +} // namespace + +// --------------------------------------------------------------------------- +// Monte Carlo RHS integrator +// --------------------------------------------------------------------------- + +TEST_CASE("OmegaHMonteCarloRHSIntegrator: sample points lie inside the domain", + "[mc_rhs_integrator]") +{ + Omega_h::Library lib; + auto target_mesh = BuildUnitSquare(lib, 0); + auto target_space = MakeP1Space(target_mesh); + + const int samples_per_element = 16; + for (const auto sampling : {pcms::MonteCarloSampling::UniformRandom}) { + pcms::OmegaHMonteCarloRHSIntegrator integrator( + target_space, samples_per_element, sampling); + const auto raw_coords = + integrator.GetIntegrationPoints().GetCoordinates().GetValues(); + auto coords_view = Kokkos::View>( + raw_coords.data_handle(), raw_coords.extent(0), raw_coords.extent(1)); + auto coords_h = + Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, coords_view); + + REQUIRE( + coords_h.extent(0) == + static_cast(target_mesh.nelems() * samples_per_element)); + for (std::size_t i = 0; i < coords_h.extent(0); ++i) { + CAPTURE(i, coords_h(i, 0), coords_h(i, 1)); + CHECK(coords_h(i, 0) >= -1e-12); + CHECK(coords_h(i, 0) <= 1.0 + 1e-12); + CHECK(coords_h(i, 1) >= -1e-12); + CHECK(coords_h(i, 1) <= 1.0 + 1e-12); + } + } +} + +TEST_CASE("OmegaHMonteCarloRHSIntegrator: constant field integrates exactly", + "[mc_rhs_integrator]") +{ + // For f = c the estimator sums to c * |domain| for any sample placement, + // because the P1 basis functions partition unity at every sample point. + Omega_h::Library lib; + auto source_mesh = BuildUnitSquare(lib, 1); + auto target_mesh = BuildUnitSquare(lib, 0); + auto source_space = MakeP1Space(source_mesh); + auto target_space = MakeP1Space(target_mesh); + + auto source_field = source_space.CreateField(); + pcms::test::SetField( + source_field, OMEGA_H_LAMBDA(pcms::Real, pcms::Real) { return 2.0; }); + + for (const auto sampling : {pcms::MonteCarloSampling::UniformRandom}) { + pcms::OmegaHMonteCarloRHSIntegrator integrator(target_space, 8, sampling); + EvaluateAndAssemble(integrator, source_space, source_field); + + PetscScalar sum = 0.0; + VecSum(integrator.GetVector(), &sum); + INFO("sampling=" << "UniformRandom"); + CHECK(static_cast(sum) == Catch::Approx(2.0).margin(1e-12)); + } +} + +// --------------------------------------------------------------------------- +// Control variate projection +// --------------------------------------------------------------------------- + +TEST_CASE("OmegaHControlVariateProjection: exact for fields in the target " + "space even with few samples", + "[mc_rhs_integrator][control_variate]") +{ + // For a globally linear field the control variate (target-space interpolant + // of the source field) equals the source field, so the sampled residual is + // identically zero and the projection is exact regardless of sample count. + Omega_h::Library lib; + auto source_mesh = BuildUnitSquare(lib, 1); + auto target_mesh = BuildUnitSquare(lib, 0); + auto source_space = MakeP1Space(source_mesh); + auto target_space = MakeP1Space(target_mesh); + + auto source_field = source_space.CreateField(); + auto target_field = target_space.CreateField(); + pcms::test::SetField( + source_field, + OMEGA_H_LAMBDA(pcms::Real x, pcms::Real y) { return 3.0 * x - y + 0.25; }); + + pcms::OmegaHControlVariateProjection projection( + source_space, target_space, /*samples_per_element=*/4, + pcms::MonteCarloSampling::UniformRandom); + projection.Apply(source_field, target_field); + + const auto values = + pcms::FlattenToRank1View(target_field.GetDOFHolderDataHost()); + const auto coords_h = Omega_h::HostRead(target_mesh.coords()); + REQUIRE(static_cast(values.size()) == target_mesh.nverts()); + for (Omega_h::LO i = 0; i < target_mesh.nverts(); ++i) { + const double expected = + 3.0 * coords_h[2 * i + 0] - coords_h[2 * i + 1] + 0.25; + CAPTURE(i, expected, values[i]); + CHECK(values[i] == Catch::Approx(expected).margin(1e-9)); + } +} + +TEST_CASE("OmegaHControlVariateProjection: reduces error vs plain Monte Carlo", + "[mc_rhs_integrator][control_variate]") +{ + // For a non-linear field the control variate does not cancel f exactly, but + // the residual is much smaller than f, so with identical sample counts and + // seeds the control-variate projection must be closer to the exact + // (intersection-quadrature) projection than the plain Monte Carlo one. + Omega_h::Library lib; + auto source_mesh = BuildUnitSquare(lib, 1); + auto target_mesh = BuildUnitSquare(lib, 0); + auto source_space = MakeP1Space(source_mesh); + auto target_space = MakeP1Space(target_mesh); + + auto source_field = source_space.CreateField(); + pcms::test::SetField( + source_field, OMEGA_H_LAMBDA(pcms::Real x, pcms::Real y) { + return x * x + x * y + 0.5 * y * y; + }); + + // Reference: exact mesh-intersection quadrature projection. + auto reference_field = target_space.CreateField(); + pcms::OmegaHConservativeProjection reference_projection(source_space, + target_space); + reference_projection.Apply(source_field, reference_field); + const auto reference = + pcms::FlattenToRank1View(reference_field.GetDOFHolderDataHost()); + + const int samples_per_element = 64; + const uint64_t seed = 20240611; + + // Plain Monte Carlo projection. + pcms::OmegaHMonteCarloRHSIntegrator mc_integrator( + target_space, samples_per_element, pcms::MonteCarloSampling::UniformRandom, + seed); + auto evaluator = source_space.CreatePointEvaluator( + pcms::EvaluationRequest::FromCoordinates( + mc_integrator.GetIntegrationPoints().GetCoordinates())); + pcms::OmegaHMassIntegrator mass_integrator(target_space); + pcms::GalerkinProjectionSolver solver(mass_integrator, mc_integrator); + const auto mc_projected = solver.Solve(*evaluator, source_field); + const auto mc_projected_h = Omega_h::HostRead(mc_projected); + + // Control-variate projection with the same sampling configuration. + auto cv_field = target_space.CreateField(); + pcms::OmegaHControlVariateProjection cv_projection( + source_space, target_space, samples_per_element, + pcms::MonteCarloSampling::UniformRandom, seed); + cv_projection.Apply(source_field, cv_field); + const auto cv_values = + pcms::FlattenToRank1View(cv_field.GetDOFHolderDataHost()); + + double mc_error = 0.0; + double cv_error = 0.0; + for (Omega_h::LO i = 0; i < target_mesh.nverts(); ++i) { + mc_error = std::max(mc_error, std::abs(mc_projected_h[i] - reference[i])); + cv_error = std::max(cv_error, std::abs(cv_values[i] - reference[i])); + } + CAPTURE(mc_error, cv_error); + CHECK(cv_error < mc_error); +} From 2e58ed90c938fd993ab278823d99325208f32503 Mon Sep 17 00:00:00 2001 From: Jacob Merson Date: Sat, 11 Jul 2026 11:11:58 -0400 Subject: [PATCH 04/15] pythonapi: expose conservative transfer; flatten DOF fill for ND numpy buffers Squashed original commits: e7488c6, 4e68c80 --- .../example_conservative_transfer.py | 305 ++++++++++++++++++ src/pcms/pythonapi/CMakeLists.txt | 8 + src/pcms/pythonapi/bind_field_base.cpp | 50 ++- src/pcms/pythonapi/bind_transfer_field.cpp | 63 ++++ src/pcms/pythonapi/pythonapi.cpp | 24 ++ 5 files changed, 444 insertions(+), 6 deletions(-) create mode 100644 examples/python-api-example/example_conservative_transfer.py diff --git a/examples/python-api-example/example_conservative_transfer.py b/examples/python-api-example/example_conservative_transfer.py new file mode 100644 index 00000000..d6e6fd54 --- /dev/null +++ b/examples/python-api-example/example_conservative_transfer.py @@ -0,0 +1,305 @@ +#!/usr/bin/env python3 +""" +Conservative field transfer between two Omega_h meshes. + +This example: + 1. Builds a source Omega_h 2D simplex mesh, writes it to disk, and reads it + back (demonstrating the mesh read path). + 2. Evaluates an analytic function onto an order-1 (linear) Lagrange field on + the source mesh. + 3. Builds, writes, and reads back a *different* target Omega_h mesh that + carries its own order-1 Lagrange field. + 4. Transfers the source field onto the target field with both the Monte + Carlo (control-variate) and the mesh-intersection conservative projection + methods, and reports the L2 / max errors and the conservation error + (preservation of the integral) for each. + 5. Demonstrates the mesh-intersection projection across element orders as + well: P1 -> P0 (projection to cell averages) and P0 -> P1, confirming the + integral is conserved in both directions. + +The mesh-intersection conservative projection supports scalar Lagrange spaces +of order 0 (piecewise constant, one DOF per element) or order 1 (piecewise +linear, one DOF per vertex) on 2D simplex meshes; source and target orders may +differ. The Monte Carlo operator requires order-1 spaces. +""" + +import argparse +import os +import gc +import tempfile + +import numpy as np +import pcms + + +# Smooth analytic field sampled onto the source mesh. +def source_function(x, y): + return np.sin(np.pi * x) * np.sin(np.pi * y) + 0.5 * x + 0.25 * y + + +def build_and_roundtrip_mesh(world, lib, nx, ny, path): + """Build a 2D box mesh, write it to disk, and read it back.""" + mesh = pcms.build_box( + world, + pcms.Family.SIMPLEX, + 1.0, 1.0, 0.0, # unit square (z = 0 -> 2D) + nx, ny, 0, + False, + ) + pcms.write_mesh_binary(path, mesh) + del mesh + gc.collect() + return pcms.read_mesh_binary(path, lib) + + +def integral_of_field(mesh, field): + """Integral of an order-1 nodal field via element-averaged quadrature.""" + values = field.get_dof_holder_data() + elem_verts = mesh.ask_elem_verts().reshape(-1, 3) + areas = mesh.ask_sizes() # element areas for a 2D simplex mesh + elem_avg = values[elem_verts].mean(axis=1) + return float(np.sum(areas * elem_avg)) + + +def integral_of_p0_field(mesh, field): + """Integral of an order-0 field: sum of per-element value * element area.""" + values = field.get_dof_holder_data() # one value per element + areas = mesh.ask_sizes() + return float(np.sum(areas * values)) + + +def evaluate_function_onto_field(space, func): + """Create a field on `space` and fill it from `func` at the DOF holders.""" + field = space.create_field() + coords = field.get_dof_holder_coordinates() + field.set_dof_holder_data(func(coords[:, 0], coords[:, 1])) + return field + + +def report(name, target_field, target_mesh, reference_coords, + reference_values, source_integral, integrate=integral_of_field): + values = target_field.get_dof_holder_data() + err = values - reference_values + l2 = np.sqrt(np.mean(err ** 2)) + linf = np.max(np.abs(err)) + target_integral = integrate(target_mesh, target_field) + cons = abs(target_integral - source_integral) + print(f" {name}") + print(f" L2 error vs analytic : {l2:.3e}") + print(f" max error vs analytic : {linf:.3e}") + print(f" target integral : {target_integral:.6f}") + print(f" conservation error |dI| : {cons:.3e}") + + +def _triangulation(mesh, vertex_coords): + """matplotlib Triangulation from an Omega_h 2D simplex mesh.""" + import matplotlib.tri as mtri + tris = mesh.ask_elem_verts().reshape(-1, 3) + return mtri.Triangulation(vertex_coords[:, 0], vertex_coords[:, 1], tris) + + +def _draw_field(ax, triang, field, order, title, vmin, vmax): + values = field.get_dof_holder_data() + if order == 1: + # Nodal (per-vertex) values -> smooth (Gouraud) shading. + art = ax.tripcolor(triang, values, shading="gouraud", + vmin=vmin, vmax=vmax) + else: + # Cell (per-element) values -> flat shading. + art = ax.tripcolor(triang, facecolors=values, shading="flat", + vmin=vmin, vmax=vmax) + ax.set_title(title, fontsize=9) + ax.set_aspect("equal") + ax.set_xticks([]) + ax.set_yticks([]) + return art + + +def make_plots(out_path, source_mesh, target_mesh, source_vertex_coords, + target_vertex_coords, panels): + """Render before/after field panels to `out_path` (PNG). + + `panels` is a list of (mesh_key, field, order, title) where mesh_key is + "source" or "target". Returns True on success, False if matplotlib is + unavailable. + """ + try: + import matplotlib + matplotlib.use("Agg") # file output, no display needed + import matplotlib.pyplot as plt + except ImportError: + print("\n[plot] matplotlib not available; skipping plots. " + "Install matplotlib to enable --plot.") + return False + + triangs = { + "source": _triangulation(source_mesh, source_vertex_coords), + "target": _triangulation(target_mesh, target_vertex_coords), + } + + # Shared color scale across panels so fields are directly comparable. + all_vals = np.concatenate([f.get_dof_holder_data() for _, f, _, _ in panels]) + vmin, vmax = float(all_vals.min()), float(all_vals.max()) + + ncols = 3 + nrows = (len(panels) + ncols - 1) // ncols + fig, axes = plt.subplots(nrows, ncols, figsize=(4 * ncols, 3.6 * nrows), + squeeze=False) + art = None + for idx, (mesh_key, field, order, title) in enumerate(panels): + ax = axes[idx // ncols][idx % ncols] + art = _draw_field(ax, triangs[mesh_key], field, order, title, vmin, vmax) + for idx in range(len(panels), nrows * ncols): + axes[idx // ncols][idx % ncols].axis("off") + + fig.colorbar(art, ax=axes, shrink=0.85, label="field value") + fig.suptitle("Conservative projection: before / after", fontsize=12) + fig.savefig(out_path, dpi=130, bbox_inches="tight") + plt.close(fig) + print(f"\n[plot] wrote {out_path}") + return True + + +def main(plot_path=None): + lib = pcms.OmegaHLibrary() + world = lib.world() + + work_dir = tempfile.mkdtemp(prefix="pcms_conservative_transfer_") + try: + print("=" * 64) + print("Conservative field transfer: MonteCarlo and mesh intersection") + print("=" * 64) + + # 1. Source mesh (read from disk) + linear Lagrange field. + source_path = os.path.join(work_dir, "source.osh") + source_mesh = build_and_roundtrip_mesh(world, lib, 12, 12, source_path) + print(f"\nSource mesh : {source_mesh.nverts()} verts, " + f"{source_mesh.nelems()} elems (read from {source_path})") + + # The conservative-projection operators require the native Omega_h + # Lagrange backend (not the default MeshFields layout). + source_space = pcms.LagrangeFunctionSpace.from_mesh( + source_mesh, 1, 1, pcms.CoordinateSystem.Cartesian, + backend=pcms.LagrangeFunctionSpace.Backend.OmegaH, + ) + source_field = evaluate_function_onto_field(source_space, + source_function) + source_integral = integral_of_field(source_mesh, source_field) + print(f"Source field: order-1 Lagrange, " + f"integral = {source_integral:.6f}") + + # 2. Target mesh (read from disk) with its own linear Lagrange field. + target_path = os.path.join(work_dir, "target.osh") + target_mesh = build_and_roundtrip_mesh(world, lib, 17, 17, target_path) + print(f"\nTarget mesh : {target_mesh.nverts()} verts, " + f"{target_mesh.nelems()} elems (read from {target_path})") + + target_space = pcms.LagrangeFunctionSpace.from_mesh( + target_mesh, 1, 1, pcms.CoordinateSystem.Cartesian, + backend=pcms.LagrangeFunctionSpace.Backend.OmegaH, + ) + + # Analytic reference values at the target DOF holders. + target_template = target_space.create_field() + target_coords = target_template.get_dof_holder_coordinates() + reference_values = source_function(target_coords[:, 0], + target_coords[:, 1]) + + print("\nTransferring source field to target field:") + + # 3a. Mesh-intersection conservative projection. + mi_target = target_space.create_field() + mesh_intersection = pcms.OmegaHConservativeProjection( + source_space, target_space + ) + mesh_intersection.apply(source_field, mi_target) + report("mesh intersection (exact quadrature)", mi_target, target_mesh, + target_coords, reference_values, source_integral) + + # 3b. Monte Carlo (control-variate) conservative projection. + mc_target = target_space.create_field() + monte_carlo = pcms.OmegaHControlVariateProjection( + source_space, + target_space, + samples_per_element=64, + sampling=pcms.MonteCarloSampling.UniformRandom, + seed=8675309, + ) + monte_carlo.apply(source_field, mc_target) + report("Monte Carlo (control variate, UniformRandom)", mc_target, target_mesh, + target_coords, reference_values, source_integral) + + # 4. Mixed-order mesh-intersection projections. The same operator works + # for any combination of order-0 (piecewise constant) and order-1 + # (piecewise linear) Lagrange spaces; the integral is conserved in + # every case. + print("\nMixed-order mesh-intersection projections:") + + # 4a. P1 source -> P0 target (project the linear field to cell averages). + target_p0_space = pcms.LagrangeFunctionSpace.from_mesh( + target_mesh, 0, 1, pcms.CoordinateSystem.Cartesian, + backend=pcms.LagrangeFunctionSpace.Backend.OmegaH, + ) + p1_to_p0 = target_p0_space.create_field() + # Analytic reference at the P0 DOF holders (element centroids). + p0_coords = p1_to_p0.get_dof_holder_coordinates() + p0_reference = source_function(p0_coords[:, 0], p0_coords[:, 1]) + proj_p1_p0 = pcms.OmegaHConservativeProjection( + source_space, target_p0_space + ) + proj_p1_p0.apply(source_field, p1_to_p0) + report("P1 -> P0 (cell averages)", p1_to_p0, target_mesh, + p0_coords, p0_reference, source_integral, + integrate=integral_of_p0_field) + + # 4b. P0 source -> P1 target. Build a P0 source field by sampling the + # analytic function at the source element centroids. + source_p0_space = pcms.LagrangeFunctionSpace.from_mesh( + source_mesh, 0, 1, pcms.CoordinateSystem.Cartesian, + backend=pcms.LagrangeFunctionSpace.Backend.OmegaH, + ) + source_p0_field = evaluate_function_onto_field(source_p0_space, + source_function) + source_p0_integral = integral_of_p0_field(source_mesh, source_p0_field) + p0_to_p1 = target_space.create_field() + proj_p0_p1 = pcms.OmegaHConservativeProjection( + source_p0_space, target_space + ) + proj_p0_p1.apply(source_p0_field, p0_to_p1) + report("P0 -> P1", p0_to_p1, target_mesh, + target_coords, reference_values, source_p0_integral) + + print("\n✓ Field transfer completed for P1<->P1, P1->P0, and P0->P1.") + + # 5. Optional before/after visualization. + if plot_path is not None: + source_vertex_coords = source_field.get_dof_holder_coordinates() + make_plots( + plot_path, source_mesh, target_mesh, + source_vertex_coords, target_coords, + panels=[ + ("source", source_field, 1, "source (P1)"), + ("target", mi_target, 1, "P1 -> P1 (mesh intersection)"), + ("target", p1_to_p0, 0, "P1 -> P0 (cell averages)"), + ("source", source_p0_field, 0, "source (P0)"), + ("target", p0_to_p1, 1, "P0 -> P1"), + ], + ) + + finally: + import shutil + shutil.rmtree(work_dir, ignore_errors=True) + #del world + #del lib + + +if __name__ == "__main__": + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument( + "--plot", nargs="?", const="conservative_transfer.png", default=None, + metavar="PATH", + help="render before/after field panels to PATH (default " + "conservative_transfer.png) using matplotlib, if installed", + ) + args = parser.parse_args() + main(plot_path=args.plot) diff --git a/src/pcms/pythonapi/CMakeLists.txt b/src/pcms/pythonapi/CMakeLists.txt index bb2c316b..f2c569cd 100644 --- a/src/pcms/pythonapi/CMakeLists.txt +++ b/src/pcms/pythonapi/CMakeLists.txt @@ -17,6 +17,14 @@ pybind11_add_module(pcms target_compile_features(pcms PUBLIC cxx_std_20) target_link_libraries(pcms PRIVATE Kokkos::kokkos pybind11::module MPI::MPI_C pcms::core pcms::transfer) +# The conservative/Monte Carlo projection bindings include transfer headers +# that transitively pull in PETSc (petscvec.h). PETSc is linked PRIVATE into +# pcms_transfer, so its include directories do not propagate; link it here so +# the bindings compile when PETSc + MeshFields are enabled. +if(PCMS_ENABLE_PETSC AND PCMS_ENABLE_MESHFIELDS) + target_link_libraries(pcms PRIVATE PETSc::PETSc meshfields::meshfields) +endif() + # if we are building as a python package from the pyproject.toml # scikit-build-core will set this variable if(SKBUILD) diff --git a/src/pcms/pythonapi/bind_field_base.cpp b/src/pcms/pythonapi/bind_field_base.cpp index ee4bff09..77be9871 100644 --- a/src/pcms/pythonapi/bind_field_base.cpp +++ b/src/pcms/pythonapi/bind_field_base.cpp @@ -8,6 +8,7 @@ #include "pcms/field/field_layout.h" #include "pcms/field/field_metadata.h" #include "pcms/field/function_space.h" +#include "pcms/discretization/discretization/omega_h.hpp" #include "pcms/field/data/simple.h" #include "pcms/field/layout/uniform_grid.h" #include "pcms/field/uniform_grid_binary_field.h" @@ -223,7 +224,12 @@ void bind_create_field_module(py::module& m) .def( "set_dof_holder_data", - [](Field& self, py::array_t arr) { + // c_style|forcecast makes pybind materialize a contiguous copy of the + // input, so non-contiguous arrays (e.g. a column slice like data[:, i]) + // are read correctly rather than by walking flat memory with the wrong + // stride. + [](Field& self, + py::array_t arr) { auto buf = arr.request(); if (buf.ndim != 1) { throw std::runtime_error("DOF holder data must be a 1D array"); @@ -293,19 +299,51 @@ void bind_create_field_module(py::module& m) py::arg("request"), "Create a reusable point evaluator from an EvaluationRequest.") .def("get_coordinate_system", &FunctionSpace::GetCoordinateSystem, - "Get the coordinate system for this function space"); + "Get the coordinate system for this function space") + .def( + "mesh", + [](const FunctionSpace& self) -> Omega_h::Mesh& { + auto disc = std::dynamic_pointer_cast( + self.GetDiscretization()); + if (!disc) { + throw std::runtime_error( + "FunctionSpace::mesh: this function space is not backed by an " + "Omega_h mesh"); + } + return disc->GetMesh(); + }, + py::return_value_policy::reference, + "The Omega_h mesh this function space is defined on (Omega_h backend " + "only)."); // Bind LagrangeFunctionSpace as a concrete FunctionSpace subtype. - py::class_(m, "LagrangeFunctionSpace") + py::class_ lagrange_space( + m, "LagrangeFunctionSpace"); + + // Backend selecting the underlying field layout. The conservative-projection + // transfer operators require the native Omega_h backend. + py::enum_(lagrange_space, "Backend") + .value("MeshFields", LagrangeFunctionSpace::Backend::MeshFields, + "MeshFields-backed layout (default when MeshFields is enabled)") + .value("OmegaH", LagrangeFunctionSpace::Backend::OmegaH, + "Native Omega_h Lagrange layout (required by the conservative and " + "Monte Carlo projection transfer operators)") + .export_values(); + + lagrange_space .def_static( "from_mesh", [](Omega_h::Mesh& mesh, int order, int num_components, - CoordinateSystem coordinate_system) { - return LagrangeFunctionSpace::FromMesh(mesh, order, num_components, - coordinate_system); + CoordinateSystem coordinate_system, std::string global_id_name, + LagrangeFunctionSpace::Backend backend) { + return LagrangeFunctionSpace::FromMesh( + mesh, order, num_components, coordinate_system, + std::move(global_id_name), backend); }, py::arg("mesh"), py::arg("order"), py::arg("num_components") = 1, py::arg("coordinate_system") = CoordinateSystem::Cartesian, + py::arg("global_id_name") = "global", + py::arg("backend") = LagrangeFunctionSpace::DefaultBackend, "Create a LagrangeFunctionSpace from an Omega_h mesh") .def_static( diff --git a/src/pcms/pythonapi/bind_transfer_field.cpp b/src/pcms/pythonapi/bind_transfer_field.cpp index 599b31f3..cc599bb0 100644 --- a/src/pcms/pythonapi/bind_transfer_field.cpp +++ b/src/pcms/pythonapi/bind_transfer_field.cpp @@ -1,6 +1,7 @@ #include #include #include +#include "pcms/configuration.h" #include "pcms/transfer/copy.h" #include "pcms/field/field.h" #include "pcms/field/function_space.h" @@ -9,6 +10,12 @@ #include "pcms/field/out_of_bounds_policy.h" #include "numpy_array_transform.h" #include "pcms/utility/types.h" +#if defined(PCMS_ENABLE_PETSC) && defined(PCMS_ENABLE_MESHFIELDS) +#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" +#include +#endif namespace py = pybind11; @@ -70,6 +77,62 @@ void bind_transfer_field_module(py::module& m) Field& target) { self.Apply(source, target); }, py::arg("source"), py::arg("target"), "Copy source field data to target field (same layout required)."); + +#if defined(PCMS_ENABLE_PETSC) && defined(PCMS_ENABLE_MESHFIELDS) + // 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_(m, "OmegaHConservativeProjection") + .def(py::init([](const FunctionSpace& source_space, + const FunctionSpace& target_space) { + return std::make_unique(source_space, + target_space); + }), + py::arg("source_space"), py::arg("target_space"), + "Construct a mesh-intersection conservative projection. Mesh " + "intersection, quadrature setup, and mass-matrix factorization happen " + "here and are cached; call apply() repeatedly.") + .def( + "apply", + [](const OmegaHConservativeProjection& self, const Field& source, + Field& target) { self.Apply(source, target); }, + py::arg("source"), py::arg("target"), + "Conservatively project the source field onto the target space using " + "exact mesh-intersection quadrature."); + + // Sampling strategy for the Monte Carlo RHS integrator. + py::enum_(m, "MonteCarloSampling") + .value("UniformRandom", MonteCarloSampling::UniformRandom, + "Independent pseudo-random uniform draws per element") + .export_values(); + + // Variance-reduced Monte Carlo Galerkin projection. The source field is + // first interpolated onto the target space as a control variate; only the + // residual is integrated stochastically, so sampling noise is small (and + // exactly zero when the source already lives in the target space). + py::class_(m, + "OmegaHControlVariateProjection") + .def(py::init([](const FunctionSpace& source_space, + const FunctionSpace& target_space, int samples_per_element, + MonteCarloSampling sampling, std::uint64_t seed) { + return std::make_unique( + source_space, target_space, samples_per_element, sampling, seed); + }), + 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), + "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.") + .def( + "apply", + [](const OmegaHControlVariateProjection& self, const Field& source, + Field& target) { self.Apply(source, target); }, + py::arg("source"), py::arg("target"), + "Project the source field onto the target space using Monte Carlo " + "integration with control-variate variance reduction."); +#endif } } // namespace pcms diff --git a/src/pcms/pythonapi/pythonapi.cpp b/src/pcms/pythonapi/pythonapi.cpp index 3e9e253e..344ae87c 100644 --- a/src/pcms/pythonapi/pythonapi.cpp +++ b/src/pcms/pythonapi/pythonapi.cpp @@ -1,5 +1,9 @@ #include +#include "pcms/configuration.h" +#ifdef PCMS_ENABLE_PETSC +#include +#endif namespace py = pybind11; @@ -33,6 +37,26 @@ void bind_mesh_utilities_module(py::module& m); } // namespace pcms PYBIND11_MODULE(pcms, m) { +#ifdef PCMS_ENABLE_PETSC + // The conservative/Monte Carlo projection solvers build PETSc objects on + // PETSC_COMM_WORLD, so PETSc must be initialized before any of them are + // constructed. Do it once at import time (PetscInitialize brings up MPI if it + // is not already running) so callers never manage PETSc state by hand. + // + // PETSc is intentionally NOT finalized via atexit: PETSc's Kokkos-backed + // objects must be torn down before Kokkos::finalize, but Kokkos is owned by + // OmegaHLibrary and an atexit hook would run after that library is destroyed. + // Leaving PETSc initialized until the process exits avoids that ordering trap + // and is harmless (the OS reclaims everything on exit). + { + PetscBool petsc_initialized = PETSC_FALSE; + PetscInitialized(&petsc_initialized); + if (!petsc_initialized) { + PetscInitializeNoArguments(); + } + } +#endif + // Bind fundamental types first (coordinate systems, etc.) pcms::bind_coordinate_system_module(m); pcms::bind_coordinate_module(m); From 28aac5195f2c1e147c98c028ee4cf1904c683f8a Mon Sep 17 00:00:00 2001 From: Jacob Merson Date: Sat, 11 Jul 2026 11:11:58 -0400 Subject: [PATCH 05/15] transfer: ad-hoc adjacency-search snapping fixes (TEMPORARY; remove after #308) Squashed original commit: 75a3045 --- src/pcms/localization/adj_search.cpp | 53 +++++++++++++++++++--------- 1 file changed, 37 insertions(+), 16 deletions(-) diff --git a/src/pcms/localization/adj_search.cpp b/src/pcms/localization/adj_search.cpp index 6f88a2ce..20e4676c 100644 --- a/src/pcms/localization/adj_search.cpp +++ b/src/pcms/localization/adj_search.cpp @@ -93,13 +93,11 @@ static Omega_h::Write locate_target_cells( Omega_h::parallel_for( nvertices_target, OMEGA_H_LAMBDA(const Omega_h::LO i) { auto source_cell_id = results(i).element_id; - if (source_cell_id < 0) - source_cell_id = Kokkos::abs(source_cell_id); - OMEGA_H_CHECK_PRINTF( - source_cell_id >= 0, - "ERROR: Source cell id not found for target %d (%f,%f)\n", i, - target_points(i, 0), target_points(i, 1)); - source_cell_ids[i] = source_cell_id; + // A negative id means the point is outside the source mesh. Until the + // localization returns the negated *closest* element for exterior + // points (PR #319), abs() is not a valid face index, so flag it as + // invalid (-1) and let the neighbor search skip the target. + source_cell_ids[i] = source_cell_id >= 0 ? source_cell_id : -1; }); } else if (dim == 3) { Kokkos::View target_points("target_points", @@ -117,12 +115,9 @@ static Omega_h::Write locate_target_cells( Omega_h::parallel_for( nvertices_target, OMEGA_H_LAMBDA(const Omega_h::LO i) { auto source_cell_id = results(i).element_id; - if (source_cell_id < 0) - source_cell_id = Kokkos::abs(source_cell_id); - OMEGA_H_CHECK_PRINTF(source_cell_id >= 0, - "ERROR: Source cell id not found for target %d\n", - i); - source_cell_ids[i] = source_cell_id; + // See the 2D branch: exterior points (negative id) are flagged invalid + // until PR #319 makes the negated-closest-element id usable. + source_cell_ids[i] = source_cell_id >= 0 ? source_cell_id : -1; }); } else { throw pcms_error("Unsupported dimension in locate_target_cells"); @@ -208,9 +203,12 @@ static void adjBasedSearchFromPoints( Omega_h::Real cutoffDistance = radii2[id]; Omega_h::LO source_cell_id = source_cell_ids[id]; - OMEGA_H_CHECK_PRINTF(source_cell_id >= 0, - "ERROR: Source cell id not found for target %d\n", - id); + // Targets outside the source mesh are flagged with an invalid cell id + // (see locate_target_cells). Skip them: no supports, no memory access. + if (source_cell_id < 0) { + nSupports[id] = 0; + return; + } const Omega_h::LO num_verts_in_dim = dim + 1; Omega_h::LO start_ptr = source_cell_id * num_verts_in_dim; @@ -382,6 +380,10 @@ SupportResults searchNeighbors(Omega_h::Mesh& source_mesh, true); } else { pcms::printInfo("INFO: Adaptive radius search... \n"); + // Cap the adaptation so targets that can never reach the required support + // count (e.g. points outside the source mesh, which are skipped and always + // report 0 supports) do not spin the loop forever. + const int max_radius_adjust_loops = 50; int r_adjust_loop = 0; while (true) { nSupports = Omega_h::Write( @@ -410,6 +412,14 @@ SupportResults searchNeighbors(Omega_h::Mesh& source_mesh, break; } + if (r_adjust_loop >= max_radius_adjust_loops) { + pcms::printInfo( + "WARNING: radius adaptation hit the %d-loop cap; some targets remain " + "outside the [%d, %d] support range\n", + max_radius_adjust_loops, min_req_support, 3 * min_req_support); + break; + } + adapt_radii(min_req_support, max_allowed_support, nvertices_target, radii2, nSupports); } @@ -471,6 +481,9 @@ SupportResults searchNeighbors(Omega_h::Mesh& mesh, search.adjBasedSearch(supports_ptr, nSupports, supports_idx, radii2, true); } else { pcms::printInfo("INFO: Adaptive radius search... \n"); + // Cap the adaptation so targets that can never reach the required support + // count do not spin the loop forever (see searchNeighborsFromPoints). + const int max_radius_adjust_loops = 50; int r_adjust_loop = 0; while (true) { // until the number of minimum support is met const auto max_radius = Omega_h::get_max(Omega_h::read(radii2)); @@ -498,6 +511,14 @@ SupportResults searchNeighbors(Omega_h::Mesh& mesh, break; } + if (r_adjust_loop >= max_radius_adjust_loops) { + pcms::printInfo( + "WARNING: radius adaptation hit the %d-loop cap; some targets remain " + "outside the [%d, %d] support range\n", + max_radius_adjust_loops, min_support, 3 * min_support); + break; + } + adapt_radii(min_support, 3 * min_support, radii2.size(), radii2, nSupports); } From 328e67430fcfc084737e577f73feaa7bb25580a7 Mon Sep 17 00:00:00 2001 From: Jacob Merson Date: Sat, 11 Jul 2026 14:57:10 -0400 Subject: [PATCH 06/15] fix petsc include guards in python bindings 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 --- src/pcms/pythonapi/pythonapi.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/pcms/pythonapi/pythonapi.cpp b/src/pcms/pythonapi/pythonapi.cpp index 344ae87c..e3491377 100644 --- a/src/pcms/pythonapi/pythonapi.cpp +++ b/src/pcms/pythonapi/pythonapi.cpp @@ -1,7 +1,7 @@ #include #include "pcms/configuration.h" -#ifdef PCMS_ENABLE_PETSC +#if defined(PCMS_ENABLE_PETSC) && defined(PCMS_ENABLE_MESHFIELDS) #include #endif @@ -37,7 +37,7 @@ void bind_mesh_utilities_module(py::module& m); } // namespace pcms PYBIND11_MODULE(pcms, m) { -#ifdef PCMS_ENABLE_PETSC +#if defined(PCMS_ENABLE_PETSC) && defined(PCMS_ENABLE_MESHFIELDS) // The conservative/Monte Carlo projection solvers build PETSc objects on // PETSC_COMM_WORLD, so PETSc must be initialized before any of them are // constructed. Do it once at import time (PetscInitialize brings up MPI if it From a69f406fa5ad4b97675b9764fa8e3b9aa270b831 Mon Sep 17 00:00:00 2001 From: Jacob Merson Date: Sat, 11 Jul 2026 20:03:53 -0400 Subject: [PATCH 07/15] transfer: index DOF holders by compact GID permutation 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. --- src/pcms/field/CMakeLists.txt | 3 + src/pcms/field/field_layout.cpp | 38 ++++++++ src/pcms/field/field_layout.h | 26 ++++- src/pcms/field/gid_permutation.cpp | 70 +++++++++++++ src/pcms/field/gid_permutation.hpp | 45 +++++++++ src/pcms/field/layout/omega_h_lagrange.cpp | 4 + .../omega_h_conservative_projection.cpp | 24 ++++- .../omega_h_conservative_projection.hpp | 2 + .../omega_h_control_variate_projection.cpp | 19 ++-- .../omega_h_form_integrator_utils.hpp | 22 +++-- .../omega_h_intersection_rhs_integrator.cpp | 9 +- src/pcms/transfer/omega_h_mass_integrator.cpp | 11 +-- .../transfer/omega_h_mc_rhs_integrator.cpp | 13 +-- .../test_mesh_intersection_field_transfer.cpp | 97 ++++++++++++++++++- 14 files changed, 333 insertions(+), 50 deletions(-) create mode 100644 src/pcms/field/field_layout.cpp create mode 100644 src/pcms/field/gid_permutation.cpp create mode 100644 src/pcms/field/gid_permutation.hpp diff --git a/src/pcms/field/CMakeLists.txt b/src/pcms/field/CMakeLists.txt index b37865be..fc15c3a6 100644 --- a/src/pcms/field/CMakeLists.txt +++ b/src/pcms/field/CMakeLists.txt @@ -5,6 +5,7 @@ set(PCMS_FIELD_HEADERS coordinate_system.h element_dispatch.h evaluation_request.h + gid_permutation.hpp field_layout.h field_metadata.h layout/empty.h @@ -18,6 +19,8 @@ set(PCMS_FIELD_HEADERS data/point_cloud.h ) set(PCMS_FIELD_SOURCES + field_layout.cpp + gid_permutation.cpp function_space/polynomial_reconstruction.cpp layout/empty.cpp layout/point_cloud.cpp diff --git a/src/pcms/field/field_layout.cpp b/src/pcms/field/field_layout.cpp new file mode 100644 index 00000000..6e785f9e --- /dev/null +++ b/src/pcms/field/field_layout.cpp @@ -0,0 +1,38 @@ +#include "pcms/field/field_layout.h" +#include "pcms/field/gid_permutation.hpp" +#include "pcms/utility/assert.h" + +namespace pcms +{ + +void FieldLayout::BuildGlobalToLocalPermutation() +{ + const auto permutation = + BuildSortedGidPermutation(GetGidsHost(), GetEntOffsets()); + global_to_local_host_ = Kokkos::View( + "global_to_local_host", permutation.size()); + for (size_t i = 0; i < permutation.size(); ++i) + global_to_local_host_(i) = permutation[i]; + + global_to_local_ = + Kokkos::View("global_to_local", permutation.size()); + Kokkos::deep_copy(global_to_local_, global_to_local_host_); +} + +Kokkos::View +FieldLayout::GetGlobalToLocalPermutationHost() const +{ + // A non-empty view means BuildGlobalToLocalPermutation() ran in the derived + // layout's constructor; an empty one means the layout never built it. + PCMS_ALWAYS_ASSERT(global_to_local_host_.size() > 0); + return global_to_local_host_; +} + +Kokkos::View +FieldLayout::GetGlobalToLocalPermutation() const +{ + PCMS_ALWAYS_ASSERT(global_to_local_.size() > 0); + return global_to_local_; +} + +} // namespace pcms diff --git a/src/pcms/field/field_layout.h b/src/pcms/field/field_layout.h index 03aef658..d01a13a0 100644 --- a/src/pcms/field/field_layout.h +++ b/src/pcms/field/field_layout.h @@ -40,14 +40,20 @@ class FieldLayout virtual Rank1View GetOwnedHost() const = 0; virtual GlobalIDView GetGidsHost() const = 0; + // Maps each local DOF holder to its contiguous active index, ordered by GID + // within each entity block. Components are not included in the permutation. + // For local GIDs [102, 7, 41, 19], the permutation is [3, 0, 2, 1]. Thus + // holder i, component c is read from values(permutation(i), c), or from + // flat_values[permutation(i) * num_components + c]. + Kokkos::View GetGlobalToLocalPermutationHost() + const; + Kokkos::View GetGlobalToLocalPermutation() + const; + // returns true if the field layout is distributed // if the field layout is distributed, the owned and global dofs are the same [[nodiscard]] virtual bool IsDistributed() const = 0; - // This class should construct the permutation arrays that are needed - // for serialization / deserialization - // - virtual EntOffsetsArray GetEntOffsets() const = 0; virtual CoordinateView GetDOFHolderCoordinates() const = 0; @@ -61,6 +67,18 @@ class FieldLayout GetDOFHolderClassificationIdsHost() const = 0; virtual ~FieldLayout() noexcept = default; + +protected: + // Builds the global-to-local permutation from GetGidsHost()/GetEntOffsets(). + // Derived layouts that support GetGlobalToLocalPermutation() must call this + // from their constructor, once their GIDs and entity offsets are available; + // building eagerly keeps first access free of the data race that lazy + // construction of the mutable cache would introduce. + void BuildGlobalToLocalPermutation(); + +private: + Kokkos::View global_to_local_host_; + Kokkos::View global_to_local_; }; } // namespace pcms diff --git a/src/pcms/field/gid_permutation.cpp b/src/pcms/field/gid_permutation.cpp new file mode 100644 index 00000000..475216f8 --- /dev/null +++ b/src/pcms/field/gid_permutation.cpp @@ -0,0 +1,70 @@ +#include "pcms/field/gid_permutation.hpp" +#include "pcms/utility/assert.h" +#include + +namespace pcms +{ + +namespace +{ + +// Sentinel permutation entry for a local holder that does not participate in an +// exchange (a non-owned/ghost holder, or an owned holder outside the coupled +// overlap region whose GID is never transmitted). It is deliberately negative +// so serialization/deserialization can skip it, and any accidental use as a +// buffer index is an obvious out-of-bounds access rather than a silent read or +// write of holder 0. +constexpr LO UnmappedPermutationIndex = -1; + +} // namespace + +std::vector MakeUnmappedPermutation(std::size_t size) +{ + return std::vector(size, UnmappedPermutationIndex); +} + +void AppendGidPermutation(GlobalIDView local_gids, + const EntOffsetsArray& ent_offsets, + const GidToIndexMaps& gid_to_index, + std::vector& permutation, bool allow_missing) +{ + for (size_t e = 0; e < ent_offsets.size() - 1; ++e) { + for (size_t i = ent_offsets[e]; i < ent_offsets[e + 1]; ++i) { + const auto entry = gid_to_index[e].find(local_gids(i)); + if (entry != gid_to_index[e].end()) { + permutation.push_back(entry->second); + continue; + } + PCMS_ALWAYS_ASSERT(allow_missing); + permutation.push_back(UnmappedPermutationIndex); + } + } +} + +std::vector BuildSortedGidPermutation(GlobalIDView gids, + const EntOffsetsArray& ent_offsets) +{ + GidToIndexMaps gid_to_index; + LO next_index = 0; + for (size_t e = 0; e < ent_offsets.size() - 1; ++e) { + std::vector sorted_gids; + sorted_gids.reserve(ent_offsets[e + 1] - ent_offsets[e]); + for (size_t i = ent_offsets[e]; i < ent_offsets[e + 1]; ++i) { + sorted_gids.push_back(gids(i)); + } + std::sort(sorted_gids.begin(), sorted_gids.end()); + PCMS_ALWAYS_ASSERT( + std::adjacent_find(sorted_gids.begin(), sorted_gids.end()) == + sorted_gids.end()); + for (const GO gid : sorted_gids) { + gid_to_index[e].emplace(gid, next_index++); + } + } + std::vector permutation; + permutation.reserve(gids.size()); + AppendGidPermutation(gids, ent_offsets, gid_to_index, permutation, + /*allow_missing=*/false); + return permutation; +} + +} // namespace pcms diff --git a/src/pcms/field/gid_permutation.hpp b/src/pcms/field/gid_permutation.hpp new file mode 100644 index 00000000..0118fd31 --- /dev/null +++ b/src/pcms/field/gid_permutation.hpp @@ -0,0 +1,45 @@ +#ifndef PCMS_FIELD_GID_PERMUTATION_HPP +#define PCMS_FIELD_GID_PERMUTATION_HPP + +#include "pcms/field/field_layout.h" +#include +#include +#include +#include + +namespace pcms +{ + +// Maps a GID to its position in an exchange buffer, grouped by entity-dimension +// block: index e holds the GIDs classified to mesh entity dimension e. +using GidToIndexMaps = std::array, ent_offsets_len - 1>; + +// Returns a permutation of `size` entries with every local holder initially +// unmapped. Unmapped entries hold a negative sentinel so that serialization +// skips them and any accidental use as a buffer index is an obvious +// out-of-bounds access rather than a silent read or write of holder 0. Callers +// overwrite the entries for holders that participate in the exchange. +std::vector MakeUnmappedPermutation(std::size_t size); + +// Appends, in entity-offset (local holder) order, the position each local GID +// maps to according to gid_to_index. When allow_missing is true, a GID absent +// from gid_to_index is left unmapped (negative sentinel): a local layout may +// contain ghost holders and holders outside the overlap region, which are never +// exchanged and which (de)serialization skips. When allow_missing is false -- +// e.g. the sorted permutation below, whose maps are built from these very GIDs +// -- a missing GID is a hard error. +void AppendGidPermutation(GlobalIDView local_gids, + const EntOffsetsArray& ent_offsets, + const GidToIndexMaps& gid_to_index, + std::vector& permutation, + bool allow_missing = false); + +// Builds a permutation that maps each local holder to a contiguous index +// ordered by GID within each entity block. Every GID must be unique within its +// block. +std::vector BuildSortedGidPermutation(GlobalIDView gids, + const EntOffsetsArray& ent_offsets); + +} // namespace pcms + +#endif // PCMS_FIELD_GID_PERMUTATION_HPP diff --git a/src/pcms/field/layout/omega_h_lagrange.cpp b/src/pcms/field/layout/omega_h_lagrange.cpp index 1b2f14a4..a325140f 100644 --- a/src/pcms/field/layout/omega_h_lagrange.cpp +++ b/src/pcms/field/layout/omega_h_lagrange.cpp @@ -146,6 +146,8 @@ OmegaHLagrangeLayout::OmegaHLagrangeLayout(Omega_h::Mesh& mesh, int order, classification_ids_host_(i) = static_cast(class_ids_host[i]); } discretization_ = std::make_shared(mesh_); + + BuildGlobalToLocalPermutation(); } OmegaHLagrangeLayout::OmegaHLagrangeLayout( @@ -189,6 +191,8 @@ OmegaHLagrangeLayout::OmegaHLagrangeLayout( classification_ids_host_(i) = static_cast(class_ids_host[i]); } discretization_ = std::make_shared(mesh_); + + BuildGlobalToLocalPermutation(); } std::shared_ptr OmegaHLagrangeLayout::GetDiscretization() diff --git a/src/pcms/transfer/omega_h_conservative_projection.cpp b/src/pcms/transfer/omega_h_conservative_projection.cpp index f5d911ec..aac5918c 100644 --- a/src/pcms/transfer/omega_h_conservative_projection.cpp +++ b/src/pcms/transfer/omega_h_conservative_projection.cpp @@ -60,6 +60,9 @@ OmegaHConservativeProjection::OmegaHConservativeProjection( target_space.GetCoordinateSystem()); solver_ = std::make_unique(mass_integrator, *rhs_integrator_); + target_values_ = Kokkos::View( + "conservative_projection_target_values", + target_layout_->GetNumOwnedDofHolder(), target_layout_->GetNumComponents()); } // Defined here so that GalerkinProjectionSolver (forward-declared in the @@ -71,11 +74,22 @@ void OmegaHConservativeProjection::Apply(const Field& source, { CheckApplyCompatible(source, target, *source_layout_, *target_layout_); - const auto target_values = solver_->Solve(*evaluator_, source); - auto target_values_h = Omega_h::HostRead(target_values); - target.SetDOFHolderDataHost(Rank2View( - target_values_h.data(), target_layout_->GetNumOwnedDofHolder(), - target_layout_->GetNumComponents())); + const auto solution = solver_->Solve(*evaluator_, source); + const auto global_to_local = target_layout_->GetGlobalToLocalPermutation(); + const int num_dof_holders = target_layout_->GetNumOwnedDofHolder(); + const int num_components = target_layout_->GetNumComponents(); + auto target_values = target_values_; + Kokkos::parallel_for( + "conservative_projection_scatter_solution", + Kokkos::RangePolicy(0, num_dof_holders), + KOKKOS_LAMBDA(int i) { + for (int c = 0; c < num_components; ++c) { + target_values(i, c) = solution[global_to_local(i) * num_components + c]; + } + }); + Kokkos::fence(); + + target.SetDOFHolderData(MakeConstRank2View(target_values_)); } } // namespace pcms diff --git a/src/pcms/transfer/omega_h_conservative_projection.hpp b/src/pcms/transfer/omega_h_conservative_projection.hpp index 114fc2a0..bf4fcb93 100644 --- a/src/pcms/transfer/omega_h_conservative_projection.hpp +++ b/src/pcms/transfer/omega_h_conservative_projection.hpp @@ -6,6 +6,7 @@ #include "pcms/field/point_evaluator.h" #include "pcms/transfer/omega_h_intersection_rhs_integrator.hpp" #include "pcms/transfer/transfer_operator.hpp" +#include #include namespace pcms @@ -45,6 +46,7 @@ class OmegaHConservativeProjection : public TransferOperator std::unique_ptr rhs_integrator_; std::unique_ptr> evaluator_; std::unique_ptr solver_; + mutable Kokkos::View target_values_; }; } // namespace pcms diff --git a/src/pcms/transfer/omega_h_control_variate_projection.cpp b/src/pcms/transfer/omega_h_control_variate_projection.cpp index c22b3f71..f30f1bc0 100644 --- a/src/pcms/transfer/omega_h_control_variate_projection.cpp +++ b/src/pcms/transfer/omega_h_control_variate_projection.cpp @@ -64,26 +64,21 @@ void OmegaHControlVariateProjection::Apply(const Field& source, // 3. Project the residual: M * delta = r. const auto delta = solver_->Solve(MakeConstRank2View(residual)); - // 4. x = g + delta. delta is indexed by global id (PETSc row), the DOF - // holder by local vertex id. + // 4. x = g + delta. delta is indexed by active PETSc row, while the field is + // indexed by local DOF holder. const auto g_nodal = control_variate_.GetDOFHolderData(); - auto g_view = Kokkos::View>( - g_nodal.data_handle(), g_nodal.extent(0)); - const auto device_gids = target_layout_->GetGids(); - const GO* gids_ptr = device_gids.data_handle(); - const int nverts = static_cast(g_view.extent(0)); + const auto global_to_local = target_layout_->GetGlobalToLocalPermutation(); + const int nverts = static_cast(g_nodal.extent(0)); - Kokkos::View result("cv_result", nverts); + Kokkos::View result("cv_result", nverts, 1); Kokkos::parallel_for( "cv_add_correction", Kokkos::RangePolicy(0, nverts), KOKKOS_LAMBDA(int i) { - result(i) = g_view(i) + delta[static_cast(gids_ptr[i])]; + result(i, 0) = g_nodal(i, 0) + delta[global_to_local(i)]; }); Kokkos::fence(); - target.SetDOFHolderData( - Rank2View(result.data(), nverts, 1)); + target.SetDOFHolderData(MakeConstRank2View(result)); } } // namespace pcms diff --git a/src/pcms/transfer/omega_h_form_integrator_utils.hpp b/src/pcms/transfer/omega_h_form_integrator_utils.hpp index 3294e45e..23f5e35d 100644 --- a/src/pcms/transfer/omega_h_form_integrator_utils.hpp +++ b/src/pcms/transfer/omega_h_form_integrator_utils.hpp @@ -282,8 +282,7 @@ struct IntegrationData // Each specialization provides, for a target element `elm` with local vertex // ids `verts` and vertex coordinates `tgt_verts`: // ndof number of local target DOFs -// Gid(...) global DOF id for local dof k (element id for P0, vertex -// id for P1) — matching OmegaHLagrangeLayout::GetGids() +// Index(...) active PETSc row for local dof k // Values(pt, ...) basis values at the (global) integration point pt template struct TargetTriBasis; @@ -293,11 +292,13 @@ struct TargetTriBasis<0> { static constexpr int ndof = 1; - KOKKOS_INLINE_FUNCTION static GO Gid(const GO* gids, int elm, - const Omega_h::Few&, - int /*k*/) + template + KOKKOS_INLINE_FUNCTION static LO Index(const Permutation& permutation, + int elm, + const Omega_h::Few&, + int /*k*/) { - return gids[elm]; + return permutation(elm); } KOKKOS_INLINE_FUNCTION static void Values( @@ -313,11 +314,12 @@ struct TargetTriBasis<1> { static constexpr int ndof = 3; - KOKKOS_INLINE_FUNCTION static GO Gid( - const GO* gids, int /*elm*/, const Omega_h::Few& verts, - int k) + template + KOKKOS_INLINE_FUNCTION static LO Index( + const Permutation& permutation, int /*elm*/, + const Omega_h::Few& verts, int k) { - return gids[verts[k]]; + return permutation(verts[k]); } // P1 basis functions are the barycentric coordinates of the target element diff --git a/src/pcms/transfer/omega_h_intersection_rhs_integrator.cpp b/src/pcms/transfer/omega_h_intersection_rhs_integrator.cpp index 10103981..84c3cdff 100644 --- a/src/pcms/transfer/omega_h_intersection_rhs_integrator.cpp +++ b/src/pcms/transfer/omega_h_intersection_rhs_integrator.cpp @@ -79,8 +79,7 @@ OmegaHIntersectionRHSIntegrator::BuildDataImpl( auto bary_coords = ip_data.bary_coords; // device view auto weights = ip_data.weights; // device view - const auto device_gids = target_layout.GetGids(); - const GO* gids_ptr = device_gids.data_handle(); + const auto global_to_local = target_layout.GetGlobalToLocalPermutation(); const Omega_h::LOs tgt2src_offsets = intersections.tgt2src_offsets; const Omega_h::LOs tgt2src_indices = intersections.tgt2src_indices; @@ -143,8 +142,8 @@ OmegaHIntersectionRHSIntegrator::BuildDataImpl( coords(global_ip, 0) = pt[0]; coords(global_ip, 1) = pt[1]; for (int k = 0; k < ndof; ++k) { - node_gids(global_ip * ndof + k) = - static_cast(Basis::Gid(gids_ptr, elm, tgt_verts, k)); + node_gids(global_ip * ndof + k) = static_cast( + Basis::Index(global_to_local, elm, tgt_verts, k)); coeffs(global_ip * ndof + k) = basis[k] * w * 2.0 * area; } ++ip_local; @@ -159,7 +158,7 @@ OmegaHIntersectionRHSIntegrator::BuildDataImpl( d.coeffs = std::move(coeffs); d.ndof_per_elem = ndof; d.num_target_dofs = - static_cast(target_layout.GetNumGlobalDofHolder()); + static_cast(target_layout.GetNumOwnedDofHolder()); return d; } diff --git a/src/pcms/transfer/omega_h_mass_integrator.cpp b/src/pcms/transfer/omega_h_mass_integrator.cpp index e5fe259d..d3f4d71d 100644 --- a/src/pcms/transfer/omega_h_mass_integrator.cpp +++ b/src/pcms/transfer/omega_h_mass_integrator.cpp @@ -29,10 +29,9 @@ OmegaHMassIntegrator::OmegaHMassIntegrator( "OmegaHMassIntegrator", "target"); Omega_h::Mesh& mesh = target_layout->GetMesh(); - const auto device_gids = target_layout->GetGids(); - const GO* gids_ptr = device_gids.data_handle(); + const auto global_to_local = target_layout->GetGlobalToLocalPermutation(); const PetscInt num_dofs = - static_cast(target_layout->GetNumGlobalDofHolder()); + static_cast(target_layout->GetNumOwnedDofHolder()); const int nelems = mesh.nelems(); if (target_layout->GetOrder() == 0) { @@ -56,7 +55,7 @@ OmegaHMassIntegrator::OmegaHMassIntegrator( basis[1] = vm[2] - vm[0]; const Omega_h::Real area = Kokkos::fabs(Omega_h::triangle_area_from_basis(basis)); - const PetscInt g = static_cast(gids_ptr[e]); + const PetscInt g = static_cast(global_to_local(e)); coo_rows(e) = g; coo_cols(e) = g; vals(e) = static_cast(area); @@ -95,8 +94,8 @@ OmegaHMassIntegrator::OmegaHMassIntegrator( 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(gids_ptr[verts[i]]); - coo_cols(idx) = static_cast(gids_ptr[verts[j]]); + coo_rows(idx) = static_cast(global_to_local(verts[i])); + coo_cols(idx) = static_cast(global_to_local(verts[j])); } } }); diff --git a/src/pcms/transfer/omega_h_mc_rhs_integrator.cpp b/src/pcms/transfer/omega_h_mc_rhs_integrator.cpp index 212539fc..e009d0a4 100644 --- a/src/pcms/transfer/omega_h_mc_rhs_integrator.cpp +++ b/src/pcms/transfer/omega_h_mc_rhs_integrator.cpp @@ -28,7 +28,8 @@ template KOKKOS_INLINE_FUNCTION void FillElementSamples( const int elm, const int samples_per_element, const Omega_h::Reals& mesh_coords, const Omega_h::LOs& faces2nodes, - const GO* gids, const Kokkos::View& coords, + const Kokkos::View& global_to_local, + const Kokkos::View& coords, const Kokkos::View& node_gids, const Kokkos::View& coeffs, const UnitSampleAt& unit_sample_at) @@ -53,7 +54,7 @@ KOKKOS_INLINE_FUNCTION void FillElementSamples( for (int k = 0; k < 3; ++k) { x += bary[k] * vert_coords[k][0]; y += bary[k] * vert_coords[k][1]; - node_gids(i * 3 + k) = static_cast(gids[verts[k]]); + node_gids(i * 3 + k) = static_cast(global_to_local(verts[k])); coeffs(i * 3 + k) = bary[k] * weight; } coords(i, 0) = x; @@ -85,8 +86,7 @@ OmegaHMonteCarloRHSIntegrator::Data OmegaHMonteCarloRHSIntegrator::BuildData( const int num_samples = nelems * samples_per_element; const auto mesh_coords = mesh.coords(); const auto faces2nodes = mesh.ask_down(Omega_h::FACE, Omega_h::VERT).ab2b; - const auto device_gids = target_layout->GetGids(); - const GO* gids_ptr = device_gids.data_handle(); + const auto global_to_local = target_layout->GetGlobalToLocalPermutation(); Data d; d.coords = @@ -107,8 +107,9 @@ OmegaHMonteCarloRHSIntegrator::Data OmegaHMonteCarloRHSIntegrator::BuildData( "mc_rhs_fill_random", Kokkos::RangePolicy(0, nelems), KOKKOS_LAMBDA(int elm) { auto gen = pool.get_state(); - FillElementSamples(elm, npe, mesh_coords, faces2nodes, gids_ptr, coords, - node_gids, coeffs, [&](int /*s*/, Real& u, Real& v) { + FillElementSamples(elm, npe, mesh_coords, faces2nodes, global_to_local, + coords, node_gids, coeffs, + [&](int /*s*/, Real& u, Real& v) { u = gen.drand(); v = gen.drand(); }); diff --git a/test/test_mesh_intersection_field_transfer.cpp b/test/test_mesh_intersection_field_transfer.cpp index a66579e7..54b16038 100644 --- a/test/test_mesh_intersection_field_transfer.cpp +++ b/test/test_mesh_intersection_field_transfer.cpp @@ -39,10 +39,11 @@ Omega_h::Mesh BuildUnitSquare(Omega_h::Library& lib, const Omega_h::LOs& ev2v) return mesh; } -pcms::LagrangeFunctionSpace MakeP1Space(Omega_h::Mesh& mesh) +pcms::LagrangeFunctionSpace MakeP1Space( + Omega_h::Mesh& mesh, const std::string& global_id_name = "global") { return pcms::LagrangeFunctionSpace::FromMesh( - mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", + mesh, 1, 1, pcms::CoordinateSystem::Cartesian, global_id_name, pcms::LagrangeFunctionSpace::Backend::OmegaH); } @@ -53,6 +54,20 @@ pcms::LagrangeFunctionSpace MakeP0Space(Omega_h::Mesh& mesh) pcms::LagrangeFunctionSpace::Backend::OmegaH); } +void AddReorderedVertexGlobalIds(Omega_h::Mesh& mesh) +{ + mesh.add_tag( + 0, "reordered_global", 1, + Omega_h::Read({2, 0, 3, 1}, "reordered_global")); +} + +void AddSparseVertexGlobalIds(Omega_h::Mesh& mesh) +{ + mesh.add_tag( + 0, "sparse_global", 1, + Omega_h::Read({102, 7, 41, 19}, "sparse_global")); +} + // Integral over the mesh of an order-0 (piecewise-constant) field whose values // are stored in element order. Exact by construction. double IntegrateP0Field(Omega_h::Mesh& mesh, @@ -216,6 +231,84 @@ TEST_CASE("OmegaHConservativeProjection conserves the integral", Catch::Approx(IntegrateP1Field(source_mesh, source)).margin(1e-12)); } +TEST_CASE("OmegaHConservativeProjection writes reordered target GIDs in local " + "DOF order", + "[transfer][mesh_intersection]") +{ + Omega_h::Library lib; + + Omega_h::Mesh source_mesh = + BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3})); + Omega_h::Mesh target_mesh = + BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3})); + AddReorderedVertexGlobalIds(target_mesh); + + auto source_space = MakeP1Space(source_mesh); + auto target_space = MakeP1Space(target_mesh, "reordered_global"); + + auto source = source_space.CreateField(); + auto target = target_space.CreateField(); + + pcms::test::SetField( + source, OMEGA_H_LAMBDA(pcms::Real x, pcms::Real y) { return x + 2.0 * y; }); + + pcms::OmegaHConservativeProjection projection(source_space, target_space); + projection.Apply(source, target); + + const auto target_values = target.GetDOFHolderDataHost(); + const auto target_coords = + target_space.GetLayout()->GetDOFHolderCoordinates().GetValues(); + const auto target_coords_h = pcms::test::CopyCoordinatesToHost( + target_coords, target_mesh.nverts(), target_mesh.dim()); + REQUIRE(static_cast(target_values.extent(0)) == + target_mesh.nverts()); + REQUIRE(target_values.extent(1) == 1); + for (Omega_h::LO i = 0; i < target_mesh.nverts(); ++i) { + const double expected = target_coords_h(i, 0) + 2.0 * target_coords_h(i, 1); + CAPTURE(i); + REQUIRE(target_values(i, 0) == Catch::Approx(expected).margin(1e-9)); + } +} + +TEST_CASE("OmegaHConservativeProjection maps sparse target GIDs to active " + "PETSc rows", + "[transfer][mesh_intersection]") +{ + Omega_h::Library lib; + + Omega_h::Mesh source_mesh = + BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3})); + Omega_h::Mesh target_mesh = + BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3})); + AddSparseVertexGlobalIds(target_mesh); + + auto source_space = MakeP1Space(source_mesh); + auto target_space = MakeP1Space(target_mesh, "sparse_global"); + + auto source = source_space.CreateField(); + auto target = target_space.CreateField(); + + pcms::test::SetField( + source, OMEGA_H_LAMBDA(pcms::Real x, pcms::Real y) { return x + 2.0 * y; }); + + pcms::OmegaHConservativeProjection projection(source_space, target_space); + projection.Apply(source, target); + + const auto target_values = target.GetDOFHolderDataHost(); + const auto target_coords = + target_space.GetLayout()->GetDOFHolderCoordinates().GetValues(); + const auto target_coords_h = pcms::test::CopyCoordinatesToHost( + target_coords, target_mesh.nverts(), target_mesh.dim()); + REQUIRE(static_cast(target_values.extent(0)) == + target_mesh.nverts()); + REQUIRE(target_values.extent(1) == 1); + for (Omega_h::LO i = 0; i < target_mesh.nverts(); ++i) { + const double expected = target_coords_h(i, 0) + 2.0 * target_coords_h(i, 1); + CAPTURE(i); + REQUIRE(target_values(i, 0) == Catch::Approx(expected).margin(1e-9)); + } +} + // P0 (piecewise-constant) source -> P1 target. A constant lives in P1, so it is // reproduced exactly; a non-constant P0 source is not reproduced but its // integral must still be conserved. From af9c175852b8eb50d25bcb44bf664f70f4debe07 Mon Sep 17 00:00:00 2001 From: Jacob Merson Date: Sat, 11 Jul 2026 20:09:40 -0400 Subject: [PATCH 08/15] coupler: separate GID-layout and field messages; skip non-exchanged holders 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. --- src/pcms/coupler/field_exchange_planner.cpp | 105 +++++++++++------- src/pcms/coupler/field_exchange_planner.h | 2 + .../coupler/field_layout_communicator.cpp | 9 +- src/pcms/coupler/field_serializer.h | 9 +- src/pcms/coupler/serializer/xgc.h | 9 +- test/CMakeLists.txt | 1 + test/test_field_exchange_planner.cpp | 64 +++++++++++ 7 files changed, 151 insertions(+), 48 deletions(-) create mode 100644 test/test_field_exchange_planner.cpp diff --git a/src/pcms/coupler/field_exchange_planner.cpp b/src/pcms/coupler/field_exchange_planner.cpp index 4bdcd574..52480d87 100644 --- a/src/pcms/coupler/field_exchange_planner.cpp +++ b/src/pcms/coupler/field_exchange_planner.cpp @@ -1,5 +1,6 @@ #include "pcms/coupler/field_exchange_planner.h" #include "partition.h" +#include "pcms/field/gid_permutation.hpp" #include "pcms/utility/assert.h" #include "pcms/utility/inclusive_scan.h" #include "pcms/utility/profile.h" @@ -45,11 +46,13 @@ static int GetMeshEntityDim(LO local_index, const EntOffsetsArray& ent_offsets) static size_t GetMessageBlockIndex( LO permutation_entry, Rank1View offsets) { - auto begin = offsets.data_handle(); - auto end = begin + offsets.size(); - auto it = std::upper_bound(begin, end, permutation_entry); - PCMS_ALWAYS_ASSERT(it != begin); - return static_cast(std::distance(begin, it - 1)); + for (size_t i = 0; i < offsets.size() - 1; ++i) { + if (permutation_entry >= offsets(i) && permutation_entry < offsets(i + 1)) { + return i; + } + } + PCMS_ALWAYS_ASSERT(false); + return 0; } static OutMsg ConstructOutMessage(const ReversePartitionMap2& reverse_partition) @@ -61,8 +64,7 @@ static OutMsg ConstructOutMessage(const ReversePartitionMap2& reverse_partition) out.dest.reserve(reverse_partition.size()); for (const auto& rank : reverse_partition) { out.dest.push_back(rank.first); - counts.push_back(rank.second.indices.size() + - rank.second.ent_offsets.size()); + counts.push_back(rank.second.indices.size()); } out.offset.resize(counts.size() + 1); out.offset[0] = 0; @@ -73,14 +75,14 @@ static OutMsg ConstructOutMessage(const ReversePartitionMap2& reverse_partition) static redev::LOs ConstructPermutation( const ReversePartitionMap2& reverse_partition, size_t num_entries, - int* length) + int& length) { PCMS_FUNCTION_TIMER; - redev::LOs permutation(num_entries); + // Holders that do not participate in this exchange (non-owned, or owned but + // outside the overlap region) stay unmapped; serialization skips them. + redev::LOs permutation = MakeUnmappedPermutation(num_entries); LO entry = 0; for (const auto& rank : reverse_partition) { - entry += ent_offsets_len; - for (unsigned e = 0; e < rank.second.ent_offsets.size() - 1; ++e) { const int start = rank.second.ent_offsets[e]; const int end = rank.second.ent_offsets[e + 1]; @@ -92,53 +94,58 @@ static redev::LOs ConstructPermutation( } } } - *length = entry; + length = entry; return permutation; } +// Parses a received GID message -- a concatenation of per-sender +// [entity-offset header | GIDs] blocks -- into the permutation mapping each +// local holder to the compact (header-excluded) buffer index of its GID. This +// owns the on-the-wire message format; the GID-to-holder matching is delegated +// to pcms::AppendGidPermutation. Holders absent from the message are left +// unmapped. Returns the total payload length (GID count excluding headers) via +// payload_length. static redev::LOs ConstructPermutation( GlobalIDView local_gids, GlobalIDView received_msg, - const EntOffsetsArray& ent_offsets) + const EntOffsetsArray& ent_offsets, int& payload_length) { PCMS_FUNCTION_TIMER; - std::array, 4> gid_to_buffer_index; + GidToIndexMaps gid_to_buffer_index; size_t offset = 0; + LO payload_offset = 0; while (true) { - GlobalIDView received_offsets( - received_msg.data_handle() + offset, ent_offsets_len); - int length = received_offsets[received_offsets.size() - 1]; - GlobalIDView received_gids( - received_msg.data_handle() + offset + ent_offsets_len, length); + // Guard the header read before dereferencing its last entry (the payload + // length), so a truncated message aborts rather than reading out of bounds. + PCMS_ALWAYS_ASSERT(offset + ent_offsets_len <= received_msg.size()); + int message_length = received_msg(offset + ent_offsets_len - 1); - PCMS_ALWAYS_ASSERT(offset + ent_offsets_len + length - 1 < + PCMS_ALWAYS_ASSERT(offset + ent_offsets_len + message_length - 1 < received_msg.size()); - for (size_t e = 0; e < received_offsets.size() - 1; ++e) { - size_t start = received_offsets[e]; - size_t end = received_offsets[e + 1]; + for (size_t e = 0; e < ent_offsets_len - 1; ++e) { + size_t start = received_msg(offset + e); + size_t end = received_msg(offset + e + 1); for (size_t i = start; i < end; ++i) { - gid_to_buffer_index[e][received_gids[i]] = offset + ent_offsets_len + i; + const GO gid = received_msg(offset + ent_offsets_len + i); + gid_to_buffer_index[e][gid] = payload_offset + i; } } - offset += length + ent_offsets_len; + payload_offset += message_length; + offset += message_length + ent_offsets_len; if (offset >= received_msg.size()) break; } redev::LOs permutation; permutation.reserve(local_gids.size()); - for (size_t e = 0; e < ent_offsets.size() - 1; ++e) { - const auto start = ent_offsets[e]; - const auto end = ent_offsets[e + 1]; - - for (size_t i = start; i < end; ++i) - permutation.push_back(gid_to_buffer_index[e][local_gids[i]]); - } + AppendGidPermutation(local_gids, ent_offsets, gid_to_buffer_index, permutation, + /*allow_missing=*/true); REDEV_ALWAYS_ASSERT(permutation.size() == local_gids.size()); + payload_length = payload_offset; return permutation; } @@ -158,8 +165,11 @@ static OutMsg ConstructOutMessage(int rank, int nproc, totInMsgs - in.srcRanks[(nAppProcs - 1) * nproc + rank]; OutMsg out; for (size_t i = 0; i < nAppProcs; ++i) { - if (senderDeg[i] > 0) + if (senderDeg[i] > 0) { + REDEV_ALWAYS_ASSERT(senderDeg[i] > ent_offsets_len); + senderDeg[i] -= ent_offsets_len; out.dest.push_back(i); + } } redev::GO sum = 0; for (auto deg : senderDeg) { @@ -203,7 +213,7 @@ static ReversePartitionMap2 BuildReversePartitionMap( auto overlap_mask_view = overlap_mask.GetMask(layout); for (LO local_index = 0; local_index < n; ++local_index) { - if (!owned[local_index]) + if (!owned(local_index)) continue; if (!overlap_mask_view[local_index]) @@ -250,7 +260,7 @@ ExchangePlan GenericFieldExchangePlanner::BuildExchangePlan( int length = 0; plan.permutation = - ConstructPermutation(reverse_partition, gids.size(), &length); + ConstructPermutation(reverse_partition, gids.size(), length); plan.msg_size = static_cast(length); return plan; @@ -268,8 +278,10 @@ ExchangePlan GenericFieldExchangePlanner::BuildReceivePlan( auto out_msg = ConstructOutMessage(rank, nproc, in_message_layout); plan.dest_ranks = std::move(out_msg.dest); plan.offsets = std::move(out_msg.offset); - plan.permutation = ConstructPermutation(gids, received_gids, ent_offsets); - plan.msg_size = received_gids.size(); + int length = 0; + plan.permutation = + ConstructPermutation(gids, received_gids, ent_offsets, length); + plan.msg_size = static_cast(length); return plan; } @@ -278,7 +290,9 @@ void GenericFieldExchangePlanner::FillGidMessage( Rank1View gid_message) const { PCMS_FUNCTION_TIMER; - PCMS_ALWAYS_ASSERT(static_cast(gid_message.size()) == plan.msg_size); + const size_t header_size = plan.dest_ranks.size() * ent_offsets_len; + PCMS_ALWAYS_ASSERT(static_cast(gid_message.size()) == + plan.msg_size + header_size); auto gids = layout.GetGidsHost(); auto owned = layout.GetOwnedHost(); @@ -294,9 +308,15 @@ void GenericFieldExchangePlanner::FillGidMessage( continue; LO perm_index = plan.permutation[local_index]; - gid_message[perm_index] = gids[local_index]; - + // Owned holders outside the overlap region carry the sentinel and have no + // slot in the message. + if (perm_index < 0) + continue; auto block_index = GetMessageBlockIndex(perm_index, offsets); + const auto gid_index = + perm_index + static_cast((block_index + 1) * ent_offsets_len); + gid_message(gid_index) = gids(local_index); + int mesh_ent_dim = GetMeshEntityDim(local_index, ent_offsets); for (size_t e = static_cast(mesh_ent_dim) + 1; e < ent_offsets_len; ++e) { @@ -306,9 +326,10 @@ void GenericFieldExchangePlanner::FillGidMessage( for (size_t block_index = 0; block_index < per_rank_offsets.size(); ++block_index) { - auto header_offset = static_cast(plan.offsets[block_index]); + auto header_offset = static_cast(plan.offsets[block_index]) + + block_index * ent_offsets_len; for (size_t e = 0; e < ent_offsets_len; ++e) { - gid_message[header_offset + e] = + gid_message(header_offset + e) = static_cast(per_rank_offsets[block_index][e]); } } diff --git a/src/pcms/coupler/field_exchange_planner.h b/src/pcms/coupler/field_exchange_planner.h index 4ea2fa54..5e5ae29f 100644 --- a/src/pcms/coupler/field_exchange_planner.h +++ b/src/pcms/coupler/field_exchange_planner.h @@ -13,6 +13,8 @@ namespace pcms struct ExchangePlan { redev::LOs dest_ranks; + // Field-payload offsets and permutation. GID-message headers are inserted + // only while exchanging the layout and are not represented in this plan. redev::LOs offsets; std::vector permutation; size_t msg_size = 0; diff --git a/src/pcms/coupler/field_layout_communicator.cpp b/src/pcms/coupler/field_layout_communicator.cpp index 713ba98c..00dabfb7 100644 --- a/src/pcms/coupler/field_layout_communicator.cpp +++ b/src/pcms/coupler/field_layout_communicator.cpp @@ -74,8 +74,13 @@ void FieldLayoutCommunicator::UpdateLayout() // overlap_mask_ is always non-null (created in constructor if not provided) plan_ = planner_->BuildExchangePlan( layout_, redev::Partition{redev_.GetPartition()}, overlap_mask_.get()); - gid_comm_.SetOutMessageLayout(plan_.dest_ranks, plan_.offsets); - std::vector gid_message(plan_.msg_size); + redev::LOs gid_offsets(plan_.offsets.size()); + for (size_t i = 0; i < plan_.offsets.size(); ++i) { + gid_offsets[i] = plan_.offsets[i] + i * ent_offsets_len; + } + gid_comm_.SetOutMessageLayout(plan_.dest_ranks, gid_offsets); + std::vector gid_message(plan_.msg_size + + plan_.dest_ranks.size() * ent_offsets_len); planner_->FillGidMessage( layout_, plan_, Rank1View(gid_message.data(), gid_message.size())); diff --git a/src/pcms/coupler/field_serializer.h b/src/pcms/coupler/field_serializer.h index 87246456..1838fb98 100644 --- a/src/pcms/coupler/field_serializer.h +++ b/src/pcms/coupler/field_serializer.h @@ -28,7 +28,9 @@ class FieldSerializer const LO num_dof = static_cast(data.extent(0)); const LO num_comp = static_cast(data.extent(1)); for (LO i = 0; i < num_dof; ++i) { - 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) { for (LO c = 0; c < num_comp; ++c) { buffer[permutation[i] * num_comp + c] = data(i, c); } @@ -48,7 +50,10 @@ class FieldSerializer Kokkos::View sorted("sorted", layout.OwnedSize()); auto owned = layout.GetOwnedHost(); for (LO i = 0; i < num_dof; ++i) { - if (owned[i]) { + // A negative permutation entry marks a holder outside the exchange (owned + // but outside the overlap region); no data was received for it, so its + // zero-initialized `sorted` slot is left as-is. + if (owned[i] && permutation[i] >= 0) { for (LO c = 0; c < num_comp; ++c) { sorted[i * num_comp + c] = buffer[permutation[i] * num_comp + c]; } diff --git a/src/pcms/coupler/serializer/xgc.h b/src/pcms/coupler/serializer/xgc.h index ce09cf74..ec16d116 100644 --- a/src/pcms/coupler/serializer/xgc.h +++ b/src/pcms/coupler/serializer/xgc.h @@ -41,7 +41,9 @@ class XGCFieldSerializer : public FieldSerializer const LO num_dof = static_cast(data.extent(0)); const LO num_comp = static_cast(data.extent(1)); for (LO i = 0; i < num_dof; ++i) { - 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) { for (LO c = 0; c < num_comp; ++c) { buffer[permutation[i] * num_comp + c] = data(i, c); } @@ -74,7 +76,10 @@ class XGCFieldSerializer : public FieldSerializer } if (rank_participates_) { for (LO i = 0; i < num_dof; ++i) { - if (owned[i]) { + // A negative permutation entry marks a holder outside the exchange + // (owned but outside the overlap region); no data was received for it, + // so it keeps its current value (pre-filled above). + if (owned[i] && permutation[i] >= 0) { for (LO c = 0; c < num_comp; ++c) { full_data[i * num_comp + c] = buffer[permutation[i] * num_comp + c]; } diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 914ab322..5efdc81e 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -385,6 +385,7 @@ if(Catch2_FOUND) test_field_evaluation.cpp test_field_interpolation.cpp test_field_copy.cpp + test_field_exchange_planner.cpp test_point_search.cpp test_mls_basis.cpp test_rbf_interp.cpp diff --git a/test/test_field_exchange_planner.cpp b/test/test_field_exchange_planner.cpp new file mode 100644 index 00000000..d4b318ac --- /dev/null +++ b/test/test_field_exchange_planner.cpp @@ -0,0 +1,64 @@ +#include + +#include +#include + +TEST_CASE("Field exchange plans exclude GID message headers", + "[field_exchange]") +{ + Kokkos::View coords("coords", 4, 2); + pcms::PointCloudLayout layout(2, coords, pcms::CoordinateSystem::Cartesian); + + // Two incoming messages. Each contains a five-entry entity-offset header + // followed by two vertex GIDs. + Kokkos::View received( + "received_gids", 2 * (pcms::ent_offsets_len + 2)); + const pcms::GO message[] = { + 0, 2, 2, 2, 2, 2, 0, 0, 2, 2, 2, 2, 3, 1, + }; + for (size_t i = 0; i < received.size(); ++i) + received(i) = message[i]; + + redev::InMessageLayout in; + in.srcRanks = {0, 7}; + in.offset = {0, 14}; + + pcms::GenericFieldExchangePlanner planner; + const auto plan = planner.BuildReceivePlan( + layout, + pcms::GlobalIDView(received.data(), received.size()), + 0, 1, in); + + REQUIRE(plan.dest_ranks == redev::LOs{0, 1}); + REQUIRE(plan.offsets == redev::LOs{0, 2, 4}); + REQUIRE(plan.permutation == std::vector{1, 3, 0, 2}); + REQUIRE(plan.msg_size == 4); +} + +TEST_CASE("GID messages insert headers around compact field payloads", + "[field_exchange]") +{ + Kokkos::View coords("coords", 4, 2); + pcms::PointCloudLayout layout(2, coords, pcms::CoordinateSystem::Cartesian); + + pcms::ExchangePlan plan; + plan.dest_ranks = {0, 1}; + plan.offsets = {0, 2, 4}; + plan.permutation = {0, 2, 1, 3}; + plan.msg_size = 4; + + Kokkos::View message( + "gid_message", plan.msg_size + 2 * pcms::ent_offsets_len); + pcms::GenericFieldExchangePlanner planner; + planner.FillGidMessage(layout, plan, + pcms::Rank1View( + message.data(), message.size())); + + const pcms::GO expected[] = { + 0, 2, 2, 2, 2, 0, 2, 0, 2, 2, 2, 2, 1, 3, + }; + for (size_t i = 0; i < message.size(); ++i) { + CAPTURE(i); + REQUIRE(message(i) == expected[i]); + } +} From d1f768d3aed1f2c28a77293dd3a627fe7013bbdf Mon Sep 17 00:00:00 2001 From: Jacob Merson Date: Sun, 12 Jul 2026 15:24:43 -0400 Subject: [PATCH 09/15] clang-format field_exchange_planner.cpp --- src/pcms/coupler/field_exchange_planner.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/pcms/coupler/field_exchange_planner.cpp b/src/pcms/coupler/field_exchange_planner.cpp index 52480d87..fb37ce29 100644 --- a/src/pcms/coupler/field_exchange_planner.cpp +++ b/src/pcms/coupler/field_exchange_planner.cpp @@ -141,7 +141,8 @@ static redev::LOs ConstructPermutation( redev::LOs permutation; permutation.reserve(local_gids.size()); - AppendGidPermutation(local_gids, ent_offsets, gid_to_buffer_index, permutation, + AppendGidPermutation(local_gids, ent_offsets, gid_to_buffer_index, + permutation, /*allow_missing=*/true); REDEV_ALWAYS_ASSERT(permutation.size() == local_gids.size()); From ab51122da430d50ba93e90d8e776d34eb9a239ad Mon Sep 17 00:00:00 2001 From: Jacob Merson Date: Sun, 12 Jul 2026 15:24:43 -0400 Subject: [PATCH 10/15] ci: add PETSc lib dir to LD_LIBRARY_PATH for python API tests 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. --- .github/workflows/cmake-test.yml | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/.github/workflows/cmake-test.yml b/.github/workflows/cmake-test.yml index 20c58f9b..d736a952 100644 --- a/.github/workflows/cmake-test.yml +++ b/.github/workflows/cmake-test.yml @@ -182,7 +182,11 @@ jobs: if [ "${{ matrix.meshfields }}" = "ON" ]; then meshfields_lib="${{ runner.temp }}/build-meshFields/install/lib:" fi - echo "LD_LIBRARY_PATH=${{ runner.temp }}/build-kokkos/install/lib:${{ runner.temp }}/build-omega_h/install/lib:${meshfields_lib}${{ runner.temp }}/build-redev/install/lib:${{ runner.temp }}/build-ADIOS2/install/lib:${{ runner.temp }}/build-perfstubs/install/lib:$LD_LIBRARY_PATH" >> $GITHUB_ENV + petsc_lib="" + if [ "${{ matrix.petsc }}" = "ON" ]; then + petsc_lib="${{ runner.temp }}/petsc/ubuntu-kokkos/lib:" + fi + echo "LD_LIBRARY_PATH=${{ runner.temp }}/build-kokkos/install/lib:${{ runner.temp }}/build-omega_h/install/lib:${meshfields_lib}${petsc_lib}${{ runner.temp }}/build-redev/install/lib:${{ runner.temp }}/build-ADIOS2/install/lib:${{ runner.temp }}/build-perfstubs/install/lib:$LD_LIBRARY_PATH" >> $GITHUB_ENV - name: clone petsc if: matrix.petsc == 'ON' From 25ab342006621b148cd11ef891c16b5e2206df8c Mon Sep 17 00:00:00 2001 From: Jacob Merson Date: Sun, 12 Jul 2026 15:50:51 -0400 Subject: [PATCH 11/15] transfer: widen View-size multiplications to size_t 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. --- src/pcms/transfer/omega_h_intersection_rhs_integrator.cpp | 7 ++++--- src/pcms/transfer/omega_h_mc_rhs_integrator.cpp | 8 ++++---- test/test_field_exchange_planner.cpp | 5 +++-- 3 files changed, 11 insertions(+), 9 deletions(-) diff --git a/src/pcms/transfer/omega_h_intersection_rhs_integrator.cpp b/src/pcms/transfer/omega_h_intersection_rhs_integrator.cpp index 84c3cdff..8fc7876e 100644 --- a/src/pcms/transfer/omega_h_intersection_rhs_integrator.cpp +++ b/src/pcms/transfer/omega_h_intersection_rhs_integrator.cpp @@ -108,9 +108,10 @@ OmegaHIntersectionRHSIntegrator::BuildDataImpl( // Pass 2: fill coords, node_gids, and coeffs on device. node_gids/coeffs hold // ndof (target DOFs per element) entries per integration point. Kokkos::View coords("rhs_coords", num_pts, 2); - Kokkos::View node_gids("rhs_node_gids", - num_pts * ndof); - Kokkos::View coeffs("rhs_coeffs", num_pts * ndof); + Kokkos::View node_gids( + "rhs_node_gids", static_cast(num_pts) * ndof); + Kokkos::View coeffs( + "rhs_coeffs", static_cast(num_pts) * ndof); Kokkos::parallel_for( "rhs_fill", nelems, KOKKOS_LAMBDA(int elm) { diff --git a/src/pcms/transfer/omega_h_mc_rhs_integrator.cpp b/src/pcms/transfer/omega_h_mc_rhs_integrator.cpp index e009d0a4..26794ab3 100644 --- a/src/pcms/transfer/omega_h_mc_rhs_integrator.cpp +++ b/src/pcms/transfer/omega_h_mc_rhs_integrator.cpp @@ -91,10 +91,10 @@ OmegaHMonteCarloRHSIntegrator::Data OmegaHMonteCarloRHSIntegrator::BuildData( Data d; d.coords = Kokkos::View("mc_rhs_coords", num_samples, 2); - d.node_gids = Kokkos::View("mc_rhs_node_gids", - num_samples * 3); - d.coeffs = - Kokkos::View("mc_rhs_coeffs", num_samples * 3); + d.node_gids = Kokkos::View( + "mc_rhs_node_gids", static_cast(num_samples) * 3); + d.coeffs = Kokkos::View( + "mc_rhs_coeffs", static_cast(num_samples) * 3); d.nverts = static_cast(mesh.nverts()); const int npe = samples_per_element; diff --git a/test/test_field_exchange_planner.cpp b/test/test_field_exchange_planner.cpp index d4b318ac..44678e47 100644 --- a/test/test_field_exchange_planner.cpp +++ b/test/test_field_exchange_planner.cpp @@ -12,7 +12,7 @@ TEST_CASE("Field exchange plans exclude GID message headers", // Two incoming messages. Each contains a five-entry entity-offset header // followed by two vertex GIDs. Kokkos::View received( - "received_gids", 2 * (pcms::ent_offsets_len + 2)); + "received_gids", 2 * static_cast(pcms::ent_offsets_len + 2)); const pcms::GO message[] = { 0, 2, 2, 2, 2, 2, 0, 0, 2, 2, 2, 2, 3, 1, }; @@ -48,7 +48,8 @@ TEST_CASE("GID messages insert headers around compact field payloads", plan.msg_size = 4; Kokkos::View message( - "gid_message", plan.msg_size + 2 * pcms::ent_offsets_len); + "gid_message", + plan.msg_size + 2 * static_cast(pcms::ent_offsets_len)); pcms::GenericFieldExchangePlanner planner; planner.FillGidMessage(layout, plan, pcms::Rank1View( From 9c1ff5c90267517082ec96e5f990ee36d60c3510 Mon Sep 17 00:00:00 2001 From: Jacob Merson Date: Sun, 12 Jul 2026 18:47:31 -0400 Subject: [PATCH 12/15] transfer: keep PETSc/MeshField backends out of public headers Forward-declare the RHS integrators --- src/pcms/transfer/CMakeLists.txt | 2 +- src/pcms/transfer/monte_carlo_sampling.hpp | 15 +++++++++++++++ .../transfer/omega_h_conservative_projection.cpp | 1 + .../transfer/omega_h_conservative_projection.hpp | 3 +-- .../omega_h_control_variate_projection.cpp | 1 + .../omega_h_control_variate_projection.hpp | 4 ++-- src/pcms/transfer/omega_h_mc_rhs_integrator.hpp | 7 +------ 7 files changed, 22 insertions(+), 11 deletions(-) create mode 100644 src/pcms/transfer/monte_carlo_sampling.hpp diff --git a/src/pcms/transfer/CMakeLists.txt b/src/pcms/transfer/CMakeLists.txt index 2de13931..dddd25a2 100644 --- a/src/pcms/transfer/CMakeLists.txt +++ b/src/pcms/transfer/CMakeLists.txt @@ -65,7 +65,7 @@ target_link_libraries(pcms_transfer PUBLIC Omega_h::omega_h) if(PCMS_ENABLE_MESHFIELDS) - target_link_libraries(pcms_transfer PUBLIC meshfields::meshfields) + target_link_libraries(pcms_transfer PRIVATE meshfields::meshfields) endif() if(PCMS_ENABLE_PETSC AND PCMS_ENABLE_MESHFIELDS) diff --git a/src/pcms/transfer/monte_carlo_sampling.hpp b/src/pcms/transfer/monte_carlo_sampling.hpp new file mode 100644 index 00000000..0e5738e9 --- /dev/null +++ b/src/pcms/transfer/monte_carlo_sampling.hpp @@ -0,0 +1,15 @@ +#ifndef PCMS_TRANSFER_MONTE_CARLO_SAMPLING_HPP +#define PCMS_TRANSFER_MONTE_CARLO_SAMPLING_HPP + +namespace pcms +{ + +// Strategy for generating sample points on each target element. +enum class MonteCarloSampling +{ + UniformRandom, // independent uniform draws per element (XorShift64 pool) +}; + +} // namespace pcms + +#endif // PCMS_TRANSFER_MONTE_CARLO_SAMPLING_HPP diff --git a/src/pcms/transfer/omega_h_conservative_projection.cpp b/src/pcms/transfer/omega_h_conservative_projection.cpp index aac5918c..0faef4d0 100644 --- a/src/pcms/transfer/omega_h_conservative_projection.cpp +++ b/src/pcms/transfer/omega_h_conservative_projection.cpp @@ -1,5 +1,6 @@ #include "pcms/transfer/omega_h_conservative_projection.hpp" #include "pcms/transfer/conservative_projection_solver.hpp" +#include "pcms/transfer/omega_h_intersection_rhs_integrator.hpp" #include "pcms/transfer/omega_h_mass_integrator.hpp" #include "pcms/utility/arrays.h" #include diff --git a/src/pcms/transfer/omega_h_conservative_projection.hpp b/src/pcms/transfer/omega_h_conservative_projection.hpp index bf4fcb93..18b7a2cb 100644 --- a/src/pcms/transfer/omega_h_conservative_projection.hpp +++ b/src/pcms/transfer/omega_h_conservative_projection.hpp @@ -4,7 +4,6 @@ #include "pcms/field/function_space.h" #include "pcms/field/layout/omega_h_lagrange.h" #include "pcms/field/point_evaluator.h" -#include "pcms/transfer/omega_h_intersection_rhs_integrator.hpp" #include "pcms/transfer/transfer_operator.hpp" #include #include @@ -12,8 +11,8 @@ namespace pcms { -// Forward declaration keeps PETSc headers out of this header. class GalerkinProjectionSolver; +class OmegaHIntersectionRHSIntegrator; // Conservative Galerkin projection between Omega_h order-1 Lagrange spaces. // diff --git a/src/pcms/transfer/omega_h_control_variate_projection.cpp b/src/pcms/transfer/omega_h_control_variate_projection.cpp index f30f1bc0..bce1a89a 100644 --- a/src/pcms/transfer/omega_h_control_variate_projection.cpp +++ b/src/pcms/transfer/omega_h_control_variate_projection.cpp @@ -1,6 +1,7 @@ #include "pcms/transfer/omega_h_control_variate_projection.hpp" #include "pcms/transfer/conservative_projection_solver.hpp" #include "pcms/transfer/omega_h_mass_integrator.hpp" +#include "pcms/transfer/omega_h_mc_rhs_integrator.hpp" #include "pcms/utility/arrays.h" #include #include diff --git a/src/pcms/transfer/omega_h_control_variate_projection.hpp b/src/pcms/transfer/omega_h_control_variate_projection.hpp index dd3c93ff..4196ea5d 100644 --- a/src/pcms/transfer/omega_h_control_variate_projection.hpp +++ b/src/pcms/transfer/omega_h_control_variate_projection.hpp @@ -6,7 +6,7 @@ #include "pcms/field/layout/omega_h_lagrange.h" #include "pcms/field/point_evaluator.h" #include "pcms/transfer/interpolator.h" -#include "pcms/transfer/omega_h_mc_rhs_integrator.hpp" +#include "pcms/transfer/monte_carlo_sampling.hpp" #include "pcms/transfer/transfer_operator.hpp" #include #include @@ -14,8 +14,8 @@ namespace pcms { -// Forward declaration keeps PETSc headers out of this header. class GalerkinProjectionSolver; +class OmegaHMonteCarloRHSIntegrator; // Variance-reduced Monte Carlo Galerkin projection using a control variate. // diff --git a/src/pcms/transfer/omega_h_mc_rhs_integrator.hpp b/src/pcms/transfer/omega_h_mc_rhs_integrator.hpp index 1fed3ef4..eb23be95 100644 --- a/src/pcms/transfer/omega_h_mc_rhs_integrator.hpp +++ b/src/pcms/transfer/omega_h_mc_rhs_integrator.hpp @@ -5,6 +5,7 @@ #include "pcms/field/layout/omega_h_lagrange.h" #include "pcms/transfer/integration_point_set.hpp" #include "pcms/transfer/linear_form_integrator.hpp" +#include "pcms/transfer/monte_carlo_sampling.hpp" #include #include #include @@ -12,12 +13,6 @@ namespace pcms { -// Strategy for generating sample points on each target element. -enum class MonteCarloSampling -{ - UniformRandom, // independent uniform draws per element (XorShift64 pool) -}; - // Monte Carlo RHS integrator for conservative L2 projection onto order-1 // Lagrange spaces on Omega_h 2D simplex meshes. // From a44f8c622f827d9881b19e3909035ff427500050 Mon Sep 17 00:00:00 2001 From: Jacob Merson Date: Mon, 13 Jul 2026 02:03:22 -0400 Subject: [PATCH 13/15] Fix GPU build for integrators We had kokkos lambdas in private functions which does not compile. --- src/pcms/transfer/integration_point_set.hpp | 4 +- .../omega_h_intersection_rhs_integrator.cpp | 64 ++++++++--------- .../omega_h_intersection_rhs_integrator.hpp | 27 ------- .../transfer/omega_h_mc_rhs_integrator.cpp | 71 ++++++------------- .../transfer/omega_h_mc_rhs_integrator.hpp | 16 ----- 5 files changed, 57 insertions(+), 125 deletions(-) diff --git a/src/pcms/transfer/integration_point_set.hpp b/src/pcms/transfer/integration_point_set.hpp index 3a91c9e1..5cd3c7fe 100644 --- a/src/pcms/transfer/integration_point_set.hpp +++ b/src/pcms/transfer/integration_point_set.hpp @@ -25,6 +25,8 @@ template class IntegrationPointSet { public: + IntegrationPointSet() = default; + IntegrationPointSet(CoordinateSystem coordinate_system, Kokkos::View coords) : coordinate_system_(coordinate_system), coords_(std::move(coords)) @@ -46,7 +48,7 @@ class IntegrationPointSet } private: - CoordinateSystem coordinate_system_; + CoordinateSystem coordinate_system_ = CoordinateSystem::Cartesian; Kokkos::View coords_; }; diff --git a/src/pcms/transfer/omega_h_intersection_rhs_integrator.cpp b/src/pcms/transfer/omega_h_intersection_rhs_integrator.cpp index 8fc7876e..83300b55 100644 --- a/src/pcms/transfer/omega_h_intersection_rhs_integrator.cpp +++ b/src/pcms/transfer/omega_h_intersection_rhs_integrator.cpp @@ -11,28 +11,29 @@ namespace pcms { -// --------------------------------------------------------------------------- -// OmegaHIntersectionRHSIntegrator::BuildData -// --------------------------------------------------------------------------- -OmegaHIntersectionRHSIntegrator::Data -OmegaHIntersectionRHSIntegrator::BuildData(const FunctionSpace& source_space, - const FunctionSpace& target_space) +namespace { - return BuildData(std::dynamic_pointer_cast( - source_space.GetLayout()), - source_space.GetCoordinateSystem(), - std::dynamic_pointer_cast( - target_space.GetLayout()), - target_space.GetCoordinateSystem()); -} -OmegaHIntersectionRHSIntegrator::Data -OmegaHIntersectionRHSIntegrator::BuildData( - std::shared_ptr source_layout, - CoordinateSystem source_coordinate_system, - std::shared_ptr target_layout, - CoordinateSystem target_coordinate_system) +// Results of the two-pass intersection kernel; the constructor moves these into +// the integrator's members. +struct Data +{ + Kokkos::View coords; + Kokkos::View node_gids; + Kokkos::View coeffs; + int ndof_per_elem = 0; + PetscInt num_target_dofs = 0; +}; + +template +Data BuildDataImpl(const OmegaHLagrangeLayout& source_layout, + const OmegaHLagrangeLayout& target_layout, int quad_order); + +Data BuildData(const std::shared_ptr& source_layout, + CoordinateSystem source_coordinate_system, + const std::shared_ptr& target_layout, + CoordinateSystem target_coordinate_system) { detail::CheckOmegaHScalarLagrangeLayout( source_coordinate_system, source_layout, "OmegaHIntersectionRHSIntegrator", @@ -54,10 +55,8 @@ OmegaHIntersectionRHSIntegrator::BuildData( } template -OmegaHIntersectionRHSIntegrator::Data -OmegaHIntersectionRHSIntegrator::BuildDataImpl( - const OmegaHLagrangeLayout& source_layout, - const OmegaHLagrangeLayout& target_layout, int quad_order) +Data BuildDataImpl(const OmegaHLagrangeLayout& source_layout, + const OmegaHLagrangeLayout& target_layout, int quad_order) { using Basis = detail::TargetTriBasis; constexpr int ndof = Basis::ndof; @@ -163,6 +162,8 @@ OmegaHIntersectionRHSIntegrator::BuildDataImpl( return d; } +} // namespace + // --------------------------------------------------------------------------- // OmegaHIntersectionRHSIntegrator constructors // --------------------------------------------------------------------------- @@ -184,18 +185,15 @@ OmegaHIntersectionRHSIntegrator::OmegaHIntersectionRHSIntegrator( CoordinateSystem source_coordinate_system, std::shared_ptr target_layout, CoordinateSystem target_coordinate_system) - : OmegaHIntersectionRHSIntegrator( - BuildData(std::move(source_layout), source_coordinate_system, - std::move(target_layout), target_coordinate_system)) { -} + Data data = BuildData(source_layout, source_coordinate_system, target_layout, + target_coordinate_system); + point_set_ = IntegrationPointSet( + CoordinateSystem::Cartesian, std::move(data.coords)); + node_gids_ = std::move(data.node_gids); + coeffs_ = std::move(data.coeffs); + ndof_per_elem_ = data.ndof_per_elem; -OmegaHIntersectionRHSIntegrator::OmegaHIntersectionRHSIntegrator(Data data) - : point_set_(CoordinateSystem::Cartesian, std::move(data.coords)), - node_gids_(std::move(data.node_gids)), - coeffs_(std::move(data.coeffs)), - ndof_per_elem_(data.ndof_per_elem) -{ const PetscInt nnz = static_cast(node_gids_.extent(0)); PetscErrorCode ierr = createSeqVec(PETSC_COMM_WORLD, data.num_target_dofs, &vec_); diff --git a/src/pcms/transfer/omega_h_intersection_rhs_integrator.hpp b/src/pcms/transfer/omega_h_intersection_rhs_integrator.hpp index 45e36205..f6a3a9d5 100644 --- a/src/pcms/transfer/omega_h_intersection_rhs_integrator.hpp +++ b/src/pcms/transfer/omega_h_intersection_rhs_integrator.hpp @@ -49,33 +49,6 @@ class OmegaHIntersectionRHSIntegrator : public LinearFormIntegrator Rank2View sampled_values) override; private: - // Internal data bundle produced by the two-pass device kernel. node_gids and - // coeffs are laid out with ndof_per_elem entries per integration point. - struct Data - { - Kokkos::View coords; // [num_pts, 2] - Kokkos::View node_gids; // [num_pts*ndof] - Kokkos::View coeffs; // [num_pts*ndof] - int ndof_per_elem = 0; // target DOFs/element - PetscInt num_target_dofs = 0; // target DOF count - }; - - static Data BuildData(const FunctionSpace& source_space, - const FunctionSpace& target_space); - static Data BuildData( - std::shared_ptr source_layout, - CoordinateSystem source_coordinate_system, - std::shared_ptr target_layout, - CoordinateSystem target_coordinate_system); - - // Order-templated two-pass builder, selected at runtime on the target order. - template - static Data BuildDataImpl(const OmegaHLagrangeLayout& source_layout, - const OmegaHLagrangeLayout& target_layout, - int quad_order); - - OmegaHIntersectionRHSIntegrator(Data data); - Vec vec_ = nullptr; IntegrationPointSet point_set_; Kokkos::View diff --git a/src/pcms/transfer/omega_h_mc_rhs_integrator.cpp b/src/pcms/transfer/omega_h_mc_rhs_integrator.cpp index 26794ab3..c411171b 100644 --- a/src/pcms/transfer/omega_h_mc_rhs_integrator.cpp +++ b/src/pcms/transfer/omega_h_mc_rhs_integrator.cpp @@ -64,14 +64,21 @@ KOKKOS_INLINE_FUNCTION void FillElementSamples( } // namespace -// --------------------------------------------------------------------------- -// OmegaHMonteCarloRHSIntegrator::BuildData -// --------------------------------------------------------------------------- -OmegaHMonteCarloRHSIntegrator::Data OmegaHMonteCarloRHSIntegrator::BuildData( - const std::shared_ptr& target_layout, - CoordinateSystem target_coordinate_system, int samples_per_element, +OmegaHMonteCarloRHSIntegrator::OmegaHMonteCarloRHSIntegrator( + const FunctionSpace& target_space, int samples_per_element, MonteCarloSampling sampling, uint64_t seed) + : OmegaHMonteCarloRHSIntegrator( + std::dynamic_pointer_cast( + target_space.GetLayout()), + target_space.GetCoordinateSystem(), samples_per_element, sampling, seed) +{ +} + +OmegaHMonteCarloRHSIntegrator::OmegaHMonteCarloRHSIntegrator( + std::shared_ptr target_layout, + CoordinateSystem target_coordinate_system, int samples_per_element, + MonteCarloSampling /*sampling*/, uint64_t seed) { detail::CheckOmegaHScalarP1Layout(target_coordinate_system, target_layout, "OmegaHMonteCarloRHSIntegrator", "target"); @@ -88,20 +95,14 @@ OmegaHMonteCarloRHSIntegrator::Data OmegaHMonteCarloRHSIntegrator::BuildData( const auto faces2nodes = mesh.ask_down(Omega_h::FACE, Omega_h::VERT).ab2b; const auto global_to_local = target_layout->GetGlobalToLocalPermutation(); - Data d; - d.coords = - Kokkos::View("mc_rhs_coords", num_samples, 2); - d.node_gids = Kokkos::View( + Kokkos::View coords("mc_rhs_coords", num_samples, + 2); + Kokkos::View node_gids( "mc_rhs_node_gids", static_cast(num_samples) * 3); - d.coeffs = Kokkos::View( + Kokkos::View coeffs( "mc_rhs_coeffs", static_cast(num_samples) * 3); - d.nverts = static_cast(mesh.nverts()); const int npe = samples_per_element; - auto coords = d.coords; - auto node_gids = d.node_gids; - auto coeffs = d.coeffs; - Kokkos::Random_XorShift64_Pool pool(seed); Kokkos::parallel_for( "mc_rhs_fill_random", Kokkos::RangePolicy(0, nelems), @@ -117,40 +118,14 @@ OmegaHMonteCarloRHSIntegrator::Data OmegaHMonteCarloRHSIntegrator::BuildData( }); Kokkos::fence(); - return d; -} - -// --------------------------------------------------------------------------- -// OmegaHMonteCarloRHSIntegrator constructors -// --------------------------------------------------------------------------- + point_set_ = IntegrationPointSet( + CoordinateSystem::Cartesian, std::move(coords)); + node_gids_ = std::move(node_gids); + coeffs_ = std::move(coeffs); -OmegaHMonteCarloRHSIntegrator::OmegaHMonteCarloRHSIntegrator( - const FunctionSpace& target_space, int samples_per_element, - MonteCarloSampling sampling, uint64_t seed) - : OmegaHMonteCarloRHSIntegrator( - std::dynamic_pointer_cast( - target_space.GetLayout()), - target_space.GetCoordinateSystem(), samples_per_element, sampling, seed) -{ -} - -OmegaHMonteCarloRHSIntegrator::OmegaHMonteCarloRHSIntegrator( - std::shared_ptr target_layout, - CoordinateSystem target_coordinate_system, int samples_per_element, - MonteCarloSampling sampling, uint64_t seed) - : OmegaHMonteCarloRHSIntegrator( - BuildData(target_layout, target_coordinate_system, samples_per_element, - sampling, seed)) -{ -} - -OmegaHMonteCarloRHSIntegrator::OmegaHMonteCarloRHSIntegrator(Data data) - : point_set_(CoordinateSystem::Cartesian, std::move(data.coords)), - node_gids_(std::move(data.node_gids)), - coeffs_(std::move(data.coeffs)) -{ const PetscInt nnz = static_cast(node_gids_.extent(0)); - PetscErrorCode ierr = createSeqVec(PETSC_COMM_WORLD, data.nverts, &vec_); + PetscErrorCode ierr = + createSeqVec(PETSC_COMM_WORLD, static_cast(mesh.nverts()), &vec_); CHKERRABORT(PETSC_COMM_WORLD, ierr); ierr = VecSetPreallocationCOO(vec_, nnz, node_gids_.data()); CHKERRABORT(PETSC_COMM_WORLD, ierr); diff --git a/src/pcms/transfer/omega_h_mc_rhs_integrator.hpp b/src/pcms/transfer/omega_h_mc_rhs_integrator.hpp index eb23be95..68fad801 100644 --- a/src/pcms/transfer/omega_h_mc_rhs_integrator.hpp +++ b/src/pcms/transfer/omega_h_mc_rhs_integrator.hpp @@ -52,22 +52,6 @@ class OmegaHMonteCarloRHSIntegrator : public LinearFormIntegrator Rank2View sampled_values) override; private: - // Internal data bundle produced by the sample-generation kernel. - struct Data - { - Kokkos::View coords; // [num_samples, 2] - Kokkos::View node_gids; // [num_samples * 3] - Kokkos::View coeffs; // [num_samples * 3] - PetscInt nverts = 0; // target DOF count - }; - - static Data BuildData( - const std::shared_ptr& target_layout, - CoordinateSystem target_coordinate_system, int samples_per_element, - MonteCarloSampling sampling, uint64_t seed); - - explicit OmegaHMonteCarloRHSIntegrator(Data data); - Vec vec_ = nullptr; IntegrationPointSet point_set_; Kokkos::View node_gids_; // COO indices From b473bb3eba32835fb891280df5f2c330e8bcbfc1 Mon Sep 17 00:00:00 2001 From: Jacob Merson Date: Mon, 13 Jul 2026 02:24:01 -0400 Subject: [PATCH 14/15] extract out more kokkos lambdas to addressable functions --- src/pcms/transfer/omega_h_mass_integrator.cpp | 85 ++++++++++++------- .../transfer/omega_h_mc_rhs_integrator.cpp | 44 ++++++---- 2 files changed, 84 insertions(+), 45 deletions(-) diff --git a/src/pcms/transfer/omega_h_mass_integrator.cpp b/src/pcms/transfer/omega_h_mass_integrator.cpp index d3f4d71d..2105a87f 100644 --- a/src/pcms/transfer/omega_h_mass_integrator.cpp +++ b/src/pcms/transfer/omega_h_mass_integrator.cpp @@ -14,6 +14,58 @@ namespace pcms { +namespace +{ + +// Fills the diagonal COO entries of a P0 (piecewise-constant) mass matrix: one +// entry per element on its own DOF, valued at the element area. +void FillP0MassCoo( + int nelems, const Omega_h::Reals& coords, const Omega_h::LOs& faces2nodes, + const Kokkos::View& global_to_local, + const Kokkos::View& coo_rows, + const Kokkos::View& coo_cols, + const Kokkos::View& vals) +{ + Kokkos::parallel_for( + "mass_p0_diag", nelems, KOKKOS_LAMBDA(int e) { + const auto verts = Omega_h::gather_verts<3>(faces2nodes, e); + const Omega_h::Matrix<2, 3> vm = + Omega_h::gather_vectors<3, 2>(coords, verts); + Omega_h::Few, 2> basis; + basis[0] = vm[1] - vm[0]; + basis[1] = vm[2] - vm[0]; + const Omega_h::Real area = + Kokkos::fabs(Omega_h::triangle_area_from_basis(basis)); + const PetscInt g = static_cast(global_to_local(e)); + coo_rows(e) = g; + coo_cols(e) = g; + vals(e) = static_cast(area); + }); +} + +// Fills the 3x3-block COO sparsity pattern of a P1 (linear) mass matrix: each +// element contributes a dense block coupling its three vertex DOFs. +void FillP1MassCooPattern( + int nelems, const Omega_h::LOs& faces2nodes, + const Kokkos::View& global_to_local, + const Kokkos::View& coo_rows, + const Kokkos::View& coo_cols) +{ + Kokkos::parallel_for( + "mass_coo_pattern", nelems, KOKKOS_LAMBDA(int e) { + const auto verts = Omega_h::gather_verts<3>(faces2nodes, e); + 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(global_to_local(verts[i])); + coo_cols(idx) = static_cast(global_to_local(verts[j])); + } + } + }); +} + +} // namespace + OmegaHMassIntegrator::OmegaHMassIntegrator(const FunctionSpace& target_space) : OmegaHMassIntegrator(std::dynamic_pointer_cast( target_space.GetLayout()), @@ -45,21 +97,8 @@ OmegaHMassIntegrator::OmegaHMassIntegrator( Kokkos::View coo_rows("mass_coo_rows", nnz); Kokkos::View coo_cols("mass_coo_cols", nnz); Kokkos::View vals("mass_vals", nnz); - Kokkos::parallel_for( - "mass_p0_diag", nelems, KOKKOS_LAMBDA(int e) { - const auto verts = Omega_h::gather_verts<3>(faces2nodes, e); - const Omega_h::Matrix<2, 3> vm = - Omega_h::gather_vectors<3, 2>(coords, verts); - Omega_h::Few, 2> basis; - basis[0] = vm[1] - vm[0]; - basis[1] = vm[2] - vm[0]; - const Omega_h::Real area = - Kokkos::fabs(Omega_h::triangle_area_from_basis(basis)); - const PetscInt g = static_cast(global_to_local(e)); - coo_rows(e) = g; - coo_cols(e) = g; - vals(e) = static_cast(area); - }); + FillP0MassCoo(nelems, coords, faces2nodes, global_to_local, coo_rows, + coo_cols, vals); PetscErrorCode ierr = createSeqAIJMat(PETSC_COMM_WORLD, num_dofs, num_dofs, 0, nullptr, &mat_); @@ -88,17 +127,8 @@ OmegaHMassIntegrator::OmegaHMassIntegrator( const PetscInt nnz = static_cast(nelems) * 9; Kokkos::View coo_rows("mass_coo_rows", nnz); Kokkos::View coo_cols("mass_coo_cols", nnz); - Kokkos::parallel_for( - "mass_coo_pattern", nelems, KOKKOS_LAMBDA(int e) { - const auto verts = Omega_h::gather_verts<3>(faces2nodes, e); - 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(global_to_local(verts[i])); - coo_cols(idx) = static_cast(global_to_local(verts[j])); - } - } - }); + FillP1MassCooPattern(nelems, faces2nodes, global_to_local, 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 @@ -124,9 +154,6 @@ Mat OmegaHMassIntegrator::GetMatrix() const noexcept return mat_; } -// --------------------------------------------------------------------------- -// Builder -// --------------------------------------------------------------------------- std::unique_ptr BuildOmegaHMassIntegrator( const FunctionSpace& target_space) diff --git a/src/pcms/transfer/omega_h_mc_rhs_integrator.cpp b/src/pcms/transfer/omega_h_mc_rhs_integrator.cpp index c411171b..9221db70 100644 --- a/src/pcms/transfer/omega_h_mc_rhs_integrator.cpp +++ b/src/pcms/transfer/omega_h_mc_rhs_integrator.cpp @@ -62,8 +62,33 @@ KOKKOS_INLINE_FUNCTION void FillElementSamples( } } -} // namespace +// Samples samples_per_element from a uniform random distribution over the +// target element, writing their coordinates, node GIDs, and coefficients. +void FillElementSamples( + int nelems, int samples_per_element, const Omega_h::Reals& mesh_coords, + const Omega_h::LOs& faces2nodes, + const Kokkos::View& global_to_local, + const Kokkos::View& coords, + const Kokkos::View& node_gids, + const Kokkos::View& coeffs, uint64_t seed) +{ + Kokkos::Random_XorShift64_Pool pool(seed); + Kokkos::parallel_for( + "mc_rhs_fill_random", Kokkos::RangePolicy(0, nelems), + KOKKOS_LAMBDA(int elm) { + auto gen = pool.get_state(); + FillElementSamples(elm, samples_per_element, mesh_coords, faces2nodes, + global_to_local, coords, node_gids, coeffs, + [&](int /*s*/, Real& u, Real& v) { + u = gen.drand(); + v = gen.drand(); + }); + pool.free_state(gen); + }); + Kokkos::fence(); +} +} // namespace OmegaHMonteCarloRHSIntegrator::OmegaHMonteCarloRHSIntegrator( const FunctionSpace& target_space, int samples_per_element, @@ -102,21 +127,8 @@ OmegaHMonteCarloRHSIntegrator::OmegaHMonteCarloRHSIntegrator( Kokkos::View coeffs( "mc_rhs_coeffs", static_cast(num_samples) * 3); - const int npe = samples_per_element; - Kokkos::Random_XorShift64_Pool pool(seed); - Kokkos::parallel_for( - "mc_rhs_fill_random", Kokkos::RangePolicy(0, nelems), - KOKKOS_LAMBDA(int elm) { - auto gen = pool.get_state(); - FillElementSamples(elm, npe, mesh_coords, faces2nodes, global_to_local, - coords, node_gids, coeffs, - [&](int /*s*/, Real& u, Real& v) { - u = gen.drand(); - v = gen.drand(); - }); - pool.free_state(gen); - }); - Kokkos::fence(); + FillElementSamples(nelems, samples_per_element, mesh_coords, faces2nodes, + global_to_local, coords, node_gids, coeffs, seed); point_set_ = IntegrationPointSet( CoordinateSystem::Cartesian, std::move(coords)); From c72ccada304341ff9677b3adcd99f0eb49473990 Mon Sep 17 00:00:00 2001 From: Jacob Merson Date: Mon, 13 Jul 2026 02:46:21 -0400 Subject: [PATCH 15/15] Petsc functions to set the sparsity patterns take CPU pointers 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. --- .../omega_h_intersection_rhs_integrator.cpp | 7 +++++-- src/pcms/transfer/omega_h_mass_integrator.cpp | 15 ++++++++++++--- src/pcms/transfer/omega_h_mc_rhs_integrator.cpp | 6 +++++- 3 files changed, 22 insertions(+), 6 deletions(-) diff --git a/src/pcms/transfer/omega_h_intersection_rhs_integrator.cpp b/src/pcms/transfer/omega_h_intersection_rhs_integrator.cpp index 83300b55..9adb5394 100644 --- a/src/pcms/transfer/omega_h_intersection_rhs_integrator.cpp +++ b/src/pcms/transfer/omega_h_intersection_rhs_integrator.cpp @@ -11,7 +11,6 @@ namespace pcms { - namespace { @@ -198,7 +197,11 @@ OmegaHIntersectionRHSIntegrator::OmegaHIntersectionRHSIntegrator( PetscErrorCode ierr = createSeqVec(PETSC_COMM_WORLD, data.num_target_dofs, &vec_); CHKERRABORT(PETSC_COMM_WORLD, ierr); - ierr = VecSetPreallocationCOO(vec_, nnz, node_gids_.data()); + // VecSetPreallocationCOO takes the COO indices on the host + // TODO: ask Todd/PETSc folks if there is a better way to do this for GPU support + auto node_gids_host = + Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, node_gids_); + ierr = VecSetPreallocationCOO(vec_, nnz, node_gids_host.data()); CHKERRABORT(PETSC_COMM_WORLD, ierr); } diff --git a/src/pcms/transfer/omega_h_mass_integrator.cpp b/src/pcms/transfer/omega_h_mass_integrator.cpp index 2105a87f..29cd90e2 100644 --- a/src/pcms/transfer/omega_h_mass_integrator.cpp +++ b/src/pcms/transfer/omega_h_mass_integrator.cpp @@ -103,7 +103,12 @@ OmegaHMassIntegrator::OmegaHMassIntegrator( PetscErrorCode ierr = createSeqAIJMat(PETSC_COMM_WORLD, num_dofs, num_dofs, 0, nullptr, &mat_); CHKERRABORT(PETSC_COMM_WORLD, ierr); - ierr = MatSetPreallocationCOO(mat_, nnz, coo_rows.data(), coo_cols.data()); + auto coo_rows_host = + Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, coo_rows); + auto coo_cols_host = + Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, coo_cols); + ierr = MatSetPreallocationCOO(mat_, nnz, coo_rows_host.data(), + coo_cols_host.data()); CHKERRABORT(PETSC_COMM_WORLD, ierr); ierr = MatSetValuesCOO(mat_, vals.data(), INSERT_VALUES); CHKERRABORT(PETSC_COMM_WORLD, ierr); @@ -136,7 +141,12 @@ OmegaHMassIntegrator::OmegaHMassIntegrator( PetscErrorCode ierr = createSeqAIJMat(PETSC_COMM_WORLD, num_dofs, num_dofs, 0, nullptr, &mat_); CHKERRABORT(PETSC_COMM_WORLD, ierr); - ierr = MatSetPreallocationCOO(mat_, nnz, coo_rows.data(), coo_cols.data()); + auto coo_rows_host = + Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, coo_rows); + auto coo_cols_host = + Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, coo_cols); + ierr = MatSetPreallocationCOO(mat_, nnz, coo_rows_host.data(), + coo_cols_host.data()); CHKERRABORT(PETSC_COMM_WORLD, ierr); ierr = MatSetValuesCOO(mat_, elm_mass_dev.data(), INSERT_VALUES); CHKERRABORT(PETSC_COMM_WORLD, ierr); @@ -154,7 +164,6 @@ Mat OmegaHMassIntegrator::GetMatrix() const noexcept return mat_; } - std::unique_ptr BuildOmegaHMassIntegrator( const FunctionSpace& target_space) { diff --git a/src/pcms/transfer/omega_h_mc_rhs_integrator.cpp b/src/pcms/transfer/omega_h_mc_rhs_integrator.cpp index 9221db70..ecd57088 100644 --- a/src/pcms/transfer/omega_h_mc_rhs_integrator.cpp +++ b/src/pcms/transfer/omega_h_mc_rhs_integrator.cpp @@ -139,7 +139,11 @@ OmegaHMonteCarloRHSIntegrator::OmegaHMonteCarloRHSIntegrator( PetscErrorCode ierr = createSeqVec(PETSC_COMM_WORLD, static_cast(mesh.nverts()), &vec_); CHKERRABORT(PETSC_COMM_WORLD, ierr); - ierr = VecSetPreallocationCOO(vec_, nnz, node_gids_.data()); + // VecSetPreallocationCOO takes the COO indices on the host + // TODO: ask Todd/PETSc folks if there is a better way to do this for GPU support + auto node_gids_host = + Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, node_gids_); + ierr = VecSetPreallocationCOO(vec_, nnz, node_gids_host.data()); CHKERRABORT(PETSC_COMM_WORLD, ierr); }