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/PackageDataImpl.hpp"
30 : #include "Evolution/DiscontinuousGalerkin/MortarTags.hpp"
31 : #include "Evolution/Systems/ScalarAdvection/FiniteDifference/Reconstructor.hpp"
32 : #include "Evolution/Systems/ScalarAdvection/FiniteDifference/Tags.hpp"
33 : #include "Evolution/Systems/ScalarAdvection/Fluxes.hpp"
34 : #include "Evolution/Systems/ScalarAdvection/Subcell/ComputeFluxes.hpp"
35 : #include "Evolution/Systems/ScalarAdvection/System.hpp"
36 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
37 : #include "Utilities/CallWithDynamicType.hpp"
38 : #include "Utilities/ErrorHandling/Assert.hpp"
39 : #include "Utilities/Gsl.hpp"
40 : #include "Utilities/MakeArray.hpp"
41 :
42 : namespace ScalarAdvection::subcell {
43 : /*!
44 : * \brief Compute the time derivative on the subcell grid using FD
45 : * reconstruction.
46 : *
47 : * The code makes the following unchecked assumptions:
48 : * - Assumes Cartesian coordinates with a diagonal Jacobian matrix
49 : * from the logical to the inertial frame
50 : * - Assumes the mesh is not moving (grid and inertial frame are the same)
51 : */
52 : template <size_t Dim>
53 1 : struct TimeDerivative {
54 : template <typename DbTagsList>
55 0 : static void apply(const gsl::not_null<db::DataBox<DbTagsList>*> box) {
56 : ASSERT((db::get<::domain::CoordinateMaps::Tags::CoordinateMap<
57 : Dim, Frame::Grid, Frame::Inertial>>(*box))
58 : .is_identity(),
59 : "Do not yet support moving mesh with DG-subcell.");
60 : // subcell is currently not supported for external boundary elements
61 : const Element<Dim>& element = db::get<domain::Tags::Element<Dim>>(*box);
62 : ASSERT(element.external_boundaries().size() == 0,
63 : "Can't have external boundaries right now with subcell. ElementID "
64 : << element.id());
65 :
66 : using evolved_vars_tags = typename System<Dim>::variables_tag::tags_list;
67 : using fluxes_tags = typename Fluxes<Dim>::return_tags;
68 :
69 : // The copy of Mesh is intentional to avoid a GCC-7 internal compiler error.
70 : const Mesh<Dim> subcell_mesh =
71 : db::get<evolution::dg::subcell::Tags::Mesh<Dim>>(*box);
72 : ASSERT(
73 : subcell_mesh == Mesh<Dim>(subcell_mesh.extents(0),
74 : subcell_mesh.basis(0),
75 : subcell_mesh.quadrature(0)),
76 : "The subcell/FD mesh must be isotropic for the FD time derivative but "
77 : "got "
78 : << subcell_mesh);
79 :
80 : const size_t num_reconstructed_pts =
81 : (subcell_mesh.extents(0) + 1) *
82 : subcell_mesh.extents().slice_away(0).product();
83 :
84 : const ScalarAdvection::fd::Reconstructor<Dim>& recons =
85 : db::get<ScalarAdvection::fd::Tags::Reconstructor<Dim>>(*box);
86 :
87 : const auto& boundary_correction =
88 : db::get<evolution::Tags::BoundaryCorrection<System<Dim>>>(*box);
89 : using derived_boundary_corrections =
90 : typename std::decay_t<decltype(boundary_correction)>::creatable_classes;
91 :
92 : // Variables to store the boundary correction terms on FD subinterfaces
93 : std::array<Variables<evolved_vars_tags>, Dim> fd_boundary_corrections{};
94 :
95 : // package the data and compute the boundary correction
96 : tmpl::for_each<derived_boundary_corrections>([&boundary_correction,
97 : &fd_boundary_corrections,
98 : &box, &element,
99 : &num_reconstructed_pts,
100 : &recons, &subcell_mesh](
101 : auto
102 : derived_correction_v) {
103 : using derived_correction =
104 : tmpl::type_from<decltype(derived_correction_v)>;
105 : if (typeid(boundary_correction) == typeid(derived_correction)) {
106 : using dg_package_field_tags =
107 : typename derived_correction::dg_package_field_tags;
108 : using dg_package_data_temporary_tags =
109 : typename derived_correction::dg_package_data_temporary_tags;
110 : using dg_package_data_argument_tags =
111 : tmpl::append<evolved_vars_tags, fluxes_tags,
112 : dg_package_data_temporary_tags>;
113 :
114 : // Variables that need to be reconstructed for dg_package_data()
115 : auto package_data_argvars_lower_face = make_array<Dim>(
116 : Variables<dg_package_data_argument_tags>(num_reconstructed_pts));
117 : auto package_data_argvars_upper_face = make_array<Dim>(
118 : Variables<dg_package_data_argument_tags>(num_reconstructed_pts));
119 :
120 : // Reconstruct the fields on interfaces
121 : call_with_dynamic_type<void, typename ScalarAdvection::fd::
122 : Reconstructor<Dim>::creatable_classes>(
123 : &recons,
124 : [&box, &package_data_argvars_lower_face,
125 : &package_data_argvars_upper_face](const auto& reconstructor) {
126 : db::apply<typename std::decay_t<decltype(
127 : *reconstructor)>::reconstruction_argument_tags>(
128 : [&package_data_argvars_lower_face,
129 : &package_data_argvars_upper_face,
130 : &reconstructor](const auto&... args) {
131 : reconstructor->reconstruct(
132 : make_not_null(&package_data_argvars_lower_face),
133 : make_not_null(&package_data_argvars_upper_face),
134 : args...);
135 : },
136 : *box);
137 : });
138 :
139 : // copy over the face values of the velocity field
140 : {
141 : using tag = ::ScalarAdvection::Tags::VelocityField<Dim>;
142 : const auto& velocity_on_face =
143 : db::get<evolution::dg::subcell::Tags::OnSubcellFaces<tag, Dim>>(
144 : *box);
145 : for (size_t d = 0; d < Dim; ++d) {
146 : get<tag>(gsl::at(package_data_argvars_lower_face, d)) =
147 : gsl::at(velocity_on_face, d);
148 : get<tag>(gsl::at(package_data_argvars_upper_face, d)) =
149 : gsl::at(velocity_on_face, d);
150 : }
151 : }
152 :
153 : // Variables to store packaged data. Allocate outside of loop to
154 : // reduce allocations
155 : Variables<dg_package_field_tags> upper_packaged_data{
156 : num_reconstructed_pts};
157 : Variables<dg_package_field_tags> lower_packaged_data{
158 : num_reconstructed_pts};
159 :
160 : for (size_t dim = 0; dim < Dim; ++dim) {
161 : auto& vars_upper_face = gsl::at(package_data_argvars_upper_face, dim);
162 : auto& vars_lower_face = gsl::at(package_data_argvars_lower_face, dim);
163 :
164 : // Compute fluxes on each faces
165 : ScalarAdvection::subcell::compute_fluxes<Dim>(
166 : make_not_null(&vars_upper_face));
167 : ScalarAdvection::subcell::compute_fluxes<Dim>(
168 : make_not_null(&vars_lower_face));
169 :
170 : // Note that we use the sign convention on the normal vectors to
171 : // be compatible with DG.
172 : tnsr::i<DataVector, Dim, Frame::Inertial> lower_outward_conormal{
173 : num_reconstructed_pts, 0.0};
174 : lower_outward_conormal.get(dim) = 1.0;
175 : tnsr::i<DataVector, Dim, Frame::Inertial> upper_outward_conormal{
176 : num_reconstructed_pts, 0.0};
177 : upper_outward_conormal.get(dim) = -1.0;
178 :
179 : // Compute the packaged data
180 : evolution::dg::Actions::detail::dg_package_data<System<Dim>>(
181 : make_not_null(&upper_packaged_data),
182 : dynamic_cast<const derived_correction&>(boundary_correction),
183 : vars_upper_face, upper_outward_conormal, {std::nullopt}, *box,
184 : typename derived_correction::dg_package_data_volume_tags{},
185 : dg_package_data_argument_tags{});
186 : evolution::dg::Actions::detail::dg_package_data<System<Dim>>(
187 : make_not_null(&lower_packaged_data),
188 : dynamic_cast<const derived_correction&>(boundary_correction),
189 : vars_lower_face, lower_outward_conormal, {std::nullopt}, *box,
190 : typename derived_correction::dg_package_data_volume_tags{},
191 : dg_package_data_argument_tags{});
192 :
193 : // Now need to check if any of our neighbors are doing DG, because
194 : // if so then we need to use whatever boundary data they sent
195 : // instead of what we computed locally.
196 : //
197 : // Note: We could check this beforehand to avoid the extra work of
198 : // reconstruction and flux computations at the boundaries.
199 : evolution::dg::subcell::correct_package_data<true>(
200 : make_not_null(&lower_packaged_data),
201 : make_not_null(&upper_packaged_data), dim, element, subcell_mesh,
202 : db::get<evolution::dg::Tags::MortarData<Dim>>(*box), 0);
203 :
204 : // Compute the corrections on the faces. We only need to compute this
205 : // once because we can just flip the normal vectors then
206 : gsl::at(fd_boundary_corrections, dim)
207 : .initialize(num_reconstructed_pts);
208 : evolution::dg::subcell::compute_boundary_terms(
209 : make_not_null(&gsl::at(fd_boundary_corrections, dim)),
210 : dynamic_cast<const derived_correction&>(boundary_correction),
211 : upper_packaged_data, lower_packaged_data, db::as_access(*box),
212 : typename derived_correction::dg_boundary_terms_volume_tags{});
213 : }
214 : }
215 : });
216 :
217 : std::array<double, Dim> one_over_delta_xi{};
218 : {
219 : const tnsr::I<DataVector, Dim, Frame::ElementLogical>&
220 : cell_centered_logical_coords =
221 : db::get<evolution::dg::subcell::Tags::Coordinates<
222 : Dim, Frame::ElementLogical>>(*box);
223 :
224 : for (size_t i = 0; i < Dim; ++i) {
225 : // Note: assumes isotropic extents
226 : gsl::at(one_over_delta_xi, i) =
227 : 1.0 / (get<0>(cell_centered_logical_coords)[1] -
228 : get<0>(cell_centered_logical_coords)[0]);
229 : }
230 : }
231 :
232 : // Now compute the actual time derivative
233 : using dt_variables_tag =
234 : db::add_tag_prefix<::Tags::dt, typename System<Dim>::variables_tag>;
235 : const size_t num_pts = subcell_mesh.number_of_grid_points();
236 : db::mutate<dt_variables_tag>(
237 : [&num_pts, &fd_boundary_corrections, &subcell_mesh, &one_over_delta_xi](
238 : const auto dt_vars_ptr,
239 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
240 : Frame::Grid>&
241 : cell_centered_logical_to_grid_inv_jacobian) {
242 : dt_vars_ptr->initialize(num_pts, 0.0);
243 : auto& dt_u =
244 : get<::Tags::dt<::ScalarAdvection::Tags::U>>(*dt_vars_ptr);
245 :
246 : for (size_t dim = 0; dim < Dim; ++dim) {
247 : Scalar<DataVector>& u_correction = get<::ScalarAdvection ::Tags::U>(
248 : gsl::at(fd_boundary_corrections, dim));
249 : evolution::dg::subcell::add_cartesian_flux_divergence(
250 : make_not_null(&get(dt_u)), gsl::at(one_over_delta_xi, dim),
251 : cell_centered_logical_to_grid_inv_jacobian.get(dim, dim),
252 : get(u_correction), subcell_mesh.extents(), dim);
253 : }
254 : },
255 : box,
256 : db::get<evolution::dg::subcell::fd::Tags::InverseJacobianLogicalToGrid<
257 : Dim>>(*box));
258 : }
259 : };
260 : } // namespace ScalarAdvection::subcell
|
