ryujin 2.1.1 revision c38afcfdec63b34961924f05408760247496a1f0
parabolic_solver_gmg_operators.h
Go to the documentation of this file.
1//
2// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
3// Copyright (C) 2023 by the ryujin authors
4//
5
6#pragma once
7
8#include <offline_data.h>
9#include <openmp.h>
10#include <simd.h>
11
12#include "../euler/hyperbolic_system.h"
13#include "parabolic_system.h"
14
15#include <deal.II/base/vectorization.h>
16#include <deal.II/lac/diagonal_matrix.h>
17#include <deal.II/matrix_free/fe_evaluation.h>
18#include <deal.II/multigrid/mg_base.h>
19#include <deal.II/multigrid/mg_transfer_matrix_free.h>
20
21/*
22 * FIXME: generalize and make these operators equation independent and
23 * refactor into ../parabolic_module_gmg_operators.h
24 */
25
26namespace ryujin
27{
28 namespace NavierStokes
29 {
31
40 template <int dim, typename Number>
42 {
43 public:
47 using vector_type = dealii::LinearAlgebra::distributed::Vector<Number>;
48
53 dealii::LinearAlgebra::distributed::BlockVector<Number>;
54
58 DiagonalMatrix() = default;
59
64 void reinit(const vector_type &lumped_mass_matrix,
65 const vector_type &density,
66 const dealii::AffineConstraints<Number> &affine_constraints)
67 {
68 diagonal.reinit(density, true);
69
70 DEAL_II_OPENMP_SIMD_PRAGMA
71 for (unsigned int i = 0;
72 i < density.get_partitioner()->locally_owned_size();
73 ++i) {
74 diagonal.local_element(i) =
75 Number(1.0) /
76 (density.local_element(i) * lumped_mass_matrix.local_element(i));
77 }
78
79 /*
80 * Fix up diagonal entries for constrained degrees of freedom due to
81 * periodic boundary conditions.
82 */
83 affine_constraints.set_zero(diagonal);
84 }
85
90 {
91 return diagonal;
92 }
93
98 {
99 return diagonal_block;
100 }
101
105 void vmult(vector_type &dst, const vector_type &src) const
106 {
107 AssertDimension(diagonal_block.size(), 0);
108 DEAL_II_OPENMP_SIMD_PRAGMA
109 for (unsigned int i = 0;
110 i < diagonal.get_partitioner()->locally_owned_size();
111 ++i)
112 dst.local_element(i) =
113 diagonal.local_element(i) * src.local_element(i);
114 }
115
119 void vmult(block_vector_type &dst, const block_vector_type &src) const
120 {
121 AssertDimension(dim, dst.n_blocks());
122 AssertDimension(dim, src.n_blocks());
123 if (diagonal_block.size() == 0) {
124 DEAL_II_OPENMP_SIMD_PRAGMA
125 for (unsigned int i = 0;
126 i < diagonal.get_partitioner()->locally_owned_size();
127 ++i)
128 for (unsigned int d = 0; d < dim; ++d)
129 dst.block(d).local_element(i) =
130 diagonal.local_element(i) * src.block(d).local_element(i);
131 } else
132 for (unsigned int d = 0; d < dim; ++d) {
133 DEAL_II_OPENMP_SIMD_PRAGMA
134 for (unsigned int i = 0;
135 i < src.block(d).get_partitioner()->locally_owned_size();
136 ++i)
137 dst.block(d).local_element(i) =
138 diagonal_block.block(d).local_element(i) *
139 src.block(d).local_element(i);
140 }
141 }
142
143 private:
144 vector_type diagonal;
145 block_vector_type diagonal_block;
146 };
147
148
155 template <int dim, typename Number, typename Number2>
156 class VelocityMatrix : public dealii::Subscriptor
157 {
158 public:
159 // FIXME: refactor
160 static constexpr unsigned int order_fe = 1;
161 static constexpr unsigned int order_quad = 2;
162
163 using vector_type = dealii::LinearAlgebra::distributed::Vector<Number>;
165 dealii::LinearAlgebra::distributed::BlockVector<Number>;
166
167 VelocityMatrix() = default;
168
170 const ParabolicSystem &parabolic_system,
171 const OfflineData<dim, Number2> &offline_data,
172 const dealii::MatrixFree<dim, Number> &matrix_free,
173 const dealii::LinearAlgebra::distributed::Vector<Number> &density,
174 const Number theta_x_tau,
175 const unsigned int level = dealii::numbers::invalid_unsigned_int)
176 {
177 parabolic_system_ = &parabolic_system;
178 offline_data_ = &offline_data;
179 matrix_free_ = &matrix_free;
180 density_ = &density;
181 theta_x_tau_ = theta_x_tau;
182 level_ = level;
183 }
184
185 void Tvmult(block_vector_type &dst, const block_vector_type &src) const
186 {
187 vmult(dst, src);
188 }
189
190 void vmult(block_vector_type &dst, const block_vector_type &src) const
191 {
192 /* Apply action of m_i rho_i V_i: */
193
194 using VA = dealii::VectorizedArray<Number>;
195 constexpr auto simd_length = VA::size();
196
197 const vector_type *lumped_mass_matrix = nullptr;
198 if constexpr (std::is_same<Number, Number2>::value) {
199 if constexpr (std::is_same<Number, float>::value) {
200 if (level_ == dealii::numbers::invalid_unsigned_int)
201 lumped_mass_matrix = &offline_data_->lumped_mass_matrix();
202 else
203 lumped_mass_matrix =
204 &offline_data_->level_lumped_mass_matrix()[level_];
205 } else {
206 Assert(level_ == dealii::numbers::invalid_unsigned_int,
207 dealii::ExcInternalError());
208 lumped_mass_matrix = &offline_data_->lumped_mass_matrix();
209 }
210 } else
211 lumped_mass_matrix =
212 &offline_data_->level_lumped_mass_matrix()[level_];
213
214 const unsigned int n_owned =
215 lumped_mass_matrix->get_partitioner()->locally_owned_size();
216 const unsigned int size_regular = n_owned / simd_length * simd_length;
217
219
221 for (unsigned int i = 0; i < size_regular; i += simd_length) {
222 const auto m_i = get_entry<VA>(*lumped_mass_matrix, i);
223 const auto rho_i = get_entry<VA>(*density_, i);
224 for (unsigned int d = 0; d < dim; ++d) {
225 const auto temp = get_entry<VA>(src.block(d), i);
226 write_entry<VA>(dst.block(d), m_i * rho_i * temp, i);
227 }
228 }
229
231
232 for (unsigned int i = size_regular; i < n_owned; ++i) {
233 const auto m_i = lumped_mass_matrix->local_element(i);
234 const auto rho_i = density_->local_element(i);
235
236 for (unsigned int d = 0; d < dim; ++d) {
237 const auto temp = src.block(d).local_element(i);
238 dst.block(d).local_element(i) = m_i * rho_i * temp;
239 }
240 }
241
242 /* Apply action of stress tensor: + theta * \sum_j B_ij V_j: */
243
244 const auto integrator = [this](const auto &data,
245 auto &dst,
246 const auto &src,
247 const auto range) {
248 dealii::FEEvaluation<dim, order_fe, order_quad, dim, Number> velocity(
249 data);
250
251 for (unsigned int cell = range.first; cell < range.second; ++cell) {
252 velocity.reinit(cell);
253 velocity.read_dof_values(src);
254 apply_local_operator(velocity);
255 velocity.distribute_local_to_global(dst);
256 }
257 };
258
259 matrix_free_->template cell_loop<block_vector_type, block_vector_type>(
260 integrator, dst, src, /* zero destination */ false);
261
262 /* (5.4a) Fix up constrained degrees of freedom: */
263
264 const auto &boundary_map =
265 level_ == dealii::numbers::invalid_unsigned_int
266 ? offline_data_->boundary_map()
267 : offline_data_->level_boundary_map()[level_];
268
269 for (auto entry : boundary_map) {
270 // [i, normal, normal_mass, boundary_mass, id, position] = entry
271 const auto i = std::get<0>(entry);
272 if (i >= n_owned)
273 continue;
274
275 const dealii::Tensor<1, dim, Number> normal = std::get<1>(entry);
276 const auto id = std::get<4>(entry);
277
278 if (id == Boundary::slip) {
279 dealii::Tensor<1, dim, Number> V_i;
280 for (unsigned int d = 0; d < dim; ++d)
281 V_i[d] = dst.block(d).local_element(i);
282
283 /* replace normal component by source */
284 V_i -= 1. * (V_i * normal) * normal;
285 for (unsigned int d = 0; d < dim; ++d) {
286 const auto src_d = src.block(d).local_element(i);
287 V_i += 1. * (src_d * normal[d]) * normal;
288 }
289
290 for (unsigned int d = 0; d < dim; ++d)
291 dst.block(d).local_element(i) = V_i[d];
292
293 } else if (id == Boundary::no_slip || id == Boundary::dirichlet) {
294
295 /* set dst to src vector: */
296 for (unsigned int d = 0; d < dim; ++d)
297 dst.block(d).local_element(i) = src.block(d).local_element(i);
298 }
299 }
300 }
301
303 std::shared_ptr<DiagonalMatrix<dim, Number>> &matrix) const
304 {
305 Assert(level_ != dealii::numbers::invalid_unsigned_int,
306 dealii::ExcNotImplemented());
307 matrix = std::make_shared<DiagonalMatrix<dim, Number>>();
308 block_vector_type &vector = matrix->get_block_vector();
309 vector.reinit(dim);
310 for (unsigned int d = 0; d < dim; ++d)
311 matrix_free_->initialize_dof_vector(vector.block(d));
312 vector.collect_sizes();
313
314 const auto &lumped_mass_matrix =
315 offline_data_->level_lumped_mass_matrix()[level_];
316
317 unsigned int dummy = 0;
318 matrix_free_->template cell_loop<block_vector_type, unsigned int>(
319 [this](
320 const auto &data, auto &dst, const auto &, const auto range) {
321 dealii::FEEvaluation<dim, order_fe, order_quad, dim, Number>
322 velocity(data);
323 dealii::FEEvaluation<dim, order_fe, order_quad, dim, Number>
324 writer(data);
325
326 for (unsigned int cell = range.first; cell < range.second;
327 ++cell) {
328 velocity.reinit(cell);
329 writer.reinit(cell);
330 for (unsigned int i = 0; i < velocity.dofs_per_cell; ++i) {
331 for (unsigned int j = 0; j < velocity.dofs_per_cell; ++j)
332 velocity.begin_dof_values()[j] =
333 dealii::VectorizedArray<Number>();
334 velocity.begin_dof_values()[i] =
335 dealii::make_vectorized_array<Number>(1.);
336 apply_local_operator(velocity);
337 writer.begin_dof_values()[i] = velocity.begin_dof_values()[i];
338 }
339 writer.distribute_local_to_global(dst);
340 }
341 },
342 vector,
343 dummy,
344 /* zero destination */ true);
345
346 const unsigned int n_owned =
347 lumped_mass_matrix.get_partitioner()->locally_owned_size();
348
350
352 for (unsigned int i = 0; i < n_owned; ++i) {
353 const auto m_i = lumped_mass_matrix.local_element(i);
354 const auto rho_i = density_->local_element(i);
355 for (unsigned int d = 0; d < dim; ++d)
356 vector.block(d).local_element(i) =
357 1. / (m_i * rho_i + vector.block(d).local_element(i));
358 }
359
361
362 const auto &boundary_map = offline_data_->level_boundary_map()[level_];
363
364 for (auto entry : boundary_map) {
365 // [i, normal, normal_mass, boundary_mass, id, position] = entry
366 const auto i = std::get<0>(entry);
367 if (i >= n_owned)
368 continue;
369
370 const dealii::Tensor<1, dim, Number> normal = std::get<1>(entry);
371 const auto id = std::get<4>(entry);
372
373 if (id == Boundary::slip) {
374 dealii::Tensor<1, dim, Number> V_i;
375 for (unsigned int d = 0; d < dim; ++d)
376 V_i[d] = vector.block(d).local_element(i);
377
378 /* replace normal component by 1. */
379 V_i -= 1. * (V_i * normal) * normal;
380 for (unsigned int d = 0; d < dim; ++d) {
381 V_i += 1. * (1. * normal[d]) * normal;
382 }
383
384 for (unsigned int d = 0; d < dim; ++d)
385 vector.block(d).local_element(i) = V_i[d];
386
387 } else if (id == Boundary::no_slip || id == Boundary::dirichlet) {
388
389 /* set dst to src vector: */
390 for (unsigned int d = 0; d < dim; ++d)
391 vector.block(d).local_element(i) = 1.;
392 }
393 }
394 }
395
396 private:
397 const ParabolicSystem *parabolic_system_;
398 const OfflineData<dim, Number2> *offline_data_;
399 const dealii::MatrixFree<dim, Number> *matrix_free_;
400 const vector_type *density_;
401 Number theta_x_tau_;
402 unsigned int level_;
403
404 template <typename Evaluator>
405 void apply_local_operator(Evaluator &velocity) const
406 {
407 const Number mu = parabolic_system_->mu();
408 const Number lambda = parabolic_system_->lambda();
409
410 velocity.evaluate(dealii::EvaluationFlags::gradients);
411
412 for (unsigned int q = 0; q < velocity.n_q_points; ++q) {
413 if constexpr (dim == 1) {
414 /* Workaround: no symmetric gradient for dim == 1: */
415 const auto gradient = velocity.get_gradient(q);
416 auto S = (4. / 3. * mu + lambda) * gradient;
417 velocity.submit_gradient(theta_x_tau_ * S, q);
418
419 } else {
420
421 const auto symmetric_gradient = velocity.get_symmetric_gradient(q);
422 const auto divergence = trace(symmetric_gradient);
423 // S = (2 mu nabla^S(v) + (lambda - 2/3*mu) div(v) Id) : nabla phi
424 auto S = 2. * mu * symmetric_gradient;
425 for (unsigned int d = 0; d < dim; ++d)
426 S[d][d] += (lambda - 2. / 3. * mu) * divergence;
427 velocity.submit_symmetric_gradient(theta_x_tau_ * S, q);
428 }
429 }
430
431 velocity.integrate(dealii::EvaluationFlags::gradients);
432 }
433 };
434
435
439 template <int dim, typename Number>
441 : public dealii::MGTransferBase<
442 dealii::LinearAlgebra::distributed::BlockVector<Number>>
443 {
444 public:
445 using scalar_type = dealii::LinearAlgebra::distributed::Vector<Number>;
447 dealii::LinearAlgebra::distributed::BlockVector<Number>;
448
450
451 void build(const dealii::DoFHandler<dim> &dof_handler,
452 const dealii::MGConstrainedDoFs &mg_constrained_dofs,
453 const dealii::MGLevelObject<dealii::MatrixFree<dim, Number>>
454 &matrix_free)
455 {
456 transfer_.initialize_constraints(mg_constrained_dofs);
457 transfer_.build(dof_handler);
458 level_matrix_free_ = &matrix_free;
459 scalar_vector.resize(matrix_free.min_level(), matrix_free.max_level());
460 for (unsigned int level = matrix_free.min_level();
461 level < matrix_free.max_level();
462 ++level)
463 matrix_free[level].initialize_dof_vector(scalar_vector[level]);
464 }
465
466 void prolongate(const unsigned int to_level,
467 vector_type &dst,
468 const vector_type &src) const override
469 {
470 for (unsigned int block = 0; block < src.n_blocks(); ++block)
471 transfer_.prolongate(to_level, dst.block(block), src.block(block));
472 }
473
474 void restrict_and_add(const unsigned int to_level,
475 vector_type &dst,
476 const vector_type &src) const override
477 {
478 for (unsigned int block = 0; block < src.n_blocks(); ++block)
479 transfer_.restrict_and_add(
480 to_level, dst.block(block), src.block(block));
481 }
482
483 template <typename Number2>
485 const dealii::DoFHandler<dim> &dof_handler,
486 dealii::MGLevelObject<scalar_type> &dst,
487 const dealii::LinearAlgebra::distributed::Vector<Number2> &src) const
488 {
489 if (dst[dst.min_level()].size() == 0)
490 for (unsigned int l = dst.min_level(); l <= dst.max_level(); ++l)
491 (*level_matrix_free_)[l].initialize_dof_vector(dst[l]);
492 transfer_.interpolate_to_mg(dof_handler, dst, src);
493 }
494
495 template <typename Number2>
496 void
497 copy_to_mg(const dealii::DoFHandler<dim> &dof_handler,
498 dealii::MGLevelObject<vector_type> &dst,
499 const dealii::LinearAlgebra::distributed::BlockVector<Number2>
500 &src) const
501 {
502 if (dst[dst.min_level()].size() == 0)
503 for (unsigned int l = dst.min_level(); l <= dst.max_level(); ++l) {
504 dst[l].reinit(src.n_blocks());
505 for (unsigned int block = 0; block < src.n_blocks(); ++block)
506 (*level_matrix_free_)[l].initialize_dof_vector(
507 dst[l].block(block));
508 dst[l].collect_sizes();
509 }
510
511 for (unsigned int block = 0; block < src.n_blocks(); ++block) {
512 transfer_.copy_to_mg(dof_handler, scalar_vector, src.block(block));
513 for (unsigned int level = dst.min_level(); level <= dst.max_level();
514 ++level)
515 dst[level].block(block).copy_locally_owned_data_from(
516 scalar_vector[level]);
517 }
518 }
519
520 template <typename Number2>
522 const dealii::DoFHandler<dim> &dof_handler,
523 dealii::LinearAlgebra::distributed::BlockVector<Number2> &dst,
524 const dealii::MGLevelObject<vector_type> &src) const
525 {
526 for (unsigned int block = 0; block < dst.n_blocks(); ++block) {
527 for (unsigned int level = src.min_level(); level <= src.max_level();
528 ++level)
529 scalar_vector[level].copy_locally_owned_data_from(
530 src[level].block(block));
531 transfer_.copy_from_mg(dof_handler, dst.block(block), scalar_vector);
532 }
533 }
534
535 private:
536 dealii::MGTransferMatrixFree<dim, Number> transfer_;
537 const dealii::MGLevelObject<dealii::MatrixFree<dim, Number>>
538 *level_matrix_free_;
539 mutable dealii::MGLevelObject<scalar_type> scalar_vector;
540 };
541
542
546 template <int dim, typename Number, typename Number2>
547 class EnergyMatrix : public dealii::Subscriptor
548 {
549 public:
550 // FIXME: refactor
551 static constexpr unsigned int order_fe = 1;
552 static constexpr unsigned int order_quad = 2;
553
554 using vector_type = dealii::LinearAlgebra::distributed::Vector<Number>;
555
556 EnergyMatrix() = default;
557
559 const OfflineData<dim, Number2> &offline_data,
560 const dealii::MatrixFree<dim, Number> &matrix_free,
561 const dealii::LinearAlgebra::distributed::Vector<Number> &density,
562 const Number time_factor,
563 const unsigned int level = dealii::numbers::invalid_unsigned_int)
564 {
565 offline_data_ = &offline_data;
566 matrix_free_ = &matrix_free;
567 density_ = &density;
568 factor_ = time_factor;
569 level_ = level;
570 }
571
572 void Tvmult(vector_type &dst, const vector_type &src) const
573 {
574 vmult(dst, src);
575 }
576
577 dealii::types::global_dof_index m() const
578 {
579 return density_->size();
580 }
581
582 Number el(const unsigned int, const unsigned int) const
583 {
584 Assert(false, dealii::ExcNotImplemented());
585 return Number();
586 }
587
588 void vmult(vector_type &dst, const vector_type &src) const
589 {
590 /* Apply action of m_i rho_i V_i: */
591
592 using VA = dealii::VectorizedArray<Number>;
593 constexpr auto simd_length = VA::size();
594
595 const vector_type *lumped_mass_matrix = nullptr;
596 if constexpr (std::is_same<Number, Number2>::value) {
597 if constexpr (std::is_same<Number, float>::value) {
598 if (level_ == dealii::numbers::invalid_unsigned_int)
599 lumped_mass_matrix = &offline_data_->lumped_mass_matrix();
600 else
601 lumped_mass_matrix =
602 &offline_data_->level_lumped_mass_matrix()[level_];
603 } else {
604 Assert(level_ == dealii::numbers::invalid_unsigned_int,
605 dealii::ExcInternalError());
606 lumped_mass_matrix = &offline_data_->lumped_mass_matrix();
607 }
608 } else
609 lumped_mass_matrix =
610 &offline_data_->level_lumped_mass_matrix()[level_];
611
612 const unsigned int n_owned =
613 lumped_mass_matrix->get_partitioner()->locally_owned_size();
614 const unsigned int size_regular = n_owned / simd_length * simd_length;
615
617
619 for (unsigned int i = 0; i < size_regular; i += simd_length) {
620 const auto m_i = get_entry<VA>(*lumped_mass_matrix, i);
621 const auto rho_i = get_entry<VA>(*density_, i);
622 const auto e_i = get_entry<VA>(src, i);
623 write_entry<VA>(dst, m_i * rho_i * e_i, i);
624 }
625
627
628 for (unsigned int i = size_regular; i < n_owned; ++i) {
629 const auto m_i = lumped_mass_matrix->local_element(i);
630 const auto rho_i = density_->local_element(i);
631 const auto e_i = src.local_element(i);
632 dst.local_element(i) = m_i * rho_i * e_i;
633 }
634
635 /* Apply action of diffusion operator \sum_j beta_ij e_j: */
636
637 const auto integrator = [this](const auto &data,
638 auto &dst,
639 const auto &src,
640 const auto range) {
641 dealii::FEEvaluation<dim, order_fe, order_quad, 1, Number> energy(
642 data);
643
644 for (unsigned int cell = range.first; cell < range.second; ++cell) {
645 energy.reinit(cell);
646 energy.read_dof_values(src);
647 apply_local_operator(energy);
648 energy.distribute_local_to_global(dst);
649 }
650 };
651
652 matrix_free_->template cell_loop<vector_type, vector_type>(
653 integrator, dst, src, /* zero destination */ false);
654
655 /* Fix up constrained degrees of freedom: */
656
657 const auto &boundary_map =
658 (level_ == dealii::numbers::invalid_unsigned_int)
659 ? offline_data_->boundary_map()
660 : offline_data_->level_boundary_map()[level_];
661
662 for (auto entry : boundary_map) {
663 const auto i = std::get<0>(entry);
664 if (i >= n_owned)
665 continue;
666
667 const auto id = std::get<4>(entry);
668 if (id == Boundary::dirichlet)
669 dst.local_element(i) = src.local_element(i);
670 }
671 }
672
674 std::shared_ptr<dealii::DiagonalMatrix<vector_type>> &matrix) const
675 {
676 Assert(level_ != dealii::numbers::invalid_unsigned_int,
677 dealii::ExcNotImplemented());
678 matrix = std::make_shared<dealii::DiagonalMatrix<vector_type>>();
679 vector_type &vector = matrix->get_vector();
680 matrix_free_->initialize_dof_vector(vector);
681
682 const vector_type &lumped_mass_matrix =
683 offline_data_->level_lumped_mass_matrix()[level_];
684
685 unsigned int dummy = 0;
686 matrix_free_->template cell_loop<vector_type, unsigned int>(
687 [this](
688 const auto &data, auto &dst, const auto &, const auto range) {
689 dealii::FEEvaluation<dim, order_fe, order_quad, 1, Number> energy(
690 data);
691 dealii::FEEvaluation<dim, order_fe, order_quad, 1, Number> writer(
692 data);
693
694 for (unsigned int cell = range.first; cell < range.second;
695 ++cell) {
696 energy.reinit(cell);
697 writer.reinit(cell);
698 for (unsigned int i = 0; i < energy.dofs_per_cell; ++i) {
699 for (unsigned int j = 0; j < energy.dofs_per_cell; ++j)
700 energy.begin_dof_values()[j] =
701 dealii::VectorizedArray<Number>();
702 energy.begin_dof_values()[i] =
703 dealii::make_vectorized_array<Number>(1.);
704 apply_local_operator(energy);
705 writer.begin_dof_values()[i] = energy.begin_dof_values()[i];
706 }
707 writer.distribute_local_to_global(dst);
708 }
709 },
710 vector,
711 dummy,
712 /* zero destination */ true);
713
714 const unsigned int n_owned =
715 lumped_mass_matrix.get_partitioner()->locally_owned_size();
716
718
720 for (unsigned int i = 0; i < n_owned; ++i) {
721 const auto m_i = lumped_mass_matrix.local_element(i);
722 const auto rho_i = density_->local_element(i);
723 vector.local_element(i) =
724 1. / (m_i * rho_i + vector.local_element(i));
725 }
726
728
729 const auto &boundary_map = offline_data_->level_boundary_map()[level_];
730
731 for (auto entry : boundary_map) {
732 const auto i = std::get<0>(entry);
733 if (i >= n_owned)
734 continue;
735
736 const auto id = std::get<4>(entry);
737 if (id == Boundary::dirichlet)
738 vector.local_element(i) = 1.;
739 }
740 }
741
742 private:
743 const OfflineData<dim, Number2> *offline_data_;
744 const dealii::MatrixFree<dim, Number> *matrix_free_;
745 const dealii::LinearAlgebra::distributed::Vector<Number> *density_;
746 Number factor_;
747 unsigned int level_;
748
749 template <typename Evaluator>
750 void apply_local_operator(Evaluator &energy) const
751 {
752 energy.evaluate(dealii::EvaluationFlags::gradients);
753 for (unsigned int q = 0; q < energy.n_q_points; ++q) {
754 energy.submit_gradient(factor_ * energy.get_gradient(q), q);
755 }
756 energy.integrate(dealii::EvaluationFlags::gradients);
757 }
758 };
759
760
764 template <int dim, typename Number>
765 class MGTransferEnergy : public dealii::MGTransferMatrixFree<dim, Number>
766 {
767 public:
768 void build(const dealii::DoFHandler<dim> &dof_handler,
769 const dealii::MGLevelObject<dealii::MatrixFree<dim, Number>>
770 &matrix_free)
771 {
772 dealii::MGTransferMatrixFree<dim, Number>::build(dof_handler);
773 level_matrix_free_ = &matrix_free;
774 }
775
776 template <typename Number2>
778 const dealii::DoFHandler<dim> &dof_handler,
779 dealii::MGLevelObject<
780 dealii::LinearAlgebra::distributed::Vector<Number>> &dst,
781 const dealii::LinearAlgebra::distributed::Vector<Number2> &src) const
782 {
783 if (dst[dst.min_level()].size() == 0)
784 for (unsigned int l = dst.min_level(); l <= dst.max_level(); ++l)
785 (*level_matrix_free_)[l].initialize_dof_vector(dst[l]);
786 dealii::MGTransferMatrixFree<dim, Number>::copy_to_mg(
787 dof_handler, dst, src);
788 }
789
790 private:
791 const dealii::MGLevelObject<dealii::MatrixFree<dim, Number>>
792 *level_matrix_free_;
793 };
794
795 } // namespace NavierStokes
796} /* namespace ryujin */
797
798#undef locally_owned_size
void vmult(vector_type &dst, const vector_type &src) const
void vmult(block_vector_type &dst, const block_vector_type &src) const
dealii::LinearAlgebra::distributed::Vector< Number > vector_type
void reinit(const vector_type &lumped_mass_matrix, const vector_type &density, const dealii::AffineConstraints< Number > &affine_constraints)
dealii::LinearAlgebra::distributed::BlockVector< Number > block_vector_type
dealii::LinearAlgebra::distributed::Vector< Number > vector_type
void compute_diagonal(std::shared_ptr< dealii::DiagonalMatrix< vector_type > > &matrix) const
void Tvmult(vector_type &dst, const vector_type &src) const
void vmult(vector_type &dst, const vector_type &src) const
Number el(const unsigned int, const unsigned int) const
dealii::types::global_dof_index m() const
void initialize(const OfflineData< dim, Number2 > &offline_data, const dealii::MatrixFree< dim, Number > &matrix_free, const dealii::LinearAlgebra::distributed::Vector< Number > &density, const Number time_factor, const unsigned int level=dealii::numbers::invalid_unsigned_int)
void build(const dealii::DoFHandler< dim > &dof_handler, const dealii::MGLevelObject< dealii::MatrixFree< dim, Number > > &matrix_free)
void copy_to_mg(const dealii::DoFHandler< dim > &dof_handler, dealii::MGLevelObject< dealii::LinearAlgebra::distributed::Vector< Number > > &dst, const dealii::LinearAlgebra::distributed::Vector< Number2 > &src) const
void prolongate(const unsigned int to_level, vector_type &dst, const vector_type &src) const override
void copy_to_mg(const dealii::DoFHandler< dim > &dof_handler, dealii::MGLevelObject< vector_type > &dst, const dealii::LinearAlgebra::distributed::BlockVector< Number2 > &src) const
void restrict_and_add(const unsigned int to_level, vector_type &dst, const vector_type &src) const override
dealii::LinearAlgebra::distributed::Vector< Number > scalar_type
dealii::LinearAlgebra::distributed::BlockVector< Number > vector_type
void interpolate_to_mg(const dealii::DoFHandler< dim > &dof_handler, dealii::MGLevelObject< scalar_type > &dst, const dealii::LinearAlgebra::distributed::Vector< Number2 > &src) const
void build(const dealii::DoFHandler< dim > &dof_handler, const dealii::MGConstrainedDoFs &mg_constrained_dofs, const dealii::MGLevelObject< dealii::MatrixFree< dim, Number > > &matrix_free)
void copy_from_mg(const dealii::DoFHandler< dim > &dof_handler, dealii::LinearAlgebra::distributed::BlockVector< Number2 > &dst, const dealii::MGLevelObject< vector_type > &src) const
void vmult(block_vector_type &dst, const block_vector_type &src) const
dealii::LinearAlgebra::distributed::BlockVector< Number > block_vector_type
void initialize(const ParabolicSystem &parabolic_system, const OfflineData< dim, Number2 > &offline_data, const dealii::MatrixFree< dim, Number > &matrix_free, const dealii::LinearAlgebra::distributed::Vector< Number > &density, const Number theta_x_tau, const unsigned int level=dealii::numbers::invalid_unsigned_int)
dealii::LinearAlgebra::distributed::Vector< Number > vector_type
void Tvmult(block_vector_type &dst, const block_vector_type &src) const
void compute_diagonal(std::shared_ptr< DiagonalMatrix< dim, Number > > &matrix) const
auto & boundary_map() const
Definition: offline_data.h:194
auto & level_lumped_mass_matrix() const
Definition: offline_data.h:250
auto & lumped_mass_matrix() const
Definition: offline_data.h:238
auto & level_boundary_map() const
Definition: offline_data.h:208
#define RYUJIN_PARALLEL_REGION_BEGIN
Definition: openmp.h:54
#define RYUJIN_OMP_FOR
Definition: openmp.h:70
#define RYUJIN_PARALLEL_REGION_END
Definition: openmp.h:63
Euler::HyperbolicSystem HyperbolicSystem