ryujin 2.1.1 revision 350e54cc11f3d21282bddcf3e6517944dbc508bf
solution_transfer.template.h
Go to the documentation of this file.
1//
2// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception or LGPL-2.1-or-later
3// Copyright (C) 2024 - 2024 by the ryujin authors
4//
5
6#pragma once
7
8#include <compile_time_options.h>
9#include <deal.II/base/exceptions.h>
10
11#include "solution_transfer.h"
12#if DEAL_II_VERSION_GTE(9, 6, 0)
14#endif
15
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>
22#endif
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>
34
35namespace ryujin
36{
37 template <typename Description, int dim, typename Number>
39 const MPIEnsemble &mpi_ensemble,
40 const OfflineData<dim, Number> &offline_data,
41 const HyperbolicSystem &hyperbolic_system,
42 const ParabolicSystem &parabolic_system,
43 const std::string &subsection /* = "/SolutionTransfer" */)
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_(&parabolic_system)
50 , handle_(dealii::numbers::invalid_unsigned_int)
51 {
52 }
53
54
55 namespace
56 {
60 template <typename state_type>
61 std::vector<char>
62 pack_state_values(const std::vector<state_type> &state_values)
63 {
64 std::vector<char> buffer(sizeof(state_type) * state_values.size());
65 std::memcpy(buffer.data(), state_values.data(), buffer.size());
66 return buffer;
67 }
68
69
73 template <typename state_type>
74 std::vector<state_type> unpack_state_values(
75 const boost::iterator_range<std::vector<char>::const_iterator>
76 &data_range)
77 {
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(),
82 &data_range[0],
83 state_values.size() * sizeof(state_type));
84 return state_values;
85 }
86 } // namespace
87
88
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)
93 -> state_type
94 {
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)) {
99 state_type result;
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);
104 }
105 return result;
106 } else {
107 return U.get_tensor(local_i);
108 }
109 }
110
111
112 template <typename Description, int dim, typename Number>
113 inline DEAL_II_ALWAYS_INLINE void
114 SolutionTransfer<Description, dim, Number>::add_tensor(
115 HyperbolicVector &U,
116 const state_type &new_U_i,
117 const dealii::types::global_dof_index global_i)
118 {
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);
122 }
123
124
125 template <typename Description, int dim, typename Number>
127 const StateVector &old_state_vector [[maybe_unused]])
128 {
129#ifdef DEBUG_OUTPUT
130 std::cout
131 << "SolutionTransfer<Description, dim, Number>::prepare_projection()"
132 << std::endl;
133#endif
134
135#if !DEAL_II_VERSION_GTE(9, 6, 0)
136 AssertThrow(
137 false,
138 dealii::ExcMessage(
139 "The SolutionTransfer class needs deal.II version 9.6.0 or newer"));
140
141#else
142 const auto &discretization = offline_data_->discretization();
143 auto &triangulation = *discretization.triangulation_; /* writable */
144
145 Assert(handle_ == dealii::numbers::invalid_unsigned_int,
146 dealii::ExcMessage(
147 "You can only add one solution per SolutionTransfer object."));
148
149 /*
150 * -----------------------------------------------------------------------
151 * Cell-level projection to parent cells and packing data:
152 * -----------------------------------------------------------------------
153 */
154
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(),
161 cell->level(),
162 cell->index(),
163 &dof_handler);
164
165 const auto &scalar_partitioner = offline_data_->scalar_partitioner();
166
167 const auto &U = std::get<0>(old_state_vector);
168 /* precomputed needs to be valid for bounds computation */
169 const auto &precomputed = std::get<1>(old_state_vector);
170
171 using Limiter = typename Description::template Limiter<dim, Number>;
172 const Limiter limiter(
173 *hyperbolic_system_, limiter_parameters_, precomputed);
174
175 /*
176 * Collect state values for packing:
177 */
178
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);
181
182 switch (status) {
183 case dealii::CellStatus::cell_will_persist:
184 [[fallthrough]];
185 case dealii::CellStatus::cell_will_be_refined: {
186 /*
187 * For both cases we need state values from the currently
188 * active cell:
189 */
190
191 Assert(dof_cell->is_active(), dealii::ExcInternalError());
192 std::vector<dealii::types::global_dof_index> dof_indices(
193 n_dofs_per_cell);
194 dof_cell->get_dof_indices(dof_indices);
195
196 std::transform(
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); });
201 } break;
202
203 case dealii::CellStatus::children_will_be_coarsened: {
204 /*
205 * We need to project values from the active child cells up to
206 * the present parent cell that will become active after
207 * coarsening.
208 */
209
210 Assert(dof_cell->has_children(), dealii::ExcInternalError());
211
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];
217
218 dealii::FEValues<dim> fe_values(
219 mapping,
220 finite_element,
221 quadrature,
222 dealii::update_values | dealii::update_JxW_values |
223 dealii::update_quadrature_points);
224
225 const auto polynomial_space =
226 dealii::internal::FEPointEvaluation::get_polynomial_space(
227 finite_element);
228
229 std::vector<dealii::Point<dim, Number>> unit_points(
230 quadrature.size());
231 /*
232 * for Number == float we need a temporary vector for the
233 * transform_points_real_to_unit_cell() function:
234 */
235 std::vector<dealii::Point<dim>> unit_points_temp(
236 std::is_same_v<Number, float> ? quadrature.size() : 0);
237
238 /* Step 1: build up right hand side by iterating over children: */
239
240 std::vector<state_type> state_values_quad(quadrature.size());
241 std::vector<state_type> local_rhs(n_dofs_per_cell);
242
243 std::vector<dealii::types::global_dof_index> dof_indices(
244 n_dofs_per_cell);
245
246 Bounds bounds;
247
248 for (unsigned int child = 0; child < dof_cell->n_children();
249 ++child) {
250 const auto child_cell = dof_cell->child(child);
251
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."));
257
258 fe_values.reinit(child_cell);
259
260 if constexpr (std::is_same_v<Number, float>) {
261 mapping.transform_points_real_to_unit_cell(
262 dof_cell,
263 fe_values.get_quadrature_points(),
264 unit_points_temp);
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; });
269 } else {
270 mapping.transform_points_real_to_unit_cell(
271 dof_cell, fe_values.get_quadrature_points(), unit_points);
272 }
273
274 child_cell->get_dof_indices(dof_indices);
275
276 /* We want a "left fold first" for the bounds: */
277 if (child == 0 &&
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);
281 const auto local_i =
282 scalar_partitioner->global_to_local(global_i);
283 bounds = limiter.projection_bounds_from_state(local_i, U_i);
284 }
285
286 for (auto &it : state_values_quad)
287 it = state_type{};
288
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);
292 const auto local_i =
293 scalar_partitioner->global_to_local(global_i);
294 const auto bounds_i =
295 limiter.projection_bounds_from_state(local_i, U_i);
296 bounds = limiter.combine_bounds(bounds, bounds_i);
297
298 for (unsigned int q = 0; q < quadrature.size(); ++q) {
299 state_values_quad[q] += U_i * fe_values.shape_value(i, q);
300 }
301 }
302
303 for (unsigned int q = 0; q < quadrature.size(); ++q)
304 state_values_quad[q] *= fe_values.JxW(q);
305
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;
310 // Evaluate 1d polynomials and their derivatives
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]);
316
317 Assert(finite_element.degree == 1, dealii::ExcNotImplemented());
318
320 /*is linear*/ true,
321 dim,
322 Number,
323 state_type>(shapes.data(),
324 n_shapes,
325 state_values_quad[q],
326 local_rhs.data(),
327 unit_points[q],
328 true);
329 }
330 }
331
332 /* Step 2: construct inverse mass matrices: */
333
334 fe_values.reinit(dof_cell);
335
336 dealii::FullMatrix<double> mass_inverse(n_dofs_per_cell,
337 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);
340
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) {
344 double sum = 0;
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;
350 }
351 lumped_mass_inverse(i) = Number(1.) / lumped_mass(i);
352 total_mass += lumped_mass(i);
353 }
354 mass_inverse.gauss_jordan();
355
356 /* Step 3: compute low-order update and P_ij matrix: */
357
358 bounds = limiter.fully_relax_bounds(bounds, total_mass);
359
360 std::vector<state_type> pij_matrix(n_dofs_per_cell *
361 n_dofs_per_cell);
362 dealii::FullMatrix<Number> lij_matrix(n_dofs_per_cell,
363 n_dofs_per_cell);
364
365 const auto kappa_inverse = Number(n_dofs_per_cell);
366 const auto kappa = Number(1.) / kappa_inverse;
367
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;
371
372 for (unsigned int j = 0; j < n_dofs_per_cell; ++j) {
373 const auto kronecker_ij = Number(i == j ? 1. : 0.);
374 const auto b_ij =
375 lumped_mass(i) * mass_inverse(i, j) - kronecker_ij;
376 const auto b_ji =
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;
381 }
382 }
383
384 /* Step 4: compute l_ij matrix and apply limited update: */
385
386 const auto n_iterations = limiter_parameters_.iterations();
387 for (unsigned int pass = 0; pass < n_iterations; ++pass) {
388
389 for (unsigned int i = 0; i < n_dofs_per_cell; ++i) {
390 const auto &U_i = state_values[i];
391
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;
396
397#ifdef DEBUG_EXPENSIVE_BOUNDS_CHECK
398 AssertThrow(
399 check,
400 dealii::ExcMessage(
401 "Error: low-order state out of bounds in "
402 "register_data_attach / children_will_be_coarsened"));
403#endif
404 }
405 }
406
407 for (unsigned int i = 0; i < n_dofs_per_cell; ++i) {
408 auto &U_i = state_values[i];
409
410 for (unsigned int j = 0; j < n_dofs_per_cell; ++j) {
411 const auto l_ij =
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;
415 P_ij -= l_ij * P_ij;
416 }
417
418#ifdef DEBUG_EXPENSIVE_BOUNDS_CHECK
419 const auto view =
420 hyperbolic_system_->template view<dim, Number>();
421 AssertThrow(
422 view.is_admissible(U_i),
423 dealii::ExcMessage(
424 "Error: inadmissible state encountered in "
425 "register_data_attach / children_will_be_coarsened"));
426#endif
427 }
428 }
429 } break;
430
431 case dealii::CellStatus::cell_invalid:
432 Assert(false, dealii::ExcInternalError());
433 __builtin_trap();
434 break;
435 }
436
437 return pack_state_values(state_values);
438 },
439 /* returns_variable_size_data =*/false);
440#endif
441 }
442
443
444 template <typename Description, int dim, typename Number>
446 StateVector &new_state_vector [[maybe_unused]])
447 {
448#ifdef DEBUG_OUTPUT
449 std::cout << "SolutionTransfer<Description, dim, Number>::project()"
450 << std::endl;
451#endif
452
453#if !DEAL_II_VERSION_GTE(9, 6, 0)
454 AssertThrow(
455 false,
456 dealii::ExcMessage(
457 "The SolutionTransfer class needs deal.II version 9.6.0 or newer"));
458
459#else
460
461 const auto &discretization = offline_data_->discretization();
462 auto &triangulation = *discretization.triangulation_; /* writable */
463
464 Assert(
465 handle_ != dealii::numbers::invalid_unsigned_int,
466 dealii::ExcMessage(
467 "Cannot project() a state vector without valid handle. "
468 "prepare_projection() or set_handle() have to be called first."));
469
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();
473
474
475 ScalarVector projected_mass;
476 projected_mass.reinit(offline_data_->scalar_partitioner());
477 HyperbolicVector projected_state;
478 projected_state.reinit(offline_data_->hyperbolic_vector_partitioner());
479
480 /*
481 * We only need to construct entries in a pik_matrix for a small subset
482 * of affected degrees of freedom for which we have to construct the
483 * entire pik_matrix first for the limiting process (in contrast to the
484 * entirely cell-local limiting done before). Let's simply use a map.
485 */
486 std::map<std::tuple<unsigned int /*i*/, unsigned int /*k*/>, state_type>
487 pik_matrix;
488 std::map<unsigned int /*i*/, Bounds> bounds_map;
489
490 ScalarVector kappa;
491 kappa.reinit(offline_data_->scalar_partitioner());
492
493 /*
494 * -----------------------------------------------------------------------
495 * Unpacking data and cell-level interpolation/projection to child cells:
496 * -----------------------------------------------------------------------
497 */
498
499 triangulation.notify_ready_to_unpack( //
500 handle_,
501 [this, &projected_mass, &projected_state](
502 const auto &cell,
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(),
508 cell->level(),
509 cell->index(),
510 &dof_handler);
511
512 /*
513 * Retrieve packed values and project onto cell:
514 */
515
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(
518 n_dofs_per_cell);
519
520 const auto state_values = unpack_state_values<state_type>(data_range);
521
522 switch (status) {
523 case dealii::CellStatus::cell_will_persist:
524 [[fallthrough]];
525 case dealii::CellStatus::children_will_be_coarsened: {
526 /*
527 * For both cases we distribute stored state_values to the
528 * projected_state and projected_mass vectors.
529 */
530
531 Assert(dof_cell->is_active(), dealii::ExcInternalError());
532 dof_cell->get_dof_indices(dof_indices);
533
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];
539
540 dealii::FEValues<dim> fe_values(mapping,
541 finite_element,
542 quadrature,
543 dealii::update_values |
544 dealii::update_JxW_values);
545
546 fe_values.reinit(dof_cell);
547
548 dealii::Vector<double> mi(n_dofs_per_cell);
549 for (unsigned int i = 0; i < n_dofs_per_cell; ++i) {
550 double sum = 0;
551 for (unsigned int q = 0; q < quadrature.size(); ++q)
552 sum += fe_values.shape_value(i, q) * fe_values.JxW(q);
553 mi(i) += sum;
554 }
555
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);
560 }
561
562 } break;
563
564 case dealii::CellStatus::cell_will_be_refined: {
565 /*
566 * We are on a (non active) cell that has been refined. Project
567 * onto the children and do a local mass projection there:
568 */
569
570 Assert(dof_cell->has_children(), dealii::ExcInternalError());
571
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];
577
578 dealii::FEValues<dim> fe_values(
579 mapping,
580 finite_element,
581 quadrature,
582 dealii::update_values | dealii::update_JxW_values |
583 dealii::update_quadrature_points);
584
585 const auto polynomial_space =
586 dealii::internal::FEPointEvaluation::get_polynomial_space(
587 finite_element);
588 std::vector<dealii::Point<dim, Number>> unit_points(
589 quadrature.size());
590 /*
591 * for Number == float we need a temporary vector for the
592 * transform_points_real_to_unit_cell() function:
593 */
594 std::vector<dealii::Point<dim>> unit_points_temp(
595 std::is_same_v<Number, float> ? quadrature.size() : 0);
596
597 dealii::FullMatrix<double> mass_inverse(n_dofs_per_cell,
598 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);
601
602 for (unsigned int child = 0; child < dof_cell->n_children();
603 ++child) {
604 const auto child_cell = dof_cell->child(child);
605
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."));
611
612 child_cell->get_dof_indices(dof_indices);
613
614 /* Step 1: build up right hand side on child cell: */
615
616 fe_values.reinit(child_cell);
617
618 if constexpr (std::is_same_v<Number, float>) {
619 mapping.transform_points_real_to_unit_cell(
620 dof_cell,
621 fe_values.get_quadrature_points(),
622 unit_points_temp);
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; });
627 } else {
628 mapping.transform_points_real_to_unit_cell(
629 dof_cell, fe_values.get_quadrature_points(), unit_points);
630 }
631
632 for (auto &it : local_rhs)
633 it = state_type{};
634
635 for (unsigned int q = 0; q < quadrature.size(); ++q) {
636 Assert(finite_element.degree == 1, dealii::ExcNotImplemented());
637 auto coefficient =
638 dealii::internal::evaluate_tensor_product_value(
639 polynomial_space,
640 make_const_array_view(state_values),
641 unit_points[q],
642 /*is linear*/ true);
643 coefficient *= fe_values.JxW(q);
644
645 for (unsigned int i = 0; i < n_dofs_per_cell; ++i)
646 local_rhs[i] += coefficient * fe_values.shape_value(i, q);
647 }
648
649 /* Step 2: solve with inverse mass matrix on child cell: */
650
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) {
655 double sum = 0;
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;
661 }
662 }
663 mass_inverse.gauss_jordan();
664
665 /* Step 3: compute high order update and write back: */
666
667 for (unsigned int i = 0; i < n_dofs_per_cell; ++i) {
668 state_type U_i;
669 for (unsigned int j = 0; j < n_dofs_per_cell; ++j) {
670 U_i += mass_inverse(i, j) * local_rhs[j];
671 }
672
673#ifdef DEBUG_EXPENSIVE_BOUNDS_CHECK
674 const auto view =
675 hyperbolic_system_->template view<dim, Number>();
676 AssertThrow(view.is_admissible(U_i),
677 dealii::ExcMessage(
678 "Error: inadmissible state encountered in "
679 "ready_to_unpack / cell_will_be_refined"));
680#endif
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);
684 }
685 } /*child*/
686
687 } break;
688
689 case dealii::CellStatus::cell_invalid:
690 Assert(false, dealii::ExcInternalError());
691 __builtin_trap();
692 break;
693 }
694 });
695
696 projected_mass.compress(dealii::VectorOperation::add);
697 projected_state.compress(dealii::VectorOperation::add);
698
699 /*
700 * -----------------------------------------------------------------------
701 * Redistribute masses to satisfy hanging-node constraints:
702 *
703 * Now redistribute the mass defect introduced by constrained degrees
704 * of freedom. This mostly affects hanging nodes neighboring a
705 * coarsened cell. Here, cell-wise mass projection might lead to a
706 * value for a constrained degree of freedom that differs from the
707 * algebraic relationship expressed by our affine constraints. Thus, we
708 * first compute the defect and then we redistribute it to all degrees
709 * of freedom on the constraint line.
710 * -----------------------------------------------------------------------
711 */
712
713 auto &new_U = std::get<0>(new_state_vector);
714
715 /*
716 * A small lambda that takes the weighted average of all degrees of
717 * freedom, and stores the result in new_U:
718 */
719 const auto update_new_state_vector = [&]() {
720 for (unsigned int local_i = 0; local_i < n_locally_owned; ++local_i) {
721
722 const auto mU_i = projected_state.get_tensor(local_i);
723 const auto m_i = projected_mass.local_element(local_i);
724
725#ifdef DEBUG_EXPENSIVE_BOUNDS_CHECK
726 const auto view = hyperbolic_system_->template view<dim, Number>();
727 AssertThrow(
728 view.is_admissible(mU_i / m_i),
729 dealii::ExcMessage("Error: inadmissible state encountered in "
730 "update_new_state_vector()"));
731#endif
732
733 new_U.write_tensor(mU_i / m_i, local_i);
734 }
735 new_U.update_ghost_values();
736 };
737
738 update_new_state_vector();
739
740 /* precomputed needs to be valid for bounds computation */
741
742 const auto &precomputed = std::get<1>(new_state_vector);
743
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;
753
754 new_U.update_ghost_values();
755
756 if constexpr (n_precomputation_cycles != 0) {
757 for (unsigned int cycle = 0; cycle < n_precomputation_cycles; ++cycle) {
758
759 SynchronizationDispatch synchronization_dispatch([&]() {
760 precomputed.update_ghost_values_start(channel++);
761 precomputed.update_ghost_values_finish();
762 });
763
765
766 auto loop = [&](auto sentinel,
767 unsigned int left,
768 unsigned int right) {
769 using T = decltype(sentinel);
770
771 /* Stored thread locally: */
772 bool thread_ready = false;
773
774 const auto view = hyperbolic_system_->template view<dim, T>();
775 view.precomputation_loop(
776 cycle,
777 [&](const unsigned int i) {
778 synchronization_dispatch.check(
779 thread_ready, i >= n_export_indices && i < n_internal);
780 },
781 sparsity_simd,
782 new_state_vector,
783 left,
784 right,
785 /*skip_constrained_dofs*/ false);
786 };
787
788 /* Parallel non-vectorized loop: */
789 loop(Number(), n_internal, n_owned);
790 /* Parallel vectorized SIMD loop: */
791 loop(VA(), 0, n_internal);
792
794 }
795 }
796 };
797
798 update_precomputed_values();
799
800 using Limiter = typename Description::template Limiter<dim, Number>;
801 const Limiter limiter(
802 *hyperbolic_system_, limiter_parameters_, precomputed);
803
804 /*
805 * Step 1: compute low-order update P_ij matrix, and bounds:
806 *
807 * We compute limiter bounds as a single value over the constraint
808 * line. This makes sense as we need to limit the update for each
809 * affected (unconstrained) degree of freedom of a constraint line with
810 * a single limiter value anyway to ensure mass conservation.
811 * (Incidentally, this avoids having to update a global, distributed
812 * bounds vector over all MPI ranks.)
813 */
814
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);
818
819 /* Only operate on a locally owned, constrained degree of freedom: */
820 if (local_i >= n_locally_owned)
821 continue;
822
823 /* The result of the mass projection: */
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;
826
827 auto &bounds = bounds_map[local_i]; /* by reference */
828 bounds = limiter.projection_bounds_from_state(local_i, U_i_star);
829
830 /* The value obtained from the affine constraints object: */
831 state_type U_i_interp;
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);
835 }
836
837 /* And redistribute low order update: */
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);
841
842 const auto bounds_k =
843 limiter.projection_bounds_from_state(local_k, U_k);
844 bounds = limiter.combine_bounds(bounds, bounds_k);
845
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;
848
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);
851 }
852 }
853
854 /* Compress vectors, recalculate unconstrained states: */
855 projected_mass.compress(dealii::VectorOperation::add);
856 projected_state.compress(dealii::VectorOperation::add);
857 kappa.compress(dealii::VectorOperation::add);
858 update_new_state_vector();
859
860 /* Redistribute ghost layer for masses and kappa: */
861 projected_mass.update_ghost_values();
862 kappa.update_ghost_values();
863
864 /* Step 2: Apply limiter: */
865
866 const auto n_iterations = limiter_parameters_.iterations();
867 for (unsigned int pass = 0; pass < n_iterations; ++pass) {
868
869 /* Update precomputed values for bounds correction: */
870 update_precomputed_values();
871
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);
875
876 /* Only operate on a locally owned, constrained degree of freedom: */
877 if (local_i >= n_locally_owned)
878 continue;
879
880 /*
881 * We are computing bounds only over a local constraint line
882 * without recombining such bounds per (unconstrained) degree of
883 * freedom globally. We avoid doing the latter because it would
884 * require a custom "VectorOperation" invoking
885 * Limiter::combine_bounds(), which we currently do not have at our
886 * disposal.
887 *
888 * As a simple workaround we simply recompute bounds for the
889 * constraint line after the low-order update and each limiter pass
890 * and recombine those into the stored value.
891 */
892
893 auto &bounds = bounds_map[local_i]; /* by reference */
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 =
899 limiter.projection_bounds_from_state(local_k, U_k);
900 bounds = limiter.combine_bounds(bounds, bounds_k);
901
902 const auto m_k = projected_mass.local_element(local_k);
903 total_mass += m_k;
904 }
905
906 auto l = Number(1.);
907
908 /* Apply relaxation: */
909 const auto relaxed_bounds =
910 limiter.fully_relax_bounds(bounds, total_mass);
911
912 /* Compute limiter values: */
913
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;
920
921 const auto &[l_k, check] = limiter.limit(relaxed_bounds, U_k, P_ik);
922#ifdef DEBUG_EXPENSIVE_BOUNDS_CHECK
923 AssertThrow(check,
924 dealii::ExcMessage("Error: low-order state out of bounds "
925 "in project / state redistribution"));
926#endif
927 l = std::min(l, l_k);
928 }
929
930 /* Apply limiter values: */
931
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);
936 mP_ik -= l * mP_ik;
937 }
938 }
939
940 /* Compress state vector, recalculate unconstrained states: */
941 projected_state.compress(dealii::VectorOperation::add);
942 update_new_state_vector();
943 }
944
945 /* Zero out constrained degrees of freedom: */
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))
949 new_U.write_tensor(state_type{}, local_i);
950 }
951 new_U.update_ghost_values();
952
953#ifdef DEBUG_SYMMETRY_CHECK
954 /*
955 * Sanity check: Final masses must agree:
956 */
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))
961 continue;
962
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,
966 dealii::ExcMessage(
967 "SolutionTransfer::projection(): something went wrong. Final "
968 "masses do not agree with those computed in OfflineData."));
969 }
970#endif
971#endif
972 }
973} // namespace ryujin
Bounds fully_relax_bounds(const Bounds &bounds, const Number &hd) const
Definition: limiter.h:235
Bounds combine_bounds(const Bounds &bounds_left, const Bounds &bounds_right) const
Definition: limiter.h:223
Bounds projection_bounds_from_state(const unsigned int i, const state_type &U_i) const
Definition: limiter.h:213
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 &parabolic_system, const std::string &subsection="/SolutionTransfer")
Vectors::ScalarVector< Number > ScalarVector
Description::template Limiter< dim, double >::Bounds Bounds
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
Definition: openmp.h:54
#define RYUJIN_PARALLEL_REGION_END
Definition: openmp.h:63
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)