Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <array>
7 : #include <cstddef>
8 : #include <cstdint>
9 : #include <optional>
10 : #include <type_traits>
11 :
12 : #include "DataStructures/DataBox/AsAccess.hpp"
13 : #include "DataStructures/DataBox/DataBox.hpp"
14 : #include "DataStructures/DataBox/MetavariablesTag.hpp"
15 : #include "DataStructures/DataBox/PrefixHelpers.hpp"
16 : #include "DataStructures/DataBox/Prefixes.hpp"
17 : #include "DataStructures/DataVector.hpp"
18 : #include "DataStructures/TaggedContainers.hpp"
19 : #include "DataStructures/Tensor/Tensor.hpp"
20 : #include "DataStructures/Variables.hpp"
21 : #include "Domain/FunctionsOfTime/FunctionOfTime.hpp"
22 : #include "Domain/FunctionsOfTime/Tags.hpp"
23 : #include "Domain/Structure/Element.hpp"
24 : #include "Domain/Tags.hpp"
25 : #include "Evolution/BoundaryCorrection.hpp"
26 : #include "Evolution/BoundaryCorrectionTags.hpp"
27 : #include "Evolution/DgSubcell/CartesianFluxDivergence.hpp"
28 : #include "Evolution/DgSubcell/ComputeBoundaryTerms.hpp"
29 : #include "Evolution/DgSubcell/CorrectPackagedData.hpp"
30 : #include "Evolution/DgSubcell/Mesh.hpp"
31 : #include "Evolution/DgSubcell/Projection.hpp"
32 : #include "Evolution/DgSubcell/ReconstructionOrder.hpp"
33 : #include "Evolution/DgSubcell/SubcellOptions.hpp"
34 : #include "Evolution/DgSubcell/Tags/CellCenteredFlux.hpp"
35 : #include "Evolution/DgSubcell/Tags/Coordinates.hpp"
36 : #include "Evolution/DgSubcell/Tags/GhostZoneInverseJacobian.hpp"
37 : #include "Evolution/DgSubcell/Tags/Jacobians.hpp"
38 : #include "Evolution/DgSubcell/Tags/Mesh.hpp"
39 : #include "Evolution/DgSubcell/Tags/OnSubcellFaces.hpp"
40 : #include "Evolution/DgSubcell/Tags/SubcellOptions.hpp"
41 : #include "Evolution/DiscontinuousGalerkin/Actions/NormalCovectorAndMagnitude.hpp"
42 : #include "Evolution/DiscontinuousGalerkin/Actions/PackageDataImpl.hpp"
43 : #include "Evolution/DiscontinuousGalerkin/MortarTags.hpp"
44 : #include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/BoundaryConditionGhostData.hpp"
45 : #include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Reconstructor.hpp"
46 : #include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Tag.hpp"
47 : #include "Evolution/Systems/GrMhd/ValenciaDivClean/Fluxes.hpp"
48 : #include "Evolution/Systems/GrMhd/ValenciaDivClean/Sources.hpp"
49 : #include "Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/ComputeFluxes.hpp"
50 : #include "Evolution/Systems/GrMhd/ValenciaDivClean/System.hpp"
51 : #include "NumericalAlgorithms/FiniteDifference/DerivativeOrder.hpp"
52 : #include "NumericalAlgorithms/FiniteDifference/HighOrderFluxCorrection.hpp"
53 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
54 : #include "NumericalAlgorithms/Spectral/Parity.hpp"
55 : #include "PointwiseFunctions/Hydro/Tags.hpp"
56 : #include "Utilities/CallWithDynamicType.hpp"
57 : #include "Utilities/ErrorHandling/Assert.hpp"
58 : #include "Utilities/Gsl.hpp"
59 : #include "Utilities/MakeWithValue.hpp"
60 : #include "Utilities/TMPL.hpp"
61 :
62 : namespace grmhd::ValenciaDivClean::subcell {
63 : /*!
64 : * \brief Compute the time derivative on the subcell grid using FD
65 : * reconstruction.
66 : */
67 1 : struct TimeDerivative {
68 : template <typename DbTagsList>
69 0 : static void apply(const gsl::not_null<db::DataBox<DbTagsList>*> box) {
70 : using metavariables =
71 : typename std::decay_t<decltype(db::get<Parallel::Tags::Metavariables>(
72 : *box))>;
73 : using evolved_vars_tag = typename System::variables_tag;
74 : using evolved_vars_tags = typename evolved_vars_tag::tags_list;
75 : using prim_tags = typename System::primitive_variables_tag::tags_list;
76 : using recons_prim_tags = tmpl::push_back<
77 : prim_tags,
78 : hydro::Tags::LorentzFactorTimesSpatialVelocity<DataVector, 3>>;
79 : using fluxes_tags = db::wrap_tags_in<::Tags::Flux, evolved_vars_tags,
80 : tmpl::size_t<3>, Frame::Inertial>;
81 :
82 : ASSERT(
83 : (db::get<::domain::CoordinateMaps::Tags::CoordinateMap<
84 : 3, Frame::Grid, Frame::Inertial>>(*box))
85 : .is_identity(),
86 : "Moving mesh is only partly implemented in ValenciaDivClean. If you "
87 : "need this look at the complete implementation in GhValenciaDivClean. "
88 : "You will at least need to update the high-order boundary correction "
89 : "code to include the right normal vectors/Jacobians.");
90 :
91 : // Arrays representing the conormals used for
92 : // `cartesian_high_order_flux_corrections`. These are not normalized
93 : // by the magnitude of the conormals, but they do include a factor of the
94 : // determinant of the Jacobian, i.e. \tilde{n}_\hat{i} = J n_\hat{i}
95 : std::array<tnsr::i<DataVector, 3, Frame::Inertial>, 3> conormal;
96 : const auto& ghost_zone_inv_jac =
97 : db::get<evolution::dg::subcell::Tags::GhostZoneInverseJacobian<3>>(
98 : *box);
99 : std::array<DirectionMap<3, tnsr::i<DataVector, 3, Frame::Inertial>>, 3>
100 : ghost_cells_conormal;
101 :
102 : const Mesh<3>& subcell_mesh =
103 : db::get<evolution::dg::subcell::Tags::Mesh<3>>(*box);
104 : const Mesh<3>& dg_mesh = db::get<domain::Tags::Mesh<3>>(*box);
105 : const size_t comp_dim =
106 : evolution::dg::subcell::fd::get_computational_dim(subcell_mesh);
107 : evolution::dg::subcell::fd::verify_subcell_extents(subcell_mesh.extents());
108 :
109 : const size_t reconstructed_num_pts =
110 : (subcell_mesh.extents(0) + 1) *
111 : subcell_mesh.extents().slice_away(0).product();
112 :
113 : const tnsr::I<DataVector, 3, Frame::ElementLogical>&
114 : cell_centered_logical_coords =
115 : db::get<evolution::dg::subcell::Tags::Coordinates<
116 : 3, Frame::ElementLogical>>(*box);
117 : std::array<double, 3> one_over_delta_xi{};
118 : for (size_t i = 0; i < 3; ++i) {
119 : // Note: assumes isotropic extents
120 : gsl::at(one_over_delta_xi, i) =
121 : 1.0 / (get<0>(cell_centered_logical_coords)[1] -
122 : get<0>(cell_centered_logical_coords)[0]);
123 : }
124 :
125 : // Inverse jacobian, to be projected on faces
126 : const auto& inv_jacobian_dg =
127 : db::get<domain::Tags::InverseJacobian<3, Frame::ElementLogical,
128 : Frame::Inertial>>(*box);
129 : const auto& det_inv_jacobian_dg = db::get<
130 : domain::Tags::DetInvJacobian<Frame::ElementLogical, Frame::Inertial>>(
131 : *box);
132 :
133 : // Velocity of the moving mesh on the DG grid, if applicable.
134 : const std::optional<tnsr::I<DataVector, 3, Frame::Inertial>>&
135 : mesh_velocity_dg = db::get<domain::Tags::MeshVelocity<3>>(*box);
136 : const std::optional<Scalar<DataVector>>& div_mesh_velocity =
137 : db::get<domain::Tags::DivMeshVelocity>(*box);
138 :
139 : const grmhd::ValenciaDivClean::fd::Reconstructor& recons =
140 : db::get<grmhd::ValenciaDivClean::fd::Tags::Reconstructor>(*box);
141 :
142 : const Element<3>& element = db::get<domain::Tags::Element<3>>(*box);
143 : const auto fd_derivative_order =
144 : db::get<evolution::dg::subcell::Tags::SubcellOptions<3>>(*box)
145 : .finite_difference_derivative_order();
146 : std::optional<std::array<std::vector<std::uint8_t>, 3>>
147 : reconstruction_order_data{};
148 : std::optional<std::array<gsl::span<std::uint8_t>, 3>>
149 : reconstruction_order{};
150 : if (static_cast<int>(fd_derivative_order) < 0) {
151 : reconstruction_order_data = make_array<3>(std::vector<std::uint8_t>(
152 : (subcell_mesh.extents(0) + 2) * subcell_mesh.extents(1) *
153 : subcell_mesh.extents(2),
154 : std::numeric_limits<std::uint8_t>::max()));
155 : reconstruction_order = std::array<gsl::span<std::uint8_t>, 3>{};
156 : for (size_t i = 0; i < 3; ++i) {
157 : gsl::at(reconstruction_order.value(), i) = gsl::make_span(
158 : gsl::at(reconstruction_order_data.value(), i).data(),
159 : gsl::at(reconstruction_order_data.value(), i).size());
160 : }
161 : }
162 :
163 : const bool element_is_interior = element.external_boundaries().empty();
164 : constexpr bool subcell_enabled_at_external_boundary =
165 : metavariables::SubcellOptions::subcell_enabled_at_external_boundary;
166 :
167 : ASSERT(element_is_interior or subcell_enabled_at_external_boundary,
168 : "Subcell time derivative is called at a boundary element while "
169 : "using subcell is disabled at external boundaries."
170 : "ElementID "
171 : << element.id());
172 :
173 : // Now package the data and compute the correction
174 : const auto& boundary_correction =
175 : db::get<evolution::Tags::BoundaryCorrection>(*box);
176 : using derived_boundary_corrections =
177 : tmpl::at<typename metavariables::factory_creation::factory_classes,
178 : evolution::BoundaryCorrection>;
179 : std::array<Variables<evolved_vars_tags>, 3> boundary_corrections{};
180 :
181 : // If the element has external boundaries and subcell is enabled for
182 : // boundary elements, compute FD ghost data with a given boundary condition.
183 : if constexpr (subcell_enabled_at_external_boundary) {
184 : if (not element.external_boundaries().empty()) {
185 : fd::BoundaryConditionGhostData::apply(box, element, recons);
186 : }
187 : }
188 :
189 : call_with_dynamic_type<void, derived_boundary_corrections>(
190 : &boundary_correction, [&](const auto* derived_correction) {
191 : using DerivedCorrection = std::decay_t<decltype(*derived_correction)>;
192 : using dg_package_data_temporary_tags =
193 : typename DerivedCorrection::dg_package_data_temporary_tags;
194 : using dg_package_data_argument_tags = tmpl::append<
195 : evolved_vars_tags, recons_prim_tags, fluxes_tags,
196 : tmpl::remove_duplicates<tmpl::push_back<
197 : dg_package_data_temporary_tags,
198 : gr::Tags::SpatialMetric<DataVector, 3>,
199 : gr::Tags::SqrtDetSpatialMetric<DataVector>,
200 : gr::Tags::InverseSpatialMetric<DataVector, 3>,
201 : evolution::dg::Actions::detail::NormalVector<3>>>>;
202 : // Computed prims and cons on face via reconstruction
203 : auto package_data_argvars_lower_face = make_array<3>(
204 : Variables<dg_package_data_argument_tags>(reconstructed_num_pts));
205 : auto package_data_argvars_upper_face = make_array<3>(
206 : Variables<dg_package_data_argument_tags>(reconstructed_num_pts));
207 : // Copy over the face values of the metric quantities.
208 : using spacetime_vars_to_copy =
209 : tmpl::list<gr::Tags::Lapse<DataVector>,
210 : gr::Tags::Shift<DataVector, 3>,
211 : gr::Tags::SpatialMetric<DataVector, 3>,
212 : gr::Tags::SqrtDetSpatialMetric<DataVector>,
213 : gr::Tags::InverseSpatialMetric<DataVector, 3>>;
214 : tmpl::for_each<spacetime_vars_to_copy>(
215 : [&package_data_argvars_lower_face,
216 : &package_data_argvars_upper_face,
217 : &spacetime_vars_on_face =
218 : db::get<evolution::dg::subcell::Tags::OnSubcellFaces<
219 : typename System::flux_spacetime_variables_tag, 3>>(*box),
220 : comp_dim](auto tag_v) {
221 : using tag = tmpl::type_from<decltype(tag_v)>;
222 : for (size_t d = 0; d < comp_dim; ++d) { // comp_dim
223 : get<tag>(gsl::at(package_data_argvars_lower_face, d)) =
224 : get<tag>(gsl::at(spacetime_vars_on_face, d));
225 : get<tag>(gsl::at(package_data_argvars_upper_face, d)) =
226 : get<tag>(gsl::at(spacetime_vars_on_face, d));
227 : }
228 : });
229 :
230 : // Reconstruct data to the face
231 : call_with_dynamic_type<void, typename grmhd::ValenciaDivClean::fd::
232 : Reconstructor::creatable_classes>(
233 : &recons, [&box, &package_data_argvars_lower_face,
234 : &package_data_argvars_upper_face,
235 : &reconstruction_order](const auto& reconstructor) {
236 : using ReconstructorType =
237 : std::decay_t<decltype(*reconstructor)>;
238 : db::apply<
239 : typename ReconstructorType::reconstruction_argument_tags>(
240 : [&package_data_argvars_lower_face,
241 : &package_data_argvars_upper_face, &reconstructor,
242 : &reconstruction_order](const auto&... args) {
243 : if constexpr (ReconstructorType::use_adaptive_order) {
244 : reconstructor->reconstruct(
245 : make_not_null(&package_data_argvars_lower_face),
246 : make_not_null(&package_data_argvars_upper_face),
247 : make_not_null(&reconstruction_order), args...);
248 : } else {
249 : (void)reconstruction_order;
250 : reconstructor->reconstruct(
251 : make_not_null(&package_data_argvars_lower_face),
252 : make_not_null(&package_data_argvars_upper_face),
253 : args...);
254 : }
255 : },
256 : *box);
257 : });
258 :
259 : using dg_package_field_tags =
260 : typename DerivedCorrection::dg_package_field_tags;
261 : // Allocated outside for loop to reduce allocations
262 : Variables<dg_package_field_tags> upper_packaged_data{
263 : reconstructed_num_pts};
264 : Variables<dg_package_field_tags> lower_packaged_data{
265 : reconstructed_num_pts};
266 :
267 : // Compute fluxes on faces
268 : for (size_t i = 0; i < comp_dim; ++i) {
269 : // Build extents of mesh shifted by half a grid cell in direction i
270 : const unsigned long& num_subcells_1d = subcell_mesh.extents(0);
271 : Index<3> face_mesh_extents = subcell_mesh.extents();
272 : face_mesh_extents[i] = num_subcells_1d + 1;
273 :
274 : auto& vars_upper_face = gsl::at(package_data_argvars_upper_face, i);
275 : auto& vars_lower_face = gsl::at(package_data_argvars_lower_face, i);
276 : grmhd::ValenciaDivClean::subcell::compute_fluxes(
277 : make_not_null(&vars_upper_face));
278 : grmhd::ValenciaDivClean::subcell::compute_fluxes(
279 : make_not_null(&vars_lower_face));
280 :
281 : // Add moving mesh corrections to the fluxes, if needed
282 : std::optional<tnsr::I<DataVector, 3, Frame::Inertial>>
283 : mesh_velocity_on_face = {};
284 : if (mesh_velocity_dg.has_value()) {
285 : // Project mesh velocity on face mesh.
286 : // Can we get away with only doing the normal component? It
287 : // is also used in the packaged data...
288 : mesh_velocity_on_face = tnsr::I<DataVector, 3, Frame::Inertial>{
289 : reconstructed_num_pts};
290 : for (size_t j = 0; j < 3; j++) {
291 : // j^th component of the velocity on the i^th directed face
292 : mesh_velocity_on_face.value().get(j) =
293 : evolution::dg::subcell::fd::project_to_faces(
294 : mesh_velocity_dg.value().get(j), dg_mesh,
295 : face_mesh_extents, i,
296 : j == 0 ? Spectral::Parity::Odd
297 : : Spectral::Parity::Even);
298 : }
299 : tmpl::for_each<evolved_vars_tags>([&vars_upper_face,
300 : &vars_lower_face,
301 : &mesh_velocity_on_face](
302 : auto tag_v) {
303 : using tag = tmpl::type_from<decltype(tag_v)>;
304 : using flux_tag =
305 : ::Tags::Flux<tag, tmpl::size_t<3>, Frame::Inertial>;
306 : using FluxTensor = typename flux_tag::type;
307 : const auto& var_upper = get<tag>(vars_upper_face);
308 : const auto& var_lower = get<tag>(vars_lower_face);
309 : auto& flux_upper = get<flux_tag>(vars_upper_face);
310 : auto& flux_lower = get<flux_tag>(vars_lower_face);
311 : for (size_t storage_index = 0; storage_index < var_upper.size();
312 : ++storage_index) {
313 : const auto tensor_index =
314 : var_upper.get_tensor_index(storage_index);
315 : for (size_t j = 0; j < 3; j++) {
316 : const auto flux_storage_index =
317 : FluxTensor::get_storage_index(prepend(tensor_index, j));
318 : flux_upper[flux_storage_index] -=
319 : mesh_velocity_on_face.value().get(j) *
320 : var_upper[storage_index];
321 : flux_lower[flux_storage_index] -=
322 : mesh_velocity_on_face.value().get(j) *
323 : var_lower[storage_index];
324 : }
325 : }
326 : });
327 : }
328 :
329 : // Normal vectors in curved spacetime normalized by inverse
330 : // spatial metric. Note that we use the sign convention on
331 : // the normal vectors to be compatible with DG.
332 : //
333 : // Note that these normal vectors are on all faces inside the DG
334 : // element since there are a bunch of subcells. We don't use the
335 : // NormalCovectorAndMagnitude tag in the DataBox right now to avoid
336 : // conflicts with the DG solver. We can explore in the future if
337 : // it's possible to reuse that allocation.
338 : //
339 : // The unnormalized normal vector is
340 : // n_j = d \xi^{\hat i}/dx^j
341 : // with "i" the current face.
342 : tnsr::i<DataVector, 3, Frame::Inertial> lower_outward_conormal_face{
343 : reconstructed_num_pts, 0.0};
344 : tnsr::i<DataVector, 3, Frame::Inertial> conormal_in_dir{
345 : subcell_mesh.extents().product(), 0.0};
346 : for (size_t j = 0; j < 3; j++) {
347 : conormal_in_dir.get(j) = evolution::dg::subcell::fd::project(
348 : inv_jacobian_dg.get(i, j), dg_mesh, subcell_mesh.extents(),
349 : (i == 0) != (j == 0) ? Spectral::Parity::Odd
350 : : Spectral::Parity::Even);
351 : lower_outward_conormal_face.get(j) =
352 : evolution::dg::subcell::fd::project_to_faces(
353 : inv_jacobian_dg.get(i, j), dg_mesh, face_mesh_extents, i,
354 : (i == 0) != (j == 0) ? Spectral::Parity::Odd
355 : : Spectral::Parity::Even);
356 : }
357 : const auto det_inv_jacobian = evolution::dg::subcell::fd::project(
358 : get(det_inv_jacobian_dg), dg_mesh, subcell_mesh.extents(),
359 : Spectral::Parity::Even);
360 : const auto det_inv_jacobian_face =
361 : evolution::dg::subcell::fd::project_to_faces(
362 : get(det_inv_jacobian_dg), dg_mesh, face_mesh_extents, i,
363 : Spectral::Parity::Even);
364 :
365 : const Scalar<DataVector> normalization{sqrt(get(dot_product(
366 : lower_outward_conormal_face, lower_outward_conormal_face,
367 : get<gr::Tags::InverseSpatialMetric<DataVector, 3>>(
368 : vars_upper_face))))};
369 : for (size_t j = 0; j < 3; j++) {
370 : lower_outward_conormal_face.get(j) =
371 : lower_outward_conormal_face.get(j) / get(normalization);
372 : }
373 :
374 : tnsr::i<DataVector, 3, Frame::Inertial> upper_outward_conormal_face{
375 : reconstructed_num_pts, 0.0};
376 : for (size_t j = 0; j < 3; j++) {
377 : upper_outward_conormal_face.get(j) =
378 : -lower_outward_conormal_face.get(j);
379 : gsl::at(conormal, i).get(j) =
380 : conormal_in_dir.get(j) / det_inv_jacobian;
381 : }
382 :
383 : for (const auto side : {Side::Lower, Side::Upper}) {
384 : const Direction<3> direction{i, side};
385 : const auto& ghost_cells_grid_coords = get<
386 : evolution::dg::subcell::Tags::Coordinates<3, Frame::Grid>>(
387 : ghost_zone_inv_jac.at(direction));
388 : const auto& ghost_cells_grid_inv_jacobian =
389 : get<evolution::dg::subcell::fd::Tags::
390 : InverseJacobianLogicalToGrid<3>>(
391 : ghost_zone_inv_jac.at(direction));
392 : const auto& ghost_cells_inertial_inv_jacobian =
393 : db::get<domain::CoordinateMaps::Tags::CoordinateMap<
394 : 3, Frame::Grid, Frame::Inertial>>(*box)
395 : .inv_jacobian(
396 : ghost_cells_grid_coords, db::get<::Tags::Time>(*box),
397 : db::get<::domain::Tags::FunctionsOfTime>(*box));
398 : auto total_inv_jacobian =
399 : make_with_value<tnsr::Ij<DataVector, 3>>(
400 : ghost_cells_inertial_inv_jacobian, 0.0);
401 : for (size_t m = 0; m < 3; m++) {
402 : for (size_t n = 0; n < 3; n++) {
403 : for (size_t j = 0; j < 3; j++) {
404 : total_inv_jacobian.get(m, n) +=
405 : ghost_cells_grid_inv_jacobian.get(m, j) *
406 : ghost_cells_inertial_inv_jacobian.get(j, n);
407 : }
408 : }
409 : }
410 : const auto ghost_cells_det_inv_jacobian =
411 : determinant(total_inv_jacobian);
412 : tnsr::i<DataVector, 3, Frame::Inertial> tmp_ghost_cells_conormal{
413 : ghost_cells_grid_coords.size()};
414 : for (size_t j = 0; j < 3; j++) {
415 : tmp_ghost_cells_conormal.get(j) =
416 : total_inv_jacobian.get(i, j) /
417 : get(ghost_cells_det_inv_jacobian);
418 : }
419 : gsl::at(ghost_cells_conormal, i)
420 : .insert_or_assign(direction, tmp_ghost_cells_conormal);
421 : }
422 : // Note: we probably should compute the normal vector in addition to
423 : // the co-vector. Not a huge issue since we'll get an FPE right now
424 : // if it's used by a Riemann solver.
425 :
426 : // Compute the packaged data
427 : using dg_package_data_projected_tags = tmpl::append<
428 : evolved_vars_tags, fluxes_tags, dg_package_data_temporary_tags,
429 : typename DerivedCorrection::dg_package_data_primitive_tags>;
430 : evolution::dg::Actions::detail::dg_package_data<System>(
431 : make_not_null(&upper_packaged_data), *derived_correction,
432 : vars_upper_face, upper_outward_conormal_face,
433 : mesh_velocity_on_face, *box,
434 : typename DerivedCorrection::dg_package_data_volume_tags{},
435 : dg_package_data_projected_tags{});
436 :
437 : evolution::dg::Actions::detail::dg_package_data<System>(
438 : make_not_null(&lower_packaged_data), *derived_correction,
439 : vars_lower_face, lower_outward_conormal_face,
440 : mesh_velocity_on_face, *box,
441 : typename DerivedCorrection::dg_package_data_volume_tags{},
442 : dg_package_data_projected_tags{});
443 :
444 : // Now need to check if any of our neighbors are doing DG,
445 : // because if so then we need to use whatever boundary data
446 : // they sent instead of what we computed locally.
447 : //
448 : // Note: We could check this beforehand to avoid the extra
449 : // work of reconstruction and flux computations at the
450 : // boundaries.
451 : evolution::dg::subcell::correct_package_data<true>(
452 : make_not_null(&lower_packaged_data),
453 : make_not_null(&upper_packaged_data), i, element, subcell_mesh,
454 : db::get<evolution::dg::Tags::MortarData<3>>(*box), 0);
455 :
456 : // Compute the corrections on the faces. We only need to
457 : // compute this once because we can just flip the normal
458 : // vectors then
459 : gsl::at(boundary_corrections, i).initialize(reconstructed_num_pts);
460 : evolution::dg::subcell::compute_boundary_terms(
461 : make_not_null(&gsl::at(boundary_corrections, i)),
462 : *derived_correction, upper_packaged_data, lower_packaged_data,
463 : db::as_access(*box),
464 : typename DerivedCorrection::dg_boundary_terms_volume_tags{});
465 : // We need to multiply by the normal vector normalization
466 : gsl::at(boundary_corrections, i) *= get(normalization);
467 : // Also multiply by determinant of Jacobian, following Eq.(34)
468 : // of 2109.11645
469 : gsl::at(boundary_corrections, i) *= 1.0 / det_inv_jacobian_face;
470 : }
471 : });
472 :
473 : // Now compute the actual time derivatives.
474 : using variables_tag = typename System::variables_tag;
475 : using dt_variables_tag = db::add_tag_prefix<::Tags::dt, variables_tag>;
476 : const gsl::not_null<typename dt_variables_tag::type*> dt_vars_ptr =
477 : db::mutate<dt_variables_tag>(
478 : [](const auto local_dt_vars_ptr) { return local_dt_vars_ptr; },
479 : box);
480 : dt_vars_ptr->initialize(subcell_mesh.number_of_grid_points());
481 :
482 : using grmhd_source_tags =
483 : tmpl::transform<ValenciaDivClean::ComputeSources::return_tags,
484 : tmpl::bind<db::remove_tag_prefix, tmpl::_1>>;
485 : sources_impl(
486 : dt_vars_ptr, *box, grmhd_source_tags{},
487 : typename grmhd::ValenciaDivClean::ComputeSources::argument_tags{});
488 :
489 : // Zero GRMHD tags that don't have sources.
490 : tmpl::for_each<typename variables_tag::tags_list>(
491 : [&dt_vars_ptr](auto evolved_var_tag_v) {
492 : using evolved_var_tag = tmpl::type_from<decltype(evolved_var_tag_v)>;
493 : using dt_tag = ::Tags::dt<evolved_var_tag>;
494 : auto& dt_var = get<dt_tag>(*dt_vars_ptr);
495 : for (size_t i = 0; i < dt_var.size(); ++i) {
496 : if constexpr (not tmpl::list_contains_v<grmhd_source_tags,
497 : evolved_var_tag>) {
498 : dt_var[i] = 0.0;
499 : }
500 : }
501 : });
502 :
503 : // Correction to source terms due to moving mesh
504 : if (div_mesh_velocity.has_value()) {
505 : const DataVector div_mesh_velocity_subcell =
506 : evolution::dg::subcell::fd::project(div_mesh_velocity.value().get(),
507 : dg_mesh, subcell_mesh.extents(),
508 : Spectral::Parity::Even);
509 : const auto& evolved_vars = db::get<evolved_vars_tag>(*box);
510 :
511 : tmpl::for_each<typename variables_tag::tags_list>(
512 : [&dt_vars_ptr, &div_mesh_velocity_subcell,
513 : &evolved_vars](auto evolved_var_tag_v) {
514 : using evolved_var_tag =
515 : tmpl::type_from<decltype(evolved_var_tag_v)>;
516 : using dt_tag = ::Tags::dt<evolved_var_tag>;
517 : auto& dt_var = get<dt_tag>(*dt_vars_ptr);
518 : const auto& evolved_var = get<evolved_var_tag>(evolved_vars);
519 : for (size_t i = 0; i < dt_var.size(); ++i) {
520 : dt_var[i] -= div_mesh_velocity_subcell * evolved_var[i];
521 : }
522 : });
523 : }
524 :
525 : std::optional<std::array<Variables<evolved_vars_tags>, 3>>
526 : high_order_corrections{};
527 : ::fd::cartesian_high_order_flux_corrections(
528 : make_not_null(&high_order_corrections),
529 :
530 : db::get<evolution::dg::subcell::Tags::CellCenteredFlux<
531 : evolved_vars_tags, 3>>(*box),
532 : boundary_corrections, fd_derivative_order,
533 : db::get<evolution::dg::subcell::Tags::GhostDataForReconstruction<3>>(
534 : *box),
535 : subcell_mesh, recons.ghost_zone_size(),
536 : reconstruction_order.value_or(std::array<gsl::span<std::uint8_t>, 3>{}),
537 : false, conormal, ghost_cells_conormal);
538 :
539 : const auto& cell_centered_det_inv_jacobian = db::get<
540 : evolution::dg::subcell::fd::Tags::DetInverseJacobianLogicalToInertial>(
541 : *box);
542 : for (size_t dim = 0; dim < comp_dim; ++dim) {
543 : const auto& boundary_correction_in_axis =
544 : high_order_corrections.has_value()
545 : ? gsl::at(high_order_corrections.value(), dim)
546 : : gsl::at(boundary_corrections, dim);
547 : const double inverse_delta = gsl::at(one_over_delta_xi, dim);
548 : tmpl::for_each<typename variables_tag::tags_list>(
549 : [&dt_vars_ptr, &boundary_correction_in_axis,
550 : &cell_centered_det_inv_jacobian, dim, inverse_delta, &subcell_mesh,
551 : comp_dim, &box](auto evolved_var_tag_v) {
552 : using evolved_var_tag =
553 : tmpl::type_from<decltype(evolved_var_tag_v)>;
554 : using dt_tag = ::Tags::dt<evolved_var_tag>;
555 : auto& dt_var = get<dt_tag>(*dt_vars_ptr);
556 : const auto& var_correction =
557 : get<evolved_var_tag>(boundary_correction_in_axis);
558 : for (size_t i = 0; i < dt_var.size(); ++i) {
559 : if (comp_dim == 3) {
560 : evolution::dg::subcell::add_cartesian_flux_divergence(
561 : make_not_null(&dt_var[i]), inverse_delta,
562 : get(cell_centered_det_inv_jacobian), var_correction[i],
563 : subcell_mesh.extents(), dim);
564 : } else {
565 : evolution::dg::subcell::add_cartoon_cartesian_flux_divergence(
566 : make_not_null(&dt_var[i]), inverse_delta,
567 : get(cell_centered_det_inv_jacobian), var_correction[i],
568 : subcell_mesh.extents(), dim,
569 : db::get<evolution::dg::subcell::Tags::Coordinates<
570 : 3, Frame::Inertial>>(*box),
571 : get<domain::Tags::ElementMap<3, Frame::Grid>>(*box),
572 : get<domain::CoordinateMaps::Tags::CoordinateMap<
573 : 3, Frame::Grid, Frame::Inertial>>(*box),
574 : get<::Tags::Time>(*box),
575 : get<domain::Tags::FunctionsOfTime>(*box));
576 : }
577 : }
578 : });
579 : }
580 :
581 : evolution::dg::subcell::store_reconstruction_order_in_databox(
582 : box, reconstruction_order);
583 : }
584 :
585 : private:
586 : template <typename DtVarsList, typename DbTagsList, typename... SourcedTags,
587 : typename... ArgsTags>
588 0 : static void sources_impl(
589 : const gsl::not_null<Variables<DtVarsList>*> dt_vars_ptr,
590 : const db::DataBox<DbTagsList>& box, tmpl::list<SourcedTags...> /*meta*/,
591 : tmpl::list<ArgsTags...> /*meta*/) {
592 : grmhd::ValenciaDivClean::ComputeSources::apply(
593 : get<::Tags::dt<SourcedTags>>(dt_vars_ptr)..., get<ArgsTags>(box)...);
594 : }
595 : };
596 : } // namespace grmhd::ValenciaDivClean::subcell
|
