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