diff --git a/demos/fel/MithraBunch.h b/demos/fel/MithraBunch.h index 626eb0a121..e99513a862 100644 --- a/demos/fel/MithraBunch.h +++ b/demos/fel/MithraBunch.h @@ -176,16 +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; Double bF = bunchInit.bF_; - Double bFi; unsigned int bmi; std::vector randomNumbers; @@ -235,7 +233,7 @@ void initializeBunchEllipsoid(BunchInitialize bunchInit, ChargeVector bunchInit, ChargeVector(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/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/FEM/LagrangeSpace.h b/src/FEM/LagrangeSpace.h index 194368d158..8ee9335b48 100644 --- a/src/FEM/LagrangeSpace.h +++ b/src/FEM/LagrangeSpace.h @@ -287,27 +287,192 @@ 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 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 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 + * 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. + * + */ 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 numElementDOFs = LagrangeSpace::numElementDOFs; + 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 cbcfe5991b..2dc4ff0ddb 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 diff --git a/src/FEM/NedelecSpace.h b/src/FEM/NedelecSpace.h index 14a28a6cb5..7d46368bb0 100644 --- a/src/FEM/NedelecSpace.h +++ b/src/FEM/NedelecSpace.h @@ -392,6 +392,149 @@ 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`. + * + */ + 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 3f9e1787fb..46c76905cc 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. diff --git a/src/FFT/Transform/PrunedCC.h b/src/FFT/Transform/PrunedCC.h index 5cf4655f41..35eef3e33c 100644 --- a/src/FFT/Transform/PrunedCC.h +++ b/src/FFT/Transform/PrunedCC.h @@ -297,11 +297,12 @@ namespace ippl { const long g0 = static_cast(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; diff --git a/src/Field/BareField.hpp b/src/Field/BareField.hpp index 14e7b1ec6b..6622813803 100644 --- a/src/Field/BareField.hpp +++ b/src/Field/BareField.hpp @@ -137,10 +137,12 @@ 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)...); + // 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)...); }; - resize(std::make_index_sequence{}); + reallocate(std::make_index_sequence{}); } template 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/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/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.h b/src/Particle/ParticleSpatialLayout.h index aca4106c36..c0fd93cd8c 100644 --- a/src/Particle/ParticleSpatialLayout.h +++ b/src/Particle/ParticleSpatialLayout.h @@ -156,8 +156,11 @@ 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_; + // 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_; // Neigbour cache locate_type neighbors_d_; // [neighborSize] cached device neighbors list diff --git a/src/Particle/ParticleSpatialLayout.hpp b/src/Particle/ParticleSpatialLayout.hpp index e0590974c6..5042a8dc7b 100644 --- a/src/Particle/ParticleSpatialLayout.hpp +++ b/src/Particle/ParticleSpatialLayout.hpp @@ -147,15 +147,17 @@ 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) { 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)); } @@ -454,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); } }); @@ -471,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); @@ -538,8 +542,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; } 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/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); } /*! 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); 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)); + } } // ============================================================