VerifyGrMhdSolution.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include "DataStructures/DataVector.hpp"
10 #include "Domain/Block.hpp"
13 #include "Domain/Mesh.hpp"
14 #include "Evolution/Systems/GrMhd/ValenciaDivClean/ConservativeFromPrimitive.hpp"
15 #include "Evolution/Systems/GrMhd/ValenciaDivClean/Fluxes.hpp"
16 #include "Evolution/Systems/GrMhd/ValenciaDivClean/Sources.hpp"
17 #include "Evolution/Systems/GrMhd/ValenciaDivClean/System.hpp"
18 #include "NumericalAlgorithms/LinearOperators/Divergence.tpp"
19 #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.tpp"
20 #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
21 #include "PointwiseFunctions/Hydro/Tags.hpp"
22 #include "Utilities/Gsl.hpp"
23 #include "Utilities/MakeArray.hpp"
25 #include "Utilities/TMPL.hpp"
26 #include "Utilities/TypeTraits.hpp"
27 
29 
30 using valencia_tags = tmpl::list<grmhd::ValenciaDivClean::Tags::TildeD,
35 
36 // compute the time derivative using a centered sixth-order stencil
37 template <typename Solution>
38 Variables<valencia_tags> numerical_dt(
39  const Solution& solution, const tnsr::I<DataVector, 3, Frame::Inertial>& x,
40  const double time, const double delta_time) {
41  std::array<double, 6> six_times{
42  {time - 3.0 * delta_time, time - 2.0 * delta_time, time - delta_time,
43  time + delta_time, time + 2.0 * delta_time, time + 3.0 * delta_time}};
44  const size_t number_of_points = get<0>(x).size();
45  auto solution_at_six_times =
46  make_array<6>(Variables<valencia_tags>(number_of_points));
47 
48  using solution_tags =
49  tmpl::list<hydro::Tags::RestMassDensity<DataVector>,
59 
60  for (size_t i = 0; i < 6; ++i) {
61  const auto vars =
62  solution.variables(x, gsl::at(six_times, i), solution_tags{});
63 
64  const auto& rest_mass_density =
65  get<hydro::Tags::RestMassDensity<DataVector>>(vars);
66  const auto& specific_internal_energy =
67  get<hydro::Tags::SpecificInternalEnergy<DataVector>>(vars);
68  const auto& spatial_velocity =
69  get<hydro::Tags::SpatialVelocity<DataVector, 3>>(vars);
70  const auto& magnetic_field =
71  get<hydro::Tags::MagneticField<DataVector, 3>>(vars);
72  const auto& divergence_cleaning_field =
73  get<hydro::Tags::DivergenceCleaningField<DataVector>>(vars);
74  const auto& lorentz_factor =
75  get<hydro::Tags::LorentzFactor<DataVector>>(vars);
76  const auto& pressure = get<hydro::Tags::Pressure<DataVector>>(vars);
77  const auto& specific_enthalpy =
78  get<hydro::Tags::SpecificEnthalpy<DataVector>>(vars);
79  const auto& spatial_metric =
80  get<gr::Tags::SpatialMetric<3, Frame::Inertial, DataVector>>(vars);
81  const auto& sqrt_det_spatial_metric =
82  get<gr::Tags::SqrtDetSpatialMetric<DataVector>>(vars);
83 
84  grmhd::ValenciaDivClean::ConservativeFromPrimitive::apply(
85  make_not_null(&get<grmhd::ValenciaDivClean::Tags::TildeD>(
86  gsl::at(solution_at_six_times, i))),
87  make_not_null(&get<grmhd::ValenciaDivClean::Tags::TildeTau>(
88  gsl::at(solution_at_six_times, i))),
89  make_not_null(&get<grmhd::ValenciaDivClean::Tags::TildeS<>>(
90  gsl::at(solution_at_six_times, i))),
91  make_not_null(&get<grmhd::ValenciaDivClean::Tags::TildeB<>>(
92  gsl::at(solution_at_six_times, i))),
93  make_not_null(&get<grmhd::ValenciaDivClean::Tags::TildePhi>(
94  gsl::at(solution_at_six_times, i))),
95  rest_mass_density, specific_internal_energy, specific_enthalpy,
96  pressure, spatial_velocity, lorentz_factor, magnetic_field,
97  sqrt_det_spatial_metric, spatial_metric, divergence_cleaning_field);
98  }
99 
100  return (-1.0 / (60.0 * delta_time)) * solution_at_six_times[0] +
101  (3.0 / (20.0 * delta_time)) * solution_at_six_times[1] +
102  (-0.75 / delta_time) * solution_at_six_times[2] +
103  (0.75 / delta_time) * solution_at_six_times[3] +
104  (-3.0 / (20.0 * delta_time)) * solution_at_six_times[4] +
105  (1.0 / (60.0 * delta_time)) * solution_at_six_times[5];
106 }
107 } // namespace VerifyGrMhdSolution_detail
108 
109 /// \ingroup TestingFrameworkGroup
110 /// \brief Determines if the given `solution` is a solution of the GRMHD
111 /// equations.
112 ///
113 /// Uses numerical derivatives to compute the solution, on the given `mesh` of
114 /// the root Element of the given `block` at the given `time` using a
115 /// sixth-order derivative in time for the given `delta_time`. The maximum
116 /// residual of the GRMHD equations must be zero within `error_tolerance`
117 template <typename Solution>
118 void verify_grmhd_solution(const Solution& solution,
119  const Block<3, Frame::Inertial>& block,
120  const Mesh<3>& mesh, const double error_tolerance,
121  const double time,
122  const double delta_time) noexcept {
123  // Set up coordinates
124  const auto x_logical = logical_coordinates(mesh);
125  const auto x = block.coordinate_map()(x_logical);
126 
127  // Evaluate analytic solution
128  using solution_tags = tmpl::list<
130  hydro::Tags::SpecificInternalEnergy<DataVector>,
131  hydro::Tags::SpatialVelocity<DataVector, 3, Frame::Inertial>,
132  hydro::Tags::MagneticField<DataVector, 3, Frame::Inertial>,
133  hydro::Tags::DivergenceCleaningField<DataVector>,
134  hydro::Tags::LorentzFactor<DataVector>, hydro::Tags::Pressure<DataVector>,
135  hydro::Tags::SpecificEnthalpy<DataVector>, gr::Tags::Lapse<DataVector>,
137  gr::Tags::SpatialMetric<3, Frame::Inertial, DataVector>,
139  gr::Tags::SqrtDetSpatialMetric<DataVector>,
140  ::Tags::deriv<gr::Tags::Lapse<DataVector>, tmpl::size_t<3>,
142  ::Tags::deriv<gr::Tags::Shift<3, Frame::Inertial, DataVector>,
143  tmpl::size_t<3>, Frame::Inertial>,
144  ::Tags::deriv<gr::Tags::SpatialMetric<3, Frame::Inertial, DataVector>,
145  tmpl::size_t<3>, Frame::Inertial>,
147  const auto vars = solution.variables(x, time, solution_tags{});
148 
149  const auto& rest_mass_density =
150  get<hydro::Tags::RestMassDensity<DataVector>>(vars);
151  const auto& specific_internal_energy =
152  get<hydro::Tags::SpecificInternalEnergy<DataVector>>(vars);
153  const auto& spatial_velocity =
154  get<hydro::Tags::SpatialVelocity<DataVector, 3>>(vars);
155  const auto& magnetic_field =
156  get<hydro::Tags::MagneticField<DataVector, 3>>(vars);
157  const auto& divergence_cleaning_field =
158  get<hydro::Tags::DivergenceCleaningField<DataVector>>(vars);
159  const auto& lorentz_factor =
160  get<hydro::Tags::LorentzFactor<DataVector>>(vars);
161  const auto& pressure = get<hydro::Tags::Pressure<DataVector>>(vars);
162  const auto& specific_enthalpy =
163  get<hydro::Tags::SpecificEnthalpy<DataVector>>(vars);
164  const auto& lapse = get<gr::Tags::Lapse<DataVector>>(vars);
165  const auto& d_lapse =
166  get<::Tags::deriv<gr::Tags::Lapse<DataVector>, tmpl::size_t<3>,
167  Frame::Inertial>>(vars);
168  const auto& shift =
169  get<gr::Tags::Shift<3, Frame::Inertial, DataVector>>(vars);
170  const auto& d_shift =
171  get<::Tags::deriv<gr::Tags::Shift<3, Frame::Inertial, DataVector>,
172  tmpl::size_t<3>, Frame::Inertial>>(vars);
173  const auto& spatial_metric =
174  get<gr::Tags::SpatialMetric<3, Frame::Inertial, DataVector>>(vars);
175  const auto& d_spatial_metric =
176  get<::Tags::deriv<gr::Tags::SpatialMetric<3, Frame::Inertial, DataVector>,
177  tmpl::size_t<3>, Frame::Inertial>>(vars);
178  const auto& inv_spatial_metric =
179  get<gr::Tags::InverseSpatialMetric<3, Frame::Inertial, DataVector>>(vars);
180  const auto& sqrt_det_spatial_metric =
181  get<gr::Tags::SqrtDetSpatialMetric<DataVector>>(vars);
182  const auto& extrinsic_curvature =
183  get<gr::Tags::ExtrinsicCurvature<3, Frame::Inertial, DataVector>>(vars);
184 
185  const size_t number_of_points = mesh.number_of_grid_points();
186  Scalar<DataVector> tilde_d(number_of_points);
187  Scalar<DataVector> tilde_tau(number_of_points);
188  tnsr::i<DataVector, 3> tilde_s(number_of_points);
189  tnsr::I<DataVector, 3> tilde_b(number_of_points);
190  Scalar<DataVector> tilde_phi(number_of_points);
191 
192  grmhd::ValenciaDivClean::ConservativeFromPrimitive::apply(
193  make_not_null(&tilde_d), make_not_null(&tilde_tau),
194  make_not_null(&tilde_s), make_not_null(&tilde_b),
195  make_not_null(&tilde_phi), rest_mass_density, specific_internal_energy,
196  specific_enthalpy, pressure, spatial_velocity, lorentz_factor,
197  magnetic_field, sqrt_det_spatial_metric, spatial_metric,
198  divergence_cleaning_field);
199 
200  using flux_tags =
201  tmpl::list<Tags::Flux<grmhd::ValenciaDivClean::Tags::TildeD,
202  tmpl::size_t<3>, Frame::Inertial>,
203  Tags::Flux<grmhd::ValenciaDivClean::Tags::TildeTau,
204  tmpl::size_t<3>, Frame::Inertial>,
206  tmpl::size_t<3>, Frame::Inertial>,
208  tmpl::size_t<3>, Frame::Inertial>,
209  Tags::Flux<grmhd::ValenciaDivClean::Tags::TildePhi,
210  tmpl::size_t<3>, Frame::Inertial>>;
211  Variables<flux_tags> fluxes(number_of_points);
212  auto& flux_tilde_d =
213  get<Tags::Flux<grmhd::ValenciaDivClean::Tags::TildeD, tmpl::size_t<3>,
214  Frame::Inertial>>(fluxes);
215  auto& flux_tilde_tau =
216  get<Tags::Flux<grmhd::ValenciaDivClean::Tags::TildeTau, tmpl::size_t<3>,
217  Frame::Inertial>>(fluxes);
218  auto& flux_tilde_s =
219  get<Tags::Flux<grmhd::ValenciaDivClean::Tags::TildeS<>, tmpl::size_t<3>,
220  Frame::Inertial>>(fluxes);
221  auto& flux_tilde_b =
222  get<Tags::Flux<grmhd::ValenciaDivClean::Tags::TildeB<>, tmpl::size_t<3>,
223  Frame::Inertial>>(fluxes);
224  auto& flux_tilde_phi =
225  get<Tags::Flux<grmhd::ValenciaDivClean::Tags::TildePhi, tmpl::size_t<3>,
226  Frame::Inertial>>(fluxes);
227 
228  grmhd::ValenciaDivClean::ComputeFluxes::apply(
229  make_not_null(&flux_tilde_d), make_not_null(&flux_tilde_tau),
230  make_not_null(&flux_tilde_s), make_not_null(&flux_tilde_b),
231  make_not_null(&flux_tilde_phi), tilde_d, tilde_tau, tilde_s, tilde_b,
232  tilde_phi, lapse, shift, sqrt_det_spatial_metric, spatial_metric,
233  inv_spatial_metric, pressure, spatial_velocity, lorentz_factor,
234  magnetic_field);
235 
236  const auto div_of_fluxes = divergence<flux_tags, 3, Frame::Inertial>(
237  fluxes, mesh, block.coordinate_map().inv_jacobian(x_logical));
238 
239  const auto& div_flux_tilde_d =
240  get<Tags::div<Tags::Flux<grmhd::ValenciaDivClean::Tags::TildeD,
241  tmpl::size_t<3>, Frame::Inertial>>>(
242  div_of_fluxes);
243  const auto& div_flux_tilde_tau =
244  get<Tags::div<Tags::Flux<grmhd::ValenciaDivClean::Tags::TildeTau,
245  tmpl::size_t<3>, Frame::Inertial>>>(
246  div_of_fluxes);
247  const auto& div_flux_tilde_s =
248  get<Tags::div<Tags::Flux<grmhd::ValenciaDivClean::Tags::TildeS<>,
249  tmpl::size_t<3>, Frame::Inertial>>>(
250  div_of_fluxes);
251  const auto& div_flux_tilde_b =
252  get<Tags::div<Tags::Flux<grmhd::ValenciaDivClean::Tags::TildeB<>,
253  tmpl::size_t<3>, Frame::Inertial>>>(
254  div_of_fluxes);
255  const auto& div_flux_tilde_phi =
256  get<Tags::div<Tags::Flux<grmhd::ValenciaDivClean::Tags::TildePhi,
257  tmpl::size_t<3>, Frame::Inertial>>>(
258  div_of_fluxes);
259 
260  Scalar<DataVector> source_tilde_tau(number_of_points);
261  tnsr::i<DataVector, 3> source_tilde_s(number_of_points);
262  tnsr::I<DataVector, 3> source_tilde_b(number_of_points);
263  Scalar<DataVector> source_tilde_phi(number_of_points);
264  grmhd::ValenciaDivClean::ComputeSources::apply(
265  make_not_null(&source_tilde_tau), make_not_null(&source_tilde_s),
266  make_not_null(&source_tilde_b), make_not_null(&source_tilde_phi), tilde_d,
267  tilde_tau, tilde_s, tilde_b, tilde_phi, spatial_velocity, magnetic_field,
268  rest_mass_density, specific_enthalpy, lorentz_factor, pressure, lapse,
269  d_lapse, d_shift, spatial_metric, d_spatial_metric, inv_spatial_metric,
270  sqrt_det_spatial_metric, extrinsic_curvature, 0.0);
271 
272  Scalar<DataVector> residual_tilde_d(number_of_points, 0.);
273  Scalar<DataVector> residual_tilde_tau(number_of_points, 0.);
274  tnsr::i<DataVector, 3> residual_tilde_s(number_of_points, 0.);
275  tnsr::I<DataVector, 3> residual_tilde_b(number_of_points, 0.);
276  Scalar<DataVector> residual_tilde_phi(number_of_points, 0.);
277 
278  Variables<VerifyGrMhdSolution_detail::valencia_tags> dt_solution =
279  VerifyGrMhdSolution_detail::numerical_dt(solution, x, time, delta_time);
280 
281  get(residual_tilde_d) =
282  get(get<grmhd::ValenciaDivClean::Tags::TildeD>(dt_solution)) +
283  get(div_flux_tilde_d);
284  get(residual_tilde_tau) =
285  get(get<grmhd::ValenciaDivClean::Tags::TildeTau>(dt_solution)) +
286  get(div_flux_tilde_tau) - get(source_tilde_tau);
287  get<0>(residual_tilde_s) =
288  get<0>(get<grmhd::ValenciaDivClean::Tags::TildeS<>>(dt_solution)) +
289  get<0>(div_flux_tilde_s) - get<0>(source_tilde_s);
290  get<1>(residual_tilde_s) =
291  get<1>(get<grmhd::ValenciaDivClean::Tags::TildeS<>>(dt_solution)) +
292  get<1>(div_flux_tilde_s) - get<1>(source_tilde_s);
293  get<2>(residual_tilde_s) =
294  get<2>(get<grmhd::ValenciaDivClean::Tags::TildeS<>>(dt_solution)) +
295  get<2>(div_flux_tilde_s) - get<2>(source_tilde_s);
296  get<0>(residual_tilde_b) =
297  get<0>(get<grmhd::ValenciaDivClean::Tags::TildeB<>>(dt_solution)) +
298  get<0>(div_flux_tilde_b) - get<0>(source_tilde_b);
299  get<1>(residual_tilde_b) =
300  get<1>(get<grmhd::ValenciaDivClean::Tags::TildeB<>>(dt_solution)) +
301  get<1>(div_flux_tilde_b) - get<1>(source_tilde_b);
302  get<2>(residual_tilde_b) =
303  get<2>(get<grmhd::ValenciaDivClean::Tags::TildeB<>>(dt_solution)) +
304  get<2>(div_flux_tilde_b) - get<2>(source_tilde_b);
305  get(residual_tilde_phi) +=
306  get(get<grmhd::ValenciaDivClean::Tags::TildePhi>(dt_solution)) +
307  get(div_flux_tilde_phi) - get(source_tilde_phi);
308 
309  Approx numerical_approx =
310  Approx::custom().epsilon(error_tolerance).scale(1.0);
311  CHECK(max(abs(get(residual_tilde_d))) == numerical_approx(0.0));
312  CHECK(max(abs(get(residual_tilde_tau))) == numerical_approx(0.0));
313  CHECK(max(abs(get<0>(residual_tilde_s))) == numerical_approx(0.0));
314  CHECK(max(abs(get<1>(residual_tilde_s))) == numerical_approx(0.0));
315  CHECK(max(abs(get<2>(residual_tilde_s))) == numerical_approx(0.0));
316  CHECK(max(abs(get<0>(residual_tilde_b))) == numerical_approx(0.0));
317  CHECK(max(abs(get<1>(residual_tilde_b))) == numerical_approx(0.0));
318  CHECK(max(abs(get<2>(residual_tilde_b))) == numerical_approx(0.0));
319  CHECK(max(abs(get(residual_tilde_phi))) == numerical_approx(0.0));
320 }
void verify_grmhd_solution(const Solution &solution, const Block< 3, Frame::Inertial > &block, const Mesh< 3 > &mesh, const double error_tolerance, const double time, const double delta_time) noexcept
Determines if the given solution is a solution of the GRMHD equations.
Definition: VerifyGrMhdSolution.hpp:118
Scalar< DataType > lorentz_factor(const tnsr::I< DataType, Dim, Fr > &spatial_velocity, const tnsr::i< DataType, Dim, Fr > &spatial_velocity_form) noexcept
Computes the Lorentz factor .
The spatial velocity of the fluid, where . Here is the spatial part of the 4-velocity of the fluid...
Definition: Tags.hpp:156
The fluid pressure .
Definition: Tags.hpp:130
The specific internal energy .
Definition: Tags.hpp:189
Defines function make_array.
Definition: Tags.hpp:47
Holds the number of grid points, basis, and quadrature in each direction of the computational grid...
Definition: Mesh.hpp:49
The densitized divergence-cleaning field .
Definition: Tags.hpp:70
Definition: Tags.hpp:57
Inverse of the spatial metric.
Definition: Tags.hpp:34
tnsr::ii< DataType, SpatialDim, Frame > spatial_metric(const tnsr::aa< DataType, SpatialDim, Frame > &spacetime_metric) noexcept
Compute spatial metric from spacetime metric.
Definition: ComputeSpacetimeQuantities.cpp:62
Scalar< DataType > specific_enthalpy(const Scalar< DataType > &rest_mass_density, const Scalar< DataType > &specific_internal_energy, const Scalar< DataType > &pressure) noexcept
Computes the relativistic specific enthalpy as: where is the specific internal energy...
The magnetic field measured by an Eulerian observer, where is the normal to the spatial hypersurfac...
Definition: Tags.hpp:86
tnsr::ii< DataVector, 3, Frame > extrinsic_curvature(const tnsr::ii< DataVector, 3, Frame > &grad_normal, const tnsr::i< DataVector, 3, Frame > &unit_normal_one_form, const tnsr::I< DataVector, 3, Frame > &unit_normal_vector) noexcept
Extrinsic curvature of a 2D Strahlkorper embedded in a 3D space.
The divergence-cleaning field .
Definition: Tags.hpp:50
The densitized energy density .
Definition: Tags.hpp:50
The densitized momentum density .
Definition: Tags.hpp:57
Defines functions logical_coordinates and interface_logical_coordinates.
The densitized magnetic field .
Definition: Tags.hpp:64
Definition: Tags.hpp:26
Prefix indicating a flux.
Definition: Prefixes.hpp:38
Defines class Variables.
The Lorentz factor , where is the spatial velocity of the fluid.
Definition: Tags.hpp:68
Definition: Tags.hpp:52
Scalar< DataType > lapse(const tnsr::I< DataType, SpatialDim, Frame > &shift, const tnsr::aa< DataType, SpatialDim, Frame > &spacetime_metric) noexcept
Compute lapse from shift and spacetime metric.
Definition: ComputeSpacetimeQuantities.cpp:117
Defines the class template Mesh.
Defines class CoordinateMap.
Defines classes for Tensor.
Defines a list of useful type aliases for tensors.
tnsr::I< DataVector, VolumeDim, Frame::Logical > logical_coordinates(const Mesh< VolumeDim > &mesh) noexcept
Compute the logical coordinates in an Element.
Definition: LogicalCoordinates.cpp:20
T size(T... args)
tnsr::I< DataType, SpatialDim, Frame > shift(const tnsr::aa< DataType, SpatialDim, Frame > &spacetime_metric, const tnsr::II< DataType, SpatialDim, Frame > &inverse_spatial_metric) noexcept
Compute shift from spacetime metric and inverse spatial metric.
Definition: ComputeSpacetimeQuantities.cpp:101
The densitized rest-mass density .
Definition: Tags.hpp:44
Wraps the template metaprogramming library used (brigand)
Definition: VerifyGrMhdSolution.hpp:28
Defines functions and classes from the GSL.
gsl::not_null< T * > make_not_null(T *ptr) noexcept
Construct a not_null from a pointer. Often this will be done as an implicit conversion, but it may be necessary to perform the conversion explicitly when type deduction is desired.
Definition: Gsl.hpp:863
The rest-mass density .
Definition: Tags.hpp:137
Definition: IndexType.hpp:44
Defines type traits, some of which are future STL type_traits header.
The specific enthalpy .
Definition: Tags.hpp:182
Defines class template Block.
Tensor< T, Symmetry<>, index_list<> > Scalar
Scalar type.
Definition: TypeAliases.hpp:21
A Block<VolumeDim> is a region of a VolumeDim-dimensional computational domain that defines the root ...
Definition: Block.hpp:42
Defines make_with_value.
constexpr T & at(std::array< T, N > &arr, Size index)
Retrieve a entry from a container, with checks in Debug mode that the index being retrieved is valid...
Definition: Gsl.hpp:124
Definition: Tags.hpp:146