ryujin 2.1.1 revision 2062b43aae59b7ee510dec4ae666e29d56c4f0be
parabolic_solver_gmg_operators.h
Go to the documentation of this file.
1//
2// SPDX-License-Identifier: MIT
3// Copyright (C) 2020 - 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 using vector_type = dealii::LinearAlgebra::distributed::Vector<Number>;
161 dealii::LinearAlgebra::distributed::BlockVector<Number>;
162
163 VelocityMatrix() = default;
164
166 const ParabolicSystem &parabolic_system,
167 const OfflineData<dim, Number2> &offline_data,
168 const dealii::MatrixFree<dim, Number> &matrix_free,
169 const dealii::LinearAlgebra::distributed::Vector<Number> &density,
170 const Number theta_x_tau,
171 const unsigned int level = dealii::numbers::invalid_unsigned_int)
172 {
173 parabolic_system_ = &parabolic_system;
174 offline_data_ = &offline_data;
175 matrix_free_ = &matrix_free;
176 density_ = &density;
177 theta_x_tau_ = theta_x_tau;
178 level_ = level;
179 }
180
181 void Tvmult(block_vector_type &dst, const block_vector_type &src) const
182 {
183 vmult(dst, src);
184 }
185
186 void vmult(block_vector_type &dst, const block_vector_type &src) const
187 {
188 /* Apply action of m_i rho_i V_i: */
189
190 using VA = dealii::VectorizedArray<Number>;
191 constexpr auto simd_length = VA::size();
192
193 const vector_type *lumped_mass_matrix = nullptr;
194 if constexpr (std::is_same<Number, Number2>::value) {
195 if constexpr (std::is_same<Number, float>::value) {
196 if (level_ == dealii::numbers::invalid_unsigned_int)
197 lumped_mass_matrix = &offline_data_->lumped_mass_matrix();
198 else
199 lumped_mass_matrix =
200 &offline_data_->level_lumped_mass_matrix()[level_];
201 } else {
202 Assert(level_ == dealii::numbers::invalid_unsigned_int,
203 dealii::ExcInternalError());
204 lumped_mass_matrix = &offline_data_->lumped_mass_matrix();
205 }
206 } else
207 lumped_mass_matrix =
208 &offline_data_->level_lumped_mass_matrix()[level_];
209
210 const unsigned int n_owned =
211 lumped_mass_matrix->get_partitioner()->locally_owned_size();
212 const unsigned int size_regular = n_owned / simd_length * simd_length;
213
215
217 for (unsigned int i = 0; i < size_regular; i += simd_length) {
218 const auto m_i = load_value<VA>(*lumped_mass_matrix, i);
219 const auto rho_i = load_value<VA>(*density_, i);
220 for (unsigned int d = 0; d < dim; ++d) {
221 const auto temp = load_value<VA>(src.block(d), i);
222 store_value<VA>(dst.block(d), m_i * rho_i * temp, i);
223 }
224 }
225
227
228 for (unsigned int i = size_regular; i < n_owned; ++i) {
229 const auto m_i = lumped_mass_matrix->local_element(i);
230 const auto rho_i = density_->local_element(i);
231
232 for (unsigned int d = 0; d < dim; ++d) {
233 const auto temp = src.block(d).local_element(i);
234 dst.block(d).local_element(i) = m_i * rho_i * temp;
235 }
236 }
237
238 /* Apply action of stress tensor: + theta * \sum_j B_ij V_j: */
239
240 const auto integrator = [this](const auto &data,
241 auto &dst,
242 const auto &src,
243 const auto range) {
244 constexpr auto order_fe = Discretization<dim>::order_finite_element;
245 constexpr auto order_quad = Discretization<dim>::order_quadrature;
246 dealii::FEEvaluation<dim, order_fe, order_quad, dim, Number> velocity(
247 data);
248
249 for (unsigned int cell = range.first; cell < range.second; ++cell) {
250 velocity.reinit(cell);
251 velocity.read_dof_values(src);
252 apply_local_operator(velocity);
253 velocity.distribute_local_to_global(dst);
254 }
255 };
256
257 matrix_free_->template cell_loop<block_vector_type, block_vector_type>(
258 integrator, dst, src, /* zero destination */ false);
259
260 /* (5.4a) Fix up constrained degrees of freedom: */
261
262 const auto &boundary_map =
263 level_ == dealii::numbers::invalid_unsigned_int
264 ? offline_data_->boundary_map()
265 : offline_data_->level_boundary_map()[level_];
266
267 for (auto entry : boundary_map) {
268 const auto i = entry.first;
269 if (i >= n_owned)
270 continue;
271
272 const dealii::Tensor<1, dim, Number> normal =
273 std::get<0>(entry.second);
274 const auto id = std::get<3>(entry.second);
275
276 if (id == Boundary::slip) {
277 dealii::Tensor<1, dim, Number> V_i;
278 for (unsigned int d = 0; d < dim; ++d)
279 V_i[d] = dst.block(d).local_element(i);
280
281 /* replace normal component by source */
282 V_i -= 1. * (V_i * normal) * normal;
283 for (unsigned int d = 0; d < dim; ++d) {
284 const auto src_d = src.block(d).local_element(i);
285 V_i += 1. * (src_d * normal[d]) * normal;
286 }
287
288 for (unsigned int d = 0; d < dim; ++d)
289 dst.block(d).local_element(i) = V_i[d];
290
291 } else if (id == Boundary::no_slip || id == Boundary::dirichlet) {
292
293 /* set dst to src vector: */
294 for (unsigned int d = 0; d < dim; ++d)
295 dst.block(d).local_element(i) = src.block(d).local_element(i);
296 }
297 }
298 }
299
301 std::shared_ptr<DiagonalMatrix<dim, Number>> &matrix) const
302 {
303 Assert(level_ != dealii::numbers::invalid_unsigned_int,
304 dealii::ExcNotImplemented());
305 matrix = std::make_shared<DiagonalMatrix<dim, Number>>();
306 block_vector_type &vector = matrix->get_block_vector();
307 vector.reinit(dim);
308 for (unsigned int d = 0; d < dim; ++d)
309 matrix_free_->initialize_dof_vector(vector.block(d));
310 vector.collect_sizes();
311
312 const auto &lumped_mass_matrix =
313 offline_data_->level_lumped_mass_matrix()[level_];
314
315 unsigned int dummy = 0;
316 matrix_free_->template cell_loop<block_vector_type, unsigned int>(
317 [this](
318 const auto &data, auto &dst, const auto &, const auto range) {
319 constexpr auto order_fe =
321 constexpr auto order_quad = Discretization<dim>::order_quadrature;
322 dealii::FEEvaluation<dim, order_fe, order_quad, dim, Number>
323 velocity(data);
324 dealii::FEEvaluation<dim, order_fe, order_quad, dim, Number>
325 writer(data);
326
327 for (unsigned int cell = range.first; cell < range.second;
328 ++cell) {
329 velocity.reinit(cell);
330 writer.reinit(cell);
331 for (unsigned int i = 0; i < velocity.dofs_per_cell; ++i) {
332 for (unsigned int j = 0; j < velocity.dofs_per_cell; ++j)
333 velocity.begin_dof_values()[j] =
334 dealii::VectorizedArray<Number>();
335 velocity.begin_dof_values()[i] =
336 dealii::make_vectorized_array<Number>(1.);
337 apply_local_operator(velocity);
338 writer.begin_dof_values()[i] = velocity.begin_dof_values()[i];
339 }
340 writer.distribute_local_to_global(dst);
341 }
342 },
343 vector,
344 dummy,
345 /* zero destination */ true);
346
347 const unsigned int n_owned =
348 lumped_mass_matrix.get_partitioner()->locally_owned_size();
349
351
353 for (unsigned int i = 0; i < n_owned; ++i) {
354 const auto m_i = lumped_mass_matrix.local_element(i);
355 const auto rho_i = density_->local_element(i);
356 for (unsigned int d = 0; d < dim; ++d)
357 vector.block(d).local_element(i) =
358 1. / (m_i * rho_i + vector.block(d).local_element(i));
359 }
360
362
363 const auto &boundary_map = offline_data_->level_boundary_map()[level_];
364
365 for (auto entry : boundary_map) {
366 const auto i = entry.first;
367 if (i >= n_owned)
368 continue;
369
370 const dealii::Tensor<1, dim, Number> normal =
371 std::get<0>(entry.second);
372 const auto id = std::get<3>(entry.second);
373
374 if (id == Boundary::slip) {
375 dealii::Tensor<1, dim, Number> V_i;
376 for (unsigned int d = 0; d < dim; ++d)
377 V_i[d] = vector.block(d).local_element(i);
378
379 /* replace normal component by 1. */
380 V_i -= 1. * (V_i * normal) * normal;
381 for (unsigned int d = 0; d < dim; ++d) {
382 V_i += 1. * (1. * normal[d]) * normal;
383 }
384
385 for (unsigned int d = 0; d < dim; ++d)
386 vector.block(d).local_element(i) = V_i[d];
387
388 } else if (id == Boundary::no_slip || id == Boundary::dirichlet) {
389
390 /* set dst to src vector: */
391 for (unsigned int d = 0; d < dim; ++d)
392 vector.block(d).local_element(i) = 1.;
393 }
394 }
395 }
396
397 private:
398 const ParabolicSystem *parabolic_system_;
399 const OfflineData<dim, Number2> *offline_data_;
400 const dealii::MatrixFree<dim, Number> *matrix_free_;
401 const vector_type *density_;
402 Number theta_x_tau_;
403 unsigned int level_;
404
405 template <typename Evaluator>
406 void apply_local_operator(Evaluator &velocity) const
407 {
408 const Number mu = parabolic_system_->mu();
409 const Number lambda = parabolic_system_->lambda();
410
411 velocity.evaluate(dealii::EvaluationFlags::gradients);
412
413 for (unsigned int q = 0; q < velocity.n_q_points; ++q) {
414 if constexpr (dim == 1) {
415 /* Workaround: no symmetric gradient for dim == 1: */
416 const auto gradient = velocity.get_gradient(q);
417 auto S = (4. / 3. * mu + lambda) * gradient;
418 velocity.submit_gradient(theta_x_tau_ * S, q);
419
420 } else {
421
422 const auto symmetric_gradient = velocity.get_symmetric_gradient(q);
423 const auto divergence = trace(symmetric_gradient);
424 // S = (2 mu nabla^S(v) + (lambda - 2/3*mu) div(v) Id) : nabla phi
425 auto S = 2. * mu * symmetric_gradient;
426 for (unsigned int d = 0; d < dim; ++d)
427 S[d][d] += (lambda - 2. / 3. * mu) * divergence;
428 velocity.submit_symmetric_gradient(theta_x_tau_ * S, q);
429 }
430 }
431
432 velocity.integrate(dealii::EvaluationFlags::gradients);
433 }
434 };
435
436
440 template <int dim, typename Number>
442 : public dealii::MGTransferBase<
443 dealii::LinearAlgebra::distributed::BlockVector<Number>>
444 {
445 public:
446 using scalar_type = dealii::LinearAlgebra::distributed::Vector<Number>;
448 dealii::LinearAlgebra::distributed::BlockVector<Number>;
449
451
452 void build(const dealii::DoFHandler<dim> &dof_handler,
453 const dealii::MGConstrainedDoFs &mg_constrained_dofs,
454 const dealii::MGLevelObject<dealii::MatrixFree<dim, Number>>
455 &matrix_free)
456 {
457 transfer_.initialize_constraints(mg_constrained_dofs);
458 transfer_.build(dof_handler);
459 level_matrix_free_ = &matrix_free;
460 scalar_vector.resize(matrix_free.min_level(), matrix_free.max_level());
461 for (unsigned int level = matrix_free.min_level();
462 level < matrix_free.max_level();
463 ++level)
464 matrix_free[level].initialize_dof_vector(scalar_vector[level]);
465 }
466
467 void prolongate(const unsigned int to_level,
468 vector_type &dst,
469 const vector_type &src) const override
470 {
471 for (unsigned int block = 0; block < src.n_blocks(); ++block)
472 transfer_.prolongate(to_level, dst.block(block), src.block(block));
473 }
474
475 void restrict_and_add(const unsigned int to_level,
476 vector_type &dst,
477 const vector_type &src) const override
478 {
479 for (unsigned int block = 0; block < src.n_blocks(); ++block)
480 transfer_.restrict_and_add(
481 to_level, dst.block(block), src.block(block));
482 }
483
484 template <typename Number2>
486 const dealii::DoFHandler<dim> &dof_handler,
487 dealii::MGLevelObject<scalar_type> &dst,
488 const dealii::LinearAlgebra::distributed::Vector<Number2> &src) const
489 {
490 if (dst[dst.min_level()].size() == 0)
491 for (unsigned int l = dst.min_level(); l <= dst.max_level(); ++l)
492 (*level_matrix_free_)[l].initialize_dof_vector(dst[l]);
493 transfer_.interpolate_to_mg(dof_handler, dst, src);
494 }
495
496 template <typename Number2>
497 void
498 copy_to_mg(const dealii::DoFHandler<dim> &dof_handler,
499 dealii::MGLevelObject<vector_type> &dst,
500 const dealii::LinearAlgebra::distributed::BlockVector<Number2>
501 &src) const
502 {
503 if (dst[dst.min_level()].size() == 0)
504 for (unsigned int l = dst.min_level(); l <= dst.max_level(); ++l) {
505 dst[l].reinit(src.n_blocks());
506 for (unsigned int block = 0; block < src.n_blocks(); ++block)
507 (*level_matrix_free_)[l].initialize_dof_vector(
508 dst[l].block(block));
509 dst[l].collect_sizes();
510 }
511
512 for (unsigned int block = 0; block < src.n_blocks(); ++block) {
513 transfer_.copy_to_mg(dof_handler, scalar_vector, src.block(block));
514 for (unsigned int level = dst.min_level(); level <= dst.max_level();
515 ++level)
516 dst[level].block(block).copy_locally_owned_data_from(
517 scalar_vector[level]);
518 }
519 }
520
521 template <typename Number2>
523 const dealii::DoFHandler<dim> &dof_handler,
524 dealii::LinearAlgebra::distributed::BlockVector<Number2> &dst,
525 const dealii::MGLevelObject<vector_type> &src) const
526 {
527 for (unsigned int block = 0; block < dst.n_blocks(); ++block) {
528 for (unsigned int level = src.min_level(); level <= src.max_level();
529 ++level)
530 scalar_vector[level].copy_locally_owned_data_from(
531 src[level].block(block));
532 transfer_.copy_from_mg(dof_handler, dst.block(block), scalar_vector);
533 }
534 }
535
536 private:
537 dealii::MGTransferMatrixFree<dim, Number> transfer_;
538 const dealii::MGLevelObject<dealii::MatrixFree<dim, Number>>
539 *level_matrix_free_;
540 mutable dealii::MGLevelObject<scalar_type> scalar_vector;
541 };
542
543
547 template <int dim, typename Number, typename Number2>
548 class EnergyMatrix : public dealii::Subscriptor
549 {
550 public:
551 using vector_type = dealii::LinearAlgebra::distributed::Vector<Number>;
552
553 EnergyMatrix() = default;
554
556 const OfflineData<dim, Number2> &offline_data,
557 const dealii::MatrixFree<dim, Number> &matrix_free,
558 const dealii::LinearAlgebra::distributed::Vector<Number> &density,
559 const Number time_factor,
560 const unsigned int level = dealii::numbers::invalid_unsigned_int)
561 {
562 offline_data_ = &offline_data;
563 matrix_free_ = &matrix_free;
564 density_ = &density;
565 factor_ = time_factor;
566 level_ = level;
567 }
568
569 void Tvmult(vector_type &dst, const vector_type &src) const
570 {
571 vmult(dst, src);
572 }
573
574 dealii::types::global_dof_index m() const
575 {
576 return density_->size();
577 }
578
579 Number el(const unsigned int, const unsigned int) const
580 {
581 Assert(false, dealii::ExcNotImplemented());
582 return Number();
583 }
584
585 void vmult(vector_type &dst, const vector_type &src) const
586 {
587 /* Apply action of m_i rho_i V_i: */
588
589 using VA = dealii::VectorizedArray<Number>;
590 constexpr auto simd_length = VA::size();
591
592 const vector_type *lumped_mass_matrix = nullptr;
593 if constexpr (std::is_same<Number, Number2>::value) {
594 if constexpr (std::is_same<Number, float>::value) {
595 if (level_ == dealii::numbers::invalid_unsigned_int)
596 lumped_mass_matrix = &offline_data_->lumped_mass_matrix();
597 else
598 lumped_mass_matrix =
599 &offline_data_->level_lumped_mass_matrix()[level_];
600 } else {
601 Assert(level_ == dealii::numbers::invalid_unsigned_int,
602 dealii::ExcInternalError());
603 lumped_mass_matrix = &offline_data_->lumped_mass_matrix();
604 }
605 } else
606 lumped_mass_matrix =
607 &offline_data_->level_lumped_mass_matrix()[level_];
608
609 const unsigned int n_owned =
610 lumped_mass_matrix->get_partitioner()->locally_owned_size();
611 const unsigned int size_regular = n_owned / simd_length * simd_length;
612
614
616 for (unsigned int i = 0; i < size_regular; i += simd_length) {
617 const auto m_i = load_value<VA>(*lumped_mass_matrix, i);
618 const auto rho_i = load_value<VA>(*density_, i);
619 const auto e_i = load_value<VA>(src, i);
620 store_value<VA>(dst, m_i * rho_i * e_i, i);
621 }
622
624
625 for (unsigned int i = size_regular; i < n_owned; ++i) {
626 const auto m_i = lumped_mass_matrix->local_element(i);
627 const auto rho_i = density_->local_element(i);
628 const auto e_i = src.local_element(i);
629 dst.local_element(i) = m_i * rho_i * e_i;
630 }
631
632 /* Apply action of diffusion operator \sum_j beta_ij e_j: */
633
634 const auto integrator = [this](const auto &data,
635 auto &dst,
636 const auto &src,
637 const auto range) {
638 constexpr auto order_fe = Discretization<dim>::order_finite_element;
639 constexpr auto order_quad = Discretization<dim>::order_quadrature;
640 dealii::FEEvaluation<dim, order_fe, order_quad, 1, Number> energy(
641 data);
642
643 for (unsigned int cell = range.first; cell < range.second; ++cell) {
644 energy.reinit(cell);
645 energy.read_dof_values(src);
646 apply_local_operator(energy);
647 energy.distribute_local_to_global(dst);
648 }
649 };
650
651 matrix_free_->template cell_loop<vector_type, vector_type>(
652 integrator, dst, src, /* zero destination */ false);
653
654 /* Fix up constrained degrees of freedom: */
655
656 const auto &boundary_map =
657 (level_ == dealii::numbers::invalid_unsigned_int)
658 ? offline_data_->boundary_map()
659 : offline_data_->level_boundary_map()[level_];
660
661 for (auto entry : boundary_map) {
662 const auto i = entry.first;
663 if (i >= n_owned)
664 continue;
665
666 const auto id = std::get<3>(entry.second);
667 if (id == Boundary::dirichlet)
668 dst.local_element(i) = src.local_element(i);
669 }
670 }
671
673 std::shared_ptr<dealii::DiagonalMatrix<vector_type>> &matrix) const
674 {
675 Assert(level_ != dealii::numbers::invalid_unsigned_int,
676 dealii::ExcNotImplemented());
677 matrix = std::make_shared<dealii::DiagonalMatrix<vector_type>>();
678 vector_type &vector = matrix->get_vector();
679 matrix_free_->initialize_dof_vector(vector);
680
681 const vector_type &lumped_mass_matrix =
682 offline_data_->level_lumped_mass_matrix()[level_];
683
684 unsigned int dummy = 0;
685 matrix_free_->template cell_loop<vector_type, unsigned int>(
686 [this](
687 const auto &data, auto &dst, const auto &, const auto range) {
688 constexpr auto order_fe =
690 constexpr auto order_quad = Discretization<dim>::order_quadrature;
691 dealii::FEEvaluation<dim, order_fe, order_quad, 1, Number> energy(
692 data);
693 dealii::FEEvaluation<dim, order_fe, order_quad, 1, Number> writer(
694 data);
695
696 for (unsigned int cell = range.first; cell < range.second;
697 ++cell) {
698 energy.reinit(cell);
699 writer.reinit(cell);
700 for (unsigned int i = 0; i < energy.dofs_per_cell; ++i) {
701 for (unsigned int j = 0; j < energy.dofs_per_cell; ++j)
702 energy.begin_dof_values()[j] =
703 dealii::VectorizedArray<Number>();
704 energy.begin_dof_values()[i] =
705 dealii::make_vectorized_array<Number>(1.);
706 apply_local_operator(energy);
707 writer.begin_dof_values()[i] = energy.begin_dof_values()[i];
708 }
709 writer.distribute_local_to_global(dst);
710 }
711 },
712 vector,
713 dummy,
714 /* zero destination */ true);
715
716 const unsigned int n_owned =
717 lumped_mass_matrix.get_partitioner()->locally_owned_size();
718
720
722 for (unsigned int i = 0; i < n_owned; ++i) {
723 const auto m_i = lumped_mass_matrix.local_element(i);
724 const auto rho_i = density_->local_element(i);
725 vector.local_element(i) =
726 1. / (m_i * rho_i + vector.local_element(i));
727 }
728
730
731 const auto &boundary_map = offline_data_->level_boundary_map()[level_];
732
733 for (auto entry : boundary_map) {
734 const auto i = entry.first;
735 if (i >= n_owned)
736 continue;
737
738 const auto id = std::get<3>(entry.second);
739 if (id == Boundary::dirichlet)
740 vector.local_element(i) = 1.;
741 }
742 }
743
744 private:
745 const OfflineData<dim, Number2> *offline_data_;
746 const dealii::MatrixFree<dim, Number> *matrix_free_;
747 const dealii::LinearAlgebra::distributed::Vector<Number> *density_;
748 Number factor_;
749 unsigned int level_;
750
751 template <typename Evaluator>
752 void apply_local_operator(Evaluator &energy) const
753 {
754 energy.evaluate(dealii::EvaluationFlags::gradients);
755 for (unsigned int q = 0; q < energy.n_q_points; ++q) {
756 energy.submit_gradient(factor_ * energy.get_gradient(q), q);
757 }
758 energy.integrate(dealii::EvaluationFlags::gradients);
759 }
760 };
761
762
766 template <int dim, typename Number>
767 class MGTransferEnergy : public dealii::MGTransferMatrixFree<dim, Number>
768 {
769 public:
770 void build(const dealii::DoFHandler<dim> &dof_handler,
771 const dealii::MGLevelObject<dealii::MatrixFree<dim, Number>>
772 &matrix_free)
773 {
774 dealii::MGTransferMatrixFree<dim, Number>::build(dof_handler);
775 level_matrix_free_ = &matrix_free;
776 }
777
778 template <typename Number2>
780 const dealii::DoFHandler<dim> &dof_handler,
781 dealii::MGLevelObject<
782 dealii::LinearAlgebra::distributed::Vector<Number>> &dst,
783 const dealii::LinearAlgebra::distributed::Vector<Number2> &src) const
784 {
785 if (dst[dst.min_level()].size() == 0)
786 for (unsigned int l = dst.min_level(); l <= dst.max_level(); ++l)
787 (*level_matrix_free_)[l].initialize_dof_vector(dst[l]);
788 dealii::MGTransferMatrixFree<dim, Number>::copy_to_mg(
789 dof_handler, dst, src);
790 }
791
792 private:
793 const dealii::MGLevelObject<dealii::MatrixFree<dim, Number>>
794 *level_matrix_free_;
795 };
796
797 } // namespace NavierStokes
798} /* namespace ryujin */
799
800#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:161
auto & level_lumped_mass_matrix() const
Definition: offline_data.h:208
auto & lumped_mass_matrix() const
Definition: offline_data.h:197
auto & level_boundary_map() const
Definition: offline_data.h:175
#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