hyper.deal
Loading...
Searching...
No Matches
vector_tools.h
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2020 by the hyper.deal authors
4//
5// This file is part of the hyper.deal library.
6//
7// The hyper.deal library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 3.0 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.MD at
12// the top level directory of hyper.deal.
13//
14// ---------------------------------------------------------------------
15
16#ifndef NDIM_INTERPOLATION
17#define NDIM_INTERPOLATION
18
19#include <hyper.deal/base/config.h>
20
21#include <deal.II/base/function.h>
22
23#include <hyper.deal/matrix_free/fe_evaluation_cell.h>
24#include <hyper.deal/matrix_free/matrix_free.h>
25
26namespace hyperdeal
27{
28 namespace VectorTools
29 {
30 namespace internal
31 {
32 template <int dim_x,
33 int dim_v,
34 int degree,
35 int n_points,
36 typename Number,
37 typename VectorizedArrayType,
38 typename VectorType,
39 typename ID>
40 VectorizedArrayType *
41 interpolate(FEEvaluation<dim_x,
42 dim_v,
43 degree,
44 n_points,
45 Number,
46 VectorizedArrayType> &phi,
47 const VectorType & src,
48 const ID cell)
49 {
50 static const int dim = dim_x + dim_v;
51
52 static const dealii::internal::EvaluatorVariant tensorproduct =
53 dealii::internal::EvaluatorVariant::evaluate_evenodd;
54
55
56 // get data and scratch
57 VectorizedArrayType *data_ptr = phi.get_data_ptr();
58
59 // get cell values
60 phi.reinit(cell);
61 phi.read_dof_values(src);
62
63 dealii::internal::FEEvaluationImplBasisChange<
64 tensorproduct,
65 dealii::internal::EvaluatorQuantity::value,
66 dim,
67 degree + 1,
68 n_points>::do_forward(1, *phi.get_shape_values(), data_ptr, data_ptr);
69
70 return data_ptr;
71 }
72 } // namespace internal
73
81 template <int degree,
82 int n_points,
83 int dim_x,
84 int dim_v,
85 typename Number,
86 typename VectorType,
87 typename VectorizedArrayType>
88 void
89 interpolate(
90 const std::shared_ptr<dealii::Function<dim_x + dim_v, Number>>
91 analytical_solution,
92 const MatrixFree<dim_x, dim_v, Number, VectorizedArrayType> &matrix_free,
93 VectorType & dst,
94 const unsigned int dof_no_x,
95 const unsigned int dof_no_v,
96 const unsigned int quad_no_x,
97 const unsigned int quad_no_v)
98 {
99 const static int dim = dim_x + dim_v;
100
101 FEEvaluation<dim_x, dim_v, degree, n_points, Number, VectorizedArrayType>
102 phi(matrix_free, dof_no_x, dof_no_v, quad_no_x, quad_no_v);
103
104 const int dummy = 0;
105
106 matrix_free.template cell_loop<VectorType, int>(
107 [&](const auto &, auto &dst, const auto &, const auto cell) mutable {
108 VectorizedArrayType *data_ptr = phi.get_data_ptr();
109
110 phi.reinit(cell);
111
112 // loop over quadrature points
113 for (unsigned int qv = 0, q = 0; qv < phi.n_q_points_v; qv++)
114 for (unsigned int qx = 0; qx < phi.n_q_points_x; qx++, q++)
115 {
116 // get reference to quadrature points
117 const auto q_point = phi.get_quadrature_point(qx, qv);
118 // loop over all lanes
119 for (unsigned int v = 0; v < phi.n_vectorization_lanes_filled();
120 v++)
121 {
122 // setup point ...
123 dealii::Point<dim, Number> p;
124 for (unsigned int d = 0; d < dim; ++d)
125 p[d] = q_point[d][v];
126 // ... evaluate function at point
127 data_ptr[q][v] = analytical_solution->value(p);
128 }
129 }
130
131 phi.set_dof_values(dst);
132 },
133 dst,
134 dummy);
135 }
136
144 template <int degree,
145 int n_points,
146 int dim_x,
147 int dim_v,
148 typename Number,
149 typename VectorType,
150 typename VectorizedArrayType>
151 std::array<Number, 2>
152 norm_and_error(
153 const std::shared_ptr<dealii::Function<dim_x + dim_v, Number>>
154 analytical_solution,
155 const MatrixFree<dim_x, dim_v, Number, VectorizedArrayType> &matrix_free,
156 const VectorType & src,
157 const unsigned int dof_no_x,
158 const unsigned int dof_no_v,
159 const unsigned int quad_no_x,
160 const unsigned int quad_no_v)
161 {
162 const static int dim = dim_x + dim_v;
163
164 FEEvaluation<dim_x, dim_v, degree, n_points, Number, VectorizedArrayType>
165 phi(matrix_free, dof_no_x, dof_no_v, quad_no_x, quad_no_v);
166
167 std::array<Number, 2> result{{0.0, 0.0}};
168
169 int dummy;
170
171 matrix_free.template cell_loop<int, VectorType>(
172 [&](
173 const auto &, int &, const VectorType &src, const auto cell) mutable {
174 const VectorizedArrayType *f_ptr =
175 VectorTools::internal::interpolate(phi, src, cell);
176
177 std::array<VectorizedArrayType, 2> temp;
178 std::fill(temp.begin(), temp.end(), 0.);
179
180 for (unsigned int qv = 0, q = 0; qv < phi.n_q_points_v; qv++)
181 for (unsigned int qx = 0; qx < phi.n_q_points_x; qx++, q++)
182 {
183 // determine exact solution at quadrature point
184 VectorizedArrayType solution_q;
185 const auto q_point = phi.get_quadrature_point(qx, qv);
186
187 for (unsigned int v = 0; v < phi.n_vectorization_lanes_filled();
188 v++)
189 {
190 dealii::Point<dim> p;
191 for (unsigned int d = 0; d < dim; ++d)
192 p[d] = q_point[d][v];
193 solution_q[v] = analytical_solution->value(p);
194 }
195
196 const auto f = f_ptr[q]; // value at quadrature point
197 const auto e = f - solution_q; // error
198 const auto JxW = phi.JxW(qx, qv);
199
200 temp[0] += f * f * JxW; // L2: norm
201 temp[1] += e * e * JxW; // L2: error
202 }
203
204 // gather results
205 for (unsigned int v = 0; v < phi.n_vectorization_lanes_filled(); v++)
206 for (unsigned int i = 0; i < temp.size(); i++)
207 result[i] += temp[i][v];
208 },
209 dummy,
210 src);
211
212 MPI_Allreduce(MPI_IN_PLACE,
213 result.data(),
214 result.size(),
215 MPI_DOUBLE, // [TODO]: mpi_type_id
216 MPI_SUM,
217 matrix_free.get_communicator());
218
219 result[0] = std::sqrt(result[0]);
220 result[1] = std::sqrt(result[1]);
221
222 return result;
223 }
224
230 template <int degree,
231 int n_points,
232 int dim_x,
233 int dim_v,
234 typename Number,
235 typename VectorizedArrayType,
236 typename Vector_Out,
237 typename Vector_In>
238 void
239 velocity_space_integration(
240 const MatrixFree<dim_x, dim_v, Number, VectorizedArrayType> &data,
241 Vector_Out & dst,
242 const Vector_In & src,
243 const unsigned int dof_no_x,
244 const unsigned int dof_no_v,
245 const unsigned int quad_no_v)
246 {
247 using MF =
248 hyperdeal::MatrixFree<dim_x, dim_v, Number, VectorizedArrayType>;
249
250 using VectorizedArrayTypeX = typename MF::VectorizedArrayTypeX;
251 using VectorizedArrayTypeV = typename MF::VectorizedArrayTypeV;
252
253 FEEvaluation<dim_x, dim_v, degree, n_points, Number, VectorizedArrayType>
254 phi(const_cast<MatrixFree<dim_x, dim_v, Number, VectorizedArrayType> &>(
255 data),
256 dof_no_x,
257 dof_no_v,
258 0 /*dummy*/,
259 0 /*dummy*/);
260
261 dealii::
262 FEEvaluation<dim_x, degree, n_points, 1, Number, VectorizedArrayTypeX>
263 phi_x(data.get_matrix_free_x(), dof_no_x, 0 /*dummy*/);
264 dealii::
265 FEEvaluation<dim_v, degree, n_points, 1, Number, VectorizedArrayTypeV>
266 phi_v(data.get_matrix_free_v(), dof_no_v, quad_no_v);
267
268 dst = 0.0; // clear destination vector
269
270 data.template cell_loop<Vector_Out, Vector_In>(
271 [&](const auto &,
272 Vector_Out & dst,
273 const Vector_In &src,
274 const auto cell) mutable {
275 phi.reinit(cell);
276 phi.read_dof_values(src);
277
278 phi_x.reinit(cell.x);
279
280 unsigned int index_v = cell.v / VectorizedArrayTypeV::size();
281 unsigned int lane_v = cell.v % VectorizedArrayTypeV::size();
282 phi_v.reinit(index_v);
283
284 // get data and scratch
285 const VectorizedArrayType *data_ptr_src = phi.get_data_ptr();
286 VectorizedArrayTypeX * data_ptr_dst = phi_x.begin_dof_values();
287
288 // loop over all x points and integrate over all v points
289 for (unsigned int qx = 0; qx < phi_x.n_q_points; ++qx)
290 {
291 VectorizedArrayType sum_v = VectorizedArrayType();
292 for (unsigned int qv = 0;
293 qv < dealii::Utilities::pow<unsigned int>(n_points, dim_v);
294 ++qv)
295 sum_v +=
296 data_ptr_src[qx +
297 dealii::Utilities::pow(n_points, dim_v) * qv] *
298 phi_v.JxW(qv)[lane_v];
299 data_ptr_dst[qx] = sum_v;
300 }
301
302 phi_x.distribute_local_to_global(dst);
303 },
304 dst,
305 src);
306
307 // collect global contributions
308 MPI_Allreduce(
309 MPI_IN_PLACE,
310 dst.get_values(),
311 dst.locally_owned_size(),
312 MPI_DOUBLE,
313 MPI_SUM,
314 data.get_matrix_free_v().get_dof_handler().get_communicator());
315 }
316
322 template <int degree,
323 int n_points,
324 int dim_x,
325 int dim_v,
326 typename Number,
327 typename VectorizedArrayType,
328 typename Vector_Out,
329 typename Vector_In>
330 void
331 coordinate_space_integration(
332 const MatrixFree<dim_x, dim_v, Number, VectorizedArrayType> &data,
333 Vector_Out & dst,
334 const Vector_In & src,
335 const unsigned int dof_no_x,
336 const unsigned int dof_no_v,
337 const unsigned int quad_no_x)
338 {
339 using MF =
340 hyperdeal::MatrixFree<dim_x, dim_v, Number, VectorizedArrayType>;
341 using VectorizedArrayTypeX = typename MF::VectorizedArrayTypeX;
342 using VectorizedArrayTypeV = typename MF::VectorizedArrayTypeV;
343 static_assert(VectorizedArrayTypeV::size() == 1);
344
345 FEEvaluation<dim_x, dim_v, degree, n_points, Number, VectorizedArrayType>
346 phi(data, dof_no_x, dof_no_v, 0 /*dummy*/, 0 /*dummy*/);
347 dealii::
348 FEEvaluation<dim_x, degree, n_points, 1, Number, VectorizedArrayTypeX>
349 phi_x(data.get_matrix_free_x(), dof_no_x, quad_no_x);
350 dealii::
351 FEEvaluation<dim_v, degree, n_points, 1, Number, VectorizedArrayTypeV>
352 phi_v(data.get_matrix_free_v(), dof_no_v, 0 /*dummy*/);
353
354 // Clear destination vector
355 dst = 0.0;
356
357 constexpr auto N_v_points = dealii::Utilities::pow(n_points, dim_v);
358
359 dealii::AlignedVector<VectorizedArrayTypeX> scratch(N_v_points);
360
361 data.template cell_loop<Vector_Out, Vector_In>(
362 [&](const auto &,
363 Vector_Out & dst,
364 const Vector_In &src,
365 const auto cell) mutable {
366 phi.reinit(cell);
367 phi.read_dof_values(src);
368 phi_x.reinit(cell.x);
369 phi_v.reinit(cell.v);
370
371 const VectorizedArrayType *data_ptr_src = phi.get_data_ptr();
372 VectorizedArrayTypeV * data_ptr_dst = phi_v.begin_dof_values();
373
374 const auto N_lanes_x = phi.n_vectorization_lanes_filled();
375
376 // reduce in x-direction in a vectorized fashion
377 for (unsigned int qv = 0; qv < phi_v.n_q_points; ++qv)
378 {
379 scratch[qv] = 0.0;
380 for (unsigned int qx = 0; qx < phi_x.n_q_points; ++qx)
381 scratch[qv] +=
382 data_ptr_src[qx + qv * N_v_points] * phi_x.JxW(qx);
383 }
384
385 // perform cross-lane reduction
386 for (unsigned int qv = 0; qv < phi_v.n_q_points; ++qv)
387 {
388 data_ptr_dst[qv] = 0.0;
389 for (unsigned int lane = 0; lane < N_lanes_x; ++lane)
390 data_ptr_dst[qv] += scratch[qv][lane];
391 }
392
393 phi_v.distribute_local_to_global(dst);
394 },
395 dst,
396 src);
397
398 // Collect global contribution
399 const auto comm =
400 data.get_matrix_free_x().get_dof_handler().get_communicator();
401 MPI_Allreduce(MPI_IN_PLACE,
402 &*dst.begin(),
403 dst.locally_owned_size(),
404 MPI_DOUBLE,
405 MPI_SUM,
406 comm);
407 }
408
409 } // namespace VectorTools
410
411
412
413} // namespace hyperdeal
414
415#endif
hyperdeal::MatrixFree
Definition matrix_free.h:40
Generated by doxygen 1.12.0