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 12 : KOKKOS_FUNCTION size_t LagrangeSpace<T, Dim, Order, ElementType, QuadratureType, FieldLHS,
110 : FieldRHS>::numGlobalDOFs() const {
111 12 : size_t num_global_dofs = 1;
112 36 : for (size_t d = 0; d < Dim; ++d) {
113 24 : num_global_dofs *= this->nr_m[d] * Order;
114 : }
115 :
116 12 : 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 : // Make local element matrix -- does not change through the element mesh
388 : // Element matrix
389 8 : Vector<Vector<T, numElementDOFs>, numElementDOFs> A_K;
390 :
391 : // 1. Compute the Galerkin element matrix A_K
392 48 : for (size_t i = 0; i < numElementDOFs; ++i) {
393 312 : for (size_t j = 0; j < numElementDOFs; ++j) {
394 272 : A_K[i][j] = 0.0;
395 544 : for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
396 272 : A_K[i][j] += w[k] * evalFunction(i, j, grad_b_q[k]);
397 : }
398 : }
399 : }
400 :
401 : // Get field data and atomic result data,
402 : // since it will be added to during the kokkos loop
403 8 : ViewType view = field.getView();
404 8 : AtomicViewType resultView = resultField.getView();
405 :
406 : // Get boundary conditions from field
407 8 : BConds<FieldLHS, Dim>& bcField = field.getFieldBC();
408 8 : FieldBC bcType = bcField[0]->getBCType();
409 :
410 : // Get domain information
411 8 : auto ldom = (field.getLayout()).getLocalNDIndex();
412 :
413 : using exec_space = typename Kokkos::View<const size_t*>::execution_space;
414 : using policy_type = Kokkos::RangePolicy<exec_space>;
415 :
416 : // start a timer
417 8 : static IpplTimings::TimerRef outer_loop = IpplTimings::getTimer("evaluateAx: outer loop");
418 8 : IpplTimings::startTimer(outer_loop);
419 :
420 : // Loop over elements to compute contributions
421 8 : Kokkos::parallel_for(
422 16 : "Loop over elements", policy_type(0, elementIndices.extent(0)),
423 152 : KOKKOS_CLASS_LAMBDA(const size_t index) {
424 0 : const size_t elementIndex = elementIndices(index);
425 136 : const Vector<size_t, numElementDOFs> local_dof = this->getLocalDOFIndices();
426 0 : const Vector<size_t, numElementDOFs> global_dofs =
427 136 : this->getGlobalDOFIndices(elementIndex);
428 136 : Vector<indices_t, numElementDOFs> global_dof_ndindices;
429 :
430 1176 : for (size_t i = 0; i < numElementDOFs; ++i) {
431 1040 : global_dof_ndindices[i] = this->getMeshVertexNDIndex(global_dofs[i]);
432 : }
433 :
434 : // local DOF indices (both i and j go from 0 to numDOFs-1 in the element)
435 : size_t i, j;
436 :
437 : // global DOF n-dimensional indices (Vector of N indices representing indices in
438 : // each dimension)
439 136 : indices_t I_nd, J_nd;
440 :
441 : // 2. Compute the contribution to resultAx = A*x with A_K
442 1176 : for (i = 0; i < numElementDOFs; ++i) {
443 0 : I_nd = global_dof_ndindices[i];
444 :
445 : // Handle boundary DOFs
446 : // If Zero Dirichlet BCs, skip this DOF
447 : // If Constant Dirichlet BCs, identity
448 1040 : if ((bcType == CONSTANT_FACE) && (this->isDOFOnBoundary(I_nd))) {
449 0 : for (unsigned d = 0; d < Dim; ++d) {
450 0 : I_nd[d] = I_nd[d] - ldom[d].first() + nghost;
451 : }
452 0 : apply(resultView, I_nd) = apply(view, I_nd);
453 0 : continue;
454 1040 : } else if ((bcType == ZERO_FACE) && (this->isDOFOnBoundary(I_nd))) {
455 0 : continue;
456 : }
457 :
458 : // get the appropriate index for the Kokkos view of the field
459 1752 : for (unsigned d = 0; d < Dim; ++d) {
460 0 : I_nd[d] = I_nd[d] - ldom[d].first() + nghost;
461 : }
462 :
463 3924 : for (j = 0; j < numElementDOFs; ++j) {
464 0 : J_nd = global_dof_ndindices[j];
465 :
466 : // Skip boundary DOFs (Zero & Constant Dirichlet BCs)
467 0 : if (((bcType == ZERO_FACE) || (bcType == CONSTANT_FACE))
468 3480 : && this->isDOFOnBoundary(J_nd)) {
469 0 : continue;
470 : }
471 :
472 : // get the appropriate index for the Kokkos view of the field
473 8040 : for (unsigned d = 0; d < Dim; ++d) {
474 0 : J_nd[d] = J_nd[d] - ldom[d].first() + nghost;
475 : }
476 :
477 2020 : apply(resultView, I_nd) += A_K[i][j] * apply(view, J_nd);
478 : }
479 : }
480 0 : });
481 8 : IpplTimings::stopTimer(outer_loop);
482 :
483 8 : if (bcType == PERIODIC_FACE) {
484 0 : resultField.accumulateHalo();
485 0 : bcField.apply(resultField);
486 0 : bcField.assignGhostToPhysical(resultField);
487 : } else {
488 8 : resultField.accumulateHalo_noghost();
489 : }
490 :
491 8 : IpplTimings::stopTimer(evalAx);
492 :
493 8 : return resultField;
494 8 : }
495 :
496 : template <typename T, unsigned Dim, unsigned Order, typename ElementType,
497 : typename QuadratureType, typename FieldLHS, typename FieldRHS>
498 : template <typename F>
499 0 : FieldLHS LagrangeSpace<T, Dim, Order, ElementType, QuadratureType, FieldLHS,
500 : FieldRHS>::evaluateAx_lower(FieldLHS& field, F& evalFunction) const {
501 0 : Inform m("");
502 :
503 : // start a timer
504 0 : static IpplTimings::TimerRef evalAx = IpplTimings::getTimer("evaluateAx");
505 0 : IpplTimings::startTimer(evalAx);
506 :
507 : // get number of ghost cells in field
508 0 : const int nghost = field.getNghost();
509 :
510 : // create a new field for result with view initialized to zero (views are initialized to
511 : // zero by default)
512 0 : FieldLHS resultField(field.get_mesh(), field.getLayout(), nghost);
513 :
514 : // List of quadrature weights
515 0 : const Vector<T, QuadratureType::numElementNodes> w =
516 0 : this->quadrature_m.getWeightsForRefElement();
517 :
518 : // List of quadrature nodes
519 0 : const Vector<point_t, QuadratureType::numElementNodes> q =
520 0 : this->quadrature_m.getIntegrationNodesForRefElement();
521 :
522 : // TODO move outside of evaluateAx (I think it is possible for other problems as well)
523 : // Gradients of the basis functions for the DOF at the quadrature nodes
524 0 : Vector<Vector<point_t, numElementDOFs>, QuadratureType::numElementNodes> grad_b_q;
525 0 : for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
526 0 : for (size_t i = 0; i < numElementDOFs; ++i) {
527 0 : grad_b_q[k][i] = this->evaluateRefElementShapeFunctionGradient(i, q[k]);
528 : }
529 : }
530 :
531 : // Make local element matrix -- does not change through the element mesh
532 : // Element matrix
533 0 : Vector<Vector<T, numElementDOFs>, numElementDOFs> A_K;
534 :
535 : // 1. Compute the Galerkin element matrix A_K
536 0 : for (size_t i = 0; i < numElementDOFs; ++i) {
537 0 : for (size_t j = 0; j < numElementDOFs; ++j) {
538 0 : A_K[i][j] = 0.0;
539 0 : for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
540 0 : A_K[i][j] += w[k] * evalFunction(i, j, grad_b_q[k]);
541 : }
542 : }
543 : }
544 :
545 : // Get field data and atomic result data,
546 : // since it will be added to during the kokkos loop
547 0 : ViewType view = field.getView();
548 0 : AtomicViewType resultView = resultField.getView();
549 :
550 : // Get boundary conditions from field
551 0 : BConds<FieldLHS, Dim>& bcField = field.getFieldBC();
552 0 : FieldBC bcType = bcField[0]->getBCType();
553 :
554 : // Get domain information
555 0 : auto ldom = (field.getLayout()).getLocalNDIndex();
556 :
557 : using exec_space = typename Kokkos::View<const size_t*>::execution_space;
558 : using policy_type = Kokkos::RangePolicy<exec_space>;
559 :
560 : // start a timer
561 0 : static IpplTimings::TimerRef outer_loop = IpplTimings::getTimer("evaluateAx: outer loop");
562 0 : IpplTimings::startTimer(outer_loop);
563 :
564 : // Loop over elements to compute contributions
565 0 : Kokkos::parallel_for(
566 0 : "Loop over elements", policy_type(0, elementIndices.extent(0)),
567 0 : KOKKOS_CLASS_LAMBDA(const size_t index) {
568 0 : const size_t elementIndex = elementIndices(index);
569 0 : const Vector<size_t, numElementDOFs> local_dof = this->getLocalDOFIndices();
570 0 : const Vector<size_t, numElementDOFs> global_dofs =
571 0 : this->getGlobalDOFIndices(elementIndex);
572 0 : Vector<indices_t, numElementDOFs> global_dof_ndindices;
573 :
574 0 : for (size_t i = 0; i < numElementDOFs; ++i) {
575 0 : global_dof_ndindices[i] = this->getMeshVertexNDIndex(global_dofs[i]);
576 : }
577 :
578 : // local DOF indices
579 : size_t i, j;
580 :
581 : // global DOF n-dimensional indices (Vector of N indices representing indices in
582 : // each dimension)
583 0 : indices_t I_nd, J_nd;
584 :
585 : // 2. Compute the contribution to resultAx = A*x with A_K
586 0 : for (i = 0; i < numElementDOFs; ++i) {
587 0 : I_nd = global_dof_ndindices[i];
588 :
589 : // Handle boundary DOFs
590 : // If Zero Dirichlet BCs, skip this DOF
591 : // If Constant Dirichlet BCs, identity
592 0 : if ((bcType == CONSTANT_FACE) && (this->isDOFOnBoundary(I_nd))) {
593 0 : for (unsigned d = 0; d < Dim; ++d) {
594 0 : I_nd[d] = I_nd[d] - ldom[d].first() + nghost;
595 : }
596 0 : apply(resultView, I_nd) = apply(view, I_nd);
597 0 : continue;
598 0 : } else if ((bcType == ZERO_FACE) && (this->isDOFOnBoundary(I_nd))) {
599 0 : continue;
600 : }
601 :
602 : // get the appropriate index for the Kokkos view of the field
603 0 : for (unsigned d = 0; d < Dim; ++d) {
604 0 : I_nd[d] = I_nd[d] - ldom[d].first() + nghost;
605 : }
606 :
607 0 : for (j = 0; j < numElementDOFs; ++j) {
608 0 : J_nd = global_dof_ndindices[j];
609 :
610 0 : if (global_dofs[i] >= global_dofs[j]) {
611 0 : continue;
612 : }
613 :
614 : // Skip boundary DOFs (Zero & Constant Dirichlet BCs)
615 0 : if (((bcType == ZERO_FACE) || (bcType == CONSTANT_FACE))
616 0 : && this->isDOFOnBoundary(J_nd)) {
617 0 : continue;
618 : }
619 :
620 : // get the appropriate index for the Kokkos view of the field
621 0 : for (unsigned d = 0; d < Dim; ++d) {
622 0 : J_nd[d] = J_nd[d] - ldom[d].first() + nghost;
623 : }
624 :
625 0 : apply(resultView, I_nd) += A_K[i][j] * apply(view, J_nd);
626 : }
627 : }
628 0 : });
629 0 : IpplTimings::stopTimer(outer_loop);
630 :
631 0 : if (bcType == PERIODIC_FACE) {
632 0 : resultField.accumulateHalo();
633 0 : bcField.apply(resultField);
634 0 : bcField.assignGhostToPhysical(resultField);
635 : } else {
636 0 : resultField.accumulateHalo_noghost();
637 : }
638 :
639 0 : IpplTimings::stopTimer(evalAx);
640 :
641 0 : return resultField;
642 0 : }
643 :
644 : template <typename T, unsigned Dim, unsigned Order, typename ElementType,
645 : typename QuadratureType, typename FieldLHS, typename FieldRHS>
646 : template <typename F>
647 0 : FieldLHS LagrangeSpace<T, Dim, Order, ElementType, QuadratureType, FieldLHS,
648 : FieldRHS>::evaluateAx_upper(FieldLHS& field, F& evalFunction) const {
649 0 : Inform m("");
650 :
651 : // start a timer
652 0 : static IpplTimings::TimerRef evalAx = IpplTimings::getTimer("evaluateAx");
653 0 : IpplTimings::startTimer(evalAx);
654 :
655 : // get number of ghost cells in field
656 0 : const int nghost = field.getNghost();
657 :
658 : // create a new field for result with view initialized to zero (views are initialized to
659 : // zero by default)
660 0 : FieldLHS resultField(field.get_mesh(), field.getLayout(), nghost);
661 :
662 : // List of quadrature weights
663 0 : const Vector<T, QuadratureType::numElementNodes> w =
664 0 : this->quadrature_m.getWeightsForRefElement();
665 :
666 : // List of quadrature nodes
667 0 : const Vector<point_t, QuadratureType::numElementNodes> q =
668 0 : this->quadrature_m.getIntegrationNodesForRefElement();
669 :
670 : // TODO move outside of evaluateAx (I think it is possible for other problems as well)
671 : // Gradients of the basis functions for the DOF at the quadrature nodes
672 0 : Vector<Vector<point_t, numElementDOFs>, QuadratureType::numElementNodes> grad_b_q;
673 0 : for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
674 0 : for (size_t i = 0; i < numElementDOFs; ++i) {
675 0 : grad_b_q[k][i] = this->evaluateRefElementShapeFunctionGradient(i, q[k]);
676 : }
677 : }
678 :
679 : // Make local element matrix -- does not change through the element mesh
680 : // Element matrix
681 0 : Vector<Vector<T, numElementDOFs>, numElementDOFs> A_K;
682 :
683 : // 1. Compute the Galerkin element matrix A_K
684 0 : for (size_t i = 0; i < numElementDOFs; ++i) {
685 0 : for (size_t j = 0; j < numElementDOFs; ++j) {
686 0 : A_K[i][j] = 0.0;
687 0 : for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
688 0 : A_K[i][j] += w[k] * evalFunction(i, j, grad_b_q[k]);
689 : }
690 : }
691 : }
692 :
693 : // Get field data and atomic result data,
694 : // since it will be added to during the kokkos loop
695 0 : ViewType view = field.getView();
696 0 : AtomicViewType resultView = resultField.getView();
697 :
698 : // Get boundary conditions from field
699 0 : BConds<FieldLHS, Dim>& bcField = field.getFieldBC();
700 0 : FieldBC bcType = bcField[0]->getBCType();
701 :
702 : // Get domain information
703 0 : auto ldom = (field.getLayout()).getLocalNDIndex();
704 :
705 : using exec_space = typename Kokkos::View<const size_t*>::execution_space;
706 : using policy_type = Kokkos::RangePolicy<exec_space>;
707 :
708 : // start a timer
709 0 : static IpplTimings::TimerRef outer_loop = IpplTimings::getTimer("evaluateAx: outer loop");
710 0 : IpplTimings::startTimer(outer_loop);
711 :
712 : // Loop over elements to compute contributions
713 0 : Kokkos::parallel_for(
714 0 : "Loop over elements", policy_type(0, elementIndices.extent(0)),
715 0 : KOKKOS_CLASS_LAMBDA(const size_t index) {
716 0 : const size_t elementIndex = elementIndices(index);
717 0 : const Vector<size_t, numElementDOFs> local_dof = this->getLocalDOFIndices();
718 0 : const Vector<size_t, numElementDOFs> global_dofs =
719 0 : this->getGlobalDOFIndices(elementIndex);
720 0 : Vector<indices_t, numElementDOFs> global_dof_ndindices;
721 :
722 0 : for (size_t i = 0; i < numElementDOFs; ++i) {
723 0 : global_dof_ndindices[i] = this->getMeshVertexNDIndex(global_dofs[i]);
724 : }
725 :
726 : // local DOF indices
727 : size_t i, j;
728 :
729 : // global DOF n-dimensional indices (Vector of N indices representing indices in
730 : // each dimension)
731 0 : indices_t I_nd, J_nd;
732 :
733 : // 2. Compute the contribution to resultAx = A*x with A_K
734 0 : for (i = 0; i < numElementDOFs; ++i) {
735 0 : I_nd = global_dof_ndindices[i];
736 :
737 : // Handle boundary DOFs
738 : // If Zero Dirichlet BCs, skip this DOF
739 : // If Constant Dirichlet BCs, identity
740 0 : if ((bcType == CONSTANT_FACE) && (this->isDOFOnBoundary(I_nd))) {
741 0 : for (unsigned d = 0; d < Dim; ++d) {
742 0 : I_nd[d] = I_nd[d] - ldom[d].first() + nghost;
743 : }
744 0 : apply(resultView, I_nd) = apply(view, I_nd);
745 0 : continue;
746 0 : } else if ((bcType == ZERO_FACE) && (this->isDOFOnBoundary(I_nd))) {
747 0 : continue;
748 : }
749 :
750 : // get the appropriate index for the Kokkos view of the field
751 0 : for (unsigned d = 0; d < Dim; ++d) {
752 0 : I_nd[d] = I_nd[d] - ldom[d].first() + nghost;
753 : }
754 :
755 0 : for (j = 0; j < numElementDOFs; ++j) {
756 0 : J_nd = global_dof_ndindices[j];
757 :
758 0 : if (global_dofs[i] <= global_dofs[j]) {
759 0 : continue;
760 : }
761 :
762 : // Skip boundary DOFs (Zero & Constant Dirichlet BCs)
763 0 : if (((bcType == ZERO_FACE) || (bcType == CONSTANT_FACE))
764 0 : && this->isDOFOnBoundary(J_nd)) {
765 0 : continue;
766 : }
767 :
768 : // get the appropriate index for the Kokkos view of the field
769 0 : for (unsigned d = 0; d < Dim; ++d) {
770 0 : J_nd[d] = J_nd[d] - ldom[d].first() + nghost;
771 : }
772 :
773 0 : apply(resultView, I_nd) += A_K[i][j] * apply(view, J_nd);
774 : }
775 : }
776 0 : });
777 0 : IpplTimings::stopTimer(outer_loop);
778 :
779 0 : if (bcType == PERIODIC_FACE) {
780 0 : resultField.accumulateHalo();
781 0 : bcField.apply(resultField);
782 0 : bcField.assignGhostToPhysical(resultField);
783 : } else {
784 0 : resultField.accumulateHalo_noghost();
785 : }
786 :
787 0 : IpplTimings::stopTimer(evalAx);
788 :
789 0 : return resultField;
790 0 : }
791 :
792 : template <typename T, unsigned Dim, unsigned Order, typename ElementType,
793 : typename QuadratureType, typename FieldLHS, typename FieldRHS>
794 : template <typename F>
795 0 : FieldLHS LagrangeSpace<T, Dim, Order, ElementType, QuadratureType, FieldLHS,
796 : FieldRHS>::evaluateAx_upperlower(FieldLHS& field, F& evalFunction) const {
797 0 : Inform m("");
798 :
799 : // start a timer
800 0 : static IpplTimings::TimerRef evalAx = IpplTimings::getTimer("evaluateAx");
801 0 : IpplTimings::startTimer(evalAx);
802 :
803 : // get number of ghost cells in field
804 0 : const int nghost = field.getNghost();
805 :
806 : // create a new field for result with view initialized to zero (views are initialized to
807 : // zero by default)
808 0 : FieldLHS resultField(field.get_mesh(), field.getLayout(), nghost);
809 :
810 : // List of quadrature weights
811 0 : const Vector<T, QuadratureType::numElementNodes> w =
812 0 : this->quadrature_m.getWeightsForRefElement();
813 :
814 : // List of quadrature nodes
815 0 : const Vector<point_t, QuadratureType::numElementNodes> q =
816 0 : this->quadrature_m.getIntegrationNodesForRefElement();
817 :
818 : // TODO move outside of evaluateAx (I think it is possible for other problems as well)
819 : // Gradients of the basis functions for the DOF at the quadrature nodes
820 0 : Vector<Vector<point_t, numElementDOFs>, QuadratureType::numElementNodes> grad_b_q;
821 0 : for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
822 0 : for (size_t i = 0; i < numElementDOFs; ++i) {
823 0 : grad_b_q[k][i] = this->evaluateRefElementShapeFunctionGradient(i, q[k]);
824 : }
825 : }
826 :
827 : // Make local element matrix -- does not change through the element mesh
828 : // Element matrix
829 0 : Vector<Vector<T, numElementDOFs>, numElementDOFs> A_K;
830 :
831 : // 1. Compute the Galerkin element matrix A_K
832 0 : for (size_t i = 0; i < numElementDOFs; ++i) {
833 0 : for (size_t j = 0; j < numElementDOFs; ++j) {
834 0 : A_K[i][j] = 0.0;
835 0 : for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
836 0 : A_K[i][j] += w[k] * evalFunction(i, j, grad_b_q[k]);
837 : }
838 : }
839 : }
840 :
841 : // Get field data and atomic result data,
842 : // since it will be added to during the kokkos loop
843 0 : ViewType view = field.getView();
844 0 : AtomicViewType resultView = resultField.getView();
845 :
846 : // Get boundary conditions from field
847 0 : BConds<FieldLHS, Dim>& bcField = field.getFieldBC();
848 0 : FieldBC bcType = bcField[0]->getBCType();
849 :
850 : // Get domain information
851 0 : auto ldom = (field.getLayout()).getLocalNDIndex();
852 :
853 : using exec_space = typename Kokkos::View<const size_t*>::execution_space;
854 : using policy_type = Kokkos::RangePolicy<exec_space>;
855 :
856 : // start a timer
857 0 : static IpplTimings::TimerRef outer_loop = IpplTimings::getTimer("evaluateAx: outer loop");
858 0 : IpplTimings::startTimer(outer_loop);
859 :
860 : // Loop over elements to compute contributions
861 0 : Kokkos::parallel_for(
862 0 : "Loop over elements", policy_type(0, elementIndices.extent(0)),
863 0 : KOKKOS_CLASS_LAMBDA(const size_t index) {
864 0 : const size_t elementIndex = elementIndices(index);
865 0 : const Vector<size_t, numElementDOFs> local_dof = this->getLocalDOFIndices();
866 0 : const Vector<size_t, numElementDOFs> global_dofs =
867 0 : this->getGlobalDOFIndices(elementIndex);
868 0 : Vector<indices_t, numElementDOFs> global_dof_ndindices;
869 :
870 0 : for (size_t i = 0; i < numElementDOFs; ++i) {
871 0 : global_dof_ndindices[i] = this->getMeshVertexNDIndex(global_dofs[i]);
872 : }
873 :
874 : // local DOF indices
875 : size_t i, j;
876 :
877 : // global DOF n-dimensional indices (Vector of N indices representing indices in
878 : // each dimension)
879 0 : indices_t I_nd, J_nd;
880 :
881 : // 2. Compute the contribution to resultAx = A*x with A_K
882 0 : for (i = 0; i < numElementDOFs; ++i) {
883 0 : I_nd = global_dof_ndindices[i];
884 :
885 : // Handle boundary DOFs
886 : // If Zero Dirichlet BCs, skip this DOF
887 : // If Constant Dirichlet BCs, identity
888 0 : if ((bcType == CONSTANT_FACE) && (this->isDOFOnBoundary(I_nd))) {
889 0 : for (unsigned d = 0; d < Dim; ++d) {
890 0 : I_nd[d] = I_nd[d] - ldom[d].first() + nghost;
891 : }
892 0 : apply(resultView, I_nd) = apply(view, I_nd);
893 0 : continue;
894 0 : } else if ((bcType == ZERO_FACE) && (this->isDOFOnBoundary(I_nd))) {
895 0 : continue;
896 : }
897 :
898 : // get the appropriate index for the Kokkos view of the field
899 0 : for (unsigned d = 0; d < Dim; ++d) {
900 0 : I_nd[d] = I_nd[d] - ldom[d].first() + nghost;
901 : }
902 :
903 0 : for (j = 0; j < numElementDOFs; ++j) {
904 0 : J_nd = global_dof_ndindices[j];
905 :
906 0 : if (global_dofs[i] == global_dofs[j]) {
907 0 : continue;
908 : }
909 :
910 : // Skip boundary DOFs (Zero & Constant Dirichlet BCs)
911 0 : if (((bcType == ZERO_FACE) || (bcType == CONSTANT_FACE))
912 0 : && this->isDOFOnBoundary(J_nd)) {
913 0 : continue;
914 : }
915 :
916 : // get the appropriate index for the Kokkos view of the field
917 0 : for (unsigned d = 0; d < Dim; ++d) {
918 0 : J_nd[d] = J_nd[d] - ldom[d].first() + nghost;
919 : }
920 :
921 0 : apply(resultView, I_nd) += A_K[i][j] * apply(view, J_nd);
922 : }
923 : }
924 0 : });
925 0 : IpplTimings::stopTimer(outer_loop);
926 :
927 0 : if (bcType == PERIODIC_FACE) {
928 0 : resultField.accumulateHalo();
929 0 : bcField.apply(resultField);
930 0 : bcField.assignGhostToPhysical(resultField);
931 : } else {
932 0 : resultField.accumulateHalo_noghost();
933 : }
934 :
935 0 : IpplTimings::stopTimer(evalAx);
936 :
937 0 : return resultField;
938 0 : }
939 :
940 : template <typename T, unsigned Dim, unsigned Order, typename ElementType,
941 : typename QuadratureType, typename FieldLHS, typename FieldRHS>
942 : template <typename F>
943 0 : FieldLHS LagrangeSpace<T, Dim, Order, ElementType, QuadratureType, FieldLHS,
944 : FieldRHS>::evaluateAx_inversediag(FieldLHS& field, F& evalFunction) const {
945 0 : Inform m("");
946 :
947 : // start a timer
948 0 : static IpplTimings::TimerRef evalAx = IpplTimings::getTimer("evaluateAx");
949 0 : IpplTimings::startTimer(evalAx);
950 :
951 : // get number of ghost cells in field
952 0 : const int nghost = field.getNghost();
953 :
954 : // create a new field for result with view initialized to zero (views are initialized to
955 : // zero by default)
956 0 : FieldLHS resultField(field.get_mesh(), field.getLayout(), nghost);
957 :
958 : // List of quadrature weights
959 0 : const Vector<T, QuadratureType::numElementNodes> w =
960 0 : this->quadrature_m.getWeightsForRefElement();
961 :
962 : // List of quadrature nodes
963 0 : const Vector<point_t, QuadratureType::numElementNodes> q =
964 0 : this->quadrature_m.getIntegrationNodesForRefElement();
965 :
966 : // TODO move outside of evaluateAx (I think it is possible for other problems as well)
967 : // Gradients of the basis functions for the DOF at the quadrature nodes
968 0 : Vector<Vector<point_t, numElementDOFs>, QuadratureType::numElementNodes> grad_b_q;
969 0 : for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
970 0 : for (size_t i = 0; i < numElementDOFs; ++i) {
971 0 : grad_b_q[k][i] = this->evaluateRefElementShapeFunctionGradient(i, q[k]);
972 : }
973 : }
974 :
975 : // Make local element matrix -- does not change through the element mesh
976 : // Element matrix
977 0 : Vector<Vector<T, numElementDOFs>, numElementDOFs> A_K;
978 :
979 : // 1. Compute the Galerkin element matrix A_K
980 0 : for (size_t i = 0; i < numElementDOFs; ++i) {
981 0 : for (size_t j = 0; j < numElementDOFs; ++j) {
982 0 : A_K[i][j] = 0.0;
983 0 : for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
984 0 : A_K[i][j] += w[k] * evalFunction(i, j, grad_b_q[k]);
985 : }
986 : }
987 : }
988 :
989 : // Get field data and atomic result data,
990 : // since it will be added to during the kokkos loop
991 0 : ViewType view = field.getView();
992 0 : AtomicViewType resultView = resultField.getView();
993 :
994 : // Get boundary conditions from field
995 0 : BConds<FieldLHS, Dim>& bcField = field.getFieldBC();
996 0 : FieldBC bcType = bcField[0]->getBCType();
997 :
998 : // Get domain information
999 0 : auto ldom = (field.getLayout()).getLocalNDIndex();
1000 :
1001 : using exec_space = typename Kokkos::View<const size_t*>::execution_space;
1002 : using policy_type = Kokkos::RangePolicy<exec_space>;
1003 :
1004 : // start a timer
1005 0 : static IpplTimings::TimerRef outer_loop = IpplTimings::getTimer("evaluateAx: outer loop");
1006 0 : IpplTimings::startTimer(outer_loop);
1007 :
1008 : // Loop over elements to compute contributions
1009 0 : Kokkos::parallel_for(
1010 0 : "Loop over elements", policy_type(0, elementIndices.extent(0)),
1011 0 : KOKKOS_CLASS_LAMBDA(const size_t index) {
1012 0 : const size_t elementIndex = elementIndices(index);
1013 0 : const Vector<size_t, numElementDOFs> local_dof = this->getLocalDOFIndices();
1014 0 : const Vector<size_t, numElementDOFs> global_dofs =
1015 0 : this->getGlobalDOFIndices(elementIndex);
1016 0 : Vector<indices_t, numElementDOFs> global_dof_ndindices;
1017 :
1018 0 : for (size_t i = 0; i < numElementDOFs; ++i) {
1019 0 : global_dof_ndindices[i] = this->getMeshVertexNDIndex(global_dofs[i]);
1020 : }
1021 :
1022 : // local DOF indices
1023 : size_t i, j;
1024 :
1025 : // global DOF n-dimensional indices (Vector of N indices representing indices in
1026 : // each dimension)
1027 0 : indices_t I_nd, J_nd;
1028 :
1029 : // 2. Compute the contribution to resultAx = A*x with A_K
1030 0 : for (i = 0; i < numElementDOFs; ++i) {
1031 0 : I_nd = global_dof_ndindices[i];
1032 :
1033 : // Handle boundary DOFs
1034 : // If Zero Dirichlet BCs, skip this DOF
1035 : // If Constant Dirichlet BCs, identity
1036 0 : if ((bcType == CONSTANT_FACE) && (this->isDOFOnBoundary(I_nd))) {
1037 0 : for (unsigned d = 0; d < Dim; ++d) {
1038 0 : I_nd[d] = I_nd[d] - ldom[d].first() + nghost;
1039 : }
1040 0 : apply(resultView, I_nd) = 1.0;
1041 0 : continue;
1042 0 : } else if ((bcType == ZERO_FACE) && (this->isDOFOnBoundary(I_nd))) {
1043 0 : continue;
1044 : }
1045 :
1046 : // get the appropriate index for the Kokkos view of the field
1047 0 : for (unsigned d = 0; d < Dim; ++d) {
1048 0 : I_nd[d] = I_nd[d] - ldom[d].first() + nghost;
1049 : }
1050 :
1051 0 : for (j = 0; j < numElementDOFs; ++j) {
1052 0 : if (global_dofs[i] == global_dofs[j]) {
1053 0 : J_nd = global_dof_ndindices[j];
1054 :
1055 : // Skip boundary DOFs (Zero & Constant Dirichlet BCs)
1056 0 : if (((bcType == ZERO_FACE) || (bcType == CONSTANT_FACE))
1057 0 : && this->isDOFOnBoundary(J_nd)) {
1058 0 : continue;
1059 : }
1060 :
1061 : // get the appropriate index for the Kokkos view of the field
1062 0 : for (unsigned d = 0; d < Dim; ++d) {
1063 0 : J_nd[d] = J_nd[d] - ldom[d].first() + nghost;
1064 : }
1065 :
1066 : // sum up all contributions of element matrix
1067 0 : apply(resultView, I_nd) += A_K[i][j];
1068 : }
1069 : }
1070 : }
1071 0 : });
1072 :
1073 0 : if (bcType == PERIODIC_FACE) {
1074 0 : resultField.accumulateHalo();
1075 0 : bcField.apply(resultField);
1076 0 : bcField.assignGhostToPhysical(resultField);
1077 : } else {
1078 0 : resultField.accumulateHalo_noghost();
1079 : }
1080 :
1081 : // apply the inverse diagonal after already summed all contributions from element matrices
1082 : using index_array_type = typename RangePolicy<Dim, exec_space>::index_array_type;
1083 0 : ippl::parallel_for("Loop over result view to apply inverse", field.getFieldRangePolicy(),
1084 0 : KOKKOS_LAMBDA(const index_array_type& args) {
1085 0 : if (apply(resultView, args) != 0.0) {
1086 0 : apply(resultView, args) = (1.0 / apply(resultView, args)) * apply(view, args);
1087 : }
1088 : });
1089 0 : IpplTimings::stopTimer(outer_loop);
1090 :
1091 0 : IpplTimings::stopTimer(evalAx);
1092 :
1093 0 : return resultField;
1094 0 : }
1095 :
1096 : template <typename T, unsigned Dim, unsigned Order, typename ElementType,
1097 : typename QuadratureType, typename FieldLHS, typename FieldRHS>
1098 : template <typename F>
1099 0 : FieldLHS LagrangeSpace<T, Dim, Order, ElementType, QuadratureType, FieldLHS,
1100 : FieldRHS>::evaluateAx_diag(FieldLHS& field, F& evalFunction) const {
1101 0 : Inform m("");
1102 :
1103 : // start a timer
1104 0 : static IpplTimings::TimerRef evalAx = IpplTimings::getTimer("evaluateAx");
1105 0 : IpplTimings::startTimer(evalAx);
1106 :
1107 : // get number of ghost cells in field
1108 0 : const int nghost = field.getNghost();
1109 :
1110 : // create a new field for result with view initialized to zero (views are initialized to
1111 : // zero by default)
1112 0 : FieldLHS resultField(field.get_mesh(), field.getLayout(), nghost);
1113 :
1114 : // List of quadrature weights
1115 0 : const Vector<T, QuadratureType::numElementNodes> w =
1116 0 : this->quadrature_m.getWeightsForRefElement();
1117 :
1118 : // List of quadrature nodes
1119 0 : const Vector<point_t, QuadratureType::numElementNodes> q =
1120 0 : this->quadrature_m.getIntegrationNodesForRefElement();
1121 :
1122 : // TODO move outside of evaluateAx (I think it is possible for other problems as well)
1123 : // Gradients of the basis functions for the DOF at the quadrature nodes
1124 0 : Vector<Vector<point_t, numElementDOFs>, QuadratureType::numElementNodes> grad_b_q;
1125 0 : for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
1126 0 : for (size_t i = 0; i < numElementDOFs; ++i) {
1127 0 : grad_b_q[k][i] = this->evaluateRefElementShapeFunctionGradient(i, q[k]);
1128 : }
1129 : }
1130 :
1131 : // Make local element matrix -- does not change through the element mesh
1132 : // Element matrix
1133 0 : Vector<Vector<T, numElementDOFs>, numElementDOFs> A_K;
1134 :
1135 : // 1. Compute the Galerkin element matrix A_K
1136 0 : for (size_t i = 0; i < numElementDOFs; ++i) {
1137 0 : for (size_t j = 0; j < numElementDOFs; ++j) {
1138 0 : A_K[i][j] = 0.0;
1139 0 : for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
1140 0 : A_K[i][j] += w[k] * evalFunction(i, j, grad_b_q[k]);
1141 : }
1142 : }
1143 : }
1144 :
1145 : // Get field data and atomic result data,
1146 : // since it will be added to during the kokkos loop
1147 0 : ViewType view = field.getView();
1148 0 : AtomicViewType resultView = resultField.getView();
1149 :
1150 : // Get boundary conditions from field
1151 0 : BConds<FieldLHS, Dim>& bcField = field.getFieldBC();
1152 0 : FieldBC bcType = bcField[0]->getBCType();
1153 :
1154 : // Get domain information
1155 0 : auto ldom = (field.getLayout()).getLocalNDIndex();
1156 :
1157 : using exec_space = typename Kokkos::View<const size_t*>::execution_space;
1158 : using policy_type = Kokkos::RangePolicy<exec_space>;
1159 :
1160 : // start a timer
1161 0 : static IpplTimings::TimerRef outer_loop = IpplTimings::getTimer("evaluateAx: outer loop");
1162 0 : IpplTimings::startTimer(outer_loop);
1163 :
1164 : // Loop over elements to compute contributions
1165 0 : Kokkos::parallel_for(
1166 0 : "Loop over elements", policy_type(0, elementIndices.extent(0)),
1167 0 : KOKKOS_CLASS_LAMBDA(const size_t index) {
1168 0 : const size_t elementIndex = elementIndices(index);
1169 0 : const Vector<size_t, numElementDOFs> local_dof = this->getLocalDOFIndices();
1170 0 : const Vector<size_t, numElementDOFs> global_dofs =
1171 0 : this->getGlobalDOFIndices(elementIndex);
1172 0 : Vector<indices_t, numElementDOFs> global_dof_ndindices;
1173 :
1174 0 : for (size_t i = 0; i < numElementDOFs; ++i) {
1175 0 : global_dof_ndindices[i] = this->getMeshVertexNDIndex(global_dofs[i]);
1176 : }
1177 :
1178 : // local DOF indices
1179 : size_t i, j;
1180 :
1181 : // global DOF n-dimensional indices (Vector of N indices representing indices in
1182 : // each dimension)
1183 0 : indices_t I_nd, J_nd;
1184 :
1185 : // 2. Compute the contribution to resultAx = A*x with A_K
1186 0 : for (i = 0; i < numElementDOFs; ++i) {
1187 0 : I_nd = global_dof_ndindices[i];
1188 :
1189 : // Handle boundary DOFs
1190 : // If Zero Dirichlet BCs, skip this DOF
1191 : // If Constant Dirichlet BCs, identity
1192 0 : if ((bcType == CONSTANT_FACE) && (this->isDOFOnBoundary(I_nd))) {
1193 0 : for (unsigned d = 0; d < Dim; ++d) {
1194 0 : I_nd[d] = I_nd[d] - ldom[d].first() + nghost;
1195 : }
1196 0 : apply(resultView, I_nd) = apply(view, I_nd);
1197 0 : continue;
1198 0 : } else if ((bcType == ZERO_FACE) && (this->isDOFOnBoundary(I_nd))) {
1199 0 : continue;
1200 : }
1201 :
1202 : // get the appropriate index for the Kokkos view of the field
1203 0 : for (unsigned d = 0; d < Dim; ++d) {
1204 0 : I_nd[d] = I_nd[d] - ldom[d].first() + nghost;
1205 : }
1206 :
1207 0 : for (j = 0; j < numElementDOFs; ++j) {
1208 0 : if (global_dofs[i] == global_dofs[j]) {
1209 0 : J_nd = global_dof_ndindices[j];
1210 :
1211 : // Skip boundary DOFs (Zero & Constant Dirichlet BCs)
1212 0 : if (((bcType == ZERO_FACE) || (bcType == CONSTANT_FACE))
1213 0 : && this->isDOFOnBoundary(J_nd)) {
1214 0 : continue;
1215 : }
1216 :
1217 : // get the appropriate index for the Kokkos view of the field
1218 0 : for (unsigned d = 0; d < Dim; ++d) {
1219 0 : J_nd[d] = J_nd[d] - ldom[d].first() + nghost;
1220 : }
1221 :
1222 0 : apply(resultView, I_nd) += A_K[i][j] * apply(view, J_nd);
1223 : }
1224 : }
1225 : }
1226 0 : });
1227 0 : IpplTimings::stopTimer(outer_loop);
1228 :
1229 0 : if (bcType == PERIODIC_FACE) {
1230 0 : resultField.accumulateHalo();
1231 0 : bcField.apply(resultField);
1232 0 : bcField.assignGhostToPhysical(resultField);
1233 : } else {
1234 0 : resultField.accumulateHalo_noghost();
1235 : }
1236 :
1237 0 : IpplTimings::stopTimer(evalAx);
1238 :
1239 0 : return resultField;
1240 0 : }
1241 :
1242 : template <typename T, unsigned Dim, unsigned Order, typename ElementType,
1243 : typename QuadratureType, typename FieldLHS, typename FieldRHS>
1244 : template <typename F>
1245 0 : FieldLHS LagrangeSpace<T, Dim, Order, ElementType, QuadratureType, FieldLHS,
1246 : FieldRHS>::evaluateAx_lift(FieldLHS& field, F& evalFunction) const {
1247 0 : Inform m("");
1248 :
1249 : // get number of ghost cells in field
1250 0 : const int nghost = field.getNghost();
1251 :
1252 : // create a new field for result with view initialized to zero (views are initialized to
1253 : // zero by default)
1254 0 : FieldLHS resultField(field.get_mesh(), field.getLayout(), nghost);
1255 :
1256 : // List of quadrature weights
1257 0 : const Vector<T, QuadratureType::numElementNodes> w =
1258 0 : this->quadrature_m.getWeightsForRefElement();
1259 :
1260 : // List of quadrature nodes
1261 0 : const Vector<point_t, QuadratureType::numElementNodes> q =
1262 0 : this->quadrature_m.getIntegrationNodesForRefElement();
1263 :
1264 : // TODO move outside of evaluateAx (I think it is possible for other problems as well)
1265 : // Gradients of the basis functions for the DOF at the quadrature nodes
1266 0 : Vector<Vector<point_t, numElementDOFs>, QuadratureType::numElementNodes> grad_b_q;
1267 0 : for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
1268 0 : for (size_t i = 0; i < numElementDOFs; ++i) {
1269 0 : grad_b_q[k][i] = this->evaluateRefElementShapeFunctionGradient(i, q[k]);
1270 : }
1271 : }
1272 :
1273 : // Make local element matrix -- does not change through the element mesh
1274 : // Element matrix
1275 0 : Vector<Vector<T, numElementDOFs>, numElementDOFs> A_K;
1276 :
1277 : // 1. Compute the Galerkin element matrix A_K
1278 0 : for (size_t i = 0; i < numElementDOFs; ++i) {
1279 0 : for (size_t j = 0; j < numElementDOFs; ++j) {
1280 0 : A_K[i][j] = 0.0;
1281 0 : for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
1282 0 : A_K[i][j] += w[k] * evalFunction(i, j, grad_b_q[k]);
1283 : }
1284 : }
1285 : }
1286 :
1287 : // Get field data and atomic result data,
1288 : // since it will be added to during the kokkos loop
1289 0 : ViewType view = field.getView();
1290 0 : AtomicViewType resultView = resultField.getView();
1291 :
1292 : // Get domain information
1293 0 : auto ldom = (field.getLayout()).getLocalNDIndex();
1294 :
1295 : using exec_space = typename Kokkos::View<const size_t*>::execution_space;
1296 : using policy_type = Kokkos::RangePolicy<exec_space>;
1297 :
1298 : // Loop over elements to compute contributions
1299 0 : Kokkos::parallel_for(
1300 0 : "Loop over elements", policy_type(0, elementIndices.extent(0)),
1301 0 : KOKKOS_CLASS_LAMBDA(const size_t index) {
1302 0 : const size_t elementIndex = elementIndices(index);
1303 0 : const Vector<size_t, numElementDOFs> local_dof = this->getLocalDOFIndices();
1304 0 : const Vector<size_t, numElementDOFs> global_dofs =
1305 0 : this->getGlobalDOFIndices(elementIndex);
1306 0 : Vector<indices_t, numElementDOFs> global_dof_ndindices;
1307 :
1308 0 : for (size_t i = 0; i < numElementDOFs; ++i) {
1309 0 : global_dof_ndindices[i] = this->getMeshVertexNDIndex(global_dofs[i]);
1310 : }
1311 :
1312 : // local DOF indices (both i and j go from 0 to numDOFs-1 in the element)
1313 : size_t i, j;
1314 :
1315 : // global DOF n-dimensional indices (Vector of N indices representing indices in
1316 : // each dimension)
1317 0 : indices_t I_nd, J_nd;
1318 :
1319 : // 2. Compute the contribution to resultAx = A*x with A_K
1320 0 : for (i = 0; i < numElementDOFs; ++i) {
1321 0 : I_nd = global_dof_ndindices[i];
1322 :
1323 : // Skip if on a row of the matrix
1324 0 : if (this->isDOFOnBoundary(I_nd)) {
1325 0 : continue;
1326 : }
1327 :
1328 : // get the appropriate index for the Kokkos view of the field
1329 0 : for (unsigned d = 0; d < Dim; ++d) {
1330 0 : I_nd[d] = I_nd[d] - ldom[d].first() + nghost;
1331 : }
1332 :
1333 0 : for (j = 0; j < numElementDOFs; ++j) {
1334 0 : J_nd = global_dof_ndindices[j];
1335 :
1336 : // Contribute to lifting only if on a boundary DOF
1337 0 : if (this->isDOFOnBoundary(J_nd)) {
1338 : // get the appropriate index for the Kokkos view of the field
1339 0 : for (unsigned d = 0; d < Dim; ++d) {
1340 0 : J_nd[d] = J_nd[d] - ldom[d].first() + nghost;
1341 : }
1342 0 : apply(resultView, I_nd) += A_K[i][j] * apply(view, J_nd);
1343 0 : continue;
1344 0 : }
1345 :
1346 : }
1347 : }
1348 0 : });
1349 0 : resultField.accumulateHalo();
1350 :
1351 0 : return resultField;
1352 0 : }
1353 :
1354 : template <typename T, unsigned Dim, unsigned Order, typename ElementType,
1355 : typename QuadratureType, typename FieldLHS, typename FieldRHS>
1356 20 : void LagrangeSpace<T, Dim, Order, ElementType, QuadratureType, FieldLHS,
1357 : FieldRHS>::evaluateLoadVector(FieldRHS& field) const {
1358 20 : Inform m("");
1359 :
1360 : // start a timer
1361 20 : static IpplTimings::TimerRef evalLoadV = IpplTimings::getTimer("evaluateLoadVector");
1362 20 : IpplTimings::startTimer(evalLoadV);
1363 :
1364 : // List of quadrature weights
1365 20 : const Vector<T, QuadratureType::numElementNodes> w =
1366 20 : this->quadrature_m.getWeightsForRefElement();
1367 :
1368 : // List of quadrature nodes
1369 20 : const Vector<point_t, QuadratureType::numElementNodes> q =
1370 20 : this->quadrature_m.getIntegrationNodesForRefElement();
1371 :
1372 20 : const indices_t zeroNdIndex = Vector<size_t, Dim>(0);
1373 :
1374 : // Evaluate the basis functions for the DOF at the quadrature nodes
1375 20 : Vector<Vector<T, numElementDOFs>, QuadratureType::numElementNodes> basis_q;
1376 648 : for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
1377 5108 : for (size_t i = 0; i < numElementDOFs; ++i) {
1378 4480 : basis_q[k][i] = this->evaluateRefElementShapeFunction(i, q[k]);
1379 : }
1380 : }
1381 :
1382 : // Absolute value of det Phi_K
1383 20 : const T absDetDPhi = Kokkos::abs(this->ref_element_m.getDeterminantOfTransformationJacobian(
1384 40 : this->getElementMeshVertexPoints(zeroNdIndex)));
1385 :
1386 : // Get domain information and ghost cells
1387 20 : auto ldom = (field.getLayout()).getLocalNDIndex();
1388 20 : const int nghost = field.getNghost();
1389 :
1390 : // Get boundary conditions from field
1391 20 : BConds<FieldRHS, Dim>& bcField = field.getFieldBC();
1392 20 : FieldBC bcType = bcField[0]->getBCType();
1393 :
1394 20 : FieldRHS temp_field(field.get_mesh(), field.getLayout(), nghost);
1395 20 : temp_field.setFieldBC(bcField);
1396 :
1397 : // Get field data and make it atomic,
1398 : // since it will be added to during the kokkos loop
1399 : // We work with a temporary field since we need to use field
1400 : // to evaluate the load vector; then we assign temp to RHS field
1401 20 : AtomicViewType atomic_view = temp_field.getView();
1402 :
1403 : using exec_space = typename Kokkos::View<const size_t*>::execution_space;
1404 : using policy_type = Kokkos::RangePolicy<exec_space>;
1405 :
1406 : // start a timer
1407 20 : static IpplTimings::TimerRef outer_loop =
1408 20 : IpplTimings::getTimer("evaluateLoadVec: outer loop");
1409 20 : IpplTimings::startTimer(outer_loop);
1410 :
1411 : // Loop over elements to compute contributions
1412 20 : Kokkos::parallel_for(
1413 40 : "Loop over elements", policy_type(0, elementIndices.extent(0)),
1414 344 : KOKKOS_CLASS_LAMBDA(size_t index) {
1415 0 : const size_t elementIndex = elementIndices(index);
1416 304 : const Vector<size_t, numElementDOFs> local_dofs = this->getLocalDOFIndices();
1417 0 : const Vector<size_t, numElementDOFs> global_dofs =
1418 304 : this->getGlobalDOFIndices(elementIndex);
1419 :
1420 : size_t i, I;
1421 :
1422 : // 1. Compute b_K
1423 4720 : for (i = 0; i < numElementDOFs; ++i) {
1424 0 : I = global_dofs[i];
1425 :
1426 : // TODO fix for higher order
1427 2208 : auto dof_ndindex_I = this->getMeshVertexNDIndex(I);
1428 :
1429 : // Skip boundary DOFs (Zero and Constant Dirichlet BCs)
1430 0 : if (((bcType == ZERO_FACE) || (bcType == CONSTANT_FACE))
1431 2208 : && (this->isDOFOnBoundary(dof_ndindex_I))) {
1432 0 : continue;
1433 : }
1434 :
1435 : // calculate the contribution of this element
1436 0 : T contrib = 0;
1437 57264 : for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
1438 0 : T val = 0;
1439 499104 : for (size_t j = 0; j < numElementDOFs; ++j) {
1440 : // get field index corresponding to this DOF
1441 0 : size_t J = global_dofs[j];
1442 442800 : auto dof_ndindex_J = this->getMeshVertexNDIndex(J);
1443 1763712 : for (unsigned d = 0; d < Dim; ++d) {
1444 0 : dof_ndindex_J[d] = dof_ndindex_J[d] - ldom[d].first() + nghost;
1445 : }
1446 :
1447 : // get field value at DOF and interpolate to q_k
1448 442800 : val += basis_q[k][j] * apply(field, dof_ndindex_J);
1449 : }
1450 :
1451 0 : contrib += w[k] * basis_q[k][i] * absDetDPhi * val;
1452 : }
1453 :
1454 : // get the appropriate index for the Kokkos view of the field
1455 3720 : for (unsigned d = 0; d < Dim; ++d) {
1456 0 : dof_ndindex_I[d] = dof_ndindex_I[d] - ldom[d].first() + nghost;
1457 : }
1458 :
1459 : // add the contribution of the element to the field
1460 960 : apply(atomic_view, dof_ndindex_I) += contrib;
1461 :
1462 : }
1463 0 : });
1464 20 : IpplTimings::stopTimer(outer_loop);
1465 :
1466 20 : temp_field.accumulateHalo();
1467 :
1468 20 : if ((bcType == PERIODIC_FACE) || (bcType == CONSTANT_FACE)) {
1469 0 : bcField.apply(temp_field);
1470 0 : bcField.assignGhostToPhysical(temp_field);
1471 : }
1472 :
1473 20 : field = temp_field;
1474 :
1475 20 : IpplTimings::stopTimer(evalLoadV);
1476 20 : }
1477 :
1478 : template <typename T, unsigned Dim, unsigned Order, typename ElementType,
1479 : typename QuadratureType, typename FieldLHS, typename FieldRHS>
1480 : template <typename F>
1481 0 : T LagrangeSpace<T, Dim, Order, ElementType, QuadratureType, FieldLHS, FieldRHS>::computeErrorL2(
1482 : const FieldLHS& u_h, const F& u_sol) const {
1483 0 : if (this->quadrature_m.getOrder() < (2 * Order + 1)) {
1484 : // throw exception
1485 0 : throw IpplException(
1486 : "LagrangeSpace::computeErrorL2()",
1487 : "Order of quadrature rule for error computation should be > 2*p + 1");
1488 : }
1489 :
1490 : // List of quadrature weights
1491 0 : const Vector<T, QuadratureType::numElementNodes> w =
1492 0 : this->quadrature_m.getWeightsForRefElement();
1493 :
1494 : // List of quadrature nodes
1495 0 : const Vector<point_t, QuadratureType::numElementNodes> q =
1496 0 : this->quadrature_m.getIntegrationNodesForRefElement();
1497 :
1498 : // Evaluate the basis functions for the DOF at the quadrature nodes
1499 0 : Vector<Vector<T, numElementDOFs>, QuadratureType::numElementNodes> basis_q;
1500 0 : for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
1501 0 : for (size_t i = 0; i < numElementDOFs; ++i) {
1502 0 : basis_q[k][i] = this->evaluateRefElementShapeFunction(i, q[k]);
1503 : }
1504 : }
1505 :
1506 0 : const indices_t zeroNdIndex = Vector<size_t, Dim>(0);
1507 :
1508 : // Absolute value of det Phi_K
1509 0 : const T absDetDPhi = Kokkos::abs(this->ref_element_m.getDeterminantOfTransformationJacobian(
1510 0 : this->getElementMeshVertexPoints(zeroNdIndex)));
1511 :
1512 : // Variable to sum the error to
1513 0 : T error = 0;
1514 :
1515 : // Get domain information and ghost cells
1516 0 : auto ldom = (u_h.getLayout()).getLocalNDIndex();
1517 0 : const int nghost = u_h.getNghost();
1518 :
1519 : using exec_space = typename Kokkos::View<const size_t*>::execution_space;
1520 : using policy_type = Kokkos::RangePolicy<exec_space>;
1521 :
1522 : // Loop over elements to compute contributions
1523 0 : Kokkos::parallel_reduce(
1524 0 : "Compute error over elements", policy_type(0, elementIndices.extent(0)),
1525 0 : KOKKOS_CLASS_LAMBDA(size_t index, double& local) {
1526 0 : const size_t elementIndex = elementIndices(index);
1527 0 : const Vector<size_t, numElementDOFs> global_dofs =
1528 0 : this->getGlobalDOFIndices(elementIndex);
1529 :
1530 : // contribution of this element to the error
1531 0 : T contrib = 0;
1532 0 : for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
1533 0 : T val_u_sol = u_sol(this->ref_element_m.localToGlobal(
1534 0 : this->getElementMeshVertexPoints(this->getElementNDIndex(elementIndex)),
1535 0 : q[k]));
1536 :
1537 0 : T val_u_h = 0;
1538 0 : for (size_t i = 0; i < numElementDOFs; ++i) {
1539 : // get field index corresponding to this DOF
1540 0 : size_t I = global_dofs[i];
1541 0 : auto dof_ndindex_I = this->getMeshVertexNDIndex(I);
1542 0 : for (unsigned d = 0; d < Dim; ++d) {
1543 0 : dof_ndindex_I[d] = dof_ndindex_I[d] - ldom[d].first() + nghost;
1544 : }
1545 :
1546 : // get field value at DOF and interpolate to q_k
1547 0 : val_u_h += basis_q[k][i] * apply(u_h, dof_ndindex_I);
1548 : }
1549 :
1550 0 : contrib += w[k] * Kokkos::pow(val_u_sol - val_u_h, 2) * absDetDPhi;
1551 : }
1552 0 : local += contrib;
1553 0 : },
1554 0 : Kokkos::Sum<double>(error));
1555 :
1556 : // MPI reduce
1557 0 : T global_error = 0.0;
1558 0 : Comm->allreduce(error, global_error, 1, std::plus<T>());
1559 :
1560 0 : return Kokkos::sqrt(global_error);
1561 0 : }
1562 :
1563 : template <typename T, unsigned Dim, unsigned Order, typename ElementType,
1564 : typename QuadratureType, typename FieldLHS, typename FieldRHS>
1565 0 : T LagrangeSpace<T, Dim, Order, ElementType, QuadratureType, FieldLHS, FieldRHS>::computeAvg(
1566 : const FieldLHS& u_h) const {
1567 0 : if (this->quadrature_m.getOrder() < (2 * Order + 1)) {
1568 : // throw exception
1569 0 : throw IpplException(
1570 : "LagrangeSpace::computeAvg()",
1571 : "Order of quadrature rule for error computation should be > 2*p + 1");
1572 : }
1573 :
1574 : // List of quadrature weights
1575 0 : const Vector<T, QuadratureType::numElementNodes> w =
1576 0 : this->quadrature_m.getWeightsForRefElement();
1577 :
1578 : // List of quadrature nodes
1579 0 : const Vector<point_t, QuadratureType::numElementNodes> q =
1580 0 : this->quadrature_m.getIntegrationNodesForRefElement();
1581 :
1582 : // Evaluate the basis functions for the DOF at the quadrature nodes
1583 0 : Vector<Vector<T, numElementDOFs>, QuadratureType::numElementNodes> basis_q;
1584 0 : for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
1585 0 : for (size_t i = 0; i < numElementDOFs; ++i) {
1586 0 : basis_q[k][i] = this->evaluateRefElementShapeFunction(i, q[k]);
1587 : }
1588 : }
1589 :
1590 0 : const indices_t zeroNdIndex = Vector<size_t, Dim>(0);
1591 :
1592 : // Absolute value of det Phi_K
1593 0 : const T absDetDPhi = Kokkos::abs(this->ref_element_m.getDeterminantOfTransformationJacobian(
1594 0 : this->getElementMeshVertexPoints(zeroNdIndex)));
1595 :
1596 : // Variable to sum the error to
1597 0 : T avg = 0;
1598 :
1599 : // Get domain information and ghost cells
1600 0 : auto ldom = (u_h.getLayout()).getLocalNDIndex();
1601 0 : const int nghost = u_h.getNghost();
1602 :
1603 : using exec_space = typename Kokkos::View<const size_t*>::execution_space;
1604 : using policy_type = Kokkos::RangePolicy<exec_space>;
1605 :
1606 : // Loop over elements to compute contributions
1607 0 : Kokkos::parallel_reduce(
1608 0 : "Compute average over elements", policy_type(0, elementIndices.extent(0)),
1609 0 : KOKKOS_CLASS_LAMBDA(size_t index, double& local) {
1610 0 : const size_t elementIndex = elementIndices(index);
1611 0 : const Vector<size_t, numElementDOFs> global_dofs =
1612 0 : this->getGlobalDOFIndices(elementIndex);
1613 :
1614 : // contribution of this element to the error
1615 0 : T contrib = 0;
1616 0 : for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
1617 0 : T val_u_h = 0;
1618 0 : for (size_t i = 0; i < numElementDOFs; ++i) {
1619 : // get field index corresponding to this DOF
1620 0 : size_t I = global_dofs[i];
1621 0 : auto dof_ndindex_I = this->getMeshVertexNDIndex(I);
1622 0 : for (unsigned d = 0; d < Dim; ++d) {
1623 0 : dof_ndindex_I[d] = dof_ndindex_I[d] - ldom[d].first() + nghost;
1624 : }
1625 :
1626 : // get field value at DOF and interpolate to q_k
1627 0 : val_u_h += basis_q[k][i] * apply(u_h, dof_ndindex_I);
1628 : }
1629 :
1630 0 : contrib += w[k] * val_u_h * absDetDPhi;
1631 : }
1632 0 : local += contrib;
1633 0 : },
1634 0 : Kokkos::Sum<double>(avg));
1635 :
1636 : // MPI reduce
1637 0 : T global_avg = 0.0;
1638 0 : Comm->allreduce(avg, global_avg, 1, std::plus<T>());
1639 :
1640 0 : return global_avg;
1641 0 : }
1642 :
1643 : } // namespace ippl
|
