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 "DataStructures/VariablesTag.hpp"
18 : #include "Domain/CoordinateMaps/CoordinateMap.hpp"
19 : #include "Domain/FaceNormal.hpp"
20 : #include "Domain/Structure/Direction.hpp"
21 : #include "Domain/Structure/DirectionMap.hpp"
22 : #include "Domain/Structure/DirectionalIdMap.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/SegmentSize.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, bool ComputeAuxiliary = false,
47 : typename BoundaryCorrection, typename EvolvedVariablesTags,
48 : typename TemporaryTags, typename... PackageDataVolumeArgs>
49 : void internal_mortar_data_impl(
50 : const gsl::not_null<
51 : DirectionMap<Dim, std::optional<Variables<tmpl::list<
52 : evolution::dg::Tags::MagnitudeOfNormal,
53 : evolution::dg::Tags::NormalCovector<Dim>>>>>*>
54 : normal_covector_and_magnitude_ptr,
55 : const gsl::not_null<
56 : DirectionalIdMap<Dim, evolution::dg::MortarDataHolder<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<EvolvedVariablesTags>& volume_evolved_vars,
62 : // The physical boundary correction may read the auxiliary variables,
63 : // so they are projected to the face for the physical pass.
64 : // No-op when there are no auxiliary variables or during the
65 : // auxiliary pass (which computes them but does not read them).
66 : const Variables<
67 : get_auxiliary_variables_or_default_t<System, tmpl::list<>>>* const
68 : volume_auxiliary_variables,
69 : const Variables<
70 : db::wrap_tags_in<::Tags::Flux, typename System::flux_variables,
71 : tmpl::size_t<Dim>, Frame::Inertial>>& volume_fluxes,
72 : const Variables<TemporaryTags>& volume_temporaries,
73 : const Variables<get_primitive_vars_tags_from_system<System>>* const
74 : volume_primitive_variables,
75 : const Element<Dim>& element, const Mesh<Dim>& volume_mesh,
76 : const DirectionalIdMap<Dim, Mesh<Dim - 1>>& mortar_meshes,
77 : const DirectionalIdMap<Dim, MortarInfo<Dim>>& mortar_infos,
78 : const domain::CoordinateMapBase<Frame::Grid, Frame::Inertial, Dim>&
79 : moving_mesh_map,
80 : const std::optional<tnsr::I<DataVector, Dim>>& volume_mesh_velocity,
81 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
82 : Frame::Inertial>& volume_inverse_jacobian,
83 : const PackageDataVolumeArgs&... package_data_volume_args) {
84 : using variables_tags = EvolvedVariablesTags;
85 : using auxiliary_variables =
86 : get_auxiliary_variables_or_default_t<System, tmpl::list<>>;
87 : using flux_variables = typename System::flux_variables;
88 : using fluxes_tags = db::wrap_tags_in<::Tags::Flux, flux_variables,
89 : tmpl::size_t<Dim>, Frame::Inertial>;
90 : using temporary_tags_for_face = tmpl::conditional_t<
91 : ComputeAuxiliary,
92 : get_dg_auxiliary_package_data_temporary_tags_or_default_t<
93 : BoundaryCorrection, tmpl::list<>>,
94 : typename BoundaryCorrection::dg_package_data_temporary_tags>;
95 : using primitive_tags_for_face = typename detail::get_primitive_vars<
96 : System::has_primitive_and_conservative_vars>::
97 : template f<BoundaryCorrection>;
98 : using mortar_tags_list =
99 : tmpl::conditional_t<ComputeAuxiliary,
100 : get_dg_auxiliary_package_field_tags_or_default_t<
101 : BoundaryCorrection, tmpl::list<>>,
102 : typename BoundaryCorrection::dg_package_field_tags>;
103 : using projected_auxiliary_vars_tags =
104 : tmpl::conditional_t<ComputeAuxiliary, tmpl::list<>, auxiliary_variables>;
105 : using dg_package_data_projected_tags =
106 : tmpl::append<variables_tags, projected_auxiliary_vars_tags, fluxes_tags,
107 : temporary_tags_for_face, primitive_tags_for_face>;
108 : using FieldsOnFace = Variables<tmpl::remove_duplicates<tmpl::push_back<
109 : tmpl::append<dg_package_data_projected_tags,
110 : detail::inverse_spatial_metric_tag<System>>,
111 : detail::OneOverNormalVectorMagnitude, detail::NormalVector<Dim>>>>;
112 : FieldsOnFace fields_on_face{};
113 : std::optional<tnsr::I<DataVector, Dim>> face_mesh_velocity{};
114 : for (const auto& [direction, neighbors_in_direction] : element.neighbors()) {
115 : const Mesh<Dim - 1> face_mesh =
116 : volume_mesh.on_interface(direction.dimension());
117 :
118 : // The face_temporaries buffer is guaranteed to be big enough because we
119 : // allocated it in ComputeTimeDerivative with the max number of grid points
120 : // over all faces. We still check anyways in Debug mode to be safe
121 : ASSERT(face_temporaries->size() >=
122 : FieldsOnFace::number_of_independent_components *
123 : face_mesh.number_of_grid_points(),
124 : "The buffer for computing fields on faces which was allocated in "
125 : "ComputeTimeDerivative is not large enough. It's size is "
126 : << face_temporaries->size() << ", but needs to be at least "
127 : << FieldsOnFace::number_of_independent_components *
128 : face_mesh.number_of_grid_points());
129 :
130 : fields_on_face.set_data_ref(face_temporaries->data(),
131 : FieldsOnFace::number_of_independent_components *
132 : face_mesh.number_of_grid_points());
133 :
134 : // We may not need to bring the volume fluxes or temporaries to the
135 : // boundary since that depends on the specific boundary correction we
136 : // are using. Silence compilers warnings about them being unused.
137 : (void)volume_fluxes;
138 : (void)volume_temporaries;
139 :
140 : // This does the following:
141 : //
142 : // 1. Use a helper function to get data onto the faces. Done either by
143 : // slicing (Gauss-Lobatto points) or interpolation (Gauss points).
144 : // This is done using the `project_contiguous_data_to_boundary` and
145 : // `project_tensors_to_boundary` functions.
146 : //
147 : // 2. Invoke the boundary correction to get the packaged data. Note
148 : // that this is done on the *face* and NOT the mortar.
149 : //
150 : // 3. Project the packaged data onto the DG mortars (these might need
151 : // re-projection onto subcell mortars later).
152 :
153 : // Perform step 1
154 : ::dg::project_contiguous_data_to_boundary(make_not_null(&fields_on_face),
155 : volume_evolved_vars, volume_mesh,
156 : direction);
157 : if constexpr (not ComputeAuxiliary and
158 : tmpl::size<auxiliary_variables>::value != 0) {
159 : ASSERT(volume_auxiliary_variables != nullptr,
160 : "The auxiliary variables must be provided to the physical pass "
161 : "when the system has auxiliary variables.");
162 : ::dg::project_tensors_to_boundary<auxiliary_variables>(
163 : make_not_null(&fields_on_face), *volume_auxiliary_variables,
164 : volume_mesh, direction);
165 : } else {
166 : (void)volume_auxiliary_variables;
167 : }
168 : if constexpr (tmpl::size<fluxes_tags>::value != 0) {
169 : ::dg::project_contiguous_data_to_boundary(make_not_null(&fields_on_face),
170 : volume_fluxes, volume_mesh,
171 : direction);
172 : }
173 : if constexpr (tmpl::size<tmpl::append<
174 : temporary_tags_for_face,
175 : detail::inverse_spatial_metric_tag<System>>>::value !=
176 : 0) {
177 : ::dg::project_tensors_to_boundary<tmpl::append<
178 : temporary_tags_for_face, detail::inverse_spatial_metric_tag<System>>>(
179 : make_not_null(&fields_on_face), volume_temporaries, volume_mesh,
180 : direction);
181 : }
182 : if constexpr (System::has_primitive_and_conservative_vars and
183 : tmpl::size<primitive_tags_for_face>::value != 0) {
184 : ASSERT(volume_primitive_variables != nullptr,
185 : "The volume primitive variables are not set even though the "
186 : "system has primitive variables.");
187 : ::dg::project_tensors_to_boundary<primitive_tags_for_face>(
188 : make_not_null(&fields_on_face), *volume_primitive_variables,
189 : volume_mesh, direction);
190 : } else {
191 : (void)volume_primitive_variables;
192 : }
193 : if (volume_mesh_velocity.has_value()) {
194 : if (not face_mesh_velocity.has_value() or
195 : (*face_mesh_velocity)[0].size() !=
196 : face_mesh.number_of_grid_points()) {
197 : face_mesh_velocity =
198 : tnsr::I<DataVector, Dim>{face_mesh.number_of_grid_points()};
199 : }
200 : ::dg::project_tensor_to_boundary(make_not_null(&*face_mesh_velocity),
201 : *volume_mesh_velocity, volume_mesh,
202 : direction);
203 : }
204 :
205 : // Normalize the normal vectors. We cache the unit normal covector For
206 : // flat geometry and static meshes.
207 : const bool mesh_is_moving = not moving_mesh_map.is_identity();
208 : if (auto& normal_covector_quantity =
209 : normal_covector_and_magnitude_ptr->at(direction);
210 : detail::has_inverse_spatial_metric_tag_v<System> or mesh_is_moving or
211 : not normal_covector_quantity.has_value()) {
212 : if (not normal_covector_quantity.has_value()) {
213 : normal_covector_quantity =
214 : Variables<tmpl::list<evolution::dg::Tags::MagnitudeOfNormal,
215 : evolution::dg::Tags::NormalCovector<Dim>>>{
216 : fields_on_face.number_of_grid_points()};
217 : }
218 : tnsr::i<DataVector, Dim> volume_unnormalized_normal_covector{};
219 :
220 : for (size_t inertial_index = 0; inertial_index < Dim; ++inertial_index) {
221 : volume_unnormalized_normal_covector.get(inertial_index)
222 : .set_data_ref(const_cast<double*>( // NOLINT
223 : volume_inverse_jacobian
224 : .get(direction.dimension(), inertial_index)
225 : .data()),
226 : volume_mesh.number_of_grid_points());
227 : }
228 : ::dg::project_tensor_to_boundary(
229 : make_not_null(&get<evolution::dg::Tags::NormalCovector<Dim>>(
230 : *normal_covector_quantity)),
231 : volume_unnormalized_normal_covector, volume_mesh, direction);
232 :
233 : if (direction.side() == Side::Lower) {
234 : for (auto& normal_covector_component :
235 : get<evolution::dg::Tags::NormalCovector<Dim>>(
236 : *normal_covector_quantity)) {
237 : normal_covector_component *= -1.0;
238 : }
239 : }
240 :
241 : detail::unit_normal_vector_and_covector_and_magnitude_impl<System>(
242 : make_not_null(&get<evolution::dg::Tags::MagnitudeOfNormal>(
243 : *normal_covector_quantity)),
244 : make_not_null(&get<evolution::dg::Tags::NormalCovector<Dim>>(
245 : *normal_covector_quantity)),
246 : make_not_null(&fields_on_face),
247 : get<evolution::dg::Tags::NormalCovector<Dim>>(
248 : *normal_covector_quantity));
249 : }
250 :
251 : // Perform step 2
252 : ASSERT(normal_covector_and_magnitude_ptr->at(direction).has_value(),
253 : "The magnitude of the normal vector and the unit normal "
254 : "covector have not been computed, even though they should "
255 : "have been. Direction: "
256 : << direction);
257 :
258 : const size_t total_face_size =
259 : face_mesh.number_of_grid_points() *
260 : Variables<mortar_tags_list>::number_of_independent_components;
261 : Variables<mortar_tags_list> packaged_data{};
262 :
263 : // If there are multiple non-conforming neighbors, we only create a single
264 : // mortar labeled by the host ElementId. This is done because the data
265 : // from all neighbors will be combined onto a single mortar as it makes no
266 : // sense to have multiple mortars between non-conforming Elements.
267 : const bool has_multiple_non_conforming_neighbors =
268 : neighbors_in_direction.size() > 1 and
269 : not neighbors_in_direction.are_conforming();
270 : if (neighbors_in_direction.size() == 1 or
271 : not neighbors_in_direction.are_conforming()) {
272 : const auto& neighbor = *neighbors_in_direction.begin();
273 : const DirectionalId<Dim> mortar_id{
274 : direction,
275 : has_multiple_non_conforming_neighbors ? element.id() : neighbor};
276 : const auto& mortar_mesh = mortar_meshes.at(mortar_id);
277 : const auto& mortar_size = mortar_infos.at(mortar_id).mortar_size();
278 :
279 : // If we only have one conforming neighbor in this direction, we may or
280 : // may not have to do any projection. If we don't have to do projection,
281 : // then we can use the local_mortar_data itself to calculate the
282 : // dg_package_data. However, if we need to project, then we have to use
283 : // the packaged_data_buffer that was passed in.
284 : if (neighbors_in_direction.are_conforming() and
285 : Spectral::needs_projection(face_mesh, mortar_mesh, mortar_size)) {
286 : // The face mesh will be assigned below along with ensuring the size of
287 : // the mortar data is correct
288 : packaged_data.set_data_ref(packaged_data_buffer->data(),
289 : total_face_size);
290 : } else {
291 : // Can use the local_mortar_data
292 : auto& local_mortar = mortar_data_ptr->at(mortar_id).local();
293 : local_mortar.face_mesh = face_mesh;
294 : local_mortar.mortar_mesh = mortar_mesh;
295 : // If this is the first time, initialize the data. If we don't do this,
296 : // then the DataVector will be non-owning which we don't want
297 : if (UNLIKELY(not local_mortar.mortar_data.has_value())) {
298 : local_mortar.mortar_data = DataVector{};
299 : }
300 :
301 : DataVector& local_mortar_data = local_mortar.mortar_data.value();
302 :
303 : // Do a destructive resize to account for potential p-refinement
304 : local_mortar_data.destructive_resize(total_face_size);
305 :
306 : packaged_data.set_data_ref(local_mortar_data.data(),
307 : local_mortar_data.size());
308 : }
309 : } else {
310 : // In this case, we have multiple conforming neighbors in this direction
311 : // so all will need to project their data which means we use the
312 : // packaged_data_buffer to calculate the dg_package_data
313 : packaged_data.set_data_ref(packaged_data_buffer->data(), total_face_size);
314 : }
315 :
316 : if constexpr (ComputeAuxiliary) {
317 : detail::dg_auxiliary_package_data<System>(
318 : make_not_null(&packaged_data), boundary_correction, fields_on_face,
319 : get<evolution::dg::Tags::NormalCovector<Dim>>(
320 : *normal_covector_and_magnitude_ptr->at(direction)),
321 : face_mesh_velocity, dg_package_data_projected_tags{},
322 : package_data_volume_args...);
323 : } else {
324 : detail::dg_package_data<System>(
325 : make_not_null(&packaged_data), boundary_correction, fields_on_face,
326 : get<evolution::dg::Tags::NormalCovector<Dim>>(
327 : *normal_covector_and_magnitude_ptr->at(direction)),
328 : face_mesh_velocity, dg_package_data_projected_tags{},
329 : package_data_volume_args...);
330 : }
331 :
332 : // Perform step 3
333 : // This will only do something if neighbors are conforming and either
334 : // a) we have multiple neighbors in this direction
335 : // or
336 : // b) the one (and only) neighbor in this direction needed projection
337 : if (neighbors_in_direction.are_conforming()) {
338 : for (const auto& neighbor : neighbors_in_direction) {
339 : const DirectionalId<Dim> mortar_id{direction, neighbor};
340 : const auto& mortar_mesh = mortar_meshes.at(mortar_id);
341 : const auto& mortar_size = mortar_infos.at(mortar_id).mortar_size();
342 :
343 : if (Spectral::needs_projection(face_mesh, mortar_mesh, mortar_size)) {
344 : auto& local_mortar = mortar_data_ptr->at(mortar_id).local();
345 : local_mortar.face_mesh = face_mesh;
346 : local_mortar.mortar_mesh = mortar_mesh;
347 : // If this is the first time, initialize the data. If we don't do
348 : // this, then the DataVector will be non-owning which we don't want
349 : if (UNLIKELY(not local_mortar.mortar_data.has_value())) {
350 : local_mortar.mortar_data = DataVector{};
351 : }
352 :
353 : DataVector& local_mortar_data = local_mortar.mortar_data.value();
354 :
355 : // Do a destructive resize to account for potential p-refinement
356 : local_mortar_data.destructive_resize(
357 : mortar_mesh.number_of_grid_points() *
358 : Variables<mortar_tags_list>::number_of_independent_components);
359 :
360 : Variables<mortar_tags_list> projected_packaged_data{
361 : local_mortar_data.data(), local_mortar_data.size()};
362 : ::dg::project_to_mortar(make_not_null(&projected_packaged_data),
363 : packaged_data, face_mesh, mortar_mesh,
364 : mortar_size);
365 : }
366 : }
367 : }
368 : }
369 : }
370 :
371 : template <typename System, size_t Dim, bool ComputeAuxiliary = false,
372 : typename BoundaryCorrection, typename DbTagsList,
373 : typename EvolvedVariablesTags, typename... PackageDataVolumeTags>
374 : void internal_mortar_data(
375 : const gsl::not_null<db::DataBox<DbTagsList>*> box,
376 : const gsl::not_null<gsl::span<double>*> face_temporaries,
377 : const gsl::not_null<gsl::span<double>*> packaged_data_buffer,
378 : const BoundaryCorrection& boundary_correction,
379 : const Variables<EvolvedVariablesTags>& evolved_variables,
380 : const Variables<
381 : db::wrap_tags_in<::Tags::Flux, typename System::flux_variables,
382 : tmpl::size_t<Dim>, Frame::Inertial>>& volume_fluxes,
383 : const Variables<
384 : typename System::compute_volume_time_derivative_terms::temporary_tags>&
385 : temporaries,
386 : const Variables<get_primitive_vars_tags_from_system<System>>* const
387 : primitive_vars,
388 : tmpl::list<PackageDataVolumeTags...> /*meta*/) {
389 : using auxiliary_variables =
390 : get_auxiliary_variables_or_default_t<System, tmpl::list<>>;
391 : const Variables<auxiliary_variables>* volume_auxiliary_variables = nullptr;
392 : if constexpr (not ComputeAuxiliary and
393 : tmpl::size<auxiliary_variables>::value != 0) {
394 : volume_auxiliary_variables =
395 : &db::get<::Tags::Variables<auxiliary_variables>>(*box);
396 : }
397 : db::mutate<evolution::dg::Tags::NormalCovectorAndMagnitude<Dim>,
398 : evolution::dg::Tags::MortarData<Dim>>(
399 : [&boundary_correction, &face_temporaries, &packaged_data_buffer,
400 : &element = db::get<domain::Tags::Element<Dim>>(*box), &evolved_variables,
401 : volume_auxiliary_variables,
402 : &logical_to_inertial_inverse_jacobian =
403 : db::get<domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
404 : Frame::Inertial>>(*box),
405 : &mesh = db::get<domain::Tags::Mesh<Dim>>(*box),
406 : &mesh_velocity = db::get<domain::Tags::MeshVelocity<Dim>>(*box),
407 : &mortar_meshes = db::get<Tags::MortarMesh<Dim>>(*box),
408 : &mortar_infos = db::get<Tags::MortarInfo<Dim>>(*box),
409 : &moving_mesh_map = db::get<domain::CoordinateMaps::Tags::CoordinateMap<
410 : Dim, Frame::Grid, Frame::Inertial>>(*box),
411 : &primitive_vars, &temporaries, &volume_fluxes](
412 : const auto normal_covector_and_magnitude_ptr,
413 : const auto mortar_data_ptr, const auto&... package_data_volume_args) {
414 : detail::internal_mortar_data_impl<System, Dim, ComputeAuxiliary>(
415 : normal_covector_and_magnitude_ptr, mortar_data_ptr,
416 : face_temporaries, packaged_data_buffer, boundary_correction,
417 : evolved_variables, volume_auxiliary_variables, volume_fluxes,
418 : temporaries, primitive_vars, element, mesh, mortar_meshes,
419 : mortar_infos, moving_mesh_map, mesh_velocity,
420 : logical_to_inertial_inverse_jacobian, package_data_volume_args...);
421 : },
422 : box, db::get<PackageDataVolumeTags>(*box)...);
423 : }
424 : } // namespace evolution::dg::Actions::detail
|
