Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <algorithm>
7 : #include <cstddef>
8 : #include <utility>
9 :
10 : #include "DataStructures/DataVector.hpp"
11 : #include "DataStructures/Index.hpp"
12 : #include "DataStructures/SliceIterator.hpp"
13 : #include "DataStructures/Tensor/Tensor.hpp"
14 : #include "DataStructures/Variables.hpp"
15 : #include "Domain/Structure/Direction.hpp"
16 : #include "Domain/Structure/DirectionalId.hpp"
17 : #include "Domain/Structure/DirectionalIdMap.hpp"
18 : #include "Domain/Structure/Element.hpp"
19 : #include "Domain/Structure/ElementId.hpp"
20 : #include "Evolution/DgSubcell/Projection.hpp"
21 : #include "Evolution/DiscontinuousGalerkin/MortarData.hpp"
22 : #include "Evolution/DiscontinuousGalerkin/MortarDataHolder.hpp"
23 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
24 : #include "Utilities/ErrorHandling/Assert.hpp"
25 : #include "Utilities/Gsl.hpp"
26 : #include "Utilities/TMPL.hpp"
27 :
28 : namespace evolution::dg::subcell {
29 : /*!
30 : * \brief Project the DG package data to the subcells. Data received from a
31 : * neighboring element doing DG is always projected, while the data we sent to
32 : * our neighbors before doing a rollback from DG to subcell is only projected if
33 : * `OverwriteInternalMortarData` is `true`.
34 : *
35 : * In order for the hybrid DG-FD/FV scheme to be conservative between elements
36 : * using DG and elements using subcell, the boundary terms must be the same on
37 : * both elements. In practice this means the boundary corrections \f$G+D\f$ must
38 : * be computed on the same grid. Consider the element doing subcell which
39 : * receives data from an element doing DG. In this case the DG element's
40 : * ingredients going into \f$G+D\f$ are projected to the subcells and then
41 : * \f$G+D\f$ are computed on the subcells. Similarly, for strict conservation
42 : * the element doing DG must first project the data it sent to the neighbor to
43 : * the subcells, then compute \f$G+D\f$ on the subcells, and finally reconstrct
44 : * \f$G+D\f$ back to the DG grid before lifting \f$G+D\f$ to the volume.
45 : *
46 : * This function updates the `packaged_data` (ingredients into \f$G+D\f$)
47 : * received by an element doing subcell by projecting the neighbor's DG data
48 : * onto the subcells. Note that this is only half of what is required for strict
49 : * conservation, the DG element must also compute \f$G+D\f$ on the subcells.
50 : * Note that we currently do not perform the other half of the correction
51 : * needed to be strictly conservative.
52 : *
53 : * If we are retaking a time step after the DG step failed then maintaining
54 : * conservation requires additional care. If `OverwriteInternalMortarData` is
55 : * `true` then the local (the element switching from DG to subcell) ingredients
56 : * into \f$G+D\f$ are projected and overwrite the data computed from
57 : * the FD reconstruction to the interface. However, even this is insufficient to
58 : * guarantee conservation. To guarantee conservation (which we do not currently
59 : * do) the correction \f$G+D\f$ must be computed on the DG grid and then
60 : * projected to the subcells.
61 : *
62 : * Note that our practical experience shows that since the DG-subcell hybrid
63 : * scheme switches to the subcell solver _before_ the local solution contains
64 : * discontinuities, strict conservation is not necessary between DG and FD/FV
65 : * regions. This was also observed with a block-adaptive finite difference AMR
66 : * code \cite CHEN2016604
67 : *
68 : * The variable `variables_to_offset_in_dg_grid` is used in cases like the
69 : * combined generalized harmonic and GRMHD system where the DG grid uses
70 : * boundary corrections for both GH and GRMHD, but the subcell grid only does
71 : * GRMHD.
72 : */
73 : template <bool OverwriteInternalMortarData, size_t Dim,
74 : typename DgPackageFieldTags>
75 1 : void correct_package_data(
76 : const gsl::not_null<Variables<DgPackageFieldTags>*> lower_packaged_data,
77 : const gsl::not_null<Variables<DgPackageFieldTags>*> upper_packaged_data,
78 : const size_t logical_dimension_to_operate_in, const Element<Dim>& element,
79 : const Mesh<Dim>& subcell_volume_mesh,
80 : const DirectionalIdMap<Dim, evolution::dg::MortarDataHolder<Dim>>&
81 : mortar_data,
82 : const size_t variables_to_offset_in_dg_grid) {
83 : const Direction<Dim> upper_direction{logical_dimension_to_operate_in,
84 : Side::Upper};
85 : const Direction<Dim> lower_direction{logical_dimension_to_operate_in,
86 : Side::Lower};
87 : const bool has_upper_neighbor = element.neighbors().contains(upper_direction);
88 : const bool has_lower_neighbor = element.neighbors().contains(lower_direction);
89 : const DirectionalId<Dim> upper_mortar_id =
90 : has_upper_neighbor
91 : ? DirectionalId<Dim>{upper_direction,
92 : *element.neighbors().at(upper_direction).begin()}
93 : : DirectionalId<Dim>{};
94 : const DirectionalId<Dim> lower_mortar_id =
95 : has_lower_neighbor
96 : ? DirectionalId<Dim>{lower_direction,
97 : *element.neighbors().at(lower_direction).begin()}
98 : : DirectionalId<Dim>{};
99 :
100 : Index<Dim> subcell_extents_with_faces = subcell_volume_mesh.extents();
101 : ++subcell_extents_with_faces[logical_dimension_to_operate_in];
102 : const Mesh<Dim - 1>& subcell_face_mesh =
103 : subcell_volume_mesh.slice_away(logical_dimension_to_operate_in);
104 :
105 : const auto project_dg_data_to_subcells =
106 : [logical_dimension_to_operate_in, &subcell_extents_with_faces,
107 : &subcell_face_mesh, variables_to_offset_in_dg_grid](
108 : const gsl::not_null<Variables<DgPackageFieldTags>*>
109 : subcell_packaged_data,
110 : const size_t subcell_index, const Mesh<Dim - 1>& neighbor_face_mesh,
111 : const DataVector& neighbor_data) {
112 : const size_t dg_variables_offset_size =
113 : variables_to_offset_in_dg_grid *
114 : neighbor_face_mesh.number_of_grid_points();
115 : const double* slice_data =
116 : // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-pointer-arithmetic)
117 : neighbor_data.data() + dg_variables_offset_size;
118 : // Warning: projected_data can't be inside the `if constexpr` since that
119 : // would lead to a dangling pointer.
120 : DataVector projected_data{};
121 : if constexpr (Dim > 1) {
122 : projected_data = DataVector{
123 : Variables<DgPackageFieldTags>::number_of_independent_components *
124 : subcell_face_mesh.number_of_grid_points()};
125 : evolution::dg::subcell::fd::detail::project_impl_with_tag_list<
126 : DgPackageFieldTags>(
127 : gsl::make_span(projected_data.data(), projected_data.size()),
128 : gsl::make_span(
129 : slice_data,
130 : Variables<
131 : DgPackageFieldTags>::number_of_independent_components *
132 : neighbor_face_mesh.number_of_grid_points()),
133 : neighbor_face_mesh, subcell_face_mesh.extents());
134 : slice_data = projected_data.data();
135 : } else {
136 : (void)subcell_face_mesh;
137 : (void)projected_data;
138 : }
139 : const size_t volume_grid_points = subcell_extents_with_faces.product();
140 : const size_t slice_grid_points =
141 : subcell_extents_with_faces
142 : .slice_away(logical_dimension_to_operate_in)
143 : .product();
144 : double* const volume_data = subcell_packaged_data->data();
145 : for (SliceIterator si(subcell_extents_with_faces,
146 : logical_dimension_to_operate_in, subcell_index);
147 : si; ++si) {
148 : for (size_t i = 0;
149 : i <
150 : Variables<DgPackageFieldTags>::number_of_independent_components;
151 : ++i) {
152 : // clang-tidy: do not use pointer arithmetic
153 : volume_data[si.volume_offset() +
154 : i * volume_grid_points] = // NOLINT
155 : slice_data[si.slice_offset() +
156 : i * slice_grid_points]; // NOLINT
157 : }
158 : }
159 : };
160 :
161 : // Project DG data to the subcells
162 : if (auto neighbor_mortar_data_it = mortar_data.find(upper_mortar_id);
163 : has_upper_neighbor and neighbor_mortar_data_it != mortar_data.end()) {
164 : if (neighbor_mortar_data_it->second.neighbor().mortar_data.has_value()) {
165 : project_dg_data_to_subcells(
166 : upper_packaged_data,
167 : subcell_extents_with_faces[logical_dimension_to_operate_in] - 1,
168 : neighbor_mortar_data_it->second.neighbor().face_mesh.value(),
169 : neighbor_mortar_data_it->second.neighbor().mortar_data.value());
170 : }
171 : if constexpr (OverwriteInternalMortarData) {
172 : if (neighbor_mortar_data_it->second.local().mortar_data.has_value()) {
173 : project_dg_data_to_subcells(
174 : lower_packaged_data,
175 : subcell_extents_with_faces[logical_dimension_to_operate_in] - 1,
176 : neighbor_mortar_data_it->second.local().face_mesh.value(),
177 : neighbor_mortar_data_it->second.local().mortar_data.value());
178 : }
179 : }
180 : }
181 : if (auto neighbor_mortar_data_it = mortar_data.find(lower_mortar_id);
182 : has_lower_neighbor and neighbor_mortar_data_it != mortar_data.end()) {
183 : if (neighbor_mortar_data_it->second.neighbor().mortar_data.has_value()) {
184 : project_dg_data_to_subcells(
185 : lower_packaged_data, 0,
186 : neighbor_mortar_data_it->second.neighbor().face_mesh.value(),
187 : neighbor_mortar_data_it->second.neighbor().mortar_data.value());
188 : }
189 : if constexpr (OverwriteInternalMortarData) {
190 : if (neighbor_mortar_data_it->second.local().mortar_data.has_value()) {
191 : project_dg_data_to_subcells(
192 : upper_packaged_data, 0,
193 : neighbor_mortar_data_it->second.local().face_mesh.value(),
194 : neighbor_mortar_data_it->second.local().mortar_data.value());
195 : }
196 : }
197 : }
198 : }
199 : } // namespace evolution::dg::subcell
|
