hyper.deal
Loading...
Searching...
No Matches
advection_operation.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 HYPERDEAL_NDIM_OPERATORS_ADVECTIONOPERATION
17#define HYPERDEAL_NDIM_OPERATORS_ADVECTIONOPERATION
18
19#include <hyper.deal/base/config.h>
20
21#include <deal.II/matrix_free/evaluation_kernels.h>
22
23#include <hyper.deal/base/dynamic_convergence_table.h>
24#include <hyper.deal/matrix_free/evaluation_kernels.h>
25#include <hyper.deal/matrix_free/fe_evaluation_cell.h>
26#include <hyper.deal/matrix_free/fe_evaluation_cell_inverse.h>
27#include <hyper.deal/matrix_free/fe_evaluation_face.h>
28#include <hyper.deal/matrix_free/matrix_free.h>
29#include <hyper.deal/matrix_free/tools.h>
30#include <hyper.deal/operators/advection/advection_operation_parameters.h>
31#include <hyper.deal/operators/advection/boundary_descriptor.h>
32
33namespace hyperdeal
34{
35 namespace advection
36 {
37 enum class AdvectionOperationEvaluationLevel
38 {
39 cell,
40 all_without_neighbor_load,
41 all
42 };
43
48 template <int dim_x,
49 int dim_v,
50 int degree,
51 int n_points,
52 typename Number,
53 typename VectorType,
54 typename VelocityField,
55 typename VectorizedArrayType>
56 class AdvectionOperation
57 {
58 public:
59 using This = AdvectionOperation<dim_x,
60 dim_v,
61 degree,
62 n_points,
63 Number,
64 VectorType,
65 VelocityField,
66 VectorizedArrayType>;
67 using VNumber = VectorizedArrayType;
68
69 static const int dim = dim_x + dim_v;
70
71 static const dealii::internal::EvaluatorVariant tensorproduct =
72 dealii::internal::EvaluatorVariant::evaluate_evenodd;
73
74 using FECellEval =
75 FEEvaluation<dim_x, dim_v, degree, n_points, Number, VNumber>;
76 using FEFaceEval =
77 FEFaceEvaluation<dim_x, dim_v, degree, n_points, Number, VNumber>;
78 using FECellEval_inv =
79 FEEvaluationInverse<dim_x, dim_v, degree, n_points, Number, VNumber>;
80
81
85 AdvectionOperation(
86 const MatrixFree<dim_x, dim_v, Number, VectorizedArrayType> &data,
87 DynamicConvergenceTable & table)
88 : data(data)
89 , table(table)
90 , do_collocation(false)
91 {}
92
97 void
98 reinit(
99 std::shared_ptr<BoundaryDescriptor<dim, Number>> boundary_descriptor,
100 std::shared_ptr<VelocityField> velocity_field,
101 const AdvectionOperationParamters additional_data)
102 {
103 this->factor_skew = additional_data.factor_skew;
104 this->boundary_descriptor = boundary_descriptor;
105 this->velocity_field = velocity_field;
106
107 AssertDimension(
108 (data.get_matrix_free_x().get_shape_info(0, 0).data[0].element_type ==
109 dealii::internal::MatrixFreeFunctions::ElementType::
110 tensor_symmetric_collocation),
111 (data.get_matrix_free_v().get_shape_info(0, 0).data[0].element_type ==
112 dealii::internal::MatrixFreeFunctions::ElementType::
113 tensor_symmetric_collocation));
114
115 const bool do_collocation =
116 data.get_matrix_free_x().get_shape_info(0, 0).data[0].element_type ==
117 dealii::internal::MatrixFreeFunctions::ElementType::
118 tensor_symmetric_collocation;
119
120 this->do_collocation = do_collocation;
121
122 // clang-format off
123 phi_cell.reset(new FEEvaluation<dim_x, dim_v, degree, n_points, Number, VNumber>(data, 0, 0, 0, 0));
124 phi_cell_inv.reset(new FEEvaluationInverse<dim_x, dim_v, degree, n_points, Number, VNumber>(data, 0, 0, 0, 0));
125 phi_cell_inv_co.reset(new FEEvaluationInverse<dim_x, dim_v, degree, degree + 1, Number, VNumber>(data, 0, 0, do_collocation ? 0 : 1, do_collocation ? 0 : 1));
126 phi_face_m.reset(new FEFaceEvaluation<dim_x, dim_v, degree, n_points, Number, VNumber>(data, true, 0, 0, 0, 0));
127 phi_face_p.reset(new FEFaceEvaluation<dim_x, dim_v, degree, n_points, Number, VNumber>(data, false, 0, 0, 0, 0));
128 // clang-format on
129 }
130
134 template <AdvectionOperationEvaluationLevel eval_level =
135 AdvectionOperationEvaluationLevel::all>
136 void
137 apply(VectorType & dst,
138 const VectorType &src,
139 const Number time,
140 Timers * timers = nullptr)
141 {
142 // set time of boundary functions
143 boundary_descriptor->set_time(time);
144
145 // TODO: also for velocity field
146
147 // loop over all cells/faces in phase-space
148 if (!data.is_ecl_supported()) // FCL
149 {
150 if (timers != nullptr)
151 timers->enter("FCL");
152
153 // advection operator
154 {
155 hyperdeal::ScopedTimerWrapper timer(timers, "advection");
156 if (timers != nullptr)
157 timers->enter("advection");
158
159 data.loop(&This::local_apply_cell,
160 &This::local_apply_face,
161 &This::local_apply_boundary,
162 this,
163 dst,
164 src,
165 MatrixFree<dim_x, dim_v, Number, VectorizedArrayType>::
166 DataAccessOnFaces::values,
167 MatrixFree<dim_x, dim_v, Number, VectorizedArrayType>::
168 DataAccessOnFaces::values,
169 timers);
170
171 if (timers != nullptr)
172 timers->leave();
173 }
174
175 // inverse-mass matrix operator
176 {
177 hyperdeal::ScopedTimerWrapper timer(timers, "mass");
178
179 data.cell_loop(&This::local_apply_inverse_mass_matrix,
180 this,
181 dst,
182 dst);
183 }
184
185 if (timers != nullptr)
186 timers->leave();
187 }
188 else // ECL
189 {
190 if (timers != nullptr)
191 timers->enter("ECL");
192
193 // advection and inverse-mass matrix operator in one go
194 data.loop_cell_centric(
195 &This::local_apply_advect_and_inverse_mass_matrix<eval_level>,
196 this,
197 dst,
198 src,
199 eval_level == AdvectionOperationEvaluationLevel::all ?
200 MatrixFree<dim_x, dim_v, Number, VectorizedArrayType>::
201 DataAccessOnFaces::values :
202 MatrixFree<dim_x, dim_v, Number, VectorizedArrayType>::
203 DataAccessOnFaces::none,
204 timers);
205
206 if (timers != nullptr)
207 timers->leave();
208 }
209 }
210
211 private:
212 using ID =
213 typename MatrixFree<dim_x, dim_v, Number, VectorizedArrayType>::ID;
214
218 template <AdvectionOperationEvaluationLevel eval_level =
219 AdvectionOperationEvaluationLevel::all>
220 void
221 local_apply_advect_and_inverse_mass_matrix(
222 const MatrixFree<dim_x, dim_v, Number, VectorizedArrayType> &data,
223 VectorType & dst,
224 const VectorType & src,
225 const ID cell)
226 {
227 (void)data;
228
229 auto &phi = *this->phi_cell;
230 auto &phi_m = *this->phi_face_m;
231 auto &phi_p = *this->phi_face_p;
232 auto &phi_inv = *this->phi_cell_inv;
233
234 // get data and scratch
235 VNumber *data_ptr = phi.get_data_ptr();
236 VNumber *data_ptr1 = phi_m.get_data_ptr();
237 VNumber *data_ptr2 = phi_p.get_data_ptr();
238 VNumber *data_ptr_inv = phi_inv.get_data_ptr();
239
240 // initialize tensor product kernels
241 const dealii::internal::EvaluatorTensorProduct<tensorproduct,
242 dim,
243 degree + 1,
244 n_points,
245 VNumber,
246 Number>
247 eval(*phi.get_shape_values(),
248 dealii::AlignedVector<Number>(),
249 dealii::AlignedVector<Number>());
250
251 const dealii::internal::EvaluatorTensorProduct<tensorproduct,
252 dim - 1,
253 degree + 1,
254 n_points,
255 VNumber,
256 Number>
257 eval_face(*phi.get_shape_values(),
258 dealii::AlignedVector<Number>(),
259 dealii::AlignedVector<Number>());
260
261 const dealii::internal::EvaluatorTensorProduct<tensorproduct,
262 dim,
263 n_points,
264 n_points,
265 VNumber,
266 Number>
267 eval_(dealii::AlignedVector<Number>(),
268 *phi.get_shape_gradients(),
269 dealii::AlignedVector<Number>());
270
271 const dealii::internal::EvaluatorTensorProduct<tensorproduct,
272 dim,
273 degree + 1,
274 n_points,
275 VNumber,
276 Number>
277 eval_inv(dealii::AlignedVector<Number>(),
278 dealii::AlignedVector<Number>(),
279 data.get_matrix_free_x()
280 .get_shape_info()
281 .data[0]
282 .inverse_shape_values_eo);
283
284 // clang-format off
285
286 // 1) advection: cell contribution
287 {
288 this->velocity_field->reinit(cell);
289
290 // load from global structure
291 phi.reinit(cell);
292 phi.read_dof_values(src);
293
294 if(do_collocation == false)
295 {
296 if(degree + 1 == n_points)
297 {
298 if (dim >= 1) eval.template values<0, true, false>(data_ptr, data_ptr);
299 if (dim >= 2) eval.template values<1, true, false>(data_ptr, data_ptr);
300 if (dim >= 3) eval.template values<2, true, false>(data_ptr, data_ptr);
301 if (dim >= 4) eval.template values<3, true, false>(data_ptr, data_ptr);
302 if (dim >= 5) eval.template values<4, true, false>(data_ptr, data_ptr);
303 if (dim >= 6) eval.template values<5, true, false>(data_ptr, data_ptr);
304 }
305 else
306 {
307 dealii::internal::FEEvaluationImplBasisChange<tensorproduct,
308 dealii::internal::EvaluatorQuantity::value,
309 dim,
310 degree + 1,
311 n_points>::
312 do_forward(1, data.get_matrix_free_x().get_shape_info().data.front().shape_values_eo,
313 data_ptr,
314 data_ptr);
315 }
316 }
317
318 // copy quadrature values into buffer
319 VNumber *buffer = phi_cell_inv->get_data_ptr();
320
321 if(eval_level != AdvectionOperationEvaluationLevel::cell)
322 for (auto i = 0u; i < dealii::Utilities::pow<unsigned int>(n_points, dim); i++)
323 buffer[i] = data_ptr[i];
324
325 // x-space
326 {
327 dealii::AlignedVector<VNumber> scratch_data_array;
328 scratch_data_array.resize_fast(dealii::Utilities::pow(n_points, dim) * dim_x);
329 VNumber *tempp = scratch_data_array.begin();
330
331 if(factor_skew != 0.0)
332 {
333 if (dim_x >= 1) eval_.template gradients<0, true, false>(buffer, tempp + dealii::Utilities::pow(n_points, dim) * 0);
334 if (dim_x >= 2) eval_.template gradients<1, true, false>(buffer, tempp + dealii::Utilities::pow(n_points, dim) * 1);
335 if (dim_x >= 3) eval_.template gradients<2, true, false>(buffer, tempp + dealii::Utilities::pow(n_points, dim) * 2);
336 }
337
338 for (auto qv = 0u, q = 0u; qv < dealii::Utilities::pow<unsigned int>(n_points, dim_v); ++qv)
339 for (auto qx = 0u; qx < dealii::Utilities::pow<unsigned int>(n_points, dim_x); ++qx, ++q)
340 {
341 VNumber grad_in[dim_x];
342 const auto vel = velocity_field->evaluate_x(q, qx, qv);
343
344 if(factor_skew != 0.0)
345 phi.template submit_value<false>(data_ptr, - factor_skew * (phi.get_gradient_x(tempp, q, qx, qv) * vel), q, qx, qv);
346
347 if (factor_skew != 1.0)
348 {
349 for (int d = 0; d < dim_x; d++)
350 grad_in[d] = (1.0-factor_skew) * buffer[q] * vel[d];
351 phi.submit_gradient_x(tempp, grad_in, q, qx, qv);
352 }
353 }
354
355 if(factor_skew != 1.0)
356 {
357 if (dim_x >= 1 && (factor_skew != 0.0))
358 eval_.template gradients<0, false, true >(tempp + dealii::Utilities::pow(n_points, dim) * 0, data_ptr);
359 else if (dim_x >= 1)
360 eval_.template gradients<0, false, false>(tempp + dealii::Utilities::pow(n_points, dim) * 0, data_ptr);
361
362 if (dim_x >= 2) eval_.template gradients<1, false, true >(tempp + dealii::Utilities::pow(n_points, dim) * 1, data_ptr);
363 if (dim_x >= 3) eval_.template gradients<2, false, true >(tempp + dealii::Utilities::pow(n_points, dim) * 2, data_ptr);
364 }
365 }
366
367 // v-space
368 {
369 dealii::AlignedVector<VNumber> scratch_data_array;
370 scratch_data_array.resize_fast(dealii::Utilities::pow(n_points, dim) * dim_v);
371 VNumber *tempp = scratch_data_array.begin();
372
373 if(factor_skew != 0.0)
374 {
375 if (dim_v >= 1) eval_.template gradients<0 + dim_x, true, false>(buffer, tempp + dealii::Utilities::pow(n_points, dim) * 0);
376 if (dim_v >= 2) eval_.template gradients<1 + dim_x, true, false>(buffer, tempp + dealii::Utilities::pow(n_points, dim) * 1);
377 if (dim_v >= 3) eval_.template gradients<2 + dim_x, true, false>(buffer, tempp + dealii::Utilities::pow(n_points, dim) * 2);
378 }
379
380 for (auto qv = 0u, q = 0u; qv < dealii::Utilities::pow<unsigned int>(n_points, dim_v); ++qv)
381 for (auto qx = 0u; qx < dealii::Utilities::pow<unsigned int>(n_points, dim_x); ++qx, ++q)
382 {
383 VNumber grad_in[dim_v];
384 const auto vel = velocity_field->evaluate_v(q, qx, qv);
385
386 if(factor_skew != 0.0)
387 phi.template submit_value<true>(data_ptr, - factor_skew * (phi.get_gradient_v(tempp, q, qx, qv) * vel), q, qx, qv);
388
389 if (factor_skew != 1.0)
390 {
391 for (int d = 0; d < dim_v; d++)
392 grad_in[d] = (1.0-factor_skew) * buffer[q] * vel[d];
393 phi.submit_gradient_v(tempp, grad_in, q, qx, qv);
394 }
395 }
396
397 if(factor_skew != 1.0)
398 {
399 if (dim_v >= 1) eval_.template gradients<0 + dim_x, false, true>(tempp + dealii::Utilities::pow(n_points, dim) * 0, data_ptr);
400 if (dim_v >= 2) eval_.template gradients<1 + dim_x, false, true>(tempp + dealii::Utilities::pow(n_points, dim) * 1, data_ptr);
401 if (dim_v >= 3) eval_.template gradients<2 + dim_x, false, true>(tempp + dealii::Utilities::pow(n_points, dim) * 2, data_ptr);
402 }
403 }
404 }
405
406 // 2) advection: faces
407 if(eval_level != AdvectionOperationEvaluationLevel::cell)
408 for (auto face = 0u; face < dim * 2; face++)
409 {
410 this->velocity_field->reinit_face(cell, face);
411
412 // load negative side from buffer
413 phi_m.reinit(cell, face);
414
415 const auto bid = data.get_faces_by_cells_boundary_id(cell, face);
416
417 // load positive side from global structure
418 if(bid == dealii::numbers::internal_face_boundary_id)
419 {
420 phi_p.reinit(cell, face);
421
422 if(eval_level == AdvectionOperationEvaluationLevel::all)
423 phi_p.read_dof_values(src);
424 }
425
426 if(do_collocation == false)
427 {
428 hyperdeal::internal::FEFaceNormalEvaluationImpl<dim_x, dim_v, n_points - 1, Number>::template interpolate_quadrature<true, false>(1, dealii::EvaluationFlags::values, data.get_matrix_free_x().get_shape_info(), /*out=*/data_ptr_inv, /*in=*/data_ptr1, face);
429
430 if(degree + 1 == n_points)
431 {
432 if (dim >= 2) eval_face.template values<0, true, false>(data_ptr2, data_ptr2);
433 if (dim >= 3) eval_face.template values<1, true, false>(data_ptr2, data_ptr2);
434 if (dim >= 4) eval_face.template values<2, true, false>(data_ptr2, data_ptr2);
435 if (dim >= 5) eval_face.template values<3, true, false>(data_ptr2, data_ptr2);
436 if (dim >= 6) eval_face.template values<4, true, false>(data_ptr2, data_ptr2);
437 }
438 else
439 {
440 dealii::internal::FEEvaluationImplBasisChange<tensorproduct,
441 dealii::internal::EvaluatorQuantity::value,
442 dim - 1,
443 degree + 1,
444 n_points>::
445 do_forward(1 ,data.get_matrix_free_x().get_shape_info().data.front().shape_values_eo,
446 data_ptr2,
447 data_ptr2);
448 }
449 }
450 else
451 {
452 phi_m.read_dof_values_from_buffer(this->phi_cell_inv->get_data_ptr());
453 }
454
455 if(bid == dealii::numbers::internal_face_boundary_id)
456 {
457 if (face < dim_x * 2)
458 {
459 for (unsigned int qv = 0, q = 0; qv < dealii::Utilities::pow<unsigned int>(n_points, dim_v); ++qv)
460 for (unsigned int qx = 0; qx < dealii::Utilities::pow<unsigned int>(n_points, dim_x - 1); ++qx, ++q)
461 {
462 const VectorizedArrayType u_minus = data_ptr1[q];
463 const VectorizedArrayType u_plus = data_ptr2[q];
464 const VectorizedArrayType normal_times_speed = velocity_field->evaluate_face_x(q, qx, qv) * phi_m.get_normal_vector_x(qx);
465 const VectorizedArrayType flux_times_normal_of_minus = 0.5 * ((u_minus + u_plus) * normal_times_speed + std::abs(normal_times_speed) * (u_minus - u_plus)) * alpha;
466
467 phi_m.template submit_value<ID::SpaceType::X>(data_ptr1, flux_times_normal_of_minus - factor_skew*u_minus*normal_times_speed, q, qx, qv);
468 }
469 }
470 else
471 {
472 for (unsigned int qv = 0, q = 0; qv < dealii::Utilities::pow<unsigned int>(n_points, dim_v - 1); ++qv)
473 for (unsigned int qx = 0; qx < dealii::Utilities::pow<unsigned int>(n_points, dim_x); ++qx, ++q)
474 {
475 const VectorizedArrayType u_minus = data_ptr1[q];
476 const VectorizedArrayType u_plus = data_ptr2[q];
477 const VectorizedArrayType normal_times_speed = velocity_field->evaluate_face_v(q, qx, qv) * phi_m.get_normal_vector_v(qv);
478 const VectorizedArrayType flux_times_normal_of_minus = 0.5 * ((u_minus + u_plus) * normal_times_speed + std::abs(normal_times_speed) * (u_minus - u_plus)) * alpha;
479
480 phi_m.template submit_value<ID::SpaceType::V>(data_ptr1, flux_times_normal_of_minus - factor_skew*u_minus*normal_times_speed, q, qx, qv);
481 }
482 }
483 }
484 else
485 {
486 const auto boundary_pair = boundary_descriptor->get_boundary(bid);
487
488 if (face < dim_x * 2)
489 {
490 for (unsigned int qv = 0, q = 0; qv < dealii::Utilities::pow<unsigned int>(n_points, dim_v); ++qv)
491 for (unsigned int qx = 0; qx < dealii::Utilities::pow<unsigned int>(n_points, dim_x - 1); ++qx, ++q)
492 {
493 const VectorizedArrayType u_minus = data_ptr1[q];
494 const VectorizedArrayType u_plus = boundary_pair.first == BoundaryType::DirichletHomogenous ?
495 (-u_minus) :
496 (-u_minus + 2.0 * hyperdeal::MatrixFreeTools::evaluate_scalar_function(phi_m.template get_quadrature_point<ID::SpaceType::X>(qx, qv), *boundary_pair.second, phi_m.n_vectorization_lanes_filled()));
497
498 const VectorizedArrayType normal_times_speed = velocity_field->evaluate_face_x(q, qx, qv) * phi_m.get_normal_vector_x(qx);
499 const VectorizedArrayType flux_times_normal_of_minus = 0.5 * ((u_minus + u_plus) * normal_times_speed + std::abs(normal_times_speed) * (u_minus - u_plus)) * alpha;
500
501 phi_m.template submit_value<ID::SpaceType::X>(data_ptr1, flux_times_normal_of_minus - factor_skew*u_minus*normal_times_speed, q, qx, qv);
502 }
503 }
504 else
505 {
506 for (unsigned int qv = 0, q = 0; qv < dealii::Utilities::pow<unsigned int>(n_points, dim_v - 1); ++qv)
507 for (unsigned int qx = 0; qx < dealii::Utilities::pow<unsigned int>(n_points, dim_x); ++qx, ++q)
508 {
509 const VectorizedArrayType u_minus = data_ptr1[q];
510 const VectorizedArrayType u_plus = boundary_pair.first == BoundaryType::DirichletHomogenous ?
511 (-u_minus) :
512 (-u_minus + 2.0 * hyperdeal::MatrixFreeTools::evaluate_scalar_function(phi_m.template get_quadrature_point<ID::SpaceType::V>(qx, qv), *boundary_pair.second, phi_m.n_vectorization_lanes_filled()));
513
514 const VectorizedArrayType normal_times_speed = velocity_field->evaluate_face_v(q, qx, qv) * phi_m.get_normal_vector_v(qv);
515 const VectorizedArrayType flux_times_normal_of_minus = 0.5 * ((u_minus + u_plus) * normal_times_speed + std::abs(normal_times_speed) * (u_minus - u_plus)) * alpha;
516
517 phi_m.template submit_value<ID::SpaceType::V>(data_ptr1, flux_times_normal_of_minus - factor_skew*u_minus*normal_times_speed, q, qx, qv);
518 }
519 }
520 }
521
522 if(do_collocation == false)
523 hyperdeal::internal::FEFaceNormalEvaluationImpl<dim_x, dim_v, n_points - 1, Number>::template interpolate_quadrature<false, true>(1, dealii::EvaluationFlags::values, data.get_matrix_free_x().get_shape_info(), /*out=*/data_ptr1, /*in=*/data_ptr, face);
524 else
525 phi_m.distribute_to_buffer(this->phi_cell->get_data_ptr());
526 }
527
528 // 3) inverse mass matrix
529 {
530 phi_inv.reinit(cell);
531
532 for (auto qv = 0u, q = 0u; qv < dealii::Utilities::pow<unsigned int>(n_points, dim_v); ++qv)
533 for (auto qx = 0u; qx < dealii::Utilities::pow<unsigned int>(n_points, dim_x); ++qx, ++q)
534 phi_inv.submit_inv(data_ptr, q, qx, qv);
535
536 if(do_collocation == false)
537 {
538 if(degree + 1 == n_points)
539 {
540 if (dim >= 6) eval_inv.template hessians<5, false, false>(data_ptr, data_ptr);
541 if (dim >= 5) eval_inv.template hessians<4, false, false>(data_ptr, data_ptr);
542 if (dim >= 4) eval_inv.template hessians<3, false, false>(data_ptr, data_ptr);
543 if (dim >= 3) eval_inv.template hessians<2, false, false>(data_ptr, data_ptr);
544 if (dim >= 2) eval_inv.template hessians<1, false, false>(data_ptr, data_ptr);
545 if (dim >= 1) eval_inv.template hessians<0, false, false>(data_ptr, data_ptr);
546 }
547 else
548 {
549 dealii::internal::FEEvaluationImplBasisChange<tensorproduct,
550 dealii::internal::EvaluatorQuantity::hessian,
551 dim,
552 degree + 1,
553 n_points>::
554 do_backward(1, data.get_matrix_free_x().get_shape_info().data.front().inverse_shape_values_eo,
555 false,
556 data_ptr,
557 data_ptr);
558 }
559 }
560
561 // write into global structure back
562 phi.set_dof_values(dst);
563 }
564
565 // clang-format on
566 }
567
568
569
573 void
574 local_apply_inverse_mass_matrix(
575 const MatrixFree<dim_x, dim_v, Number, VectorizedArrayType> &data,
576 VectorType & dst,
577 const VectorType & src,
578 const ID cell)
579 {
580 (void)data;
581
582 auto &phi_inv = *this->phi_cell_inv_co;
583
584 // get data and scratch
585 VectorizedArrayType *data_ptr = phi_inv.get_data_ptr();
586
587 // load from global structure
588 phi_inv.reinit(cell);
589 phi_inv.read_dof_values(src);
590
591 // initialize tensor product kernels
592 const dealii::internal::EvaluatorTensorProduct<tensorproduct,
593 dim,
594 degree + 1,
595 degree + 1,
596 VectorizedArrayType,
597 Number>
598 eval_inv(dealii::AlignedVector<Number>(),
599 dealii::AlignedVector<Number>(),
600 *phi_inv.get_inverse_shape());
601
602 // clang-format off
603
604 if(do_collocation == false)
605 {
606 if (dim >= 1) eval_inv.template hessians<0, true, false>(data_ptr, data_ptr);
607 if (dim >= 2) eval_inv.template hessians<1, true, false>(data_ptr, data_ptr);
608 if (dim >= 3) eval_inv.template hessians<2, true, false>(data_ptr, data_ptr);
609 if (dim >= 4) eval_inv.template hessians<3, true, false>(data_ptr, data_ptr);
610 if (dim >= 5) eval_inv.template hessians<4, true, false>(data_ptr, data_ptr);
611 if (dim >= 6) eval_inv.template hessians<5, true, false>(data_ptr, data_ptr);
612 }
613
614 for (auto qv = 0u, q = 0u; qv < dealii::Utilities::pow<unsigned int>(degree + 1, dim_v); ++qv)
615 for (auto qx = 0u; qx < dealii::Utilities::pow<unsigned int>(degree + 1, dim_x); ++qx, ++q)
616 phi_inv.submit_inv(data_ptr, q, qx, qv);
617
618 if(do_collocation == false)
619 {
620 if (dim >= 6) eval_inv.template hessians<5, false, false>(data_ptr, data_ptr);
621 if (dim >= 5) eval_inv.template hessians<4, false, false>(data_ptr, data_ptr);
622 if (dim >= 4) eval_inv.template hessians<3, false, false>(data_ptr, data_ptr);
623 if (dim >= 3) eval_inv.template hessians<2, false, false>(data_ptr, data_ptr);
624 if (dim >= 2) eval_inv.template hessians<1, false, false>(data_ptr, data_ptr);
625 if (dim >= 1) eval_inv.template hessians<0, false, false>(data_ptr, data_ptr);
626 }
627
628 // clang-format on
629
630 // write into global structure back
631 phi_inv.set_dof_values(dst);
632 }
633
634
635
639 void
640 local_apply_cell(
641 const MatrixFree<dim_x, dim_v, Number, VectorizedArrayType> &data,
642 VectorType & dst,
643 const VectorType & src,
644 const ID cell)
645 {
646 (void)data;
647
648 auto &phi = *this->phi_cell;
649
650 // get data and scratch
651 VNumber *data_ptr = phi.get_data_ptr();
652
653 // initialize tensor product kernels
654 const dealii::internal::EvaluatorTensorProduct<tensorproduct,
655 dim,
656 degree + 1,
657 n_points,
658 VNumber,
659 Number>
660 eval(*phi.get_shape_values(),
661 dealii::AlignedVector<Number>(),
662 dealii::AlignedVector<Number>());
663
664 const dealii::internal::EvaluatorTensorProduct<tensorproduct,
665 dim,
666 n_points,
667 n_points,
668 VNumber,
669 Number>
670 eval_(dealii::AlignedVector<Number>(),
671 *phi.get_shape_gradients(),
672 dealii::AlignedVector<Number>());
673
674 // clang-format off
675
676 this->velocity_field->reinit(cell);
677
678 // load from global structure
679 phi.reinit(cell);
680 phi.read_dof_values(src);
681
682 if(do_collocation == false)
683 {
684 if(degree + 1 == n_points)
685 {
686 if (dim >= 1) eval.template values<0, true, false>(data_ptr, data_ptr);
687 if (dim >= 2) eval.template values<1, true, false>(data_ptr, data_ptr);
688 if (dim >= 3) eval.template values<2, true, false>(data_ptr, data_ptr);
689 if (dim >= 4) eval.template values<3, true, false>(data_ptr, data_ptr);
690 if (dim >= 5) eval.template values<4, true, false>(data_ptr, data_ptr);
691 if (dim >= 6) eval.template values<5, true, false>(data_ptr, data_ptr);
692 }
693 else
694 {
695 dealii::internal::FEEvaluationImplBasisChange<tensorproduct,
696 dealii::internal::EvaluatorQuantity::value,
697 dim,
698 degree + 1,
699 n_points>::
700 do_forward(1, data.get_matrix_free_x().get_shape_info().data.front().shape_values_eo,
701 data_ptr,
702 data_ptr);
703 }
704 }
705
706 // copy quadrature values into buffer
707 VNumber *buffer = phi_cell_inv->get_data_ptr();
708 for (auto i = 0u; i < dealii::Utilities::pow<unsigned int>(n_points, dim); i++)
709 buffer[i] = data_ptr[i];
710
711 // x-space
712 {
713 dealii::AlignedVector<VNumber> scratch_data_array;
714 scratch_data_array.resize_fast(dealii::Utilities::pow(n_points, dim) * dim_x);
715 VNumber *tempp = scratch_data_array.begin();
716
717 if(factor_skew != 0.0)
718 {
719 if (dim_x >= 1) eval_.template gradients<0, true, false>(buffer, tempp + dealii::Utilities::pow(n_points, dim) * 0);
720 if (dim_x >= 2) eval_.template gradients<1, true, false>(buffer, tempp + dealii::Utilities::pow(n_points, dim) * 1);
721 if (dim_x >= 3) eval_.template gradients<2, true, false>(buffer, tempp + dealii::Utilities::pow(n_points, dim) * 2);
722 }
723
724 for (auto qv = 0u, q = 0u; qv < dealii::Utilities::pow<unsigned int>(n_points, dim_v); ++qv)
725 for (auto qx = 0u; qx < dealii::Utilities::pow<unsigned int>(n_points, dim_x); ++qx, ++q)
726 {
727 VNumber grad_in[dim_x];
728 const auto vel = velocity_field->evaluate_x(q, qx, qv);
729
730 if(factor_skew != 0.0)
731 phi.template submit_value<false>(data_ptr, - factor_skew * (phi.get_gradient_x(tempp, q, qx, qv) * vel), q, qx, qv);
732
733 if (factor_skew != 1.0)
734 {
735 for (int d = 0; d < dim_x; d++)
736 grad_in[d] = (1.0-factor_skew) * buffer[q] * vel[d];
737 phi.submit_gradient_x(tempp, grad_in, q, qx, qv);
738 }
739 }
740
741 if(factor_skew != 1.0)
742 {
743 if (dim_x >= 1 && (factor_skew != 0.0))
744 eval_.template gradients<0, false, true >(tempp + dealii::Utilities::pow(n_points, dim) * 0, data_ptr);
745 else if (dim_x >= 1)
746 eval_.template gradients<0, false, false>(tempp + dealii::Utilities::pow(n_points, dim) * 0, data_ptr);
747 if (dim_x >= 2) eval_.template gradients<1, false, true >(tempp + dealii::Utilities::pow(n_points, dim) * 1, data_ptr);
748 if (dim_x >= 3) eval_.template gradients<2, false, true >(tempp + dealii::Utilities::pow(n_points, dim) * 2, data_ptr);
749 }
750 }
751 // v-space
752 {
753 dealii::AlignedVector<VNumber> scratch_data_array;
754 scratch_data_array.resize_fast(dealii::Utilities::pow(n_points, dim) * dim_v);
755 VNumber *tempp = scratch_data_array.begin();
756
757 if(factor_skew != 0.0)
758 {
759 if (dim_v >= 1) eval_.template gradients<0 + dim_x, true, false>(buffer, tempp + dealii::Utilities::pow(n_points, dim) * 0);
760 if (dim_v >= 2) eval_.template gradients<1 + dim_x, true, false>(buffer, tempp + dealii::Utilities::pow(n_points, dim) * 1);
761 if (dim_v >= 3) eval_.template gradients<2 + dim_x, true, false>(buffer, tempp + dealii::Utilities::pow(n_points, dim) * 2);
762 }
763
764 for (auto qv = 0u, q = 0u; qv < dealii::Utilities::pow<unsigned int>(n_points, dim_v); ++qv)
765 for (auto qx = 0u; qx < dealii::Utilities::pow<unsigned int>(n_points, dim_x); ++qx, ++q)
766 {
767 VNumber grad_in[dim_v];
768 const auto vel = velocity_field->evaluate_v(q, qx, qv);
769
770 if(factor_skew != 0.0)
771 phi.template submit_value<true>(data_ptr, - factor_skew * (phi.get_gradient_v(tempp, q, qx, qv) * vel), q, qx, qv);
772
773 if (factor_skew != 1.0)
774 {
775 for (int d = 0; d < dim_v; d++)
776 grad_in[d] = (1.0-factor_skew) * buffer[q] * vel[d];
777 phi.submit_gradient_v(tempp, grad_in, q, qx, qv);
778 }
779 }
780
781 if(factor_skew != 1.0)
782 {
783 if (dim_v >= 1) eval_.template gradients<0 + dim_x, false, true>(tempp + dealii::Utilities::pow(n_points, dim) * 0, data_ptr);
784 if (dim_v >= 2) eval_.template gradients<1 + dim_x, false, true>(tempp + dealii::Utilities::pow(n_points, dim) * 1, data_ptr);
785 if (dim_v >= 3) eval_.template gradients<2 + dim_x, false, true>(tempp + dealii::Utilities::pow(n_points, dim) * 2, data_ptr);
786 }
787 }
788
789 if(do_collocation == false)
790 {
791 if(degree + 1 == n_points)
792 {
793 if(dim >= 6) eval.template values<5, false, false>(data_ptr, data_ptr);
794 if(dim >= 5) eval.template values<4, false, false>(data_ptr, data_ptr);
795 if(dim >= 4) eval.template values<3, false, false>(data_ptr, data_ptr);
796 if(dim >= 3) eval.template values<2, false, false>(data_ptr, data_ptr);
797 if(dim >= 2) eval.template values<1, false, false>(data_ptr, data_ptr);
798 if(dim >= 1) eval.template values<0, false, false>(data_ptr, data_ptr);
799 }
800 else
801 {
802 dealii::internal::FEEvaluationImplBasisChange<tensorproduct,
803 dealii::internal::EvaluatorQuantity::value,
804 dim,
805 degree + 1,
806 n_points>::
807 do_backward(1, data.get_matrix_free_x().get_shape_info().data.front().shape_values_eo,
808 false,
809 data_ptr,
810 data_ptr);
811 }
812 }
813
814 // clang-format on
815
816 // write into global structure back
817 phi.set_dof_values(dst);
818 }
819
820
821
825 void
826 local_apply_face(
827 const MatrixFree<dim_x, dim_v, Number, VectorizedArrayType> &data,
828 VectorType & dst,
829 const VectorType & src,
830 const ID face)
831 {
832 (void)data;
833
834 auto &phi_m = *this->phi_face_m;
835 auto &phi_p = *this->phi_face_p;
836
837 const dealii::internal::EvaluatorTensorProduct<tensorproduct,
838 dim - 1,
839 degree + 1,
840 n_points,
841 VNumber,
842 Number>
843 eval1(*phi_m.get_shape_values(),
844 dealii::AlignedVector<Number>(),
845 dealii::AlignedVector<Number>());
846 const dealii::internal::EvaluatorTensorProduct<tensorproduct,
847 dim - 1,
848 degree + 1,
849 n_points,
850 VNumber,
851 Number>
852 eval2(*phi_p.get_shape_values(),
853 dealii::AlignedVector<Number>(),
854 dealii::AlignedVector<Number>());
855
856 this->velocity_field->reinit_face(face);
857
858 // get data and scratch
859 VNumber *data_ptr1 = phi_m.get_data_ptr();
860 VNumber *data_ptr2 = phi_p.get_data_ptr();
861
862 // load from global structure
863 phi_m.reinit(face);
864 phi_p.reinit(face);
865
866 // clang-format off
867
868 phi_m.read_dof_values(src);
869
870 if(do_collocation == false)
871 {
872 if(degree + 1 == n_points)
873 {
874 if (dim >= 2) eval1.template values<0, true, false>(data_ptr1, data_ptr1);
875 if (dim >= 3) eval1.template values<1, true, false>(data_ptr1, data_ptr1);
876 if (dim >= 4) eval1.template values<2, true, false>(data_ptr1, data_ptr1);
877 if (dim >= 5) eval1.template values<3, true, false>(data_ptr1, data_ptr1);
878 if (dim >= 6) eval1.template values<4, true, false>(data_ptr1, data_ptr1);
879 }
880 else
881 {
882 dealii::internal::FEEvaluationImplBasisChange<tensorproduct,
883 dealii::internal::EvaluatorQuantity::value,
884 dim - 1,
885 degree + 1,
886 n_points>::
887 do_forward(1, data.get_matrix_free_x().get_shape_info().data.front().shape_values_eo,
888 data_ptr1,
889 data_ptr1);
890 }
891 }
892
893 phi_p.read_dof_values(src);
894
895 if(do_collocation == false)
896 {
897 if(degree + 1 == n_points)
898 {
899 if (dim >= 2) eval2.template values<0, true, false>(data_ptr2, data_ptr2);
900 if (dim >= 3) eval2.template values<1, true, false>(data_ptr2, data_ptr2);
901 if (dim >= 4) eval2.template values<2, true, false>(data_ptr2, data_ptr2);
902 if (dim >= 5) eval2.template values<3, true, false>(data_ptr2, data_ptr2);
903 if (dim >= 6) eval2.template values<4, true, false>(data_ptr2, data_ptr2);
904 }
905 else
906 {
907 dealii::internal::FEEvaluationImplBasisChange<tensorproduct,
908 dealii::internal::EvaluatorQuantity::value,
909 dim - 1,
910 degree + 1,
911 n_points>::
912 do_forward(1, data.get_matrix_free_x().get_shape_info().data.front().shape_values_eo,
913 data_ptr2,
914 data_ptr2);
915 }
916 }
917
918 if (face.type == ID::SpaceType::X)
919 {
920 for (unsigned int qv = 0, q = 0; qv < dealii::Utilities::pow<unsigned int>(n_points, dim_v); ++qv)
921 for (unsigned int qx = 0; qx < dealii::Utilities::pow<unsigned int>(n_points, dim_x - 1); ++qx, ++q)
922 {
923 const VectorizedArrayType u_minus = data_ptr1[q];
924 const VectorizedArrayType u_plus = data_ptr2[q];
925 const VectorizedArrayType normal_times_speed = velocity_field->evaluate_face_x(q, qx, qv) * phi_m.get_normal_vector_x(qx);
926 const VectorizedArrayType flux_times_normal_of_minus = 0.5 * ((u_minus + u_plus) * normal_times_speed + std::abs(normal_times_speed) * (u_minus - u_plus)) * alpha;
927
928 phi_m.template submit_value<ID::SpaceType::X>(data_ptr1, +flux_times_normal_of_minus - factor_skew*u_minus*normal_times_speed, q, qx, qv);
929 phi_p.template submit_value<ID::SpaceType::X>(data_ptr2, -flux_times_normal_of_minus + factor_skew*u_plus*normal_times_speed, q, qx, qv);
930 }
931 }
932 else
933 {
934 for (unsigned int qv = 0, q = 0; qv < dealii::Utilities::pow<unsigned int>(n_points, dim_v - 1); ++qv)
935 for (unsigned int qx = 0; qx < dealii::Utilities::pow<unsigned int>(n_points, dim_x); ++qx, ++q)
936 {
937 const VectorizedArrayType u_minus = data_ptr1[q];
938 const VectorizedArrayType u_plus = data_ptr2[q];
939 const VectorizedArrayType normal_times_speed = velocity_field->evaluate_face_v(q, qx, qv) * phi_m.get_normal_vector_v(qv);
940 const VectorizedArrayType flux_times_normal_of_minus = 0.5 * ((u_minus + u_plus) * normal_times_speed + std::abs(normal_times_speed) * (u_minus - u_plus)) * alpha;
941
942 phi_m.template submit_value<ID::SpaceType::V>(data_ptr1, +flux_times_normal_of_minus - factor_skew*u_minus*normal_times_speed, q, qx, qv);
943 phi_p.template submit_value<ID::SpaceType::V>(data_ptr2, -flux_times_normal_of_minus + factor_skew*u_plus*normal_times_speed, q, qx, qv);
944 }
945 }
946
947 if(do_collocation == false)
948 {
949 if(degree + 1 == n_points)
950 {
951 if (dim >= 6) eval1.template values<4, false, false>(data_ptr1, data_ptr1);
952 if (dim >= 5) eval1.template values<3, false, false>(data_ptr1, data_ptr1);
953 if (dim >= 4) eval1.template values<2, false, false>(data_ptr1, data_ptr1);
954 if (dim >= 3) eval1.template values<1, false, false>(data_ptr1, data_ptr1);
955 if (dim >= 2) eval1.template values<0, false, false>(data_ptr1, data_ptr1);
956 }
957 else
958 {
959 dealii::internal::FEEvaluationImplBasisChange<tensorproduct,
960 dealii::internal::EvaluatorQuantity::value,
961 dim - 1,
962 degree + 1,
963 n_points>::
964 do_backward(1, data.get_matrix_free_x().get_shape_info().data.front().shape_values_eo, false,
965 data_ptr1,
966 data_ptr1);
967 }
968 }
969
970 // write into global structure back
971 phi_m.distribute_local_to_global(dst);
972
973 if(do_collocation == false)
974 {
975 if(degree + 1 == n_points)
976 {
977 if (dim >= 6) eval2.template values<4, false, false>(data_ptr2, data_ptr2);
978 if (dim >= 5) eval2.template values<3, false, false>(data_ptr2, data_ptr2);
979 if (dim >= 4) eval2.template values<2, false, false>(data_ptr2, data_ptr2);
980 if (dim >= 3) eval2.template values<1, false, false>(data_ptr2, data_ptr2);
981 if (dim >= 2) eval2.template values<0, false, false>(data_ptr2, data_ptr2);
982 }
983 else
984 {
985 dealii::internal::FEEvaluationImplBasisChange<tensorproduct,
986 dealii::internal::EvaluatorQuantity::value,
987 dim - 1,
988 degree + 1,
989 n_points>::
990 do_backward(1, data.get_matrix_free_x().get_shape_info().data.front().shape_values_eo, false,
991 data_ptr2,
992 data_ptr2);
993 }
994 }
995
996 // clang-format on
997
998 // write into global structure back
999 phi_p.distribute_local_to_global(dst);
1000 }
1001
1002
1003
1009 void
1010 local_apply_boundary(
1011 const MatrixFree<dim_x, dim_v, Number, VectorizedArrayType> &data,
1012 VectorType & dst,
1013 const VectorType & src,
1014 const ID face)
1015 {
1016 (void)data;
1017
1018 const auto bid = data.get_boundary_id(face);
1019
1020 Assert(bid != dealii::numbers::internal_face_boundary_id,
1021 dealii::StandardExceptions::ExcInternalError());
1022
1023 const auto boundary_pair = boundary_descriptor->get_boundary(bid);
1024
1025 auto &phi_m = *this->phi_face_m;
1026
1027 const dealii::internal::EvaluatorTensorProduct<tensorproduct,
1028 dim - 1,
1029 degree + 1,
1030 n_points,
1031 VNumber,
1032 Number>
1033 eval1(*phi_m.get_shape_values(),
1034 dealii::AlignedVector<Number>(),
1035 dealii::AlignedVector<Number>());
1036
1037 this->velocity_field->reinit_face(face);
1038
1039 // get data and scratch
1040 VNumber *data_ptr1 = phi_m.get_data_ptr();
1041
1042 // load from global structure
1043 phi_m.reinit(face);
1044
1045 // clang-format off
1046
1047 phi_m.read_dof_values(src);
1048
1049 if(do_collocation == false)
1050 {
1051 if(degree + 1 == n_points)
1052 {
1053 if (dim >= 2) eval1.template values<0, true, false>(data_ptr1, data_ptr1);
1054 if (dim >= 3) eval1.template values<1, true, false>(data_ptr1, data_ptr1);
1055 if (dim >= 4) eval1.template values<2, true, false>(data_ptr1, data_ptr1);
1056 if (dim >= 5) eval1.template values<3, true, false>(data_ptr1, data_ptr1);
1057 if (dim >= 6) eval1.template values<4, true, false>(data_ptr1, data_ptr1);
1058 }
1059 else
1060 {
1061 dealii::internal::FEEvaluationImplBasisChange<tensorproduct,
1062 dealii::internal::EvaluatorQuantity::value,
1063 dim - 1,
1064 degree + 1,
1065 n_points>::
1066 do_forward(1, data.get_matrix_free_x().get_shape_info().data.front().shape_values_eo,
1067 data_ptr1,
1068 data_ptr1);
1069 }
1070 }
1071
1072 if (face.type == ID::SpaceType::X)
1073 {
1074 for (unsigned int qv = 0, q = 0; qv < dealii::Utilities::pow<unsigned int>(n_points, dim_v); ++qv)
1075 for (unsigned int qx = 0; qx < dealii::Utilities::pow<unsigned int>(n_points, dim_x - 1); ++qx, ++q)
1076 {
1077 const VectorizedArrayType u_minus = data_ptr1[q];
1078 const VectorizedArrayType u_plus = (boundary_pair.first == BoundaryType::DirichletHomogenous) ?
1079 (-u_minus) :
1080 (-u_minus + 2.0 * hyperdeal::MatrixFreeTools::evaluate_scalar_function(phi_m.template get_quadrature_point<ID::SpaceType::X>(qx, qv), *boundary_pair.second, phi_m.n_vectorization_lanes_filled()));
1081
1082 const VectorizedArrayType normal_times_speed = velocity_field->evaluate_face_x(q, qx, qv) * phi_m.get_normal_vector_x(qx);
1083 const VectorizedArrayType flux_times_normal_of_minus = 0.5 * ((u_minus + u_plus) * normal_times_speed + std::abs(normal_times_speed) * (u_minus - u_plus)) * alpha;
1084
1085 phi_m.template submit_value<ID::SpaceType::X>(data_ptr1, flux_times_normal_of_minus - factor_skew*u_minus*normal_times_speed, q, qx, qv);
1086 }
1087 }
1088 else
1089 {
1090 for (unsigned int qv = 0, q = 0; qv < dealii::Utilities::pow<unsigned int>(n_points, dim_v - 1); ++qv)
1091 for (unsigned int qx = 0; qx < dealii::Utilities::pow<unsigned int>(n_points, dim_x); ++qx, ++q)
1092 {
1093 const VectorizedArrayType u_minus = data_ptr1[q];
1094 const VectorizedArrayType u_plus = (boundary_pair.first == BoundaryType::DirichletHomogenous) ?
1095 (-u_minus) :
1096 (-u_minus + 2.0 * hyperdeal::MatrixFreeTools::evaluate_scalar_function(phi_m.template get_quadrature_point<ID::SpaceType::V>(qx, qv), *boundary_pair.second, phi_m.n_vectorization_lanes_filled()));
1097
1098 const VectorizedArrayType normal_times_speed = velocity_field->evaluate_face_v(q, qx, qv) * phi_m.get_normal_vector_v(qv);
1099 const VectorizedArrayType flux_times_normal_of_minus = 0.5 * ((u_minus + u_plus) * normal_times_speed + std::abs(normal_times_speed) * (u_minus - u_plus)) * alpha;
1100
1101 phi_m.template submit_value<ID::SpaceType::V>(data_ptr1, flux_times_normal_of_minus - factor_skew*u_minus*normal_times_speed, q, qx, qv);
1102 }
1103 }
1104
1105 if(do_collocation == false)
1106 {
1107 if(degree + 1 == n_points)
1108 {
1109 if (dim >= 6) eval1.template values<4, false, false>(data_ptr1, data_ptr1);
1110 if (dim >= 5) eval1.template values<3, false, false>(data_ptr1, data_ptr1);
1111 if (dim >= 4) eval1.template values<2, false, false>(data_ptr1, data_ptr1);
1112 if (dim >= 3) eval1.template values<1, false, false>(data_ptr1, data_ptr1);
1113 if (dim >= 2) eval1.template values<0, false, false>(data_ptr1, data_ptr1);
1114 }
1115 else
1116 {
1117 dealii::internal::FEEvaluationImplBasisChange<tensorproduct,
1118 dealii::internal::EvaluatorQuantity::value,
1119 dim - 1,
1120 degree + 1,
1121 n_points>::
1122 do_backward(1, data.get_matrix_free_x().get_shape_info().data.front().shape_values_eo,
1123 false,
1124 data_ptr1,
1125 data_ptr1);
1126 }
1127 }
1128
1129 // write into global structure back
1130 phi_m.distribute_local_to_global(dst);
1131 }
1132
1133 const MatrixFree<dim_x, dim_v, Number, VectorizedArrayType> &data;
1134 DynamicConvergenceTable & table;
1135
1136 // clang-format off
1137 std::shared_ptr<FEEvaluation<dim_x, dim_v, degree, n_points, Number, VNumber>> phi_cell;
1138 std::shared_ptr<FEEvaluationInverse<dim_x, dim_v, degree, n_points, Number, VNumber>> phi_cell_inv;
1139 std::shared_ptr<FEEvaluationInverse<dim_x, dim_v, degree, degree + 1, Number, VNumber>> phi_cell_inv_co;
1140 std::shared_ptr<FEFaceEvaluation<dim_x, dim_v, degree, n_points, Number, VNumber>> phi_face_m;
1141 std::shared_ptr<FEFaceEvaluation<dim_x, dim_v, degree, n_points, Number, VNumber>> phi_face_p;
1142 // clang-format on
1143
1144 std::shared_ptr<BoundaryDescriptor<dim, Number>> boundary_descriptor;
1145 std::shared_ptr<VelocityField> velocity_field;
1146
1147 bool do_collocation;
1148
1149 const double alpha = 1.0;
1150
1151 // skew factor: conservative (skew=0) and convective (skew=1)
1152 double factor_skew = 0.0;
1153 };
1154 } // namespace advection
1155} // namespace hyperdeal
1156
1157
1158#endif
hyperdeal::DynamicConvergenceTable
Definition dynamic_convergence_table.h:33
hyperdeal::FEEvaluationInverse
Definition fe_evaluation_cell_inverse.h:42
hyperdeal::FEEvaluation
Definition fe_evaluation_cell.h:47
hyperdeal::FEFaceEvaluation
Definition fe_evaluation_face.h:47
hyperdeal::MatrixFree
Definition matrix_free.h:40
hyperdeal::MatrixFree::get_faces_by_cells_boundary_id
dealii::types::boundary_id get_faces_by_cells_boundary_id(const TensorID &macro_cell, const unsigned int face_number) const
Definition matrix_free.templates.h:1781
hyperdeal::MatrixFree::get_matrix_free_x
const dealii::MatrixFree< dim_x, Number, VectorizedArrayTypeX > & get_matrix_free_x() const
Definition matrix_free.templates.h:1848
hyperdeal::ScopedTimerWrapper
Definition timers.h:270
hyperdeal::Timers
Definition timers.h:97
hyperdeal::advection::AdvectionOperation
Definition advection_operation.h:57
hyperdeal::advection::AdvectionOperation::AdvectionOperation
AdvectionOperation(const MatrixFree< dim_x, dim_v, Number, VectorizedArrayType > &data, DynamicConvergenceTable &table)
Definition advection_operation.h:85
hyperdeal::advection::AdvectionOperation::apply
void apply(VectorType &dst, const VectorType &src, const Number time, Timers *timers=nullptr)
Definition advection_operation.h:137
hyperdeal::advection::AdvectionOperation::reinit
void reinit(std::shared_ptr< BoundaryDescriptor< dim, Number > > boundary_descriptor, std::shared_ptr< VelocityField > velocity_field, const AdvectionOperationParamters additional_data)
Definition advection_operation.h:98
hyperdeal::TensorID
Definition id.h:27
hyperdeal::advection::AdvectionOperationParamters
Definition advection_operation_parameters.h:30
hyperdeal::advection::BoundaryDescriptor
Definition boundary_descriptor.h:43
hyperdeal::internal::FEFaceNormalEvaluationImpl
Definition evaluation_kernels.h:222
Generated by doxygen 1.12.0