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 : template<class View, class IVec, std::size_t... Is>
36 : KOKKOS_INLINE_FUNCTION
37 56000 : auto view_ptr_impl(View& v, const IVec& I, std::index_sequence<Is...>)
38 : -> decltype(&v(I[Is]...)) {
39 112000 : return &v(I[Is]...);
40 : }
41 :
42 : template<int D, class View, class IVec>
43 : KOKKOS_INLINE_FUNCTION
44 56000 : auto view_ptr(View& v, const IVec& I)
45 : -> decltype(view_ptr_impl(v, I, std::make_index_sequence<D>{})) {
46 56000 : return view_ptr_impl(v, I, std::make_index_sequence<D>{});
47 : }
48 :
49 : /**
50 : * @brief Assemble a P1 FEM load vector (RHS) from particle attributes.
51 : *
52 : * For each particle position x, locate the owning element (ND index e_nd) and
53 : * reference coordinate xi. Deposit the particle attribute value into the
54 : * element's nodal DOFs using P1 Lagrange shape functions evaluated at xi.
55 : *
56 : * @tparam AttribIn Particle attribute type with getView()(p) -> scalar
57 : * @tparam Field ippl::Field with rank=Dim nodal coefficients (RHS)
58 : * @tparam PosAttrib Particle position attribute with getView()(p) -> Vector<T,Dim>
59 : * @tparam Space Lagrange space providing element/DOF/topology queries
60 : * @tparam policy_type Kokkos execution policy (defaults to Field::execution_space)
61 : */
62 : template <typename AttribIn, typename Field, typename PosAttrib, typename Space,
63 : typename policy_type = Kokkos::RangePolicy<typename Field::execution_space>>
64 12 : inline void assemble_rhs_from_particles(const AttribIn& attrib, Field& f,
65 : const PosAttrib& pp, const Space& space,
66 : policy_type iteration_policy)
67 : {
68 12 : constexpr unsigned Dim = Field::dim;
69 : using T = typename Field::value_type;
70 : using view_type = typename Field::view_type;
71 : using mesh_type = typename Field::Mesh_t;
72 :
73 12 : static IpplTimings::TimerRef t = IpplTimings::getTimer("assemble_rhs_from_particles(P1)");
74 :
75 12 : IpplTimings::startTimer(t);
76 :
77 12 : view_type view = f.getView();
78 :
79 : // Mesh / layout (for locating + indexing into the field view)
80 12 : const mesh_type& mesh = f.get_mesh();
81 :
82 12 : const ippl::FieldLayout<Dim>& layout = f.getLayout();
83 12 : const ippl::NDIndex<Dim>& lDom = layout.getLocalNDIndex();
84 12 : const int nghost = f.getNghost();
85 :
86 : // Particle attribute/device views
87 12 : auto d_attr = attrib.getView(); // scalar weight per particle (e.g. charge)
88 12 : auto d_pos = pp.getView(); // positions (Vector<T,Dim>) per particle
89 :
90 12 : Kokkos::parallel_for("assemble_rhs_from_particles_P1", iteration_policy,
91 12024 : KOKKOS_LAMBDA(const size_t p) {
92 12000 : const Vector<T,Dim> x = d_pos(p);
93 12000 : const T val = d_attr(p);
94 :
95 12000 : Vector<size_t,Dim> e_nd;
96 12000 : Vector<T,Dim> xi;
97 12000 : locate_element_nd_and_xi<T,Dim>(mesh, x, e_nd, xi);
98 : // Convert to the element's linear index
99 12000 : const size_t e_lin = space.getElementIndex(e_nd);
100 :
101 : // DOFs for this element
102 12000 : const auto dofs = space.getGlobalDOFIndices(e_lin);
103 :
104 : // Deposit into each vertex/DOF
105 68000 : for (size_t a = 0; a < dofs.dim; ++a) {
106 56000 : const size_t local = space.getLocalDOFIndex(e_lin, dofs[a]);
107 56000 : const T w = space.evaluateRefElementShapeFunction(local, xi);
108 :
109 56000 : const auto v_nd = space.getMeshVertexNDIndex(dofs[a]); // ND coords (global, vertex-centered)
110 56000 : ippl::Vector<size_t,Dim> I; // indices into view
111 :
112 192000 : for (unsigned d = 0; d < Dim; ++d) {
113 136000 : I[d] = static_cast<size_t>(v_nd[d] - lDom.first()[d] + nghost);
114 : }
115 56000 : Kokkos::atomic_add(view_ptr<Dim>(view, I), val * w);
116 : }
117 12000 : }
118 : );
119 :
120 12 : static IpplTimings::TimerRef accumulateHaloTimer = IpplTimings::getTimer("accumulateHalo");
121 12 : IpplTimings::startTimer(accumulateHaloTimer);
122 12 : f.accumulateHalo();
123 12 : IpplTimings::stopTimer(accumulateHaloTimer);
124 :
125 12 : }
126 :
127 : template<class View, class IVec, std::size_t... Is>
128 : KOKKOS_INLINE_FUNCTION
129 56000 : decltype(auto) view_ref_impl(View& v, const IVec& I, std::index_sequence<Is...>) {
130 112000 : return v(I[Is]...);
131 : }
132 :
133 : template<int D, class View, class IVec>
134 : KOKKOS_INLINE_FUNCTION
135 56000 : decltype(auto) view_ref(View& v, const IVec& I) {
136 56000 : return view_ref_impl(v, I, std::make_index_sequence<D>{});
137 : }
138 :
139 : /**
140 : * @brief Interpolate a P1 FEM field to particle positions.
141 : *
142 : * For each particle position x, locate the owning element (ND index e_nd) and
143 : * reference coordinate xi. Evaluate P1 Lagrange shape functions at xi to
144 : * combine nodal coefficients and write u(x) to the particle attribute.
145 : *
146 : * @tparam AttribOut Particle attribute type with getView()(p) -> scalar
147 : * @tparam Field ippl::Field with rank=Dim nodal coefficients
148 : * @tparam PosAttrib Particle position attribute with getView()(p) -> Vector<T,Dim>
149 : * @tparam Space Lagrange space providing element/DOF/topology queries
150 : * @tparam policy_type Kokkos execution policy (defaults to Field::execution_space)
151 : */
152 : template <typename AttribOut, typename Field, typename PosAttrib, typename Space,
153 : typename policy_type = Kokkos::RangePolicy<typename Field::execution_space>>
154 12 : inline void interpolate_to_diracs(AttribOut& attrib_out,
155 : const Field& coeffs,
156 : const PosAttrib& pp,
157 : const Space& space,
158 : policy_type iteration_policy)
159 : {
160 12 : constexpr unsigned Dim = Field::dim;
161 : using T = typename AttribOut::value_type;
162 : using field_value_type = typename Field::value_type;
163 : using view_type = typename Field::view_type;
164 : using mesh_type = typename Field::Mesh_t;
165 :
166 12 : static IpplTimings::TimerRef timer =
167 12 : IpplTimings::getTimer("interpolate_field_to_particles(P1)");
168 12 : IpplTimings::startTimer(timer);
169 :
170 12 : view_type view = coeffs.getView();
171 12 : const mesh_type& M = coeffs.get_mesh();
172 :
173 :
174 12 : const ippl::FieldLayout<Dim>& layout = coeffs.getLayout();
175 12 : const ippl::NDIndex<Dim>& lDom = layout.getLocalNDIndex();
176 12 : const int nghost = coeffs.getNghost();
177 :
178 : // Particle device views
179 12 : auto d_pos = pp.getView();
180 12 : auto d_out = attrib_out.getView();
181 :
182 12 : Kokkos::parallel_for("interpolate_to_diracs_P1", iteration_policy,
183 12024 : KOKKOS_LAMBDA(const size_t p) {
184 :
185 18000 : const Vector<T,Dim> x = d_pos(p);
186 :
187 12000 : ippl::Vector<size_t,Dim> e_nd;
188 12000 : ippl::Vector<T,Dim> xi;
189 12000 : locate_element_nd_and_xi<T,Dim>(M, x, e_nd, xi);
190 12000 : const size_t e_lin = space.getElementIndex(e_nd);
191 :
192 12000 : const auto dofs = space.getGlobalDOFIndices(e_lin);
193 :
194 12000 : field_value_type up = field_value_type(0);
195 :
196 68000 : for (size_t a = 0; a < dofs.dim; ++a) {
197 56000 : const size_t local = space.getLocalDOFIndex(e_lin, dofs[a]);
198 56000 : const field_value_type w = space.evaluateRefElementShapeFunction(local, xi);
199 :
200 56000 : const auto v_nd = space.getMeshVertexNDIndex(dofs[a]);
201 56000 : ippl::Vector<size_t,Dim> I;
202 192000 : for (unsigned d = 0; d < Dim; ++d) {
203 136000 : I[d] = static_cast<size_t>(v_nd[d] - lDom.first()[d] + nghost);
204 : }
205 :
206 56000 : up += view_ref<Dim>(view, I) * w;
207 : }
208 12000 : d_out(p) = static_cast<T>(up);
209 12000 : });
210 12 : }
211 :
212 : } // namespace ippl
213 : #endif
|
