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 <optional>
9 : #include <type_traits>
10 :
11 : #include "DataStructures/DataBox/AsAccess.hpp"
12 : #include "DataStructures/DataBox/DataBox.hpp"
13 : #include "DataStructures/DataBox/MetavariablesTag.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/Tags.hpp"
20 : #include "Domain/Structure/Element.hpp"
21 : #include "Domain/Tags.hpp"
22 : #include "Evolution/BoundaryCorrection.hpp"
23 : #include "Evolution/BoundaryCorrectionTags.hpp"
24 : #include "Evolution/DgSubcell/CartesianFluxDivergence.hpp"
25 : #include "Evolution/DgSubcell/ComputeBoundaryTerms.hpp"
26 : #include "Evolution/DgSubcell/CorrectPackagedData.hpp"
27 : #include "Evolution/DgSubcell/Tags/Coordinates.hpp"
28 : #include "Evolution/DgSubcell/Tags/Jacobians.hpp"
29 : #include "Evolution/DgSubcell/Tags/Mesh.hpp"
30 : #include "Evolution/DgSubcell/Tags/OnSubcellFaces.hpp"
31 : #include "Evolution/DiscontinuousGalerkin/Actions/NormalCovectorAndMagnitude.hpp"
32 : #include "Evolution/DiscontinuousGalerkin/Actions/PackageDataImpl.hpp"
33 : #include "Evolution/DiscontinuousGalerkin/MortarTags.hpp"
34 : #include "Evolution/Systems/NewtonianMhd/FiniteDifference/Reconstructor.hpp"
35 : #include "Evolution/Systems/NewtonianMhd/FiniteDifference/Tag.hpp"
36 : #include "Evolution/Systems/NewtonianMhd/Fluxes.hpp"
37 : #include "Evolution/Systems/NewtonianMhd/Subcell/ComputeFluxes.hpp"
38 : #include "Evolution/Systems/NewtonianMhd/System.hpp"
39 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
40 : #include "PointwiseFunctions/Hydro/Tags.hpp"
41 : #include "Utilities/CallWithDynamicType.hpp"
42 : #include "Utilities/ErrorHandling/Assert.hpp"
43 : #include "Utilities/Gsl.hpp"
44 : #include "Utilities/TMPL.hpp"
45 :
46 : namespace NewtonianMhd::subcell {
47 : /*!
48 : * \brief Compute the time derivative on the subcell grid using FD
49 : * reconstruction.
50 : *
51 : * The code makes the following unchecked assumptions:
52 : * - Assumes Cartesian coordinates with a diagonal Jacobian matrix
53 : * from the logical to the inertial frame
54 : * - Assumes the mesh is not moving (grid and inertial frame are the same)
55 : */
56 : template <size_t Dim>
57 1 : struct TimeDerivative {
58 : template <typename DbTagsList>
59 0 : static void apply(const gsl::not_null<db::DataBox<DbTagsList>*> box) {
60 : using metavariables =
61 : typename std::decay_t<decltype(db::get<Parallel::Tags::Metavariables>(
62 : *box))>;
63 : using system = typename metavariables::system;
64 : using evolved_vars_tag = typename system::variables_tag;
65 : using evolved_vars_tags = typename evolved_vars_tag::tags_list;
66 : using prim_tags = typename system::primitive_variables_tag::tags_list;
67 : using fluxes_tags = db::wrap_tags_in<::Tags::Flux, evolved_vars_tags,
68 : tmpl::size_t<Dim>, Frame::Inertial>;
69 :
70 : ASSERT((db::get<::domain::CoordinateMaps::Tags::CoordinateMap<
71 : Dim, Frame::Grid, Frame::Inertial>>(*box))
72 : .is_identity(),
73 : "Do not yet support moving mesh with DG-subcell.");
74 :
75 : // The copy of Mesh is intentional to avoid a GCC-7 internal compiler error.
76 : const Mesh<Dim> subcell_mesh =
77 : db::get<evolution::dg::subcell::Tags::Mesh<Dim>>(*box);
78 : ASSERT(
79 : subcell_mesh == Mesh<Dim>(subcell_mesh.extents(0),
80 : subcell_mesh.basis(0),
81 : subcell_mesh.quadrature(0)),
82 : "The subcell/FD mesh must be isotropic for the FD time derivative but "
83 : "got "
84 : << subcell_mesh);
85 : const size_t num_pts = subcell_mesh.number_of_grid_points();
86 : const size_t reconstructed_num_pts =
87 : (subcell_mesh.extents(0) + 1) *
88 : subcell_mesh.extents().slice_away(0).product();
89 :
90 : const tnsr::I<DataVector, Dim, Frame::ElementLogical>&
91 : cell_centered_logical_coords =
92 : db::get<evolution::dg::subcell::Tags::Coordinates<
93 : Dim, Frame::ElementLogical>>(*box);
94 : std::array<double, Dim> one_over_delta_xi{};
95 : for (size_t i = 0; i < Dim; ++i) {
96 : // Note: assumes isotropic extents
97 : gsl::at(one_over_delta_xi, i) =
98 : 1.0 / (get<0>(cell_centered_logical_coords)[1] -
99 : get<0>(cell_centered_logical_coords)[0]);
100 : }
101 :
102 : const NewtonianMhd::fd::Reconstructor<Dim>& recons =
103 : db::get<NewtonianMhd::fd::Tags::Reconstructor<Dim>>(*box);
104 :
105 : const Element<Dim>& element = db::get<domain::Tags::Element<Dim>>(*box);
106 : ASSERT(element.external_boundaries().size() == 0,
107 : "Can't have external boundaries right now with subcell. ElementID "
108 : << element.id());
109 :
110 : // Now package the data and compute the correction
111 : const auto& boundary_correction =
112 : db::get<evolution::Tags::BoundaryCorrection>(*box);
113 : using derived_boundary_corrections =
114 : tmpl::at<typename metavariables::factory_creation::factory_classes,
115 : evolution::BoundaryCorrection>;
116 : std::array<Variables<evolved_vars_tags>, Dim> boundary_corrections{};
117 : tmpl::for_each<derived_boundary_corrections>([&boundary_correction,
118 : &reconstructed_num_pts,
119 : &recons, &box, &element,
120 : &subcell_mesh,
121 : &boundary_corrections](
122 : auto
123 : derived_correction_v) {
124 : using DerivedCorrection = tmpl::type_from<decltype(derived_correction_v)>;
125 : if (typeid(boundary_correction) == typeid(DerivedCorrection)) {
126 : using dg_package_data_temporary_tags =
127 : typename DerivedCorrection::dg_package_data_temporary_tags;
128 : using dg_package_data_argument_tags =
129 : tmpl::append<evolved_vars_tags, prim_tags, fluxes_tags,
130 : dg_package_data_temporary_tags>;
131 : // Computed prims and cons on face via reconstruction
132 : auto package_data_argvars_lower_face = make_array<Dim>(
133 : Variables<dg_package_data_argument_tags>(reconstructed_num_pts));
134 : auto package_data_argvars_upper_face = make_array<Dim>(
135 : Variables<dg_package_data_argument_tags>(reconstructed_num_pts));
136 :
137 : // Reconstruct data to the face
138 : call_with_dynamic_type<void, typename NewtonianMhd::fd::Reconstructor<
139 : Dim>::creatable_classes>(
140 : &recons,
141 : [&box, &package_data_argvars_lower_face,
142 : &package_data_argvars_upper_face](const auto& reconstructor) {
143 : db::apply<typename std::decay_t<
144 : decltype(*reconstructor)>::reconstruction_argument_tags>(
145 : [&package_data_argvars_lower_face,
146 : &package_data_argvars_upper_face,
147 : &reconstructor](const auto&... args) {
148 : reconstructor->reconstruct(
149 : make_not_null(&package_data_argvars_lower_face),
150 : make_not_null(&package_data_argvars_upper_face),
151 : args...);
152 : },
153 : *box);
154 : });
155 :
156 : using dg_package_field_tags =
157 : typename DerivedCorrection::dg_package_field_tags;
158 : // Allocated outside for loop to reduce allocations
159 : Variables<dg_package_field_tags> upper_packaged_data{
160 : reconstructed_num_pts};
161 : Variables<dg_package_field_tags> lower_packaged_data{
162 : reconstructed_num_pts};
163 :
164 : // Compute fluxes on faces
165 : for (size_t i = 0; i < Dim; ++i) {
166 : auto& vars_upper_face = gsl::at(package_data_argvars_upper_face, i);
167 : auto& vars_lower_face = gsl::at(package_data_argvars_lower_face, i);
168 : NewtonianMhd::subcell::compute_fluxes<Dim>(
169 : make_not_null(&vars_upper_face));
170 : NewtonianMhd::subcell::compute_fluxes<Dim>(
171 : make_not_null(&vars_lower_face));
172 :
173 : tnsr::i<DataVector, Dim, Frame::Inertial> lower_outward_conormal{
174 : reconstructed_num_pts, 0.0};
175 : lower_outward_conormal.get(i) = 1.0;
176 :
177 : tnsr::i<DataVector, Dim, Frame::Inertial> upper_outward_conormal{
178 : reconstructed_num_pts, 0.0};
179 : upper_outward_conormal.get(i) = -1.0;
180 :
181 : // Compute the packaged data
182 : using dg_package_data_projected_tags = tmpl::append<
183 : evolved_vars_tags, fluxes_tags, dg_package_data_temporary_tags,
184 : typename DerivedCorrection::dg_package_data_primitive_tags>;
185 : evolution::dg::Actions::detail::dg_package_data<system>(
186 : make_not_null(&upper_packaged_data),
187 : dynamic_cast<const DerivedCorrection&>(boundary_correction),
188 : vars_upper_face, upper_outward_conormal, {std::nullopt}, *box,
189 : typename DerivedCorrection::dg_package_data_volume_tags{},
190 : dg_package_data_projected_tags{});
191 :
192 : evolution::dg::Actions::detail::dg_package_data<system>(
193 : make_not_null(&lower_packaged_data),
194 : dynamic_cast<const DerivedCorrection&>(boundary_correction),
195 : vars_lower_face, lower_outward_conormal, {std::nullopt}, *box,
196 : typename DerivedCorrection::dg_package_data_volume_tags{},
197 : dg_package_data_projected_tags{});
198 :
199 : // Now need to check if any of our neighbors are doing DG,
200 : // because if so then we need to use whatever boundary data
201 : // they sent instead of what we computed locally.
202 : //
203 : // Note: We could check this beforehand to avoid the extra
204 : // work of reconstruction and flux computations at the
205 : // boundaries.
206 : evolution::dg::subcell::correct_package_data<true>(
207 : make_not_null(&lower_packaged_data),
208 : make_not_null(&upper_packaged_data), i, element, subcell_mesh,
209 : db::get<evolution::dg::Tags::MortarData<Dim>>(*box), 0);
210 :
211 : // Compute the corrections on the faces. We only need to
212 : // compute this once because we can just flip the normal
213 : // vectors then
214 : gsl::at(boundary_corrections, i).initialize(reconstructed_num_pts);
215 : evolution::dg::subcell::compute_boundary_terms(
216 : make_not_null(&gsl::at(boundary_corrections, i)),
217 : dynamic_cast<const DerivedCorrection&>(boundary_correction),
218 : upper_packaged_data, lower_packaged_data, db::as_access(*box),
219 : typename DerivedCorrection::dg_boundary_terms_volume_tags{});
220 : }
221 : }
222 : });
223 :
224 : // Now compute the actual time derivatives.
225 : using dt_variables_tag = db::add_tag_prefix<::Tags::dt, evolved_vars_tag>;
226 : using source_argument_tags = tmpl::list<
227 : Tags::MassDensityCons, Tags::MomentumDensity<Dim>, Tags::EnergyDensity,
228 : hydro::Tags::SpatialVelocity<DataVector, Dim>,
229 : hydro::Tags::Pressure<DataVector>,
230 : hydro::Tags::SpecificInternalEnergy<DataVector>,
231 : hydro::Tags::EquationOfState<false, 2>,
232 : evolution::dg::subcell::Tags::Coordinates<Dim, Frame::Inertial>,
233 : ::Tags::Time, NewtonianMhd::Tags::SourceTerm<Dim>>;
234 : db::mutate_apply<tmpl::list<dt_variables_tag>, source_argument_tags>(
235 : [&num_pts, &boundary_corrections, &subcell_mesh, &one_over_delta_xi,
236 : &cell_centered_logical_to_grid_inv_jacobian =
237 : db::get<evolution::dg::subcell::fd::Tags::
238 : InverseJacobianLogicalToGrid<Dim>>(*box)](
239 : const auto dt_vars_ptr, const Scalar<DataVector>& mass_density_cons,
240 : const tnsr::I<DataVector, Dim>& momentum_density,
241 : const Scalar<DataVector>& energy_density,
242 : const tnsr::I<DataVector, Dim>& velocity,
243 : const Scalar<DataVector>& pressure,
244 : const Scalar<DataVector>& specific_internal_energy,
245 : const EquationsOfState::EquationOfState<false, 2>& eos,
246 : const tnsr::I<DataVector, Dim>& coords, const double time,
247 : const Sources::Source<Dim>& source) {
248 : dt_vars_ptr->initialize(num_pts, 0.0);
249 : using MassDensityCons = NewtonianMhd::Tags::MassDensityCons;
250 : using MomentumDensity = NewtonianMhd::Tags::MomentumDensity<Dim>;
251 : using EnergyDensity = NewtonianMhd::Tags::EnergyDensity;
252 :
253 : auto& dt_mass = get<::Tags::dt<MassDensityCons>>(*dt_vars_ptr);
254 : auto& dt_momentum = get<::Tags::dt<MomentumDensity>>(*dt_vars_ptr);
255 : auto& dt_energy = get<::Tags::dt<EnergyDensity>>(*dt_vars_ptr);
256 :
257 : const auto eos_2d = eos.promote_to_2d_eos();
258 : source(make_not_null(&dt_mass), make_not_null(&dt_momentum),
259 : make_not_null(&dt_energy), mass_density_cons, momentum_density,
260 : energy_density, velocity, pressure, specific_internal_energy,
261 : *eos_2d, coords, time);
262 :
263 : for (size_t dim = 0; dim < Dim; ++dim) {
264 : Scalar<DataVector>& mass_density_correction =
265 : get<MassDensityCons>(gsl::at(boundary_corrections, dim));
266 : Scalar<DataVector>& energy_density_correction =
267 : get<EnergyDensity>(gsl::at(boundary_corrections, dim));
268 : tnsr::I<DataVector, Dim, Frame::Inertial>&
269 : momentum_density_correction =
270 : get<MomentumDensity>(gsl::at(boundary_corrections, dim));
271 : evolution::dg::subcell::add_cartesian_flux_divergence(
272 : make_not_null(&get(dt_mass)), gsl::at(one_over_delta_xi, dim),
273 : cell_centered_logical_to_grid_inv_jacobian.get(dim, dim),
274 : get(mass_density_correction), subcell_mesh.extents(), dim);
275 : evolution::dg::subcell::add_cartesian_flux_divergence(
276 : make_not_null(&get(dt_energy)), gsl::at(one_over_delta_xi, dim),
277 : cell_centered_logical_to_grid_inv_jacobian.get(dim, dim),
278 : get(energy_density_correction), subcell_mesh.extents(), dim);
279 : for (size_t d = 0; d < Dim; ++d) {
280 : evolution::dg::subcell::add_cartesian_flux_divergence(
281 : make_not_null(&dt_momentum.get(d)),
282 : gsl::at(one_over_delta_xi, dim),
283 : cell_centered_logical_to_grid_inv_jacobian.get(dim, dim),
284 : momentum_density_correction.get(d), subcell_mesh.extents(),
285 : dim);
286 : }
287 : }
288 : },
289 : box);
290 : }
291 :
292 : private:
293 : template <typename SourceTerm, typename... SourceTermArgs,
294 : typename... SourcedVars>
295 0 : static void sources_impl(std::tuple<gsl::not_null<Scalar<DataVector>*>,
296 : gsl::not_null<tnsr::I<DataVector, Dim>*>,
297 : gsl::not_null<Scalar<DataVector>*>>
298 : dt_vars,
299 : tmpl::list<SourcedVars...> /*meta*/,
300 : const SourceTerm& source,
301 : const SourceTermArgs&... source_term_args) {
302 : using dt_vars_list =
303 : tmpl::list<Tags::MassDensityCons, Tags::MomentumDensity<Dim>,
304 : Tags::EnergyDensity>;
305 :
306 : source.apply(
307 : std::get<tmpl::index_of<dt_vars_list, SourcedVars>::value>(dt_vars)...,
308 : source_term_args...);
309 : }
310 : };
311 : } // namespace NewtonianMhd::subcell
|
