Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cstddef>
7 : #include <optional>
8 : #include <type_traits>
9 : #include <utility>
10 :
11 : #include "DataStructures/DataBox/DataBox.hpp"
12 : #include "DataStructures/DataBox/PrefixHelpers.hpp"
13 : #include "DataStructures/DataBox/Prefixes.hpp"
14 : #include "DataStructures/DataVector.hpp"
15 : #include "DataStructures/Tensor/Tensor.hpp"
16 : #include "DataStructures/Variables.hpp"
17 : #include "Domain/CoordinateMaps/CoordinateMap.hpp"
18 : #include "Domain/FaceNormal.hpp"
19 : #include "Domain/Structure/Direction.hpp"
20 : #include "Domain/Structure/DirectionMap.hpp"
21 : #include "Domain/Structure/DirectionalIdMap.hpp"
22 : #include "Domain/Structure/Element.hpp"
23 : #include "Domain/Structure/ElementId.hpp"
24 : #include "Domain/Tags.hpp"
25 : #include "Domain/TagsTimeDependent.hpp"
26 : #include "Evolution/DiscontinuousGalerkin/Actions/ComputeTimeDerivativeHelpers.hpp"
27 : #include "Evolution/DiscontinuousGalerkin/Actions/NormalCovectorAndMagnitude.hpp"
28 : #include "Evolution/DiscontinuousGalerkin/Actions/PackageDataImpl.hpp"
29 : #include "Evolution/DiscontinuousGalerkin/MortarTags.hpp"
30 : #include "NumericalAlgorithms/DiscontinuousGalerkin/MortarHelpers.hpp"
31 : #include "NumericalAlgorithms/DiscontinuousGalerkin/ProjectToBoundary.hpp"
32 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
33 : #include "NumericalAlgorithms/Spectral/SegmentSize.hpp"
34 : #include "Time/TimeStepId.hpp"
35 : #include "Utilities/Gsl.hpp"
36 : #include "Utilities/TMPL.hpp"
37 :
38 : /// \cond
39 : namespace Tags {
40 : struct TimeStepId;
41 : } // namespace Tags
42 : /// \endcond
43 :
44 : namespace evolution::dg::Actions::detail {
45 : template <typename System, size_t Dim, typename BoundaryCorrection,
46 : typename TemporaryTags, typename... PackageDataVolumeArgs>
47 : void internal_mortar_data_impl(
48 : const gsl::not_null<
49 : DirectionMap<Dim, std::optional<Variables<tmpl::list<
50 : evolution::dg::Tags::MagnitudeOfNormal,
51 : evolution::dg::Tags::NormalCovector<Dim>>>>>*>
52 : normal_covector_and_magnitude_ptr,
53 : const gsl::not_null<
54 : DirectionalIdMap<Dim, evolution::dg::MortarDataHolder<Dim>>*>
55 : mortar_data_ptr,
56 : const gsl::not_null<gsl::span<double>*> face_temporaries,
57 : const gsl::not_null<gsl::span<double>*> packaged_data_buffer,
58 : const BoundaryCorrection& boundary_correction,
59 : const Variables<typename System::variables_tag::tags_list>&
60 : volume_evolved_vars,
61 : const Variables<
62 : db::wrap_tags_in<::Tags::Flux, typename System::flux_variables,
63 : tmpl::size_t<Dim>, Frame::Inertial>>& volume_fluxes,
64 : const Variables<TemporaryTags>& volume_temporaries,
65 : const Variables<get_primitive_vars_tags_from_system<System>>* const
66 : volume_primitive_variables,
67 : const Element<Dim>& element, const Mesh<Dim>& volume_mesh,
68 : const DirectionalIdMap<Dim, Mesh<Dim - 1>>& mortar_meshes,
69 : const DirectionalIdMap<Dim, MortarInfo<Dim>>& mortar_infos,
70 : const domain::CoordinateMapBase<Frame::Grid, Frame::Inertial, Dim>&
71 : moving_mesh_map,
72 : const std::optional<tnsr::I<DataVector, Dim>>& volume_mesh_velocity,
73 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
74 : Frame::Inertial>& volume_inverse_jacobian,
75 : const PackageDataVolumeArgs&... package_data_volume_args) {
76 : using variables_tags = typename System::variables_tag::tags_list;
77 : using flux_variables = typename System::flux_variables;
78 : using fluxes_tags = db::wrap_tags_in<::Tags::Flux, flux_variables,
79 : tmpl::size_t<Dim>, Frame::Inertial>;
80 : using temporary_tags_for_face =
81 : typename BoundaryCorrection::dg_package_data_temporary_tags;
82 : using primitive_tags_for_face = typename detail::get_primitive_vars<
83 : System::has_primitive_and_conservative_vars>::
84 : template f<BoundaryCorrection>;
85 : using mortar_tags_list = typename BoundaryCorrection::dg_package_field_tags;
86 :
87 : using dg_package_data_projected_tags =
88 : tmpl::append<variables_tags, fluxes_tags, temporary_tags_for_face,
89 : primitive_tags_for_face>;
90 : using FieldsOnFace = Variables<tmpl::remove_duplicates<tmpl::push_back<
91 : tmpl::append<dg_package_data_projected_tags,
92 : detail::inverse_spatial_metric_tag<System>>,
93 : detail::OneOverNormalVectorMagnitude, detail::NormalVector<Dim>>>>;
94 : FieldsOnFace fields_on_face{};
95 : std::optional<tnsr::I<DataVector, Dim>> face_mesh_velocity{};
96 : for (const auto& [direction, neighbors_in_direction] : element.neighbors()) {
97 : const Mesh<Dim - 1> face_mesh =
98 : volume_mesh.slice_away(direction.dimension());
99 :
100 : // The face_temporaries buffer is guaranteed to be big enough because we
101 : // allocated it in ComputeTimeDerivative with the max number of grid points
102 : // over all faces. We still check anyways in Debug mode to be safe
103 : ASSERT(face_temporaries->size() >=
104 : FieldsOnFace::number_of_independent_components *
105 : face_mesh.number_of_grid_points(),
106 : "The buffer for computing fields on faces which was allocated in "
107 : "ComputeTimeDerivative is not large enough. It's size is "
108 : << face_temporaries->size() << ", but needs to be at least "
109 : << FieldsOnFace::number_of_independent_components *
110 : face_mesh.number_of_grid_points());
111 :
112 : fields_on_face.set_data_ref(face_temporaries->data(),
113 : FieldsOnFace::number_of_independent_components *
114 : face_mesh.number_of_grid_points());
115 :
116 : // We may not need to bring the volume fluxes or temporaries to the
117 : // boundary since that depends on the specific boundary correction we
118 : // are using. Silence compilers warnings about them being unused.
119 : (void)volume_fluxes;
120 : (void)volume_temporaries;
121 :
122 : // This does the following:
123 : //
124 : // 1. Use a helper function to get data onto the faces. Done either by
125 : // slicing (Gauss-Lobatto points) or interpolation (Gauss points).
126 : // This is done using the `project_contiguous_data_to_boundary` and
127 : // `project_tensors_to_boundary` functions.
128 : //
129 : // 2. Invoke the boundary correction to get the packaged data. Note
130 : // that this is done on the *face* and NOT the mortar.
131 : //
132 : // 3. Project the packaged data onto the DG mortars (these might need
133 : // re-projection onto subcell mortars later).
134 :
135 : // Perform step 1
136 : ::dg::project_contiguous_data_to_boundary(make_not_null(&fields_on_face),
137 : volume_evolved_vars, volume_mesh,
138 : direction);
139 : if constexpr (tmpl::size<fluxes_tags>::value != 0) {
140 : ::dg::project_contiguous_data_to_boundary(make_not_null(&fields_on_face),
141 : volume_fluxes, volume_mesh,
142 : direction);
143 : }
144 : if constexpr (tmpl::size<tmpl::append<
145 : temporary_tags_for_face,
146 : detail::inverse_spatial_metric_tag<System>>>::value !=
147 : 0) {
148 : ::dg::project_tensors_to_boundary<tmpl::append<
149 : temporary_tags_for_face, detail::inverse_spatial_metric_tag<System>>>(
150 : make_not_null(&fields_on_face), volume_temporaries, volume_mesh,
151 : direction);
152 : }
153 : if constexpr (System::has_primitive_and_conservative_vars and
154 : tmpl::size<primitive_tags_for_face>::value != 0) {
155 : ASSERT(volume_primitive_variables != nullptr,
156 : "The volume primitive variables are not set even though the "
157 : "system has primitive variables.");
158 : ::dg::project_tensors_to_boundary<primitive_tags_for_face>(
159 : make_not_null(&fields_on_face), *volume_primitive_variables,
160 : volume_mesh, direction);
161 : } else {
162 : (void)volume_primitive_variables;
163 : }
164 : if (volume_mesh_velocity.has_value()) {
165 : if (not face_mesh_velocity.has_value() or
166 : (*face_mesh_velocity)[0].size() !=
167 : face_mesh.number_of_grid_points()) {
168 : face_mesh_velocity =
169 : tnsr::I<DataVector, Dim>{face_mesh.number_of_grid_points()};
170 : }
171 : ::dg::project_tensor_to_boundary(make_not_null(&*face_mesh_velocity),
172 : *volume_mesh_velocity, volume_mesh,
173 : direction);
174 : }
175 :
176 : // Normalize the normal vectors. We cache the unit normal covector For
177 : // flat geometry and static meshes.
178 : const bool mesh_is_moving = not moving_mesh_map.is_identity();
179 : if (auto& normal_covector_quantity =
180 : normal_covector_and_magnitude_ptr->at(direction);
181 : detail::has_inverse_spatial_metric_tag_v<System> or mesh_is_moving or
182 : not normal_covector_quantity.has_value()) {
183 : if (not normal_covector_quantity.has_value()) {
184 : normal_covector_quantity =
185 : Variables<tmpl::list<evolution::dg::Tags::MagnitudeOfNormal,
186 : evolution::dg::Tags::NormalCovector<Dim>>>{
187 : fields_on_face.number_of_grid_points()};
188 : }
189 : tnsr::i<DataVector, Dim> volume_unnormalized_normal_covector{};
190 :
191 : for (size_t inertial_index = 0; inertial_index < Dim; ++inertial_index) {
192 : volume_unnormalized_normal_covector.get(inertial_index)
193 : .set_data_ref(const_cast<double*>( // NOLINT
194 : volume_inverse_jacobian
195 : .get(direction.dimension(), inertial_index)
196 : .data()),
197 : volume_mesh.number_of_grid_points());
198 : }
199 : ::dg::project_tensor_to_boundary(
200 : make_not_null(&get<evolution::dg::Tags::NormalCovector<Dim>>(
201 : *normal_covector_quantity)),
202 : volume_unnormalized_normal_covector, volume_mesh, direction);
203 :
204 : if (direction.side() == Side::Lower) {
205 : for (auto& normal_covector_component :
206 : get<evolution::dg::Tags::NormalCovector<Dim>>(
207 : *normal_covector_quantity)) {
208 : normal_covector_component *= -1.0;
209 : }
210 : }
211 :
212 : detail::unit_normal_vector_and_covector_and_magnitude_impl<System>(
213 : make_not_null(&get<evolution::dg::Tags::MagnitudeOfNormal>(
214 : *normal_covector_quantity)),
215 : make_not_null(&get<evolution::dg::Tags::NormalCovector<Dim>>(
216 : *normal_covector_quantity)),
217 : make_not_null(&fields_on_face),
218 : get<evolution::dg::Tags::NormalCovector<Dim>>(
219 : *normal_covector_quantity));
220 : }
221 :
222 : // Perform step 2
223 : ASSERT(normal_covector_and_magnitude_ptr->at(direction).has_value(),
224 : "The magnitude of the normal vector and the unit normal "
225 : "covector have not been computed, even though they should "
226 : "have been. Direction: "
227 : << direction);
228 :
229 : const size_t total_face_size =
230 : face_mesh.number_of_grid_points() *
231 : Variables<mortar_tags_list>::number_of_independent_components;
232 : Variables<mortar_tags_list> packaged_data{};
233 :
234 : // This is the case where we only have one neighbor in this direction, so we
235 : // may or may not have to do any projection. If we don't have to do
236 : // projection, then we can use the local_mortar_data itself to calculate the
237 : // dg_package_data. However, if we need to project, then we hae to use the
238 : // packaged_data_buffer that was passed in.
239 : if (neighbors_in_direction.size() == 1) {
240 : const auto& neighbor = *neighbors_in_direction.begin();
241 : const auto mortar_id = DirectionalId<Dim>{direction, neighbor};
242 : const auto& mortar_mesh = mortar_meshes.at(mortar_id);
243 : const auto& mortar_size = mortar_infos.at(mortar_id).mortar_size();
244 :
245 : // Have to use packaged_data_buffer
246 : if (Spectral::needs_projection(face_mesh, mortar_mesh, mortar_size)) {
247 : // The face mesh will be assigned below along with ensuring the size of
248 : // the mortar data is correct
249 : packaged_data.set_data_ref(packaged_data_buffer->data(),
250 : total_face_size);
251 : } else {
252 : // Can use the local_mortar_data
253 : auto& local_mortar = mortar_data_ptr->at(mortar_id).local();
254 : local_mortar.face_mesh = face_mesh;
255 : local_mortar.mortar_mesh = mortar_mesh;
256 : // If this is the first time, initialize the data. If we don't do this,
257 : // then the DataVector will be non-owning which we don't want
258 : if (UNLIKELY(not local_mortar.mortar_data.has_value())) {
259 : local_mortar.mortar_data = DataVector{};
260 : }
261 :
262 : DataVector& local_mortar_data = local_mortar.mortar_data.value();
263 :
264 : // Do a destructive resize to account for potential p-refinement
265 : local_mortar_data.destructive_resize(total_face_size);
266 :
267 : packaged_data.set_data_ref(local_mortar_data.data(),
268 : local_mortar_data.size());
269 : }
270 : } else {
271 : // In this case, we have multiple neighbors in this direction so all will
272 : // need to project their data which means we use the
273 : // packaged_data_buffer to calculate the dg_package_data
274 : packaged_data.set_data_ref(packaged_data_buffer->data(), total_face_size);
275 : }
276 :
277 : detail::dg_package_data<System>(
278 : make_not_null(&packaged_data), boundary_correction, fields_on_face,
279 : get<evolution::dg::Tags::NormalCovector<Dim>>(
280 : *normal_covector_and_magnitude_ptr->at(direction)),
281 : face_mesh_velocity, dg_package_data_projected_tags{},
282 : package_data_volume_args...);
283 :
284 : // Perform step 3
285 : // This will only do something if
286 : // a) we have multiple neighbors in this direction
287 : // or
288 : // b) the one (and only) neighbor in this direction needed projection
289 : for (const auto& neighbor : neighbors_in_direction) {
290 : const DirectionalId<Dim> mortar_id{direction, neighbor};
291 : const auto& mortar_mesh = mortar_meshes.at(mortar_id);
292 : const auto& mortar_size = mortar_infos.at(mortar_id).mortar_size();
293 :
294 : if (Spectral::needs_projection(face_mesh, mortar_mesh, mortar_size)) {
295 : auto& local_mortar = mortar_data_ptr->at(mortar_id).local();
296 : local_mortar.face_mesh = face_mesh;
297 : local_mortar.mortar_mesh = mortar_mesh;
298 : // If this is the first time, initialize the data. If we don't do this,
299 : // then the DataVector will be non-owning which we don't want
300 : if (UNLIKELY(not local_mortar.mortar_data.has_value())) {
301 : local_mortar.mortar_data = DataVector{};
302 : }
303 :
304 : DataVector& local_mortar_data = local_mortar.mortar_data.value();
305 :
306 : // Do a destructive resize to account for potential p-refinement
307 : local_mortar_data.destructive_resize(
308 : mortar_mesh.number_of_grid_points() *
309 : Variables<mortar_tags_list>::number_of_independent_components);
310 :
311 : Variables<mortar_tags_list> projected_packaged_data{
312 : local_mortar_data.data(), local_mortar_data.size()};
313 : ::dg::project_to_mortar(make_not_null(&projected_packaged_data),
314 : packaged_data, face_mesh, mortar_mesh,
315 : mortar_size);
316 : }
317 : }
318 : }
319 : }
320 :
321 : template <typename System, size_t Dim, typename BoundaryCorrection,
322 : typename DbTagsList, typename... PackageDataVolumeTags>
323 : void internal_mortar_data(
324 : const gsl::not_null<db::DataBox<DbTagsList>*> box,
325 : const gsl::not_null<gsl::span<double>*> face_temporaries,
326 : const gsl::not_null<gsl::span<double>*> packaged_data_buffer,
327 : const BoundaryCorrection& boundary_correction,
328 : const Variables<typename System::variables_tag::tags_list>
329 : evolved_variables,
330 : const Variables<
331 : db::wrap_tags_in<::Tags::Flux, typename System::flux_variables,
332 : tmpl::size_t<Dim>, Frame::Inertial>>& volume_fluxes,
333 : const Variables<
334 : typename System::compute_volume_time_derivative_terms::temporary_tags>&
335 : temporaries,
336 : const Variables<get_primitive_vars_tags_from_system<System>>* const
337 : primitive_vars,
338 : tmpl::list<PackageDataVolumeTags...> /*meta*/) {
339 : db::mutate<evolution::dg::Tags::NormalCovectorAndMagnitude<Dim>,
340 : evolution::dg::Tags::MortarData<Dim>>(
341 : [&boundary_correction, &face_temporaries, &packaged_data_buffer,
342 : &element = db::get<domain::Tags::Element<Dim>>(*box), &evolved_variables,
343 : &logical_to_inertial_inverse_jacobian =
344 : db::get<domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
345 : Frame::Inertial>>(*box),
346 : &mesh = db::get<domain::Tags::Mesh<Dim>>(*box),
347 : &mesh_velocity = db::get<domain::Tags::MeshVelocity<Dim>>(*box),
348 : &mortar_meshes = db::get<Tags::MortarMesh<Dim>>(*box),
349 : &mortar_infos = db::get<Tags::MortarInfo<Dim>>(*box),
350 : &moving_mesh_map = db::get<domain::CoordinateMaps::Tags::CoordinateMap<
351 : Dim, Frame::Grid, Frame::Inertial>>(*box),
352 : &primitive_vars, &temporaries, &volume_fluxes](
353 : const auto normal_covector_and_magnitude_ptr,
354 : const auto mortar_data_ptr, const auto&... package_data_volume_args) {
355 : detail::internal_mortar_data_impl<System>(
356 : normal_covector_and_magnitude_ptr, mortar_data_ptr,
357 : face_temporaries, packaged_data_buffer, boundary_correction,
358 : evolved_variables, volume_fluxes, temporaries, primitive_vars,
359 : element, mesh, mortar_meshes, mortar_infos, moving_mesh_map,
360 : mesh_velocity, logical_to_inertial_inverse_jacobian,
361 : package_data_volume_args...);
362 : },
363 : box, db::get<PackageDataVolumeTags>(*box)...);
364 : }
365 : } // namespace evolution::dg::Actions::detail
|
