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