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