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/PrefixHelpers.hpp"
14 : #include "DataStructures/DataBox/Prefixes.hpp"
15 : #include "DataStructures/DataVector.hpp"
16 : #include "DataStructures/Tensor/Tensor.hpp"
17 : #include "DataStructures/Variables.hpp"
18 : #include "Domain/CoordinateMaps/Tags.hpp"
19 : #include "Domain/Structure/Element.hpp"
20 : #include "Domain/Tags.hpp"
21 : #include "Evolution/BoundaryCorrectionTags.hpp"
22 : #include "Evolution/DgSubcell/CartesianFluxDivergence.hpp"
23 : #include "Evolution/DgSubcell/ComputeBoundaryTerms.hpp"
24 : #include "Evolution/DgSubcell/CorrectPackagedData.hpp"
25 : #include "Evolution/DgSubcell/Tags/Coordinates.hpp"
26 : #include "Evolution/DgSubcell/Tags/Jacobians.hpp"
27 : #include "Evolution/DgSubcell/Tags/Mesh.hpp"
28 : #include "Evolution/DgSubcell/Tags/OnSubcellFaces.hpp"
29 : #include "Evolution/DiscontinuousGalerkin/Actions/NormalCovectorAndMagnitude.hpp"
30 : #include "Evolution/DiscontinuousGalerkin/Actions/PackageDataImpl.hpp"
31 : #include "Evolution/DiscontinuousGalerkin/MortarTags.hpp"
32 : #include "Evolution/Systems/NewtonianEuler/BoundaryCorrections/BoundaryCorrection.hpp"
33 : #include "Evolution/Systems/NewtonianEuler/FiniteDifference/Reconstructor.hpp"
34 : #include "Evolution/Systems/NewtonianEuler/FiniteDifference/Tag.hpp"
35 : #include "Evolution/Systems/NewtonianEuler/Fluxes.hpp"
36 : #include "Evolution/Systems/NewtonianEuler/Subcell/ComputeFluxes.hpp"
37 : #include "Evolution/Systems/NewtonianEuler/System.hpp"
38 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
39 : #include "Parallel/Tags/Metavariables.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<system>>(*box);
112 : using derived_boundary_corrections =
113 : typename std::decay_t<decltype(boundary_correction)>::creatable_classes;
114 : std::array<Variables<evolved_vars_tags>, Dim> boundary_corrections{};
115 : tmpl::for_each<derived_boundary_corrections>([&boundary_correction,
116 : &reconstructed_num_pts,
117 : &recons, &box, &element,
118 : &subcell_mesh,
119 : &boundary_corrections](
120 : auto
121 : derived_correction_v) {
122 : using DerivedCorrection = tmpl::type_from<decltype(derived_correction_v)>;
123 : if (typeid(boundary_correction) == typeid(DerivedCorrection)) {
124 : using dg_package_data_temporary_tags =
125 : typename DerivedCorrection::dg_package_data_temporary_tags;
126 : using dg_package_data_argument_tags =
127 : tmpl::append<evolved_vars_tags, prim_tags, fluxes_tags,
128 : dg_package_data_temporary_tags>;
129 : // Computed prims and cons on face via reconstruction
130 : auto package_data_argvars_lower_face = make_array<Dim>(
131 : Variables<dg_package_data_argument_tags>(reconstructed_num_pts));
132 : auto package_data_argvars_upper_face = make_array<Dim>(
133 : Variables<dg_package_data_argument_tags>(reconstructed_num_pts));
134 :
135 : // Reconstruct data to the face
136 : call_with_dynamic_type<void, typename NewtonianEuler::fd::Reconstructor<
137 : Dim>::creatable_classes>(
138 : &recons,
139 : [&box, &package_data_argvars_lower_face,
140 : &package_data_argvars_upper_face](const auto& reconstructor) {
141 : db::apply<typename std::decay_t<decltype(
142 : *reconstructor)>::reconstruction_argument_tags>(
143 : [&package_data_argvars_lower_face,
144 : &package_data_argvars_upper_face,
145 : &reconstructor](const auto&... args) {
146 : reconstructor->reconstruct(
147 : make_not_null(&package_data_argvars_lower_face),
148 : make_not_null(&package_data_argvars_upper_face),
149 : args...);
150 : },
151 : *box);
152 : });
153 :
154 : using dg_package_field_tags =
155 : typename DerivedCorrection::dg_package_field_tags;
156 : // Allocated outside for loop to reduce allocations
157 : Variables<dg_package_field_tags> upper_packaged_data{
158 : reconstructed_num_pts};
159 : Variables<dg_package_field_tags> lower_packaged_data{
160 : reconstructed_num_pts};
161 :
162 : // Compute fluxes on faces
163 : for (size_t i = 0; i < Dim; ++i) {
164 : auto& vars_upper_face = gsl::at(package_data_argvars_upper_face, i);
165 : auto& vars_lower_face = gsl::at(package_data_argvars_lower_face, i);
166 : NewtonianEuler::subcell::compute_fluxes<Dim>(
167 : make_not_null(&vars_upper_face));
168 : NewtonianEuler::subcell::compute_fluxes<Dim>(
169 : make_not_null(&vars_lower_face));
170 :
171 : tnsr::i<DataVector, Dim, Frame::Inertial> lower_outward_conormal{
172 : reconstructed_num_pts, 0.0};
173 : lower_outward_conormal.get(i) = 1.0;
174 :
175 : tnsr::i<DataVector, Dim, Frame::Inertial> upper_outward_conormal{
176 : reconstructed_num_pts, 0.0};
177 : upper_outward_conormal.get(i) = -1.0;
178 :
179 : // Compute the packaged data
180 : using dg_package_data_projected_tags = tmpl::append<
181 : evolved_vars_tags, fluxes_tags, dg_package_data_temporary_tags,
182 : typename DerivedCorrection::dg_package_data_primitive_tags>;
183 : evolution::dg::Actions::detail::dg_package_data<system>(
184 : make_not_null(&upper_packaged_data),
185 : dynamic_cast<const DerivedCorrection&>(boundary_correction),
186 : vars_upper_face, upper_outward_conormal, {std::nullopt}, *box,
187 : typename DerivedCorrection::dg_package_data_volume_tags{},
188 : dg_package_data_projected_tags{});
189 :
190 : evolution::dg::Actions::detail::dg_package_data<system>(
191 : make_not_null(&lower_packaged_data),
192 : dynamic_cast<const DerivedCorrection&>(boundary_correction),
193 : vars_lower_face, lower_outward_conormal, {std::nullopt}, *box,
194 : typename DerivedCorrection::dg_package_data_volume_tags{},
195 : dg_package_data_projected_tags{});
196 :
197 : // Now need to check if any of our neighbors are doing DG,
198 : // because if so then we need to use whatever boundary data
199 : // they sent instead of what we computed locally.
200 : //
201 : // Note: We could check this beforehand to avoid the extra
202 : // work of reconstruction and flux computations at the
203 : // boundaries.
204 : evolution::dg::subcell::correct_package_data<true>(
205 : make_not_null(&lower_packaged_data),
206 : make_not_null(&upper_packaged_data), i, element, subcell_mesh,
207 : db::get<evolution::dg::Tags::MortarData<Dim>>(*box), 0);
208 :
209 : // Compute the corrections on the faces. We only need to
210 : // compute this once because we can just flip the normal
211 : // vectors then
212 : gsl::at(boundary_corrections, i).initialize(reconstructed_num_pts);
213 : evolution::dg::subcell::compute_boundary_terms(
214 : make_not_null(&gsl::at(boundary_corrections, i)),
215 : dynamic_cast<const DerivedCorrection&>(boundary_correction),
216 : upper_packaged_data, lower_packaged_data, db::as_access(*box),
217 : typename DerivedCorrection::dg_boundary_terms_volume_tags{});
218 : }
219 : }
220 : });
221 :
222 : // Now compute the actual time derivatives.
223 : using dt_variables_tag = db::add_tag_prefix<::Tags::dt, evolved_vars_tag>;
224 : using source_argument_tags = tmpl::list<
225 : Tags::MassDensityCons, Tags::MomentumDensity<Dim>, Tags::EnergyDensity,
226 : hydro::Tags::SpatialVelocity<DataVector, Dim>,
227 : hydro::Tags::Pressure<DataVector>,
228 : hydro::Tags::SpecificInternalEnergy<DataVector>,
229 : hydro::Tags::EquationOfState<false, 2>,
230 : evolution::dg::subcell::Tags::Coordinates<Dim, Frame::Inertial>,
231 : ::Tags::Time, NewtonianEuler::Tags::SourceTerm<Dim>>;
232 : db::mutate_apply<tmpl::list<dt_variables_tag>, source_argument_tags>(
233 : [&num_pts, &boundary_corrections, &subcell_mesh, &one_over_delta_xi,
234 : &cell_centered_logical_to_grid_inv_jacobian =
235 : db::get<evolution::dg::subcell::fd::Tags::
236 : InverseJacobianLogicalToGrid<Dim>>(*box)](
237 : const auto dt_vars_ptr, const Scalar<DataVector>& mass_density_cons,
238 : const tnsr::I<DataVector, Dim>& momentum_density,
239 : const Scalar<DataVector>& energy_density,
240 : const tnsr::I<DataVector, Dim>& velocity,
241 : const Scalar<DataVector>& pressure,
242 : const Scalar<DataVector>& specific_internal_energy,
243 : const EquationsOfState::EquationOfState<false, 2>& eos,
244 : const tnsr::I<DataVector, Dim>& coords, const double time,
245 : const Sources::Source<Dim>& source) {
246 : dt_vars_ptr->initialize(num_pts, 0.0);
247 : using MassDensityCons = NewtonianEuler::Tags::MassDensityCons;
248 : using MomentumDensity = NewtonianEuler::Tags::MomentumDensity<Dim>;
249 : using EnergyDensity = NewtonianEuler::Tags::EnergyDensity;
250 :
251 : auto& dt_mass = get<::Tags::dt<MassDensityCons>>(*dt_vars_ptr);
252 : auto& dt_momentum = get<::Tags::dt<MomentumDensity>>(*dt_vars_ptr);
253 : auto& dt_energy = get<::Tags::dt<EnergyDensity>>(*dt_vars_ptr);
254 :
255 : const auto eos_2d = eos.promote_to_2d_eos();
256 : source(make_not_null(&dt_mass), make_not_null(&dt_momentum),
257 : make_not_null(&dt_energy), mass_density_cons, momentum_density,
258 : energy_density, velocity, pressure, specific_internal_energy,
259 : *eos_2d, coords, time);
260 :
261 : for (size_t dim = 0; dim < Dim; ++dim) {
262 : Scalar<DataVector>& mass_density_correction =
263 : get<MassDensityCons>(gsl::at(boundary_corrections, dim));
264 : Scalar<DataVector>& energy_density_correction =
265 : get<EnergyDensity>(gsl::at(boundary_corrections, dim));
266 : tnsr::I<DataVector, Dim, Frame::Inertial>&
267 : momentum_density_correction =
268 : get<MomentumDensity>(gsl::at(boundary_corrections, dim));
269 : evolution::dg::subcell::add_cartesian_flux_divergence(
270 : make_not_null(&get(dt_mass)), gsl::at(one_over_delta_xi, dim),
271 : cell_centered_logical_to_grid_inv_jacobian.get(dim, dim),
272 : get(mass_density_correction), subcell_mesh.extents(), dim);
273 : evolution::dg::subcell::add_cartesian_flux_divergence(
274 : make_not_null(&get(dt_energy)), gsl::at(one_over_delta_xi, dim),
275 : cell_centered_logical_to_grid_inv_jacobian.get(dim, dim),
276 : get(energy_density_correction), subcell_mesh.extents(), dim);
277 : for (size_t d = 0; d < Dim; ++d) {
278 : evolution::dg::subcell::add_cartesian_flux_divergence(
279 : make_not_null(&dt_momentum.get(d)),
280 : gsl::at(one_over_delta_xi, dim),
281 : cell_centered_logical_to_grid_inv_jacobian.get(dim, dim),
282 : momentum_density_correction.get(d), subcell_mesh.extents(),
283 : dim);
284 : }
285 : }
286 : },
287 : box);
288 : }
289 :
290 : private:
291 : template <typename SourceTerm, typename... SourceTermArgs,
292 : typename... SourcedVars>
293 0 : static void sources_impl(std::tuple<gsl::not_null<Scalar<DataVector>*>,
294 : gsl::not_null<tnsr::I<DataVector, Dim>*>,
295 : gsl::not_null<Scalar<DataVector>*>>
296 : dt_vars,
297 : tmpl::list<SourcedVars...> /*meta*/,
298 : const SourceTerm& source,
299 : const SourceTermArgs&... source_term_args) {
300 : using dt_vars_list =
301 : tmpl::list<Tags::MassDensityCons, Tags::MomentumDensity<Dim>,
302 : Tags::EnergyDensity>;
303 :
304 : source.apply(
305 : std::get<tmpl::index_of<dt_vars_list, SourcedVars>::value>(dt_vars)...,
306 : source_term_args...);
307 : }
308 : };
309 : } // namespace NewtonianEuler::subcell
|