ryujin 2.5.0 revision c64172211cf21344cb9b9b0f3a1266aa6400d392
Loading...
Searching...
No Matches
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 - 2026 by the ryujin authors
4//
5
6#pragma once
7
8#include <ryujin/base/computing_timer.h>
9#include <ryujin/discretization/tensor_product_point_kernels.h>
10#include <ryujin/typed/solution_transfer.h>
11
12#include <deal.II/base/exceptions.h>
13#include <deal.II/distributed/tria.h>
14#include <deal.II/dofs/dof_accessor.h>
15#include <deal.II/dofs/dof_tools.h>
16#include <deal.II/grid/cell_status.h>
17#include <deal.II/grid/tria_accessor.h>
18#include <deal.II/grid/tria_iterator.h>
19#include <deal.II/lac/block_vector.h>
20#include <deal.II/lac/la_parallel_block_vector.h>
21#include <deal.II/lac/la_parallel_vector.h>
22#include <deal.II/lac/vector.h>
23#include <deal.II/matrix_free/fe_point_evaluation.h>
24
25namespace ryujin
26{
27 template <typename HyperbolicDescription, int dim, typename Number>
28 SolutionTransfer<HyperbolicDescription, dim, Number>::SolutionTransfer(
29 const MPIEnsemble &mpi_ensemble,
30 Discretization<dim> &discretization,
31 const OfflineData<dim, Number> &offline_data,
32 const HyperbolicSystem &hyperbolic_system,
33 const Interface::ParabolicSystem &parabolic_system,
34 const std::string &subsection /* = "/SolutionTransfer" */)
35 : Interface::SolutionTransfer<dim, Number>(subsection)
36 , limiter_(hyperbolic_system, subsection + "/mass transfer limiter")
37 , mpi_ensemble_(mpi_ensemble)
38 , discretization_(&discretization)
39 , offline_data_(&offline_data)
40 , hyperbolic_system_(&hyperbolic_system)
41 , parabolic_system_(&parabolic_system)
42 {
43 }
44
45
46 namespace
47 {
52 template <typename state_type, typename Number>
53 std::vector<char>
54 pack_state_values(const std::vector<state_type> &state_values,
55 const std::vector<Number> &dynamic_values)
56 {
57 const auto n_state_bytes = sizeof(state_type) * state_values.size();
58 const auto n_dynamic_bytes = sizeof(Number) * dynamic_values.size();
59 std::vector<char> buffer(n_state_bytes + n_dynamic_bytes);
60 std::memcpy(buffer.data(), state_values.data(), n_state_bytes);
61
62 if (n_dynamic_bytes > 0)
63 std::memcpy(buffer.data() + n_state_bytes,
64 dynamic_values.data(),
65 n_dynamic_bytes);
66
67 return buffer;
68 }
69
70
76 template <typename state_type, typename Number>
77 std::pair<std::vector<state_type>, std::vector<Number>> unpack_state_values(
78 const boost::iterator_range<std::vector<char>::const_iterator>
79 &data_range,
80 const unsigned int n_state_values)
81 {
82 const std::size_t n_bytes = data_range.size();
83 const std::size_t n_state_bytes = sizeof(state_type) * n_state_values;
84 Assert(n_bytes >= n_state_bytes, dealii::ExcInternalError());
85 Assert((n_bytes - n_state_bytes) % sizeof(Number) == 0,
86 dealii::ExcInternalError());
87
88 std::vector<state_type> state_values(n_state_values);
89 std::memcpy(state_values.data(), &data_range[0], n_state_bytes);
90
91 std::vector<Number> dynamic_values((n_bytes - n_state_bytes) /
92 sizeof(Number));
93
94 if (!dynamic_values.empty())
95 std::memcpy(dynamic_values.data(),
96 &data_range[0] + n_state_bytes,
97 dynamic_values.size() * sizeof(Number));
98
99 return {std::move(state_values), std::move(dynamic_values)};
100 }
101 } // namespace
102
103
104 template <typename HyperbolicDescription, int dim, typename Number>
105 void SolutionTransfer<HyperbolicDescription, dim, Number>::prepare_projection(
106 const StateVector<Number> &old_state_vector [[maybe_unused]])
107 {
108#ifdef DEBUG_OUTPUT
109 std::cout << "SolutionTransfer<HyperbolicDescription, dim, "
110 "Number>::prepare_projection()"
111 << std::endl;
112#endif
113
114 copy_to_memory_space<dealii::MemorySpace::Host>(
115 old_state_vector,
116 StateVectorPart::state,
117 StateVectorPart::precomputed,
118 StateVectorPart::dynamic_state);
119
120 auto &triangulation = discretization_->triangulation();
121
122 Assert(this->handle_ == dealii::numbers::invalid_unsigned_int,
123 dealii::ExcMessage(
124 "You can only add one solution per SolutionTransfer object."));
125
126 for (const auto &U_dynamic_k : old_state_vector.dynamic_state)
127 U_dynamic_k.deal_ii_vector().update_ghost_values();
128
129 /*
130 * -----------------------------------------------------------------------
131 * Cell-level projection to parent cells and packing data:
132 * -----------------------------------------------------------------------
133 */
134
135 this->handle_ = triangulation.register_data_attach(
136 [this, &old_state_vector](const auto cell,
137 const dealii::CellStatus status) {
138 const auto &dof_handler = offline_data_->dof_handler();
139 const auto dof_cell = typename dealii::DoFHandler<dim>::cell_iterator(
140 &cell->get_triangulation(),
141 cell->level(),
142 cell->index(),
143 &dof_handler);
144
145 const auto &scalar_partitioner = offline_data_->scalar_partitioner();
146
147 const auto &U = old_state_vector.state;
148 const auto &U_dynamic = old_state_vector.dynamic_state;
149 /* precomputed needs to be valid for bounds computation */
150 const auto precomputed_view =
151 old_state_vector.precomputed
152 .template view<n_precomputed_values>();
153
154 const auto limiter_view = limiter_.template view<dim, Number>();
155
156 /*
157 * Collect state values for packing:
158 */
159
160 const auto n_dofs_per_cell = dof_cell->get_fe().n_dofs_per_cell();
161 std::vector<state_type> state_values(n_dofs_per_cell);
162
163 switch (status) {
164 case dealii::CellStatus::cell_will_persist:
165 [[fallthrough]];
166 case dealii::CellStatus::cell_will_be_refined: {
167 /*
168 * For both cases we need state values from the currently
169 * active cell:
170 */
171
172 Assert(dof_cell->is_active(), dealii::ExcInternalError());
173 std::vector<dealii::types::global_dof_index> dof_indices(
174 n_dofs_per_cell);
175 dof_cell->get_dof_indices(dof_indices);
176
177 std::transform(
178 std::begin(dof_indices),
179 std::end(dof_indices),
180 std::begin(state_values),
181 [&](const auto global_i) { return read_tensor(U, global_i); });
182 } break;
183
184 case dealii::CellStatus::children_will_be_coarsened: {
185 /*
186 * We need to project values from the active child cells up to
187 * the present parent cell that will become active after
188 * coarsening.
189 */
190
191 Assert(dof_cell->has_children(), dealii::ExcInternalError());
192
193 const auto &discretization = offline_data_->discretization();
194 const auto index = dof_cell->active_fe_index();
195 const auto &finite_element = discretization.finite_element()[index];
196 const auto &mapping = discretization.mapping()[index];
197 const auto &quadrature = discretization.quadrature()[index];
198
199 dealii::FEValues<dim> fe_values(
200 mapping,
201 finite_element,
202 quadrature,
203 dealii::update_values | dealii::update_JxW_values |
204 dealii::update_quadrature_points);
205
206 const auto polynomial_space =
207 dealii::internal::FEPointEvaluation::get_polynomial_space(
208 finite_element);
209
210 std::vector<dealii::Point<dim, Number>> unit_points(
211 quadrature.size());
212 /*
213 * for Number == float we need a temporary vector for the
214 * transform_points_real_to_unit_cell() function:
215 */
216 std::vector<dealii::Point<dim>> unit_points_temp(
217 std::is_same_v<Number, float> ? quadrature.size() : 0);
218
219 /* Step 1: build up right hand side by iterating over children: */
220
221 std::vector<state_type> state_values_quad(quadrature.size());
222 std::vector<state_type> local_rhs(n_dofs_per_cell);
223
224 std::vector<dealii::types::global_dof_index> dof_indices(
225 n_dofs_per_cell);
226
227 Bounds bounds;
228
229 for (unsigned int child = 0; child < dof_cell->n_children();
230 ++child) {
231 const auto child_cell = dof_cell->child(child);
232
233 Assert(child_cell->is_active(), dealii::ExcInternalError());
234 Assert(dof_cell->active_fe_index() ==
235 child_cell->active_fe_index(),
236 dealii::ExcMessage("SolutionTransfer: projection between "
237 "different FE space is not set up."));
238
239 fe_values.reinit(child_cell);
240
241 if constexpr (std::is_same_v<Number, float>) {
242 mapping.transform_points_real_to_unit_cell(
243 dof_cell,
244 fe_values.get_quadrature_points(),
245 unit_points_temp);
246 std::transform(std::begin(unit_points_temp),
247 std::end(unit_points_temp),
248 std::begin(unit_points),
249 [](const auto &x) { return x; });
250 } else {
251 mapping.transform_points_real_to_unit_cell(
252 dof_cell, fe_values.get_quadrature_points(), unit_points);
253 }
254
255 child_cell->get_dof_indices(dof_indices);
256
257 /* We want a "left fold first" for the bounds: */
258 if (child == 0 &&
259 std::begin(dof_indices) != std::end(dof_indices)) {
260 const auto global_i = dof_indices[0];
261 const auto U_i = read_tensor(U, global_i);
262 const auto local_i =
263 scalar_partitioner->global_to_local(global_i);
264 bounds = limiter_view.projection_bounds_from_state(
265 precomputed_view, local_i, U_i);
266 }
267
268 for (auto &it : state_values_quad)
269 it = state_type{};
270
271 for (unsigned int i = 0; i < n_dofs_per_cell; ++i) {
272 const auto global_i = dof_indices[i];
273 const auto U_i = read_tensor(U, global_i);
274 const auto local_i =
275 scalar_partitioner->global_to_local(global_i);
276 const auto bounds_i = limiter_view.projection_bounds_from_state(
277 precomputed_view, local_i, U_i);
278 bounds = limiter_view.combine_bounds(bounds, bounds_i);
279
280 for (unsigned int q = 0; q < quadrature.size(); ++q) {
281 state_values_quad[q] += U_i * fe_values.shape_value(i, q);
282 }
283 }
284
285 for (unsigned int q = 0; q < quadrature.size(); ++q)
286 state_values_quad[q] *= fe_values.JxW(q);
287
288 for (unsigned int q = 0; q < quadrature.size(); ++q) {
289 const unsigned int n_shapes = polynomial_space.size();
290 AssertIndexRange(n_shapes, 10);
291 dealii::ndarray<Number, 10, 2, dim> shapes;
292 // Evaluate 1d polynomials and their derivatives
293 std::array<Number, dim> point;
294 for (unsigned int d = 0; d < dim; ++d)
295 point[d] = unit_points[q][d];
296 for (unsigned int i = 0; i < n_shapes; ++i)
297 polynomial_space[i].values_of_array(point, 1, &shapes[i][0]);
298
299 Assert(finite_element.degree == 1, dealii::ExcNotImplemented());
300
301 ryujin::internal::integrate_tensor_product_value<
302 /*is linear*/ true,
303 dim,
304 Number,
305 state_type>(shapes.data(),
306 n_shapes,
307 state_values_quad[q],
308 local_rhs.data(),
309 unit_points[q],
310 true);
311 }
312 }
313
314 /* Step 2: construct inverse mass matrices: */
315
316 fe_values.reinit(dof_cell);
317
318 dealii::FullMatrix<double> mass_inverse(n_dofs_per_cell,
319 n_dofs_per_cell);
320 dealii::Vector<double> lumped_mass(n_dofs_per_cell);
321 dealii::Vector<double> lumped_mass_inverse(n_dofs_per_cell);
322
323 auto total_mass = Number(0.);
324 for (unsigned int i = 0; i < n_dofs_per_cell; ++i) {
325 for (unsigned int j = 0; j < n_dofs_per_cell; ++j) {
326 double sum = 0;
327 for (unsigned int q = 0; q < quadrature.size(); ++q)
328 sum += fe_values.shape_value(i, q) *
329 fe_values.shape_value(j, q) * fe_values.JxW(q);
330 mass_inverse(i, j) = sum;
331 lumped_mass(i) += sum;
332 }
333 lumped_mass_inverse(i) = Number(1.) / lumped_mass(i);
334 total_mass += lumped_mass(i);
335 }
336 mass_inverse.gauss_jordan();
337
338 /* Step 3: compute low-order update and P_ij matrix: */
339
340 bounds = limiter_view.fully_relax_bounds(bounds, total_mass);
341
342 std::vector<state_type> pij_matrix(n_dofs_per_cell *
343 n_dofs_per_cell);
344 dealii::FullMatrix<Number> lij_matrix(n_dofs_per_cell,
345 n_dofs_per_cell);
346
347 const auto kappa_inverse = Number(n_dofs_per_cell);
348 const auto kappa = Number(1.) / kappa_inverse;
349
350 for (unsigned int i = 0; i < n_dofs_per_cell; ++i) {
351 const state_type U_i = lumped_mass_inverse(i) * local_rhs[i];
352 state_values[i] = U_i;
353
354 for (unsigned int j = 0; j < n_dofs_per_cell; ++j) {
355 const auto kronecker_ij = Number(i == j ? 1. : 0.);
356 const auto b_ij =
357 lumped_mass(i) * mass_inverse(i, j) - kronecker_ij;
358 const auto b_ji =
359 lumped_mass(j) * mass_inverse(i, j) - kronecker_ij;
360 const auto P_ij = kappa_inverse * lumped_mass_inverse(i) *
361 (b_ij * local_rhs[j] - b_ji * local_rhs[i]);
362 pij_matrix[n_dofs_per_cell * i + j] = P_ij;
363 }
364 }
365
366 /* Step 4: compute l_ij matrix and apply limited update: */
367
368 const auto n_iterations = limiter_view.iterations();
369 for (unsigned int pass = 0; pass < n_iterations; ++pass) {
370
371 for (unsigned int i = 0; i < n_dofs_per_cell; ++i) {
372 const auto &U_i = state_values[i];
373
374 for (unsigned int j = 0; j < n_dofs_per_cell; ++j) {
375 const auto &P_ij = pij_matrix[n_dofs_per_cell * i + j];
376 const auto &[l_ij, check] =
377 limiter_view.limit(bounds, U_i, P_ij);
378 lij_matrix(i, j) = l_ij;
379 }
380 }
381
382 for (unsigned int i = 0; i < n_dofs_per_cell; ++i) {
383 auto &U_i = state_values[i];
384
385 for (unsigned int j = 0; j < n_dofs_per_cell; ++j) {
386 const auto l_ij =
387 std::min(lij_matrix(i, j), lij_matrix(j, i));
388 auto &P_ij = pij_matrix[n_dofs_per_cell * i + j];
389 U_i += kappa * l_ij * P_ij;
390 P_ij -= l_ij * P_ij;
391 }
392
393#ifdef DEBUG_EXPENSIVE_BOUNDS_CHECK
394 const auto view =
395 hyperbolic_system_->template view<dim, Number>();
396 AssertThrow(
397 view.is_admissible(U_i),
398 dealii::ExcMessage(
399 "Error: inadmissible state encountered in "
400 "register_data_attach / children_will_be_coarsened"));
401#endif
402 }
403 }
404 } break;
405
406 case dealii::CellStatus::cell_invalid:
407 Assert(false, dealii::ExcInternalError());
408 __builtin_trap();
409 break;
410 }
411
412 /*
413 * Collect dynamic state values for packing. We perform a plain
414 * Lagrange interpolation (without limiting):
415 */
416
417 std::vector<Number> dynamic_values;
418
419 if (!U_dynamic.empty()) {
420 const auto dof_cell_cg =
421 typename dealii::DoFHandler<dim>::cell_iterator(
422 &cell->get_triangulation(),
423 cell->level(),
424 cell->index(),
425 &offline_data_->dof_handler_cg());
426
427 const auto n_dofs_per_cell_cg =
428 dof_cell_cg->get_fe().n_dofs_per_cell();
429 dealii::Vector<Number> local_values(n_dofs_per_cell_cg);
430
431 for (const auto &U_k : U_dynamic) {
432 dof_cell_cg->get_interpolated_dof_values(U_k.deal_ii_vector(),
433 local_values);
434 dynamic_values.insert(std::end(dynamic_values),
435 std::begin(local_values),
436 std::end(local_values));
437 }
438 }
439
440 return pack_state_values(state_values, dynamic_values);
441 },
442 /* returns_variable_size_data =*/false);
443 }
444
445
446 template <typename HyperbolicDescription, int dim, typename Number>
447 void SolutionTransfer<HyperbolicDescription, dim, Number>::project(
448 StateVector<Number> &new_state_vector [[maybe_unused]])
449 {
450#ifdef DEBUG_OUTPUT
451 std::cout
452 << "SolutionTransfer<HyperbolicDescription, dim, Number>::project()"
453 << std::endl;
454#endif
455
456 move_to_memory_space<dealii::MemorySpace::Host>(
457 new_state_vector,
458 StateVectorPart::state,
459 StateVectorPart::precomputed,
460 StateVectorPart::dynamic_state);
461
462 auto &triangulation = discretization_->triangulation();
463
464 Assert(
465 this->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 ScalarHostVector<Number> projected_mass;
475 projected_mass.reinit(offline_data_->scalar_partitioner());
476 MultiComponentVector<Number> projected_state;
477 projected_state.reinit_with_vector_partitioner(
478 offline_data_->hyperbolic_vector_partitioner(), problem_dimension);
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 ScalarHostVector<Number> 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 /*
500 * The interpolation of the dynamic state writes into the ghost range.
501 * Make sure that the vectors are not in ghosted state:
502 */
503
504 auto &U_dynamic = new_state_vector.dynamic_state;
505 for (auto &U_dynamic_k : U_dynamic)
506 U_dynamic_k.deal_ii_vector().zero_out_ghost_values();
507
508 triangulation.notify_ready_to_unpack( //
509 this->handle_,
510 [this, &projected_mass, &projected_state, &U_dynamic](
511 const auto &cell,
512 const dealii::CellStatus status,
513 const auto &data_range) {
514 const auto &dof_handler = offline_data_->dof_handler();
515 const auto dof_cell = typename dealii::DoFHandler<dim>::cell_iterator(
516 &cell->get_triangulation(),
517 cell->level(),
518 cell->index(),
519 &dof_handler);
520
521 /*
522 * Retrieve packed values and project onto cell:
523 */
524
525 const auto n_dofs_per_cell = dof_cell->get_fe().n_dofs_per_cell();
526 std::vector<dealii::types::global_dof_index> dof_indices(
527 n_dofs_per_cell);
528
529 const auto [state_values, dynamic_values] =
530 unpack_state_values<state_type, Number>(data_range,
531 n_dofs_per_cell);
532
533 Assert(dynamic_values.empty() == U_dynamic.empty(),
534 dealii::ExcMessage("The attached data does not match the "
535 "dynamic state of the state vector."));
536
537 /*
538 * Interpolate the dynamic state: set_dof_values_by_interpolation()
539 * sets the values of an active cell, or prolongs the values onto
540 * the children of a cell that has been refined.
541 */
542
543 if (!U_dynamic.empty()) {
544 const auto dof_cell_cg =
545 typename dealii::DoFHandler<dim>::cell_iterator(
546 &cell->get_triangulation(),
547 cell->level(),
548 cell->index(),
549 &offline_data_->dof_handler_cg());
550
551 const auto n_dofs_per_cell_cg =
552 dof_cell_cg->get_fe().n_dofs_per_cell();
553
554 Assert(dynamic_values.size() ==
555 U_dynamic.size() * n_dofs_per_cell_cg,
556 dealii::ExcMessage(
557 "The attached data has " +
558 std::to_string(dynamic_values.size()) +
559 " dynamic state values for this cell, but " +
560 std::to_string(U_dynamic.size()) +
561 " dynamic state components with " +
562 std::to_string(n_dofs_per_cell_cg) +
563 " degrees of freedom per cell require " +
564 std::to_string(U_dynamic.size() * n_dofs_per_cell_cg) +
565 "."));
566
567 dealii::Vector<Number> local_values(n_dofs_per_cell_cg);
568 for (std::size_t k = 0; k < U_dynamic.size(); ++k) {
569 std::copy_n(std::begin(dynamic_values) + k * n_dofs_per_cell_cg,
570 n_dofs_per_cell_cg,
571 std::begin(local_values));
572 dof_cell_cg->set_dof_values_by_interpolation(
573 local_values, U_dynamic[k].deal_ii_vector());
574 }
575 }
576
577 switch (status) {
578 case dealii::CellStatus::cell_will_persist:
579 [[fallthrough]];
580 case dealii::CellStatus::children_will_be_coarsened: {
581 /*
582 * For both cases we distribute stored state_values to the
583 * projected_state and projected_mass vectors.
584 */
585
586 Assert(dof_cell->is_active(), dealii::ExcInternalError());
587 dof_cell->get_dof_indices(dof_indices);
588
589 const auto &discretization = offline_data_->discretization();
590 const auto index = dof_cell->active_fe_index();
591 const auto &finite_element = discretization.finite_element()[index];
592 const auto &mapping = discretization.mapping()[index];
593 const auto &quadrature = discretization.quadrature()[index];
594
595 dealii::FEValues<dim> fe_values(mapping,
596 finite_element,
597 quadrature,
598 dealii::update_values |
599 dealii::update_JxW_values);
600
601 fe_values.reinit(dof_cell);
602
603 dealii::Vector<double> mi(n_dofs_per_cell);
604 for (unsigned int i = 0; i < n_dofs_per_cell; ++i) {
605 double sum = 0;
606 for (unsigned int q = 0; q < quadrature.size(); ++q)
607 sum += fe_values.shape_value(i, q) * fe_values.JxW(q);
608 mi(i) += sum;
609 }
610
611 for (unsigned int i = 0; i < n_dofs_per_cell; ++i) {
612 const auto global_i = dof_indices[i];
613 add_tensor(projected_state, mi(i) * state_values[i], global_i);
614 projected_mass(global_i) += mi(i);
615 }
616
617 } break;
618
619 case dealii::CellStatus::cell_will_be_refined: {
620 /*
621 * We are on a (non active) cell that has been refined. Project
622 * onto the children and do a local mass projection there:
623 */
624
625 Assert(dof_cell->has_children(), dealii::ExcInternalError());
626
627 const auto &discretization = offline_data_->discretization();
628 const auto index = dof_cell->active_fe_index();
629 const auto &finite_element = discretization.finite_element()[index];
630 const auto &mapping = discretization.mapping()[index];
631 const auto &quadrature = discretization.quadrature()[index];
632
633 dealii::FEValues<dim> fe_values(
634 mapping,
635 finite_element,
636 quadrature,
637 dealii::update_values | dealii::update_JxW_values |
638 dealii::update_quadrature_points);
639
640 const auto polynomial_space =
641 dealii::internal::FEPointEvaluation::get_polynomial_space(
642 finite_element);
643 std::vector<dealii::Point<dim, Number>> unit_points(
644 quadrature.size());
645 /*
646 * for Number == float we need a temporary vector for the
647 * transform_points_real_to_unit_cell() function:
648 */
649 std::vector<dealii::Point<dim>> unit_points_temp(
650 std::is_same_v<Number, float> ? quadrature.size() : 0);
651
652 dealii::FullMatrix<double> mass_inverse(n_dofs_per_cell,
653 n_dofs_per_cell);
654 dealii::Vector<double> lumped_mass(n_dofs_per_cell);
655 std::vector<state_type> local_rhs(n_dofs_per_cell);
656
657 for (unsigned int child = 0; child < dof_cell->n_children();
658 ++child) {
659 const auto child_cell = dof_cell->child(child);
660
661 Assert(child_cell->is_active(), dealii::ExcInternalError());
662 Assert(dof_cell->active_fe_index() ==
663 child_cell->active_fe_index(),
664 dealii::ExcMessage("SolutionTransfer: projection between "
665 "different FE space is not set up."));
666
667 child_cell->get_dof_indices(dof_indices);
668
669 /* Step 1: build up right hand side on child cell: */
670
671 fe_values.reinit(child_cell);
672
673 if constexpr (std::is_same_v<Number, float>) {
674 mapping.transform_points_real_to_unit_cell(
675 dof_cell,
676 fe_values.get_quadrature_points(),
677 unit_points_temp);
678 std::transform(std::begin(unit_points_temp),
679 std::end(unit_points_temp),
680 std::begin(unit_points),
681 [](const auto &x) { return x; });
682 } else {
683 mapping.transform_points_real_to_unit_cell(
684 dof_cell, fe_values.get_quadrature_points(), unit_points);
685 }
686
687 for (auto &it : local_rhs)
688 it = state_type{};
689
690 for (unsigned int q = 0; q < quadrature.size(); ++q) {
691 Assert(finite_element.degree == 1, dealii::ExcNotImplemented());
692 auto coefficient =
693 dealii::internal::evaluate_tensor_product_value(
694 polynomial_space,
695 make_const_array_view(state_values),
696 unit_points[q],
697 /*is linear*/ true);
698 coefficient *= fe_values.JxW(q);
699
700 for (unsigned int i = 0; i < n_dofs_per_cell; ++i)
701 local_rhs[i] += coefficient * fe_values.shape_value(i, q);
702 }
703
704 /* Step 2: solve with inverse mass matrix on child cell: */
705
706 mass_inverse = Number(0.);
707 lumped_mass = Number(0.);
708 for (unsigned int i = 0; i < n_dofs_per_cell; ++i) {
709 for (unsigned int j = 0; j < n_dofs_per_cell; ++j) {
710 double sum = 0;
711 for (unsigned int q = 0; q < quadrature.size(); ++q)
712 sum += fe_values.shape_value(i, q) *
713 fe_values.shape_value(j, q) * fe_values.JxW(q);
714 mass_inverse(i, j) = sum;
715 lumped_mass(i) += sum;
716 }
717 }
718 mass_inverse.gauss_jordan();
719
720 /* Step 3: compute high order update and write back: */
721
722 for (unsigned int i = 0; i < n_dofs_per_cell; ++i) {
723 state_type U_i;
724 for (unsigned int j = 0; j < n_dofs_per_cell; ++j) {
725 U_i += mass_inverse(i, j) * local_rhs[j];
726 }
727
728#ifdef DEBUG_EXPENSIVE_BOUNDS_CHECK
729 const auto view =
730 hyperbolic_system_->template view<dim, Number>();
731 AssertThrow(view.is_admissible(U_i),
732 dealii::ExcMessage(
733 "Error: inadmissible state encountered in "
734 "ready_to_unpack / cell_will_be_refined"));
735#endif
736 const auto global_i = dof_indices[i];
737 add_tensor(projected_state, lumped_mass(i) * U_i, global_i);
738 projected_mass(global_i) += lumped_mass(i);
739 }
740 } /*child*/
741
742 } break;
743
744 case dealii::CellStatus::cell_invalid:
745 Assert(false, dealii::ExcInternalError());
746 __builtin_trap();
747 break;
748 }
749 });
750
751 /*
752 * Discard the values written into the ghost range, distribute the
753 * interpolated dynamic state to constrained degrees of freedom and
754 * update the ghost range:
755 */
756 for (auto &U_dynamic_k : U_dynamic) {
757 auto &vector = U_dynamic_k.deal_ii_vector();
758 vector.compress(dealii::VectorOperation::insert);
759 offline_data_->affine_constraints_cg().distribute(vector);
760 vector.update_ghost_values();
761 }
762
763 const auto projected_state_view =
764 projected_state.template view<problem_dimension>();
765
766 projected_mass.compress(dealii::VectorOperation::add);
767 projected_state_view.compress(dealii::VectorOperation::add);
768
769 /*
770 * -----------------------------------------------------------------------
771 * Redistribute masses to satisfy hanging-node constraints:
772 *
773 * Now redistribute the mass defect introduced by constrained degrees
774 * of freedom. This mostly affects hanging nodes neighboring a
775 * coarsened cell. Here, cell-wise mass projection might lead to a
776 * value for a constrained degree of freedom that differs from the
777 * algebraic relationship expressed by our affine constraints. Thus, we
778 * first compute the defect and then we redistribute it to all degrees
779 * of freedom on the constraint line.
780 * -----------------------------------------------------------------------
781 */
782
783 const auto new_U_view =
784 new_state_vector.state.template view<problem_dimension>();
785
786 /*
787 * A small lambda that takes the weighted average of all degrees of
788 * freedom, and stores the result in new_U:
789 */
790 const auto update_new_state_vector = [&]() {
791 for (unsigned int local_i = 0; local_i < n_locally_owned; ++local_i) {
792
793 const auto mU_i = projected_state_view.read_tensor(local_i);
794 const auto m_i = projected_mass.local_element(local_i);
795
796#ifdef DEBUG_EXPENSIVE_BOUNDS_CHECK
797 const auto view = hyperbolic_system_->template view<dim, Number>();
798 AssertThrow(
799 view.is_admissible(mU_i / m_i),
800 dealii::ExcMessage("Error: inadmissible state encountered in "
801 "update_new_state_vector()"));
802#endif
803
804 new_U_view.write_tensor(mU_i / m_i, local_i);
805 }
806 new_U_view.update_ghost_values();
807 };
808
809 update_new_state_vector();
810
811 const auto precomputed_view =
812 new_state_vector.precomputed.template view<n_precomputed_values>();
813
814 /* The limiter requires valid precomputed values. Therefore, update: */
815 const auto update_precomputed_values = [&]() {
816 new_U_view.update_ghost_values();
817 hyperbolic_system_->fill_precomputed_values(
818 *offline_data_, new_state_vector, /*skip_constrainted_dofs*/ false);
819 precomputed_view.update_ghost_values();
820 };
821
822 update_precomputed_values();
823
824 const auto limiter_view = limiter_.template view<dim, Number>();
825
826 /*
827 * Step 1: compute low-order update P_ij matrix, and bounds:
828 *
829 * We compute limiter bounds as a single value over the constraint
830 * line. This makes sense as we need to limit the update for each
831 * affected (unconstrained) degree of freedom of a constraint line with
832 * a single limiter value anyway to ensure mass conservation.
833 * (Incidentally, this avoids having to update a global, distributed
834 * bounds vector over all MPI ranks.)
835 */
836
837 for (const auto &line : affine_constraints.get_lines()) {
838 const auto global_i = line.index;
839 const auto local_i = scalar_partitioner->global_to_local(global_i);
840
841 /* Only operate on a locally owned, constrained degree of freedom: */
842 if (local_i >= n_locally_owned)
843 continue;
844
845 /* The result of the mass projection: */
846 const auto m_i_star = projected_mass.local_element(local_i);
847 const auto U_i_star =
848 projected_state_view.read_tensor(local_i) / m_i_star;
849
850 auto &bounds = bounds_map[local_i]; /* by reference */
851 bounds = limiter_view.projection_bounds_from_state(
852 precomputed_view, local_i, U_i_star);
853
854 /* The value obtained from the affine constraints object: */
855 state_type U_i_interp;
856 for (const auto &[global_k, c_k] : line.entries) {
857 const auto local_k = scalar_partitioner->global_to_local(global_k);
858 U_i_interp += c_k * new_U_view.read_tensor(local_k);
859 }
860
861 /* And redistribute low order update: */
862 for (const auto &[global_k, c_k] : line.entries) {
863 const auto local_k = scalar_partitioner->global_to_local(global_k);
864 const auto U_k = new_U_view.read_tensor(local_k);
865
866 const auto bounds_k = limiter_view.projection_bounds_from_state(
867 precomputed_view, local_k, U_k);
868 bounds = limiter_view.combine_bounds(bounds, bounds_k);
869
870 projected_state_view.add_tensor(c_k * m_i_star * U_i_star, local_k);
871 projected_mass.local_element(local_k) += c_k * m_i_star;
872
873 kappa.local_element(local_k) += Number(1.);
874 pik_matrix[{local_i, local_k}] = c_k * m_i_star * (U_k - U_i_interp);
875 }
876 }
877
878 /* Compress vectors, recalculate unconstrained states: */
879 projected_mass.compress(dealii::VectorOperation::add);
880 projected_state_view.compress(dealii::VectorOperation::add);
881 kappa.compress(dealii::VectorOperation::add);
882 update_new_state_vector();
883
884 /* Redistribute ghost layer for masses and kappa: */
885 projected_mass.update_ghost_values();
886 kappa.update_ghost_values();
887
888 /* Step 2: Apply limiter: */
889
890 const auto n_iterations = limiter_view.iterations();
891 for (unsigned int pass = 0; pass < n_iterations; ++pass) {
892
893 /* Update precomputed values for bounds correction: */
894 update_precomputed_values();
895
896 for (const auto &line : affine_constraints.get_lines()) {
897 const auto global_i = line.index;
898 const auto local_i = scalar_partitioner->global_to_local(global_i);
899
900 /* Only operate on a locally owned, constrained degree of freedom: */
901 if (local_i >= n_locally_owned)
902 continue;
903
904 /*
905 * We are computing bounds only over a local constraint line
906 * without recombining such bounds per (unconstrained) degree of
907 * freedom globally. We avoid doing the latter because it would
908 * require a custom "VectorOperation" invoking
909 * LimiterView::combine_bounds(), which we currently do not have at our
910 * disposal.
911 *
912 * As a simple workaround we simply recompute bounds for the
913 * constraint line after the low-order update and each limiter pass
914 * and recombine those into the stored value.
915 */
916
917 auto &bounds = bounds_map[local_i]; /* by reference */
918 auto total_mass = Number(0.);
919 for (const auto &[global_k, c_k] : line.entries) {
920 const auto local_k = scalar_partitioner->global_to_local(global_k);
921 const auto U_k = new_U_view.read_tensor(local_k);
922 const auto bounds_k = limiter_view.projection_bounds_from_state(
923 precomputed_view, local_k, U_k);
924 bounds = limiter_view.combine_bounds(bounds, bounds_k);
925
926 const auto m_k = projected_mass.local_element(local_k);
927 total_mass += m_k;
928 }
929
930 auto l = Number(1.);
931
932 /* Apply relaxation: */
933 const auto relaxed_bounds =
934 limiter_view.fully_relax_bounds(bounds, total_mass);
935
936 /* Compute limiter values: */
937
938 for (const auto &[global_k, c_k] : line.entries) {
939 const auto local_k = scalar_partitioner->global_to_local(global_k);
940 const auto kappa_k = kappa.local_element(local_k);
941 const auto m_k = projected_mass.local_element(local_k);
942 const auto U_k = new_U_view.read_tensor(local_k);
943 const auto P_ik = pik_matrix[{local_i, local_k}] * kappa_k / m_k;
944
945 const auto &[l_k, check] =
946 limiter_view.limit(relaxed_bounds, U_k, P_ik);
947 l = std::min(l, l_k);
948 }
949
950 /* Apply limiter values: */
951
952 for (const auto &[global_k, c_k] : line.entries) {
953 const auto local_k = scalar_partitioner->global_to_local(global_k);
954 auto &mP_ik = pik_matrix[{local_i, local_k}];
955 projected_state_view.add_tensor(l * mP_ik, local_k);
956 mP_ik -= l * mP_ik;
957 }
958 }
959
960 /* Compress state vector, recalculate unconstrained states: */
961 projected_state_view.compress(dealii::VectorOperation::add);
962 update_new_state_vector();
963 }
964
965 /* Zero out constrained degrees of freedom: */
966 for (unsigned int local_i = 0; local_i < n_locally_owned; ++local_i) {
967 const auto global_i = scalar_partitioner->local_to_global(local_i);
968 if (affine_constraints.is_constrained(global_i))
969 new_U_view.write_tensor(state_type{}, local_i);
970 }
971 new_U_view.update_ghost_values();
972
973#ifdef DEBUG_SYMMETRY_CHECK
974 /*
975 * Sanity check: Final masses must agree:
976 */
977 const auto &lumped_mass_matrix = offline_data_->lumped_mass_matrix();
978 for (unsigned int local_i = 0; local_i < n_locally_owned; ++local_i) {
979 const auto global_i = scalar_partitioner->local_to_global(local_i);
980 if (affine_constraints.is_constrained(global_i))
981 continue;
982
983 const auto m_i = projected_mass.local_element(local_i);
984 const auto m_i_reference = lumped_mass_matrix.view().read_entry(local_i);
985 Assert(std::abs(m_i - m_i_reference) < 1.e-10,
986 dealii::ExcMessage(
987 "SolutionTransfer::projection(): something went wrong. Final "
988 "masses do not agree with those computed in OfflineData."));
989 }
990#endif
991 }
992
993
994 template <typename HyperbolicDescription, int dim, typename Number>
995 inline DEAL_II_ALWAYS_INLINE auto
996 SolutionTransfer<HyperbolicDescription, dim, Number>::read_tensor(
997 const MultiComponentVector<Number> &U,
998 const dealii::types::global_dof_index global_i) -> state_type
999 {
1000 const auto &scalar_partitioner = offline_data_->scalar_partitioner();
1001 const auto &affine_constraints = offline_data_->affine_constraints();
1002 const auto local_i = scalar_partitioner->global_to_local(global_i);
1003 const auto U_view = U.template view<problem_dimension>();
1004 if (affine_constraints.is_constrained(global_i)) {
1005 state_type result;
1006 const auto &line = *affine_constraints.get_constraint_entries(global_i);
1007 for (const auto &[global_k, c_k] : line) {
1008 const auto local_k = scalar_partitioner->global_to_local(global_k);
1009 result += c_k * U_view.read_tensor(local_k);
1010 }
1011 return result;
1012 } else {
1013 return U_view.read_tensor(local_i);
1014 }
1015 }
1016
1017
1018 template <typename HyperbolicDescription, int dim, typename Number>
1019 inline DEAL_II_ALWAYS_INLINE void
1020 SolutionTransfer<HyperbolicDescription, dim, Number>::add_tensor(
1021 MultiComponentVector<Number> &U,
1022 const state_type &new_U_i,
1023 const dealii::types::global_dof_index global_i)
1024 {
1025 const auto &scalar_partitioner = offline_data_->scalar_partitioner();
1026 const auto local_i = scalar_partitioner->global_to_local(global_i);
1027 U.template view<problem_dimension>().add_tensor(new_U_i, local_i);
1028 }
1029} // namespace ryujin
ryujin::MultiComponentVector::reinit_with_vector_partitioner
void reinit_with_vector_partitioner(const std::shared_ptr< const dealii::Utilities::MPI::Partitioner > &vector_partitioner, const unsigned int n_comp, const TransferPolicy transfer_policy=TransferPolicy::explicit_transfers)
ryujin::SolutionTransfer::HyperbolicSystem
typename HyperbolicDescription::HyperbolicSystem HyperbolicSystem
Definition solution_transfer.h:49
ryujin::SolutionTransfer::Bounds
typename LimiterView::Bounds Bounds
Definition solution_transfer.h:72
ryujin::SolutionTransfer::project
void project(StateVector< Number > &new_state_vector) override
Definition solution_transfer.template.h:447
ryujin::SolutionTransfer::state_type
typename View::state_type state_type
Definition solution_transfer.h:62
ryujin::SolutionTransfer::SolutionTransfer
SolutionTransfer(const MPIEnsemble &mpi_ensemble, Discretization< dim > &discretization, const OfflineData< dim, Number > &offline_data, const HyperbolicSystem &hyperbolic_system, const Interface::ParabolicSystem &parabolic_system, const std::string &subsection="/SolutionTransfer")
Definition solution_transfer.template.h:28
ryujin::SolutionTransfer::prepare_projection
void prepare_projection(const StateVector< Number > &old_state_vector) override
Definition solution_transfer.template.h:105
computing_timer.h
ryujin::ScalarHostVector
dealii::LinearAlgebra::distributed::Vector< Number > ScalarHostVector
Definition state_vector.h:33
ryujin::StateVectorPart::dynamic_state
@ dynamic_state
ryujin::StateVectorPart::precomputed
@ precomputed
ryujin::StateVectorPart::state
@ state
ryujin::internal::integrate_tensor_product_value
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)
Definition tensor_product_point_kernels.h:108
tensor_product_point_kernels.h
solution_transfer.h
Generated by doxygen 1.9.8