diff --git a/CMakeLists.txt b/CMakeLists.txt index 9d407a53..40be9ed6 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -49,6 +49,15 @@ option(PCMS_ENABLE_Python "Enable pcms Python api" OFF) option(PCMS_ENABLE_PRINT "PCMS print statements enabled" ON) +# Per-thread on-stack ring-buffer capacities for the adjacency-based +# intersection/adjacency BFS (queue_visited.hpp). Defaults are sized for 3D +# tetrahedral meshes on the serial/OpenMP backends. Lower them for device +# (GPU) builds, where each search thread's stack must hold these arrays. +set(PCMS_INTERSECTION_QUEUE_SIZE 1024 CACHE STRING + "Max BFS queue size per target element in the intersection search") +set(PCMS_INTERSECTION_TRACK_SIZE 2048 CACHE STRING + "Max visited-set size per target element in the intersection search") + option(PETSC_LINK_STATIC "Use pkg-config --static results for PETSc" ${_pcms_link_petsc_static_default}) diff --git a/src/pcms/configuration.h.in b/src/pcms/configuration.h.in index 51859ec0..24585146 100644 --- a/src/pcms/configuration.h.in +++ b/src/pcms/configuration.h.in @@ -7,3 +7,6 @@ #cmakedefine PCMS_ENABLE_Fortran #cmakedefine PCMS_ENABLE_MESHFIELDS #cmakedefine PCMS_ENABLE_PETSC + +#cmakedefine PCMS_INTERSECTION_QUEUE_SIZE @PCMS_INTERSECTION_QUEUE_SIZE@ +#cmakedefine PCMS_INTERSECTION_TRACK_SIZE @PCMS_INTERSECTION_TRACK_SIZE@ diff --git a/src/pcms/field/layout/omega_h_lagrange.cpp b/src/pcms/field/layout/omega_h_lagrange.cpp index a325140f..c9d7bdb9 100644 --- a/src/pcms/field/layout/omega_h_lagrange.cpp +++ b/src/pcms/field/layout/omega_h_lagrange.cpp @@ -171,6 +171,7 @@ OmegaHLagrangeLayout::OmegaHLagrangeLayout( owned_ = BuildOwned(mesh_, entity_dim, owned_mask); owned_host_ = Kokkos::View("owned_host", owned_.size()); + Kokkos::deep_copy(owned_host_, owned_); class_ids_ = Omega_h::Read( mesh_.get_array(entity_dim, "class_id")); diff --git a/src/pcms/localization/queue_visited.hpp b/src/pcms/localization/queue_visited.hpp index 8d768eec..3db919e8 100644 --- a/src/pcms/localization/queue_visited.hpp +++ b/src/pcms/localization/queue_visited.hpp @@ -7,8 +7,23 @@ #include #include #include -#define MAX_SIZE_QUEUE 500 -#define MAX_SIZE_TRACK 800 +#include "pcms/configuration.h" +// Per-thread BFS ring-buffer capacities for the adjacency-based intersection +// search. Sized for 3D: a target tetrahedron can overlap many more source +// elements than a 2D triangle, and the visited set additionally holds the +// non-intersecting neighbors probed along the way. These are on-stack arrays, +// so raising them raises per-thread stack usage; lower them for device (GPU) +// builds via the PCMS_INTERSECTION_{QUEUE,TRACK}_SIZE CMake cache variables, +// which flow in through pcms/configuration.h. The fallbacks below apply only if +// the configured header is unavailable (e.g. a header-only consumer). +#ifndef PCMS_INTERSECTION_QUEUE_SIZE +#define PCMS_INTERSECTION_QUEUE_SIZE 1024 +#endif +#ifndef PCMS_INTERSECTION_TRACK_SIZE +#define PCMS_INTERSECTION_TRACK_SIZE 2048 +#endif +#define MAX_SIZE_QUEUE PCMS_INTERSECTION_QUEUE_SIZE +#define MAX_SIZE_TRACK PCMS_INTERSECTION_TRACK_SIZE namespace pcms { diff --git a/src/pcms/transfer/mass_matrix_integrator.hpp b/src/pcms/transfer/mass_matrix_integrator.hpp index 7d36c975..8eabf475 100644 --- a/src/pcms/transfer/mass_matrix_integrator.hpp +++ b/src/pcms/transfer/mass_matrix_integrator.hpp @@ -17,16 +17,19 @@ template class MassMatrixIntegrator : public MeshField::Integrator { public: + // Linear simplex: numNodes = spatial dim + 1 (3 for triangles, 4 for tets). + static constexpr int numNodes = FieldElement::MeshEntDim + 1; + MassMatrixIntegrator(Omega_h::Mesh& mesh_in, FieldElement& fe_in, int order = 2) : mesh(mesh_in), fe(fe_in), - subMatrixSize(3 * 3), // FIXME remove hard coded size - elmMassMatrix("elmMassMatrix", mesh_in.nelems() * 3 * 3), + subMatrixSize(numNodes * numNodes), + elmMassMatrix("elmMassMatrix", mesh_in.nelems() * numNodes * numNodes), Integrator(order) { Kokkos::deep_copy(elmMassMatrix, 0); - assert(mesh.dim() == 2); // TODO support 1d,2d,3d + assert(mesh.dim() == 2 || mesh.dim() == 3); assert(mesh.family() == OMEGA_H_SIMPLEX); } void atPoints(Kokkos::View p, @@ -58,12 +61,18 @@ class MassMatrixIntegrator : public MeshField::Integrator } const auto N = shapeFn.getValues(localCoord); const auto wPt = w(pt); - const auto dVPt = dV(pt); + // Use the unsigned volume element: MeshField returns a signed + // Jacobian determinant, which is negative for tetrahedra whose vertex + // ordering has negative orientation. A mass matrix integrates against + // the positive volume measure, so take the magnitude (a no-op in 2D + // where the differential area is already positive). + const auto dVPt = Kokkos::fabs(dV(pt)); // printf("Shape Functions: %f, %f, %f \n", N[0], N[1], N[2]); // printf("wPt, dVPt: %f, %f \n", wPt, dVPt); for (auto i = 0; i < N.size(); i++) { for (auto j = 0; j < N.size(); j++) { - massMatrix(elm * subMat + i * 3 + j) += N[i] * N[j] * wPt * dVPt; + massMatrix(elm * subMat + i * numNodes + j) += + N[i] * N[j] * wPt * dVPt; } } } diff --git a/src/pcms/transfer/mesh_intersection.cpp b/src/pcms/transfer/mesh_intersection.cpp index 564ae350..06f26d20 100644 --- a/src/pcms/transfer/mesh_intersection.cpp +++ b/src/pcms/transfer/mesh_intersection.cpp @@ -4,51 +4,66 @@ namespace pcms { +namespace +{ +// Construct the source-mesh containing-element search appropriate to the +// spatial dimension: a uniform 20^Dim background grid over the source mesh. +template +auto MakeGridPointSearch(Omega_h::Mesh& source_mesh) +{ + if constexpr (Dim == 3) { + return pcms::GridPointSearch3D(source_mesh, 20, 20, 20); + } else { + return pcms::GridPointSearch2D(source_mesh, 20, 20); + } +} +} // namespace + +template void FindIntersections::adjBasedIntersectSearch( const Omega_h::LOs& tgt2src_offsets, Omega_h::Write& nIntersections, Omega_h::Write& tgt2src_indices, bool is_count_only) { - + // Element entity dimension equals the spatial dimension (FACE for 2D, REGION + // for 3D); measures are triangle areas (2D) or tet volumes (3D). 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; - const auto& src_elem_areas = measure_elements_real(&source_mesh_); - const auto& tgt_elem_areas = measure_elements_real(&target_mesh_); + const auto& tgt_elems2nodes = target_mesh_.ask_down(Dim, Omega_h::VERT).ab2b; + const auto& src_elems2nodes = source_mesh_.ask_down(Dim, Omega_h::VERT).ab2b; + const auto& src_elem_measures = measure_elements_real(&source_mesh_); + const auto& tgt_elem_measures = measure_elements_real(&target_mesh_); const auto& t2t = source_mesh_.ask_dual(); // gives connected element neighbors const auto& t2tt = t2t.a2ab; const auto& tt2t = t2t.ab2b; - const auto flat_centroids = - pcms::get_entity_centroids(target_mesh_, Omega_h::FACE); + const auto flat_centroids = pcms::get_entity_centroids(target_mesh_, Dim); // Convert layout_right 1D Omega_h array to 2D Kokkos view with correct layout - auto centroids = ConvertCoordsTo2D(flat_centroids, target_mesh_.nfaces(), 2); + auto centroids = + ConvertCoordsTo2D(flat_centroids, target_mesh_.nelems(), Dim); - pcms::GridPointSearch2D search_cell(source_mesh_, 20, 20); + auto search_cell = MakeGridPointSearch(source_mesh_); auto results = search_cell(centroids); auto owning_cell_ids = search_cell.GetOwningElementIds(results); - auto nfaces_target = target_mesh_.nfaces(); + auto nelems_target = target_mesh_.nelems(); Omega_h::parallel_for( - nfaces_target, + nelems_target, OMEGA_H_LAMBDA(const Omega_h::LO id) { Queue queue; Track visited; auto current_cell_id = owning_cell_ids(id); - auto current_tgt_elm_area = tgt_elem_areas[id]; + auto current_tgt_elm_measure = tgt_elem_measures[id]; OMEGA_H_CHECK_PRINTF(current_cell_id >= 0, "ERROR: source cell id not found for given target " - "centroid %d (%f, %f)\n", - id, centroids(id, 0), centroids(id, 1)); + "centroid %d\n", + id); auto tgt_elm_vert_coords = - get_vert_coords_of_elem(tgt_coords, tgt_faces2nodes, id); + get_vert_coords_of_elem(tgt_coords, tgt_elems2nodes, id); Omega_h::LO start_counter; if (!is_count_only) { @@ -76,22 +91,42 @@ void FindIntersections::adjBasedIntersectSearch( auto neighborElmId = tt2t[i]; if (visited.notVisited(neighborElmId)) { - visited.push_back(neighborElmId); - auto elm_vert_coords = get_vert_coords_of_elem( - src_coords, src_faces2nodes, neighborElmId); - r3d::Polytope<2> intersection; + // If the visited buffer is full, skip this neighbor so the BFS + // terminates. Without this, an unrecorded neighbor stays "not + // visited" and is re-queued forever (infinite loop). Mirrors the + // guard in adj_search.cpp. Raise MAX_SIZE_TRACK if this fires. + if (!visited.push_back(neighborElmId)) { + printf("ERROR: visited buffer full (MAX_SIZE_TRACK=%d) for " + "target %d; some intersections may be missed\n", + MAX_SIZE_TRACK, id); + continue; + } + auto elm_vert_coords = get_vert_coords_of_elem( + src_coords, src_elems2nodes, neighborElmId); + r3d::Polytope intersection; r3d::intersect_simplices(intersection, tgt_elm_vert_coords, elm_vert_coords); - auto intersected_area = r3d::measure(intersection); - auto current_src_elm_area = src_elem_areas[neighborElmId]; + // Take the magnitude: r3d::measure is signed by the orientation of + // the target simplex used to initialize the polytope, which can be + // negative for tetrahedra. This mirrors the fabs applied in the + // sub-simplex decomposition and mass assembly; without it a + // negatively-oriented target element would reject all of its real + // overlaps and break conservation. + auto intersected_measure = Kokkos::fabs(r3d::measure(intersection)); + auto current_src_elm_measure = src_elem_measures[neighborElmId]; auto scale = - Kokkos::fmax(current_tgt_elm_area, current_src_elm_area); + Kokkos::fmax(current_tgt_elm_measure, current_src_elm_measure); auto eps = Kokkos::fmax(abs_tol, rel_tol * scale); - if (intersection.nverts >= 3 && intersected_area >= eps) { + // A valid intersection is a non-degenerate simplex-simplex overlap: + // at least Dim+1 vertices (a polygon in 2D, a polyhedron in 3D). + if (intersection.nverts >= Dim + 1 && intersected_measure >= eps) { count++; OMEGA_H_CHECK_PRINTF( - count < 500, "WARNING: count exceeds 500 for target %d", id); + count < kMaxIntersectionsPerTarget, + "intersection count for target %d reached the cap %d; raise " + "kMaxIntersectionsPerTarget/MAX_SIZE_QUEUE", + id, kMaxIntersectionsPerTarget); queue.push_back(neighborElmId); @@ -113,20 +148,32 @@ void FindIntersections::adjBasedIntersectSearch( }, // end of lambda "count the number of intersections for each target element"); } -IntersectionResults intersectTargets(Omega_h::Mesh& source_mesh, - Omega_h::Mesh& target_mesh) + +// Explicit instantiations for the supported spatial dimensions. +template void FindIntersections::adjBasedIntersectSearch<2>( + const Omega_h::LOs&, Omega_h::Write&, + Omega_h::Write&, bool); +template void FindIntersections::adjBasedIntersectSearch<3>( + const Omega_h::LOs&, Omega_h::Write&, + Omega_h::Write&, bool); + +namespace +{ +template +IntersectionResults intersectTargetsImpl(Omega_h::Mesh& source_mesh, + Omega_h::Mesh& target_mesh) { FindIntersections intersect(source_mesh, target_mesh); - auto nfaces_target = target_mesh.nfaces(); + auto nelems_target = target_mesh.nelems(); Omega_h::Write nIntersections( - nfaces_target, 0, "number of intersections in each target vertex"); + nelems_target, 0, "number of intersections in each target element"); Omega_h::Write tgt2src_indices; - intersect.adjBasedIntersectSearch(Omega_h::LOs(), nIntersections, - tgt2src_indices, true); + intersect.adjBasedIntersectSearch(Omega_h::LOs(), nIntersections, + tgt2src_indices, true); Kokkos::fence(); auto tgt2src_offsets = Omega_h::offset_scan(Omega_h::Read(nIntersections), @@ -139,9 +186,20 @@ IntersectionResults intersectTargets(Omega_h::Mesh& source_mesh, ntotal_intersections, 0, "indices of the source elements that intersect the given target element"); - intersect.adjBasedIntersectSearch(tgt2src_offsets, nIntersections, - tgt2src_indices, false); + intersect.adjBasedIntersectSearch(tgt2src_offsets, nIntersections, + tgt2src_indices, false); return {.tgt2src_offsets = tgt2src_offsets, .tgt2src_indices = Omega_h::read(tgt2src_indices)}; } +} // namespace + +IntersectionResults intersectTargets(Omega_h::Mesh& source_mesh, + Omega_h::Mesh& target_mesh) +{ + OMEGA_H_CHECK(source_mesh.dim() == target_mesh.dim()); + if (source_mesh.dim() == 3) { + return intersectTargetsImpl<3>(source_mesh, target_mesh); + } + return intersectTargetsImpl<2>(source_mesh, target_mesh); +} } // namespace pcms diff --git a/src/pcms/transfer/mesh_intersection.hpp b/src/pcms/transfer/mesh_intersection.hpp index 84a52114..f72b5c86 100644 --- a/src/pcms/transfer/mesh_intersection.hpp +++ b/src/pcms/transfer/mesh_intersection.hpp @@ -15,19 +15,29 @@ namespace pcms constexpr static double abs_tol = 1e-18; /// abs tolerance constexpr static double rel_tol = 1e-12; /// rel tolerance -[[nodiscard]] OMEGA_H_INLINE r3d::Few, 3> +// Upper bound on intersecting source elements recorded per target element. +// Must not exceed the BFS queue capacity (MAX_SIZE_QUEUE); raised for 3D where +// a single target tet can overlap many source tets. +constexpr static int kMaxIntersectionsPerTarget = MAX_SIZE_QUEUE; + +// Gather the vertex coordinates of a simplex element (triangle for Dim==2, +// tetrahedron for Dim==3) into an r3d simplex, ready for +// r3d::intersect_simplices. +template +[[nodiscard]] OMEGA_H_INLINE r3d::Few, Dim + 1> get_vert_coords_of_elem(const Omega_h::Reals& coords, - const Omega_h::LOs& faces2nodes, const int id) + const Omega_h::LOs& elems2nodes, const int id) { - const auto elm_verts = Omega_h::gather_verts<3>(faces2nodes, id); + const auto elm_verts = Omega_h::gather_verts(elems2nodes, id); - const Omega_h::Matrix<2, 3> elm_vert_coords = - Omega_h::gather_vectors<3, 2>(coords, elm_verts); + const Omega_h::Matrix elm_vert_coords = + Omega_h::gather_vectors(coords, elm_verts); - r3d::Few, 3> r3d_vector; - for (int i = 0; i < 3; ++i) { - r3d_vector[i][0] = elm_vert_coords[i][0]; - r3d_vector[i][1] = elm_vert_coords[i][1]; + r3d::Few, Dim + 1> r3d_vector; + for (int i = 0; i < Dim + 1; ++i) { + for (int d = 0; d < Dim; ++d) { + r3d_vector[i][d] = elm_vert_coords[i][d]; + } } return r3d_vector; @@ -78,11 +88,13 @@ class FindIntersections * @param is_count_only If true, only counts intersections; if false, also * fills tgt2src_indices. * - * @note This method assumes 2D linear triangles and uses - * `r3d::intersect_simplices` for geometric intersection. + * @note Templated on spatial dimension `Dim`: linear triangles (Dim==2) or + * linear tetrahedra (Dim==3), using `r3d::intersect_simplices` for geometric + * intersection. * * @see r3d::intersect_simplices, intersectTargets */ + template void adjBasedIntersectSearch(const Omega_h::LOs& tgt2src_offsets, Omega_h::Write& nIntersections, Omega_h::Write& tgt2src_indices, diff --git a/src/pcms/transfer/omega_h_form_integrator_utils.hpp b/src/pcms/transfer/omega_h_form_integrator_utils.hpp index 23f5e35d..ee2b5194 100644 --- a/src/pcms/transfer/omega_h_form_integrator_utils.hpp +++ b/src/pcms/transfer/omega_h_form_integrator_utils.hpp @@ -13,9 +13,10 @@ 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( +// Shared checks for a scalar Cartesian Lagrange space on a simplex mesh +// (triangles in 2D, tetrahedra in 3D), independent of order. Order is validated +// separately by the callers below. +inline void CheckOmegaHScalarSimplexLayout( CoordinateSystem coordinate_system, const std::shared_ptr& layout, const char* context, const char* role) @@ -33,12 +34,13 @@ inline void CheckOmegaHScalarSimplex2DLayout( " 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.dim() != 2 && mesh.dim() != 3) { + throw pcms_error(std::string(context) + ": " + role + + " mesh must be 2D or 3D"); } if (mesh.family() != OMEGA_H_SIMPLEX) { throw pcms_error(std::string(context) + ": " + role + - " mesh must be a simplex (triangle) mesh"); + " mesh must be a simplex (triangle/tetrahedron) mesh"); } } @@ -48,7 +50,7 @@ inline void CheckOmegaHScalarP1Layout( const std::shared_ptr& layout, const char* context, const char* role) { - CheckOmegaHScalarSimplex2DLayout(coordinate_system, layout, context, role); + CheckOmegaHScalarSimplexLayout(coordinate_system, layout, context, role); if (layout->GetOrder() != 1) { throw pcms_error(std::string(context) + ": " + role + " space must be order-1"); @@ -63,7 +65,7 @@ inline void CheckOmegaHScalarLagrangeLayout( const std::shared_ptr& layout, const char* context, const char* role) { - CheckOmegaHScalarSimplex2DLayout(coordinate_system, layout, context, role); + CheckOmegaHScalarSimplexLayout(coordinate_system, layout, context, role); const int order = layout->GetOrder(); if (order != 0 && order != 1) { throw pcms_error(std::string(context) + ": " + role + @@ -71,28 +73,23 @@ inline void CheckOmegaHScalarLagrangeLayout( } } -[[nodiscard]] OMEGA_H_INLINE Omega_h::Vector<2> GlobalFromBarycentric( - const MeshField::Vector3& barycentric_coord, - const Omega_h::Few, 3>& verts_coord) +// Map barycentric coordinates on a simplex (Dim+1 barycentric components) to +// the global Cartesian point, given the simplex's Dim+1 vertex coordinates. +template +[[nodiscard]] OMEGA_H_INLINE Omega_h::Vector GlobalFromBarycentric( + const Omega_h::Vector& barycentric_coord, + const Omega_h::Few, Dim + 1>& 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]; + Omega_h::Vector real_coords; + for (int d = 0; d < Dim; ++d) { + real_coords[d] = 0.0; } - 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]; + for (int i = 0; i < Dim + 1; ++i) { + for (int d = 0; d < Dim; ++d) { + real_coords[d] += barycentric_coord[i] * verts_coord[i][d]; + } } - return Omega_h::barycentric_from_global<2, 2>(point, omegah_vector); + return real_coords; } [[nodiscard]] OMEGA_H_INLINE int RemoveDuplicateVerticesAndFixLinks( @@ -188,22 +185,24 @@ inline void CheckOmegaHScalarLagrangeLayout( return new_n; } -template -OMEGA_H_INLINE void ForEachIntersectionSubtriangle( +// 2D: fan the clipped intersection polygon into triangles anchored at vertex 0, +// invoking op(sub_triangle, src_elm, area) for each non-degenerate piece. +template +OMEGA_H_INLINE void ForEachIntersectionSubtriangleImpl( 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) + const Omega_h::LOs& tgt_elems2nodes, const Omega_h::LOs& src_elems2nodes, + SimplexOp&& op) { auto tgt_elm_vert_coords = - get_vert_coords_of_elem(tgt_coords, tgt_faces2nodes, elm); + get_vert_coords_of_elem<2>(tgt_coords, tgt_elems2nodes, 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); + get_vert_coords_of_elem<2>(src_coords, src_elems2nodes, 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); @@ -231,40 +230,201 @@ OMEGA_H_INLINE void ForEachIntersectionSubtriangle( continue; } - op(tri_coords, tgt_elm_vert_coords, src_elm_vert_coords, current_src_elm, - area); + op(tri_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. +// Walk each face of a clipped r3d polyhedron exactly once and fan it into +// triangles, invoking op(v0, v1, v2) for each triangle (v0 is the face's anchor +// vertex, so a face with k vertices yields k-2 triangles). Every vertex of an +// r3d clipped Polytope<3> has exactly three face-neighbors (pnbrs); marking +// each directed edge as it is consumed guarantees every face is emitted once. +// This mirrors the edge-marking traversal buried inside r3d::reduce, which r3d +// does not expose for reuse, so the traversal is reproduced here once and +// shared. Vertices arrive as r3d::Vector<3> (indexable [0..2]). +template +OMEGA_H_INLINE void ForEachPolytopeFaceTriangle(const r3d::Polytope<3>& poly, + TriangleOp&& op) +{ + // emarks[v][p] == 1 once the directed edge (v, pnbr p) has been consumed. + int emarks[r3d::Polytope<3>::max_verts][3] = {{}}; + for (int vstart = 0; vstart < poly.nverts; ++vstart) { + for (int pstart = 0; pstart < 3; ++pstart) { + if (emarks[vstart][pstart]) { + continue; + } + int pnext = pstart; + int vcur = vstart; + emarks[vcur][pnext] = 1; + int vnext = poly.verts[vcur].pnbrs[pnext]; + const auto face_v0 = poly.verts[vcur].pos; + + // Move to the second edge of this face. + int np = 0; + for (np = 0; np < 3; ++np) { + if (poly.verts[vnext].pnbrs[np] == vcur) { + break; + } + } + vcur = vnext; + pnext = (np + 1) % 3; + emarks[vcur][pnext] = 1; + vnext = poly.verts[vcur].pnbrs[pnext]; + + // Fan the face into triangles anchored at face_v0. + while (vnext != vstart) { + op(face_v0, poly.verts[vnext].pos, poly.verts[vcur].pos); + + // Advance around the face. + for (np = 0; np < 3; ++np) { + if (poly.verts[vnext].pnbrs[np] == vcur) { + break; + } + } + vcur = vnext; + pnext = (np + 1) % 3; + emarks[vcur][pnext] = 1; + vnext = poly.verts[vcur].pnbrs[pnext]; + } + } + } +} + +// 3D: star-decompose the clipped intersection polyhedron into tetrahedra from +// its centroid. The centroid lies strictly inside the convex intersection, so +// lifting each boundary-face triangle (enumerated by +// ForEachPolytopeFaceTriangle) to the centroid tiles the polyhedron without +// overlap; op(sub_tet, src_elm, volume) fires for each non-degenerate piece. +template +OMEGA_H_INLINE void ForEachIntersectionSubtetImpl( + const int elm, const IntersectionResults& intersection, + const Omega_h::Reals& tgt_coords, const Omega_h::Reals& src_coords, + const Omega_h::LOs& tgt_elems2nodes, const Omega_h::LOs& src_elems2nodes, + SimplexOp&& op) +{ + auto tgt_elm_vert_coords = + get_vert_coords_of_elem<3>(tgt_coords, tgt_elems2nodes, 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<3>(src_coords, src_elems2nodes, current_src_elm); + r3d::Polytope<3> poly; + r3d::intersect_simplices(poly, tgt_elm_vert_coords, src_elm_vert_coords); + if (poly.nverts < 4) { + continue; + } + const double poly_vol = Kokkos::fabs(r3d::measure(poly)); + const double eps_vol = abs_tol + rel_tol * poly_vol; + + // Centroid of the (convex) intersection polyhedron: interior apex. + Omega_h::Vector<3> apex = {0.0, 0.0, 0.0}; + for (int v = 0; v < poly.nverts; ++v) { + apex[0] += poly.verts[v].pos[0]; + apex[1] += poly.verts[v].pos[1]; + apex[2] += poly.verts[v].pos[2]; + } + apex[0] /= poly.nverts; + apex[1] /= poly.nverts; + apex[2] /= poly.nverts; + + // Lift each boundary-face triangle to the interior centroid to form a tet. + ForEachPolytopeFaceTriangle(poly, [&](const r3d::Vector<3>& a, + const r3d::Vector<3>& b, + const r3d::Vector<3>& c) { + Omega_h::Few, 4> tet_coords; + tet_coords[0] = apex; + tet_coords[1] = {a[0], a[1], a[2]}; + tet_coords[2] = {b[0], b[1], b[2]}; + tet_coords[3] = {c[0], c[1], c[2]}; + + Omega_h::Few, 3> basis; + basis[0] = tet_coords[1] - tet_coords[0]; + basis[1] = tet_coords[2] - tet_coords[0]; + basis[2] = tet_coords[3] - tet_coords[0]; + + const Omega_h::Real vol = + Kokkos::fabs(Omega_h::tet_volume_from_basis(basis)); + if (vol > eps_vol) { + op(tet_coords, current_src_elm, vol); + } + }); + } +} + +// Dimension-generic driver over the sub-simplices (triangles in 2D, tets in 3D) +// that tile each target element's intersection with the source mesh. Invokes +// op(sub_simplex_coords, src_elm, measure) for every non-degenerate piece. +template +OMEGA_H_INLINE void ForEachIntersectionSubsimplex( + const int elm, const IntersectionResults& intersection, + const Omega_h::Reals& tgt_coords, const Omega_h::Reals& src_coords, + const Omega_h::LOs& tgt_elems2nodes, const Omega_h::LOs& src_elems2nodes, + SimplexOp&& op) +{ + if constexpr (Dim == 3) { + ForEachIntersectionSubtetImpl(elm, intersection, tgt_coords, src_coords, + tgt_elems2nodes, src_elems2nodes, op); + } else { + ForEachIntersectionSubtriangleImpl(elm, intersection, tgt_coords, + src_coords, tgt_elems2nodes, + src_elems2nodes, op); + } +} + +// Maps spatial dimension to the MeshField simplex topology whose reference +// quadrature rules we use (triangle in 2D, tetrahedron in 3D). +template +struct SimplexTopology; +template <> +struct SimplexTopology<2> +{ + static constexpr MeshField::Mesh_Topology value = MeshField::Triangle; +}; +template <> +struct SimplexTopology<3> +{ + static constexpr MeshField::Mesh_Topology value = MeshField::Tetrahedron; +}; + +// Barycentric integration points and weights for a reference simplex (triangle +// in 2D, tetrahedron in 3D), taken from MeshField's predefined 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. +// Kokkos views once at construction. Each barycentric point has Dim+1 +// components. // // 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. +template struct IntegrationData { - Kokkos::View bary_coords; // barycentric coordinates - Kokkos::View weights; // quadrature weights + Kokkos::View + bary_coords; // barycentric coordinates + Kokkos::View weights; // quadrature weights explicit IntegrationData(int order) { - auto ip_vec = MeshField::getIntegrationPoints(order); + auto ip_vec = + MeshField::getIntegrationPoints::value>(order); const std::size_t num_ip = ip_vec.size(); - bary_coords = Kokkos::View("bary_coords", num_ip); + 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; + for (int d = 0; d < Dim + 1; ++d) { + bary_coords_host(i, d) = ip_vec[i].param[d]; + } weights_host(i) = ip_vec[i].weight; } Kokkos::deep_copy(bary_coords, bary_coords_host); @@ -274,9 +434,10 @@ struct IntegrationData int size() const { return bary_coords.extent(0); } }; -// Target Lagrange basis on a triangle, parameterized by element order, for the +// Target Lagrange basis on a simplex (triangle in 2D, tetrahedron in 3D), +// parameterized by spatial dimension and 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 +// DOF; order 1 is the Dim+1 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 @@ -284,40 +445,39 @@ struct IntegrationData // ndof number of local target DOFs // Index(...) active PETSc row for local dof k // Values(pt, ...) basis values at the (global) integration point pt -template -struct TargetTriBasis; +template +struct TargetSimplexBasis; -template <> -struct TargetTriBasis<0> +template +struct TargetSimplexBasis { static constexpr int ndof = 1; template - KOKKOS_INLINE_FUNCTION static LO Index(const Permutation& permutation, - int elm, - const Omega_h::Few&, - int /*k*/) + KOKKOS_INLINE_FUNCTION static LO Index( + const Permutation& permutation, int elm, + const Omega_h::Few&, int /*k*/) { return permutation(elm); } KOKKOS_INLINE_FUNCTION static void Values( - const Omega_h::Vector<2>&, const Omega_h::Few, 3>&, - Omega_h::Real out[ndof]) + const Omega_h::Vector&, + const Omega_h::Few, Dim + 1>&, Omega_h::Real out[ndof]) { out[0] = 1.0; } }; -template <> -struct TargetTriBasis<1> +template +struct TargetSimplexBasis { - static constexpr int ndof = 3; + static constexpr int ndof = Dim + 1; template KOKKOS_INLINE_FUNCTION static LO Index( const Permutation& permutation, int /*elm*/, - const Omega_h::Few& verts, int k) + const Omega_h::Few& verts, int k) { return permutation(verts[k]); } @@ -325,14 +485,14 @@ struct TargetTriBasis<1> // 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, + const Omega_h::Vector& pt, + const Omega_h::Few, Dim + 1>& 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]; + const auto bary = Omega_h::barycentric_from_global(pt, tgt_verts); + for (int i = 0; i < ndof; ++i) { + out[i] = bary[i]; + } } }; diff --git a/src/pcms/transfer/omega_h_intersection_rhs_integrator.cpp b/src/pcms/transfer/omega_h_intersection_rhs_integrator.cpp index ae9b293f..f3c16df1 100644 --- a/src/pcms/transfer/omega_h_intersection_rhs_integrator.cpp +++ b/src/pcms/transfer/omega_h_intersection_rhs_integrator.cpp @@ -25,7 +25,7 @@ struct Data PetscInt num_target_dofs = 0; }; -template +template Data BuildDataImpl(const OmegaHLagrangeLayout& source_layout, const OmegaHLagrangeLayout& target_layout, int quad_order); @@ -41,24 +41,39 @@ Data BuildData(const std::shared_ptr& source_layout, target_coordinate_system, target_layout, "OmegaHIntersectionRHSIntegrator", "target"); + const int dim = target_layout->GetMesh().dim(); + if (source_layout->GetMesh().dim() != dim) { + throw pcms_error("OmegaHIntersectionRHSIntegrator: source and target mesh " + "dimensions differ"); + } + // The integrand f_src * phi_target has polynomial degree source_order + - // target_order on each intersection subtriangle; integrate it exactly (with a + // target_order on each intersection sub-simplex; 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); + if (dim == 3) { + return BuildDataImpl<3, TgtOrder>(*source_layout, *target_layout, + quad_order); + } + return BuildDataImpl<2, TgtOrder>(*source_layout, *target_layout, + quad_order); }); } -template +template Data BuildDataImpl(const OmegaHLagrangeLayout& source_layout, const OmegaHLagrangeLayout& target_layout, int quad_order) { - using Basis = detail::TargetTriBasis; + using Basis = detail::TargetSimplexBasis; constexpr int ndof = Basis::ndof; + // Reference-to-physical Jacobian factor for a simplex: the reference simplex + // measure is 1/Dim! (1/2 in 2D, 1/6 in 3D), so a physical sub-simplex of + // measure `m` scales the reference quadrature weights by Dim! * m. + constexpr Omega_h::Real ref_factor = (Dim == 3) ? 6.0 : 2.0; Omega_h::Mesh& source_mesh = source_layout.GetMesh(); Omega_h::Mesh& target_mesh = target_layout.GetMesh(); @@ -66,13 +81,11 @@ Data BuildDataImpl(const OmegaHLagrangeLayout& source_layout, 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& tgt_elems2nodes = target_mesh.ask_down(Dim, 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; + const auto& src_elems2nodes = source_mesh.ask_down(Dim, Omega_h::VERT).ab2b; - detail::IntegrationData ip_data(quad_order); + 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 @@ -88,12 +101,10 @@ Data BuildDataImpl(const OmegaHLagrangeLayout& source_layout, Kokkos::parallel_for( "rhs_count", nelems, KOKKOS_LAMBDA(int elm) { int count = 0; - detail::ForEachIntersectionSubtriangle( + detail::ForEachIntersectionSubsimplex( 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, + tgt_elems2nodes, src_elems2nodes, + [&](const Omega_h::Few, Dim + 1>&, int, Omega_h::Real) { count += npts; }); ip_counts[elm] = count; }); @@ -105,7 +116,7 @@ Data BuildDataImpl(const OmegaHLagrangeLayout& source_layout, // 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 coords("rhs_coords", num_pts, Dim); Kokkos::View node_gids( "rhs_node_gids", static_cast(num_pts) * ndof); Kokkos::View coeffs( @@ -113,37 +124,45 @@ Data BuildDataImpl(const OmegaHLagrangeLayout& source_layout, 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]}; + const auto tgt_verts = + Omega_h::gather_verts(tgt_elems2nodes, elm); + const Omega_h::Matrix tgt_vert_mat = + Omega_h::gather_vectors(tgt_coords, tgt_verts); + Omega_h::Few, Dim + 1> tgt_omh; + for (int i = 0; i < Dim + 1; ++i) { + for (int d = 0; d < Dim; ++d) { + tgt_omh[i][d] = tgt_vert_mat[i][d]; + } + } int ip_local = 0; const int offset = ip_offsets[elm]; - detail::ForEachIntersectionSubtriangle( + detail::ForEachIntersectionSubsimplex( 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) { + tgt_elems2nodes, src_elems2nodes, + [&](const Omega_h::Few, Dim + 1>& sub, int, + Omega_h::Real measure) { for (int ip_idx = 0; ip_idx < npts; ++ip_idx) { - const auto bary = bary_coords(ip_idx); + Omega_h::Vector bary; + for (int d = 0; d < Dim + 1; ++d) { + bary[d] = bary_coords(ip_idx, d); + } const double w = weights(ip_idx); - const auto pt = detail::GlobalFromBarycentric(bary, tri); + const auto pt = detail::GlobalFromBarycentric(bary, sub); 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 d = 0; d < Dim; ++d) { + coords(global_ip, d) = pt[d]; + } for (int k = 0; k < ndof; ++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; + coeffs(global_ip * ndof + k) = + basis[k] * w * ref_factor * measure; } ++ip_local; } diff --git a/src/pcms/transfer/omega_h_mass_integrator.cpp b/src/pcms/transfer/omega_h_mass_integrator.cpp index 29cd90e2..99edb7eb 100644 --- a/src/pcms/transfer/omega_h_mass_integrator.cpp +++ b/src/pcms/transfer/omega_h_mass_integrator.cpp @@ -17,10 +17,24 @@ namespace pcms namespace { +// Measure (area in 2D, volume in 3D) of a simplex from its vertex-difference +// basis. +template +KOKKOS_INLINE_FUNCTION Omega_h::Real SimplexMeasure( + const Omega_h::Few, Dim>& basis) +{ + if constexpr (Dim == 3) { + return Kokkos::fabs(Omega_h::tet_volume_from_basis(basis)); + } else { + return Kokkos::fabs(Omega_h::triangle_area_from_basis(basis)); + } +} + // 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. +// entry per element on its own DOF, valued at the element measure. +template void FillP0MassCoo( - int nelems, const Omega_h::Reals& coords, const Omega_h::LOs& faces2nodes, + int nelems, const Omega_h::Reals& coords, const Omega_h::LOs& elems2nodes, const Kokkos::View& global_to_local, const Kokkos::View& coo_rows, const Kokkos::View& coo_cols, @@ -28,35 +42,37 @@ void FillP0MassCoo( { 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 auto verts = Omega_h::gather_verts(elems2nodes, e); + const Omega_h::Matrix vm = + Omega_h::gather_vectors(coords, verts); + Omega_h::Few, Dim> basis; + for (int d = 0; d < Dim; ++d) { + basis[d] = vm[d + 1] - vm[0]; + } + const Omega_h::Real measure = SimplexMeasure(basis); const PetscInt g = static_cast(global_to_local(e)); coo_rows(e) = g; coo_cols(e) = g; - vals(e) = static_cast(area); + vals(e) = static_cast(measure); }); } -// Fills the 3x3-block COO sparsity pattern of a P1 (linear) mass matrix: each -// element contributes a dense block coupling its three vertex DOFs. +// Fills the (Dim+1)x(Dim+1)-block COO sparsity pattern of a P1 (linear) mass +// matrix: each element contributes a dense block coupling its vertex DOFs. +template void FillP1MassCooPattern( - int nelems, const Omega_h::LOs& faces2nodes, + int nelems, const Omega_h::LOs& elems2nodes, const Kokkos::View& global_to_local, const Kokkos::View& coo_rows, const Kokkos::View& coo_cols) { + constexpr int nv = Dim + 1; 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; + const auto verts = Omega_h::gather_verts(elems2nodes, e); + for (int i = 0; i < nv; ++i) { + for (int j = 0; j < nv; ++j) { + const int idx = e * (nv * nv) + i * nv + j; coo_rows(idx) = static_cast(global_to_local(verts[i])); coo_cols(idx) = static_cast(global_to_local(verts[j])); } @@ -64,92 +80,115 @@ void FillP1MassCooPattern( }); } -} // namespace - -OmegaHMassIntegrator::OmegaHMassIntegrator(const FunctionSpace& target_space) - : OmegaHMassIntegrator(std::dynamic_pointer_cast( - target_space.GetLayout()), - target_space.GetCoordinateSystem()) +// Assembles the target-space mass matrix for spatial dimension Dim (triangles +// for Dim==2, tetrahedra for Dim==3) and returns the owned PETSc matrix. +template +Mat BuildOmegaHMassMatrixImpl(Omega_h::Mesh& mesh, + const OmegaHLagrangeLayout& target_layout) { -} - -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 global_to_local = target_layout->GetGlobalToLocalPermutation(); + const auto global_to_local = target_layout.GetGlobalToLocalPermutation(); const PetscInt num_dofs = - static_cast(target_layout->GetNumOwnedDofHolder()); + static_cast(target_layout.GetNumOwnedDofHolder()); const int nelems = mesh.nelems(); + Mat mat = nullptr; - if (target_layout->GetOrder() == 0) { + 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 + // the mass matrix is diagonal with M_ee = measure(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 auto& elems2nodes = mesh.ask_down(Dim, 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); - FillP0MassCoo(nelems, coords, faces2nodes, global_to_local, coo_rows, - coo_cols, vals); + FillP0MassCoo(nelems, coords, elems2nodes, global_to_local, coo_rows, + coo_cols, vals); PetscErrorCode ierr = - createSeqAIJMat(PETSC_COMM_WORLD, num_dofs, num_dofs, 0, nullptr, &mat_); + createSeqAIJMat(PETSC_COMM_WORLD, num_dofs, num_dofs, 0, nullptr, &mat); CHKERRABORT(PETSC_COMM_WORLD, ierr); 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(), + ierr = MatSetPreallocationCOO(mat, nnz, coo_rows_host.data(), coo_cols_host.data()); CHKERRABORT(PETSC_COMM_WORLD, ierr); - ierr = MatSetValuesCOO(mat_, vals.data(), INSERT_VALUES); + ierr = MatSetValuesCOO(mat, vals.data(), INSERT_VALUES); CHKERRABORT(PETSC_COMM_WORLD, ierr); - return; + return mat; } - // 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); + Kokkos::View elm_mass_dev; + if constexpr (Dim == 3) { + const auto [shp, map] = MeshField::Omegah::getTetrahedronElement<1>(mesh); + MeshField::FieldElement coordFe(mesh.nelems(), coordField, shp, map); + elm_mass_dev = buildElementMassMatrix(mesh, coordFe); + } else { + const auto [shp, map] = MeshField::Omegah::getTriangleElement<1>(mesh); + MeshField::FieldElement coordFe(mesh.nelems(), coordField, shp, map); + 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; + // Build COO sparsity pattern on device: each element contributes a + // (Dim+1)x(Dim+1) block. + const auto& elems2nodes = mesh.ask_down(Dim, Omega_h::VERT).ab2b; - const PetscInt nnz = static_cast(nelems) * 9; + constexpr int nv = Dim + 1; + const PetscInt nnz = static_cast(nelems) * (nv * nv); Kokkos::View coo_rows("mass_coo_rows", nnz); Kokkos::View coo_cols("mass_coo_cols", nnz); - FillP1MassCooPattern(nelems, faces2nodes, global_to_local, coo_rows, - coo_cols); + FillP1MassCooPattern(nelems, elems2nodes, 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 // it can be passed directly to MatSetValuesCOO — no host copy needed. PetscErrorCode ierr = - createSeqAIJMat(PETSC_COMM_WORLD, num_dofs, num_dofs, 0, nullptr, &mat_); + createSeqAIJMat(PETSC_COMM_WORLD, num_dofs, num_dofs, 0, nullptr, &mat); CHKERRABORT(PETSC_COMM_WORLD, ierr); 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(), + 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); + ierr = MatSetValuesCOO(mat, elm_mass_dev.data(), INSERT_VALUES); CHKERRABORT(PETSC_COMM_WORLD, ierr); + return mat; +} + +} // namespace + +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(); + if (mesh.dim() == 3) { + mat_ = BuildOmegaHMassMatrixImpl<3>(mesh, *target_layout); + } else { + mat_ = BuildOmegaHMassMatrixImpl<2>(mesh, *target_layout); + } } OmegaHMassIntegrator::~OmegaHMassIntegrator() diff --git a/src/pcms/transfer/omega_h_mc_rhs_integrator.cpp b/src/pcms/transfer/omega_h_mc_rhs_integrator.cpp index c57ce9f1..05016379 100644 --- a/src/pcms/transfer/omega_h_mc_rhs_integrator.cpp +++ b/src/pcms/transfer/omega_h_mc_rhs_integrator.cpp @@ -12,9 +12,9 @@ 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 +// Maps a uniform point on the unit square to uniform barycentric coordinates on +// a triangle. 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) { @@ -22,51 +22,96 @@ KOKKOS_INLINE_FUNCTION Omega_h::Vector<3> UniformTriangleBarycentric(Real u, return {1.0 - s, s * (1.0 - v), s * v}; } +// Maps three uniform draws on the unit cube to uniform barycentric coordinates +// on a tetrahedron via the cut-and-fold method (Rocchini & Cignoni, +// "Generating Random Points in a Tetrahedron", J. Graphics Tools 2000). +KOKKOS_INLINE_FUNCTION Omega_h::Vector<4> UniformTetBarycentric(Real s, Real t, + Real u) +{ + if (s + t > 1.0) { // fold the cube into a prism + s = 1.0 - s; + t = 1.0 - t; + } + if (t + u > 1.0) { // fold the prism into a tetrahedron + const Real tmp = u; + u = 1.0 - s - t; + t = 1.0 - tmp; + } else if (s + t + u > 1.0) { + const Real tmp = u; + u = s + t + u - 1.0; + s = 1.0 - t - tmp; + } + const Real a = 1.0 - s - t - u; + return {a, s, t, u}; +} + +template +KOKKOS_INLINE_FUNCTION Omega_h::Vector UniformSimplexBarycentric( + const Real r[Dim]) +{ + if constexpr (Dim == 3) { + return UniformTetBarycentric(r[0], r[1], r[2]); + } else { + return UniformTriangleBarycentric(r[0], r[1]); + } +} + // 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 +// unit_sample_at(s, r) fills the s-th sample's Dim unit-hypercube draws. +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 Omega_h::Reals& mesh_coords, const Omega_h::LOs& elems2nodes, const Kokkos::View& global_to_local, 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; + constexpr int nv = Dim + 1; + const auto verts = Omega_h::gather_verts(elems2nodes, elm); + const auto vert_coords = Omega_h::gather_vectors(mesh_coords, verts); + Omega_h::Few, Dim> basis; + for (int d = 0; d < Dim; ++d) { + basis[d] = vert_coords[d + 1] - vert_coords[0]; + } + Real measure; + if constexpr (Dim == 3) { + measure = Kokkos::fabs(Omega_h::tet_volume_from_basis(basis)); + } else { + measure = Kokkos::fabs(Omega_h::triangle_area_from_basis(basis)); + } + const Real weight = measure / 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); + Real r[Dim]; + unit_sample_at(s, r); + const auto bary = UniformSimplexBarycentric(r); 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(global_to_local(verts[k])); - coeffs(i * 3 + k) = bary[k] * weight; + Omega_h::Vector x; + for (int d = 0; d < Dim; ++d) { + x[d] = 0.0; + } + for (int k = 0; k < nv; ++k) { + for (int d = 0; d < Dim; ++d) { + x[d] += bary[k] * vert_coords[k][d]; + } + node_gids(i * nv + k) = static_cast(global_to_local(verts[k])); + coeffs(i * nv + k) = bary[k] * weight; + } + for (int d = 0; d < Dim; ++d) { + coords(i, d) = x[d]; } - coords(i, 0) = x; - coords(i, 1) = y; } } -// Samples samples_per_element from a uniform random distribution over the +// Samples samples_per_element from a uniform random distribution over each // target element, writing their coordinates, node GIDs, and coefficients. +template void FillElementSamples( int nelems, int samples_per_element, const Omega_h::Reals& mesh_coords, - const Omega_h::LOs& faces2nodes, + const Omega_h::LOs& elems2nodes, const Kokkos::View& global_to_local, const Kokkos::View& coords, const Kokkos::View& node_gids, @@ -77,12 +122,13 @@ void FillElementSamples( "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(); - }); + FillElementSamples(elm, samples_per_element, mesh_coords, + elems2nodes, global_to_local, coords, node_gids, + coeffs, [&](int /*s*/, Real r[Dim]) { + for (int d = 0; d < Dim; ++d) { + r[d] = gen.drand(); + } + }); pool.free_state(gen); }); Kokkos::fence(); @@ -113,22 +159,29 @@ OmegaHMonteCarloRHSIntegrator::OmegaHMonteCarloRHSIntegrator( } Omega_h::Mesh& mesh = target_layout->GetMesh(); + const int dim = mesh.dim(); + nbary_ = dim + 1; 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 elems2nodes = mesh.ask_down(dim, Omega_h::VERT).ab2b; const auto global_to_local = target_layout->GetGlobalToLocalPermutation(); Kokkos::View coords("mc_rhs_coords", num_samples, - 2); + dim); Kokkos::View node_gids( - "mc_rhs_node_gids", static_cast(num_samples) * 3); + "mc_rhs_node_gids", static_cast(num_samples) * nbary_); Kokkos::View coeffs( - "mc_rhs_coeffs", static_cast(num_samples) * 3); + "mc_rhs_coeffs", static_cast(num_samples) * nbary_); - FillElementSamples(nelems, samples_per_element, mesh_coords, faces2nodes, - global_to_local, coords, node_gids, coeffs, seed); + if (dim == 3) { + FillElementSamples<3>(nelems, samples_per_element, mesh_coords, elems2nodes, + global_to_local, coords, node_gids, coeffs, seed); + } else { + FillElementSamples<2>(nelems, samples_per_element, mesh_coords, elems2nodes, + global_to_local, coords, node_gids, coeffs, seed); + } coords_ = std::move(coords); node_gids_ = std::move(node_gids); @@ -173,8 +226,9 @@ Vec OmegaHMonteCarloRHSIntegrator::GetVector() const noexcept void OmegaHMonteCarloRHSIntegrator::Assemble( Rank2View sampled_values) { + const int nbary = nbary_; const std::size_t num_samples = - static_cast(node_gids_.extent(0) / 3); + static_cast(node_gids_.extent(0) / nbary); PCMS_ALWAYS_ASSERT(static_cast(sampled_values.extent(0)) == num_samples); PCMS_ALWAYS_ASSERT(sampled_values.extent(1) >= 1); @@ -187,13 +241,14 @@ void OmegaHMonteCarloRHSIntegrator::Assemble( sampled_values.data_handle(), sampled_values.extent(0), sampled_values.extent(1)); Kokkos::View coo_vals("mc_rhs_coo_vals", - num_samples * 3); + num_samples * nbary); 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; + for (int j = 0; j < nbary; ++j) { + coo_vals(i * nbary + j) = + static_cast(coeffs(i * nbary + j)) * f; } }); diff --git a/src/pcms/transfer/omega_h_mc_rhs_integrator.hpp b/src/pcms/transfer/omega_h_mc_rhs_integrator.hpp index c9b533b4..5d8946eb 100644 --- a/src/pcms/transfer/omega_h_mc_rhs_integrator.hpp +++ b/src/pcms/transfer/omega_h_mc_rhs_integrator.hpp @@ -14,7 +14,8 @@ namespace pcms { // Monte Carlo RHS integrator for conservative L2 projection onto order-1 -// Lagrange spaces on Omega_h 2D simplex meshes. +// Lagrange spaces on Omega_h simplex meshes (triangles in 2D, tetrahedra in +// 3D). // // Instead of intersection-based quadrature, each load-vector entry // b_j = \int phi_j f dx @@ -57,6 +58,7 @@ class OmegaHMonteCarloRHSIntegrator : public LinearFormIntegrator coords_; // [num_pts][dim] sample coordinates Kokkos::View node_gids_; // COO indices Kokkos::View coeffs_; // basis * |T| / N + int nbary_ = 3; // barycentric DOFs per sample (Dim+1): 3 in 2D, 4 in 3D }; } // namespace pcms diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 5efdc81e..06a24e7a 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -413,7 +413,8 @@ if(Catch2_FOUND) 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) + test_mesh_intersection_field_transfer.cpp + test_omega_h_3d_conservative_projection.cpp) endif() add_executable(unit_tests ${PCMS_UNIT_TEST_SOURCES}) diff --git a/test/field_test_utils.h b/test/field_test_utils.h index 5aaac78e..2ddd2f45 100644 --- a/test/field_test_utils.h +++ b/test/field_test_utils.h @@ -3,6 +3,7 @@ #include #include +#include "pcms/configuration.h" #include "pcms/field/field.h" #include "pcms/field/field_data.h" #include "pcms/field/field_evaluator_factory.h" @@ -11,10 +12,22 @@ #include "pcms/field/out_of_bounds_policy.h" #include "pcms/coupler/field_serializer.h" #include "pcms/field/coordinate_system.h" -#include "pcms/localization/adj_search.hpp" #include "pcms/utility/arrays.h" #include "pcms/utility/memory_spaces.h" +#ifdef PCMS_ENABLE_OMEGA_H +#include +#include +#include +#include "pcms/field/function_space/lagrange.h" +#include "pcms/localization/adj_search.hpp" +#endif +#if defined(PCMS_ENABLE_PETSC) && defined(PCMS_ENABLE_MESHFIELDS) +#include "pcms/transfer/linear_form_integrator.hpp" +#endif #include +#include +#include +#include #include #include @@ -36,6 +49,62 @@ inline std::vector StandardEvalCoords2D() return {0.1, 0.2, 0.5, 0.5, 0.7, 0.3, 0.9, 0.1, 0.2, 0.8}; } +// Points strictly outside a unit [0,1]^2 box mesh. +inline std::vector StandardOutsideCoords2D() +{ + return {-0.5, 0.5, 1.5, 0.5, 0.5, -0.5, 0.5, 1.5}; +} + +#ifdef PCMS_ENABLE_OMEGA_H +// 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. +inline 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; +} + +// Convenience overload: diagonal=0 splits the square along vertices (1,3), +// diagonal=1 along (0,2). +inline Omega_h::Mesh BuildUnitSquare(Omega_h::Library& lib, int diagonal) +{ + return BuildUnitSquare(lib, (diagonal == 0) + ? Omega_h::LOs({0, 1, 3, 1, 2, 3}) + : Omega_h::LOs({0, 1, 2, 0, 2, 3})); +} + +inline std::shared_ptr MakeP1Space( + Omega_h::Mesh& mesh, const std::string& global_id_name = "global") +{ + return LagrangeFunctionSpace::FromMesh( + mesh, 1, 1, CoordinateSystem::Cartesian, global_id_name, + LagrangeFunctionSpace::Backend::OmegaH); +} + +inline std::shared_ptr MakeP0Space(Omega_h::Mesh& mesh) +{ + return LagrangeFunctionSpace::FromMesh( + mesh, 0, 1, CoordinateSystem::Cartesian, "global", + LagrangeFunctionSpace::Backend::OmegaH); +} + inline bool AreArraysEqualUnordered( const Omega_h::HostRead& array1, const Omega_h::HostRead& array2, int start, int end) @@ -76,6 +145,40 @@ inline std::vector CopyOmegaHRealsToVector(const Omega_h::Reals& coords) coords_read.data() + coords_read.size()); } +inline double IntegrateP0Field(Omega_h::Mesh& mesh, const Field& field) +{ + const auto values = FlattenToRank1View(field.GetDOFHolderDataHost()); + const auto measures = Omega_h::measure_elements_real(&mesh); + const auto measures_h = Omega_h::HostRead(measures); + + double integral = 0.0; + for (Omega_h::LO e = 0; e < mesh.nelems(); ++e) { + integral += measures_h[e] * values[e]; + } + return integral; +} + +inline double IntegrateP1Field(Omega_h::Mesh& mesh, const Field& field) +{ + const auto values = FlattenToRank1View(field.GetDOFHolderDataHost()); + const auto measures = Omega_h::measure_elements_real(&mesh); + const auto measures_h = Omega_h::HostRead(measures); + const auto elem_verts_h = + Omega_h::HostRead(mesh.ask_elem_verts()); + const int verts_per_elem = mesh.dim() + 1; + + double integral = 0.0; + for (Omega_h::LO e = 0; e < mesh.nelems(); ++e) { + double avg = 0.0; + for (int k = 0; k < verts_per_elem; ++k) { + avg += values[elem_verts_h[verts_per_elem * e + k]]; + } + integral += measures_h[e] * (avg / verts_per_elem); + } + return integral; +} +#endif // PCMS_ENABLE_OMEGA_H + // Copy coordinates from device memory to a host view. // This handles potential layout mismatches between host and device memory // spaces. @@ -86,14 +189,8 @@ inline Kokkos::View CopyCoordinatesToHost( auto coords_view = Kokkos::View("coords_view", nents, dim); auto coords_view_device = - Kokkos::create_mirror(DeviceMemorySpace(), coords_view); - Kokkos::parallel_for( - "copy_coords_to_host_view", Kokkos::RangePolicy<>(0, nents), - KOKKOS_LAMBDA(int i) { - for (int d = 0; d < dim; ++d) { - coords_view_device(i, d) = coords_device(i, d); - } - }); + Kokkos::create_mirror_view(DeviceMemorySpace(), coords_view); + ConvertMismatchLayoutView2D(coords_view_device, coords_device); Kokkos::deep_copy(coords_view, coords_view_device); return coords_view; } @@ -122,44 +219,49 @@ inline std::vector EvaluateReferenceFunction(const std::vector& pts, expected_host.data() + expected_host.extent(0)); } -template +template struct SetFieldFunctor { - Kokkos::View data; - Kokkos::View coords; + DataView data; + CoordsView coords; Func f; - SetFieldFunctor(Kokkos::View data_, - Kokkos::View coords_, Func f_) + SetFieldFunctor(DataView data_, CoordsView coords_, Func f_) : data(data_), coords(coords_), f(f_) { } KOKKOS_INLINE_FUNCTION - void operator()(int i) const { data(i) = f(coords(i, 0), coords(i, 1)); } + void operator()(int i) const + { + if constexpr (std::is_invocable_v) { + data(i) = f(coords(i, 0), coords(i, 1), coords(i, 2)); + } else { + data(i) = f(coords(i, 0), coords(i, 1)); + } + } }; -// Set scalar DOF data by sampling func at each DOF-holder coordinate. +// Set scalar DOF data by sampling func at each DOF-holder coordinate. The +// arity of func selects the spatial dimension: func(x, y) for 2D layouts, +// func(x, y, z) for 3D. template inline void SetField(const FieldLayout& layout, FieldData& field, Func func) { using MemorySpace = typename ExecutionSpace::memory_space; + static_assert(std::is_invocable_v || + std::is_invocable_v, + "SetField requires func(x, y) or func(x, y, z)"); + auto dof_coords = layout.GetDOFHolderCoordinates().GetValues(); int n = static_cast(dof_coords.extent(0)); - Kokkos::View coords_device("coords_device", n); - Kokkos::parallel_for( - "field_test_utils_copy_coords", Kokkos::RangePolicy(0, n), - KOKKOS_LAMBDA(int i) { - coords_device(i, 0) = dof_coords(i, 0); - coords_device(i, 1) = dof_coords(i, 1); - }); Kokkos::View data_device("data_device", n); Kokkos::parallel_for("field_test_utils_set_field", Kokkos::RangePolicy(0, n), - SetFieldFunctor{data_device, coords_device, func}); + SetFieldFunctor{data_device, dof_coords, func}); auto data_host = Kokkos::create_mirror_view_and_copy(HostMemorySpace(), data_device); @@ -253,12 +355,8 @@ void CheckEvaluation(const PointEvaluator& evaluator, { int n = static_cast(pts.size()) / 2; - Kokkos::View out_device("out_device", n); - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - Rank2View out(out_device.data(), n, 1); - // Rank2View out(eval.data(), n, 1); - evaluator.Evaluate(field, out); + Kokkos::View out_device("out_device", n, 1); + evaluator.Evaluate(field, MakeRank2View(out_device)); auto out_host = Kokkos::create_mirror_view_and_copy(HostMemorySpace(), out_device); @@ -267,8 +365,8 @@ void CheckEvaluation(const PointEvaluator& evaluator, for (int i = 0; i < n; ++i) { INFO("Point " << i << " (" << pts[2 * static_cast(i)] << ", " << pts[2 * static_cast(i) + 1] << ")" - << " got=" << out_host(i) << " expected=" << expected[i]); - REQUIRE(out_host(i) == Catch::Approx(expected[i]).margin(abs_tol)); + << " got=" << out_host(i, 0) << " expected=" << expected[i]); + REQUIRE(out_host(i, 0) == Catch::Approx(expected[i]).margin(abs_tol)); } } @@ -295,16 +393,13 @@ inline void CheckFillMode(const PointEvaluator& evaluator, { int n = static_cast(outside_pts.size()) / 2; - Kokkos::View out_device("out_device", n); - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - Rank2View out(out_device.data(), n, 1); - evaluator.Evaluate(field, out); + Kokkos::View out_device("out_device", n, 1); + evaluator.Evaluate(field, MakeRank2View(out_device)); auto out_host = Kokkos::create_mirror_view_and_copy(HostMemorySpace(), out_device); for (int i = 0; i < n; ++i) { - REQUIRE(out_host(i) == fill_value); + REQUIRE(out_host(i, 0) == fill_value); } } @@ -339,11 +434,8 @@ void CheckEvaluationWithFill(const Factory& factory, const Field& field, auto evaluator = factory->template CreatePointEvaluator( EvaluationRequest::FromCoordinates(device_coords.coordinate_view, policy)); - Kokkos::View out_device("out_device", n); - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - Rank2View out(out_device.data(), n, 1); - evaluator->Evaluate(field, out); + Kokkos::View out_device("out_device", n, 1); + evaluator->Evaluate(field, MakeRank2View(out_device)); auto out_host = Kokkos::create_mirror_view_and_copy(HostMemorySpace(), out_device); @@ -352,16 +444,34 @@ void CheckEvaluationWithFill(const Factory& factory, const Field& field, for (int i = 0; i < n; ++i) { INFO("Point " << i << " (" << pts[2 * static_cast(i)] << ", " << pts[2 * static_cast(i) + 1] << ")" - << " got=" << out_host(i)); + << " got=" << out_host(i, 0)); if (is_inside[i]) { INFO(" expected=" << expected[i]); - REQUIRE(out_host(i) == Catch::Approx(expected[i]).margin(abs_tol)); + REQUIRE(out_host(i, 0) == Catch::Approx(expected[i]).margin(abs_tol)); } else { - REQUIRE(out_host(i) == fill_value); + REQUIRE(out_host(i, 0) == fill_value); } } } +#if defined(PCMS_ENABLE_PETSC) && defined(PCMS_ENABLE_MESHFIELDS) +// Evaluates source_field at the integrator's sample points and assembles the +// load vector. +inline void EvaluateAndAssemble( + LinearFormIntegrator& integrator, + const std::shared_ptr& source_space, + const Field& source_field) +{ + const auto& pts = integrator.GetIntegrationPoints(); + const std::size_t npts = pts.GetValues().extent(0); + auto evaluator = source_space->CreatePointEvaluator( + EvaluationRequest::FromCoordinates(pts)); + Kokkos::View sampled("sampled", npts, 1); + evaluator->Evaluate(source_field, MakeRank2View(sampled)); + integrator.Assemble(MakeConstRank2View(sampled)); +} +#endif // PCMS_ENABLE_PETSC && PCMS_ENABLE_MESHFIELDS + } // namespace pcms::test #endif // PCMS_TEST_FIELD_TEST_UTILS_H diff --git a/test/test_eqdsk.cpp b/test/test_eqdsk.cpp index d546fbbe..c51f46b1 100644 --- a/test/test_eqdsk.cpp +++ b/test/test_eqdsk.cpp @@ -7,6 +7,7 @@ #include #include #include +#include "field_test_utils.h" #include #include @@ -212,42 +213,24 @@ TEST_CASE("EQDSKData with SplineFunctionSpace") const int num_eval_points = 3; - auto eval_coords_host = Kokkos::View( - "eval_coords_host", num_eval_points, 2); - for (size_t i = 0; i < num_eval_points; ++i) { - eval_coords_host(i, 0) = eval_coords[2 * i]; // R - eval_coords_host(i, 1) = eval_coords[2 * i + 1]; // Z - } - - auto eval_coords_device = - Kokkos::View("eval_coords_device", - num_eval_points, 2); - pcms::DeepCopyMismatchLayouts(eval_coords_device, eval_coords_host); - - auto coords_view = pcms::MakeRank2View(eval_coords_device); - auto coord_view = pcms::CoordinateView{ - CoordinateSystem::Cartesian, coords_view}; - auto eval_request = pcms::EvaluationRequest::FromCoordinates(coord_view); + auto device_coords = pcms::test::CreateDeviceCoordinateView( + eval_coords, CoordinateSystem::Cartesian); + auto eval_request = + pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view); auto evaluator = spline_space->CreatePointEvaluator(eval_request); - auto eval_results_1d = Kokkos::View( - "eval_results", num_eval_points); - - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - pcms::Rank2View - eval_results(eval_results_1d.data(), num_eval_points, 1); - - evaluator->Evaluate(psi_field, eval_results); + Kokkos::View eval_results( + "eval_results", num_eval_points, 1); + evaluator->Evaluate(psi_field, pcms::MakeRank2View(eval_results)); auto results_host = Kokkos::create_mirror_view_and_copy( - pcms::HostMemorySpace(), eval_results_1d); + pcms::HostMemorySpace(), eval_results); // Verify that all results are finite for (int i = 0; i < num_eval_points; ++i) { - REQUIRE(std::isfinite(results_host(i))); + REQUIRE(std::isfinite(results_host(i, 0))); } } } diff --git a/test/test_field_evaluation.cpp b/test/test_field_evaluation.cpp index ff238e11..85af38b4 100644 --- a/test/test_field_evaluation.cpp +++ b/test/test_field_evaluation.cpp @@ -34,10 +34,10 @@ TEST_CASE("evaluate linear 2d omega_h_field") pcms::test::SetField( field.GetData(), *factory->GetLayout(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); pcms::test::CheckEvaluation( factory, field, pcms::test::StandardEvalCoords2D(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); } #ifdef PCMS_ENABLE_MESHFIELDS @@ -50,35 +50,11 @@ TEST_CASE("evaluate quadratic 2d meshfields_field") mesh, 2, 1, pcms::CoordinateSystem::Cartesian, "global", pcms::LagrangeFunctionSpace::Backend::MeshFields); - // Quadratic DOF holders span vertices and edge midpoints; set them inline. - const auto nverts = mesh.nents(0); - const auto nedges = mesh.nents(1); - auto mesh_coords = mesh.coords(); - auto edge_verts = mesh.ask_verts_of(1); - - Omega_h::Write test_f(nverts + nedges); - Omega_h::parallel_for( - nverts, OMEGA_H_LAMBDA(int i) { - test_f[i] = sin_f(mesh_coords[2 * static_cast(i)], - mesh_coords[2 * static_cast(i) + 1]); - }); - Omega_h::parallel_for( - nedges, OMEGA_H_LAMBDA(int i) { - auto ep = Omega_h::gather_verts<2>(edge_verts, i); - Real cx = (mesh_coords[2 * static_cast(ep[0])] + - mesh_coords[2 * static_cast(ep[1])]) / - 2; - Real cy = (mesh_coords[2 * static_cast(ep[0]) + 1] + - mesh_coords[2 * static_cast(ep[1]) + 1]) / - 2; - test_f[nverts + i] = sin_f(cx, cy); - }); - - Omega_h::HostWrite test_f_host(test_f); + // Quadratic DOF holders span vertices and edge midpoints; the layout's DOF + // coordinates cover both, so SetField samples sin_f at every holder. auto field = factory->CreateFunction(); - field.GetData().SetDOFHolderDataHost( - pcms::Rank2View(test_f_host.data(), - test_f_host.size(), 1)); + pcms::test::SetField( + field, OMEGA_H_LAMBDA(Real x, Real y) { return sin_f(x, y); }); pcms::test::CheckEvaluation( factory, field, kEvalCoords, diff --git a/test/test_field_exchange_planner.cpp b/test/test_field_exchange_planner.cpp index 61b5bc54..99e6181b 100644 --- a/test/test_field_exchange_planner.cpp +++ b/test/test_field_exchange_planner.cpp @@ -51,9 +51,7 @@ TEST_CASE("GID messages insert headers around compact field payloads", "gid_message", plan.msg_size + 2 * static_cast(pcms::ent_offsets_len)); pcms::GenericFieldExchangePlanner planner; - planner.FillGidMessage(layout, plan, - pcms::Rank1View( - message.data(), message.size())); + planner.FillGidMessage(layout, plan, pcms::MakeRank1View(message)); const pcms::GO expected[] = { 0, 2, 2, 2, 2, 0, 2, 0, 2, 2, 2, 2, 1, 3, diff --git a/test/test_field_interpolation.cpp b/test/test_field_interpolation.cpp index bc88ab97..f8e2432d 100644 --- a/test/test_field_interpolation.cpp +++ b/test/test_field_interpolation.cpp @@ -55,35 +55,11 @@ TEST_CASE("interpolate quadratic 2d meshfields_field") auto factory2 = pcms::LagrangeFunctionSpace::FromMesh( mesh, 2, 1, pcms::CoordinateSystem::Cartesian, "global", pcms::LagrangeFunctionSpace::Backend::MeshFields); - auto layout = factory2->GetLayout(); - const auto nverts = mesh.nents(0); - const auto nedges = mesh.nents(1); - auto mesh_coords = mesh.coords(); - auto edge_verts = mesh.ask_verts_of(1); - Omega_h::Write test_f(nverts + nedges); - Omega_h::parallel_for( - nverts, OMEGA_H_LAMBDA(int i) { - Real x = mesh_coords[2 * i + 0]; - Real y = mesh_coords[2 * i + 1]; - test_f[i] = interpolation_linear_f(x, y); - }); - Omega_h::parallel_for( - nedges, OMEGA_H_LAMBDA(int i) { - auto endpoints = Omega_h::gather_verts<2>(edge_verts, i); - Real x0 = mesh_coords[2 * endpoints[0] + 0]; - Real y0 = mesh_coords[2 * endpoints[0] + 1]; - Real x1 = mesh_coords[2 * endpoints[1] + 0]; - Real y1 = mesh_coords[2 * endpoints[1] + 1]; - Real cx = (x0 + x1) / 2; - Real cy = (y0 + y1) / 2; - test_f[nverts + i] = interpolation_linear_f(cx, cy); - }); - - Omega_h::HostWrite test_f_host(test_f); auto field = factory2->CreateFunction(); auto interpolated = factory2->CreateFunction(); - field.SetDOFHolderDataHost(pcms::Rank2View( - test_f_host.data(), test_f_host.size(), 1)); + pcms::test::SetField( + field, + OMEGA_H_LAMBDA(Real x, Real y) { return interpolation_linear_f(x, y); }); pcms::Interpolator interp(*factory2, *factory2); interp.Apply(field, interpolated); diff --git a/test/test_interpolation_class.cpp b/test/test_interpolation_class.cpp index 7a34af4d..878f6288 100644 --- a/test/test_interpolation_class.cpp +++ b/test/test_interpolation_class.cpp @@ -10,36 +10,10 @@ #include #include #include +#include "field_test_utils.h" #include #include -#include - -bool areArraysEqualUnordered(const Omega_h::HostRead& array1, - const Omega_h::HostRead& array2, - int start, int end) -{ - // Ensure the indices are valid - assert(start >= 0 && end <= array1.size() && start <= end); - assert(start >= 0 && end <= array2.size() && start <= end); - - // Use frequency maps to count occurrences of each value - std::unordered_map freq1, freq2; - - for (int i = start; i < end; ++i) { - freq1[array1[i]]++; - freq2[array2[i]]++; - } - - // Compare the frequency maps - if (freq1 != freq2) { - pcms::printError("[ERROR] Arrays differ in the range [%d, %d)\n", start, - end); - return false; - } - - return true; -} void translate_mesh(Omega_h::Mesh* mesh, Omega_h::Vector<2> translation_vector) { @@ -135,28 +109,13 @@ TEST_CASE("Test MLSMeshInterpolation") auto mls_single = pcms::MLSMeshInterpolation(source_mesh, 0.12, 15, 3, true, 0.0, 5.0); - auto source_points_reals = - pcms::get_entity_centroids(source_mesh, Omega_h::FACE); - auto source_points_host = - Omega_h::HostRead(source_points_reals); - auto source_points_host_write = - Omega_h::HostWrite(source_points_host.size()); - for (int i = 0; i < source_points_host.size(); i++) { - source_points_host_write[i] = source_points_host[i]; - } - auto source_points_view = pcms::Rank1View( - source_points_host_write.data(), source_points_host_write.size()); - - auto target_points_reals = source_mesh.coords(); - auto target_points_host = - Omega_h::HostRead(target_points_reals); - auto target_points_host_write = - Omega_h::HostWrite(target_points_host.size()); - for (int i = 0; i < target_points_host.size(); i++) { - target_points_host_write[i] = target_points_host[i]; - } - auto target_points_view = pcms::Rank1View( - target_points_host_write.data(), target_points_host_write.size()); + auto source_points_vec = pcms::test::CopyOmegaHRealsToVector( + pcms::get_entity_centroids(source_mesh, Omega_h::FACE)); + auto source_points_view = pcms::make_array_view(source_points_vec); + + auto target_points_vec = + pcms::test::CopyOmegaHRealsToVector(source_mesh.coords()); + auto target_points_view = pcms::make_array_view(target_points_vec); REQUIRE(source_mesh.dim() == 2); pcms::printInfo("Point cloud based search...\n"); auto point_mls = pcms::MLSPointCloudInterpolation( @@ -173,14 +132,11 @@ TEST_CASE("Test MLSMeshInterpolation") source_mesh.nverts()); Omega_h::HostWrite exact_values_at_nodes(source_sinxcosy_node); - pcms::Rank1View sourceArrayView( - source_data_host_write.data(), source_data_host_write.size()); - pcms::Rank1View interpolatedArrayView( - interpolated_data_hwrite.data(), interpolated_data_hwrite.size()); - pcms::Rank1View - point_cloud_interpolatedArrayView( - point_cloud_interpolated_data_hwrite.data(), - point_cloud_interpolated_data_hwrite.size()); + auto sourceArrayView = pcms::make_array_view(source_data_host_write); + auto interpolatedArrayView = + pcms::make_array_view(interpolated_data_hwrite); + auto point_cloud_interpolatedArrayView = + pcms::make_array_view(point_cloud_interpolated_data_hwrite); OMEGA_H_CHECK_PRINTF(sourceArrayView.size() == mls_single.getSourceSize(), "Source size mismatch: %zu vs %zu\n", @@ -222,33 +178,8 @@ TEST_CASE("Test MLSMeshInterpolation") ///*****************************// auto mesh_based_supports = mls_single.getSupports(); auto point_cloud_based_supports = point_mls.getSupports(); - auto mesh_based_support_ptr_host = - Omega_h::HostRead(mesh_based_supports.supports_ptr); - auto mesh_based_support_idx_host = - Omega_h::HostRead(mesh_based_supports.supports_idx); - auto point_cloud_based_support_ptr_host = - Omega_h::HostRead(point_cloud_based_supports.supports_ptr); - auto point_cloud_based_support_idx_host = - Omega_h::HostRead(point_cloud_based_supports.supports_idx); - - REQUIRE(point_cloud_based_support_idx_host.size() == - mesh_based_support_idx_host.size()); - REQUIRE(point_cloud_based_support_ptr_host.size() == - mesh_based_support_ptr_host.size()); - for (int i = 0; i < mesh_based_support_ptr_host.size(); i++) { - REQUIRE(point_cloud_based_support_ptr_host[i] == - mesh_based_support_ptr_host[i]); - } - - for (int i = 0; i < mesh_based_support_ptr_host.size() - 1; i++) { - auto start = mesh_based_support_ptr_host[i]; - auto end = mesh_based_support_ptr_host[i + 1]; - - bool isEqual = - areArraysEqualUnordered(mesh_based_support_idx_host, - point_cloud_based_support_idx_host, start, end); - REQUIRE(isEqual); - } + pcms::test::CheckSupportResultsEquivalent(point_cloud_based_supports, + mesh_based_supports); // Check if the point cloud interpolation is same as the MLS interpolation pcms::printDebugInfo("Interpolated data size: %d\n", @@ -289,10 +220,9 @@ TEST_CASE("Test MLSMeshInterpolation") Omega_h::HostWrite interpolated_data_hwrite( mls_double.getTargetSize()); - pcms::Rank1View sourceArrayView( - source_data_host_write.data(), source_data_host_write.size()); - pcms::Rank1View interpolatedArrayView( - interpolated_data_hwrite.data(), interpolated_data_hwrite.size()); + auto sourceArrayView = pcms::make_array_view(source_data_host_write); + auto interpolatedArrayView = + pcms::make_array_view(interpolated_data_hwrite); mls_double.eval(sourceArrayView, interpolatedArrayView); @@ -326,10 +256,8 @@ TEST_CASE("MLSPointCloudInterpolation honors provided dimension in eval") } auto target_points = source_points; - auto source_points_view = pcms::Rank1View( - source_points.data(), source_points.size()); - auto target_points_view = pcms::Rank1View( - target_points.data(), target_points.size()); + auto source_points_view = pcms::make_array_view(source_points); + auto target_points_view = pcms::make_array_view(target_points); // Degree-1 polynomial that depends on z to catch accidental 2D behavior. auto source_values = Omega_h::HostWrite(27); @@ -341,10 +269,8 @@ TEST_CASE("MLSPointCloudInterpolation honors provided dimension in eval") } auto output_values = Omega_h::HostWrite(27, "output_values"); - auto source_values_view = pcms::Rank1View( - source_values.data(), source_values.size()); - auto output_values_view = pcms::Rank1View( - output_values.data(), output_values.size()); + auto source_values_view = pcms::make_array_view(source_values); + auto output_values_view = pcms::make_array_view(output_values); auto mls = pcms::MLSPointCloudInterpolation( source_points_view, target_points_view, 3, 2.5, 10, 1, true, 0.0, 5.0); diff --git a/test/test_interpolation_on_ltx_mesh.cpp b/test/test_interpolation_on_ltx_mesh.cpp index a348acae..141ebfa8 100644 --- a/test/test_interpolation_on_ltx_mesh.cpp +++ b/test/test_interpolation_on_ltx_mesh.cpp @@ -83,16 +83,13 @@ TEST_CASE("Test Interpolation on LTX Mesh", "[interpolation]") printf("[INFO] Degas2 Mesh loaded from %s with %d elements\n", degas2_mesh_filename.c_str(), degas2_num_elems); const auto degas2_mesh_centroids_view = - pcms::Rank1View( - degas2_mesh_centroids_host.data(), degas2_mesh_centroids_host.size()); + pcms::make_array_view(degas2_mesh_centroids_host); auto xgc_mesh_points = read_xgc_mesh_nodes(ltx_mesh_base_filename + ".node"); const int xgc_num_nodes = xgc_mesh_points.size() / 2; printf("[INFO] XGC Mesh loaded from %s with %d points\n", ltx_mesh_base_filename.c_str(), xgc_num_nodes); - const auto xgc_mesh_points_view = - pcms::Rank1View( - xgc_mesh_points.data(), xgc_mesh_points.size()); + const auto xgc_mesh_points_view = pcms::make_array_view(xgc_mesh_points); auto xgc_to_degas2_interpolator = pcms::MLSPointCloudInterpolation( xgc_mesh_points_view, degas2_mesh_centroids_view, 2, 0.000001, 10, 1, true, @@ -120,18 +117,13 @@ TEST_CASE("Test Interpolation on LTX Mesh", "[interpolation]") // ------------------ First Interpolation ------------------ // const auto density_at_xgc_nodes_view = - pcms::Rank1View(density_at_xgc_nodes.data(), - density_at_xgc_nodes.size()); - const auto temp_at_xgc_nodes_view = - pcms::Rank1View(temp_at_xgc_nodes.data(), - temp_at_xgc_nodes.size()); + pcms::make_array_view(density_at_xgc_nodes); + const auto temp_at_xgc_nodes_view = pcms::make_array_view(temp_at_xgc_nodes); const auto density_at_degas2_centroids_view = - pcms::Rank1View( - density_at_degas2_centroids.data(), density_at_degas2_centroids.size()); + pcms::make_array_view(density_at_degas2_centroids); const auto temp_at_degas2_centroids_view = - pcms::Rank1View( - temp_at_degas2_centroids.data(), temp_at_degas2_centroids.size()); + pcms::make_array_view(temp_at_degas2_centroids); Omega_h::HostWrite interpolated_xgc_density(degas2_num_elems); Omega_h::HostWrite interpolated_xgc_temp(degas2_num_elems); @@ -139,17 +131,13 @@ TEST_CASE("Test Interpolation on LTX Mesh", "[interpolation]") Omega_h::HostWrite interpolated_degas2_temp(xgc_num_nodes); const auto interpolated_xgc_density_view = - pcms::Rank1View( - interpolated_xgc_density.data(), interpolated_xgc_density.size()); + pcms::make_array_view(interpolated_xgc_density); const auto interpolated_xgc_temp_view = - pcms::Rank1View( - interpolated_xgc_temp.data(), interpolated_xgc_temp.size()); + pcms::make_array_view(interpolated_xgc_temp); const auto interpolated_degas2_density_view = - pcms::Rank1View( - interpolated_degas2_density.data(), interpolated_degas2_density.size()); + pcms::make_array_view(interpolated_degas2_density); const auto interpolated_degas2_temp_view = - pcms::Rank1View( - interpolated_degas2_temp.data(), interpolated_degas2_temp.size()); + pcms::make_array_view(interpolated_degas2_temp); xgc_to_degas2_interpolator.eval(density_at_xgc_nodes_view, interpolated_xgc_density_view); @@ -174,21 +162,13 @@ TEST_CASE("Test Interpolation on LTX Mesh", "[interpolation]") degas2_num_elems); const auto interpolated_back_density_at_xgc_nodes_view = - pcms::Rank1View( - interpolated_back_density_at_xgc_nodes.data(), - interpolated_back_density_at_xgc_nodes.size()); + pcms::make_array_view(interpolated_back_density_at_xgc_nodes); const auto interpolated_back_temp_at_xgc_nodes_view = - pcms::Rank1View( - interpolated_back_temp_at_xgc_nodes.data(), - interpolated_back_temp_at_xgc_nodes.size()); + pcms::make_array_view(interpolated_back_temp_at_xgc_nodes); const auto interpolated_back_density_at_degas2_centroids_view = - pcms::Rank1View( - interpolated_back_density_at_degas2_centroids.data(), - interpolated_back_density_at_degas2_centroids.size()); + pcms::make_array_view(interpolated_back_density_at_degas2_centroids); const auto interpolated_back_temp_at_degas2_centroids_view = - pcms::Rank1View( - interpolated_back_temp_at_degas2_centroids.data(), - interpolated_back_temp_at_degas2_centroids.size()); + pcms::make_array_view(interpolated_back_temp_at_degas2_centroids); degas2_to_xgc_interpolator.eval(interpolated_xgc_density_view, interpolated_back_density_at_xgc_nodes_view); diff --git a/test/test_intersections.cpp b/test/test_intersections.cpp index 02598f0a..36ccb766 100644 --- a/test/test_intersections.cpp +++ b/test/test_intersections.cpp @@ -111,7 +111,7 @@ TEST_CASE("Mesh intersection test with source and target", "[intersection]") Omega_h::parallel_for( ntgt, OMEGA_H_LAMBDA(int t) { auto tgt_vert_coords = - get_vert_coords_of_elem(tgt_coords, tgt_faces2verts, t); + get_vert_coords_of_elem<2>(tgt_coords, tgt_faces2verts, t); int start = intersection.tgt2src_offsets[t]; int end = intersection.tgt2src_offsets[t + 1]; @@ -120,7 +120,7 @@ TEST_CASE("Mesh intersection test with source and target", "[intersection]") for (int i = start; i < end; ++i) { int sid = intersection.tgt2src_indices[i]; auto src_vert_coords = - get_vert_coords_of_elem(src_coords, src_faces2verts, sid); + get_vert_coords_of_elem<2>(src_coords, src_faces2verts, sid); r3d::Polytope<2> poly; r3d::intersect_simplices(poly, tgt_vert_coords, src_vert_coords); diff --git a/test/test_localization_factory.cpp b/test/test_localization_factory.cpp index 571ac175..d65e8fc0 100644 --- a/test/test_localization_factory.cpp +++ b/test/test_localization_factory.cpp @@ -14,6 +14,7 @@ #include "pcms/localization/point_cloud_localization.h" #include "pcms/discretization/discretization/omega_h.hpp" #include "pcms/utility/arrays.h" +#include "pcms/utility/omega_h_array_utils.h" #include "field_test_utils.h" namespace @@ -46,17 +47,7 @@ TEST_CASE( auto target_device = pcms::test::CreateDeviceCoordinateView( target_coords, pcms::CoordinateSystem::Cartesian, dim); - auto coords_read = Omega_h::HostRead(source_coords); - - // auto coords_dev = Kokkos::create_mirror_view_and_copy( - // Kokkos::DefaultExecutionSpace{}, coords_host); - auto coords_dev = Kokkos::View( - "point_cloud_coords", mesh.nverts(), dim); - auto coords_host = Kokkos::create_mirror(pcms::HostMemorySpace(), coords_dev); - for (int i = 0; i < mesh.nverts(); ++i) - for (int d = 0; d < dim; ++d) - coords_host(i, d) = coords_read[i * dim + d]; - Kokkos::deep_copy(coords_dev, coords_host); + auto coords_dev = pcms::ConvertCoordsTo2D(source_coords, mesh.nverts(), dim); auto layout = std::make_shared( dim, coords_dev, pcms::CoordinateSystem::Cartesian); diff --git a/test/test_mesh_intersection_field_transfer.cpp b/test/test_mesh_intersection_field_transfer.cpp index fbb34ddb..a627ab81 100644 --- a/test/test_mesh_intersection_field_transfer.cpp +++ b/test/test_mesh_intersection_field_transfer.cpp @@ -7,53 +7,14 @@ #include #include +#include #include "field_test_utils.h" -#include #include namespace { -// 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; -} - -std::shared_ptr MakeP1Space( - Omega_h::Mesh& mesh, const std::string& global_id_name = "global") -{ - return pcms::LagrangeFunctionSpace::FromMesh( - mesh, 1, 1, pcms::CoordinateSystem::Cartesian, global_id_name, - pcms::LagrangeFunctionSpace::Backend::OmegaH); -} - -std::shared_ptr MakeP0Space(Omega_h::Mesh& mesh) -{ - return pcms::LagrangeFunctionSpace::FromMesh( - mesh, 0, 1, pcms::CoordinateSystem::Cartesian, "global", - pcms::LagrangeFunctionSpace::Backend::OmegaH); -} - void AddReorderedVertexGlobalIds(Omega_h::Mesh& mesh) { mesh.add_tag( @@ -68,64 +29,6 @@ void AddSparseVertexGlobalIds(Omega_h::Mesh& mesh) 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, - const pcms::Field& field) -{ - const auto values = pcms::FlattenToRank1View(field.GetDOFHolderDataHost()); - const auto elem_areas = Omega_h::measure_elements_real(&mesh); - 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; -} - -// 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); - - 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 avg = (values[v0] + values[v1] + values[v2]) / 3.0; - integral += elem_areas_h[e] * avg; - } - return integral; -} - } // namespace // Fields that already live in the target order-1 space (constants and affine @@ -139,12 +42,12 @@ TEST_CASE("OmegaHConservativeProjection reproduces constant and linear fields", // 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})); + pcms::test::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})); + pcms::test::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_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh); auto source = source_space->CreateFunction(); auto target = target_space->CreateFunction(); @@ -166,8 +69,9 @@ TEST_CASE("OmegaHConservativeProjection reproduces constant and linear fields", 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(IntegrateP1Field(source_mesh, source)).margin(1e-10)); + REQUIRE(pcms::test::IntegrateP1Field(target_mesh, target) == + Catch::Approx(pcms::test::IntegrateP1Field(source_mesh, source)) + .margin(1e-10)); } SECTION("linear field is reproduced on target vertices and conserved") @@ -179,34 +83,31 @@ TEST_CASE("OmegaHConservativeProjection reproduces constant and linear fields", const auto target_values = pcms::FlattenToRank1View(target.GetDOFHolderDataHost()); - const auto tgt_coords_h = - Omega_h::HostRead(target_mesh.coords()); + const auto tgt_coords_h = pcms::test::CopyCoordinatesToHost( + pcms::MakeConstRank2View(target_mesh.coords(), 2), target_mesh.nverts(), + 2); 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]; + const double expected = tgt_coords_h(i, 0) + tgt_coords_h(i, 1); REQUIRE(target_values[i] == Catch::Approx(expected).margin(1e-9)); } - REQUIRE(IntegrateP1Field(target_mesh, target) == - Catch::Approx(IntegrateP1Field(source_mesh, source)).margin(1e-9)); + REQUIRE(pcms::test::IntegrateP1Field(target_mesh, target) == + Catch::Approx(pcms::test::IntegrateP1Field(source_mesh, source)) + .margin(1e-9)); } } -// 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::Mesh source_mesh = - BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3})); + pcms::test::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})); + pcms::test::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_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh); auto source = source_space->CreateFunction(); auto target = target_space->CreateFunction(); @@ -218,8 +119,9 @@ TEST_CASE("OmegaHConservativeProjection conserves the integral", 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)); + REQUIRE(pcms::test::IntegrateP1Field(target_mesh, target) == + Catch::Approx(pcms::test::IntegrateP1Field(source_mesh, source)) + .margin(1e-12)); // Second Apply with a different source field — verifies cached state is // reused correctly and is not stale. @@ -227,8 +129,9 @@ TEST_CASE("OmegaHConservativeProjection conserves the integral", 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)); + REQUIRE(pcms::test::IntegrateP1Field(target_mesh, target) == + Catch::Approx(pcms::test::IntegrateP1Field(source_mesh, source)) + .margin(1e-12)); } TEST_CASE("OmegaHConservativeProjection writes reordered target GIDs in local " @@ -238,13 +141,13 @@ TEST_CASE("OmegaHConservativeProjection writes reordered target GIDs in local " Omega_h::Library lib; Omega_h::Mesh source_mesh = - BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3})); + pcms::test::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})); + pcms::test::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_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh, "reordered_global"); auto source = source_space->CreateFunction(); auto target = target_space->CreateFunction(); @@ -277,13 +180,13 @@ TEST_CASE("OmegaHConservativeProjection maps sparse target GIDs to active " Omega_h::Library lib; Omega_h::Mesh source_mesh = - BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3})); + pcms::test::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})); + pcms::test::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_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh, "sparse_global"); auto source = source_space->CreateFunction(); auto target = target_space->CreateFunction(); @@ -318,12 +221,12 @@ TEST_CASE("OmegaHConservativeProjection P0 source to P1 target", Omega_h::Library lib; Omega_h::Mesh source_mesh = - BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3})); + pcms::test::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})); + pcms::test::BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3})); - auto source_space = MakeP0Space(source_mesh); - auto target_space = MakeP1Space(target_mesh); + auto source_space = pcms::test::MakeP0Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh); auto source = source_space->CreateFunction(); auto target = target_space->CreateFunction(); @@ -345,8 +248,9 @@ TEST_CASE("OmegaHConservativeProjection P0 source to P1 target", 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)); + REQUIRE(pcms::test::IntegrateP1Field(target_mesh, target) == + Catch::Approx(pcms::test::IntegrateP0Field(source_mesh, source)) + .margin(1e-10)); } SECTION("non-constant P0 source conserves the integral") @@ -357,8 +261,9 @@ TEST_CASE("OmegaHConservativeProjection P0 source to P1 target", projection.Apply(source, target); - REQUIRE(IntegrateP1Field(target_mesh, target) == - Catch::Approx(IntegrateP0Field(source_mesh, source)).margin(1e-10)); + REQUIRE(pcms::test::IntegrateP1Field(target_mesh, target) == + Catch::Approx(pcms::test::IntegrateP0Field(source_mesh, source)) + .margin(1e-10)); } } @@ -372,12 +277,12 @@ TEST_CASE("OmegaHConservativeProjection P1 source to P0 target", Omega_h::Library lib; Omega_h::Mesh source_mesh = - BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3})); + pcms::test::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})); + pcms::test::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_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP0Space(target_mesh); auto source = source_space->CreateFunction(); auto target = target_space->CreateFunction(); @@ -399,8 +304,9 @@ TEST_CASE("OmegaHConservativeProjection P1 source to P0 target", 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)); + REQUIRE(pcms::test::IntegrateP0Field(target_mesh, target) == + Catch::Approx(pcms::test::IntegrateP1Field(source_mesh, source)) + .margin(1e-10)); } SECTION("linear field projects to cell average and conserves the integral") @@ -414,13 +320,16 @@ TEST_CASE("OmegaHConservativeProjection P1 source to P0 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]); + const auto centroids_h = Omega_h::HostRead( + pcms::get_entity_centroids(target_mesh, target_mesh.dim())); + REQUIRE(static_cast(target_values.size()) == + target_mesh.nelems()); + for (Omega_h::LO e = 0; e < target_mesh.nelems(); ++e) { + const double expected = f(centroids_h[2 * e + 0], centroids_h[2 * 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)); + REQUIRE(pcms::test::IntegrateP0Field(target_mesh, target) == + Catch::Approx(pcms::test::IntegrateP1Field(source_mesh, source)) + .margin(1e-9)); } } diff --git a/test/test_omega_h_3d_conservative_projection.cpp b/test/test_omega_h_3d_conservative_projection.cpp new file mode 100644 index 00000000..d2783ae3 --- /dev/null +++ b/test/test_omega_h_3d_conservative_projection.cpp @@ -0,0 +1,227 @@ +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include "field_test_utils.h" + +#include + +namespace +{ + +// Build a tetrahedral box mesh of [0,1]^3 at the given resolution. build_box +// produces a classified simplex mesh, which is what the Omega_h Lagrange +// function space requires. +Omega_h::Mesh BuildUnitCube(Omega_h::Library& lib, int n) +{ + return Omega_h::build_box(lib.world(), OMEGA_H_SIMPLEX, 1.0, 1.0, 1.0, n, n, + n); +} + +} // namespace + +TEST_CASE("OmegaHConservativeProjection (3D tets) reproduces constant and " + "linear fields", + "[transfer][mesh_intersection][3d]") +{ + Omega_h::Library lib; + + // Two independent tessellations of the same unit cube. + Omega_h::Mesh source_mesh = BuildUnitCube(lib, 1); + Omega_h::Mesh target_mesh = BuildUnitCube(lib, 2); + + auto source_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh); + + auto source = source_space->CreateFunction(); + auto target = target_space->CreateFunction(); + + pcms::OmegaHConservativeProjection projection(*source_space, *target_space); + + SECTION("constant field is preserved and conserved") + { + const double c = 2.0; + pcms::test::SetField( + source, KOKKOS_LAMBDA(pcms::Real, 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-9)); + } + REQUIRE(pcms::test::IntegrateP1Field(target_mesh, target) == + Catch::Approx(pcms::test::IntegrateP1Field(source_mesh, source)) + .margin(1e-9)); + } + + SECTION("linear field is reproduced on target vertices and conserved") + { + pcms::test::SetField( + source, KOKKOS_LAMBDA(pcms::Real x, pcms::Real y, pcms::Real z) { + return 1.0 + x + 2.0 * y + 3.0 * z; + }); + + projection.Apply(source, target); + + const auto target_values = + pcms::FlattenToRank1View(target.GetDOFHolderDataHost()); + const auto tgt_coords_h = pcms::test::CopyCoordinatesToHost( + pcms::MakeConstRank2View(target_mesh.coords(), 3), target_mesh.nverts(), + 3); + for (Omega_h::LO i = 0; i < target_mesh.nverts(); ++i) { + const double expected = 1.0 + tgt_coords_h(i, 0) + + 2.0 * tgt_coords_h(i, 1) + + 3.0 * tgt_coords_h(i, 2); + REQUIRE(target_values[i] == Catch::Approx(expected).margin(1e-8)); + } + REQUIRE(pcms::test::IntegrateP1Field(target_mesh, target) == + Catch::Approx(pcms::test::IntegrateP1Field(source_mesh, source)) + .margin(1e-8)); + } +} + +TEST_CASE("OmegaHConservativeProjection (3D tets) conserves the integral for a " + "P0 target", + "[transfer][mesh_intersection][3d]") +{ + Omega_h::Library lib; + + Omega_h::Mesh source_mesh = BuildUnitCube(lib, 2); + Omega_h::Mesh target_mesh = BuildUnitCube(lib, 1); + + auto source_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP0Space(target_mesh); + + auto source = source_space->CreateFunction(); + auto target = target_space->CreateFunction(); + + pcms::OmegaHConservativeProjection projection(*source_space, *target_space); + + pcms::test::SetField( + source, KOKKOS_LAMBDA(pcms::Real x, pcms::Real y, pcms::Real z) { + return 1.0 + x + 2.0 * y + 3.0 * z; + }); + + projection.Apply(source, target); + + REQUIRE(pcms::test::IntegrateP0Field(target_mesh, target) == + Catch::Approx(pcms::test::IntegrateP1Field(source_mesh, source)) + .margin(1e-8)); +} + +TEST_CASE("Copy transfer (3D tets) reproduces the source field", + "[transfer][copy][3d]") +{ + Omega_h::Library lib; + Omega_h::Mesh mesh = BuildUnitCube(lib, 2); + auto space = pcms::test::MakeP1Space(mesh); + + auto source = space->CreateFunction(); + auto target = space->CreateFunction(); + + pcms::test::SetField( + source, KOKKOS_LAMBDA(pcms::Real x, pcms::Real y, pcms::Real z) { + return 1.0 + x + 2.0 * y + 3.0 * z; + }); + + pcms::Copy copy(*space, *space); + copy.Apply(source, target); + + const auto sv = pcms::FlattenToRank1View(source.GetDOFHolderDataHost()); + const auto tv = pcms::FlattenToRank1View(target.GetDOFHolderDataHost()); + REQUIRE(tv.size() == sv.size()); + for (std::size_t i = 0; i < tv.size(); ++i) { + REQUIRE(tv[i] == Catch::Approx(sv[i])); + } +} + +// Point interpolation evaluates the source field at each target DOF site. A +// linear field is reproduced exactly at the target vertices in 3D. +TEST_CASE("Interpolation transfer (3D tets) reproduces a linear field", + "[transfer][interpolation][3d]") +{ + Omega_h::Library lib; + Omega_h::Mesh source_mesh = BuildUnitCube(lib, 2); + Omega_h::Mesh target_mesh = BuildUnitCube(lib, 3); + + auto source_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh); + + auto source = source_space->CreateFunction(); + auto target = target_space->CreateFunction(); + + pcms::test::SetField( + source, KOKKOS_LAMBDA(pcms::Real x, pcms::Real y, pcms::Real z) { + return 1.0 + x + 2.0 * y + 3.0 * z; + }); + + pcms::Interpolator interp(*source_space, *target_space); + interp.Apply(source, target); + + const auto target_values = + pcms::FlattenToRank1View(target.GetDOFHolderDataHost()); + const auto tgt_coords_h = pcms::test::CopyCoordinatesToHost( + pcms::MakeConstRank2View(target_mesh.coords(), 3), target_mesh.nverts(), 3); + for (Omega_h::LO i = 0; i < target_mesh.nverts(); ++i) { + const double expected = 1.0 + tgt_coords_h(i, 0) + + 2.0 * tgt_coords_h(i, 1) + 3.0 * tgt_coords_h(i, 2); + REQUIRE(target_values[i] == Catch::Approx(expected).margin(1e-8)); + } +} + +// The Monte-Carlo/control-variate projection uses the source field interpolated +// onto the target space as a control variate, so a field already representable +// in the target P1 space (an affine function) is reproduced exactly and its +// integral conserved, even with very few stochastic samples. +TEST_CASE("OmegaHControlVariateProjection (3D tets) is exact for target-space " + "fields", + "[transfer][monte_carlo][3d]") +{ + Omega_h::Library lib; + + Omega_h::Mesh source_mesh = BuildUnitCube(lib, 1); + Omega_h::Mesh target_mesh = BuildUnitCube(lib, 2); + + auto source_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh); + + auto source = source_space->CreateFunction(); + auto target = target_space->CreateFunction(); + + pcms::test::SetField( + source, KOKKOS_LAMBDA(pcms::Real x, pcms::Real y, pcms::Real z) { + return 1.0 + x + 2.0 * y + 3.0 * z; + }); + + pcms::OmegaHControlVariateProjection projection( + *source_space, *target_space, /*samples_per_element=*/8, + pcms::MonteCarloSampling::UniformRandom, /*seed=*/12345); + projection.Apply(source, target); + + const auto target_values = + pcms::FlattenToRank1View(target.GetDOFHolderDataHost()); + const auto tgt_coords_h = pcms::test::CopyCoordinatesToHost( + pcms::MakeConstRank2View(target_mesh.coords(), 3), target_mesh.nverts(), 3); + for (Omega_h::LO i = 0; i < target_mesh.nverts(); ++i) { + const double expected = 1.0 + tgt_coords_h(i, 0) + + 2.0 * tgt_coords_h(i, 1) + 3.0 * tgt_coords_h(i, 2); + REQUIRE(target_values[i] == Catch::Approx(expected).margin(1e-8)); + } + REQUIRE(pcms::test::IntegrateP1Field(target_mesh, target) == + Catch::Approx(pcms::test::IntegrateP1Field(source_mesh, source)) + .margin(1e-8)); +} diff --git a/test/test_omega_h_form_integrator_utils.cpp b/test/test_omega_h_form_integrator_utils.cpp index 91cedac7..785fd19a 100644 --- a/test/test_omega_h_form_integrator_utils.cpp +++ b/test/test_omega_h_form_integrator_utils.cpp @@ -2,6 +2,10 @@ #include #include +#include +#include +#include + namespace { @@ -78,3 +82,117 @@ TEST_CASE( REQUIRE(n == 3); REQUIRE(poly.nverts == 3); } + +// --------------------------------------------------------------------------- +// ForEachPolytopeFaceTriangle: the 3D face-walk that ForEachIntersectionSubtet +// relies on to star-decompose a clipped intersection polyhedron. +// --------------------------------------------------------------------------- + +namespace +{ + +double TetVolume(const r3d::Vector<3>& a, const r3d::Vector<3>& b, + const r3d::Vector<3>& c, const r3d::Vector<3>& d) +{ + const double bx = b[0] - a[0], by = b[1] - a[1], bz = b[2] - a[2]; + const double cx = c[0] - a[0], cy = c[1] - a[1], cz = c[2] - a[2]; + const double dx = d[0] - a[0], dy = d[1] - a[1], dz = d[2] - a[2]; + const double triple = bx * (cy * dz - cz * dy) - by * (cx * dz - cz * dx) + + bz * (cx * dy - cy * dx); + return std::abs(triple) / 6.0; +} + +// Star-decompose `poly` from its centroid exactly as ForEachIntersectionSubtet +// does, summing the sub-tet volumes and counting the emitted face triangles. +std::pair DecomposeFromCentroid(const r3d::Polytope<3>& poly) +{ + r3d::Vector<3> apex; + apex[0] = apex[1] = apex[2] = 0.0; + for (int v = 0; v < poly.nverts; ++v) { + apex[0] += poly.verts[v].pos[0]; + apex[1] += poly.verts[v].pos[1]; + apex[2] += poly.verts[v].pos[2]; + } + apex[0] /= poly.nverts; + apex[1] /= poly.nverts; + apex[2] /= poly.nverts; + + double total_vol = 0.0; + int ntri = 0; + pcms::detail::ForEachPolytopeFaceTriangle(poly, [&](const r3d::Vector<3>& a, + const r3d::Vector<3>& b, + const r3d::Vector<3>& c) { + total_vol += TetVolume(apex, a, b, c); + ++ntri; + }); + return {total_vol, ntri}; +} + +r3d::Few, 4> MakeTetVerts( + std::initializer_list> pts) +{ + r3d::Few, 4> verts; + int i = 0; + for (const auto& p : pts) { + verts[i][0] = p[0]; + verts[i][1] = p[1]; + verts[i][2] = p[2]; + ++i; + } + return verts; +} + +} // namespace + +TEST_CASE("ForEachPolytopeFaceTriangle: a tetrahedron yields 4 face triangles " + "that tile its volume", + "[form_integrator_utils]") +{ + // A tetrahedron has 4 triangular faces, so the walk must emit exactly 4 + // triangles (= 2*(nverts-2) for nverts==4), and the centroid star-tiling must + // reproduce the tet's volume (1/6 for the reference tet). + r3d::Polytope<3> poly; + r3d::init(poly, MakeTetVerts({{{0.0, 0.0, 0.0}}, + {{1.0, 0.0, 0.0}}, + {{0.0, 1.0, 0.0}}, + {{0.0, 0.0, 1.0}}})); + + const auto [vol, ntri] = DecomposeFromCentroid(poly); + REQUIRE(ntri == 4); + REQUIRE(ntri == 2 * (poly.nverts - 2)); + REQUIRE(vol == Catch::Approx(1.0 / 6.0).epsilon(1e-12)); + REQUIRE(vol == Catch::Approx(std::abs(r3d::measure(poly))).epsilon(1e-12)); +} + +TEST_CASE( + "ForEachPolytopeFaceTriangle: a clipped polyhedron with quad faces is " + "tiled exactly", + "[form_integrator_utils]") +{ + // Clip the reference tet with the plane x <= 0.5. This truncates the corner + // at (1,0,0), producing a polyhedron with a quadrilateral face, so at least + // one face must fan into more than one triangle. The removed piece is a tet + // similar to the original at scale 0.5 (volume (1/6)*0.5^3 = 1/48), leaving + // 1/6 - 1/48 = 7/48. + r3d::Polytope<3> poly; + r3d::init(poly, MakeTetVerts({{{0.0, 0.0, 0.0}}, + {{1.0, 0.0, 0.0}}, + {{0.0, 1.0, 0.0}}, + {{0.0, 0.0, 1.0}}})); + r3d::Few, 1> planes; + planes[0].n[0] = -1.0; // keep -x + 0.5 >= 0, i.e. x <= 0.5 + planes[0].n[1] = 0.0; + planes[0].n[2] = 0.0; + planes[0].d = 0.5; + r3d::clip(poly, planes); + + REQUIRE(poly.nverts > 4); // truncation added vertices / a quad face + const double measure = std::abs(r3d::measure(poly)); + REQUIRE(measure == Catch::Approx(7.0 / 48.0).epsilon(1e-12)); + + const auto [vol, ntri] = DecomposeFromCentroid(poly); + // Every face walked exactly once <=> the genus-0 fan invariant holds. + REQUIRE(ntri == 2 * (poly.nverts - 2)); + REQUIRE(ntri > 4); // a quad face fans into 2+ triangles + REQUIRE(vol == Catch::Approx(measure).epsilon(1e-12)); +} diff --git a/test/test_omega_h_intersection_rhs_integrator.cpp b/test/test_omega_h_intersection_rhs_integrator.cpp index 18ea11bb..f5d667d9 100644 --- a/test/test_omega_h_intersection_rhs_integrator.cpp +++ b/test/test_omega_h_intersection_rhs_integrator.cpp @@ -21,37 +21,6 @@ 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 @@ -100,17 +69,7 @@ void CheckAssembledLoadMatches( const pcms::Field& source_field, const std::unordered_map& expected) { - const auto pts = integrator.GetIntegrationPoints(); - const std::size_t npts = pts.GetValues().extent(0); - - auto evaluator = source_space->CreatePointEvaluator( - pcms::EvaluationRequest::FromCoordinates(pts)); - - Kokkos::View sampled("sampled", npts, - 1); - evaluator->Evaluate(source_field, pcms::MakeRank2View(sampled)); - - integrator.Assemble(pcms::MakeConstRank2View(sampled)); + pcms::test::EvaluateAndAssemble(integrator, source_space, source_field); Vec vec = integrator.GetVector(); const PetscScalar* vec_array = nullptr; @@ -135,26 +94,19 @@ TEST_CASE("OmegaHIntersectionRHSIntegrator: integration points lie inside " "[rhs_integrator]") { Omega_h::Library lib; - auto source_mesh = BuildUnitSquare(lib, 1); - auto target_mesh = BuildUnitSquare(lib, 0); + auto source_mesh = pcms::test::BuildUnitSquare(lib, 1); + auto target_mesh = pcms::test::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_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh); auto integrator = pcms::BuildOmegaHConservativeRHSIntegrator(*source_space, *target_space); const auto raw_coords = integrator->GetIntegrationPoints().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); + auto raw_coords_host = pcms::test::CopyCoordinatesToHost( + raw_coords, static_cast(raw_coords.extent(0)), + static_cast(raw_coords.extent(1))); const std::size_t n = raw_coords_host.extent(0); REQUIRE(n > 0); @@ -172,32 +124,18 @@ TEST_CASE("OmegaHIntersectionRHSIntegrator: zero source field gives zero " "[rhs_integrator]") { Omega_h::Library lib; - auto source_mesh = BuildUnitSquare(lib, 1); - auto target_mesh = BuildUnitSquare(lib, 0); + auto source_mesh = pcms::test::BuildUnitSquare(lib, 1); + auto target_mesh = pcms::test::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_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh); auto source_field = source_space->CreateFunction(); // 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.GetValues().extent(0); - - auto evaluator = source_space->CreatePointEvaluator( - pcms::EvaluationRequest::FromCoordinates(pts)); - - Kokkos::View sampled("sampled", npts, - 1); - evaluator->Evaluate(source_field, pcms::MakeRank2View(sampled)); - - integrator->Assemble(pcms::MakeConstRank2View(sampled)); + pcms::test::EvaluateAndAssemble(*integrator, source_space, source_field); Vec vec = integrator->GetVector(); PetscReal norm = 0.0; @@ -214,15 +152,11 @@ TEST_CASE("OmegaHIntersectionRHSIntegrator: constant field matches " // 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_mesh = pcms::test::BuildUnitSquare(lib, 1); + auto target_mesh = pcms::test::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_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh); const double c = 2.0; auto source_field = source_space->CreateFunction(); @@ -251,15 +185,11 @@ TEST_CASE( // 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_mesh = pcms::test::BuildUnitSquare(lib, 1); + auto target_mesh = pcms::test::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_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh); auto source_field = source_space->CreateFunction(); pcms::test::SetField( @@ -268,11 +198,11 @@ TEST_CASE( const auto target_layout = std::dynamic_pointer_cast( target_space->GetLayout()); - const auto tgt_coords_h = - Omega_h::HostRead(target_mesh.coords()); + const auto tgt_coords_h = pcms::test::CopyCoordinatesToHost( + pcms::MakeConstRank2View(target_mesh.coords(), 2), target_mesh.nverts(), 2); 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]; + g[i] = tgt_coords_h(i, 0) + tgt_coords_h(i, 1); } const auto expected = ExpectedLoadByGid(target_mesh, *target_layout, g); @@ -285,26 +215,22 @@ TEST_CASE("OmegaHIntersectionRHSIntegrator: rejects invalid layouts", "[rhs_integrator]") { Omega_h::Library lib; - auto source_mesh = BuildUnitSquare(lib, 1); - auto target_mesh = BuildUnitSquare(lib, 0); + auto source_mesh = pcms::test::BuildUnitSquare(lib, 1); + auto target_mesh = pcms::test::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); + auto target_space = pcms::test::MakeP1Space(target_mesh); 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 source_space = pcms::test::MakeP1Space(source_mesh); auto target_space = pcms::LagrangeFunctionSpace::FromMesh( target_mesh, 1, 2, pcms::CoordinateSystem::Cartesian, "global", pcms::LagrangeFunctionSpace::Backend::OmegaH); @@ -317,18 +243,14 @@ TEST_CASE("OmegaHIntersectionRHSIntegrator: rejects invalid layouts", 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); + auto target_space = pcms::test::MakeP1Space(target_mesh); 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 source_space = pcms::test::MakeP1Space(source_mesh); auto target_space = pcms::LagrangeFunctionSpace::FromMesh( target_mesh, 1, 1, pcms::CoordinateSystem::Cylindrical, "global", pcms::LagrangeFunctionSpace::Backend::OmegaH); diff --git a/test/test_omega_h_lagrange_field.cpp b/test/test_omega_h_lagrange_field.cpp index c17d45b1..0dd5ffc7 100644 --- a/test/test_omega_h_lagrange_field.cpp +++ b/test/test_omega_h_lagrange_field.cpp @@ -155,10 +155,10 @@ TEST_CASE("OmegaHLagrangeField order-1: linear function evaluation") pcms::test::SetField( field.GetData(), *factory->GetLayout(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); pcms::test::CheckEvaluation( factory, field, pcms::test::StandardEvalCoords2D(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); } // The same linear evaluation test run on the MeshFields-backed order-1 field @@ -173,10 +173,10 @@ TEST_CASE("MeshFieldsAdapter order-1: linear function evaluation (shared util)") pcms::test::SetField( field.GetData(), *factory->GetLayout(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); pcms::test::CheckEvaluation( factory, field, pcms::test::StandardEvalCoords2D(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); } TEST_CASE("OmegaHLagrangeField order-1: out-of-bounds FILL mode") @@ -189,10 +189,10 @@ TEST_CASE("OmegaHLagrangeField order-1: out-of-bounds FILL mode") pcms::test::SetField( field.GetData(), *factory->GetLayout(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); Real fill_value = -999.0; - std::vector outside{-0.5, 0.5, 1.5, 0.5, 0.5, -0.5, 0.5, 1.5}; + auto outside = pcms::test::StandardOutsideCoords2D(); pcms::test::CheckFillMode(factory, field, fill_value, outside); } @@ -206,7 +206,7 @@ TEST_CASE("OmegaHLagrangeField order-1: serialize / deserialize round-trip") pcms::test::SetField( field.GetData(), *factory->GetLayout(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); pcms::test::CheckSerializeDeserialize(*factory->GetLayout(), field.GetData()); } @@ -270,24 +270,9 @@ TEST_CASE("OmegaHLagrangeField order-0: constant field evaluation") 1); field.GetData().SetDOFHolderDataHost(view); - auto pts = pcms::test::StandardEvalCoords2D(); - int n = static_cast(pts.size()) / 2; - auto device_coords = - pcms::test::CreateDeviceCoordinateView(pts, factory->GetCoordinateSystem()); - auto evaluator = factory->CreatePointEvaluator( - pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); - - Kokkos::View eval_device("eval", n); - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - auto out = pcms::Rank2View( - eval_device.data(), n, 1); - evaluator->Evaluate(field, out); - auto eval_host = - Kokkos::create_mirror_view_and_copy(pcms::HostMemorySpace(), eval_device); - - for (int i = 0; i < n; ++i) - REQUIRE(eval_host(i) == Catch::Approx(kValue)); + pcms::test::CheckEvaluation( + factory, field, pcms::test::StandardEvalCoords2D(), + OMEGA_H_LAMBDA(Real, Real) { return kValue; }); } TEST_CASE("OmegaHLagrangeField order-0: out-of-bounds FILL mode") @@ -357,7 +342,8 @@ TEST_CASE("OmegaHLagrangeField: field valid after layout destruction") } // factory goes out of scope; field keeps layout alive pcms::test::SetField( - *field, OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + *field, + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); // Just verify data was set correctly (no evaluator needed for this lifetime // test) auto data = field->GetDOFHolderDataHost(); diff --git a/test/test_omega_h_mass_integrator.cpp b/test/test_omega_h_mass_integrator.cpp index c9013f9d..c49a2e7d 100644 --- a/test/test_omega_h_mass_integrator.cpp +++ b/test/test_omega_h_mass_integrator.cpp @@ -8,6 +8,7 @@ #include #include #include +#include "field_test_utils.h" #include #include #include @@ -23,33 +24,6 @@ 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) { @@ -98,11 +72,9 @@ TEST_CASE( "[mass_integrator]") { Omega_h::Library lib; - auto mesh = BuildUnitSquare(lib); + auto mesh = pcms::test::BuildUnitSquare(lib, 0); - auto space = pcms::LagrangeFunctionSpace::FromMesh( - mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", - pcms::LagrangeFunctionSpace::Backend::OmegaH); + auto space = pcms::test::MakeP1Space(mesh); const auto ref = BuildReferenceMassMap(mesh, *space); @@ -128,11 +100,9 @@ TEST_CASE("OmegaHMassIntegrator: row sums match lumped mass", // 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 mesh = pcms::test::BuildUnitSquare(lib, 0); - auto space = pcms::LagrangeFunctionSpace::FromMesh( - mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", - pcms::LagrangeFunctionSpace::Backend::OmegaH); + auto space = pcms::test::MakeP1Space(mesh); const auto layout = std::dynamic_pointer_cast( @@ -169,7 +139,7 @@ TEST_CASE("OmegaHMassIntegrator: row sums match lumped mass", TEST_CASE("OmegaHMassIntegrator: rejects invalid layouts", "[mass_integrator]") { Omega_h::Library lib; - auto mesh = BuildUnitSquare(lib); + auto mesh = pcms::test::BuildUnitSquare(lib, 0); SECTION("multi-component space throws") { diff --git a/test/test_omega_h_mc_rhs_integrator.cpp b/test/test_omega_h_mc_rhs_integrator.cpp index 43172fce..4cf8a651 100644 --- a/test/test_omega_h_mc_rhs_integrator.cpp +++ b/test/test_omega_h_mc_rhs_integrator.cpp @@ -15,69 +15,6 @@ #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; -} - -std::shared_ptr 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 std::shared_ptr& source_space, - const pcms::Field& source_field) -{ - const auto& pts = integrator.GetIntegrationPoints(); - const std::size_t npts = pts.GetValues().extent(0); - auto evaluator = source_space->CreatePointEvaluator( - pcms::EvaluationRequest::FromCoordinates(pts)); - Kokkos::View sampled("sampled", npts, - 1); - evaluator->Evaluate(source_field, pcms::MakeRank2View(sampled)); - integrator.Assemble(pcms::MakeConstRank2View(sampled)); -} - -} // namespace - // --------------------------------------------------------------------------- // Monte Carlo RHS integrator // --------------------------------------------------------------------------- @@ -86,20 +23,17 @@ 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); + auto target_mesh = pcms::test::BuildUnitSquare(lib, 0); + auto target_space = pcms::test::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().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); + auto coords_h = pcms::test::CopyCoordinatesToHost( + raw_coords, static_cast(raw_coords.extent(0)), + static_cast(raw_coords.extent(1))); REQUIRE( coords_h.extent(0) == @@ -120,10 +54,10 @@ TEST_CASE("OmegaHMonteCarloRHSIntegrator: constant field integrates exactly", // 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_mesh = pcms::test::BuildUnitSquare(lib, 1); + auto target_mesh = pcms::test::BuildUnitSquare(lib, 0); + auto source_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh); auto source_field = source_space->CreateFunction(); pcms::test::SetField( @@ -131,7 +65,7 @@ TEST_CASE("OmegaHMonteCarloRHSIntegrator: constant field integrates exactly", for (const auto sampling : {pcms::MonteCarloSampling::UniformRandom}) { pcms::OmegaHMonteCarloRHSIntegrator integrator(*target_space, 8, sampling); - EvaluateAndAssemble(integrator, source_space, source_field); + pcms::test::EvaluateAndAssemble(integrator, source_space, source_field); PetscScalar sum = 0.0; VecSum(integrator.GetVector(), &sum); @@ -152,10 +86,10 @@ TEST_CASE("OmegaHControlVariateProjection: exact for fields in the target " // 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_mesh = pcms::test::BuildUnitSquare(lib, 1); + auto target_mesh = pcms::test::BuildUnitSquare(lib, 0); + auto source_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh); auto source_field = source_space->CreateFunction(); auto target_field = target_space->CreateFunction(); @@ -170,11 +104,11 @@ TEST_CASE("OmegaHControlVariateProjection: exact for fields in the target " const auto values = pcms::FlattenToRank1View(target_field.GetDOFHolderDataHost()); - const auto coords_h = Omega_h::HostRead(target_mesh.coords()); + const auto coords_h = pcms::test::CopyCoordinatesToHost( + pcms::MakeConstRank2View(target_mesh.coords(), 2), target_mesh.nverts(), 2); 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; + const double expected = 3.0 * coords_h(i, 0) - coords_h(i, 1) + 0.25; CAPTURE(i, expected, values[i]); CHECK(values[i] == Catch::Approx(expected).margin(1e-9)); } @@ -188,10 +122,10 @@ TEST_CASE("OmegaHControlVariateProjection: reduces error vs plain Monte Carlo", // 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_mesh = pcms::test::BuildUnitSquare(lib, 1); + auto target_mesh = pcms::test::BuildUnitSquare(lib, 0); + auto source_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh); auto source_field = source_space->CreateFunction(); pcms::test::SetField( diff --git a/test/test_point_evaluator.cpp b/test/test_point_evaluator.cpp index dcc57087..e9e235a9 100644 --- a/test/test_point_evaluator.cpp +++ b/test/test_point_evaluator.cpp @@ -34,7 +34,7 @@ TEST_CASE("PointEvaluator: OmegaH order-1 linear evaluation") auto field_data = factory->CreateFunction(); pcms::test::SetField( field_data.GetData(), *factory->GetLayout(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); auto pts = pcms::test::StandardEvalCoords2D(); int n = static_cast(pts.size()) / 2; @@ -44,7 +44,7 @@ TEST_CASE("PointEvaluator: OmegaH order-1 linear evaluation") pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); pcms::test::CheckEvaluation( *evaluator, field_data, pts, - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); } // ============================================================================ @@ -67,45 +67,24 @@ TEST_CASE("PointEvaluator: same evaluator reused for two FieldData objects") // field_a: linear_f; field_b: constant 42 pcms::test::SetField( field_a.GetData(), *factory->GetLayout(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); pcms::test::SetField( field_b.GetData(), *factory->GetLayout(), OMEGA_H_LAMBDA(Real, Real) { return Real(42); }); auto pts = pcms::test::StandardEvalCoords2D(); - int n = static_cast(pts.size()) / 2; auto device_coords = pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian); - // Create the PointEvaluator once + // Create the PointEvaluator once and reuse it for both fields. auto evaluator = factory->CreatePointEvaluator( pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); - Kokkos::View out_a_device("out_a", n); - Kokkos::View out_b_device("out_b", n); - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - auto view_a = pcms::Rank2View( - out_a_device.data(), n, 1); - auto view_b = pcms::Rank2View( - out_b_device.data(), n, 1); - - // Evaluate field_a then field_b with the same evaluator - evaluator->Evaluate(field_a, view_a); - evaluator->Evaluate(field_b, view_b); - - auto out_a_host = - Kokkos::create_mirror_view_and_copy(pcms::HostMemorySpace(), out_a_device); - auto out_b_host = - Kokkos::create_mirror_view_and_copy(pcms::HostMemorySpace(), out_b_device); - - for (int i = 0; i < n; ++i) { - Real x = pts[2 * static_cast(i)], - y = pts[2 * static_cast(i) + 1]; - REQUIRE(out_a_host(i) == - Catch::Approx(pcms::test::linear_f(x, y)).margin(1e-10)); - REQUIRE(out_b_host(i) == Catch::Approx(42.0).margin(1e-10)); - } + pcms::test::CheckEvaluation( + *evaluator, field_a, pts, + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); + pcms::test::CheckEvaluation( + *evaluator, field_b, pts, OMEGA_H_LAMBDA(Real, Real) { return Real(42); }); } // ============================================================================ @@ -125,11 +104,10 @@ TEST_CASE("PointEvaluator: OmegaH order-1 out-of-bounds fill") auto field_data = factory->CreateFunction(); pcms::test::SetField( field_data.GetData(), *factory->GetLayout(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); // Points clearly outside [0,1]^2 - const std::vector outside_pts = {-0.5, 0.5, 1.5, 0.5, - 0.5, -0.5, 0.5, 1.5}; + const auto outside_pts = pcms::test::StandardOutsideCoords2D(); auto device_coords = pcms::test::CreateDeviceCoordinateView( outside_pts, CoordinateSystem::Cartesian); pcms::OutOfBoundsPolicy policy{pcms::OutOfBoundsMode::FILL, -999.0}; @@ -158,7 +136,7 @@ TEST_CASE("PointEvaluator: UniformGrid order-1 linear evaluation") auto field_data = factory->CreateFunction(); pcms::test::SetField( field_data.GetData(), *factory->GetLayout(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); auto pts = pcms::test::StandardEvalCoords2D(); int n = static_cast(pts.size()) / 2; auto device_coords = @@ -167,7 +145,8 @@ TEST_CASE("PointEvaluator: UniformGrid order-1 linear evaluation") pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); pcms::test::CheckEvaluation( *evaluator, field_data, pts, - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }, 1e-8); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }, + 1e-8); } TEST_CASE("PointEvaluator: SplineFunctionSpace uniform-grid evaluation") @@ -184,7 +163,7 @@ TEST_CASE("PointEvaluator: SplineFunctionSpace uniform-grid evaluation") auto field_data = factory->CreateFunction(); pcms::test::SetField( field_data.GetData(), *factory->GetLayout(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); auto pts = pcms::test::StandardEvalCoords2D(); int n = static_cast(pts.size()) / 2; auto device_coords = @@ -193,7 +172,8 @@ TEST_CASE("PointEvaluator: SplineFunctionSpace uniform-grid evaluation") pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); pcms::test::CheckEvaluation( *evaluator, field_data, pts, - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }, 1e-8); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }, + 1e-8); } // ============================================================================ @@ -308,7 +288,7 @@ TEST_CASE("PointEvaluator: MeshFields order-1 linear evaluation") auto field_data = factory->CreateFunction(); pcms::test::SetField( field_data.GetData(), *factory->GetLayout(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); auto pts = pcms::test::StandardEvalCoords2D(); int n = static_cast(pts.size()) / 2; @@ -318,7 +298,7 @@ TEST_CASE("PointEvaluator: MeshFields order-1 linear evaluation") pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); pcms::test::CheckEvaluation( *evaluator, field_data, pts, - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); } TEST_CASE("PointEvaluator: MeshFields out-of-bounds fill") @@ -334,10 +314,9 @@ TEST_CASE("PointEvaluator: MeshFields out-of-bounds fill") auto field_data = factory->CreateFunction(); pcms::test::SetField( field_data.GetData(), *factory->GetLayout(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); - const std::vector outside_pts = {-0.5, 0.5, 1.5, 0.5, - 0.5, -0.5, 0.5, 1.5}; + const auto outside_pts = pcms::test::StandardOutsideCoords2D(); auto device_coords = pcms::test::CreateDeviceCoordinateView( outside_pts, CoordinateSystem::Cartesian); pcms::OutOfBoundsPolicy policy{pcms::OutOfBoundsMode::FILL, -999.0}; @@ -362,43 +341,24 @@ TEST_CASE( auto field_b = factory->CreateFunction(); pcms::test::SetField( field_a.GetData(), *factory->GetLayout(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); pcms::test::SetField( field_b.GetData(), *factory->GetLayout(), OMEGA_H_LAMBDA(Real, Real) { return Real(42); }); auto pts = pcms::test::StandardEvalCoords2D(); - int n = static_cast(pts.size()) / 2; auto device_coords = pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian); + // Create the PointEvaluator once and reuse it for both fields. auto evaluator = factory->CreatePointEvaluator( pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); - Kokkos::View out_a_device("out_a", n); - Kokkos::View out_b_device("out_b", n); - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - auto view_a = pcms::Rank2View( - out_a_device.data(), n, 1); - auto view_b = pcms::Rank2View( - out_b_device.data(), n, 1); - - evaluator->Evaluate(field_a, view_a); - evaluator->Evaluate(field_b, view_b); - - auto out_a_host = - Kokkos::create_mirror_view_and_copy(pcms::HostMemorySpace(), out_a_device); - auto out_b_host = - Kokkos::create_mirror_view_and_copy(pcms::HostMemorySpace(), out_b_device); - - for (int i = 0; i < n; ++i) { - Real x = pts[2 * static_cast(i)], - y = pts[2 * static_cast(i) + 1]; - REQUIRE(out_a_host(i) == - Catch::Approx(pcms::test::linear_f(x, y)).margin(1e-10)); - REQUIRE(out_b_host(i) == Catch::Approx(42.0).margin(1e-10)); - } + pcms::test::CheckEvaluation( + *evaluator, field_a, pts, + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); + pcms::test::CheckEvaluation( + *evaluator, field_b, pts, OMEGA_H_LAMBDA(Real, Real) { return Real(42); }); } TEST_CASE("LagrangeFunctionSpace: MeshFields rejects multi-component fields") diff --git a/test/test_polynomial_reconstruction_mls_evaluation.cpp b/test/test_polynomial_reconstruction_mls_evaluation.cpp index 888a83b1..415d4968 100644 --- a/test/test_polynomial_reconstruction_mls_evaluation.cpp +++ b/test/test_polynomial_reconstruction_mls_evaluation.cpp @@ -77,12 +77,6 @@ pcms::MLSOptions DefaultTestOptions3D() return opts; } -// Interior query points for a unit box — same as StandardEvalCoords2D. -std::vector QueryPoints() -{ - return pcms::test::StandardEvalCoords2D(); -} - pcms::MLSOptions SweepTestOptions(unsigned degree, pcms::RadialBasisFunction basis) { @@ -122,7 +116,7 @@ void CheckPolynomialReproduction(unsigned degree, auto field = fs->CreateFunction(); pcms::test::SetField(field.GetData(), *fs->GetLayout(), func); - auto pts = QueryPoints(); + auto pts = pcms::test::StandardEvalCoords2D(); auto device_coords = pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian); auto evaluator = fs->CreatePointEvaluator( @@ -149,7 +143,9 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: reproduces " CheckPolynomialReproduction( 0, basis, OMEGA_H_LAMBDA(Real, Real) { return 3.14; }, 5e-3); CheckPolynomialReproduction( - 1, basis, OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }, 5e-3); + 1, basis, + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }, + 5e-3); CheckPolynomialReproduction( 2, basis, OMEGA_H_LAMBDA(Real x, Real y) { return x * x + x * y + 2.0 * y * y; }, @@ -175,43 +171,26 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: same PointEvaluator " pcms::test::SetField( field_a.GetData(), *fs->GetLayout(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); const Real cval = 7.0; pcms::test::SetField( field_b.GetData(), *fs->GetLayout(), OMEGA_H_LAMBDA(Real, Real) { return cval; }); - auto pts = QueryPoints(); - int n = static_cast(pts.size()) / 2; + auto pts = pcms::test::StandardEvalCoords2D(); auto device_coords = pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian); + // Create the PointEvaluator once and reuse it for both fields. auto evaluator = fs->CreatePointEvaluator( pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); - Kokkos::View out_a_device("out_a", n); - Kokkos::View out_b_device("out_b", n); - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - Rank2View view_a(out_a_device.data(), - n, 1); - Rank2View view_b(out_b_device.data(), - n, 1); - - evaluator->Evaluate(field_a, view_a); - evaluator->Evaluate(field_b, view_b); - - auto out_a_host = - Kokkos::create_mirror_view_and_copy(HostMemorySpace(), out_a_device); - auto out_b_host = - Kokkos::create_mirror_view_and_copy(HostMemorySpace(), out_b_device); - - for (int i = 0; i < n; ++i) { - Real x = pts[2 * static_cast(i)], - y = pts[2 * static_cast(i) + 1]; - REQUIRE(out_a_host(i) == - Catch::Approx(pcms::test::linear_f(x, y)).margin(5e-3)); - REQUIRE(out_b_host(i) == Catch::Approx(cval).margin(5e-3)); - } + pcms::test::CheckEvaluation( + *evaluator, field_a, pts, + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }, + 5e-3); + pcms::test::CheckEvaluation( + *evaluator, field_b, pts, OMEGA_H_LAMBDA(Real, Real) { return cval; }, + 5e-3); } // ============================================================================ @@ -228,7 +207,7 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: Evaluate throws for " coords_view, CoordinateSystem::Cartesian); auto field = fs->CreateFunction(); - auto pts = QueryPoints(); + auto pts = pcms::test::StandardEvalCoords2D(); int n = static_cast(pts.size()) / 2; auto device_coords = pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian); @@ -236,13 +215,8 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: Evaluate throws for " pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); // Two-component output — must throw - Kokkos::View out_device("out", - static_cast(n) * 2); - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - auto out_view = - Rank2View(out_device.data(), n, 2); - REQUIRE_THROWS(evaluator->Evaluate(field, out_view)); + Kokkos::View out_device("out", n, 2); + REQUIRE_THROWS(evaluator->Evaluate(field, pcms::MakeRank2View(out_device))); } // ============================================================================ @@ -263,24 +237,20 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: default MLSOptions — " field.GetData(), *fs->GetLayout(), OMEGA_H_LAMBDA(Real, Real) { return Real(1.0); }); - auto pts = QueryPoints(); + auto pts = pcms::test::StandardEvalCoords2D(); int n = static_cast(pts.size()) / 2; auto device_coords = pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian); auto evaluator = fs->CreatePointEvaluator( pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); - Kokkos::View out_device("out", n); - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - auto out_view = - Rank2View(out_device.data(), n, 1); + Kokkos::View out_device("out", n, 1); // Just verify it runs without error and returns finite values - REQUIRE_NOTHROW(evaluator->Evaluate(field, out_view)); + REQUIRE_NOTHROW(evaluator->Evaluate(field, pcms::MakeRank2View(out_device))); auto out_host = Kokkos::create_mirror_view_and_copy(HostMemorySpace(), out_device); for (int i = 0; i < n; ++i) - REQUIRE(std::isfinite(out_host(i))); + REQUIRE(std::isfinite(out_host(i, 0))); } TEST_CASE("PolynomialReconstructionFunctionSpace MLS: CreatePointEvaluator " @@ -293,7 +263,7 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: CreatePointEvaluator " auto fs = pcms::PolynomialReconstructionFunctionSpace::Create( coords_view, CoordinateSystem::Cylindrical, DefaultTestOptions()); - auto pts = QueryPoints(); + auto pts = pcms::test::StandardEvalCoords2D(); auto device_coords = pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian); REQUIRE_THROWS(fs->CreatePointEvaluator( @@ -310,7 +280,7 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: CreatePointEvaluator " auto fs = pcms::PolynomialReconstructionFunctionSpace::Create( coords_view, CoordinateSystem::Cylindrical, DefaultTestOptions()); - auto pts = QueryPoints(); + auto pts = pcms::test::StandardEvalCoords2D(); auto device_coords = pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cylindrical); REQUIRE_THROWS(fs->CreatePointEvaluator( @@ -342,16 +312,12 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: radius option is " pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian); auto evaluator = fs->CreatePointEvaluator( pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); - Kokkos::View out_device("out", 1); - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - Rank2View out_view(out_device.data(), - 1, 1); - evaluator->Evaluate(field, out_view); + Kokkos::View out_device("out", 1, 1); + evaluator->Evaluate(field, pcms::MakeRank2View(out_device)); auto out_host = Kokkos::create_mirror_view_and_copy(HostMemorySpace(), out_device); - REQUIRE(out_host(0) == Catch::Approx(1.0).margin(1e-8)); + REQUIRE(out_host(0, 0) == Catch::Approx(1.0).margin(1e-8)); } TEST_CASE("PolynomialReconstructionFunctionSpace MLS: 3D point clouds preserve " @@ -375,34 +341,24 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: 3D point clouds preserve " } auto field = fs->CreateFunction(); - std::vector dof_values(27); - for (int i = 0; i < 27; ++i) { - const Real x = src[3 * i + 0]; - const Real y = src[3 * i + 1]; - const Real z = src[3 * i + 2]; - dof_values[i] = x + 2.0 * y + 3.0 * z; - } - field.GetData().SetDOFHolderDataHost(Rank2View( - dof_values.data(), static_cast(dof_values.size()), 1)); + pcms::test::SetField( + field, + KOKKOS_LAMBDA(Real x, Real y, Real z) { return x + 2.0 * y + 3.0 * z; }); std::vector pts{0.5, 0.5, 0.25, 0.5, 0.5, 0.75}; auto device_coords = pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian, 3); auto evaluator = fs->CreatePointEvaluator( pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); - Kokkos::View out_device("out", 2); - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - Rank2View out_view(out_device.data(), - 2, 1); - evaluator->Evaluate(field, out_view); + Kokkos::View out_device("out", 2, 1); + evaluator->Evaluate(field, pcms::MakeRank2View(out_device)); auto out_host = Kokkos::create_mirror_view_and_copy(HostMemorySpace(), out_device); // This 3D case is a small, low-resolution support cloud with MLS weights // built from a finite-radius neighborhood rather than exact nodal lookup, so // it is checked with a slightly looser tolerance than the denser 2D tests. - REQUIRE(out_host(0) == Catch::Approx(2.25).margin(1e-2)); - REQUIRE(out_host(1) == Catch::Approx(3.75).margin(1e-2)); - REQUIRE(out_host(1) - out_host(0) == Catch::Approx(1.5).margin(1e-2)); + REQUIRE(out_host(0, 0) == Catch::Approx(2.25).margin(1e-2)); + REQUIRE(out_host(1, 0) == Catch::Approx(3.75).margin(1e-2)); + REQUIRE(out_host(1, 0) - out_host(0, 0) == Catch::Approx(1.5).margin(1e-2)); } diff --git a/test/test_spr_meshfields.cpp b/test/test_spr_meshfields.cpp index 42a6acf0..7ad4470a 100644 --- a/test/test_spr_meshfields.cpp +++ b/test/test_spr_meshfields.cpp @@ -3,6 +3,7 @@ #include #include #include +#include #include #include #include @@ -127,22 +128,7 @@ TEST_CASE("meshfields_spr_test") const auto& ntargets = mesh.nverts(); - Write source_coordinates( - dim * nfaces, 0, "stores coordinates of cell centroid of each tri element"); - - const auto& faces2nodes = mesh.ask_down(FACE, VERT).ab2b; - - Kokkos::parallel_for( - "calculate the centroid in each tri element", nfaces, - OMEGA_H_LAMBDA(const LO id) { - const auto current_el_verts = gather_verts<3>(faces2nodes, id); - const Omega_h::Few, 3> current_el_vert_coords = - gather_vectors<3, 2>(target_coordinates, current_el_verts); - auto centroid = average(current_el_vert_coords); - int index = 2 * id; - source_coordinates[index] = centroid[0]; - source_coordinates[index + 1] = centroid[1]; - }); + const auto source_coordinates = pcms::get_entity_centroids(mesh, FACE); pcms::Points source_points; diff --git a/test/test_svd_serial.cpp b/test/test_svd_serial.cpp index bdc23290..1774e29b 100644 --- a/test/test_svd_serial.cpp +++ b/test/test_svd_serial.cpp @@ -38,7 +38,7 @@ TEST_CASE("test_serial_svd") Kokkos::deep_copy(A_data, host_A_data); - Kokkos::View rhs_data("Device rhs data", row, column); + Kokkos::View rhs_data("Device rhs data", row); auto host_rhs_data = Kokkos::create_mirror_view(rhs_data); host_rhs_data(0) = 1.28571; @@ -54,7 +54,7 @@ TEST_CASE("test_serial_svd") { Kokkos::View result("result", row, row); - Kokkos::View transpose_expected("result", row, row); + Kokkos::View transpose_expected("transpose_expected", row, row); Kokkos::deep_copy(result, 0.0); Kokkos::deep_copy(transpose_expected, 0.0); team_policy tp(1, Kokkos::AUTO); diff --git a/test/test_uniform_grid_field.cpp b/test/test_uniform_grid_field.cpp index 632bd4db..699ec26a 100644 --- a/test/test_uniform_grid_field.cpp +++ b/test/test_uniform_grid_field.cpp @@ -72,7 +72,7 @@ void VerifyUniformGridFieldValues( int vertex_id = j * (grid.divisions[0] + 1) + i; pcms::Real x = ug_coords.GetValues()(vertex_id, 0); pcms::Real y = ug_coords.GetValues()(vertex_id, 1); - pcms::Real expected = x + 2.0 * y; + pcms::Real expected = pcms::test::linear_f(x, y); pcms::Real actual = ug_field_data[vertex_id]; REQUIRE(std::abs(expected - actual) <= 1e-10); } @@ -151,21 +151,16 @@ TEST_CASE("UniformGrid order-0 field creation and evaluation") auto evaluator = eval_factory.CreatePointEvaluator( pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); - Kokkos::View results_host("results_host", - 4); - Kokkos::View results_device( - "results_device", 4); - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - pcms::Rank2View out( - results_device.data(), 4, 1); - evaluator->Evaluate(field, out); - Kokkos::deep_copy(results_host, results_device); - - REQUIRE(results_host(0) == Catch::Approx(1.0)); - REQUIRE(results_host(1) == Catch::Approx(2.0)); - REQUIRE(results_host(2) == Catch::Approx(3.0)); - REQUIRE(results_host(3) == Catch::Approx(4.0)); + Kokkos::View results_device( + "results_device", 4, 1); + evaluator->Evaluate(field, pcms::MakeRank2View(results_device)); + auto results_host = Kokkos::create_mirror_view_and_copy( + pcms::HostMemorySpace(), results_device); + + REQUIRE(results_host(0, 0) == Catch::Approx(1.0)); + REQUIRE(results_host(1, 0) == Catch::Approx(2.0)); + REQUIRE(results_host(2, 0) == Catch::Approx(3.0)); + REQUIRE(results_host(3, 0) == Catch::Approx(4.0)); } TEST_CASE("UniformGrid field data operations", "[uniform_grid_field]") @@ -236,26 +231,21 @@ TEST_CASE("UniformGrid field evaluation - piecewise constant") auto evaluator = eval_factory.CreatePointEvaluator( pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); - Kokkos::View results_host("results_host", - 4); - Kokkos::View results_device( - "results_device", 4); - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - pcms::Rank2View out( - results_device.data(), 4, 1); - evaluator->Evaluate(field, out); - Kokkos::deep_copy(results_host, results_device); + Kokkos::View results_device( + "results_device", 4, 1); + evaluator->Evaluate(field, pcms::MakeRank2View(results_device)); + auto results_host = Kokkos::create_mirror_view_and_copy( + pcms::HostMemorySpace(), results_device); // Check results - interpolated from vertices // Cell 0 center (2.5, 2.5): avg of v0,v1,v3,v4 = (1.0+1.5+2.0+2.5)/4 = 1.75 // Cell 1 center (7.5, 2.5): avg of v1,v2,v4,v5 = (1.5+2.0+2.5+3.0)/4 = 2.25 // Cell 2 center (2.5, 7.5): avg of v3,v4,v6,v7 = (2.0+2.5+3.0+3.5)/4 = 2.75 // Cell 3 center (7.5, 7.5): avg of v4,v5,v7,v8 = (2.5+3.0+3.5+4.0)/4 = 3.25 - REQUIRE(std::abs(results_host(0) - 1.75) < 1e-10); - REQUIRE(std::abs(results_host(1) - 2.25) < 1e-10); - REQUIRE(std::abs(results_host(2) - 2.75) < 1e-10); - REQUIRE(std::abs(results_host(3) - 3.25) < 1e-10); + REQUIRE(std::abs(results_host(0, 0) - 1.75) < 1e-10); + REQUIRE(std::abs(results_host(1, 0) - 2.25) < 1e-10); + REQUIRE(std::abs(results_host(2, 0) - 2.75) < 1e-10); + REQUIRE(std::abs(results_host(3, 0) - 3.25) < 1e-10); } TEST_CASE("UniformGrid field serialization") @@ -328,8 +318,9 @@ TEST_CASE("Transfer from OmegaH field to UniformGrid field") mesh, 1, 1, pcms::CoordinateSystem::Cartesian); auto omega_h_field = omega_h_factory->CreateFunction(); pcms::test::SetField( - omega_h_field, - OMEGA_H_LAMBDA(pcms::Real x, pcms::Real y) { return x + 2.0 * y; }); + omega_h_field, OMEGA_H_LAMBDA(pcms::Real x, pcms::Real y) { + return pcms::test::linear_f(x, y); + }); pcms::UniformGrid<2> grid; grid.edge_length = {1.0, 1.0}; @@ -354,7 +345,7 @@ TEST_CASE("Transfer from OmegaH field to UniformGrid field") for (int i = 0; i < num_ug_nodes; ++i) { pcms::Real x = ug_coords_host(i, 0); pcms::Real y = ug_coords_host(i, 1); - pcms::Real expected = x + 2.0 * y; + pcms::Real expected = pcms::test::linear_f(x, y); REQUIRE(std::abs(transferred_data[i] - expected) < 1e-6); } } @@ -641,8 +632,9 @@ TEST_CASE("UniformGrid workflow") mesh, 1, 1, pcms::CoordinateSystem::Cartesian); auto omega_h_field = omega_h_factory->CreateFunction(); pcms::test::SetField( - omega_h_field, - OMEGA_H_LAMBDA(pcms::Real x, pcms::Real y) { return x + 2.0 * y; }); + omega_h_field, OMEGA_H_LAMBDA(pcms::Real x, pcms::Real y) { + return pcms::test::linear_f(x, y); + }); auto ug_factory = pcms::LagrangeFunctionSpace::FromUniformGrid( grid, 1, pcms::CoordinateSystem::Cartesian); @@ -658,23 +650,11 @@ TEST_CASE("UniformGrid workflow") auto ug_coords_host_view = pcms::test::CopyCoordinatesToHost(ug_coords_device_view, 25, 2); - auto ug_field_data_device = ug_field.GetDOFHolderData(); - Kokkos::View ug_field_data_device_view( - "", 25); - Kokkos::parallel_for( - "CopyFieldDataToView", 25, KOKKOS_LAMBDA(int i) { - ug_field_data_device_view(i) = ug_field_data_device(i, 0); - }); - auto ug_field_data_host_view = - Kokkos::View("", 25); - Kokkos::deep_copy(ug_field_data_host_view, ug_field_data_device_view); - - pcms::Rank2View ug_coords( - ug_coords_host_view.data(), 25, 2); pcms::CoordinateView ug_coords_view( - pcms::CoordinateSystem::Cartesian, ug_coords); - pcms::Rank1View ug_field_data( - ug_field_data_host_view.data(), 25); + pcms::CoordinateSystem::Cartesian, + pcms::MakeConstRank2View(ug_coords_host_view)); + const auto ug_field_data = + pcms::FlattenToRank1View(ug_field.GetDOFHolderDataHost()); VerifyUniformGridFieldValues(grid, ug_coords_view, ug_field_data); VerifyMaskFieldValues(grid, mask_field); diff --git a/test/test_xgc_reverse_classification.cpp b/test/test_xgc_reverse_classification.cpp index 4b87fcb9..5603b52b 100644 --- a/test/test_xgc_reverse_classification.cpp +++ b/test/test_xgc_reverse_classification.cpp @@ -39,7 +39,7 @@ TEST_CASE("reverse classification") } auto vec = rc.Serialize(); pcms::ReverseClassificationVertex rc_deserialized; - pcms::Rank1View av{vec.data(), vec.size()}; + auto av = pcms::make_array_view(vec); rc_deserialized.Deserialize(av); REQUIRE(rc_deserialized == rc); }