8#include <compile_time_options.h>
9#include <deal.II/base/exceptions.h>
12#if DEAL_II_VERSION_GTE(9, 6, 0)
16#include <deal.II/base/config.h>
17#include <deal.II/distributed/tria.h>
18#include <deal.II/dofs/dof_accessor.h>
19#include <deal.II/dofs/dof_tools.h>
20#if DEAL_II_VERSION_GTE(9, 6, 0)
21#include <deal.II/grid/cell_status.h>
23#include <deal.II/grid/tria_accessor.h>
24#include <deal.II/grid/tria_iterator.h>
25#include <deal.II/lac/block_vector.h>
26#include <deal.II/lac/la_parallel_block_vector.h>
27#include <deal.II/lac/la_parallel_vector.h>
28#include <deal.II/lac/petsc_block_vector.h>
29#include <deal.II/lac/petsc_vector.h>
30#include <deal.II/lac/trilinos_parallel_block_vector.h>
31#include <deal.II/lac/trilinos_vector.h>
32#include <deal.II/lac/vector.h>
33#include <deal.II/matrix_free/fe_point_evaluation.h>
37 template <
typename Description,
int dim,
typename Number>
43 const std::string &subsection )
44 : ParameterAcceptor(subsection)
45 , limiter_parameters_(subsection +
"/mass transfer limiter")
46 , mpi_ensemble_(mpi_ensemble)
47 , offline_data_(&offline_data)
48 , hyperbolic_system_(&hyperbolic_system)
49 , parabolic_system_(¶bolic_system)
50 , handle_(dealii::numbers::invalid_unsigned_int)
60 template <
typename state_type>
62 pack_state_values(
const std::vector<state_type> &state_values)
64 std::vector<char> buffer(
sizeof(state_type) * state_values.size());
65 std::memcpy(buffer.data(), state_values.data(), buffer.size());
73 template <
typename state_type>
74 std::vector<state_type> unpack_state_values(
75 const boost::iterator_range<std::vector<char>::const_iterator>
78 const std::size_t n_bytes = data_range.size();
79 Assert(n_bytes %
sizeof(
state_type) == 0, dealii::ExcInternalError());
80 std::vector<state_type> state_values(n_bytes /
sizeof(
state_type));
81 std::memcpy(state_values.data(),
89 template <
typename Description,
int dim,
typename Number>
90 inline DEAL_II_ALWAYS_INLINE
auto
91 SolutionTransfer<Description, dim, Number>::get_tensor(
92 const HyperbolicVector &U,
const dealii::types::global_dof_index global_i)
95 const auto &scalar_partitioner = offline_data_->scalar_partitioner();
96 const auto &affine_constraints = offline_data_->affine_constraints();
97 const auto local_i = scalar_partitioner->global_to_local(global_i);
98 if (affine_constraints.is_constrained(global_i)) {
100 const auto &line = *affine_constraints.get_constraint_entries(global_i);
101 for (
const auto &[global_k, c_k] : line) {
102 const auto local_k = scalar_partitioner->global_to_local(global_k);
103 result += c_k * U.get_tensor(local_k);
107 return U.get_tensor(local_i);
112 template <
typename Description,
int dim,
typename Number>
113 inline DEAL_II_ALWAYS_INLINE
void
114 SolutionTransfer<Description, dim, Number>::add_tensor(
116 const state_type &new_U_i,
117 const dealii::types::global_dof_index global_i)
119 const auto &scalar_partitioner = offline_data_->scalar_partitioner();
120 const auto local_i = scalar_partitioner->global_to_local(global_i);
121 U.add_tensor(new_U_i, local_i);
125 template <
typename Description,
int dim,
typename Number>
127 const StateVector &old_state_vector [[maybe_unused]])
131 <<
"SolutionTransfer<Description, dim, Number>::prepare_projection()"
135#if !DEAL_II_VERSION_GTE(9, 6, 0)
139 "The SolutionTransfer class needs deal.II version 9.6.0 or newer"));
142 const auto &discretization = offline_data_->discretization();
143 auto &triangulation = *discretization.triangulation_;
145 Assert(handle_ == dealii::numbers::invalid_unsigned_int,
147 "You can only add one solution per SolutionTransfer object."));
155 handle_ = triangulation.register_data_attach(
156 [
this, &old_state_vector](
const auto cell,
157 const dealii::CellStatus status) {
158 const auto &dof_handler = offline_data_->dof_handler();
159 const auto dof_cell =
typename dealii::DoFHandler<dim>::cell_iterator(
160 &cell->get_triangulation(),
165 const auto &scalar_partitioner = offline_data_->scalar_partitioner();
167 const auto &U = std::get<0>(old_state_vector);
169 const auto &precomputed = std::get<1>(old_state_vector);
173 *hyperbolic_system_, limiter_parameters_, precomputed);
179 const auto n_dofs_per_cell = dof_cell->get_fe().n_dofs_per_cell();
180 std::vector<state_type> state_values(n_dofs_per_cell);
183 case dealii::CellStatus::cell_will_persist:
185 case dealii::CellStatus::cell_will_be_refined: {
191 Assert(dof_cell->is_active(), dealii::ExcInternalError());
192 std::vector<dealii::types::global_dof_index> dof_indices(
194 dof_cell->get_dof_indices(dof_indices);
197 std::begin(dof_indices),
198 std::end(dof_indices),
199 std::begin(state_values),
200 [&](
const auto global_i) {
return get_tensor(U, global_i); });
203 case dealii::CellStatus::children_will_be_coarsened: {
210 Assert(dof_cell->has_children(), dealii::ExcInternalError());
212 const auto &discretization = offline_data_->discretization();
213 const auto index = dof_cell->active_fe_index();
214 const auto &finite_element = discretization.finite_element()[index];
215 const auto &mapping = discretization.mapping()[index];
216 const auto &quadrature = discretization.quadrature()[index];
218 dealii::FEValues<dim> fe_values(
222 dealii::update_values | dealii::update_JxW_values |
223 dealii::update_quadrature_points);
225 const auto polynomial_space =
226 dealii::internal::FEPointEvaluation::get_polynomial_space(
229 std::vector<dealii::Point<dim, Number>> unit_points(
235 std::vector<dealii::Point<dim>> unit_points_temp(
236 std::is_same_v<Number, float> ? quadrature.size() : 0);
240 std::vector<state_type> state_values_quad(quadrature.size());
241 std::vector<state_type> local_rhs(n_dofs_per_cell);
243 std::vector<dealii::types::global_dof_index> dof_indices(
248 for (
unsigned int child = 0; child < dof_cell->n_children();
250 const auto child_cell = dof_cell->child(child);
252 Assert(child_cell->is_active(), dealii::ExcInternalError());
253 Assert(dof_cell->active_fe_index() ==
254 child_cell->active_fe_index(),
255 dealii::ExcMessage(
"SolutionTransfer: projection between "
256 "different FE space is not set up."));
258 fe_values.reinit(child_cell);
260 if constexpr (std::is_same_v<Number, float>) {
261 mapping.transform_points_real_to_unit_cell(
263 fe_values.get_quadrature_points(),
265 std::transform(std::begin(unit_points_temp),
266 std::end(unit_points_temp),
267 std::begin(unit_points),
268 [](
const auto &x) {
return x; });
270 mapping.transform_points_real_to_unit_cell(
271 dof_cell, fe_values.get_quadrature_points(), unit_points);
274 child_cell->get_dof_indices(dof_indices);
278 std::begin(dof_indices) != std::end(dof_indices)) {
279 const auto global_i = dof_indices[0];
280 const auto U_i = get_tensor(U, global_i);
282 scalar_partitioner->global_to_local(global_i);
286 for (
auto &it : state_values_quad)
289 for (
unsigned int i = 0; i < n_dofs_per_cell; ++i) {
290 const auto global_i = dof_indices[i];
291 const auto U_i = get_tensor(U, global_i);
293 scalar_partitioner->global_to_local(global_i);
294 const auto bounds_i =
298 for (
unsigned int q = 0; q < quadrature.size(); ++q) {
299 state_values_quad[q] += U_i * fe_values.shape_value(i, q);
303 for (
unsigned int q = 0; q < quadrature.size(); ++q)
304 state_values_quad[q] *= fe_values.JxW(q);
306 for (
unsigned int q = 0; q < quadrature.size(); ++q) {
307 const unsigned int n_shapes = polynomial_space.size();
308 AssertIndexRange(n_shapes, 10);
309 dealii::ndarray<Number, 10, 2, dim> shapes;
311 std::array<Number, dim> point;
312 for (
unsigned int d = 0; d < dim; ++d)
313 point[d] = unit_points[q][d];
314 for (
unsigned int i = 0; i < n_shapes; ++i)
315 polynomial_space[i].values_of_array(point, 1, &shapes[i][0]);
317 Assert(finite_element.degree == 1, dealii::ExcNotImplemented());
325 state_values_quad[q],
334 fe_values.reinit(dof_cell);
336 dealii::FullMatrix<double> mass_inverse(n_dofs_per_cell,
338 dealii::Vector<double> lumped_mass(n_dofs_per_cell);
339 dealii::Vector<double> lumped_mass_inverse(n_dofs_per_cell);
341 auto total_mass = Number(0.);
342 for (
unsigned int i = 0; i < n_dofs_per_cell; ++i) {
343 for (
unsigned int j = 0; j < n_dofs_per_cell; ++j) {
345 for (
unsigned int q = 0; q < quadrature.size(); ++q)
346 sum += fe_values.shape_value(i, q) *
347 fe_values.shape_value(j, q) * fe_values.JxW(q);
348 mass_inverse(i, j) = sum;
349 lumped_mass(i) += sum;
351 lumped_mass_inverse(i) = Number(1.) / lumped_mass(i);
352 total_mass += lumped_mass(i);
354 mass_inverse.gauss_jordan();
360 std::vector<state_type> pij_matrix(n_dofs_per_cell *
362 dealii::FullMatrix<Number> lij_matrix(n_dofs_per_cell,
365 const auto kappa_inverse = Number(n_dofs_per_cell);
366 const auto kappa = Number(1.) / kappa_inverse;
368 for (
unsigned int i = 0; i < n_dofs_per_cell; ++i) {
369 const state_type U_i = lumped_mass_inverse(i) * local_rhs[i];
370 state_values[i] = U_i;
372 for (
unsigned int j = 0; j < n_dofs_per_cell; ++j) {
373 const auto kronecker_ij = Number(i == j ? 1. : 0.);
375 lumped_mass(i) * mass_inverse(i, j) - kronecker_ij;
377 lumped_mass(j) * mass_inverse(i, j) - kronecker_ij;
378 const auto P_ij = kappa_inverse * lumped_mass_inverse(i) *
379 (b_ij * local_rhs[j] - b_ji * local_rhs[i]);
380 pij_matrix[n_dofs_per_cell * i + j] = P_ij;
386 const auto n_iterations = limiter_parameters_.iterations();
387 for (
unsigned int pass = 0; pass < n_iterations; ++pass) {
389 for (
unsigned int i = 0; i < n_dofs_per_cell; ++i) {
390 const auto &U_i = state_values[i];
392 for (
unsigned int j = 0; j < n_dofs_per_cell; ++j) {
393 const auto &P_ij = pij_matrix[n_dofs_per_cell * i + j];
394 const auto &[l_ij, check] = limiter.
limit(bounds, U_i, P_ij);
395 lij_matrix(i, j) = l_ij;
397#ifdef DEBUG_EXPENSIVE_BOUNDS_CHECK
401 "Error: low-order state out of bounds in "
402 "register_data_attach / children_will_be_coarsened"));
407 for (
unsigned int i = 0; i < n_dofs_per_cell; ++i) {
408 auto &U_i = state_values[i];
410 for (
unsigned int j = 0; j < n_dofs_per_cell; ++j) {
412 std::min(lij_matrix(i, j), lij_matrix(j, i));
413 auto &P_ij = pij_matrix[n_dofs_per_cell * i + j];
414 U_i += kappa * l_ij * P_ij;
418#ifdef DEBUG_EXPENSIVE_BOUNDS_CHECK
420 hyperbolic_system_->template view<dim, Number>();
422 view.is_admissible(U_i),
424 "Error: inadmissible state encountered in "
425 "register_data_attach / children_will_be_coarsened"));
431 case dealii::CellStatus::cell_invalid:
432 Assert(
false, dealii::ExcInternalError());
437 return pack_state_values(state_values);
444 template <
typename Description,
int dim,
typename Number>
449 std::cout <<
"SolutionTransfer<Description, dim, Number>::project()"
453#if !DEAL_II_VERSION_GTE(9, 6, 0)
457 "The SolutionTransfer class needs deal.II version 9.6.0 or newer"));
461 const auto &discretization = offline_data_->discretization();
462 auto &triangulation = *discretization.triangulation_;
465 handle_ != dealii::numbers::invalid_unsigned_int,
467 "Cannot project() a state vector without valid handle. "
468 "prepare_projection() or set_handle() have to be called first."));
470 const auto &scalar_partitioner = offline_data_->scalar_partitioner();
471 const auto &affine_constraints = offline_data_->affine_constraints();
472 const auto n_locally_owned = offline_data_->n_locally_owned();
476 projected_mass.reinit(offline_data_->scalar_partitioner());
478 projected_state.reinit(offline_data_->hyperbolic_vector_partitioner());
486 std::map<std::tuple<
unsigned int ,
unsigned int >,
state_type>
488 std::map<
unsigned int ,
Bounds> bounds_map;
491 kappa.reinit(offline_data_->scalar_partitioner());
499 triangulation.notify_ready_to_unpack(
501 [
this, &projected_mass, &projected_state](
503 const dealii::CellStatus status,
504 const auto &data_range) {
505 const auto &dof_handler = offline_data_->dof_handler();
506 const auto dof_cell =
typename dealii::DoFHandler<dim>::cell_iterator(
507 &cell->get_triangulation(),
516 const auto n_dofs_per_cell = dof_cell->get_fe().n_dofs_per_cell();
517 std::vector<dealii::types::global_dof_index> dof_indices(
520 const auto state_values = unpack_state_values<state_type>(data_range);
523 case dealii::CellStatus::cell_will_persist:
525 case dealii::CellStatus::children_will_be_coarsened: {
531 Assert(dof_cell->is_active(), dealii::ExcInternalError());
532 dof_cell->get_dof_indices(dof_indices);
534 const auto &discretization = offline_data_->discretization();
535 const auto index = dof_cell->active_fe_index();
536 const auto &finite_element = discretization.finite_element()[index];
537 const auto &mapping = discretization.mapping()[index];
538 const auto &quadrature = discretization.quadrature()[index];
540 dealii::FEValues<dim> fe_values(mapping,
543 dealii::update_values |
544 dealii::update_JxW_values);
546 fe_values.reinit(dof_cell);
548 dealii::Vector<double> mi(n_dofs_per_cell);
549 for (
unsigned int i = 0; i < n_dofs_per_cell; ++i) {
551 for (
unsigned int q = 0; q < quadrature.size(); ++q)
552 sum += fe_values.shape_value(i, q) * fe_values.JxW(q);
556 for (
unsigned int i = 0; i < n_dofs_per_cell; ++i) {
557 const auto global_i = dof_indices[i];
558 add_tensor(projected_state, mi(i) * state_values[i], global_i);
559 projected_mass(global_i) += mi(i);
564 case dealii::CellStatus::cell_will_be_refined: {
570 Assert(dof_cell->has_children(), dealii::ExcInternalError());
572 const auto &discretization = offline_data_->discretization();
573 const auto index = dof_cell->active_fe_index();
574 const auto &finite_element = discretization.finite_element()[index];
575 const auto &mapping = discretization.mapping()[index];
576 const auto &quadrature = discretization.quadrature()[index];
578 dealii::FEValues<dim> fe_values(
582 dealii::update_values | dealii::update_JxW_values |
583 dealii::update_quadrature_points);
585 const auto polynomial_space =
586 dealii::internal::FEPointEvaluation::get_polynomial_space(
588 std::vector<dealii::Point<dim, Number>> unit_points(
594 std::vector<dealii::Point<dim>> unit_points_temp(
595 std::is_same_v<Number, float> ? quadrature.size() : 0);
597 dealii::FullMatrix<double> mass_inverse(n_dofs_per_cell,
599 dealii::Vector<double> lumped_mass(n_dofs_per_cell);
600 std::vector<state_type> local_rhs(n_dofs_per_cell);
602 for (
unsigned int child = 0; child < dof_cell->n_children();
604 const auto child_cell = dof_cell->child(child);
606 Assert(child_cell->is_active(), dealii::ExcInternalError());
607 Assert(dof_cell->active_fe_index() ==
608 child_cell->active_fe_index(),
609 dealii::ExcMessage(
"SolutionTransfer: projection between "
610 "different FE space is not set up."));
612 child_cell->get_dof_indices(dof_indices);
616 fe_values.reinit(child_cell);
618 if constexpr (std::is_same_v<Number, float>) {
619 mapping.transform_points_real_to_unit_cell(
621 fe_values.get_quadrature_points(),
623 std::transform(std::begin(unit_points_temp),
624 std::end(unit_points_temp),
625 std::begin(unit_points),
626 [](
const auto &x) {
return x; });
628 mapping.transform_points_real_to_unit_cell(
629 dof_cell, fe_values.get_quadrature_points(), unit_points);
632 for (
auto &it : local_rhs)
635 for (
unsigned int q = 0; q < quadrature.size(); ++q) {
636 Assert(finite_element.degree == 1, dealii::ExcNotImplemented());
638 dealii::internal::evaluate_tensor_product_value(
640 make_const_array_view(state_values),
643 coefficient *= fe_values.JxW(q);
645 for (
unsigned int i = 0; i < n_dofs_per_cell; ++i)
646 local_rhs[i] += coefficient * fe_values.shape_value(i, q);
651 mass_inverse = Number(0.);
652 lumped_mass = Number(0.);
653 for (
unsigned int i = 0; i < n_dofs_per_cell; ++i) {
654 for (
unsigned int j = 0; j < n_dofs_per_cell; ++j) {
656 for (
unsigned int q = 0; q < quadrature.size(); ++q)
657 sum += fe_values.shape_value(i, q) *
658 fe_values.shape_value(j, q) * fe_values.JxW(q);
659 mass_inverse(i, j) = sum;
660 lumped_mass(i) += sum;
663 mass_inverse.gauss_jordan();
667 for (
unsigned int i = 0; i < n_dofs_per_cell; ++i) {
669 for (
unsigned int j = 0; j < n_dofs_per_cell; ++j) {
670 U_i += mass_inverse(i, j) * local_rhs[j];
673#ifdef DEBUG_EXPENSIVE_BOUNDS_CHECK
675 hyperbolic_system_->template view<dim, Number>();
676 AssertThrow(view.is_admissible(U_i),
678 "Error: inadmissible state encountered in "
679 "ready_to_unpack / cell_will_be_refined"));
681 const auto global_i = dof_indices[i];
682 add_tensor(projected_state, lumped_mass(i) * U_i, global_i);
683 projected_mass(global_i) += lumped_mass(i);
689 case dealii::CellStatus::cell_invalid:
690 Assert(
false, dealii::ExcInternalError());
713 auto &new_U = std::get<0>(new_state_vector);
719 const auto update_new_state_vector = [&]() {
720 for (
unsigned int local_i = 0; local_i < n_locally_owned; ++local_i) {
722 const auto mU_i = projected_state.get_tensor(local_i);
723 const auto m_i = projected_mass.local_element(local_i);
725#ifdef DEBUG_EXPENSIVE_BOUNDS_CHECK
726 const auto view = hyperbolic_system_->template view<dim, Number>();
728 view.is_admissible(mU_i / m_i),
729 dealii::ExcMessage(
"Error: inadmissible state encountered in "
730 "update_new_state_vector()"));
733 new_U.write_tensor(mU_i / m_i, local_i);
735 new_U.update_ghost_values();
738 update_new_state_vector();
742 const auto &precomputed = std::get<1>(new_state_vector);
744 const auto update_precomputed_values = [&]() {
745 const unsigned int n_export_indices = offline_data_->n_export_indices();
746 const unsigned int n_internal = offline_data_->n_locally_internal();
747 const unsigned int n_owned = offline_data_->n_locally_owned();
748 const auto &sparsity_simd = offline_data_->sparsity_pattern_simd();
749 unsigned int channel = 10;
750 using VA = dealii::VectorizedArray<Number>;
751 static constexpr auto n_precomputation_cycles =
752 View::n_precomputation_cycles;
754 new_U.update_ghost_values();
756 if constexpr (n_precomputation_cycles != 0) {
757 for (
unsigned int cycle = 0; cycle < n_precomputation_cycles; ++cycle) {
760 precomputed.update_ghost_values_start(channel++);
761 precomputed.update_ghost_values_finish();
766 auto loop = [&](
auto sentinel,
768 unsigned int right) {
769 using T =
decltype(sentinel);
772 bool thread_ready =
false;
774 const auto view = hyperbolic_system_->template view<dim, T>();
775 view.precomputation_loop(
777 [&](
const unsigned int i) {
778 synchronization_dispatch.check(
779 thread_ready, i >= n_export_indices && i < n_internal);
789 loop(Number(), n_internal, n_owned);
791 loop(VA(), 0, n_internal);
798 update_precomputed_values();
802 *hyperbolic_system_, limiter_parameters_, precomputed);
815 for (
const auto &line : affine_constraints.get_lines()) {
816 const auto global_i = line.index;
817 const auto local_i = scalar_partitioner->global_to_local(global_i);
820 if (local_i >= n_locally_owned)
824 const auto m_i_star = projected_mass.local_element(local_i);
825 const auto U_i_star = projected_state.get_tensor(local_i) / m_i_star;
827 auto &bounds = bounds_map[local_i];
832 for (
const auto &[global_k, c_k] : line.entries) {
833 const auto local_k = scalar_partitioner->global_to_local(global_k);
834 U_i_interp += c_k * new_U.get_tensor(local_k);
838 for (
const auto &[global_k, c_k] : line.entries) {
839 const auto local_k = scalar_partitioner->global_to_local(global_k);
840 const auto U_k = new_U.get_tensor(local_k);
842 const auto bounds_k =
846 projected_state.add_tensor(c_k * m_i_star * U_i_star, local_k);
847 projected_mass.local_element(local_k) += c_k * m_i_star;
849 kappa.local_element(local_k) += Number(1.);
850 pik_matrix[{local_i, local_k}] = c_k * m_i_star * (U_k - U_i_interp);
858 update_new_state_vector();
861 projected_mass.update_ghost_values();
862 kappa.update_ghost_values();
866 const auto n_iterations = limiter_parameters_.iterations();
867 for (
unsigned int pass = 0; pass < n_iterations; ++pass) {
870 update_precomputed_values();
872 for (
const auto &line : affine_constraints.get_lines()) {
873 const auto global_i = line.index;
874 const auto local_i = scalar_partitioner->global_to_local(global_i);
877 if (local_i >= n_locally_owned)
893 auto &bounds = bounds_map[local_i];
894 auto total_mass = Number(0.);
895 for (
const auto &[global_k, c_k] : line.entries) {
896 const auto local_k = scalar_partitioner->global_to_local(global_k);
897 const auto U_k = new_U.get_tensor(local_k);
898 const auto bounds_k =
902 const auto m_k = projected_mass.local_element(local_k);
909 const auto relaxed_bounds =
914 for (
const auto &[global_k, c_k] : line.entries) {
915 const auto local_k = scalar_partitioner->global_to_local(global_k);
916 const auto kappa_k = kappa.local_element(local_k);
917 const auto m_k = projected_mass.local_element(local_k);
918 const auto U_k = new_U.get_tensor(local_k);
919 const auto P_ik = pik_matrix[{local_i, local_k}] * kappa_k / m_k;
921 const auto &[l_k, check] = limiter.
limit(relaxed_bounds, U_k, P_ik);
922#ifdef DEBUG_EXPENSIVE_BOUNDS_CHECK
924 dealii::ExcMessage(
"Error: low-order state out of bounds "
925 "in project / state redistribution"));
927 l = std::min(l, l_k);
932 for (
const auto &[global_k, c_k] : line.entries) {
933 const auto local_k = scalar_partitioner->global_to_local(global_k);
934 auto &mP_ik = pik_matrix[{local_i, local_k}];
935 projected_state.add_tensor(l * mP_ik, local_k);
942 update_new_state_vector();
946 for (
unsigned int local_i = 0; local_i < n_locally_owned; ++local_i) {
947 const auto global_i = scalar_partitioner->local_to_global(local_i);
948 if (affine_constraints.is_constrained(global_i))
951 new_U.update_ghost_values();
953#ifdef DEBUG_SYMMETRY_CHECK
957 const auto &lumped_mass_matrix = offline_data_->lumped_mass_matrix();
958 for (
unsigned int local_i = 0; local_i < n_locally_owned; ++local_i) {
959 const auto global_i = scalar_partitioner->local_to_global(local_i);
960 if (affine_constraints.is_constrained(global_i))
963 const auto m_i = projected_mass.local_element(local_i);
964 const auto m_i_reference = lumped_mass_matrix.local_element(local_i);
965 Assert(std::abs(m_i - m_i_reference) < 1.e-10,
967 "SolutionTransfer::projection(): something went wrong. Final "
968 "masses do not agree with those computed in OfflineData."));
Bounds fully_relax_bounds(const Bounds &bounds, const Number &hd) const
Bounds combine_bounds(const Bounds &bounds_left, const Bounds &bounds_right) const
Bounds projection_bounds_from_state(const unsigned int i, const state_type &U_i) const
std::tuple< Number, bool > limit(const Bounds &bounds, const state_type &U, const state_type &P, const Number t_min=Number(0.), const Number t_max=Number(1.)) const
typename View::HyperbolicVector HyperbolicVector
SolutionTransfer(const MPIEnsemble &mpi_ensemble, const OfflineData< dim, Number > &offline_data, const HyperbolicSystem &hyperbolic_system, const ParabolicSystem ¶bolic_system, const std::string &subsection="/SolutionTransfer")
Vectors::ScalarVector< Number > ScalarVector
Description::template Limiter< dim, double >::Bounds Bounds
typename View::state_type state_type
typename Description::ParabolicSystem ParabolicSystem
void prepare_projection(const StateVector &old_state_vector)
typename Description::HyperbolicSystem HyperbolicSystem
void project(StateVector &new_state_vector)
typename View::StateVector StateVector
#define RYUJIN_PARALLEL_REGION_BEGIN
#define RYUJIN_PARALLEL_REGION_END
void integrate_tensor_product_value(const dealii::ndarray< Number, 2, dim > *shapes, const unsigned int n_shapes, const Number2 &value, Number2 *values, const dealii::Point< dim, Number > &p, const bool do_add)
DEAL_II_ALWAYS_INLINE FT add(const FT &flux_left_ij, const FT &flux_right_ij)