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