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 : // If there are multiple non-conforming neighbors, we only create a single
235 : // mortar labeled by the host ElementId. This is done because the data
236 : // from all neighbors will be combined onto a single mortar as it makes no
237 : // sense to have multiple mortars between non-conforming Elements.
238 : const bool has_multiple_non_conforming_neighbors =
239 : neighbors_in_direction.size() > 1 and
240 : not neighbors_in_direction.are_conforming();
241 : if (neighbors_in_direction.size() == 1 or
242 : not neighbors_in_direction.are_conforming()) {
243 : const auto& neighbor = *neighbors_in_direction.begin();
244 : const DirectionalId<Dim> mortar_id{
245 : direction,
246 : has_multiple_non_conforming_neighbors ? element.id() : neighbor};
247 : const auto& mortar_mesh = mortar_meshes.at(mortar_id);
248 : const auto& mortar_size = mortar_infos.at(mortar_id).mortar_size();
249 :
250 : // If we only have one conforming neighbor in this direction, we may or
251 : // may not have to do any projection. If we don't have to do projection,
252 : // then we can use the local_mortar_data itself to calculate the
253 : // dg_package_data. However, if we need to project, then we have to use
254 : // the packaged_data_buffer that was passed in.
255 : if (neighbors_in_direction.are_conforming() and
256 : Spectral::needs_projection(face_mesh, mortar_mesh, mortar_size)) {
257 : // The face mesh will be assigned below along with ensuring the size of
258 : // the mortar data is correct
259 : packaged_data.set_data_ref(packaged_data_buffer->data(),
260 : total_face_size);
261 : } else {
262 : // Can use the local_mortar_data
263 : auto& local_mortar = mortar_data_ptr->at(mortar_id).local();
264 : local_mortar.face_mesh = face_mesh;
265 : local_mortar.mortar_mesh = mortar_mesh;
266 : // If this is the first time, initialize the data. If we don't do this,
267 : // then the DataVector will be non-owning which we don't want
268 : if (UNLIKELY(not local_mortar.mortar_data.has_value())) {
269 : local_mortar.mortar_data = DataVector{};
270 : }
271 :
272 : DataVector& local_mortar_data = local_mortar.mortar_data.value();
273 :
274 : // Do a destructive resize to account for potential p-refinement
275 : local_mortar_data.destructive_resize(total_face_size);
276 :
277 : packaged_data.set_data_ref(local_mortar_data.data(),
278 : local_mortar_data.size());
279 : }
280 : } else {
281 : // In this case, we have multiple conforming neighbors in this direction
282 : // so all will need to project their data which means we use the
283 : // packaged_data_buffer to calculate the dg_package_data
284 : packaged_data.set_data_ref(packaged_data_buffer->data(), total_face_size);
285 : }
286 :
287 : detail::dg_package_data<System>(
288 : make_not_null(&packaged_data), boundary_correction, fields_on_face,
289 : get<evolution::dg::Tags::NormalCovector<Dim>>(
290 : *normal_covector_and_magnitude_ptr->at(direction)),
291 : face_mesh_velocity, dg_package_data_projected_tags{},
292 : package_data_volume_args...);
293 :
294 : // Perform step 3
295 : // This will only do something if neighbors are conforming and either
296 : // a) we have multiple neighbors in this direction
297 : // or
298 : // b) the one (and only) neighbor in this direction needed projection
299 : if (neighbors_in_direction.are_conforming()) {
300 : for (const auto& neighbor : neighbors_in_direction) {
301 : const DirectionalId<Dim> mortar_id{direction, neighbor};
302 : const auto& mortar_mesh = mortar_meshes.at(mortar_id);
303 : const auto& mortar_size = mortar_infos.at(mortar_id).mortar_size();
304 :
305 : if (Spectral::needs_projection(face_mesh, mortar_mesh, mortar_size)) {
306 : auto& local_mortar = mortar_data_ptr->at(mortar_id).local();
307 : local_mortar.face_mesh = face_mesh;
308 : local_mortar.mortar_mesh = mortar_mesh;
309 : // If this is the first time, initialize the data. If we don't do
310 : // this, then the DataVector will be non-owning which we don't want
311 : if (UNLIKELY(not local_mortar.mortar_data.has_value())) {
312 : local_mortar.mortar_data = DataVector{};
313 : }
314 :
315 : DataVector& local_mortar_data = local_mortar.mortar_data.value();
316 :
317 : // Do a destructive resize to account for potential p-refinement
318 : local_mortar_data.destructive_resize(
319 : mortar_mesh.number_of_grid_points() *
320 : Variables<mortar_tags_list>::number_of_independent_components);
321 :
322 : Variables<mortar_tags_list> projected_packaged_data{
323 : local_mortar_data.data(), local_mortar_data.size()};
324 : ::dg::project_to_mortar(make_not_null(&projected_packaged_data),
325 : packaged_data, face_mesh, mortar_mesh,
326 : mortar_size);
327 : }
328 : }
329 : }
330 : }
331 : }
332 :
333 : template <typename System, size_t Dim, typename BoundaryCorrection,
334 : typename DbTagsList, typename... PackageDataVolumeTags>
335 : void internal_mortar_data(
336 : const gsl::not_null<db::DataBox<DbTagsList>*> box,
337 : const gsl::not_null<gsl::span<double>*> face_temporaries,
338 : const gsl::not_null<gsl::span<double>*> packaged_data_buffer,
339 : const BoundaryCorrection& boundary_correction,
340 : const Variables<typename System::variables_tag::tags_list>
341 : evolved_variables,
342 : const Variables<
343 : db::wrap_tags_in<::Tags::Flux, typename System::flux_variables,
344 : tmpl::size_t<Dim>, Frame::Inertial>>& volume_fluxes,
345 : const Variables<
346 : typename System::compute_volume_time_derivative_terms::temporary_tags>&
347 : temporaries,
348 : const Variables<get_primitive_vars_tags_from_system<System>>* const
349 : primitive_vars,
350 : tmpl::list<PackageDataVolumeTags...> /*meta*/) {
351 : db::mutate<evolution::dg::Tags::NormalCovectorAndMagnitude<Dim>,
352 : evolution::dg::Tags::MortarData<Dim>>(
353 : [&boundary_correction, &face_temporaries, &packaged_data_buffer,
354 : &element = db::get<domain::Tags::Element<Dim>>(*box), &evolved_variables,
355 : &logical_to_inertial_inverse_jacobian =
356 : db::get<domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
357 : Frame::Inertial>>(*box),
358 : &mesh = db::get<domain::Tags::Mesh<Dim>>(*box),
359 : &mesh_velocity = db::get<domain::Tags::MeshVelocity<Dim>>(*box),
360 : &mortar_meshes = db::get<Tags::MortarMesh<Dim>>(*box),
361 : &mortar_infos = db::get<Tags::MortarInfo<Dim>>(*box),
362 : &moving_mesh_map = db::get<domain::CoordinateMaps::Tags::CoordinateMap<
363 : Dim, Frame::Grid, Frame::Inertial>>(*box),
364 : &primitive_vars, &temporaries, &volume_fluxes](
365 : const auto normal_covector_and_magnitude_ptr,
366 : const auto mortar_data_ptr, const auto&... package_data_volume_args) {
367 : detail::internal_mortar_data_impl<System>(
368 : normal_covector_and_magnitude_ptr, mortar_data_ptr,
369 : face_temporaries, packaged_data_buffer, boundary_correction,
370 : evolved_variables, volume_fluxes, temporaries, primitive_vars,
371 : element, mesh, mortar_meshes, mortar_infos, moving_mesh_map,
372 : mesh_velocity, logical_to_inertial_inverse_jacobian,
373 : package_data_volume_args...);
374 : },
375 : box, db::get<PackageDataVolumeTags>(*box)...);
376 : }
377 : } // namespace evolution::dg::Actions::detail
|
