Branch data 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 : 198 : 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 [ + - ]: 198 : 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 [ + - ]: 198 : initializeElementIndices(layout);
19 : 198 : }
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 : 198 : void LagrangeSpace<T, Dim, Order, ElementType, QuadratureType, FieldLHS,
52 : : FieldRHS>::initializeElementIndices(const Layout_t& layout) {
53 [ + - ]: 198 : const auto& ldom = layout.getLocalNDIndex();
54 : 198 : int npoints = ldom.size();
55 [ + - ]: 198 : auto first = ldom.first();
56 [ + - ]: 198 : auto last = ldom.last();
57 [ # # ]: 198 : ippl::Vector<double, Dim> bounds;
58 : :
59 [ + + ]: 612 : for (size_t d = 0; d < Dim; ++d) {
60 : 414 : bounds[d] = this->nr_m[d] - 1;
61 : : }
62 : :
63 : 198 : 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 [ + - ]: 198 : Kokkos::View<size_t*> points("npoints", npoints);
68 [ + - + - ]: 198 : Kokkos::parallel_reduce(
69 : 0 : "ComputePoints", npoints,
70 [ + - + - : 8118 : 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 [ + - ]: 396 : Kokkos::Sum<int>(upperBoundaryPoints));
86 [ + - + - ]: 198 : Kokkos::fence();
87 : :
88 : : // The elementIndices will be the same array as computed above,
89 : : // with the tagged upper boundary points removed.
90 : 198 : int elementsPerRank = npoints - upperBoundaryPoints;
91 [ + - + - ]: 198 : elementIndices = Kokkos::View<size_t*>("i", elementsPerRank);
92 [ + - ]: 198 : Kokkos::View<size_t> index("index");
93 : :
94 [ + - + - ]: 198 : Kokkos::parallel_for(
95 [ + - + - : 15840 : "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 : 198 : }
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 : 516 : 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 : 516 : const Vector<size_t, this->numElementDOFs> global_dofs =
130 [ + - ]: 516 : this->getGlobalDOFIndices(elementIndex);
131 : :
132 [ + - ]: 516 : ippl::Vector<size_t, this->numElementDOFs> dof_mapping;
133 : : if (Dim == 1) {
134 : 12 : dof_mapping = {0, 1};
135 : : } else if (Dim == 2) {
136 : 72 : dof_mapping = {0, 1, 3, 2};
137 : : } else if (Dim == 3) {
138 : 432 : 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 [ + + ]: 3616 : for (size_t i = 0; i < dof_mapping.dim; ++i) {
144 [ + + ]: 3268 : if (global_dofs[dof_mapping[i]] == globalDOFIndex) {
145 : 168 : return dof_mapping[i];
146 : : }
147 : : }
148 : 348 : 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 : 516 : }
152 : :
153 : : template <typename T, unsigned Dim, unsigned Order, typename ElementType,
154 : : typename QuadratureType, typename FieldLHS, typename FieldRHS>
155 : : KOKKOS_FUNCTION size_t
156 : 168 : LagrangeSpace<T, Dim, Order, ElementType, QuadratureType, FieldLHS,
157 : : FieldRHS>::getGlobalDOFIndex(const size_t& elementIndex,
158 : : const size_t& localDOFIndex) const {
159 [ + - ]: 168 : const auto global_dofs = this->getGlobalDOFIndices(elementIndex);
160 : :
161 : 336 : return global_dofs[localDOFIndex];
162 : 168 : }
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 : 54 : LagrangeSpace<T, Dim, Order, ElementType, QuadratureType, FieldLHS,
169 : : FieldRHS>::getLocalDOFIndices() const {
170 : 54 : Vector<size_t, numElementDOFs> localDOFs;
171 : :
172 [ + + ]: 178 : for (size_t dof = 0; dof < numElementDOFs; ++dof) {
173 : 124 : localDOFs[dof] = dof;
174 : : }
175 : :
176 : 54 : 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 : 738 : LagrangeSpace<T, Dim, Order, ElementType, QuadratureType, FieldLHS,
184 : : FieldRHS>::getGlobalDOFIndices(const size_t& elementIndex) const {
185 : 738 : Vector<size_t, this->numElementDOFs> globalDOFs(0);
186 : :
187 : : // get element pos
188 [ + - ]: 738 : indices_t elementPos = this->getElementNDIndex(elementIndex);
189 : :
190 : : // Compute the vector to multiply the ndindex with
191 [ + - ]: 738 : ippl::Vector<size_t, Dim> vec(1);
192 [ + + ]: 1968 : for (size_t d = 1; d < dim; ++d) {
193 [ + + ]: 3022 : for (size_t d2 = d; d2 < Dim; ++d2) {
194 : 1792 : vec[d2] *= this->nr_m[d - 1];
195 : : }
196 : : }
197 : 738 : vec *= Order; // Multiply each dimension by the order
198 : 738 : size_t smallestGlobalDOF = elementPos.dot(vec);
199 : :
200 : : // Add the vertex DOFs
201 : 738 : globalDOFs[0] = smallestGlobalDOF;
202 : 738 : globalDOFs[1] = smallestGlobalDOF + Order;
203 : :
204 : : if (Dim >= 2) {
205 : 668 : globalDOFs[2] = globalDOFs[1] + this->nr_m[1] * Order;
206 : 668 : globalDOFs[3] = globalDOFs[0] + this->nr_m[1] * Order;
207 : : }
208 : : if (Dim >= 3) {
209 : 562 : globalDOFs[4] = globalDOFs[0] + this->nr_m[1] * this->nr_m[2] * Order;
210 : 562 : globalDOFs[5] = globalDOFs[1] + this->nr_m[1] * this->nr_m[2] * Order;
211 : 562 : globalDOFs[6] = globalDOFs[2] + this->nr_m[1] * this->nr_m[2] * Order;
212 : 562 : 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 : 738 : return globalDOFs;
253 : 738 : }
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 : 131300 : 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 [ - + ]: 131300 : assert(localDOF < this->numElementDOFs
270 : : && "The local vertex index is invalid"); // TODO assumes 1st order Lagrange
271 : :
272 [ + - ]: 131300 : 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 [ + - ]: 131300 : 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 : 131300 : T product = 1;
281 : :
282 [ + + ]: 521800 : for (size_t d = 0; d < Dim; d++) {
283 [ + + ]: 390500 : if (localPoint[d] < ref_element_point[d]) {
284 : 195250 : product *= localPoint[d];
285 : : } else {
286 : 195250 : product *= 1.0 - localPoint[d];
287 : : }
288 : : }
289 : :
290 : 131300 : return product;
291 : 131300 : }
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 : 131300 : 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 [ - + ]: 131300 : assert(localDOF < this->numElementDOFs && "The local vertex index is invalid");
307 : :
308 [ + - ]: 131300 : assert(this->ref_element_m.isPointInRefElement(localPoint)
[ # # # # ]
309 : : && "Point is not in reference element");
310 : :
311 : : // Get the local dof nd_index
312 [ + - ]: 131300 : const vertex_points_t local_vertex_points = this->ref_element_m.getLocalVertices();
313 : :
314 : 131300 : const point_t& local_vertex_point = local_vertex_points[localDOF];
315 : :
316 [ + - ]: 131300 : 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 [ + + ]: 521800 : for (size_t d = 0; d < Dim; d++) {
323 : : // The variable that accumulates the product of the shape functions.
324 : 390500 : T product = 1;
325 : :
326 [ + + ]: 1555400 : for (size_t d2 = 0; d2 < Dim; d2++) {
327 [ + + ]: 1164900 : if (d2 == d) {
328 [ + + ]: 390500 : if (localPoint[d] < local_vertex_point[d]) {
329 : 195250 : product *= 1;
330 : : } else {
331 : 195250 : product *= -1;
332 : : }
333 : : } else {
334 [ + + ]: 774400 : if (localPoint[d2] < local_vertex_point[d2]) {
335 : 387200 : product *= localPoint[d2];
336 : : } else {
337 : 387200 : product *= 1.0 - localPoint[d2];
338 : : }
339 : : }
340 : : }
341 : :
342 : 390500 : gradient[d] = product;
343 : : }
344 : :
345 : 131300 : return gradient;
346 : 131300 : }
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 : 10 : FieldLHS LagrangeSpace<T, Dim, Order, ElementType, QuadratureType, FieldLHS,
356 : : FieldRHS>::evaluateAx(FieldLHS& field, F& evalFunction) const {
357 [ + - ]: 10 : Inform m("");
358 : :
359 : : // start a timer
360 [ + + + - : 10 : static IpplTimings::TimerRef evalAx = IpplTimings::getTimer("evaluateAx");
+ - - - ]
361 [ + - ]: 10 : IpplTimings::startTimer(evalAx);
362 : :
363 : : // get number of ghost cells in field
364 : 10 : 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 [ + - + - ]: 10 : FieldLHS resultField(field.get_mesh(), field.getLayout(), nghost);
369 : :
370 : : // List of quadrature weights
371 : 10 : const Vector<T, QuadratureType::numElementNodes> w =
372 [ + - ]: 10 : this->quadrature_m.getWeightsForRefElement();
373 : :
374 : : // List of quadrature nodes
375 : 10 : const Vector<point_t, QuadratureType::numElementNodes> q =
376 [ + - ]: 10 : 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 [ + - ]: 10 : Vector<Vector<point_t, this->numElementDOFs>, QuadratureType::numElementNodes> grad_b_q;
381 [ + + ]: 20 : for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
382 [ + + ]: 30 : for (size_t i = 0; i < this->numElementDOFs; ++i) {
383 [ + - ]: 20 : 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 [ + - ]: 10 : ViewType view = field.getView();
390 [ + - ]: 10 : AtomicViewType resultView = resultField.getView();
391 : :
392 : : // Get boundary conditions from field
393 : 10 : BConds<FieldLHS, Dim>& bcField = field.getFieldBC();
394 [ + - ]: 10 : FieldBC bcType = bcField[0]->getBCType();
395 : :
396 : : // Get domain information
397 [ + - + - ]: 10 : 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 [ + + + - : 10 : static IpplTimings::TimerRef outer_loop = IpplTimings::getTimer("evaluateAx: outer loop");
+ - - - ]
404 [ + - ]: 10 : IpplTimings::startTimer(outer_loop);
405 : :
406 : : // Loop over elements to compute contributions
407 [ + - + - ]: 10 : Kokkos::parallel_for(
408 [ + - ]: 20 : "Loop over elements", policy_type(0, elementIndices.extent(0)),
409 [ + - + - : 60 : KOKKOS_CLASS_LAMBDA(const size_t index) {
+ - - - -
- - - - -
- - ]
410 : 0 : const size_t elementIndex = elementIndices(index);
411 [ # # ][ # # : 40 : const Vector<size_t, this->numElementDOFs> local_dof = this->getLocalDOFIndices();
# # # # ]
[ - - + -
- - + - -
- - - ]
412 : 0 : const Vector<size_t, this->numElementDOFs> global_dofs =
413 [ # # ][ # # : 40 : this->getGlobalDOFIndices(elementIndex);
# # # # ]
[ - - + -
- - + - -
- - - ]
414 [ # # ][ # # : 40 : Vector<indices_t, this->numElementDOFs> global_dof_ndindices;
# # # # ]
[ - - + -
- - + - -
- - - ]
415 : :
416 [ # # ][ # # : 120 : for (size_t i = 0; i < this->numElementDOFs; ++i) {
# # # # ]
[ - - + +
- - + + -
- - - ]
417 [ # # ][ # # : 80 : 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 [ # # ][ # # : 40 : Vector<Vector<T, this->numElementDOFs>, this->numElementDOFs> A_K;
# # # # ]
[ - - + -
- - + - -
- - - ]
425 : :
426 : : // 1. Compute the Galerkin element matrix A_K
427 [ # # ][ # # : 120 : for (i = 0; i < this->numElementDOFs; ++i) {
# # # # ]
[ - - + +
- - + + -
- - - ]
428 [ # # ][ # # : 240 : for (j = 0; j < this->numElementDOFs; ++j) {
# # # # ]
[ - - + +
- - + + -
- - - ]
429 : 0 : A_K[i][j] = 0.0;
430 [ # # ][ # # : 320 : for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
# # # # ]
[ - - + +
- - + + -
- - - ]
431 [ # # ][ # # : 160 : 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 [ # # # # ]: 40 : indices_t I_nd, J_nd;
[ - - - -
+ - + - -
- - - + -
+ - - - -
- - - -
- ]
439 : :
440 : : // 2. Compute the contribution to resultAx = A*x with A_K
441 [ # # ][ # # : 120 : for (i = 0; i < this->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 [ # # # # : 80 : 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 [ # # # # : 80 : } 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 [ # # ][ # # : 120 : for (unsigned d = 0; d < Dim; ++d) {
# # # # ]
[ - - + +
- - + + -
- - - ]
459 : 0 : I_nd[d] = I_nd[d] - ldom[d].first() + nghost;
460 : : }
461 : :
462 [ # # ][ # # : 180 : for (j = 0; j < this->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 [ # # # # : 120 : && this->isDOFOnBoundary(J_nd)) {
# # ][ # #
# # # # #
# ][ # # #
# # # # #
# # # # #
# # # # #
# # ][ - -
- - - - -
+ + + + +
- - - - -
- - + + +
+ + - - -
- - - - -
- - - - ]
468 : 0 : continue;
469 : : }
470 : :
471 : : // get the appropriate index for the Kokkos view of the field
472 [ # # ][ # # : 200 : for (unsigned d = 0; d < Dim; ++d) {
# # # # ]
[ - - + +
- - + + -
- - - ]
473 : 0 : J_nd[d] = J_nd[d] - ldom[d].first() + nghost;
474 : : }
475 : :
476 [ # # # # : 100 : apply(resultView, I_nd) += A_K[i][j] * apply(view, J_nd);
# # ][ # #
# # # # #
# # # # #
# # # # #
# ][ - - -
- - - + -
+ - + - -
- - - - -
+ - + - +
- - - - -
- - - - -
- - - ]
477 : : }
478 : : }
479 : 0 : });
480 [ + - ]: 10 : IpplTimings::stopTimer(outer_loop);
481 : :
482 [ - + ]: 10 : if (bcType == PERIODIC_FACE) {
483 [ # # ]: 0 : resultField.accumulateHalo();
484 [ # # ]: 0 : bcField.apply(resultField);
485 [ # # ]: 0 : bcField.assignGhostToPhysical(resultField);
486 : : } else {
487 [ + - ]: 10 : resultField.accumulateHalo_noghost();
488 : : }
489 : :
490 [ + - ]: 10 : IpplTimings::stopTimer(evalAx);
491 : :
492 : 10 : return resultField;
493 : 10 : }
494 : :
495 : : template <typename T, unsigned Dim, unsigned Order, typename ElementType,
496 : : typename QuadratureType, typename FieldLHS, typename FieldRHS>
497 : : template <typename F>
498 : 0 : FieldLHS LagrangeSpace<T, Dim, Order, ElementType, QuadratureType, FieldLHS,
499 : : FieldRHS>::evaluateAx_lift(FieldLHS& field, F& evalFunction) const {
500 [ # # ]: 0 : Inform m("");
501 : :
502 : : // get number of ghost cells in field
503 : 0 : const int nghost = field.getNghost();
504 : :
505 : : // create a new field for result with view initialized to zero (views are initialized to
506 : : // zero by default)
507 [ # # # # ]: 0 : FieldLHS resultField(field.get_mesh(), field.getLayout(), nghost);
508 : :
509 : : // List of quadrature weights
510 : 0 : const Vector<T, QuadratureType::numElementNodes> w =
511 [ # # ]: 0 : this->quadrature_m.getWeightsForRefElement();
512 : :
513 : : // List of quadrature nodes
514 : 0 : const Vector<point_t, QuadratureType::numElementNodes> q =
515 [ # # ]: 0 : this->quadrature_m.getIntegrationNodesForRefElement();
516 : :
517 : : // TODO move outside of evaluateAx (I think it is possible for other problems as well)
518 : : // Gradients of the basis functions for the DOF at the quadrature nodes
519 [ # # ]: 0 : Vector<Vector<point_t, this->numElementDOFs>, QuadratureType::numElementNodes> grad_b_q;
520 [ # # ]: 0 : for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
521 [ # # ]: 0 : for (size_t i = 0; i < this->numElementDOFs; ++i) {
522 [ # # ]: 0 : grad_b_q[k][i] = this->evaluateRefElementShapeFunctionGradient(i, q[k]);
523 : : }
524 : : }
525 : :
526 : : // Get field data and atomic result data,
527 : : // since it will be added to during the kokkos loop
528 [ # # ]: 0 : ViewType view = field.getView();
529 [ # # ]: 0 : AtomicViewType resultView = resultField.getView();
530 : :
531 : : // Get domain information
532 [ # # # # ]: 0 : auto ldom = (field.getLayout()).getLocalNDIndex();
533 : :
534 : : using exec_space = typename Kokkos::View<const size_t*>::execution_space;
535 : : using policy_type = Kokkos::RangePolicy<exec_space>;
536 : :
537 : : // Loop over elements to compute contributions
538 [ # # # # ]: 0 : Kokkos::parallel_for(
539 [ # # ]: 0 : "Loop over elements", policy_type(0, elementIndices.extent(0)),
540 [ # # # # : 0 : KOKKOS_CLASS_LAMBDA(const size_t index) {
# # # # #
# # # # #
# # ]
541 : 0 : const size_t elementIndex = elementIndices(index);
542 [ # # ][ # # : 0 : const Vector<size_t, this->numElementDOFs> local_dof = this->getLocalDOFIndices();
# # # # ]
543 : 0 : const Vector<size_t, this->numElementDOFs> global_dofs =
544 [ # # ][ # # : 0 : this->getGlobalDOFIndices(elementIndex);
# # # # ]
545 [ # # ][ # # : 0 : Vector<indices_t, this->numElementDOFs> global_dof_ndindices;
# # # # ]
546 : :
547 [ # # ][ # # : 0 : for (size_t i = 0; i < this->numElementDOFs; ++i) {
# # # # ]
548 [ # # ][ # # : 0 : global_dof_ndindices[i] = this->getMeshVertexNDIndex(global_dofs[i]);
# # # # ]
549 : : }
550 : :
551 : : // local DOF indices
552 : : size_t i, j;
553 : :
554 : : // Element matrix
555 [ # # ][ # # : 0 : Vector<Vector<T, this->numElementDOFs>, this->numElementDOFs> A_K;
# # # # ]
556 : :
557 : : // 1. Compute the Galerkin element matrix A_K
558 [ # # ][ # # : 0 : for (i = 0; i < this->numElementDOFs; ++i) {
# # # # ]
559 [ # # ][ # # : 0 : for (j = 0; j < this->numElementDOFs; ++j) {
# # # # ]
560 : 0 : A_K[i][j] = 0.0;
561 [ # # ][ # # : 0 : for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
# # # # ]
562 [ # # ][ # # : 0 : A_K[i][j] += w[k] * evalFunction(i, j, grad_b_q[k]);
# # # # ]
563 : : }
564 : : }
565 : : }
566 : :
567 : : // global DOF n-dimensional indices (Vector of N indices representing indices in
568 : : // each dimension)
569 [ # # # # ]: 0 : indices_t I_nd, J_nd;
570 : :
571 : : // 2. Compute the contribution to resultAx = A*x with A_K
572 [ # # ][ # # : 0 : for (i = 0; i < this->numElementDOFs; ++i) {
# # # # ]
573 : 0 : I_nd = global_dof_ndindices[i];
574 : :
575 : : // Skip if on a row of the matrix
576 [ # # ]: 0 : if (this->isDOFOnBoundary(I_nd)) {
[ # # # # ]
[ # # # #
# # # # ]
577 : 0 : continue;
578 : : }
579 : :
580 : : // get the appropriate index for the Kokkos view of the field
581 [ # # ][ # # : 0 : for (unsigned d = 0; d < Dim; ++d) {
# # # # ]
582 : 0 : I_nd[d] = I_nd[d] - ldom[d].first() + nghost;
583 : : }
584 : :
585 [ # # ][ # # : 0 : for (j = 0; j < this->numElementDOFs; ++j) {
# # # # ]
586 : 0 : J_nd = global_dof_ndindices[j];
587 : :
588 : : // Contribute to lifting only if on a boundary DOF
589 [ # # ]: 0 : if (this->isDOFOnBoundary(J_nd)) {
[ # # # # ]
[ # # # #
# # # # ]
590 : : // get the appropriate index for the Kokkos view of the field
591 [ # # ][ # # : 0 : for (unsigned d = 0; d < Dim; ++d) {
# # # # ]
592 : 0 : J_nd[d] = J_nd[d] - ldom[d].first() + nghost;
593 : : }
594 [ # # # # : 0 : apply(resultView, I_nd) += A_K[i][j] * apply(view, J_nd);
# # ][ # #
# # # # #
# # # # #
# # # # #
# ]
595 : 0 : continue;
596 : 0 : }
597 : :
598 : : }
599 : : }
600 : 0 : });
601 [ # # ]: 0 : resultField.accumulateHalo();
602 : :
603 : 0 : return resultField;
604 : 0 : }
605 : :
606 : : template <typename T, unsigned Dim, unsigned Order, typename ElementType,
607 : : typename QuadratureType, typename FieldLHS, typename FieldRHS>
608 : 2 : void LagrangeSpace<T, Dim, Order, ElementType, QuadratureType, FieldLHS,
609 : : FieldRHS>::evaluateLoadVector(FieldRHS& field) const {
610 [ + - ]: 2 : Inform m("");
611 : :
612 : : // start a timer
613 [ + - + - : 2 : static IpplTimings::TimerRef evalLoadV = IpplTimings::getTimer("evaluateLoadVector");
+ - - - ]
614 [ + - ]: 2 : IpplTimings::startTimer(evalLoadV);
615 : :
616 : : // List of quadrature weights
617 : 2 : const Vector<T, QuadratureType::numElementNodes> w =
618 [ + - ]: 2 : this->quadrature_m.getWeightsForRefElement();
619 : :
620 : : // List of quadrature nodes
621 : 2 : const Vector<point_t, QuadratureType::numElementNodes> q =
622 [ + - ]: 2 : this->quadrature_m.getIntegrationNodesForRefElement();
623 : :
624 [ + - ]: 2 : const indices_t zeroNdIndex = Vector<size_t, Dim>(0);
625 : :
626 : : // Evaluate the basis functions for the DOF at the quadrature nodes
627 [ + - ]: 2 : Vector<Vector<T, this->numElementDOFs>, QuadratureType::numElementNodes> basis_q;
628 [ + + ]: 12 : for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
629 [ + + ]: 30 : for (size_t i = 0; i < this->numElementDOFs; ++i) {
630 [ + - ]: 20 : basis_q[k][i] = this->evaluateRefElementShapeFunction(i, q[k]);
631 : : }
632 : : }
633 : :
634 : : // Absolute value of det Phi_K
635 [ + - ]: 2 : const T absDetDPhi = Kokkos::abs(this->ref_element_m.getDeterminantOfTransformationJacobian(
636 [ + - ]: 4 : this->getElementMeshVertexPoints(zeroNdIndex)));
637 : :
638 : : // Get domain information and ghost cells
639 [ + - + - ]: 2 : auto ldom = (field.getLayout()).getLocalNDIndex();
640 : 2 : const int nghost = field.getNghost();
641 : :
642 : : // Get boundary conditions from field
643 : 2 : BConds<FieldRHS, Dim>& bcField = field.getFieldBC();
644 [ + - ]: 2 : FieldBC bcType = bcField[0]->getBCType();
645 : :
646 [ + - + - ]: 2 : FieldRHS temp_field(field.get_mesh(), field.getLayout(), nghost);
647 [ + - ]: 2 : temp_field.setFieldBC(bcField);
648 : :
649 : : // Get field data and make it atomic,
650 : : // since it will be added to during the kokkos loop
651 : : // We work with a temporary field since we need to use field
652 : : // to evaluate the load vector; then we assign temp to RHS field
653 [ + - ]: 2 : AtomicViewType atomic_view = temp_field.getView();
654 : :
655 : : using exec_space = typename Kokkos::View<const size_t*>::execution_space;
656 : : using policy_type = Kokkos::RangePolicy<exec_space>;
657 : :
658 : : // start a timer
659 [ + - + - : 2 : static IpplTimings::TimerRef outer_loop =
- - ]
660 [ + - ]: 2 : IpplTimings::getTimer("evaluateLoadVec: outer loop");
661 [ + - ]: 2 : IpplTimings::startTimer(outer_loop);
662 : :
663 : : // Loop over elements to compute contributions
664 [ + - + - ]: 2 : Kokkos::parallel_for(
665 [ + - ]: 4 : "Loop over elements", policy_type(0, elementIndices.extent(0)),
666 [ + - + - : 12 : KOKKOS_CLASS_LAMBDA(size_t index) {
+ - - - -
- - - -
- ]
667 : 0 : const size_t elementIndex = elementIndices(index);
668 [ # # ][ # # : 8 : const Vector<size_t, this->numElementDOFs> local_dofs = this->getLocalDOFIndices();
# # # # ]
[ - - - -
- - - - -
- + - - -
+ - - - -
- - - -
- ][ # # #
# # # #
# ]
669 : 0 : const Vector<size_t, this->numElementDOFs> global_dofs =
670 [ # # ][ # # : 8 : this->getGlobalDOFIndices(elementIndex);
# # # # ]
[ - - - -
- - - - -
- + - - -
+ - - - -
- - - -
- ][ # # #
# # # #
# ]
671 : :
672 : : size_t i, I;
673 : :
674 : : // 1. Compute b_K
675 [ # # # # ]: 40 : for (i = 0; i < this->numElementDOFs; ++i) {
[ # # # #
# # # # #
# # # ][ -
- - - - -
- - - - -
- - - - -
- - - - +
+ + + - -
- - + + +
+ - - - -
- - - - -
- - - - -
- - ][ # #
# # # # #
# # # # #
# # # # ]
676 : 0 : I = global_dofs[i];
677 : :
678 : : // TODO fix for higher order
679 [ # # ][ # # : 16 : auto dof_ndindex_I = this->getMeshVertexNDIndex(I);
# # # # ]
[ - - - -
- - - - -
- + - - -
+ - - - -
- - - -
- ][ # # #
# # # #
# ]
680 : :
681 : : // Skip boundary DOFs (Zero and Constant Dirichlet BCs)
682 [ # # ][ # # : 0 : if (((bcType == ZERO_FACE) || (bcType == CONSTANT_FACE))
# # # # ]
[ # # # #
# # # # #
# # # # #
# # # # #
# # # #
# ][ # # #
# # # #
# ]
683 [ # # # # : 16 : && (this->isDOFOnBoundary(dof_ndindex_I))) {
# # ][ # #
# # # # #
# ][ # # #
# # # # #
# # # # #
# # # # #
# # ][ - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - +
+ + + + -
- - - - -
- + + + +
+ - - - -
- - - - -
- - - - -
- - - - -
- - - -
- ][ # # #
# # # # #
# # # # #
# # # # #
# # # # #
# ]
684 : 0 : continue;
685 : : }
686 : :
687 : : // calculate the contribution of this element
688 : 0 : T contrib = 0;
689 [ # # ][ # # : 72 : for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
# # # # ]
[ - - - -
- - - - -
- + + - -
+ + - - -
- - - -
- ][ # # #
# # # #
# ]
690 : 0 : T val = 0;
691 [ # # ][ # # : 180 : for (size_t j = 0; j < this->numElementDOFs; ++j) {
# # # # ]
[ - - - -
- - - - -
- + + - -
+ + - - -
- - - -
- ][ # # #
# # # #
# ]
692 : : // get field index corresponding to this DOF
693 : 0 : size_t J = global_dofs[j];
694 [ # # ][ # # : 120 : auto dof_ndindex_J = this->getMeshVertexNDIndex(J);
# # # # ]
[ - - - -
- - - - -
- + - - -
+ - - - -
- - - -
- ][ # # #
# # # #
# ]
695 [ # # ][ # # : 240 : for (unsigned d = 0; d < Dim; ++d) {
# # # # ]
[ - - - -
- - - - -
- + + - -
+ + - - -
- - - -
- ][ # # #
# # # #
# ]
696 : 0 : dof_ndindex_J[d] = dof_ndindex_J[d] - ldom[d].first() + nghost;
697 : : }
698 : :
699 : : // get field value at DOF and interpolate to q_k
700 [ # # ][ # # : 120 : val += basis_q[k][j] * apply(field, dof_ndindex_J);
# # # # ]
[ - - - -
- - - - -
- + - - -
+ - - - -
- - - -
- ][ # # #
# # # #
# ]
701 : : }
702 : :
703 : 0 : contrib += w[k] * basis_q[k][i] * absDetDPhi * val;
704 : : }
705 : :
706 : : // get the appropriate index for the Kokkos view of the field
707 [ # # ][ # # : 24 : for (unsigned d = 0; d < Dim; ++d) {
# # # # ]
[ - - - -
- - - - -
- + + - -
+ + - - -
- - - -
- ][ # # #
# # # #
# ]
708 : 0 : dof_ndindex_I[d] = dof_ndindex_I[d] - ldom[d].first() + nghost;
709 : : }
710 : :
711 : : // add the contribution of the element to the field
712 [ # # # # ]: 12 : apply(atomic_view, dof_ndindex_I) += contrib;
[ # # # #
# # # # #
# # # ][ -
- - - - -
- - - - -
- - - - -
- - - - +
- + - - -
- - + - +
- - - - -
- - - - -
- - - - -
- - ][ # #
# # # # #
# # # # #
# # # # ]
713 : :
714 : : }
715 : 0 : });
716 [ + - ]: 2 : IpplTimings::stopTimer(outer_loop);
717 : :
718 [ + - ]: 2 : temp_field.accumulateHalo();
719 : :
720 [ + - - + ]: 2 : if ((bcType == PERIODIC_FACE) || (bcType == CONSTANT_FACE)) {
721 [ # # ]: 0 : bcField.apply(temp_field);
722 [ # # ]: 0 : bcField.assignGhostToPhysical(temp_field);
723 : : }
724 : :
725 [ + - ]: 2 : field = temp_field;
726 : :
727 [ + - ]: 2 : IpplTimings::stopTimer(evalLoadV);
728 : 2 : }
729 : :
730 : : template <typename T, unsigned Dim, unsigned Order, typename ElementType,
731 : : typename QuadratureType, typename FieldLHS, typename FieldRHS>
732 : : template <typename F>
733 : 0 : T LagrangeSpace<T, Dim, Order, ElementType, QuadratureType, FieldLHS, FieldRHS>::computeErrorL2(
734 : : const FieldLHS& u_h, const F& u_sol) const {
735 [ # # ]: 0 : if (this->quadrature_m.getOrder() < (2 * Order + 1)) {
736 : : // throw exception
737 [ # # # # : 0 : throw IpplException(
# # ]
738 : : "LagrangeSpace::computeErrorL2()",
739 : : "Order of quadrature rule for error computation should be > 2*p + 1");
740 : : }
741 : :
742 : : // List of quadrature weights
743 : 0 : const Vector<T, QuadratureType::numElementNodes> w =
744 [ # # ]: 0 : this->quadrature_m.getWeightsForRefElement();
745 : :
746 : : // List of quadrature nodes
747 : 0 : const Vector<point_t, QuadratureType::numElementNodes> q =
748 [ # # ]: 0 : this->quadrature_m.getIntegrationNodesForRefElement();
749 : :
750 : : // Evaluate the basis functions for the DOF at the quadrature nodes
751 [ # # ]: 0 : Vector<Vector<T, this->numElementDOFs>, QuadratureType::numElementNodes> basis_q;
752 [ # # ]: 0 : for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
753 [ # # ]: 0 : for (size_t i = 0; i < this->numElementDOFs; ++i) {
754 [ # # ]: 0 : basis_q[k][i] = this->evaluateRefElementShapeFunction(i, q[k]);
755 : : }
756 : : }
757 : :
758 [ # # ]: 0 : const indices_t zeroNdIndex = Vector<size_t, Dim>(0);
759 : :
760 : : // Absolute value of det Phi_K
761 [ # # ]: 0 : const T absDetDPhi = Kokkos::abs(this->ref_element_m.getDeterminantOfTransformationJacobian(
762 [ # # ]: 0 : this->getElementMeshVertexPoints(zeroNdIndex)));
763 : :
764 : : // Variable to sum the error to
765 : 0 : T error = 0;
766 : :
767 : : // Get domain information and ghost cells
768 [ # # # # ]: 0 : auto ldom = (u_h.getLayout()).getLocalNDIndex();
769 : 0 : const int nghost = u_h.getNghost();
770 : :
771 : : using exec_space = typename Kokkos::View<const size_t*>::execution_space;
772 : : using policy_type = Kokkos::RangePolicy<exec_space>;
773 : :
774 : : // Loop over elements to compute contributions
775 [ # # # # ]: 0 : Kokkos::parallel_reduce(
776 [ # # ]: 0 : "Compute error over elements", policy_type(0, elementIndices.extent(0)),
777 [ # # # # : 0 : KOKKOS_CLASS_LAMBDA(size_t index, double& local) {
# # # # #
# # # #
# ]
778 : 0 : const size_t elementIndex = elementIndices(index);
779 : 0 : const Vector<size_t, this->numElementDOFs> global_dofs =
780 [ # # ][ # # : 0 : this->getGlobalDOFIndices(elementIndex);
# # # # ]
781 : :
782 : : // contribution of this element to the error
783 : 0 : T contrib = 0;
784 [ # # ][ # # : 0 : for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
# # # # ]
785 [ # # ]: 0 : T val_u_sol = u_sol(this->ref_element_m.localToGlobal(
[ # # # # ]
[ # # # #
# # # # #
# # # ][ #
# # # #
# ]
786 [ # # # # ]: 0 : this->getElementMeshVertexPoints(this->getElementNDIndex(elementIndex)),
[ # # # #
# # # # #
# # # ]
787 : 0 : q[k]));
788 : :
789 : 0 : T val_u_h = 0;
790 [ # # ][ # # : 0 : for (size_t i = 0; i < this->numElementDOFs; ++i) {
# # # # ]
791 : : // get field index corresponding to this DOF
792 : 0 : size_t I = global_dofs[i];
793 [ # # ][ # # : 0 : auto dof_ndindex_I = this->getMeshVertexNDIndex(I);
# # # # ]
794 [ # # ][ # # : 0 : for (unsigned d = 0; d < Dim; ++d) {
# # # # ]
795 : 0 : dof_ndindex_I[d] = dof_ndindex_I[d] - ldom[d].first() + nghost;
796 : : }
797 : :
798 : : // get field value at DOF and interpolate to q_k
799 [ # # ][ # # : 0 : val_u_h += basis_q[k][i] * apply(u_h, dof_ndindex_I);
# # # # ]
800 : : }
801 : :
802 : 0 : contrib += w[k] * Kokkos::pow(val_u_sol - val_u_h, 2) * absDetDPhi;
803 : : }
804 : 0 : local += contrib;
805 : 0 : },
806 [ # # ]: 0 : Kokkos::Sum<double>(error));
807 : :
808 : : // MPI reduce
809 : 0 : T global_error = 0.0;
810 [ # # ]: 0 : Comm->allreduce(error, global_error, 1, std::plus<T>());
811 : :
812 : 0 : return Kokkos::sqrt(global_error);
813 : 0 : }
814 : :
815 : : template <typename T, unsigned Dim, unsigned Order, typename ElementType,
816 : : typename QuadratureType, typename FieldLHS, typename FieldRHS>
817 : 0 : T LagrangeSpace<T, Dim, Order, ElementType, QuadratureType, FieldLHS, FieldRHS>::computeAvg(
818 : : const FieldLHS& u_h) const {
819 [ # # ]: 0 : if (this->quadrature_m.getOrder() < (2 * Order + 1)) {
820 : : // throw exception
821 [ # # # # : 0 : throw IpplException(
# # ]
822 : : "LagrangeSpace::computeAvg()",
823 : : "Order of quadrature rule for error computation should be > 2*p + 1");
824 : : }
825 : :
826 : : // List of quadrature weights
827 : 0 : const Vector<T, QuadratureType::numElementNodes> w =
828 [ # # ]: 0 : this->quadrature_m.getWeightsForRefElement();
829 : :
830 : : // List of quadrature nodes
831 : 0 : const Vector<point_t, QuadratureType::numElementNodes> q =
832 [ # # ]: 0 : this->quadrature_m.getIntegrationNodesForRefElement();
833 : :
834 : : // Evaluate the basis functions for the DOF at the quadrature nodes
835 [ # # ]: 0 : Vector<Vector<T, this->numElementDOFs>, QuadratureType::numElementNodes> basis_q;
836 [ # # ]: 0 : for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
837 [ # # ]: 0 : for (size_t i = 0; i < this->numElementDOFs; ++i) {
838 [ # # ]: 0 : basis_q[k][i] = this->evaluateRefElementShapeFunction(i, q[k]);
839 : : }
840 : : }
841 : :
842 [ # # ]: 0 : const indices_t zeroNdIndex = Vector<size_t, Dim>(0);
843 : :
844 : : // Absolute value of det Phi_K
845 [ # # ]: 0 : const T absDetDPhi = Kokkos::abs(this->ref_element_m.getDeterminantOfTransformationJacobian(
846 [ # # ]: 0 : this->getElementMeshVertexPoints(zeroNdIndex)));
847 : :
848 : : // Variable to sum the error to
849 : 0 : T avg = 0;
850 : :
851 : : // Get domain information and ghost cells
852 [ # # # # ]: 0 : auto ldom = (u_h.getLayout()).getLocalNDIndex();
853 : 0 : const int nghost = u_h.getNghost();
854 : :
855 : : using exec_space = typename Kokkos::View<const size_t*>::execution_space;
856 : : using policy_type = Kokkos::RangePolicy<exec_space>;
857 : :
858 : : // Loop over elements to compute contributions
859 [ # # # # ]: 0 : Kokkos::parallel_reduce(
860 [ # # ]: 0 : "Compute average over elements", policy_type(0, elementIndices.extent(0)),
861 [ # # # # : 0 : KOKKOS_CLASS_LAMBDA(size_t index, double& local) {
# # # # #
# # # ]
862 : 0 : const size_t elementIndex = elementIndices(index);
863 : 0 : const Vector<size_t, this->numElementDOFs> global_dofs =
864 [ # # # # : 0 : this->getGlobalDOFIndices(elementIndex);
# # ]
865 : :
866 : : // contribution of this element to the error
867 : 0 : T contrib = 0;
868 [ # # # # : 0 : for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
# # ]
869 : 0 : T val_u_h = 0;
870 [ # # # # : 0 : for (size_t i = 0; i < this->numElementDOFs; ++i) {
# # ]
871 : : // get field index corresponding to this DOF
872 : 0 : size_t I = global_dofs[i];
873 [ # # # # : 0 : auto dof_ndindex_I = this->getMeshVertexNDIndex(I);
# # ]
874 [ # # # # : 0 : for (unsigned d = 0; d < Dim; ++d) {
# # ]
875 : 0 : dof_ndindex_I[d] = dof_ndindex_I[d] - ldom[d].first() + nghost;
876 : : }
877 : :
878 : : // get field value at DOF and interpolate to q_k
879 [ # # # # : 0 : val_u_h += basis_q[k][i] * apply(u_h, dof_ndindex_I);
# # ]
880 : : }
881 : :
882 : 0 : contrib += w[k] * val_u_h * absDetDPhi;
883 : : }
884 : 0 : local += contrib;
885 : 0 : },
886 [ # # ]: 0 : Kokkos::Sum<double>(avg));
887 : :
888 : : // MPI reduce
889 : 0 : T global_avg = 0.0;
890 [ # # ]: 0 : Comm->allreduce(avg, global_avg, 1, std::plus<T>());
891 : :
892 : 0 : return global_avg;
893 : 0 : }
894 : :
895 : : } // namespace ippl
|
