From 9d9387c57ba88a88f463982c31cf86331632c99e Mon Sep 17 00:00:00 2001 From: Andreas Adelmann Date: Sat, 15 Aug 2026 16:04:32 +0200 Subject: [PATCH 1/9] Fix LagrangeSpace device captures for CUDA - replace class captures with a lightweight device mirror - add device-side mesh, DOF, and boundary helpers - capture field views directly in Kokkos kernels - eliminate GH200 host-only destructor warnings --- src/FEM/LagrangeSpace.h | 171 +++++++++++++++++++++- src/FEM/LagrangeSpace.hpp | 293 ++++++++++++++++++++++++++++---------- 2 files changed, 387 insertions(+), 77 deletions(-) diff --git a/src/FEM/LagrangeSpace.h b/src/FEM/LagrangeSpace.h index b744cd49b6..50c231da23 100644 --- a/src/FEM/LagrangeSpace.h +++ b/src/FEM/LagrangeSpace.h @@ -287,27 +287,194 @@ namespace ippl { /////////////////////////////////////////////////////////////////////// /// Device struct for copies ////////////////////////////////////////// /////////////////////////////////////////////////////////////////////// + /** + * @brief Device-copyable snapshot of the Lagrange-space geometry and indexing rules. + * + * @details + * `LagrangeSpace` is a host-side owner. In particular, it contains fields and other + * objects whose lifetime management and destructors are not device-callable. Capturing + * the complete space (or `this`) in a `KOKKOS_CLASS_LAMBDA` therefore makes the kernel + * closure contain host-only state. CUDA compilers can then diagnose calls to host-only + * constructors or destructors while generating device code. + * + * `DeviceStruct` defines the architectural boundary between that host-side owner and + * device kernels. It contains only the immutable, device-copyable values needed to + * reproduce mesh indexing, element geometry, degree-of-freedom mappings, and reference + * element evaluations. A kernel obtains a snapshot with getDeviceMirror() and captures + * that snapshot by value in a `KOKKOS_LAMBDA`. Kokkos views such as `elementIndices` and + * field data views are captured separately; they deliberately are not owned by this + * structure. + * + * More specifically, the mirror performs the following tasks: + * + * - It snapshots the mesh vertex counts (`nr_m`), mesh spacing (`hr_m`), physical origin + * (`origin_m`), and reference element (`ref_element_m`) that device code needs for + * geometric calculations. The compile-time element/DOF counts are reused without + * adding runtime storage. + * - It reconstructs the deterministic operations that previously required access to the + * parent class: flattened/N-dimensional element-index conversion, element-vertex + * enumeration, physical vertex-coordinate construction, global-DOF lookup, boundary + * detection, and reference-element shape-function evaluation. + * - It gives a `KOKKOS_LAMBDA` a self-contained, read-only description of the space. The + * host creates the structure before launching a kernel; capture-by-value places it in + * the Kokkos kernel closure, which Kokkos makes available in the selected execution + * space. getDeviceMirror() itself performs no field allocation or field-data copy. + * - It intentionally excludes owning or host-oriented state such as `resultField`, field + * and layout objects, MPI decomposition objects, and the `elementIndices` allocation. + * A kernel captures only the Kokkos views and scalar values needed for that invocation, + * separately from this geometry/indexing snapshot. + * - It removes the need to capture `this`. Consequently, constructing and destroying the + * device closure never requires the host-only lifetime operations of `LagrangeSpace` or + * `Field`, which is the source of the GH200 CUDA diagnostic addressed by this design. + * + * "Lightweight" therefore means that this is neither a second owning finite-element space + * nor an automatically synchronized copy of one. It is a small, non-owning value snapshot + * containing only kernel-invariant geometry and algorithms. Changes to the host object are + * not reflected in an existing mirror; the host must create another mirror before the next + * kernel that observes the changed state. + * + * When a new device kernel needs additional `LagrangeSpace` functionality, add the + * smallest required device-copyable state and a `KOKKOS_FUNCTION` helper here instead of + * capturing the parent object. All members must remain safe to copy into a device closure + * and must not introduce host-only ownership or lifetime management. A new snapshot must + * be created after changing any mirrored mesh or reference-element state on the host. + * During a kernel invocation the snapshot is read-only and may be shared by all threads. + * + * Element and vertex indices use the same flattened ordering as `LagrangeSpace`: + * dimension zero varies fastest. An element's N-dimensional index identifies its lower + * mesh vertex. + */ struct DeviceStruct { // members we need to copy for the following functions: // works since numElementDOFs in LagrangeSpace is static constexpr static constexpr unsigned numElementDOFs = LagrangeSpace::numElementDOFs; - Vector nr_m; - ElementType ref_element_m; + static constexpr unsigned numElementVertices = LagrangeSpace::numElementVertices; + using indices_list_t = Vector; + using vertex_points_t = Vector; + + Vector nr_m; ///< Number of mesh vertices in each dimension. + Vector hr_m; ///< Uniform mesh spacing in each dimension. + Vector origin_m; ///< Physical coordinate of mesh vertex index zero. + ElementType ref_element_m; ///< Device-copyable reference-element description. // these are the functions needed for interpolation to the space KOKKOS_FUNCTION indices_t getMeshVertexNDIndex(const size_t& vertex_index) const; + /** + * @brief Convert a flattened element index to its N-dimensional mesh index. + * + * The returned index denotes the lower mesh vertex of the element. Dimension zero + * is the fastest-varying dimension in the flattened representation. + * + * @param element_index Zero-based flattened element index in the interval + * `[0, product(nr_m[d] - 1))`. + * @return N-dimensional element index in which component `d` is in the interval + * `[0, nr_m[d] - 1)`. + * + * @pre Every dimension contains at least two mesh vertices. + * @see getElementIndex() + */ + KOKKOS_FUNCTION indices_t getElementNDIndex(const size_t& element_index) const; + + /** + * @brief Flatten an N-dimensional element index. + * + * This is the inverse of getElementNDIndex() for valid element indices and uses a + * dimension-zero-fastest ordering. + * + * @param element_nd_index N-dimensional index of an element's lower mesh vertex. + * Each component `d` must be in `[0, nr_m[d] - 1)`. + * @return Zero-based flattened element index. + * + * @see getElementNDIndex() + */ + KOKKOS_FUNCTION size_t getElementIndex(const indices_t& element_nd_index) const; + + /** + * @brief Return the mesh indices of every vertex belonging to an element. + * + * Vertex ordering follows the tensor-product binary convention used throughout the + * finite-element implementation: bit `d` of a local vertex number selects the lower + * (`0`) or upper (`1`) vertex in dimension `d`. + * + * @param element_nd_index N-dimensional index of the element's lower mesh vertex. + * @return Fixed-size list of the element vertex indices in local vertex order. + */ + KOKKOS_FUNCTION indices_list_t + getElementMeshVertexNDIndices(const indices_t& element_nd_index) const; + + /** + * @brief Return the physical coordinates of every vertex belonging to an element. + * + * Coordinates are computed from the mirrored structured-mesh geometry as + * `origin_m[d] + vertex_index[d] * hr_m[d]`. The returned points use the same local + * vertex ordering as getElementMeshVertexNDIndices(). + * + * @param element_nd_index N-dimensional index of the element's lower mesh vertex. + * @return Fixed-size list of physical vertex coordinates in local vertex order. + * + * @see getElementMeshVertexNDIndices() + */ + KOKKOS_FUNCTION vertex_points_t + getElementMeshVertexPoints(const indices_t& element_nd_index) const; + KOKKOS_FUNCTION size_t getLocalDOFIndex(const indices_t& elementNDIndex, const size_t& globalDOFIndex) const; + + /** + * @brief Return the global degree-of-freedom indices of a flattened element. + * + * This convenience overload first converts the flattened element index with + * getElementNDIndex() and then applies the existing N-dimensional DOF mapping. It + * allows kernels that iterate over the flattened `elementIndices` view to remain + * entirely within the device-safe interface. + * + * @param elementIndex Zero-based flattened element index. + * @return Global DOF indices in local element-DOF order. + * + * @see getElementNDIndex() + * @see getGlobalDOFIndices(const indices_t&) const + */ + KOKKOS_FUNCTION Vector getGlobalDOFIndices( + const size_t& elementIndex) const; KOKKOS_FUNCTION Vector getGlobalDOFIndices( const indices_t& elementNDIndex) const; + /** + * @brief Determine whether a global degree of freedom lies on the mesh boundary. + * + * A DOF is a boundary DOF when at least one index component is on the lower boundary + * (`0`) or upper boundary (`nr_m[d] - 1`). The predicate is used inside assembly + * kernels to apply or skip Dirichlet boundary contributions without accessing the + * host-side `LagrangeSpace` object. + * + * @param ndindex N-dimensional global mesh/DOF index. + * @return `true` if any component lies on a domain boundary; otherwise `false`. + */ + KOKKOS_FUNCTION bool isDOFOnBoundary(const indices_t& ndindex) const; + KOKKOS_FUNCTION T evaluateRefElementShapeFunction(const size_t& localDOF, const point_t& localPoint) const; KOKKOS_FUNCTION point_t evaluateRefElementShapeFunctionGradient( const size_t& localDOF, const point_t& localPoint) const; }; + /** + * @brief Create the device-safe snapshot captured by LagrangeSpace kernels. + * + * Copies the mesh extents, spacing, origin, and reference element into a non-owning value + * that can be captured by `KOKKOS_LAMBDA`. The returned value supplies device-side mesh + * indexing, physical-coordinate, DOF-mapping, boundary, and reference-element operations + * without retaining a pointer or reference to the parent `LagrangeSpace`. Views containing + * the element partition and field data are intentionally captured separately by each + * kernel. + * + * @return Independent, device-copyable snapshot of the current geometric and indexing + * state. + * + * @note Recreate the mirror after changing the corresponding host-side mesh or reference + * element state. The returned object does not synchronize later host changes. + */ DeviceStruct getDeviceMirror() const; private: diff --git a/src/FEM/LagrangeSpace.hpp b/src/FEM/LagrangeSpace.hpp index 3eeb5205f5..5823d17bde 100644 --- a/src/FEM/LagrangeSpace.hpp +++ b/src/FEM/LagrangeSpace.hpp @@ -72,9 +72,10 @@ namespace ippl { // while tagging upper boundary points such that they can be removed after. Kokkos::View points("npoints", npoints); Kokkos::View is_boundary("is_boundary", npoints); + const DeviceStruct space = getDeviceMirror(); Kokkos::parallel_reduce( "ComputePoints", npoints, - KOKKOS_CLASS_LAMBDA(const int i, int& local) { + KOKKOS_LAMBDA(const int i, int& local) { int idx = i; indices_t val; bool isBoundary = false; @@ -87,7 +88,7 @@ namespace ippl { } } is_boundary(i) = isBoundary; - points(i) = this->getElementIndex(val); + points(i) = space.getElementIndex(val); local += isBoundary; }, Kokkos::Sum(upperBoundaryPoints)); @@ -100,11 +101,12 @@ namespace ippl { Kokkos::View index("index"); if (elementsPerRank > 0) { + const auto elementIndicesView = elementIndices; Kokkos::parallel_for( - "CompactElementIndices", npoints, KOKKOS_CLASS_LAMBDA(const int i) { + "CompactElementIndices", npoints, KOKKOS_LAMBDA(const int i) { if (!is_boundary(i)) { - const size_t idx = Kokkos::atomic_fetch_add(&index(), 1); - elementIndices(idx) = points(i); + const size_t idx = Kokkos::atomic_fetch_add(&index(), 1); + elementIndicesView(idx) = points(i); } }); } @@ -441,17 +443,20 @@ namespace ippl { // start a timer IpplTimings::startTimer(evalAx_outer); + const DeviceStruct space = getDeviceMirror(); + const auto elementIndicesView = elementIndices; + // Loop over elements to compute contributions Kokkos::parallel_for( "Loop over elements", policy_type(0, elementIndices.extent(0)), - KOKKOS_CLASS_LAMBDA(const size_t index) { - const size_t elementIndex = elementIndices(index); + KOKKOS_LAMBDA(const size_t index) { + const size_t elementIndex = elementIndicesView(index); const Vector global_dofs = - this->LagrangeSpace::getGlobalDOFIndices(elementIndex); + space.getGlobalDOFIndices(elementIndex); Vector global_dof_ndindices; for (size_t i = 0; i < numElementDOFs; ++i) { - global_dof_ndindices[i] = this->getMeshVertexNDIndex(global_dofs[i]); + global_dof_ndindices[i] = space.getMeshVertexNDIndex(global_dofs[i]); } // local DOF indices (both i and j go from 0 to numDOFs-1 in the element) @@ -468,13 +473,13 @@ namespace ippl { // Handle boundary DOFs // If Zero Dirichlet BCs, skip this DOF // If Constant Dirichlet BCs, identity - if ((bcType == CONSTANT_FACE) && (this->isDOFOnBoundary(I_nd))) { + if ((bcType == CONSTANT_FACE) && (space.isDOFOnBoundary(I_nd))) { for (unsigned d = 0; d < Dim; ++d) { I_nd[d] = I_nd[d] - ldom[d].first() + nghost; } apply(resultView, I_nd) = apply(view, I_nd); continue; - } else if ((bcType == ZERO_FACE) && (this->isDOFOnBoundary(I_nd))) { + } else if ((bcType == ZERO_FACE) && (space.isDOFOnBoundary(I_nd))) { continue; } @@ -488,7 +493,7 @@ namespace ippl { // Skip boundary DOFs (Zero & Constant Dirichlet BCs) if (((bcType == ZERO_FACE) || (bcType == CONSTANT_FACE)) - && this->isDOFOnBoundary(J_nd)) { + && space.isDOFOnBoundary(J_nd)) { continue; } @@ -588,17 +593,20 @@ namespace ippl { using exec_space = typename Kokkos::View::execution_space; using policy_type = Kokkos::RangePolicy; + const DeviceStruct space = getDeviceMirror(); + const auto elementIndicesView = elementIndices; + // Loop over elements to compute contributions Kokkos::parallel_for( "Loop over elements", policy_type(0, elementIndices.extent(0)), - KOKKOS_CLASS_LAMBDA(const size_t index) { - const size_t elementIndex = elementIndices(index); + KOKKOS_LAMBDA(const size_t index) { + const size_t elementIndex = elementIndicesView(index); const Vector global_dofs = - this->LagrangeSpace::getGlobalDOFIndices(elementIndex); + space.getGlobalDOFIndices(elementIndex); Vector global_dof_ndindices; for (size_t i = 0; i < numElementDOFs; ++i) { - global_dof_ndindices[i] = this->getMeshVertexNDIndex(global_dofs[i]); + global_dof_ndindices[i] = space.getMeshVertexNDIndex(global_dofs[i]); } // local DOF indices @@ -615,13 +623,13 @@ namespace ippl { // Handle boundary DOFs // If Zero Dirichlet BCs, skip this DOF // If Constant Dirichlet BCs, identity - if ((bcType == CONSTANT_FACE) && (this->isDOFOnBoundary(I_nd))) { + if ((bcType == CONSTANT_FACE) && (space.isDOFOnBoundary(I_nd))) { for (unsigned d = 0; d < Dim; ++d) { I_nd[d] = I_nd[d] - ldom[d].first() + nghost; } apply(resultView, I_nd) = apply(view, I_nd); continue; - } else if ((bcType == ZERO_FACE) && (this->isDOFOnBoundary(I_nd))) { + } else if ((bcType == ZERO_FACE) && (space.isDOFOnBoundary(I_nd))) { continue; } @@ -639,7 +647,7 @@ namespace ippl { // Skip boundary DOFs (Zero & Constant Dirichlet BCs) if (((bcType == ZERO_FACE) || (bcType == CONSTANT_FACE)) - && this->isDOFOnBoundary(J_nd)) { + && space.isDOFOnBoundary(J_nd)) { continue; } @@ -732,17 +740,20 @@ namespace ippl { using exec_space = typename Kokkos::View::execution_space; using policy_type = Kokkos::RangePolicy; + const DeviceStruct space = getDeviceMirror(); + const auto elementIndicesView = elementIndices; + // Loop over elements to compute contributions Kokkos::parallel_for( "Loop over elements", policy_type(0, elementIndices.extent(0)), - KOKKOS_CLASS_LAMBDA(const size_t index) { - const size_t elementIndex = elementIndices(index); + KOKKOS_LAMBDA(const size_t index) { + const size_t elementIndex = elementIndicesView(index); const Vector global_dofs = - this->LagrangeSpace::getGlobalDOFIndices(elementIndex); + space.getGlobalDOFIndices(elementIndex); Vector global_dof_ndindices; for (size_t i = 0; i < numElementDOFs; ++i) { - global_dof_ndindices[i] = this->getMeshVertexNDIndex(global_dofs[i]); + global_dof_ndindices[i] = space.getMeshVertexNDIndex(global_dofs[i]); } // local DOF indices @@ -759,13 +770,13 @@ namespace ippl { // Handle boundary DOFs // If Zero Dirichlet BCs, skip this DOF // If Constant Dirichlet BCs, identity - if ((bcType == CONSTANT_FACE) && (this->isDOFOnBoundary(I_nd))) { + if ((bcType == CONSTANT_FACE) && (space.isDOFOnBoundary(I_nd))) { for (unsigned d = 0; d < Dim; ++d) { I_nd[d] = I_nd[d] - ldom[d].first() + nghost; } apply(resultView, I_nd) = apply(view, I_nd); continue; - } else if ((bcType == ZERO_FACE) && (this->isDOFOnBoundary(I_nd))) { + } else if ((bcType == ZERO_FACE) && (space.isDOFOnBoundary(I_nd))) { continue; } @@ -783,7 +794,7 @@ namespace ippl { // Skip boundary DOFs (Zero & Constant Dirichlet BCs) if (((bcType == ZERO_FACE) || (bcType == CONSTANT_FACE)) - && this->isDOFOnBoundary(J_nd)) { + && space.isDOFOnBoundary(J_nd)) { continue; } @@ -877,17 +888,20 @@ namespace ippl { using exec_space = typename Kokkos::View::execution_space; using policy_type = Kokkos::RangePolicy; + const DeviceStruct space = getDeviceMirror(); + const auto elementIndicesView = elementIndices; + // Loop over elements to compute contributions Kokkos::parallel_for( "Loop over elements", policy_type(0, elementIndices.extent(0)), - KOKKOS_CLASS_LAMBDA(const size_t index) { - const size_t elementIndex = elementIndices(index); + KOKKOS_LAMBDA(const size_t index) { + const size_t elementIndex = elementIndicesView(index); const Vector global_dofs = - this->LagrangeSpace::getGlobalDOFIndices(elementIndex); + space.getGlobalDOFIndices(elementIndex); Vector global_dof_ndindices; for (size_t i = 0; i < numElementDOFs; ++i) { - global_dof_ndindices[i] = this->getMeshVertexNDIndex(global_dofs[i]); + global_dof_ndindices[i] = space.getMeshVertexNDIndex(global_dofs[i]); } // local DOF indices @@ -904,13 +918,13 @@ namespace ippl { // Handle boundary DOFs // If Zero Dirichlet BCs, skip this DOF // If Constant Dirichlet BCs, identity - if ((bcType == CONSTANT_FACE) && (this->isDOFOnBoundary(I_nd))) { + if ((bcType == CONSTANT_FACE) && (space.isDOFOnBoundary(I_nd))) { for (unsigned d = 0; d < Dim; ++d) { I_nd[d] = I_nd[d] - ldom[d].first() + nghost; } apply(resultView, I_nd) = apply(view, I_nd); continue; - } else if ((bcType == ZERO_FACE) && (this->isDOFOnBoundary(I_nd))) { + } else if ((bcType == ZERO_FACE) && (space.isDOFOnBoundary(I_nd))) { continue; } @@ -924,7 +938,7 @@ namespace ippl { // Skip boundary DOFs (Zero & Constant Dirichlet BCs) if (((bcType == ZERO_FACE) || (bcType == CONSTANT_FACE)) - && this->isDOFOnBoundary(J_nd)) { + && space.isDOFOnBoundary(J_nd)) { continue; } @@ -1016,17 +1030,20 @@ namespace ippl { using exec_space = typename Kokkos::View::execution_space; using policy_type = Kokkos::RangePolicy; + const DeviceStruct space = getDeviceMirror(); + const auto elementIndicesView = elementIndices; + // Loop over elements to compute contributions Kokkos::parallel_for( "Loop over elements", policy_type(0, elementIndices.extent(0)), - KOKKOS_CLASS_LAMBDA(const size_t index) { - const size_t elementIndex = elementIndices(index); + KOKKOS_LAMBDA(const size_t index) { + const size_t elementIndex = elementIndicesView(index); const Vector global_dofs = - this->LagrangeSpace::getGlobalDOFIndices(elementIndex); + space.getGlobalDOFIndices(elementIndex); Vector global_dof_ndindices; for (size_t i = 0; i < numElementDOFs; ++i) { - global_dof_ndindices[i] = this->getMeshVertexNDIndex(global_dofs[i]); + global_dof_ndindices[i] = space.getMeshVertexNDIndex(global_dofs[i]); } // local DOF indices @@ -1043,13 +1060,13 @@ namespace ippl { // Handle boundary DOFs // If Zero Dirichlet BCs, skip this DOF // If Constant Dirichlet BCs, identity - if ((bcType == CONSTANT_FACE) && (this->isDOFOnBoundary(I_nd))) { + if ((bcType == CONSTANT_FACE) && (space.isDOFOnBoundary(I_nd))) { for (unsigned d = 0; d < Dim; ++d) { I_nd[d] = I_nd[d] - ldom[d].first() + nghost; } apply(resultView, I_nd) = 1.0; continue; - } else if ((bcType == ZERO_FACE) && (this->isDOFOnBoundary(I_nd))) { + } else if ((bcType == ZERO_FACE) && (space.isDOFOnBoundary(I_nd))) { continue; } @@ -1147,17 +1164,20 @@ namespace ippl { using exec_space = typename Kokkos::View::execution_space; using policy_type = Kokkos::RangePolicy; + const DeviceStruct space = getDeviceMirror(); + const auto elementIndicesView = elementIndices; + // Loop over elements to compute contributions Kokkos::parallel_for( "Loop over elements", policy_type(0, elementIndices.extent(0)), - KOKKOS_CLASS_LAMBDA(const size_t index) { - const size_t elementIndex = elementIndices(index); + KOKKOS_LAMBDA(const size_t index) { + const size_t elementIndex = elementIndicesView(index); const Vector global_dofs = - this->LagrangeSpace::getGlobalDOFIndices(elementIndex); + space.getGlobalDOFIndices(elementIndex); Vector global_dof_ndindices; for (size_t i = 0; i < numElementDOFs; ++i) { - global_dof_ndindices[i] = this->getMeshVertexNDIndex(global_dofs[i]); + global_dof_ndindices[i] = space.getMeshVertexNDIndex(global_dofs[i]); } // local DOF indices @@ -1174,13 +1194,13 @@ namespace ippl { // Handle boundary DOFs // If Zero Dirichlet BCs, skip this DOF // If Constant Dirichlet BCs, identity - if ((bcType == CONSTANT_FACE) && (this->isDOFOnBoundary(I_nd))) { + if ((bcType == CONSTANT_FACE) && (space.isDOFOnBoundary(I_nd))) { for (unsigned d = 0; d < Dim; ++d) { I_nd[d] = I_nd[d] - ldom[d].first() + nghost; } apply(resultView, I_nd) = apply(view, I_nd); continue; - } else if ((bcType == ZERO_FACE) && (this->isDOFOnBoundary(I_nd))) { + } else if ((bcType == ZERO_FACE) && (space.isDOFOnBoundary(I_nd))) { continue; } @@ -1267,17 +1287,20 @@ namespace ippl { using exec_space = typename Kokkos::View::execution_space; using policy_type = Kokkos::RangePolicy; + const DeviceStruct space = getDeviceMirror(); + const auto elementIndicesView = elementIndices; + // Loop over elements to compute contributions Kokkos::parallel_for( "Loop over elements", policy_type(0, elementIndices.extent(0)), - KOKKOS_CLASS_LAMBDA(const size_t index) { - const size_t elementIndex = elementIndices(index); + KOKKOS_LAMBDA(const size_t index) { + const size_t elementIndex = elementIndicesView(index); const Vector global_dofs = - this->LagrangeSpace::getGlobalDOFIndices(elementIndex); + space.getGlobalDOFIndices(elementIndex); Vector global_dof_ndindices; for (size_t i = 0; i < numElementDOFs; ++i) { - global_dof_ndindices[i] = this->getMeshVertexNDIndex(global_dofs[i]); + global_dof_ndindices[i] = space.getMeshVertexNDIndex(global_dofs[i]); } // local DOF indices (both i and j go from 0 to numDOFs-1 in the element) @@ -1292,7 +1315,7 @@ namespace ippl { I_nd = global_dof_ndindices[i]; // Skip if on a row of the matrix - if (this->isDOFOnBoundary(I_nd)) { + if (space.isDOFOnBoundary(I_nd)) { continue; } @@ -1305,7 +1328,7 @@ namespace ippl { J_nd = global_dof_ndindices[j]; // Contribute to lifting only if on a boundary DOF - if (this->isDOFOnBoundary(J_nd)) { + if (space.isDOFOnBoundary(J_nd)) { // get the appropriate index for the Kokkos view of the field for (unsigned d = 0; d < Dim; ++d) { J_nd[d] = J_nd[d] - ldom[d].first() + nghost; @@ -1374,17 +1397,21 @@ namespace ippl { // We work with a temporary field since we need to use field // to evaluate the load vector; then we assign temp to RHS field AtomicViewType atomic_view = temp_field.getView(); + const auto fieldView = field.getView(); using exec_space = typename Kokkos::View::execution_space; using policy_type = Kokkos::RangePolicy; + const DeviceStruct space = getDeviceMirror(); + const auto elementIndicesView = elementIndices; + // Loop over elements to compute contributions Kokkos::parallel_for( "Loop over elements", policy_type(0, elementIndices.extent(0)), - KOKKOS_CLASS_LAMBDA(size_t index) { - const size_t elementIndex = elementIndices(index); + KOKKOS_LAMBDA(size_t index) { + const size_t elementIndex = elementIndicesView(index); const Vector global_dofs = - this->LagrangeSpace::getGlobalDOFIndices(elementIndex); + space.getGlobalDOFIndices(elementIndex); size_t i, I; @@ -1393,11 +1420,11 @@ namespace ippl { I = global_dofs[i]; // TODO fix for higher order - auto dof_ndindex_I = this->getMeshVertexNDIndex(I); + auto dof_ndindex_I = space.getMeshVertexNDIndex(I); // Skip boundary DOFs (Zero and Constant Dirichlet BCs) if (((bcType == ZERO_FACE) || (bcType == CONSTANT_FACE)) - && (this->isDOFOnBoundary(dof_ndindex_I))) { + && (space.isDOFOnBoundary(dof_ndindex_I))) { continue; } @@ -1408,13 +1435,13 @@ namespace ippl { for (size_t j = 0; j < numElementDOFs; ++j) { // get field index corresponding to this DOF size_t J = global_dofs[j]; - auto dof_ndindex_J = this->getMeshVertexNDIndex(J); + auto dof_ndindex_J = space.getMeshVertexNDIndex(J); for (unsigned d = 0; d < Dim; ++d) { dof_ndindex_J[d] = dof_ndindex_J[d] - ldom[d].first() + nghost; } // get field value at DOF and interpolate to q_k - val += basis_q[k][j] * apply(field, dof_ndindex_J); + val += basis_q[k][j] * apply(fieldView, dof_ndindex_J); } contrib += w[k] * basis_q[k][i] * absDetDPhi * val; @@ -1483,13 +1510,16 @@ namespace ippl { using exec_space = typename Kokkos::View::execution_space; using policy_type = Kokkos::RangePolicy; + const DeviceStruct space = getDeviceMirror(); + const auto elementIndicesView = elementIndices; + // Loop over elements to compute contributions Kokkos::parallel_for( "Loop over elements", policy_type(0, elementIndices.extent(0)), - KOKKOS_CLASS_LAMBDA(size_t index) { - const size_t elementIndex = elementIndices(index); + KOKKOS_LAMBDA(size_t index) { + const size_t elementIndex = elementIndicesView(index); const Vector global_dofs = - this->LagrangeSpace::getGlobalDOFIndices(elementIndex); + space.getGlobalDOFIndices(elementIndex); size_t i, I; @@ -1498,7 +1528,7 @@ namespace ippl { I = global_dofs[i]; // TODO fix for higher order - auto dof_ndindex_I = this->getMeshVertexNDIndex(I); + auto dof_ndindex_I = space.getMeshVertexNDIndex(I); // calculate the contribution of this element T contrib = 0; @@ -1559,36 +1589,40 @@ namespace ippl { // Get domain information and ghost cells auto ldom = (u_h.getLayout()).getLocalNDIndex(); const int nghost = u_h.getNghost(); + const auto uView = u_h.getView(); using exec_space = typename Kokkos::View::execution_space; using policy_type = Kokkos::RangePolicy; + const DeviceStruct space = getDeviceMirror(); + const auto elementIndicesView = elementIndices; + // Loop over elements to compute contributions Kokkos::parallel_reduce( "Compute error over elements", policy_type(0, elementIndices.extent(0)), - KOKKOS_CLASS_LAMBDA(size_t index, double& local) { - const size_t elementIndex = elementIndices(index); + KOKKOS_LAMBDA(size_t index, double& local) { + const size_t elementIndex = elementIndicesView(index); const Vector global_dofs = - this->LagrangeSpace::getGlobalDOFIndices(elementIndex); + space.getGlobalDOFIndices(elementIndex); // contribution of this element to the error T contrib = 0; for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) { - T val_u_sol = u_sol(this->ref_element_m.localToGlobal( - this->getElementMeshVertexPoints(this->getElementNDIndex(elementIndex)), + T val_u_sol = u_sol(space.ref_element_m.localToGlobal( + space.getElementMeshVertexPoints(space.getElementNDIndex(elementIndex)), q[k])); T val_u_h = 0; for (size_t i = 0; i < numElementDOFs; ++i) { // get field index corresponding to this DOF size_t I = global_dofs[i]; - auto dof_ndindex_I = this->getMeshVertexNDIndex(I); + auto dof_ndindex_I = space.getMeshVertexNDIndex(I); for (unsigned d = 0; d < Dim; ++d) { dof_ndindex_I[d] = dof_ndindex_I[d] - ldom[d].first() + nghost; } // get field value at DOF and interpolate to q_k - val_u_h += basis_q[k][i] * apply(u_h, dof_ndindex_I); + val_u_h += basis_q[k][i] * apply(uView, dof_ndindex_I); } contrib += w[k] * Kokkos::pow(val_u_sol - val_u_h, 2) * absDetDPhi; @@ -1643,17 +1677,21 @@ namespace ippl { // Get domain information and ghost cells auto ldom = (u_h.getLayout()).getLocalNDIndex(); const int nghost = u_h.getNghost(); + const auto uView = u_h.getView(); using exec_space = typename Kokkos::View::execution_space; using policy_type = Kokkos::RangePolicy; + const DeviceStruct space = getDeviceMirror(); + const auto elementIndicesView = elementIndices; + // Loop over elements to compute contributions Kokkos::parallel_reduce( "Compute average over elements", policy_type(0, elementIndices.extent(0)), - KOKKOS_CLASS_LAMBDA(size_t index, double& local) { - const size_t elementIndex = elementIndices(index); + KOKKOS_LAMBDA(size_t index, double& local) { + const size_t elementIndex = elementIndicesView(index); const Vector global_dofs = - this->LagrangeSpace::getGlobalDOFIndices(elementIndex); + space.getGlobalDOFIndices(elementIndex); // contribution of this element to the error T contrib = 0; @@ -1662,13 +1700,13 @@ namespace ippl { for (size_t i = 0; i < numElementDOFs; ++i) { // get field index corresponding to this DOF size_t I = global_dofs[i]; - auto dof_ndindex_I = this->getMeshVertexNDIndex(I); + auto dof_ndindex_I = space.getMeshVertexNDIndex(I); for (unsigned d = 0; d < Dim; ++d) { dof_ndindex_I[d] = dof_ndindex_I[d] - ldom[d].first() + nghost; } // get field value at DOF and interpolate to q_k - val_u_h += basis_q[k][i] * apply(u_h, dof_ndindex_I); + val_u_h += basis_q[k][i] * apply(uView, dof_ndindex_I); } contrib += w[k] * val_u_h * absDetDPhi; @@ -1696,7 +1734,9 @@ namespace ippl { LagrangeSpace:: getDeviceMirror() const { DeviceStruct space_mirror; - space_mirror.nr_m = this->nr_m; + space_mirror.nr_m = this->nr_m; + space_mirror.hr_m = this->hr_m; + space_mirror.origin_m = this->origin_m; space_mirror.ref_element_m = this->ref_element_m; return space_mirror; } @@ -1706,6 +1746,87 @@ namespace ippl { // evaluateRefElementShapeFunction, and getMeshVertexNDIndex from the // parent class FiniteElementSpace get propagated here. + template + KOKKOS_FUNCTION typename LagrangeSpace::indices_t + LagrangeSpace::DeviceStruct::getElementNDIndex(const size_t& element_index) const { + size_t index = element_index; + indices_t element_nd_index; + const Vector cells_per_dim = nr_m - 1; + + size_t remaining_number_of_cells = 1; + for (const size_t num_cells : cells_per_dim) { + remaining_number_of_cells *= num_cells; + } + + for (int d = Dim - 1; d >= 0; --d) { + remaining_number_of_cells /= cells_per_dim[d]; + element_nd_index[d] = index / remaining_number_of_cells; + index -= element_nd_index[d] * remaining_number_of_cells; + } + + return element_nd_index; + } + + template + KOKKOS_FUNCTION size_t + LagrangeSpace::DeviceStruct::getElementIndex( + const indices_t& element_nd_index) const { + size_t element_index = 0; + const Vector cells_per_dim = nr_m - 1; + size_t remaining_number_of_cells = 1; + + for (unsigned d = 0; d < Dim; ++d) { + element_index += element_nd_index[d] * remaining_number_of_cells; + remaining_number_of_cells *= cells_per_dim[d]; + } + + return element_index; + } + + template + KOKKOS_FUNCTION typename LagrangeSpace::DeviceStruct::indices_list_t + LagrangeSpace::DeviceStruct::getElementMeshVertexNDIndices( + const indices_t& element_nd_index) const { + indices_list_t vertex_nd_indices; + + for (size_t i = 0; i < numElementVertices; ++i) { + vertex_nd_indices[i] = element_nd_index; + for (size_t d = 0; d < Dim; ++d) { + vertex_nd_indices[i][d] += (i >> d) & 1; + } + } + + return vertex_nd_indices; + } + + template + KOKKOS_FUNCTION typename LagrangeSpace::DeviceStruct::vertex_points_t + LagrangeSpace::DeviceStruct::getElementMeshVertexPoints( + const indices_t& element_nd_index) const { + vertex_points_t vertex_points; + const indices_list_t vertex_nd_indices = + getElementMeshVertexNDIndices(element_nd_index); + + for (size_t i = 0; i < numElementVertices; ++i) { + for (size_t d = 0; d < Dim; ++d) { + vertex_points[i][d] = vertex_nd_indices[i][d] * hr_m[d] + origin_m[d]; + } + } + + return vertex_points; + } + template KOKKOS_FUNCTION size_t @@ -1736,6 +1857,15 @@ namespace ippl { return 0; } + template + KOKKOS_FUNCTION Vector::DeviceStruct::numElementDOFs> + LagrangeSpace:: + DeviceStruct::getGlobalDOFIndices(const size_t& elementIndex) const { + return getGlobalDOFIndices(getElementNDIndex(elementIndex)); + } + template KOKKOS_FUNCTION Vector + KOKKOS_FUNCTION bool + LagrangeSpace::DeviceStruct::isDOFOnBoundary(const indices_t& ndindex) const { + for (size_t d = 0; d < Dim; ++d) { + if (ndindex[d] <= 0 || ndindex[d] >= nr_m[d] - 1) { + return true; + } + } + return false; + } + template KOKKOS_FUNCTION T From 22095a0f82ee1f137d16d34ebcb7b6b083267e8a Mon Sep 17 00:00:00 2001 From: Andreas Adelmann Date: Sun, 16 Aug 2026 12:11:04 +0200 Subject: [PATCH 2/9] Kokkos::pair for device-compatible subview Replace std::pair ranges passed to Kokkos::subview with Kokkos::pair in binning, particle sorting, and particle spatial layout code. Kokkos::pair is device-compatible, avoiding the NVHPC GH200 host/device warnings emitted for std::pair without changing the selected ranges. Validated with OpenMP and CUDA builds and unit tests on 1, 2, and 4 MPI ranks. The only four-rank failure is the pre-existing NedelecSpace partitioning limitation. LandauDamping FFT and FEM results remain consistent with the previous baseline. --- src/Interpolation/Binning.h | 10 +++++----- src/Particle/ParticleSort.h | 7 +++++-- src/Particle/ParticleSpatialLayout.hpp | 17 +++++++++++------ 3 files changed, 21 insertions(+), 13 deletions(-) diff --git a/src/Interpolation/Binning.h b/src/Interpolation/Binning.h index c67929c39f..dbaf75d869 100644 --- a/src/Interpolation/Binning.h +++ b/src/Interpolation/Binning.h @@ -148,7 +148,7 @@ namespace ippl { { auto offsets_zero = Kokkos::subview( - bin_offsets, std::make_pair(size_t(0), n_bins + 1)); + bin_offsets, Kokkos::pair{size_t(0), n_bins + 1}); Kokkos::deep_copy(ExecSpace(), offsets_zero, typename OffsetViewType::value_type(0)); } @@ -194,10 +194,10 @@ namespace ippl { IpplTimings::startTimer(offsetTimer); if (n_particles > 0) { - auto cursor_sub = Kokkos::subview(cursor, - std::make_pair(size_t(0), n_bins)); - auto offsets_sub = Kokkos::subview(bin_offsets, - std::make_pair(size_t(0), n_bins)); + auto cursor_sub = Kokkos::subview( + cursor, Kokkos::pair{size_t(0), n_bins}); + auto offsets_sub = Kokkos::subview( + bin_offsets, Kokkos::pair{size_t(0), n_bins}); Kokkos::deep_copy(ExecSpace(), cursor_sub, offsets_sub); Kokkos::parallel_for( diff --git a/src/Particle/ParticleSort.h b/src/Particle/ParticleSort.h index f5eaa3336f..4dfcadf458 100644 --- a/src/Particle/ParticleSort.h +++ b/src/Particle/ParticleSort.h @@ -274,8 +274,11 @@ namespace ippl { auto& bufs = ippl::detail::getDefaultBinSortBuffers(); bufs.ensureCapacity(n, /*n_bins_p1=*/1); - Kokkos::deep_copy(Kokkos::subview(bufs.permute(), std::make_pair(size_t(0), n)), - Kokkos::subview(permute_host, std::make_pair(size_t(0), n))); + Kokkos::deep_copy( + Kokkos::subview(bufs.permute(), + Kokkos::pair{size_t(0), n}), + Kokkos::subview(permute_host, + Kokkos::pair{size_t(0), n})); return bufs.permute(); } diff --git a/src/Particle/ParticleSpatialLayout.hpp b/src/Particle/ParticleSpatialLayout.hpp index e0590974c6..6106161917 100644 --- a/src/Particle/ParticleSpatialLayout.hpp +++ b/src/Particle/ParticleSpatialLayout.hpp @@ -154,8 +154,10 @@ namespace ippl { if (nDest > 0) { Kokkos::deep_copy( position_execution_space{}, - Kokkos::subview(destRanks_h_, std::make_pair(size_t(0), size_t(nDest))), - Kokkos::subview(destRanks_d_, std::make_pair(size_t(0), size_t(nDest)))); + Kokkos::subview(destRanks_h_, + Kokkos::pair{size_t(0), size_t(nDest)}), + Kokkos::subview(destRanks_d_, + Kokkos::pair{size_t(0), size_t(nDest)})); } // counts + offsets @@ -203,8 +205,9 @@ namespace ippl { if (count == 0) continue; const size_type begin = static_cast(sendOffsets_h_(rank)); - auto ids_sub = - Kokkos::subview(sendIds_d_, std::make_pair((size_t)begin, (size_t)(begin + count))); + auto ids_sub = Kokkos::subview( + sendIds_d_, Kokkos::pair{static_cast(begin), + static_cast(begin + count)}); requests.push_back(pc.sendToRank(rank, tag, ids_sub)); } @@ -538,8 +541,10 @@ namespace ippl { } Kokkos::deep_copy( - Kokkos::subview(neighbors_d_, std::make_pair(size_t(0), size_t(neighborSize))), - Kokkos::subview(neighbors_h, std::make_pair(size_t(0), size_t(neighborSize)))); + Kokkos::subview(neighbors_d_, + Kokkos::pair{size_t(0), size_t(neighborSize)}), + Kokkos::subview(neighbors_h, + Kokkos::pair{size_t(0), size_t(neighborSize)})); neighbors_dirty_ = false; } From e0e88e188d37f65db1cec0e2c53650be98f291cd Mon Sep 17 00:00:00 2001 From: Andreas Adelmann Date: Sun, 16 Aug 2026 15:21:23 +0200 Subject: [PATCH 3/9] Capture field views in differential expressions Field differential expressions were storing their input as a complete Field object. When an expression was copied into a Kokkos kernel closure, that also copied host-owned layout, boundary-condition, and lifetime state. NVHPC consequently generated device destruction paths for Field and reported calls to the host-only Field destructor. Keep the Field template parameter for mesh, dimension, and value-type traits, but store only its Kokkos view in the expression object. Construct gradient, divergence, Laplacian, curl, Hessian, Poisson, and triangular Laplacian expressions from getView(), and size their Expression metadata from the view rather than the complete Field. This establishes the intended host/device boundary: halo exchange and boundary-condition application remain in the host wrapper, while the deferred device expression retains only the field data handle and its finite-difference coefficients. Include meta_div and meta_poisson even though the CDash warning limit hid those instantiations, since they had the same ownership pattern. This removes the GH200 host-only Field destructor diagnostics without device-annotating Field lifetime operations or suppressing compiler warnings. Validation: full GCC/OpenMP and NVHPC CUDA/A100 builds completed successfully. Unit tests were exercised at 1, 2, and 4 MPI ranks on both backends. All affected tests pass; the four-rank runs retain only the established NedelecSpace fixture limitation, which supports three partitions (40/41 tests). --- src/Expression/IpplOperations.h | 44 +++++++++++++++++++---------- src/Field/FieldOperations.hpp | 10 +++---- src/PoissonSolvers/LaplaceHelpers.h | 44 ++++++++++++++++++----------- 3 files changed, 62 insertions(+), 36 deletions(-) diff --git a/src/Expression/IpplOperations.h b/src/Expression/IpplOperations.h index 1e02a6b8a2..1324be8f15 100644 --- a/src/Expression/IpplOperations.h +++ b/src/Expression/IpplOperations.h @@ -304,12 +304,14 @@ namespace ippl { struct meta_grad : public Expression< meta_grad, - sizeof(E) + sizeof(typename E::Mesh_t::vector_type[E::Mesh_t::Dimension])> { + sizeof(typename E::view_type) + + sizeof(typename E::Mesh_t::vector_type[E::Mesh_t::Dimension])> { constexpr static unsigned dim = E::dim; using value_type = typename E::value_type; KOKKOS_FUNCTION - meta_grad(const E& u, const typename E::Mesh_t::vector_type vectors[]) + meta_grad(const typename E::view_type& u, + const typename E::Mesh_t::vector_type vectors[]) : u_m(u) { for (unsigned d = 0; d < E::Mesh_t::Dimension; d++) { vectors_m[d] = vectors[d]; @@ -348,7 +350,8 @@ namespace ippl { private: using Mesh_t = typename E::Mesh_t; using vector_type = typename Mesh_t::vector_type; - const E u_m; + using view_type = typename E::view_type; + const view_type u_m; vector_type vectors_m[dim]; }; } // namespace detail @@ -362,11 +365,13 @@ namespace ippl { struct meta_div : public Expression< meta_div, - sizeof(E) + sizeof(typename E::Mesh_t::vector_type[E::Mesh_t::Dimension])> { + sizeof(typename E::view_type) + + sizeof(typename E::Mesh_t::vector_type[E::Mesh_t::Dimension])> { constexpr static unsigned dim = E::dim; KOKKOS_FUNCTION - meta_div(const E& u, const typename E::Mesh_t::vector_type vectors[]) + meta_div(const typename E::view_type& u, + const typename E::Mesh_t::vector_type vectors[]) : u_m(u) { for (unsigned d = 0; d < E::Mesh_t::Dimension; d++) { vectors_m[d] = vectors[d]; @@ -404,7 +409,8 @@ namespace ippl { private: using Mesh_t = typename E::Mesh_t; using vector_type = typename Mesh_t::vector_type; - const E u_m; + using view_type = typename E::view_type; + const view_type u_m; vector_type vectors_m[dim]; }; @@ -414,12 +420,14 @@ namespace ippl { template struct meta_laplace : public Expression, - sizeof(E) + sizeof(typename E::Mesh_t::vector_type)> { + sizeof(typename E::view_type) + + sizeof(typename E::Mesh_t::vector_type)> { constexpr static unsigned dim = E::dim; using value_type = typename E::value_type; KOKKOS_FUNCTION - meta_laplace(const E& u, const typename E::Mesh_t::vector_type& hvector) + meta_laplace(const typename E::view_type& u, + const typename E::Mesh_t::vector_type& hvector) : u_m(u) , hvector_m(hvector) {} @@ -456,7 +464,8 @@ namespace ippl { private: using Mesh_t = typename E::Mesh_t; using vector_type = typename Mesh_t::vector_type; - const E u_m; + using view_type = typename E::view_type; + const view_type u_m; const vector_type hvector_m; }; } // namespace detail @@ -469,11 +478,13 @@ namespace ippl { template struct meta_curl : public Expression, - sizeof(E) + 4 * sizeof(typename E::Mesh_t::vector_type)> { + sizeof(typename E::view_type) + + 4 * sizeof(typename E::Mesh_t::vector_type)> { constexpr static unsigned dim = E::dim; KOKKOS_FUNCTION - meta_curl(const E& u, const typename E::Mesh_t::vector_type& xvector, + meta_curl(const typename E::view_type& u, + const typename E::Mesh_t::vector_type& xvector, const typename E::Mesh_t::vector_type& yvector, const typename E::Mesh_t::vector_type& zvector, const typename E::Mesh_t::vector_type& hvector) @@ -501,7 +512,8 @@ namespace ippl { private: using Mesh_t = typename E::Mesh_t; using vector_type = typename Mesh_t::vector_type; - const E u_m; + using view_type = typename E::view_type; + const view_type u_m; const vector_type xvector_m; const vector_type yvector_m; const vector_type zvector_m; @@ -517,13 +529,14 @@ namespace ippl { template struct meta_hess : public Expression, - sizeof(E) + sizeof(typename E::view_type) + sizeof(typename E::Mesh_t::vector_type[E::Mesh_t::Dimension]) + sizeof(typename E::Mesh_t::vector_type)> { constexpr static unsigned dim = E::dim; KOKKOS_FUNCTION - meta_hess(const E& u, const typename E::Mesh_t::vector_type vectors[], + meta_hess(const typename E::view_type& u, + const typename E::Mesh_t::vector_type vectors[], const typename E::Mesh_t::vector_type& hvector) : u_m(u) , hvector_m(hvector) { @@ -546,8 +559,9 @@ namespace ippl { using Mesh_t = typename E::Mesh_t; using vector_type = typename Mesh_t::vector_type; using matrix_type = typename Mesh_t::matrix_type; + using view_type = typename E::view_type; - const E u_m; + const view_type u_m; vector_type vectors_m[dim]; const vector_type hvector_m; diff --git a/src/Field/FieldOperations.hpp b/src/Field/FieldOperations.hpp index c2b8dc2d0e..409723a0e8 100644 --- a/src/Field/FieldOperations.hpp +++ b/src/Field/FieldOperations.hpp @@ -25,7 +25,7 @@ namespace ippl { vectors[d] = 0; vectors[d][d] = 0.5 / mesh.getMeshSpacing(d); } - return detail::meta_grad(u, vectors); + return detail::meta_grad(u.getView(), vectors); } /*! @@ -49,7 +49,7 @@ namespace ippl { vectors[d] = 0; vectors[d][d] = 0.5 / mesh.getMeshSpacing(d); } - return detail::meta_div(u, vectors); + return detail::meta_div(u.getView(), vectors); } /*! @@ -70,7 +70,7 @@ namespace ippl { for (unsigned d = 0; d < Dim; d++) { hvector[d] = 1.0 / std::pow(mesh.getMeshSpacing(d), 2); } - return detail::meta_laplace(u, hvector); + return detail::meta_laplace(u.getView(), hvector); } /*! @@ -95,7 +95,7 @@ namespace ippl { zvector[2] = 1.0; typename mesh_type::vector_type hvector(0); hvector = mesh.getMeshSpacing(); - return detail::meta_curl(u, xvector, yvector, zvector, hvector); + return detail::meta_curl(u.getView(), xvector, yvector, zvector, hvector); } /*! @@ -121,6 +121,6 @@ namespace ippl { } auto hvector = mesh.getMeshSpacing(); - return detail::meta_hess(u, vectors, hvector); + return detail::meta_hess(u.getView(), vectors, hvector); } } // namespace ippl diff --git a/src/PoissonSolvers/LaplaceHelpers.h b/src/PoissonSolvers/LaplaceHelpers.h index 65c8a83f2c..1aeb71539a 100644 --- a/src/PoissonSolvers/LaplaceHelpers.h +++ b/src/PoissonSolvers/LaplaceHelpers.h @@ -8,11 +8,11 @@ namespace ippl { namespace detail { // Implements the poisson matrix acting on a d dimensional field template - struct meta_poisson : public Expression, sizeof(E)> { + struct meta_poisson : public Expression, sizeof(typename E::view_type)> { constexpr static unsigned dim = E::dim; KOKKOS_FUNCTION - meta_poisson(const E& u) + meta_poisson(const typename E::view_type& u) : u_m(u) {} template @@ -36,20 +36,23 @@ namespace ippl { } private: - const E u_m; + using view_type = typename E::view_type; + const view_type u_m; }; template struct meta_lower_laplace : public Expression, - sizeof(E) + sizeof(typename E::Mesh_t::vector_type) + sizeof(typename E::view_type) + + sizeof(typename E::Mesh_t::vector_type) + 2 * sizeof(typename E::Layout_t::NDIndex_t) + sizeof(unsigned)> { constexpr static unsigned dim = E::dim; using value_type = typename E::value_type; KOKKOS_FUNCTION - meta_lower_laplace(const E& u, const typename E::Mesh_t::vector_type& hvector, + meta_lower_laplace(const typename E::view_type& u, + const typename E::Mesh_t::vector_type& hvector, unsigned nghosts, const typename E::Layout_t::NDIndex_t& ldom, const typename E::Layout_t::NDIndex_t& domain) : u_m(u) @@ -92,8 +95,9 @@ namespace ippl { using Mesh_t = typename E::Mesh_t; using Layout_t = typename E::Layout_t; using vector_type = typename Mesh_t::vector_type; + using view_type = typename E::view_type; using domain_type = typename Layout_t::NDIndex_t; - const E u_m; + const view_type u_m; const vector_type hvector_m; const unsigned nghosts_m; const domain_type ldom_m; @@ -103,14 +107,16 @@ namespace ippl { template struct meta_upper_laplace : public Expression, - sizeof(E) + sizeof(typename E::Mesh_t::vector_type) + sizeof(typename E::view_type) + + sizeof(typename E::Mesh_t::vector_type) + 2 * sizeof(typename E::Layout_t::NDIndex_t) + sizeof(unsigned)> { constexpr static unsigned dim = E::dim; using value_type = typename E::value_type; KOKKOS_FUNCTION - meta_upper_laplace(const E& u, const typename E::Mesh_t::vector_type& hvector, + meta_upper_laplace(const typename E::view_type& u, + const typename E::Mesh_t::vector_type& hvector, unsigned nghosts, const typename E::Layout_t::NDIndex_t& ldom, const typename E::Layout_t::NDIndex_t& domain) : u_m(u) @@ -153,8 +159,9 @@ namespace ippl { using Mesh_t = typename E::Mesh_t; using Layout_t = typename E::Layout_t; using vector_type = typename Mesh_t::vector_type; + using view_type = typename E::view_type; using domain_type = typename Layout_t::NDIndex_t; - const E u_m; + const view_type u_m; const vector_type hvector_m; const unsigned nghosts_m; const domain_type ldom_m; @@ -164,12 +171,14 @@ namespace ippl { template struct meta_upper_and_lower_laplace : public Expression, - sizeof(E) + sizeof(typename E::Mesh_t::vector_type)> { + sizeof(typename E::view_type) + + sizeof(typename E::Mesh_t::vector_type)> { constexpr static unsigned dim = E::dim; using value_type = typename E::value_type; KOKKOS_FUNCTION - meta_upper_and_lower_laplace(const E& u, const typename E::Mesh_t::vector_type& hvector) + meta_upper_and_lower_laplace(const typename E::view_type& u, + const typename E::Mesh_t::vector_type& hvector) : u_m(u) , hvector_m(hvector) {} @@ -194,7 +203,8 @@ namespace ippl { private: using vector_type = typename E::Mesh_t::vector_type; - const E u_m; + using view_type = typename E::view_type; + const view_type u_m; const vector_type hvector_m; }; } // namespace detail @@ -211,7 +221,7 @@ namespace ippl { BConds& bcField = u.getFieldBC(); bcField.apply(u); - return detail::meta_poisson(u); + return detail::meta_poisson(u.getView()); } /*! @@ -247,7 +257,8 @@ namespace ippl { unsigned nghosts = u.getNghost(); const auto& ldom = layout.getLocalNDIndex(); const auto& domain = layout.getDomain(); - return detail::meta_lower_laplace(u, hvector, nghosts, ldom, domain); + return detail::meta_lower_laplace(u.getView(), hvector, nghosts, ldom, + domain); } /*! @@ -283,7 +294,8 @@ namespace ippl { unsigned nghosts = u.getNghost(); const auto& ldom = layout.getLocalNDIndex(); const auto& domain = layout.getDomain(); - return detail::meta_upper_laplace(u, hvector, nghosts, ldom, domain); + return detail::meta_upper_laplace(u.getView(), hvector, nghosts, ldom, + domain); } /*! @@ -315,7 +327,7 @@ namespace ippl { for (unsigned d = 0; d < Dim; d++) { hvector[d] = 1.0 / Kokkos::pow(mesh.getMeshSpacing(d), 2); } - return detail::meta_upper_and_lower_laplace(u, hvector); + return detail::meta_upper_and_lower_laplace(u.getView(), hvector); } /*! From 4f8ac1c526c515a90b947b4856779d6918e46ff6 Mon Sep 17 00:00:00 2001 From: Andreas Adelmann Date: Sun, 16 Aug 2026 15:21:59 +0200 Subject: [PATCH 4/9] Initialize Archive vector deserialize temporaries Archive reconstructs each Vector component by copying sizeof(T) bytes from the receive buffer into a local T before assigning it to the destination view. GCC cannot prove that the device-compatible copyBytes loop initializes every byte and emits -Wmaybe-uninitialized for each template instantiation. Value-initialize the local in both vector deserialize overloads, including the offset variant. copyBytes still overwrites the complete object representation, so serialized data and communication semantics are unchanged; initialization only provides a defined starting state that is visible to compiler data-flow analysis. This removes the OpenMP CDash warning family without compiler switches, diagnostic suppression, or a host-only memcpy dependency. Validation: full GCC/OpenMP and NVHPC CUDA/A100 builds completed successfully. Unit tests were exercised at 1, 2, and 4 MPI ranks on both backends. All affected tests pass; the four-rank runs retain only the established NedelecSpace fixture limitation, which supports three partitions (40/41 tests). --- src/Communicate/Archive.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Communicate/Archive.hpp b/src/Communicate/Archive.hpp index 2240e05264..a719a38af4 100644 --- a/src/Communicate/Archive.hpp +++ b/src/Communicate/Archive.hpp @@ -337,7 +337,7 @@ namespace ippl { "Archive::deserialize()", mdrange_t({0, 0}, {(long int)nrecvs, Dim}), KOKKOS_LAMBDA(const size_type i, const size_t d) { const char* src = base + (Dim * i + d) * size + readpos; - T value; + T value{}; char* dst = reinterpret_cast(&value); copyBytes(dst, src, size); view.data()[i](d) = value; @@ -394,7 +394,7 @@ namespace ippl { "Archive::deserialize(offset, vector)", mdrange_t({0, 0}, {(long int)nrecvs, Dim}), KOKKOS_LAMBDA(const size_type i, const size_t d) { const char* src = base + (Dim * i + d) * size + readpos; - T value; + T value{}; char* dst = reinterpret_cast(&value); copyBytes(dst, src, size); view.data()[offset + i](d) = value; From 4f652c02ebb027bbd5474eab8081aedfeab6c7a3 Mon Sep 17 00:00:00 2001 From: Andreas Adelmann Date: Mon, 17 Aug 2026 11:27:38 +0200 Subject: [PATCH 5/9] Capture a lightweight device mirror in Nedelec kernels NedelecSpace is a host-side owner whose FieldLayout member carries MPI decomposition metadata and host-only lifetime operations. Capturing the complete object in KOKKOS_CLASS_LAMBDA forces CUDA to form device closure construction and destruction paths for that host-owned state, which is why the GH200 build diagnoses FieldLayout destruction from device code. Introduce a documented DeviceStruct snapshot containing only immutable, device-copyable mesh geometry and the reference element. Reproduce the element-index conversion, vertex geometry, edge-DOF mapping, ghosted FEMVector indexing, boundary classification, and basis/curl evaluation needed inside kernels, and expose getDeviceMirror() as the explicit host-to-device architectural boundary. Change every affected assembly, reconstruction, error, and point-compaction kernel to capture the mirror and required Kokkos views by value with KOKKOS_LAMBDA. The device closure no longer owns or references FieldLayout, MPI state, quadrature ownership, or the parent NedelecSpace object. The detailed header documentation records the mirror lifetime, synchronization contract, exclusions, index ordering, and extension rules. Validated with warning-free clean OpenMP and CUDA sm_80 builds. The complete unit suite passes 41/41 at one and two ranks on both backends; at four ranks, only the established Nedelec fixture whose three-cell domain cannot be split into four partitions fails. --- src/FEM/NedelecSpace.h | 146 ++++++++++++++ src/FEM/NedelecSpace.hpp | 413 +++++++++++++++++++++++++++++++++++---- 2 files changed, 522 insertions(+), 37 deletions(-) diff --git a/src/FEM/NedelecSpace.h b/src/FEM/NedelecSpace.h index f99745c3d1..af455733c4 100644 --- a/src/FEM/NedelecSpace.h +++ b/src/FEM/NedelecSpace.h @@ -392,6 +392,152 @@ namespace ippl { */ KOKKOS_FUNCTION int getBoundarySide(const size_t& dofIdx) const; + /** + * @brief Non-owning, device-copyable snapshot used by Nedelec assembly kernels. + * + * `NedelecSpace` is a host-side owner. In addition to the finite-element geometry, it + * stores `layout_m`, which owns MPI-decomposition metadata and has host-only lifetime + * operations. Capturing the complete space with `KOKKOS_CLASS_LAMBDA` therefore places + * host-owned state in the kernel closure. CUDA compilers must generate construction and + * destruction paths for that closure and can consequently diagnose an illegal device + * call to `FieldLayout::~FieldLayout`, even though the kernel never reads the layout. + * + * `DeviceStruct` is the explicit boundary between that host owner and device execution. + * It snapshots only the immutable values required by Nedelec kernels: mesh vertex counts, + * mesh spacing, physical origin, and the device-copyable reference element. From those + * values it reproduces flattened/N-dimensional element conversion, element geometry, + * edge-degree-of-freedom mappings, FEMVector indices, boundary detection, and Nedelec + * basis/curl evaluation. + * + * A host operation obtains a fresh value with getDeviceMirror() and captures it by value + * in `KOKKOS_LAMBDA`. Kokkos views holding element indices, coefficients, and output data + * are captured separately. The mirror performs no allocation, MPI communication, field + * copy, or synchronization, and it deliberately excludes `layout_m`, mesh and quadrature + * references, FEMVector ownership, and every other host-oriented object. + * + * "Lightweight" means this is a read-only value snapshot, not another finite-element + * space and not an automatically synchronized replica. If the host mesh or reference + * element changes, callers must obtain a new mirror before launching a kernel. When a + * future kernel needs more Nedelec functionality, add only the smallest device-copyable + * state and algorithm here; do not restore capture of the parent `NedelecSpace`. + * + * Element flattening uses dimension zero as the fastest-varying dimension. Nedelec DOFs + * retain the edge ordering used by the host class, so global and ghosted FEMVector index + * calculations remain identical on host and device. + */ + struct DeviceStruct { + static constexpr unsigned numElementDOFs = NedelecSpace::numElementDOFs; + static constexpr unsigned numElementVertices = NedelecSpace::numElementVertices; + + using indices_list_t = Vector; + using vertex_points_t = Vector; + + Vector nr_m; ///< Number of mesh vertices in each dimension. + Vector hr_m; ///< Uniform mesh spacing in each dimension. + Vector origin_m; ///< Physical coordinate of mesh vertex index zero. + ElementType ref_element_m; ///< Device-copyable reference-element description. + + /** + * @brief Convert a flattened element index to its N-dimensional mesh index. + * @param elementIndex Zero-based flattened element index. + * @return N-dimensional index of the element's lower mesh vertex. + */ + KOKKOS_FUNCTION indices_t getElementNDIndex(const size_t& elementIndex) const; + + /** + * @brief Flatten an N-dimensional element index. + * @param elementIndex Index of the element's lower mesh vertex. + * @return Zero-based element index with dimension zero varying fastest. + */ + KOKKOS_FUNCTION size_t getElementIndex(const indices_t& elementIndex) const; + + /** + * @brief Return the mesh indices of all vertices belonging to an element. + * @param elementIndex Index of the element's lower mesh vertex. + * @return Vertex indices in the reference element's binary local ordering. + */ + KOKKOS_FUNCTION indices_list_t + getElementMeshVertexNDIndices(const indices_t& elementIndex) const; + + /** + * @brief Return physical coordinates of all vertices belonging to an element. + * @param elementIndex Index of the element's lower mesh vertex. + * @return Physical vertex coordinates in local reference-element order. + */ + KOKKOS_FUNCTION vertex_points_t + getElementMeshVertexPoints(const indices_t& elementIndex) const; + + /** + * @brief Map a flattened element to its global edge DOF indices. + * @param elementIndex Zero-based flattened element index. + * @return Global edge indices in local Nedelec basis order. + */ + KOKKOS_FUNCTION Vector getGlobalDOFIndices( + const size_t& elementIndex) const; + + /** + * @brief Map an N-dimensional element to its global edge DOF indices. + * @param elementIndex Index of the element's lower mesh vertex. + * @return Global edge indices in local Nedelec basis order. + */ + KOKKOS_FUNCTION Vector getGlobalDOFIndices( + const indices_t& elementIndex) const; + + /** + * @brief Map a flattened element to ghosted FEMVector storage indices. + * @param elementIndex Zero-based flattened element index. + * @param ldom Local mesh domain owned by the current MPI rank. + * @return Indices into the rank-local FEMVector including its ghost layer. + */ + KOKKOS_FUNCTION Vector getFEMVectorDOFIndices( + const size_t& elementIndex, NDIndex ldom) const; + + /** + * @brief Map an N-dimensional element to ghosted FEMVector storage indices. + * @param elementIndex Index of the element's lower mesh vertex. + * @param ldom Local mesh domain owned by the current MPI rank. + * @return Indices into the rank-local FEMVector including its ghost layer. + */ + KOKKOS_FUNCTION Vector getFEMVectorDOFIndices( + indices_t elementIndex, NDIndex ldom) const; + + /** + * @brief Evaluate a local Nedelec basis function in the reference element. + * @param localDOF Local edge-basis index. + * @param localPoint Point in reference-element coordinates. + * @return Vector value of the requested basis function. + */ + KOKKOS_FUNCTION point_t evaluateRefElementShapeFunction( + const size_t& localDOF, const point_t& localPoint) const; + + /** + * @brief Evaluate the curl of a local Nedelec basis function. + * @param localDOF Local edge-basis index. + * @param localPoint Point in reference-element coordinates. + * @return Curl in the representation used by the host NedelecSpace. + */ + KOKKOS_FUNCTION point_t evaluateRefElementShapeFunctionCurl( + const size_t& localDOF, const point_t& localPoint) const; + + /** + * @brief Test whether a global edge DOF lies on the physical mesh boundary. + * @param dofIdx Global Nedelec edge index. + * @return `true` when the edge belongs to any boundary face. + */ + KOKKOS_FUNCTION bool isDOFOnBoundary(const size_t& dofIdx) const; + }; + + /** + * @brief Create the device-safe snapshot captured by Nedelec kernels. + * + * The returned value is independent of the host object's lifetime and contains no + * `FieldLayout` or MPI ownership. It is valid until a kernel needs mesh or reference- + * element state newer than the state copied by this call. + * + * @return Device-copyable geometry and edge-indexing snapshot. + */ + DeviceStruct getDeviceMirror() const; + private: /** * @brief Implementation of the \c NedelecSpace::createFEMVector diff --git a/src/FEM/NedelecSpace.hpp b/src/FEM/NedelecSpace.hpp index 5fefc0f1fb..71d565704e 100644 --- a/src/FEM/NedelecSpace.hpp +++ b/src/FEM/NedelecSpace.hpp @@ -89,9 +89,10 @@ namespace ippl { // while tagging upper boundary points such that they can be removed after. Kokkos::View points("npoints", npoints); Kokkos::View is_boundary("is_boundary", npoints); + const DeviceStruct space = getDeviceMirror(); Kokkos::parallel_reduce( "ComputePoints", npoints, - KOKKOS_CLASS_LAMBDA(const int i, int& local) { + KOKKOS_LAMBDA(const int i, int& local) { int idx = i; indices_t val; bool isBoundary = false; @@ -104,7 +105,7 @@ namespace ippl { } } is_boundary(i) = isBoundary; - points(i) = this->getElementIndex(val); + points(i) = space.getElementIndex(val); local += isBoundary; }, Kokkos::Sum(upperBoundaryPoints)); @@ -114,14 +115,15 @@ namespace ippl { // with the tagged upper boundary points removed. int elementsPerRank = npoints - upperBoundaryPoints; elementIndices = Kokkos::View("i", elementsPerRank); + auto elementIndicesView = elementIndices; Kokkos::View index("index"); if (elementsPerRank > 0) { Kokkos::parallel_for( - "CompactElementIndices", npoints, KOKKOS_CLASS_LAMBDA(const int i) { + "CompactElementIndices", npoints, KOKKOS_LAMBDA(const int i) { if (!is_boundary(i)) { const size_t idx = Kokkos::atomic_fetch_add(&index(), 1); - elementIndices(idx) = points(i); + elementIndicesView(idx) = points(i); } }); } @@ -348,6 +350,331 @@ namespace ippl { } + /////////////////////////////////////////////////////////////////////// + /// Device-safe kernel snapshot /////////////////////////////////////// + /////////////////////////////////////////////////////////////////////// + + template + typename NedelecSpace::DeviceStruct + NedelecSpace::getDeviceMirror() const { + DeviceStruct space; + space.nr_m = this->nr_m; + space.hr_m = this->hr_m; + space.origin_m = this->origin_m; + space.ref_element_m = this->ref_element_m; + return space; + } + + template + KOKKOS_FUNCTION typename NedelecSpace::indices_t + NedelecSpace::DeviceStruct::getElementNDIndex(const size_t& elementIndex) const { + size_t index = elementIndex; + indices_t elementNDIndex; + const Vector cellsPerDim = nr_m - 1; + + size_t remainingCells = 1; + for (const size_t cells : cellsPerDim) { + remainingCells *= cells; + } + + for (int d = Dim - 1; d >= 0; --d) { + remainingCells /= cellsPerDim[d]; + elementNDIndex[d] = index / remainingCells; + index -= elementNDIndex[d] * remainingCells; + } + return elementNDIndex; + } + + template + KOKKOS_FUNCTION size_t + NedelecSpace::DeviceStruct::getElementIndex(const indices_t& elementIndex) const { + size_t flatIndex = 0; + const Vector cellsPerDim = nr_m - 1; + size_t stride = 1; + + for (unsigned d = 0; d < Dim; ++d) { + flatIndex += elementIndex[d] * stride; + stride *= cellsPerDim[d]; + } + return flatIndex; + } + + template + KOKKOS_FUNCTION typename NedelecSpace::DeviceStruct::indices_list_t + NedelecSpace::DeviceStruct::getElementMeshVertexNDIndices( + const indices_t& elementIndex) const { + indices_list_t vertices; + for (size_t i = 0; i < numElementVertices; ++i) { + vertices[i] = elementIndex; + for (size_t d = 0; d < Dim; ++d) { + vertices[i][d] += (i >> d) & 1; + } + } + return vertices; + } + + template + KOKKOS_FUNCTION typename NedelecSpace::DeviceStruct::vertex_points_t + NedelecSpace::DeviceStruct::getElementMeshVertexPoints( + const indices_t& elementIndex) const { + vertex_points_t points; + const indices_list_t vertices = getElementMeshVertexNDIndices(elementIndex); + for (size_t i = 0; i < numElementVertices; ++i) { + for (size_t d = 0; d < Dim; ++d) { + points[i][d] = vertices[i][d] * hr_m[d] + origin_m[d]; + } + } + return points; + } + + template + KOKKOS_FUNCTION Vector::DeviceStruct::numElementDOFs> + NedelecSpace::DeviceStruct::getGlobalDOFIndices( + const indices_t& elementIndex) const { + Vector globalDOFs(0); + Vector stride(1); + + if constexpr (Dim == 2) { + const size_t nx = nr_m[0]; + stride(1) = 2 * nx - 1; + } else if constexpr (Dim == 3) { + const size_t nx = nr_m[0]; + const size_t ny = nr_m[1]; + stride(1) = 2 * nx - 1; + stride(2) = 3 * nx * ny - nx - ny; + } + + const size_t nx = nr_m[0]; + globalDOFs(0) = stride.dot(elementIndex); + globalDOFs(1) = globalDOFs(0) + nx - 1; + globalDOFs(2) = globalDOFs(1) + nx; + globalDOFs(3) = globalDOFs(1) + 1; + + if constexpr (Dim == 3) { + const size_t ny = nr_m[1]; + globalDOFs(4) = stride(2) * elementIndex(2) + 2 * nx * ny - nx - ny + + elementIndex(1) * nx + elementIndex(0); + globalDOFs(5) = globalDOFs(4) + 1; + globalDOFs(6) = globalDOFs(4) + nx + 1; + globalDOFs(7) = globalDOFs(4) + nx; + globalDOFs(8) = globalDOFs(0) + 3 * nx * ny - nx - ny; + globalDOFs(9) = globalDOFs(8) + nx - 1; + globalDOFs(10) = globalDOFs(9) + nx; + globalDOFs(11) = globalDOFs(9) + 1; + } + return globalDOFs; + } + + template + KOKKOS_FUNCTION Vector::DeviceStruct::numElementDOFs> + NedelecSpace::DeviceStruct::getGlobalDOFIndices( + const size_t& elementIndex) const { + return getGlobalDOFIndices(getElementNDIndex(elementIndex)); + } + + template + KOKKOS_FUNCTION Vector::DeviceStruct::numElementDOFs> + NedelecSpace::DeviceStruct::getFEMVectorDOFIndices(indices_t elementIndex, + NDIndex ldom) const { + Vector vectorDOFs(0); + + elementIndex -= ldom.first(); + elementIndex += 1; + + indices_t extent = ldom.last() - ldom.first(); + extent += 3; // Include the final owned point and one ghost point on each side. + + Vector stride(1); + if constexpr (Dim == 2) { + const size_t nx = extent[0]; + stride(1) = 2 * nx - 1; + } else if constexpr (Dim == 3) { + const size_t nx = extent[0]; + const size_t ny = extent[1]; + stride(1) = 2 * nx - 1; + stride(2) = 3 * nx * ny - nx - ny; + } + + const size_t nx = extent[0]; + vectorDOFs(0) = stride.dot(elementIndex); + vectorDOFs(1) = vectorDOFs(0) + nx - 1; + vectorDOFs(2) = vectorDOFs(1) + nx; + vectorDOFs(3) = vectorDOFs(1) + 1; + + if constexpr (Dim == 3) { + const size_t ny = extent[1]; + vectorDOFs(4) = stride(2) * elementIndex(2) + 2 * nx * ny - nx - ny + + elementIndex(1) * nx + elementIndex(0); + vectorDOFs(5) = vectorDOFs(4) + 1; + vectorDOFs(6) = vectorDOFs(4) + nx + 1; + vectorDOFs(7) = vectorDOFs(4) + nx; + vectorDOFs(8) = vectorDOFs(0) + 3 * nx * ny - nx - ny; + vectorDOFs(9) = vectorDOFs(8) + nx - 1; + vectorDOFs(10) = vectorDOFs(9) + nx; + vectorDOFs(11) = vectorDOFs(9) + 1; + } + return vectorDOFs; + } + + template + KOKKOS_FUNCTION Vector::DeviceStruct::numElementDOFs> + NedelecSpace::DeviceStruct::getFEMVectorDOFIndices( + const size_t& elementIndex, NDIndex ldom) const { + return getFEMVectorDOFIndices(getElementNDIndex(elementIndex), ldom); + } + + template + KOKKOS_FUNCTION typename NedelecSpace::point_t + NedelecSpace::DeviceStruct::evaluateRefElementShapeFunction( + const size_t& localDOF, const point_t& localPoint) const { + assert(localDOF < numElementDOFs && "The local edge index is invalid"); + assert(ref_element_m.isPointInRefElement(localPoint) + && "Point is not in reference element"); + + point_t result(0); + if constexpr (Dim == 2) { + const T x = localPoint(0); + const T y = localPoint(1); + switch (localDOF) { + case 0: result(0) = 1 - y; break; + case 1: result(1) = 1 - x; break; + case 2: result(0) = y; break; + case 3: result(1) = x; break; + } + } else if constexpr (Dim == 3) { + const T x = localPoint(0); + const T y = localPoint(1); + const T z = localPoint(2); + switch (localDOF) { + case 0: result(0) = y * z - y - z + 1; break; + case 1: result(1) = x * z - x - z + 1; break; + case 2: result(0) = y * (1 - z); break; + case 3: result(1) = x * (1 - z); break; + case 4: result(2) = x * y - x - y + 1; break; + case 5: result(2) = x * (1 - y); break; + case 6: result(2) = x * y; break; + case 7: result(2) = y * (1 - x); break; + case 8: result(0) = z * (1 - y); break; + case 9: result(1) = z * (1 - x); break; + case 10: result(0) = y * z; break; + case 11: result(1) = x * z; break; + } + } + return result; + } + + template + KOKKOS_FUNCTION typename NedelecSpace::point_t + NedelecSpace::DeviceStruct::evaluateRefElementShapeFunctionCurl( + const size_t& localDOF, const point_t& localPoint) const { + point_t result(0); + if constexpr (Dim == 2) { + switch (localDOF) { + case 0: result(0) = 1; break; + case 1: result(0) = -1; break; + case 2: result(0) = -1; break; + case 3: result(0) = 1; break; + } + } else { + const T x = localPoint(0); + const T y = localPoint(1); + const T z = localPoint(2); + switch (localDOF) { + case 0: result(0) = 0; result(1) = -1 + y; result(2) = 1 - z; break; + case 1: result(0) = 1 - x; result(1) = 0; result(2) = -1 + z; break; + case 2: result(0) = 0; result(1) = -y; result(2) = -1 + z; break; + case 3: result(0) = x; result(1) = 0; result(2) = 1 - z; break; + case 4: result(0) = -1 + x; result(1) = 1 - y; result(2) = 0; break; + case 5: result(0) = -x; result(1) = -1 + y; result(2) = 0; break; + case 6: result(0) = x; result(1) = -y; result(2) = 0; break; + case 7: result(0) = 1 - x; result(1) = y; result(2) = 0; break; + case 8: result(0) = 0; result(1) = 1 - y; result(2) = z; break; + case 9: result(0) = -1 + x; result(1) = 0; result(2) = -z; break; + case 10: result(0) = 0; result(1) = y; result(2) = -z; break; + case 11: result(0) = -x; result(1) = 0; result(2) = z; break; + } + } + return result; + } + + template + KOKKOS_FUNCTION bool + NedelecSpace::DeviceStruct::isDOFOnBoundary(const size_t& dofIdx) const { + bool onBoundary = false; + if constexpr (Dim == 2) { + const size_t nx = nr_m[0]; + const size_t ny = nr_m[1]; + onBoundary = onBoundary || dofIdx < nx - 1; + onBoundary = onBoundary || dofIdx > nx * (ny - 1) + ny * (nx - 1) - nx; + onBoundary = onBoundary + || (dofIdx >= nx - 1 && (dofIdx - (nx - 1)) % (2 * nx - 1) == 0); + onBoundary = onBoundary + || (dofIdx >= 2 * nx - 2 + && (dofIdx - 2 * nx + 2) % (2 * nx - 1) == 0); + } else if constexpr (Dim == 3) { + const size_t nx = nr_m[0]; + const size_t ny = nr_m[1]; + const size_t nz = nr_m[2]; + const size_t planeSize = nx * (ny - 1) + ny * (nx - 1) + nx * ny; + const size_t zOffset = dofIdx / planeSize; + + if (dofIdx - planeSize * zOffset >= nx * (ny - 1) + ny * (nx - 1)) { + const size_t offset = dofIdx - planeSize * zOffset + - (nx * (ny - 1) + ny * (nx - 1)); + const size_t yOffset = offset / nx; + const size_t xOffset = offset % nx; + onBoundary = onBoundary || yOffset == 0 || yOffset == ny - 1; + onBoundary = onBoundary || xOffset == 0 || xOffset == nx - 1; + } else { + onBoundary = onBoundary || zOffset == 0 || zOffset == nz - 1; + const size_t offset = dofIdx - planeSize * zOffset; + const size_t yOffset = offset / (2 * nx - 1); + size_t xOffset = offset - (2 * nx - 1) * yOffset; + + if (xOffset < nx - 1) { + onBoundary = onBoundary || yOffset == 0 || yOffset == ny - 1; + } else { + xOffset -= nx - 1; + onBoundary = onBoundary || xOffset == 0 || xOffset == nx - 1; + } + } + } + return onBoundary; + } + + template KOKKOS_FUNCTION typename NedelecSpace::point_t @@ -438,19 +765,22 @@ namespace ippl { IpplTimings::TimerRef timerAxLoop = IpplTimings::getTimer("Ax Loop"); IpplTimings::startTimer(timerAxLoop); + const DeviceStruct space = getDeviceMirror(); + auto elementIndicesView = elementIndices; + // Loop over elements to compute contributions Kokkos::parallel_for( "Loop over elements", policy_type(0, elementIndices.extent(0)), - KOKKOS_CLASS_LAMBDA(const size_t index) { - const size_t elementIndex = elementIndices(index); + KOKKOS_LAMBDA(const size_t index) { + const size_t elementIndex = elementIndicesView(index); // Here we now retrieve the global DOF indices and their // position inside of the FEMVector const Vector global_dofs = - this->NedelecSpace::getGlobalDOFIndices(elementIndex); + space.getGlobalDOFIndices(elementIndex); const Vector vectorIndices = - this->getFEMVectorDOFIndices(elementIndex, ldom); + space.getFEMVectorDOFIndices(elementIndex, ldom); // local DOF indices @@ -464,7 +794,7 @@ namespace ippl { I = global_dofs[i]; // Skip boundary DOFs (Zero Dirichlet BCs) - if (this->isDOFOnBoundary(I)) { + if (space.isDOFOnBoundary(I)) { continue; } @@ -472,7 +802,7 @@ namespace ippl { J = global_dofs[j]; // Skip boundary DOFs (Zero Dirichlet BCs) - if (this->isDOFOnBoundary(J)) { + if (space.isDOFOnBoundary(J)) { continue; } @@ -547,22 +877,25 @@ namespace ippl { using exec_space = typename Kokkos::View::execution_space; using policy_type = Kokkos::RangePolicy; + const DeviceStruct space = getDeviceMirror(); + auto elementIndicesView = elementIndices; + // Loop over elements to compute contributions Kokkos::parallel_for( "Loop over elements", policy_type(0, elementIndices.extent(0)), - KOKKOS_CLASS_LAMBDA(size_t index) { - const size_t elementIndex = elementIndices(index); + KOKKOS_LAMBDA(size_t index) { + const size_t elementIndex = elementIndicesView(index); const Vector global_dofs = - this->NedelecSpace::getGlobalDOFIndices(elementIndex); + space.getGlobalDOFIndices(elementIndex); const Vector vectorIndices = - this->getFEMVectorDOFIndices(elementIndex, ldom); + space.getFEMVectorDOFIndices(elementIndex, ldom); size_t i; for (i = 0; i < numElementDOFs; ++i) { size_t I = global_dofs[i]; - if (this->isDOFOnBoundary(I)) { + if (space.isDOFOnBoundary(I)) { continue; } @@ -645,18 +978,20 @@ namespace ippl { using exec_space = typename Kokkos::View::execution_space; using policy_type = Kokkos::RangePolicy; - + + const DeviceStruct space = getDeviceMirror(); + auto elementIndicesView = elementIndices; // Loop over elements to compute contributions Kokkos::parallel_for( "Loop over elements", policy_type(0, elementIndices.extent(0)), - KOKKOS_CLASS_LAMBDA(size_t index) { - const size_t elementIndex = elementIndices(index); + KOKKOS_LAMBDA(size_t index) { + const size_t elementIndex = elementIndicesView(index); const Vector global_dofs = - this->NedelecSpace::getGlobalDOFIndices(elementIndex); + space.getGlobalDOFIndices(elementIndex); const Vector vectorIndices = - this->getFEMVectorDOFIndices(elementIndex, ldom); + space.getFEMVectorDOFIndices(elementIndex, ldom); size_t i, I; @@ -665,7 +1000,7 @@ namespace ippl { I = global_dofs[i]; - if (this->isDOFOnBoundary(I)) { + if (space.isDOFOnBoundary(I)) { continue; } @@ -673,9 +1008,10 @@ namespace ippl { T contrib = 0; for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) { // Get the global position of the quadrature point - point_t pos = this->ref_element_m.localToGlobal( - this->getElementMeshVertexPoints(this->getElementNDIndex(elementIndex)), - q[k]); + point_t pos = space.ref_element_m.localToGlobal( + space.getElementMeshVertexPoints( + space.getElementNDIndex(elementIndex)), + q[k]); // evaluate the rhs function at this global position point_t interpolatedVal = f(pos); @@ -845,13 +1181,14 @@ namespace ippl { auto coefView = coef.getView(); Kokkos::View outView("reconstructed Func values at points", positions.extent(0)); + const DeviceStruct space = getDeviceMirror(); Kokkos::parallel_for("reconstructToPoints", positions.extent(0), - KOKKOS_CLASS_LAMBDA(size_t i) { + KOKKOS_LAMBDA(size_t i) { // get the current position and for it figure out to which // element it belongs point_t pos = positions<:i:>; - indices_t elemIdx = ((pos - this->origin_m) / domainSize) * gextent; + indices_t elemIdx = ((pos - space.origin_m) / domainSize) * gextent; // next up we have to handle the case of when a position that @@ -868,12 +1205,12 @@ namespace ippl { // get correct indices const Vector vectorIndices = - this->getFEMVectorDOFIndices(elemIdx, ldom); + space.getFEMVectorDOFIndices(elemIdx, ldom); // figure out position inside of the reference element - point_t locPos = pos - (elemIdx * this->hr_m + this->origin_m); - locPos /= this->hr_m; + point_t locPos = pos - (elemIdx * space.hr_m + space.origin_m); + locPos /= space.hr_m; // because of numerical instabilities it might happen then when // a point is on an edge this becomes marginally larger that 1 @@ -889,7 +1226,7 @@ namespace ippl { // basis functions. point_t val(0); for (size_t j = 0; j < numElementDOFs; ++j) { - point_t funcVal = this->evaluateRefElementShapeFunction(j, locPos); + point_t funcVal = space.evaluateRefElementShapeFunction(j, locPos); val += funcVal*coefView(vectorIndices<:j:>); } outView(i) = val; @@ -949,18 +1286,20 @@ namespace ippl { using policy_type = Kokkos::RangePolicy; auto view = u_h.getView(); + const DeviceStruct space = getDeviceMirror(); + auto elementIndicesView = elementIndices; // Loop over elements to compute contributions Kokkos::parallel_reduce("Compute error over elements", policy_type(0, elementIndices.extent(0)), - KOKKOS_CLASS_LAMBDA(size_t index, double& local) { - const size_t elementIndex = elementIndices(index); + KOKKOS_LAMBDA(size_t index, double& local) { + const size_t elementIndex = elementIndicesView(index); const Vector global_dofs = - this->NedelecSpace::getGlobalDOFIndices(elementIndex); + space.getGlobalDOFIndices(elementIndex); const Vector vectorIndices = - this->getFEMVectorDOFIndices(elementIndex, ldom); + space.getFEMVectorDOFIndices(elementIndex, ldom); // contribution of this element to the error @@ -968,9 +1307,9 @@ namespace ippl { for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) { // Evaluate the analystical solution at the global position // of the quadrature point - point_t val_u_sol = u_sol(this->ref_element_m.localToGlobal( - this->getElementMeshVertexPoints(this->getElementNDIndex(elementIndex)), - q[k])); + point_t val_u_sol = u_sol(space.ref_element_m.localToGlobal( + space.getElementMeshVertexPoints(space.getElementNDIndex(elementIndex)), + q[k])); // Here we now reconstruct the solution given the basis // functions. From ee4985bc3146783bb52ca9b25d5a0a837a892ea6 Mon Sep 17 00:00:00 2001 From: Andreas Adelmann Date: Mon, 17 Aug 2026 11:28:14 +0200 Subject: [PATCH 6/9] Make warning-sensitive host paths explicit Resolve the remaining OpenMP and common-host diagnostics by expressing initialization, dimensional specialization, and index-domain conversions directly in the code rather than suppressing compiler warnings. Value-initialize the FEL Charge aggregate so every member, including the electron charge field reported by GCC, has a defined state. Initialize the bunching factor before the optional modulation branch and scope the derived bunching amplitude and electron count to the branches that consume them, eliminating maybe-uninitialized data flow without changing the generated bunch distribution. Mark the pruned complex FFT third-dimension extents as maybe_unused because they are intentionally present for the shared two- and three-dimensional setup but are consumed only by the three-dimensional FFTW calls. Convert the particle scatter policy endpoint to the hash-view extent type before comparison so the bounds check uses one unsigned index domain. Place the full three-dimensional ORB corner-migration test inside the else branch of its if constexpr guard. Non-three-dimensional typed fixtures now compile only the GTEST_SKIP path, preventing unreachable rank-check code from being instantiated while preserving the complete three-dimensional regression. Validated in warning-free clean OpenMP and CUDA sm_80 builds and in the full unit suite at one, two, and four MPI ranks on both backends. All supported configurations pass; the only four-rank exception remains the unrelated Nedelec fixture with a domain that supports at most three partitions. --- demos/fel/MithraBunch.h | 9 ++- src/FFT/Transform/PrunedCC.h | 11 ++- src/Particle/ParticleAttrib.hpp | 4 +- .../Particle/ParticleUpdateNonuniform.cpp | 73 ++++++++++--------- 4 files changed, 52 insertions(+), 45 deletions(-) diff --git a/demos/fel/MithraBunch.h b/demos/fel/MithraBunch.h index bb828326d6..16dfa11069 100644 --- a/demos/fel/MithraBunch.h +++ b/demos/fel/MithraBunch.h @@ -176,14 +176,14 @@ void initializeBunchEllipsoid(BunchInitialize bunchInit, ChargeVector charge; + Charge charge{}; charge.q = bunchInit.cloudCharge_ / Np; FieldVector gb = bunchInit.initialGamma_ * bunchInit.betaVector_; FieldVector r(0.0); FieldVector t(0.0); Double t0; //, g; Double zmin = 1e100; - Double Ne, bF, bFi; + Double bF = 0.0; unsigned int bmi; std::vector randomNumbers; @@ -233,7 +233,7 @@ void initializeBunchEllipsoid(BunchInitialize bunchInit, ChargeVector bunchInit, ChargeVector(gDomFull[0].length()); const long g1 = static_cast(gDomFull[1].length()); - const long g2 = + [[maybe_unused]] const long g2 = (Dim == 3) ? static_cast(gDomFull[Dim == 3 ? 2 : 0].length()) : 1L; const long m0 = static_cast(modes[0]); const long m1 = static_cast(modes[1]); - const long m2 = (Dim == 3) ? static_cast(modes[Dim == 3 ? 2 : 0]) : 1L; + [[maybe_unused]] const long m2 = + (Dim == 3) ? static_cast(modes[Dim == 3 ? 2 : 0]) : 1L; const int lf0 = localFirst[0]; const int lf1 = localFirst[1]; const int lf2 = (Dim == 3) ? localFirst[Dim == 3 ? 2 : 0] : 0; @@ -425,10 +426,12 @@ namespace ippl { (Dim == 3) ? static_cast(owned[Dim == 3 ? 2 : 0].length()) : 1L; const long g0 = static_cast(gDomFull[0].length()); const long g1 = static_cast(gDomFull[1].length()); - const long g2 = (Dim == 3) ? static_cast(gDomFull[Dim == 3 ? 2 : 0].length()) : 1L; + [[maybe_unused]] const long g2 = + (Dim == 3) ? static_cast(gDomFull[Dim == 3 ? 2 : 0].length()) : 1L; const long m0 = static_cast(modes[0]); const long m1 = static_cast(modes[1]); - const long m2 = (Dim == 3) ? static_cast(modes[Dim == 3 ? 2 : 0]) : 1L; + [[maybe_unused]] const long m2 = + (Dim == 3) ? static_cast(modes[Dim == 3 ? 2 : 0]) : 1L; const int lf0 = localFirst[0]; const int lf1 = localFirst[1]; const int lf2 = (Dim == 3) ? localFirst[Dim == 3 ? 2 : 0] : 0; diff --git a/src/Particle/ParticleAttrib.hpp b/src/Particle/ParticleAttrib.hpp index 15cb155d9c..2ab12b0bdb 100644 --- a/src/Particle/ParticleAttrib.hpp +++ b/src/Particle/ParticleAttrib.hpp @@ -157,7 +157,9 @@ namespace ippl { // using policy_type = Kokkos::RangePolicy; const bool useHashView = hash_array.extent(0) > 0; - if (useHashView && (iteration_policy.end() > hash_array.extent(0))) { + const auto policyEnd = static_cast( + iteration_policy.end()); + if (useHashView && (policyEnd > hash_array.extent(0))) { Inform m("scatter"); m << "Hash array was passed to scatter, but size does not match iteration policy." << endl; diff --git a/unit_tests/Particle/ParticleUpdateNonuniform.cpp b/unit_tests/Particle/ParticleUpdateNonuniform.cpp index 331e577c4c..61b3aed1d5 100644 --- a/unit_tests/Particle/ParticleUpdateNonuniform.cpp +++ b/unit_tests/Particle/ParticleUpdateNonuniform.cpp @@ -877,49 +877,50 @@ TYPED_TEST(TestParticleUpdateORB, SuccessiveDisplacementsAcrossOrbBoundaries) { TYPED_TEST(TestParticleUpdateORB, ThreeDCornerMigrationAfterOrb) { if constexpr (TestFixture::Dim != 3) { GTEST_SKIP() << "3-D specific test"; - } - REQUIRE_RANKS(PREF_RANKS); - using T = typename TestFixture::T; + } else { + REQUIRE_RANKS(PREF_RANKS); + using T = typename TestFixture::T; - bool ok = this->orbGaussian(T(0.2)); - if (!ok) - GTEST_SKIP(); - this->rebuildPlayout(); + bool ok = this->orbGaussian(T(0.2)); + if (!ok) + GTEST_SKIP(); + this->rebuildPlayout(); - auto bunch = this->makeBunch(); + auto bunch = this->makeBunch(); - // Seed all particles near the origin - bunch->create(64); - { - auto R_host = bunch->R.getHostMirror(); - auto Q_host = bunch->Q.getHostMirror(); - for (size_t i = 0; i < bunch->getLocalNum(); ++i) { - for (unsigned d = 0; d < 3; d++) - R_host(i)[d] = T(0.02) * this->domain[d]; - Q_host(i) = T(1); + // Seed all particles near the origin + bunch->create(64); + { + auto R_host = bunch->R.getHostMirror(); + auto Q_host = bunch->Q.getHostMirror(); + for (size_t i = 0; i < bunch->getLocalNum(); ++i) { + for (unsigned d = 0; d < 3; d++) + R_host(i)[d] = T(0.02) * this->domain[d]; + Q_host(i) = T(1); + } + Kokkos::deep_copy(bunch->R.getView(), R_host); + Kokkos::deep_copy(bunch->Q.getView(), Q_host); } - Kokkos::deep_copy(bunch->R.getView(), R_host); - Kokkos::deep_copy(bunch->Q.getView(), Q_host); - } - const size_t before = this->totalParticles(*bunch); - bunch->update(); // settle near origin + const size_t before = this->totalParticles(*bunch); + bunch->update(); // settle near origin - // Jump diagonally to the opposite corner (wraps periodically) - { - auto R_host = bunch->R.getHostMirror(); - Kokkos::deep_copy(R_host, bunch->R.getView()); - for (size_t i = 0; i < bunch->getLocalNum(); ++i) - for (unsigned d = 0; d < 3; d++) { - T np = R_host(i)[d] + T(0.98) * this->domain[d]; - R_host(i)[d] = this->periodicWrap(np, this->domain[d]); - } - Kokkos::deep_copy(bunch->R.getView(), R_host); - } + // Jump diagonally to the opposite corner (wraps periodically) + { + auto R_host = bunch->R.getHostMirror(); + Kokkos::deep_copy(R_host, bunch->R.getView()); + for (size_t i = 0; i < bunch->getLocalNum(); ++i) + for (unsigned d = 0; d < 3; d++) { + T np = R_host(i)[d] + T(0.98) * this->domain[d]; + R_host(i)[d] = this->periodicWrap(np, this->domain[d]); + } + Kokkos::deep_copy(bunch->R.getView(), R_host); + } - bunch->update(); - EXPECT_EQ(before, this->totalParticles(*bunch)); - EXPECT_EQ(0u, this->countMisplaced(*bunch)); + bunch->update(); + EXPECT_EQ(before, this->totalParticles(*bunch)); + EXPECT_EQ(0u, this->countMisplaced(*bunch)); + } } // ============================================================ From e6e2387d0b29fc15f967f1ea6ce4780f29971520 Mon Sep 17 00:00:00 2001 From: Andreas Adelmann Date: Mon, 17 Aug 2026 11:28:45 +0200 Subject: [PATCH 7/9] Avoid deprecated Kokkos View conversion paths Remove the Kokkos 5 legacy View diagnostics that were hidden behind the CDash warning display cap. These warnings came from rank-zero scalar deep copies, View remapping during overwrite-only resize operations, direct access to complex View storage, and complex-valued Kokkos subviews. Represent particle destination and overlap counters as one-element rank-one views. ParticleSpatialLayout keeps a persistent host mirror for its counter, while the overlap layout mirrors its final send count explicitly. This preserves the atomic counter behavior without instantiating the deprecated rank-zero View conversion machinery. Use realloc when RegionLayout and BareField replace every element of their allocations, avoiding unnecessary remap kernels and accurately documenting that updateLayout does not preserve field values. Serialize and deserialize scalar values through mdspan element access and local byte-copy temporaries, so complex values no longer require legacy View pointer conversions. Build halo slices with submdspan and construct pack/unpack range policies from the slice extents, keeping complex halo exchange on the supported Kokkos 5 mdspan path. Capture nested overlap interaction lambdas by value, and give the truncated-Green test particle type a uniquely named namespace so CUDA-generated closure visibility matches the captured functor type. Validated by warning-free clean OpenMP and CUDA sm_80 builds, including focused NUFFT and truncated-Green interaction recompiles. Full unit tests pass 41/41 at one and two ranks on both backends and 40/41 at four ranks, where only the pre-existing three-partition Nedelec fixture is unsupported. --- src/Communicate/Archive.hpp | 9 +++-- src/Field/BareField.hpp | 10 ++++-- src/Field/HaloCells.hpp | 36 ++++++++++++++----- src/Particle/ParticleSpatialLayout.h | 6 ++-- src/Particle/ParticleSpatialLayout.hpp | 11 +++--- src/Particle/ParticleSpatialOverlapLayout.hpp | 22 ++++++------ src/Region/RegionLayout.hpp | 6 ++-- .../TestTruncatedGreenParticleInteraction.cpp | 6 ++-- 8 files changed, 70 insertions(+), 36 deletions(-) diff --git a/src/Communicate/Archive.hpp b/src/Communicate/Archive.hpp index a719a38af4..4f14609015 100644 --- a/src/Communicate/Archive.hpp +++ b/src/Communicate/Archive.hpp @@ -198,9 +198,11 @@ namespace ippl { size_t size = sizeof(T); auto base = bufferData(); auto writepos = writepos_m; + auto viewSpan = view.to_mdspan(); Kokkos::parallel_for( "Archive::serialize()", policy_type(0, nsends), KOKKOS_LAMBDA(const size_type i) { - const char* src = reinterpret_cast(view.data() + i); + const T value = viewSpan(i); + const char* src = reinterpret_cast(&value); char* dst = base + i * size + writepos; copyBytes(dst, src, size); }); @@ -303,11 +305,14 @@ namespace ippl { } auto base = bufferData(); auto readpos = readpos_m; + auto viewSpan = view.to_mdspan(); Kokkos::parallel_for( "Archive::deserialize()", policy_type(0, nrecvs), KOKKOS_LAMBDA(const size_type i) { const char* src = base + i * size + readpos; - char* dst = reinterpret_cast(view.data() + i); + T value{}; + char* dst = reinterpret_cast(&value); copyBytes(dst, src, size); + viewSpan(i) = value; }); // Wait for deserialization kernel to complete // (as with serialization kernels) diff --git a/src/Field/BareField.hpp b/src/Field/BareField.hpp index 14e7b1ec6b..327b138db5 100644 --- a/src/Field/BareField.hpp +++ b/src/Field/BareField.hpp @@ -137,10 +137,14 @@ namespace ippl { void BareField::setup() { owned_m = layout_m->getLocalNDIndex(); - auto resize = [&](const std::index_sequence&) { - this->resize((owned_m[Idx].length() + 2 * nghost_m)...); + // A new layout establishes new local storage; updateLayout does not + // redistribute or preserve field values. Allocate directly instead of + // asking Kokkos::resize to remap the previous allocation through its + // deprecated legacy subview path. + auto reallocate = [&](const std::index_sequence&) { + Kokkos::realloc(dview_m, (owned_m[Idx].length() + 2 * nghost_m)...); }; - resize(std::make_index_sequence{}); + reallocate(std::make_index_sequence{}); } template diff --git a/src/Field/HaloCells.hpp b/src/Field/HaloCells.hpp index 571f12ca71..3b00cdb49a 100644 --- a/src/Field/HaloCells.hpp +++ b/src/Field/HaloCells.hpp @@ -254,13 +254,21 @@ namespace ippl { Kokkos::realloc(buffer, size * overalloc); } - using index_array_type = - typename RangePolicy::index_array_type; + using exec_space = typename view_type::execution_space; + using range_type = RangePolicy; + using index_type = typename range_type::index_type; + using index_array_type = typename range_type::index_array_type; using buffer_view_type = typename databuffer_type::view_type; using functor_type = HaloPackFunctor; - ippl::parallel_for("HaloCells::pack()", getRangePolicy(subview), - functor_type{subview, buffer}); + + Kokkos::Array begin{}; + Kokkos::Array end{}; + for (unsigned d = 0; d < Dim; ++d) { + end[d] = static_cast(subview.extent(d)); + } + auto policy = createRangePolicy(begin, end); + ippl::parallel_for("HaloCells::pack()", policy, functor_type{subview, buffer}); Kokkos::fence(); } @@ -275,11 +283,20 @@ namespace ippl { // https://stackoverflow.com/questions/3735398/operator-as-template-parameter Op op; - using index_array_type = - typename RangePolicy::index_array_type; + using exec_space = typename view_type::execution_space; + using range_type = RangePolicy; + using index_type = typename range_type::index_type; + using index_array_type = typename range_type::index_array_type; using functor_type = HaloUnpackFunctor; - ippl::parallel_for("HaloCells::unpack()", getRangePolicy(subview), + + Kokkos::Array begin{}; + Kokkos::Array end{}; + for (unsigned d = 0; d < Dim; ++d) { + end[d] = static_cast(subview.extent(d)); + } + auto policy = createRangePolicy(begin, end); + ippl::parallel_for("HaloCells::unpack()", policy, functor_type{subview, buffer, op}); Kokkos::fence(); } @@ -288,8 +305,9 @@ namespace ippl { auto HaloCells::makeSubview(const view_type& view, const bound_type& intersect) { auto makeSub = [&](const std::index_sequence&) { - return Kokkos::subview(view, - Kokkos::make_pair(intersect.lo[Idx], intersect.hi[Idx])...); + return Kokkos::submdspan( + view.to_mdspan(), + Kokkos::make_pair(intersect.lo[Idx], intersect.hi[Idx])...); }; return makeSub(std::make_index_sequence{}); } diff --git a/src/Particle/ParticleSpatialLayout.h b/src/Particle/ParticleSpatialLayout.h index aca4106c36..cfab632862 100644 --- a/src/Particle/ParticleSpatialLayout.h +++ b/src/Particle/ParticleSpatialLayout.h @@ -156,8 +156,10 @@ namespace ippl { locate_type destRanks_d_; // [nRanks] (compacted list) bool_type leaving_d_; // [capacity >= max nLocal seen] mask - // Single scalar on device to count destinations - Kokkos::View nDest_d_; + // One-element views used to count destinations. A rank-1 view avoids + // Kokkos' deprecated legacy conversion path for rank-0 View deep copies. + Kokkos::View nDest_d_; + Kokkos::View nDest_h_; // Neigbour cache locate_type neighbors_d_; // [neighborSize] cached device neighbors list diff --git a/src/Particle/ParticleSpatialLayout.hpp b/src/Particle/ParticleSpatialLayout.hpp index 6106161917..5042a8dc7b 100644 --- a/src/Particle/ParticleSpatialLayout.hpp +++ b/src/Particle/ParticleSpatialLayout.hpp @@ -147,8 +147,8 @@ namespace ippl { const size_type nInvalid = locateParticlesPacked(pc); // Copy metadata to host - size_type nDest = 0; - Kokkos::deep_copy(position_execution_space{}, nDest, nDest_d_); + Kokkos::deep_copy(nDest_h_, nDest_d_); + const size_type nDest = nDest_h_(0); // destRanks prefix if (nDest > 0) { @@ -457,7 +457,7 @@ namespace ippl { if ((size_type)r == myRank) return; if (rankSendCount_d(r) > 0) { - const size_type idx = Kokkos::atomic_fetch_add(&nDest_d(), size_type(1)); + const size_type idx = Kokkos::atomic_fetch_add(&nDest_d(0), size_type(1)); destRanks_d(idx) = static_cast(r); } }); @@ -474,8 +474,9 @@ namespace ippl { Kokkos::realloc(destRanks_d_, nRanks); Kokkos::realloc(recvCounts_d_, nRanks); - // scalar counter - nDest_d_ = Kokkos::View("nDest_d"); + // One-element device counter and its persistent host mirror + Kokkos::realloc(nDest_d_, 1); + Kokkos::realloc(nDest_h_, 1); // Host mirrors Kokkos::realloc(rankSendCount_h_, nRanks); diff --git a/src/Particle/ParticleSpatialOverlapLayout.hpp b/src/Particle/ParticleSpatialOverlapLayout.hpp index 42f1a26191..f8c4771080 100644 --- a/src/Particle/ParticleSpatialOverlapLayout.hpp +++ b/src/Particle/ParticleSpatialOverlapLayout.hpp @@ -530,13 +530,13 @@ namespace ippl { Kokkos::fence(); // Step 2. Fill remaining ranks - Kokkos::View counter("counter"); - Kokkos::deep_copy(counter, 0); + Kokkos::View counter("counter", 1); + Kokkos::deep_copy(counter, size_type{0}); Kokkos::fence(); Kokkos::parallel_for( "fill_remaining", policy_type(0, total_ranks), KOKKOS_LAMBDA(const size_t& i) { if (is_remaining(i)) { - const size_type idx = Kokkos::atomic_fetch_inc(&counter()); + const size_type idx = Kokkos::atomic_fetch_inc(&counter(0)); nonNeighborRanks(idx) = i; } }); @@ -764,20 +764,20 @@ namespace ippl { Kokkos::fence(); /* compute the ranks to send to and the number of ranks to send to*/ - Kokkos::View rankSends( - "Number of Ranks we need to send to"); + Kokkos::View rankSends( + "Number of Ranks we need to send to", 1); Kokkos::parallel_for( "Calculate sends", policy_type(0, nSends_dview.extent(0)), KOKKOS_LAMBDA(const size_t rank) { if (nSends_dview(rank) != 0) { - size_type index = Kokkos::atomic_fetch_inc(&rankSends()); + size_type index = Kokkos::atomic_fetch_inc(&rankSends(0)); sends_dview(index) = rank; } }); Kokkos::fence(); - size_type temp; - Kokkos::deep_copy(temp, rankSends); - Kokkos::fence(); + const auto rankSendsHost = + Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), rankSends); + const size_type temp = rankSendsHost(0); return {invalidCount, temp}; } @@ -1135,7 +1135,7 @@ namespace ippl { /* iterate over all cell neighbors */ Kokkos::parallel_for( - Kokkos::TeamThreadRange(team, numCellNeighbors), [&](const size_t& n) { + Kokkos::TeamThreadRange(team, numCellNeighbors), [=](const size_t& n) { const auto neighborCellIdx = cellNeighbors[n]; const auto neighborCellParticleOffset = cellStartingIdx(neighborCellIdx); const auto numNeighborCellParticles = cellParticleCount(neighborCellIdx); @@ -1145,7 +1145,7 @@ namespace ippl { */ Kokkos::parallel_for(Kokkos::ThreadVectorMDRange, team_t>( team, numCellParticles, numNeighborCellParticles), - [&](const size_t& i, const size_t& j) { + [=](const size_t& i, const size_t& j) { const auto particleIdx = cellParticleOffset + i; const auto neighborIdx = neighborCellParticleOffset + j; diff --git a/src/Region/RegionLayout.hpp b/src/Region/RegionLayout.hpp index d074f6825e..fd6ee7fd33 100644 --- a/src/Region/RegionLayout.hpp +++ b/src/Region/RegionLayout.hpp @@ -103,8 +103,10 @@ namespace ippl { using domain_type = typename FieldLayout::host_mirror_type; const domain_type& ldomains = fl.getHostLocalDomains(); - Kokkos::resize(hLocalRegions_m, ldomains.size()); - Kokkos::resize(dLocalRegions_m, ldomains.size()); + // Every entry is replaced below, so preserving the previous allocation contents via + // Kokkos::resize only performs an unnecessary View remap. Reallocate instead. + Kokkos::realloc(hLocalRegions_m, ldomains.size()); + Kokkos::realloc(dLocalRegions_m, ldomains.size()); using size_type = typename domain_type::size_type; for (size_type i = 0; i < ldomains.size(); ++i) { diff --git a/test/particle/TestTruncatedGreenParticleInteraction.cpp b/test/particle/TestTruncatedGreenParticleInteraction.cpp index f41d6984c7..76dc6149c1 100644 --- a/test/particle/TestTruncatedGreenParticleInteraction.cpp +++ b/test/particle/TestTruncatedGreenParticleInteraction.cpp @@ -21,7 +21,7 @@ #include "Particle/ParticleSpatialOverlapLayout.h" -namespace { +namespace test_truncated_green_particle_interaction { constexpr unsigned Dim = 3; using Scalar_t = double; using Mesh_t = ippl::UniformCartesian; @@ -110,7 +110,9 @@ namespace { ippl::Comm->allreduce(normLocal, norm, 1, std::plus()); return norm; } -} // namespace +} // namespace test_truncated_green_particle_interaction + +using namespace test_truncated_green_particle_interaction; int main(int argc, char* argv[]) { ippl::initialize(argc, argv); From e6b99634683a1b25420ea7d1171591db730c2130 Mon Sep 17 00:00:00 2001 From: Andreas Adelmann Date: Mon, 17 Aug 2026 21:36:51 +0200 Subject: [PATCH 8/9] Removed the obsolete outer Ne and bFi declarations from demos/fel/MithraBunch.h --- demos/fel/MithraBunch.h | 2 -- 1 file changed, 2 deletions(-) diff --git a/demos/fel/MithraBunch.h b/demos/fel/MithraBunch.h index ef3b84191d..e99513a862 100644 --- a/demos/fel/MithraBunch.h +++ b/demos/fel/MithraBunch.h @@ -183,9 +183,7 @@ void initializeBunchEllipsoid(BunchInitialize bunchInit, ChargeVector t(0.0); Double t0; //, g; Double zmin = 1e100; - Double Ne; Double bF = bunchInit.bF_; - Double bFi; unsigned int bmi; std::vector randomNumbers; From 4dc1f68c932dc01715aee564287beb544545d120 Mon Sep 17 00:00:00 2001 From: Andreas Adelmann Date: Tue, 1 Sep 2026 17:58:08 +0200 Subject: [PATCH 9/9] Address comments from Sonali --- src/FEM/LagrangeSpace.h | 20 +++++++++----------- src/FEM/NedelecSpace.h | 5 +---- src/Field/BareField.hpp | 6 ++---- src/Particle/ParticleSpatialLayout.h | 5 +++-- 4 files changed, 15 insertions(+), 21 deletions(-) diff --git a/src/FEM/LagrangeSpace.h b/src/FEM/LagrangeSpace.h index 50c231da23..75ac3d18cd 100644 --- a/src/FEM/LagrangeSpace.h +++ b/src/FEM/LagrangeSpace.h @@ -319,13 +319,14 @@ namespace ippl { * host creates the structure before launching a kernel; capture-by-value places it in * the Kokkos kernel closure, which Kokkos makes available in the selected execution * space. getDeviceMirror() itself performs no field allocation or field-data copy. - * - It intentionally excludes owning or host-oriented state such as `resultField`, field - * and layout objects, MPI decomposition objects, and the `elementIndices` allocation. - * A kernel captures only the Kokkos views and scalar values needed for that invocation, - * separately from this geometry/indexing snapshot. + * - It does not contain field objects, field layouts, MPI decomposition state, or other + * host-owned infrastructure. Before each kernel launch, the host extracts the required + * device-accessible Kokkos views—such as the `resultField` data view and + * `elementIndices`—and captures those view handles separately from this geometry and + * indexing snapshot. * - It removes the need to capture `this`. Consequently, constructing and destroying the * device closure never requires the host-only lifetime operations of `LagrangeSpace` or - * `Field`, which is the source of the GH200 CUDA diagnostic addressed by this design. + * `Field`, which would cause issues when using GPUs. * * "Lightweight" therefore means that this is neither a second owning finite-element space * nor an automatically synchronized copy of one. It is a small, non-owning value snapshot @@ -340,17 +341,14 @@ namespace ippl { * be created after changing any mirrored mesh or reference-element state on the host. * During a kernel invocation the snapshot is read-only and may be shared by all threads. * - * Element and vertex indices use the same flattened ordering as `LagrangeSpace`: - * dimension zero varies fastest. An element's N-dimensional index identifies its lower - * mesh vertex. */ struct DeviceStruct { // members we need to copy for the following functions: // works since numElementDOFs in LagrangeSpace is static constexpr - static constexpr unsigned numElementDOFs = LagrangeSpace::numElementDOFs; + static constexpr unsigned numElementDOFs = LagrangeSpace::numElementDOFs; static constexpr unsigned numElementVertices = LagrangeSpace::numElementVertices; - using indices_list_t = Vector; - using vertex_points_t = Vector; + using indices_list_t = Vector; + using vertex_points_t = Vector; Vector nr_m; ///< Number of mesh vertices in each dimension. Vector hr_m; ///< Uniform mesh spacing in each dimension. diff --git a/src/FEM/NedelecSpace.h b/src/FEM/NedelecSpace.h index af455733c4..552f8bd50d 100644 --- a/src/FEM/NedelecSpace.h +++ b/src/FEM/NedelecSpace.h @@ -421,12 +421,9 @@ namespace ippl { * future kernel needs more Nedelec functionality, add only the smallest device-copyable * state and algorithm here; do not restore capture of the parent `NedelecSpace`. * - * Element flattening uses dimension zero as the fastest-varying dimension. Nedelec DOFs - * retain the edge ordering used by the host class, so global and ghosted FEMVector index - * calculations remain identical on host and device. */ struct DeviceStruct { - static constexpr unsigned numElementDOFs = NedelecSpace::numElementDOFs; + static constexpr unsigned numElementDOFs = NedelecSpace::numElementDOFs; static constexpr unsigned numElementVertices = NedelecSpace::numElementVertices; using indices_list_t = Vector; diff --git a/src/Field/BareField.hpp b/src/Field/BareField.hpp index 327b138db5..6622813803 100644 --- a/src/Field/BareField.hpp +++ b/src/Field/BareField.hpp @@ -137,10 +137,8 @@ namespace ippl { void BareField::setup() { owned_m = layout_m->getLocalNDIndex(); - // A new layout establishes new local storage; updateLayout does not - // redistribute or preserve field values. Allocate directly instead of - // asking Kokkos::resize to remap the previous allocation through its - // deprecated legacy subview path. + // Reallocate field to correspond to sizes of new updated layout. + // Does not preserve field data. auto reallocate = [&](const std::index_sequence&) { Kokkos::realloc(dview_m, (owned_m[Idx].length() + 2 * nghost_m)...); }; diff --git a/src/Particle/ParticleSpatialLayout.h b/src/Particle/ParticleSpatialLayout.h index cfab632862..c0fd93cd8c 100644 --- a/src/Particle/ParticleSpatialLayout.h +++ b/src/Particle/ParticleSpatialLayout.h @@ -156,8 +156,9 @@ namespace ippl { locate_type destRanks_d_; // [nRanks] (compacted list) bool_type leaving_d_; // [capacity >= max nLocal seen] mask - // One-element views used to count destinations. A rank-1 view avoids - // Kokkos' deprecated legacy conversion path for rank-0 View deep copies. + // Note: A rank-1 view avoids Kokkos' deprecated legacy conversion + // for rank-0 View deep copies. + // Rank-1 views used to count destinations. Kokkos::View nDest_d_; Kokkos::View nDest_h_;