Line data Source code
1 : #ifndef IPPL_FEMINTERPOLATE_H
2 : #define IPPL_FEMINTERPOLATE_H
3 :
4 :
5 : namespace ippl {
6 :
7 : /**
8 : * @brief Mapping from global position to element ND index and
9 : * reference coordinates (xi ∈ [0,1)^Dim) on a UniformCartesian mesh.
10 : *
11 : * Assumes the input x is strictly inside the computational domain so that
12 : * for each dimension d: 0 ≤ (x[d]-origin[d])/h[d] < nr[d]-1.
13 : */
14 : template <typename T, unsigned Dim, class Mesh>
15 : KOKKOS_INLINE_FUNCTION void
16 24024 : locate_element_nd_and_xi(const Mesh& mesh,
17 : const ippl::Vector<T,Dim>& x,
18 : ippl::Vector<size_t,Dim>& e_nd,
19 : ippl::Vector<T,Dim>& xi) {
20 :
21 24024 : const auto nr = mesh.getGridsize(); // vertices per axis
22 24024 : const auto h = mesh.getMeshSpacing();
23 24024 : const auto org = mesh.getOrigin();
24 :
25 :
26 72072 : for (unsigned d = 0; d < Dim; ++d) {
27 48048 : const T s = (x[d] - org[d]) / h[d]; // To cell units
28 48048 : const size_t e = static_cast<size_t>(std::floor(s));
29 48048 : e_nd[d] = e;
30 48048 : xi[d] = s - static_cast<T>(e);
31 : }
32 24024 : }
33 :
34 : /**
35 : * @brief Assemble a P1 FEM load vector (RHS) from particle attributes.
36 : *
37 : * For each particle position x, locate the owning element (ND index e_nd) and
38 : * reference coordinate xi. Deposit the particle attribute value into the
39 : * element's nodal DOFs using P1 Lagrange shape functions evaluated at xi.
40 : *
41 : * @tparam AttribIn Particle attribute type with getView()(p) -> scalar
42 : * @tparam Field ippl::Field with rank=Dim nodal coefficients (RHS)
43 : * @tparam PosAttrib Particle position attribute with getView()(p) -> Vector<T,Dim>
44 : * @tparam Space Lagrange space providing element/DOF/topology queries
45 : * @tparam policy_type Kokkos execution policy (defaults to Field::execution_space)
46 : */
47 : template <typename AttribIn, typename Field, typename PosAttrib, typename Space,
48 : typename policy_type = Kokkos::RangePolicy<typename Field::execution_space>>
49 12 : inline void assemble_rhs_from_particles(const AttribIn& attrib, Field& f,
50 : const PosAttrib& pp, const Space& space,
51 : policy_type iteration_policy)
52 : {
53 12 : constexpr unsigned Dim = Field::dim;
54 : using T = typename Field::value_type;
55 : using view_type = typename Field::view_type;
56 : using mesh_type = typename Field::Mesh_t;
57 :
58 12 : static IpplTimings::TimerRef t = IpplTimings::getTimer("assemble_rhs_from_particles(P1)");
59 :
60 12 : IpplTimings::startTimer(t);
61 :
62 12 : view_type view = f.getView();
63 :
64 : // Mesh / layout (for locating + indexing into the field view)
65 12 : const mesh_type& mesh = f.get_mesh();
66 :
67 12 : const ippl::FieldLayout<Dim>& layout = f.getLayout();
68 12 : const ippl::NDIndex<Dim>& lDom = layout.getLocalNDIndex();
69 12 : const int nghost = f.getNghost();
70 :
71 : // Particle attribute/device views
72 12 : auto d_attr = attrib.getView(); // scalar weight per particle (e.g. charge)
73 12 : auto d_pos = pp.getView(); // positions (Vector<T,Dim>) per particle
74 :
75 12 : Kokkos::parallel_for("assemble_rhs_from_particles_P1", iteration_policy,
76 12024 : KOKKOS_LAMBDA(const size_t p) {
77 12000 : const Vector<T,Dim> x = d_pos(p);
78 12000 : const T val = d_attr(p);
79 :
80 12000 : Vector<size_t,Dim> e_nd;
81 12000 : Vector<T,Dim> xi;
82 12000 : locate_element_nd_and_xi<T,Dim>(mesh, x, e_nd, xi);
83 : // Convert to the element's linear index
84 12000 : const size_t e_lin = space.getElementIndex(e_nd);
85 :
86 : // DOFs for this element
87 12000 : const auto dofs = space.getGlobalDOFIndices(e_lin);
88 :
89 : // Deposit into each vertex/DOF
90 68000 : for (size_t a = 0; a < dofs.dim; ++a) {
91 56000 : const size_t local = space.getLocalDOFIndex(e_lin, dofs[a]);
92 56000 : const T w = space.evaluateRefElementShapeFunction(local, xi);
93 :
94 56000 : const auto v_nd = space.getMeshVertexNDIndex(dofs[a]); // ND coords (global, vertex-centered)
95 56000 : ippl::Vector<size_t,Dim> I; // indices into view
96 :
97 192000 : for (unsigned d = 0; d < Dim; ++d) {
98 136000 : I[d] = static_cast<size_t>(v_nd[d] - lDom.first()[d] + nghost);
99 : }
100 : if constexpr (Dim == 1) {
101 16000 : Kokkos::atomic_add(&view(I[0]), val * w);
102 : } else if constexpr (Dim == 2) {
103 32000 : Kokkos::atomic_add(&view(I[0], I[1]), val * w);
104 : } else if constexpr (Dim == 3) {
105 64000 : Kokkos::atomic_add(&view(I[0], I[1], I[2]), val * w);
106 : }
107 : }
108 12000 : }
109 : );
110 :
111 12 : static IpplTimings::TimerRef accumulateHaloTimer = IpplTimings::getTimer("accumulateHalo");
112 12 : IpplTimings::startTimer(accumulateHaloTimer);
113 12 : f.accumulateHalo();
114 12 : IpplTimings::stopTimer(accumulateHaloTimer);
115 :
116 12 : }
117 :
118 : /**
119 : * @brief Interpolate a P1 FEM field to particle positions.
120 : *
121 : * For each particle position x, locate the owning element (ND index e_nd) and
122 : * reference coordinate xi. Evaluate P1 Lagrange shape functions at xi to
123 : * combine nodal coefficients and write u(x) to the particle attribute.
124 : *
125 : * @tparam AttribOut Particle attribute type with getView()(p) -> scalar
126 : * @tparam Field ippl::Field with rank=Dim nodal coefficients
127 : * @tparam PosAttrib Particle position attribute with getView()(p) -> Vector<T,Dim>
128 : * @tparam Space Lagrange space providing element/DOF/topology queries
129 : * @tparam policy_type Kokkos execution policy (defaults to Field::execution_space)
130 : */
131 : template <typename AttribOut, typename Field, typename PosAttrib, typename Space,
132 : typename policy_type = Kokkos::RangePolicy<typename Field::execution_space>>
133 12 : inline void interpolate_to_diracs(AttribOut& attrib_out,
134 : const Field& coeffs,
135 : const PosAttrib& pp,
136 : const Space& space,
137 : policy_type iteration_policy)
138 : {
139 12 : constexpr unsigned Dim = Field::dim;
140 : using T = typename AttribOut::value_type;
141 : using field_value_type = typename Field::value_type;
142 : using view_type = typename Field::view_type;
143 : using mesh_type = typename Field::Mesh_t;
144 :
145 12 : static IpplTimings::TimerRef timer =
146 12 : IpplTimings::getTimer("interpolate_field_to_particles(P1)");
147 12 : IpplTimings::startTimer(timer);
148 :
149 12 : view_type view = coeffs.getView();
150 12 : const mesh_type& M = coeffs.get_mesh();
151 :
152 :
153 12 : const ippl::FieldLayout<Dim>& layout = coeffs.getLayout();
154 12 : const ippl::NDIndex<Dim>& lDom = layout.getLocalNDIndex();
155 12 : const int nghost = coeffs.getNghost();
156 :
157 : // Particle device views
158 12 : auto d_pos = pp.getView();
159 12 : auto d_out = attrib_out.getView();
160 :
161 : // Small device helper to read a rank-Dim view at ND indices
162 56000 : auto read_view = KOKKOS_LAMBDA(const view_type& v,
163 : const ippl::Vector<size_t,Dim>& I) -> field_value_type {
164 :
165 16000 : if constexpr (Dim == 1) return v(I[0]);
166 32000 : if constexpr (Dim == 2) return v(I[0], I[1]);
167 64000 : if constexpr (Dim == 3) return v(I[0], I[1], I[2]);
168 : };
169 :
170 :
171 12 : Kokkos::parallel_for("interpolate_to_diracs_P1", iteration_policy,
172 12024 : KOKKOS_LAMBDA(const size_t p) {
173 :
174 18000 : const Vector<T,Dim> x = d_pos(p);
175 :
176 12000 : ippl::Vector<size_t,Dim> e_nd;
177 12000 : ippl::Vector<T,Dim> xi;
178 12000 : locate_element_nd_and_xi<T,Dim>(M, x, e_nd, xi);
179 12000 : const size_t e_lin = space.getElementIndex(e_nd);
180 :
181 12000 : const auto dofs = space.getGlobalDOFIndices(e_lin);
182 :
183 12000 : field_value_type up = field_value_type(0);
184 :
185 68000 : for (size_t a = 0; a < dofs.dim; ++a) {
186 56000 : const size_t local = space.getLocalDOFIndex(e_lin, dofs[a]);
187 56000 : const field_value_type w = space.evaluateRefElementShapeFunction(local, xi);
188 :
189 56000 : const auto v_nd = space.getMeshVertexNDIndex(dofs[a]);
190 56000 : ippl::Vector<size_t,Dim> I;
191 192000 : for (unsigned d = 0; d < Dim; ++d) {
192 136000 : I[d] = static_cast<size_t>(v_nd[d] - lDom.first()[d] + nghost);
193 : }
194 :
195 56000 : up += read_view(view, I) * w;
196 : }
197 12000 : d_out(p) = static_cast<T>(up);
198 12000 : });
199 12 : }
200 :
201 : } // namespace ippl
202 : #endif
|
