Line data Source code
1 :
2 : namespace ippl {
3 :
4 : // LagrangeSpace constructor, which calls the FiniteElementSpace constructor,
5 : // and decomposes the elements among ranks according to layout.
6 : template <typename T, unsigned Dim, unsigned Order, typename ElementType,
7 : typename QuadratureType, typename FieldLHS, typename FieldRHS>
8 396 : LagrangeSpace<T, Dim, Order, ElementType, QuadratureType, FieldLHS, FieldRHS>::LagrangeSpace(
9 : UniformCartesian<T, Dim>& mesh, ElementType& ref_element, const QuadratureType& quadrature,
10 : const Layout_t& layout)
11 : : FiniteElementSpace<T, Dim, getLagrangeNumElementDOFs(Dim, Order), ElementType,
12 396 : QuadratureType, FieldLHS, FieldRHS>(mesh, ref_element, quadrature) {
13 : // Assert that the dimension is either 1, 2 or 3.
14 : static_assert(Dim >= 1 && Dim <= 3,
15 : "Finite Element space only supports 1D, 2D and 3D meshes");
16 :
17 : // Initialize the elementIndices view
18 396 : initializeElementIndices(layout);
19 396 : }
20 :
21 : // LagrangeSpace constructor, which calls the FiniteElementSpace constructor.
22 : template <typename T, unsigned Dim, unsigned Order, typename ElementType,
23 : typename QuadratureType, typename FieldLHS, typename FieldRHS>
24 : LagrangeSpace<T, Dim, Order, ElementType, QuadratureType, FieldLHS, FieldRHS>::LagrangeSpace(
25 : UniformCartesian<T, Dim>& mesh, ElementType& ref_element, const QuadratureType& quadrature)
26 : : FiniteElementSpace<T, Dim, getLagrangeNumElementDOFs(Dim, Order), ElementType,
27 : QuadratureType, FieldLHS, FieldRHS>(mesh, ref_element, quadrature) {
28 : // Assert that the dimension is either 1, 2 or 3.
29 : static_assert(Dim >= 1 && Dim <= 3,
30 : "Finite Element space only supports 1D, 2D and 3D meshes");
31 : }
32 :
33 : // LagrangeSpace initializer, to be made available to the FEMPoissonSolver
34 : // such that we can call it from setRhs.
35 : // Sets the correct mesh ad decomposes the elements among ranks according to layout.
36 : template <typename T, unsigned Dim, unsigned Order, typename ElementType,
37 : typename QuadratureType, typename FieldLHS, typename FieldRHS>
38 0 : void LagrangeSpace<T, Dim, Order, ElementType, QuadratureType, FieldLHS, FieldRHS>::initialize(
39 : UniformCartesian<T, Dim>& mesh, const Layout_t& layout)
40 : {
41 : FiniteElementSpace<T, Dim, getLagrangeNumElementDOFs(Dim, Order), ElementType,
42 0 : QuadratureType, FieldLHS, FieldRHS>::setMesh(mesh);
43 :
44 : // Initialize the elementIndices view
45 0 : initializeElementIndices(layout);
46 0 : }
47 :
48 : // Initialize element indices Kokkos View by distributing elements among MPI ranks.
49 : template <typename T, unsigned Dim, unsigned Order, typename ElementType,
50 : typename QuadratureType, typename FieldLHS, typename FieldRHS>
51 396 : void LagrangeSpace<T, Dim, Order, ElementType, QuadratureType, FieldLHS,
52 : FieldRHS>::initializeElementIndices(const Layout_t& layout) {
53 396 : const auto& ldom = layout.getLocalNDIndex();
54 396 : int npoints = ldom.size();
55 396 : auto first = ldom.first();
56 396 : auto last = ldom.last();
57 396 : ippl::Vector<double, Dim> bounds;
58 :
59 1224 : for (size_t d = 0; d < Dim; ++d) {
60 828 : bounds[d] = this->nr_m[d] - 1;
61 : }
62 :
63 396 : int upperBoundaryPoints = -1;
64 :
65 : // We iterate over the local domain points, getting the corresponding elements,
66 : // while tagging upper boundary points such that they can be removed after.
67 396 : Kokkos::View<size_t*> points("npoints", npoints);
68 396 : Kokkos::parallel_reduce(
69 0 : "ComputePoints", npoints,
70 8514 : KOKKOS_CLASS_LAMBDA(const int i, int& local) {
71 7722 : int idx = i;
72 7722 : indices_t val;
73 7722 : bool isBoundary = false;
74 29070 : for (unsigned int d = 0; d < Dim; ++d) {
75 21348 : int range = last[d] - first[d] + 1;
76 21348 : val[d] = first[d] + (idx % range);
77 21348 : idx /= range;
78 21348 : if (val[d] == bounds[d]) {
79 4716 : isBoundary = true;
80 : }
81 : }
82 7722 : points(i) = (!isBoundary) * (this->getElementIndex(val));
83 7722 : local += isBoundary;
84 7722 : },
85 792 : Kokkos::Sum<int>(upperBoundaryPoints));
86 396 : Kokkos::fence();
87 :
88 : // The elementIndices will be the same array as computed above,
89 : // with the tagged upper boundary points removed.
90 396 : int elementsPerRank = npoints - upperBoundaryPoints;
91 396 : elementIndices = Kokkos::View<size_t*>("i", elementsPerRank);
92 396 : Kokkos::View<size_t> index("index");
93 :
94 396 : Kokkos::parallel_for(
95 16236 : "RemoveNaNs", npoints, KOKKOS_CLASS_LAMBDA(const int i) {
96 15444 : if ((points(i) != 0) || (i == 0)) {
97 7848 : const size_t idx = Kokkos::atomic_fetch_add(&index(), 1);
98 11772 : elementIndices(idx) = points(i);
99 : }
100 : });
101 396 : }
102 :
103 : ///////////////////////////////////////////////////////////////////////
104 : /// Degree of Freedom operations //////////////////////////////////////
105 : ///////////////////////////////////////////////////////////////////////
106 :
107 : template <typename T, unsigned Dim, unsigned Order, typename ElementType,
108 : typename QuadratureType, typename FieldLHS, typename FieldRHS>
109 24 : KOKKOS_FUNCTION size_t LagrangeSpace<T, Dim, Order, ElementType, QuadratureType, FieldLHS,
110 : FieldRHS>::numGlobalDOFs() const {
111 24 : size_t num_global_dofs = 1;
112 72 : for (size_t d = 0; d < Dim; ++d) {
113 48 : num_global_dofs *= this->nr_m[d] * Order;
114 : }
115 :
116 24 : return num_global_dofs;
117 : }
118 :
119 : template <typename T, unsigned Dim, unsigned Order, typename ElementType,
120 : typename QuadratureType, typename FieldLHS, typename FieldRHS>
121 : KOKKOS_FUNCTION
122 1032 : size_t LagrangeSpace<T, Dim, Order, ElementType, QuadratureType, FieldLHS, FieldRHS>::getLocalDOFIndex
123 : (const size_t& elementIndex, const size_t& globalDOFIndex) const {
124 : static_assert(Dim == 1 || Dim == 2 || Dim == 3, "Dim must be 1, 2 or 3");
125 : // TODO fix not order independent, only works for order 1
126 : static_assert(Order == 1, "Only order 1 is supported at the moment");
127 :
128 : // Get all the global DOFs for the element
129 1032 : const Vector<size_t, numElementDOFs> global_dofs =
130 1032 : this->getGlobalDOFIndices(elementIndex);
131 :
132 1032 : ippl::Vector<size_t, numElementDOFs> dof_mapping;
133 : if (Dim == 1) {
134 24 : dof_mapping = {0, 1};
135 : } else if (Dim == 2) {
136 144 : dof_mapping = {0, 1, 3, 2};
137 : } else if (Dim == 3) {
138 864 : dof_mapping = {0, 1, 3, 2, 4, 5, 7, 6};
139 : }
140 :
141 : // Find the global DOF in the vector and return the local DOF index
142 : // TODO this can be done faster since the global DOFs are sorted
143 7232 : for (size_t i = 0; i < dof_mapping.dim; ++i) {
144 6536 : if (global_dofs[dof_mapping[i]] == globalDOFIndex) {
145 336 : return dof_mapping[i];
146 : }
147 : }
148 696 : return std::numeric_limits<size_t>::quiet_NaN();
149 : //throw IpplException("LagrangeSpace::getLocalDOFIndex()",
150 : // "FEM Lagrange Space: Global DOF not found in specified element");
151 1032 : }
152 :
153 : template <typename T, unsigned Dim, unsigned Order, typename ElementType,
154 : typename QuadratureType, typename FieldLHS, typename FieldRHS>
155 : KOKKOS_FUNCTION size_t
156 336 : LagrangeSpace<T, Dim, Order, ElementType, QuadratureType, FieldLHS,
157 : FieldRHS>::getGlobalDOFIndex(const size_t& elementIndex,
158 : const size_t& localDOFIndex) const {
159 336 : const auto global_dofs = this->getGlobalDOFIndices(elementIndex);
160 :
161 672 : return global_dofs[localDOFIndex];
162 336 : }
163 :
164 : template <typename T, unsigned Dim, unsigned Order, typename ElementType,
165 : typename QuadratureType, typename FieldLHS, typename FieldRHS>
166 : KOKKOS_FUNCTION Vector<size_t, LagrangeSpace<T, Dim, Order, ElementType, QuadratureType,
167 : FieldLHS, FieldRHS>::numElementDOFs>
168 452 : LagrangeSpace<T, Dim, Order, ElementType, QuadratureType, FieldLHS,
169 : FieldRHS>::getLocalDOFIndices() const {
170 452 : Vector<size_t, numElementDOFs> localDOFs;
171 :
172 3756 : for (size_t dof = 0; dof < numElementDOFs; ++dof) {
173 3304 : localDOFs[dof] = dof;
174 : }
175 :
176 452 : return localDOFs;
177 : }
178 :
179 : template <typename T, unsigned Dim, unsigned Order, typename ElementType,
180 : typename QuadratureType, typename FieldLHS, typename FieldRHS>
181 : KOKKOS_FUNCTION Vector<size_t, LagrangeSpace<T, Dim, Order, ElementType, QuadratureType,
182 : FieldLHS, FieldRHS>::numElementDOFs>
183 1820 : LagrangeSpace<T, Dim, Order, ElementType, QuadratureType, FieldLHS,
184 : FieldRHS>::getGlobalDOFIndices(const size_t& elementIndex) const {
185 1820 : Vector<size_t, numElementDOFs> globalDOFs(0);
186 :
187 : // get element pos
188 1820 : indices_t elementPos = this->getElementNDIndex(elementIndex);
189 :
190 : // Compute the vector to multiply the ndindex with
191 1820 : ippl::Vector<size_t, Dim> vec(1);
192 5080 : for (size_t d = 1; d < dim; ++d) {
193 8028 : for (size_t d2 = d; d2 < Dim; ++d2) {
194 4768 : vec[d2] *= this->nr_m[d - 1];
195 : }
196 : }
197 1820 : vec *= Order; // Multiply each dimension by the order
198 1820 : size_t smallestGlobalDOF = elementPos.dot(vec);
199 :
200 : // Add the vertex DOFs
201 1820 : globalDOFs[0] = smallestGlobalDOF;
202 1820 : globalDOFs[1] = smallestGlobalDOF + Order;
203 :
204 : if (Dim >= 2) {
205 1752 : globalDOFs[2] = globalDOFs[1] + this->nr_m[1] * Order;
206 1752 : globalDOFs[3] = globalDOFs[0] + this->nr_m[1] * Order;
207 : }
208 : if (Dim >= 3) {
209 1508 : globalDOFs[4] = globalDOFs[0] + this->nr_m[1] * this->nr_m[2] * Order;
210 1508 : globalDOFs[5] = globalDOFs[1] + this->nr_m[1] * this->nr_m[2] * Order;
211 1508 : globalDOFs[6] = globalDOFs[2] + this->nr_m[1] * this->nr_m[2] * Order;
212 1508 : globalDOFs[7] = globalDOFs[3] + this->nr_m[1] * this->nr_m[2] * Order;
213 : }
214 :
215 : if (Order > 1) {
216 : // If the order is greater than 1, there are edge and face DOFs, otherwise the work is
217 : // done
218 :
219 : // Add the edge DOFs
220 : if (Dim >= 2) {
221 : for (size_t i = 0; i < Order - 1; ++i) {
222 : globalDOFs[8 + i] = globalDOFs[0] + i + 1;
223 : globalDOFs[8 + Order - 1 + i] = globalDOFs[1] + (i + 1) * this->nr_m[1];
224 : globalDOFs[8 + 2 * (Order - 1) + i] = globalDOFs[2] - (i + 1);
225 : globalDOFs[8 + 3 * (Order - 1) + i] = globalDOFs[3] - (i + 1) * this->nr_m[1];
226 : }
227 : }
228 : if (Dim >= 3) {
229 : // TODO
230 : }
231 :
232 : // Add the face DOFs
233 : if (Dim >= 2) {
234 : for (size_t i = 0; i < Order - 1; ++i) {
235 : for (size_t j = 0; j < Order - 1; ++j) {
236 : // TODO CHECK
237 : globalDOFs[8 + 4 * (Order - 1) + i * (Order - 1) + j] =
238 : globalDOFs[0] + (i + 1) + (j + 1) * this->nr_m[1];
239 : globalDOFs[8 + 4 * (Order - 1) + (Order - 1) * (Order - 1) + i * (Order - 1)
240 : + j] = globalDOFs[1] + (i + 1) + (j + 1) * this->nr_m[1];
241 : globalDOFs[8 + 4 * (Order - 1) + 2 * (Order - 1) * (Order - 1)
242 : + i * (Order - 1) + j] =
243 : globalDOFs[2] - (i + 1) + (j + 1) * this->nr_m[1];
244 : globalDOFs[8 + 4 * (Order - 1) + 3 * (Order - 1) * (Order - 1)
245 : + i * (Order - 1) + j] =
246 : globalDOFs[3] - (i + 1) + (j + 1) * this->nr_m[1];
247 : }
248 : }
249 : }
250 : }
251 :
252 1820 : return globalDOFs;
253 1820 : }
254 :
255 : ///////////////////////////////////////////////////////////////////////
256 : /// Basis functions and gradients /////////////////////////////////////
257 : ///////////////////////////////////////////////////////////////////////
258 :
259 : template <typename T, unsigned Dim, unsigned Order, typename ElementType,
260 : typename QuadratureType, typename FieldLHS, typename FieldRHS>
261 : KOKKOS_FUNCTION T
262 267040 : LagrangeSpace<T, Dim, Order, ElementType, QuadratureType, FieldLHS, FieldRHS>::
263 : evaluateRefElementShapeFunction(
264 : const size_t& localDOF,
265 : const LagrangeSpace<T, Dim, Order, ElementType, QuadratureType, FieldLHS,
266 : FieldRHS>::point_t& localPoint) const {
267 : static_assert(Order == 1, "Only order 1 is supported at the moment");
268 : // Assert that the local vertex index is valid.
269 267040 : assert(localDOF < numElementDOFs
270 : && "The local vertex index is invalid"); // TODO assumes 1st order Lagrange
271 :
272 267040 : assert(this->ref_element_m.isPointInRefElement(localPoint)
273 : && "Point is not in reference element");
274 :
275 : // Get the local vertex indices for the local vertex index.
276 : // TODO fix not order independent, only works for order 1
277 267040 : const point_t ref_element_point = this->ref_element_m.getLocalVertices()[localDOF];
278 :
279 : // The variable that accumulates the product of the shape functions.
280 267040 : T product = 1;
281 :
282 1060944 : for (size_t d = 0; d < Dim; d++) {
283 793904 : if (localPoint[d] < ref_element_point[d]) {
284 396952 : product *= localPoint[d];
285 : } else {
286 396952 : product *= 1.0 - localPoint[d];
287 : }
288 : }
289 :
290 267040 : return product;
291 267040 : }
292 :
293 : template <typename T, unsigned Dim, unsigned Order, typename ElementType,
294 : typename QuadratureType, typename FieldLHS, typename FieldRHS>
295 : KOKKOS_FUNCTION typename LagrangeSpace<T, Dim, Order, ElementType, QuadratureType, FieldLHS,
296 : FieldRHS>::point_t
297 262600 : LagrangeSpace<T, Dim, Order, ElementType, QuadratureType, FieldLHS, FieldRHS>::
298 : evaluateRefElementShapeFunctionGradient(
299 : const size_t& localDOF,
300 : const LagrangeSpace<T, Dim, Order, ElementType, QuadratureType, FieldLHS,
301 : FieldRHS>::point_t& localPoint) const {
302 : // TODO fix not order independent, only works for order 1
303 : static_assert(Order == 1 && "Only order 1 is supported at the moment");
304 :
305 : // Assert that the local vertex index is valid.
306 262600 : assert(localDOF < numElementDOFs && "The local vertex index is invalid");
307 :
308 262600 : assert(this->ref_element_m.isPointInRefElement(localPoint)
309 : && "Point is not in reference element");
310 :
311 : // Get the local dof nd_index
312 262600 : const vertex_points_t local_vertex_points = this->ref_element_m.getLocalVertices();
313 :
314 262600 : const point_t& local_vertex_point = local_vertex_points[localDOF];
315 :
316 262600 : point_t gradient(1);
317 :
318 : // To construct the gradient we need to loop over the dimensions and multiply the
319 : // shape functions in each dimension except the current one. The one of the current
320 : // dimension is replaced by the derivative of the shape function in that dimension,
321 : // which is either 1 or -1.
322 1043664 : for (size_t d = 0; d < Dim; d++) {
323 : // The variable that accumulates the product of the shape functions.
324 781064 : T product = 1;
325 :
326 3111120 : for (size_t d2 = 0; d2 < Dim; d2++) {
327 2330056 : if (d2 == d) {
328 781064 : if (localPoint[d] < local_vertex_point[d]) {
329 390532 : product *= 1;
330 : } else {
331 390532 : product *= -1;
332 : }
333 : } else {
334 1548992 : if (localPoint[d2] < local_vertex_point[d2]) {
335 774496 : product *= localPoint[d2];
336 : } else {
337 774496 : product *= 1.0 - localPoint[d2];
338 : }
339 : }
340 : }
341 :
342 781064 : gradient[d] = product;
343 : }
344 :
345 262600 : return gradient;
346 262600 : }
347 :
348 : ///////////////////////////////////////////////////////////////////////
349 : /// Assembly operations ///////////////////////////////////////////////
350 : ///////////////////////////////////////////////////////////////////////
351 :
352 : template <typename T, unsigned Dim, unsigned Order, typename ElementType,
353 : typename QuadratureType, typename FieldLHS, typename FieldRHS>
354 : template <typename F>
355 8 : FieldLHS LagrangeSpace<T, Dim, Order, ElementType, QuadratureType, FieldLHS,
356 : FieldRHS>::evaluateAx(FieldLHS& field, F& evalFunction) const {
357 8 : Inform m("");
358 :
359 : // start a timer
360 8 : static IpplTimings::TimerRef evalAx = IpplTimings::getTimer("evaluateAx");
361 8 : IpplTimings::startTimer(evalAx);
362 :
363 : // get number of ghost cells in field
364 8 : const int nghost = field.getNghost();
365 :
366 : // create a new field for result with view initialized to zero (views are initialized to
367 : // zero by default)
368 8 : FieldLHS resultField(field.get_mesh(), field.getLayout(), nghost);
369 :
370 : // List of quadrature weights
371 8 : const Vector<T, QuadratureType::numElementNodes> w =
372 8 : this->quadrature_m.getWeightsForRefElement();
373 :
374 : // List of quadrature nodes
375 8 : const Vector<point_t, QuadratureType::numElementNodes> q =
376 8 : this->quadrature_m.getIntegrationNodesForRefElement();
377 :
378 : // TODO move outside of evaluateAx (I think it is possible for other problems as well)
379 : // Gradients of the basis functions for the DOF at the quadrature nodes
380 8 : Vector<Vector<point_t, numElementDOFs>, QuadratureType::numElementNodes> grad_b_q;
381 16 : for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
382 48 : for (size_t i = 0; i < numElementDOFs; ++i) {
383 40 : grad_b_q[k][i] = this->evaluateRefElementShapeFunctionGradient(i, q[k]);
384 : }
385 : }
386 :
387 : // Get field data and atomic result data,
388 : // since it will be added to during the kokkos loop
389 8 : ViewType view = field.getView();
390 8 : AtomicViewType resultView = resultField.getView();
391 :
392 : // Get boundary conditions from field
393 8 : BConds<FieldLHS, Dim>& bcField = field.getFieldBC();
394 8 : FieldBC bcType = bcField[0]->getBCType();
395 :
396 : // Get domain information
397 8 : auto ldom = (field.getLayout()).getLocalNDIndex();
398 :
399 : using exec_space = typename Kokkos::View<const size_t*>::execution_space;
400 : using policy_type = Kokkos::RangePolicy<exec_space>;
401 :
402 : // start a timer
403 8 : static IpplTimings::TimerRef outer_loop = IpplTimings::getTimer("evaluateAx: outer loop");
404 8 : IpplTimings::startTimer(outer_loop);
405 :
406 : // Loop over elements to compute contributions
407 8 : Kokkos::parallel_for(
408 16 : "Loop over elements", policy_type(0, elementIndices.extent(0)),
409 152 : KOKKOS_CLASS_LAMBDA(const size_t index) {
410 0 : const size_t elementIndex = elementIndices(index);
411 136 : const Vector<size_t, numElementDOFs> local_dof = this->getLocalDOFIndices();
412 0 : const Vector<size_t, numElementDOFs> global_dofs =
413 136 : this->getGlobalDOFIndices(elementIndex);
414 136 : Vector<indices_t, numElementDOFs> global_dof_ndindices;
415 :
416 1176 : for (size_t i = 0; i < numElementDOFs; ++i) {
417 1040 : global_dof_ndindices[i] = this->getMeshVertexNDIndex(global_dofs[i]);
418 : }
419 :
420 : // local DOF indices
421 : size_t i, j;
422 :
423 : // Element matrix
424 136 : Vector<Vector<T, numElementDOFs>, numElementDOFs> A_K;
425 :
426 : // 1. Compute the Galerkin element matrix A_K
427 1176 : for (i = 0; i < numElementDOFs; ++i) {
428 9264 : for (j = 0; j < numElementDOFs; ++j) {
429 0 : A_K[i][j] = 0.0;
430 16448 : for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
431 8224 : A_K[i][j] += w[k] * evalFunction(i, j, grad_b_q[k]);
432 : }
433 : }
434 : }
435 :
436 : // global DOF n-dimensional indices (Vector of N indices representing indices in
437 : // each dimension)
438 136 : indices_t I_nd, J_nd;
439 :
440 : // 2. Compute the contribution to resultAx = A*x with A_K
441 1176 : for (i = 0; i < numElementDOFs; ++i) {
442 0 : I_nd = global_dof_ndindices[i];
443 :
444 : // Handle boundary DOFs
445 : // If Zero Dirichlet BCs, skip this DOF
446 : // If Constant Dirichlet BCs, identity
447 1040 : if ((bcType == CONSTANT_FACE) && (this->isDOFOnBoundary(I_nd))) {
448 0 : for (unsigned d = 0; d < Dim; ++d) {
449 0 : I_nd[d] = I_nd[d] - ldom[d].first() + nghost;
450 : }
451 0 : apply(resultView, I_nd) = apply(view, I_nd);
452 0 : continue;
453 1040 : } else if ((bcType == ZERO_FACE) && (this->isDOFOnBoundary(I_nd))) {
454 0 : continue;
455 : }
456 :
457 : // get the appropriate index for the Kokkos view of the field
458 1752 : for (unsigned d = 0; d < Dim; ++d) {
459 0 : I_nd[d] = I_nd[d] - ldom[d].first() + nghost;
460 : }
461 :
462 3924 : for (j = 0; j < numElementDOFs; ++j) {
463 0 : J_nd = global_dof_ndindices[j];
464 :
465 : // Skip boundary DOFs (Zero & Constant Dirichlet BCs)
466 0 : if (((bcType == ZERO_FACE) || (bcType == CONSTANT_FACE))
467 3480 : && this->isDOFOnBoundary(J_nd)) {
468 0 : continue;
469 : }
470 :
471 : // get the appropriate index for the Kokkos view of the field
472 8040 : for (unsigned d = 0; d < Dim; ++d) {
473 0 : J_nd[d] = J_nd[d] - ldom[d].first() + nghost;
474 : }
475 :
476 2020 : apply(resultView, I_nd) += A_K[i][j] * apply(view, J_nd);
477 : }
478 : }
479 0 : });
480 8 : IpplTimings::stopTimer(outer_loop);
481 :
482 8 : if (bcType == PERIODIC_FACE) {
483 0 : resultField.accumulateHalo();
484 0 : bcField.apply(resultField);
485 0 : bcField.assignGhostToPhysical(resultField);
486 : } else {
487 8 : resultField.accumulateHalo_noghost();
488 : }
489 :
490 8 : IpplTimings::stopTimer(evalAx);
491 :
492 8 : return resultField;
493 8 : }
494 :
495 : template <typename T, unsigned Dim, unsigned Order, typename ElementType,
496 : typename QuadratureType, typename FieldLHS, typename FieldRHS>
497 : template <typename F>
498 0 : FieldLHS LagrangeSpace<T, Dim, Order, ElementType, QuadratureType, FieldLHS,
499 : FieldRHS>::evaluateAx_lift(FieldLHS& field, F& evalFunction) const {
500 0 : Inform m("");
501 :
502 : // get number of ghost cells in field
503 0 : const int nghost = field.getNghost();
504 :
505 : // create a new field for result with view initialized to zero (views are initialized to
506 : // zero by default)
507 0 : FieldLHS resultField(field.get_mesh(), field.getLayout(), nghost);
508 :
509 : // List of quadrature weights
510 0 : const Vector<T, QuadratureType::numElementNodes> w =
511 0 : this->quadrature_m.getWeightsForRefElement();
512 :
513 : // List of quadrature nodes
514 0 : const Vector<point_t, QuadratureType::numElementNodes> q =
515 0 : this->quadrature_m.getIntegrationNodesForRefElement();
516 :
517 : // TODO move outside of evaluateAx (I think it is possible for other problems as well)
518 : // Gradients of the basis functions for the DOF at the quadrature nodes
519 0 : Vector<Vector<point_t, numElementDOFs>, QuadratureType::numElementNodes> grad_b_q;
520 0 : for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
521 0 : for (size_t i = 0; i < numElementDOFs; ++i) {
522 0 : grad_b_q[k][i] = this->evaluateRefElementShapeFunctionGradient(i, q[k]);
523 : }
524 : }
525 :
526 : // Get field data and atomic result data,
527 : // since it will be added to during the kokkos loop
528 0 : ViewType view = field.getView();
529 0 : AtomicViewType resultView = resultField.getView();
530 :
531 : // Get domain information
532 0 : auto ldom = (field.getLayout()).getLocalNDIndex();
533 :
534 : using exec_space = typename Kokkos::View<const size_t*>::execution_space;
535 : using policy_type = Kokkos::RangePolicy<exec_space>;
536 :
537 : // Loop over elements to compute contributions
538 0 : Kokkos::parallel_for(
539 0 : "Loop over elements", policy_type(0, elementIndices.extent(0)),
540 0 : KOKKOS_CLASS_LAMBDA(const size_t index) {
541 0 : const size_t elementIndex = elementIndices(index);
542 0 : const Vector<size_t, numElementDOFs> local_dof = this->getLocalDOFIndices();
543 0 : const Vector<size_t, numElementDOFs> global_dofs =
544 0 : this->getGlobalDOFIndices(elementIndex);
545 0 : Vector<indices_t, numElementDOFs> global_dof_ndindices;
546 :
547 0 : for (size_t i = 0; i < numElementDOFs; ++i) {
548 0 : global_dof_ndindices[i] = this->getMeshVertexNDIndex(global_dofs[i]);
549 : }
550 :
551 : // local DOF indices
552 : size_t i, j;
553 :
554 : // Element matrix
555 0 : Vector<Vector<T, numElementDOFs>, numElementDOFs> A_K;
556 :
557 : // 1. Compute the Galerkin element matrix A_K
558 0 : for (i = 0; i < numElementDOFs; ++i) {
559 0 : for (j = 0; j < numElementDOFs; ++j) {
560 0 : A_K[i][j] = 0.0;
561 0 : for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
562 0 : A_K[i][j] += w[k] * evalFunction(i, j, grad_b_q[k]);
563 : }
564 : }
565 : }
566 :
567 : // global DOF n-dimensional indices (Vector of N indices representing indices in
568 : // each dimension)
569 0 : indices_t I_nd, J_nd;
570 :
571 : // 2. Compute the contribution to resultAx = A*x with A_K
572 0 : for (i = 0; i < numElementDOFs; ++i) {
573 0 : I_nd = global_dof_ndindices[i];
574 :
575 : // Skip if on a row of the matrix
576 0 : if (this->isDOFOnBoundary(I_nd)) {
577 0 : continue;
578 : }
579 :
580 : // get the appropriate index for the Kokkos view of the field
581 0 : for (unsigned d = 0; d < Dim; ++d) {
582 0 : I_nd[d] = I_nd[d] - ldom[d].first() + nghost;
583 : }
584 :
585 0 : for (j = 0; j < numElementDOFs; ++j) {
586 0 : J_nd = global_dof_ndindices[j];
587 :
588 : // Contribute to lifting only if on a boundary DOF
589 0 : if (this->isDOFOnBoundary(J_nd)) {
590 : // get the appropriate index for the Kokkos view of the field
591 0 : for (unsigned d = 0; d < Dim; ++d) {
592 0 : J_nd[d] = J_nd[d] - ldom[d].first() + nghost;
593 : }
594 0 : apply(resultView, I_nd) += A_K[i][j] * apply(view, J_nd);
595 0 : continue;
596 0 : }
597 :
598 : }
599 : }
600 0 : });
601 0 : resultField.accumulateHalo();
602 :
603 0 : return resultField;
604 0 : }
605 :
606 : template <typename T, unsigned Dim, unsigned Order, typename ElementType,
607 : typename QuadratureType, typename FieldLHS, typename FieldRHS>
608 20 : void LagrangeSpace<T, Dim, Order, ElementType, QuadratureType, FieldLHS,
609 : FieldRHS>::evaluateLoadVector(FieldRHS& field) const {
610 20 : Inform m("");
611 :
612 : // start a timer
613 20 : static IpplTimings::TimerRef evalLoadV = IpplTimings::getTimer("evaluateLoadVector");
614 20 : IpplTimings::startTimer(evalLoadV);
615 :
616 : // List of quadrature weights
617 20 : const Vector<T, QuadratureType::numElementNodes> w =
618 20 : this->quadrature_m.getWeightsForRefElement();
619 :
620 : // List of quadrature nodes
621 20 : const Vector<point_t, QuadratureType::numElementNodes> q =
622 20 : this->quadrature_m.getIntegrationNodesForRefElement();
623 :
624 20 : const indices_t zeroNdIndex = Vector<size_t, Dim>(0);
625 :
626 : // Evaluate the basis functions for the DOF at the quadrature nodes
627 20 : Vector<Vector<T, numElementDOFs>, QuadratureType::numElementNodes> basis_q;
628 648 : for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
629 5108 : for (size_t i = 0; i < numElementDOFs; ++i) {
630 4480 : basis_q[k][i] = this->evaluateRefElementShapeFunction(i, q[k]);
631 : }
632 : }
633 :
634 : // Absolute value of det Phi_K
635 20 : const T absDetDPhi = Kokkos::abs(this->ref_element_m.getDeterminantOfTransformationJacobian(
636 40 : this->getElementMeshVertexPoints(zeroNdIndex)));
637 :
638 : // Get domain information and ghost cells
639 20 : auto ldom = (field.getLayout()).getLocalNDIndex();
640 20 : const int nghost = field.getNghost();
641 :
642 : // Get boundary conditions from field
643 20 : BConds<FieldRHS, Dim>& bcField = field.getFieldBC();
644 20 : FieldBC bcType = bcField[0]->getBCType();
645 :
646 20 : FieldRHS temp_field(field.get_mesh(), field.getLayout(), nghost);
647 20 : temp_field.setFieldBC(bcField);
648 :
649 : // Get field data and make it atomic,
650 : // since it will be added to during the kokkos loop
651 : // We work with a temporary field since we need to use field
652 : // to evaluate the load vector; then we assign temp to RHS field
653 20 : AtomicViewType atomic_view = temp_field.getView();
654 :
655 : using exec_space = typename Kokkos::View<const size_t*>::execution_space;
656 : using policy_type = Kokkos::RangePolicy<exec_space>;
657 :
658 : // start a timer
659 20 : static IpplTimings::TimerRef outer_loop =
660 20 : IpplTimings::getTimer("evaluateLoadVec: outer loop");
661 20 : IpplTimings::startTimer(outer_loop);
662 :
663 : // Loop over elements to compute contributions
664 20 : Kokkos::parallel_for(
665 40 : "Loop over elements", policy_type(0, elementIndices.extent(0)),
666 344 : KOKKOS_CLASS_LAMBDA(size_t index) {
667 0 : const size_t elementIndex = elementIndices(index);
668 304 : const Vector<size_t, numElementDOFs> local_dofs = this->getLocalDOFIndices();
669 0 : const Vector<size_t, numElementDOFs> global_dofs =
670 304 : this->getGlobalDOFIndices(elementIndex);
671 :
672 : size_t i, I;
673 :
674 : // 1. Compute b_K
675 4720 : for (i = 0; i < numElementDOFs; ++i) {
676 0 : I = global_dofs[i];
677 :
678 : // TODO fix for higher order
679 2208 : auto dof_ndindex_I = this->getMeshVertexNDIndex(I);
680 :
681 : // Skip boundary DOFs (Zero and Constant Dirichlet BCs)
682 0 : if (((bcType == ZERO_FACE) || (bcType == CONSTANT_FACE))
683 2208 : && (this->isDOFOnBoundary(dof_ndindex_I))) {
684 0 : continue;
685 : }
686 :
687 : // calculate the contribution of this element
688 0 : T contrib = 0;
689 57264 : for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
690 0 : T val = 0;
691 499104 : for (size_t j = 0; j < numElementDOFs; ++j) {
692 : // get field index corresponding to this DOF
693 0 : size_t J = global_dofs[j];
694 442800 : auto dof_ndindex_J = this->getMeshVertexNDIndex(J);
695 1763712 : for (unsigned d = 0; d < Dim; ++d) {
696 0 : dof_ndindex_J[d] = dof_ndindex_J[d] - ldom[d].first() + nghost;
697 : }
698 :
699 : // get field value at DOF and interpolate to q_k
700 442800 : val += basis_q[k][j] * apply(field, dof_ndindex_J);
701 : }
702 :
703 0 : contrib += w[k] * basis_q[k][i] * absDetDPhi * val;
704 : }
705 :
706 : // get the appropriate index for the Kokkos view of the field
707 3720 : for (unsigned d = 0; d < Dim; ++d) {
708 0 : dof_ndindex_I[d] = dof_ndindex_I[d] - ldom[d].first() + nghost;
709 : }
710 :
711 : // add the contribution of the element to the field
712 960 : apply(atomic_view, dof_ndindex_I) += contrib;
713 :
714 : }
715 0 : });
716 20 : IpplTimings::stopTimer(outer_loop);
717 :
718 20 : temp_field.accumulateHalo();
719 :
720 20 : if ((bcType == PERIODIC_FACE) || (bcType == CONSTANT_FACE)) {
721 0 : bcField.apply(temp_field);
722 0 : bcField.assignGhostToPhysical(temp_field);
723 : }
724 :
725 20 : field = temp_field;
726 :
727 20 : IpplTimings::stopTimer(evalLoadV);
728 20 : }
729 :
730 : template <typename T, unsigned Dim, unsigned Order, typename ElementType,
731 : typename QuadratureType, typename FieldLHS, typename FieldRHS>
732 : template <typename F>
733 0 : T LagrangeSpace<T, Dim, Order, ElementType, QuadratureType, FieldLHS, FieldRHS>::computeErrorL2(
734 : const FieldLHS& u_h, const F& u_sol) const {
735 0 : if (this->quadrature_m.getOrder() < (2 * Order + 1)) {
736 : // throw exception
737 0 : throw IpplException(
738 : "LagrangeSpace::computeErrorL2()",
739 : "Order of quadrature rule for error computation should be > 2*p + 1");
740 : }
741 :
742 : // List of quadrature weights
743 0 : const Vector<T, QuadratureType::numElementNodes> w =
744 0 : this->quadrature_m.getWeightsForRefElement();
745 :
746 : // List of quadrature nodes
747 0 : const Vector<point_t, QuadratureType::numElementNodes> q =
748 0 : this->quadrature_m.getIntegrationNodesForRefElement();
749 :
750 : // Evaluate the basis functions for the DOF at the quadrature nodes
751 0 : Vector<Vector<T, numElementDOFs>, QuadratureType::numElementNodes> basis_q;
752 0 : for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
753 0 : for (size_t i = 0; i < numElementDOFs; ++i) {
754 0 : basis_q[k][i] = this->evaluateRefElementShapeFunction(i, q[k]);
755 : }
756 : }
757 :
758 0 : const indices_t zeroNdIndex = Vector<size_t, Dim>(0);
759 :
760 : // Absolute value of det Phi_K
761 0 : const T absDetDPhi = Kokkos::abs(this->ref_element_m.getDeterminantOfTransformationJacobian(
762 0 : this->getElementMeshVertexPoints(zeroNdIndex)));
763 :
764 : // Variable to sum the error to
765 0 : T error = 0;
766 :
767 : // Get domain information and ghost cells
768 0 : auto ldom = (u_h.getLayout()).getLocalNDIndex();
769 0 : const int nghost = u_h.getNghost();
770 :
771 : using exec_space = typename Kokkos::View<const size_t*>::execution_space;
772 : using policy_type = Kokkos::RangePolicy<exec_space>;
773 :
774 : // Loop over elements to compute contributions
775 0 : Kokkos::parallel_reduce(
776 0 : "Compute error over elements", policy_type(0, elementIndices.extent(0)),
777 0 : KOKKOS_CLASS_LAMBDA(size_t index, double& local) {
778 0 : const size_t elementIndex = elementIndices(index);
779 0 : const Vector<size_t, numElementDOFs> global_dofs =
780 0 : this->getGlobalDOFIndices(elementIndex);
781 :
782 : // contribution of this element to the error
783 0 : T contrib = 0;
784 0 : for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
785 0 : T val_u_sol = u_sol(this->ref_element_m.localToGlobal(
786 0 : this->getElementMeshVertexPoints(this->getElementNDIndex(elementIndex)),
787 0 : q[k]));
788 :
789 0 : T val_u_h = 0;
790 0 : for (size_t i = 0; i < numElementDOFs; ++i) {
791 : // get field index corresponding to this DOF
792 0 : size_t I = global_dofs[i];
793 0 : auto dof_ndindex_I = this->getMeshVertexNDIndex(I);
794 0 : for (unsigned d = 0; d < Dim; ++d) {
795 0 : dof_ndindex_I[d] = dof_ndindex_I[d] - ldom[d].first() + nghost;
796 : }
797 :
798 : // get field value at DOF and interpolate to q_k
799 0 : val_u_h += basis_q[k][i] * apply(u_h, dof_ndindex_I);
800 : }
801 :
802 0 : contrib += w[k] * Kokkos::pow(val_u_sol - val_u_h, 2) * absDetDPhi;
803 : }
804 0 : local += contrib;
805 0 : },
806 0 : Kokkos::Sum<double>(error));
807 :
808 : // MPI reduce
809 0 : T global_error = 0.0;
810 0 : Comm->allreduce(error, global_error, 1, std::plus<T>());
811 :
812 0 : return Kokkos::sqrt(global_error);
813 0 : }
814 :
815 : template <typename T, unsigned Dim, unsigned Order, typename ElementType,
816 : typename QuadratureType, typename FieldLHS, typename FieldRHS>
817 0 : T LagrangeSpace<T, Dim, Order, ElementType, QuadratureType, FieldLHS, FieldRHS>::computeAvg(
818 : const FieldLHS& u_h) const {
819 0 : if (this->quadrature_m.getOrder() < (2 * Order + 1)) {
820 : // throw exception
821 0 : throw IpplException(
822 : "LagrangeSpace::computeAvg()",
823 : "Order of quadrature rule for error computation should be > 2*p + 1");
824 : }
825 :
826 : // List of quadrature weights
827 0 : const Vector<T, QuadratureType::numElementNodes> w =
828 0 : this->quadrature_m.getWeightsForRefElement();
829 :
830 : // List of quadrature nodes
831 0 : const Vector<point_t, QuadratureType::numElementNodes> q =
832 0 : this->quadrature_m.getIntegrationNodesForRefElement();
833 :
834 : // Evaluate the basis functions for the DOF at the quadrature nodes
835 0 : Vector<Vector<T, numElementDOFs>, QuadratureType::numElementNodes> basis_q;
836 0 : for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
837 0 : for (size_t i = 0; i < numElementDOFs; ++i) {
838 0 : basis_q[k][i] = this->evaluateRefElementShapeFunction(i, q[k]);
839 : }
840 : }
841 :
842 0 : const indices_t zeroNdIndex = Vector<size_t, Dim>(0);
843 :
844 : // Absolute value of det Phi_K
845 0 : const T absDetDPhi = Kokkos::abs(this->ref_element_m.getDeterminantOfTransformationJacobian(
846 0 : this->getElementMeshVertexPoints(zeroNdIndex)));
847 :
848 : // Variable to sum the error to
849 0 : T avg = 0;
850 :
851 : // Get domain information and ghost cells
852 0 : auto ldom = (u_h.getLayout()).getLocalNDIndex();
853 0 : const int nghost = u_h.getNghost();
854 :
855 : using exec_space = typename Kokkos::View<const size_t*>::execution_space;
856 : using policy_type = Kokkos::RangePolicy<exec_space>;
857 :
858 : // Loop over elements to compute contributions
859 0 : Kokkos::parallel_reduce(
860 0 : "Compute average over elements", policy_type(0, elementIndices.extent(0)),
861 0 : KOKKOS_CLASS_LAMBDA(size_t index, double& local) {
862 0 : const size_t elementIndex = elementIndices(index);
863 0 : const Vector<size_t, numElementDOFs> global_dofs =
864 0 : this->getGlobalDOFIndices(elementIndex);
865 :
866 : // contribution of this element to the error
867 0 : T contrib = 0;
868 0 : for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
869 0 : T val_u_h = 0;
870 0 : for (size_t i = 0; i < numElementDOFs; ++i) {
871 : // get field index corresponding to this DOF
872 0 : size_t I = global_dofs[i];
873 0 : auto dof_ndindex_I = this->getMeshVertexNDIndex(I);
874 0 : for (unsigned d = 0; d < Dim; ++d) {
875 0 : dof_ndindex_I[d] = dof_ndindex_I[d] - ldom[d].first() + nghost;
876 : }
877 :
878 : // get field value at DOF and interpolate to q_k
879 0 : val_u_h += basis_q[k][i] * apply(u_h, dof_ndindex_I);
880 : }
881 :
882 0 : contrib += w[k] * val_u_h * absDetDPhi;
883 : }
884 0 : local += contrib;
885 0 : },
886 0 : Kokkos::Sum<double>(avg));
887 :
888 : // MPI reduce
889 0 : T global_avg = 0.0;
890 0 : Comm->allreduce(avg, global_avg, 1, std::plus<T>());
891 :
892 0 : return global_avg;
893 0 : }
894 :
895 : } // namespace ippl
|
