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