8#include <compile_time_options.h>
15#include "../euler/hyperbolic_system.h"
18#include <deal.II/base/vectorization.h>
19#include <deal.II/lac/diagonal_matrix.h>
20#include <deal.II/matrix_free/fe_evaluation.h>
21#include <deal.II/matrix_free/tools.h>
22#include <deal.II/multigrid/mg_base.h>
23#include <deal.II/multigrid/mg_transfer_matrix_free.h>
32 namespace NavierStokes
44 template <
int dim,
typename Number>
51 using vector_type = dealii::LinearAlgebra::distributed::Vector<Number>;
57 dealii::LinearAlgebra::distributed::BlockVector<Number>;
68 template <
typename Vector>
69 void reinit(
const Vector &lumped_mass_matrix,
71 const dealii::AffineConstraints<Number> &affine_constraints)
73 diagonal.reinit(density,
true);
75 const auto n_owned = density.get_partitioner()->locally_owned_size();
77 const auto body_invert = [&](
auto sentinel,
const unsigned int i) {
78 using T =
decltype(sentinel);
81 if constexpr (std::is_same_v<Vector, vector_type>) {
82 m_i = read_entry<T>(lumped_mass_matrix, i);
84 m_i = lumped_mass_matrix.template read_entry<T>(i);
87 const auto rho_i = read_entry<T>(density, i);
88 write_entry<T>(diagonal, Number(1.0) / (rho_i * m_i), i);
91 cpu_simd_loop<Number>(
"", body_invert, 0, n_owned, n_owned);
97 affine_constraints.set_zero(diagonal);
113 return diagonal_block;
121 AssertDimension(diagonal_block.size(), 0);
122 DEAL_II_OPENMP_SIMD_PRAGMA
123 for (
unsigned int i = 0;
124 i < diagonal.get_partitioner()->locally_owned_size();
126 dst.local_element(i) =
127 diagonal.local_element(i) * src.local_element(i);
135 AssertDimension(dim, dst.n_blocks());
136 AssertDimension(dim, src.n_blocks());
137 if (diagonal_block.size() == 0) {
138 DEAL_II_OPENMP_SIMD_PRAGMA
139 for (
unsigned int i = 0;
140 i < diagonal.get_partitioner()->locally_owned_size();
142 for (
unsigned int d = 0; d < dim; ++d)
143 dst.block(d).local_element(i) =
144 diagonal.local_element(i) * src.block(d).local_element(i);
146 for (
unsigned int d = 0; d < dim; ++d) {
147 DEAL_II_OPENMP_SIMD_PRAGMA
148 for (
unsigned int i = 0;
149 i < src.block(d).get_partitioner()->locally_owned_size();
151 dst.block(d).local_element(i) =
152 diagonal_block.block(d).local_element(i) *
153 src.block(d).local_element(i);
169 template <
int dim,
typename Number,
typename Number2>
177 using vector_type = dealii::LinearAlgebra::distributed::Vector<Number>;
179 dealii::LinearAlgebra::distributed::BlockVector<Number>;
186 const dealii::MatrixFree<dim, Number> &matrix_free,
187 const dealii::LinearAlgebra::distributed::Vector<Number> &density,
188 const Number theta_x_tau,
189 const unsigned int level = dealii::numbers::invalid_unsigned_int)
191 parabolic_system_ = ¶bolic_system;
192 offline_data_ = &offline_data;
193 matrix_free_ = &matrix_free;
195 theta_x_tau_ = theta_x_tau;
209 const auto get_lumped_mass = [&](
auto sentinel,
unsigned int i) {
210 using T =
decltype(sentinel);
211 if constexpr (std::is_same_v<Number, Number2>) {
212 if constexpr (std::is_same_v<Number, float>) {
213 if (level_ == dealii::numbers::invalid_unsigned_int) {
216 return lumped.template read_entry<T>(i);
218 const auto &level_lumped =
220 return read_entry<T>(level_lumped, i);
223 Assert(level_ == dealii::numbers::invalid_unsigned_int,
224 dealii::ExcInternalError());
227 return lumped.template read_entry<T>(i);
230 const auto &level_lumped =
232 return read_entry<T>(level_lumped, i);
236 const unsigned int n_owned =
237 dst.block(0).get_partitioner()->locally_owned_size();
239 const auto body_mass = [&](
auto sentinel,
unsigned int i) {
240 using T =
decltype(sentinel);
242 const auto m_i = get_lumped_mass(T(), i);
244 const auto rho_i = read_entry<T>(*density_, i);
245 for (
unsigned int d = 0; d < dim; ++d) {
246 const auto temp = read_entry<T>(src.block(d), i);
247 write_entry<T>(dst.block(d), m_i * rho_i * temp, i);
251 cpu_simd_loop<Number>(
"", body_mass, 0, n_owned, n_owned);
255 const auto integrator = [
this](
const auto &data,
259 dealii::FEEvaluation<dim, order_fe, order_quad, dim, Number> velocity(
262 for (
unsigned int cell = range.first; cell < range.second; ++cell) {
263 velocity.reinit(cell);
264 velocity.read_dof_values(src);
265 apply_local_operator(velocity);
266 velocity.distribute_local_to_global(dst);
270 matrix_free_->template cell_loop<block_vector_type, block_vector_type>(
271 integrator, dst, src,
false);
275 const auto &boundary_map =
276 level_ == dealii::numbers::invalid_unsigned_int
280 for (
auto entry : boundary_map) {
282 const auto i = std::get<0>(entry);
286 const dealii::Tensor<1, dim, Number> normal = std::get<1>(entry);
287 const auto id = std::get<4>(entry);
290 dealii::Tensor<1, dim, Number> V_i;
291 for (
unsigned int d = 0; d < dim; ++d)
292 V_i[d] = dst.block(d).local_element(i);
295 V_i -= 1. * (V_i * normal) * normal;
296 for (
unsigned int d = 0; d < dim; ++d) {
297 const auto src_d = src.block(d).local_element(i);
298 V_i += 1. * (src_d * normal[d]) * normal;
301 for (
unsigned int d = 0; d < dim; ++d)
302 dst.block(d).local_element(i) = V_i[d];
307 for (
unsigned int d = 0; d < dim; ++d)
308 dst.block(d).local_element(i) = src.block(d).local_element(i);
316 Assert(level_ != dealii::numbers::invalid_unsigned_int,
317 dealii::ExcNotImplemented());
318 matrix = std::make_shared<DiagonalMatrix<dim, Number>>();
321 for (
unsigned int d = 0; d < dim; ++d)
322 matrix_free_->initialize_dof_vector(vector.block(d));
323 vector.collect_sizes();
325 dealii::MatrixFreeTools::compute_diagonal(
328 &VelocityMatrix::template apply_local_operator<
329 dealii::FEEvaluation<dim, -1, 0, dim, Number>>,
332 const auto &lumped_mass_matrix =
334 const unsigned int n_owned =
335 lumped_mass_matrix.get_partitioner()->locally_owned_size();
337 const auto body_invert = [&](
auto sentinel,
const unsigned int i) {
338 using T =
decltype(sentinel);
339 const auto m_i = lumped_mass_matrix.local_element(i);
340 const auto rho_i = density_->local_element(i);
341 for (
unsigned int d = 0; d < dim; ++d)
344 (m_i * rho_i + read_entry<T>(vector.block(d), i)),
348 cpu_simd_loop<Number>(
"", body_invert, 0, n_owned, n_owned);
352 for (
auto entry : boundary_map) {
354 const auto i = std::get<0>(entry);
358 const dealii::Tensor<1, dim, Number> normal = std::get<1>(entry);
359 const auto id = std::get<4>(entry);
362 dealii::Tensor<1, dim, Number> V_i;
363 for (
unsigned int d = 0; d < dim; ++d)
364 V_i[d] = vector.block(d).local_element(i);
367 V_i -= 1. * (V_i * normal) * normal;
368 for (
unsigned int d = 0; d < dim; ++d) {
369 V_i += 1. * (1. * normal[d]) * normal;
372 for (
unsigned int d = 0; d < dim; ++d)
373 vector.block(d).local_element(i) = V_i[d];
378 for (
unsigned int d = 0; d < dim; ++d)
379 vector.block(d).local_element(i) = 1.;
387 const dealii::MatrixFree<dim, Number> *matrix_free_;
392 template <
typename Evaluator>
393 void apply_local_operator(Evaluator &velocity)
const
395 const Number mu = parabolic_system_->mu();
396 const Number lambda = parabolic_system_->lambda();
398 velocity.evaluate(dealii::EvaluationFlags::gradients);
400 for (
const unsigned int q : velocity.quadrature_point_indices()) {
401 if constexpr (dim == 1) {
403 const auto gradient = velocity.get_gradient(q);
404 auto S = (4. / 3. * mu + lambda) * gradient;
405 velocity.submit_gradient(theta_x_tau_ * S, q);
409 const auto symmetric_gradient = velocity.get_symmetric_gradient(q);
410 const auto divergence = trace(symmetric_gradient);
412 auto S = 2. * mu * symmetric_gradient;
413 for (
unsigned int d = 0; d < dim; ++d)
414 S[d][d] += (lambda - 2. / 3. * mu) * divergence;
415 velocity.submit_symmetric_gradient(theta_x_tau_ * S, q);
419 velocity.integrate(dealii::EvaluationFlags::gradients);
427 template <
int dim,
typename Number>
429 :
public dealii::MGTransferBase<
430 dealii::LinearAlgebra::distributed::BlockVector<Number>>
433 using scalar_type = dealii::LinearAlgebra::distributed::Vector<Number>;
435 dealii::LinearAlgebra::distributed::BlockVector<Number>;
439 void build(
const dealii::DoFHandler<dim> &dof_handler,
440 const dealii::MGConstrainedDoFs &mg_constrained_dofs,
441 const dealii::MGLevelObject<dealii::MatrixFree<dim, Number>>
444 transfer_.initialize_constraints(mg_constrained_dofs);
445 transfer_.build(dof_handler);
446 level_matrix_free_ = &matrix_free;
447 scalar_vector.resize(matrix_free.min_level(), matrix_free.max_level());
448 for (
unsigned int level = matrix_free.min_level();
449 level <= matrix_free.max_level();
451 matrix_free[level].initialize_dof_vector(scalar_vector[level]);
458 for (
unsigned int block = 0; block < src.n_blocks(); ++block)
459 transfer_.prolongate(to_level, dst.block(block), src.block(block));
466 for (
unsigned int block = 0; block < src.n_blocks(); ++block)
467 transfer_.restrict_and_add(
468 to_level, dst.block(block), src.block(block));
471 template <
typename Number2>
473 const dealii::DoFHandler<dim> &dof_handler,
474 dealii::MGLevelObject<scalar_type> &dst,
475 const dealii::LinearAlgebra::distributed::Vector<Number2> &src)
const
477 if (dst[dst.min_level()].size() == 0)
478 for (
unsigned int l = dst.min_level(); l <= dst.max_level(); ++l)
479 (*level_matrix_free_)[l].initialize_dof_vector(dst[l]);
480 transfer_.interpolate_to_mg(dof_handler, dst, src);
483 template <
typename Number2>
486 dealii::MGLevelObject<vector_type> &dst,
487 const dealii::LinearAlgebra::distributed::BlockVector<Number2>
490 if (dst[dst.min_level()].size() == 0)
491 for (
unsigned int l = dst.min_level(); l <= dst.max_level(); ++l) {
492 dst[l].reinit(src.n_blocks());
493 for (
unsigned int block = 0; block < src.n_blocks(); ++block)
494 (*level_matrix_free_)[l].initialize_dof_vector(
495 dst[l].block(block));
496 dst[l].collect_sizes();
499 for (
unsigned int block = 0; block < src.n_blocks(); ++block) {
500 transfer_.copy_to_mg(dof_handler, scalar_vector, src.block(block));
501 for (
unsigned int level = dst.min_level(); level <= dst.max_level();
503 dst[level].block(block).copy_locally_owned_data_from(
504 scalar_vector[level]);
508 template <
typename Number2>
510 const dealii::DoFHandler<dim> &dof_handler,
511 dealii::LinearAlgebra::distributed::BlockVector<Number2> &dst,
512 const dealii::MGLevelObject<vector_type> &src)
const
514 for (
unsigned int block = 0; block < dst.n_blocks(); ++block) {
515 for (
unsigned int level = src.min_level(); level <= src.max_level();
517 scalar_vector[level].copy_locally_owned_data_from(
518 src[level].block(block));
519 transfer_.copy_from_mg(dof_handler, dst.block(block), scalar_vector);
524 dealii::MGTransferMatrixFree<dim, Number> transfer_;
525 const dealii::MGLevelObject<dealii::MatrixFree<dim, Number>>
527 mutable dealii::MGLevelObject<scalar_type> scalar_vector;
534 template <
int dim,
typename Number,
typename Number2>
542 using vector_type = dealii::LinearAlgebra::distributed::Vector<Number>;
548 const dealii::MatrixFree<dim, Number> &matrix_free,
549 const dealii::LinearAlgebra::distributed::Vector<Number> &density,
550 const Number time_factor,
551 const unsigned int level = dealii::numbers::invalid_unsigned_int)
553 offline_data_ = &offline_data;
554 matrix_free_ = &matrix_free;
556 factor_ = time_factor;
565 dealii::types::global_dof_index
m()
const
567 return density_->size();
570 Number
el(
const unsigned int,
const unsigned int)
const
572 Assert(
false, dealii::ExcNotImplemented());
581 const auto get_lumped_mass = [&](
auto sentinel,
unsigned int i) {
582 using T =
decltype(sentinel);
583 if constexpr (std::is_same_v<Number, Number2>) {
584 if constexpr (std::is_same_v<Number, float>) {
585 if (level_ == dealii::numbers::invalid_unsigned_int) {
588 return lumped.template read_entry<T>(i);
590 const auto &level_lumped =
592 return read_entry<T>(level_lumped, i);
595 Assert(level_ == dealii::numbers::invalid_unsigned_int,
596 dealii::ExcInternalError());
599 return lumped.template read_entry<T>(i);
602 const auto &level_lumped =
604 return read_entry<T>(level_lumped, i);
608 const unsigned int n_owned =
609 dst.get_partitioner()->locally_owned_size();
611 const auto body_mass = [&](
auto sentinel,
unsigned int i) {
612 using T =
decltype(sentinel);
613 const auto m_i = get_lumped_mass(T(), i);
614 const auto rho_i = read_entry<T>(*density_, i);
615 const auto e_i = read_entry<T>(src, i);
616 write_entry<T>(dst, m_i * rho_i * e_i, i);
619 cpu_simd_loop<Number>(
"", body_mass, 0, n_owned, n_owned);
623 const auto integrator = [
this](
const auto &data,
627 dealii::FEEvaluation<dim, order_fe, order_quad, 1, Number> energy(
630 for (
unsigned int cell = range.first; cell < range.second; ++cell) {
632 energy.read_dof_values(src);
633 apply_local_operator(energy);
634 energy.distribute_local_to_global(dst);
638 matrix_free_->template cell_loop<vector_type, vector_type>(
639 integrator, dst, src,
false);
643 const auto &boundary_map =
644 (level_ == dealii::numbers::invalid_unsigned_int)
648 for (
auto entry : boundary_map) {
649 const auto i = std::get<0>(entry);
653 const auto id = std::get<4>(entry);
655 dst.local_element(i) = src.local_element(i);
660 std::shared_ptr<dealii::DiagonalMatrix<vector_type>> &matrix)
const
662 Assert(level_ != dealii::numbers::invalid_unsigned_int,
663 dealii::ExcNotImplemented());
664 matrix = std::make_shared<dealii::DiagonalMatrix<vector_type>>();
666 matrix_free_->initialize_dof_vector(vector);
671 dealii::MatrixFreeTools::compute_diagonal(
674 &EnergyMatrix::template apply_local_operator<
675 dealii::FEEvaluation<dim, -1, 0, dim, Number>>,
678 const unsigned int n_owned =
679 lumped_mass_matrix.get_partitioner()->locally_owned_size();
681 const auto body_invert = [&](
auto sentinel,
const unsigned int i) {
682 using T =
decltype(sentinel);
684 const auto m_i = read_entry<T>(lumped_mass_matrix, i);
685 const auto rho_i = read_entry<T>(*density_, i);
687 vector, Number(1.) / (m_i * rho_i + read_entry<T>(vector, i)), i);
689 cpu_simd_loop<Number>(
"", body_invert, 0, n_owned, n_owned);
693 for (
auto entry : boundary_map) {
694 const auto i = std::get<0>(entry);
698 const auto id = std::get<4>(entry);
700 vector.local_element(i) = 1.;
706 const dealii::MatrixFree<dim, Number> *matrix_free_;
707 const dealii::LinearAlgebra::distributed::Vector<Number> *density_;
711 template <
typename Evaluator>
712 void apply_local_operator(Evaluator &energy)
const
714 energy.evaluate(dealii::EvaluationFlags::gradients);
715 for (
unsigned int q = 0; q < energy.n_q_points; ++q) {
716 energy.submit_gradient(factor_ * energy.get_gradient(q), q);
718 energy.integrate(dealii::EvaluationFlags::gradients);
726 template <
int dim,
typename Number>
730 void build(
const dealii::DoFHandler<dim> &dof_handler,
731 const dealii::MGLevelObject<dealii::MatrixFree<dim, Number>>
740 const auto n_levels = dof_handler.get_triangulation().n_global_levels();
741 std::vector<std::shared_ptr<const dealii::Utilities::MPI::Partitioner>>
742 partitioners(n_levels);
744 for (
unsigned int level = matrix_free.min_level();
745 level <= matrix_free.max_level();
747 dealii::LinearAlgebra::distributed::Vector<Number> vector;
748 matrix_free[level].initialize_dof_vector(vector);
749 partitioners[level] = vector.get_partitioner();
752 dealii::MGTransferMatrixFree<dim, Number>::build(dof_handler,
760#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
vector_type & get_vector()
dealii::LinearAlgebra::distributed::BlockVector< Number > block_vector_type
block_vector_type & get_block_vector()
void reinit(const Vector &lumped_mass_matrix, const vector_type &density, const dealii::AffineConstraints< Number > &affine_constraints)
static constexpr unsigned int order_quad
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
static constexpr unsigned int order_fe
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 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
MGTransferVelocity()=default
void vmult(block_vector_type &dst, const block_vector_type &src) const
static constexpr unsigned int order_quad
dealii::LinearAlgebra::distributed::BlockVector< Number > block_vector_type
void initialize(const ParabolicSystem ¶bolic_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
static constexpr unsigned int order_fe
void Tvmult(block_vector_type &dst, const block_vector_type &src) const
void compute_diagonal(std::shared_ptr< DiagonalMatrix< dim, Number > > &matrix) const
const auto & boundary_map() const
const auto & level_lumped_mass_matrix() const
const auto & level_boundary_map() const
const auto & lumped_mass_matrix() const
DEAL_II_ALWAYS_INLINE void write_entry(V &vector, const T &values, unsigned int i)