FirstOrderEllipticSolutionsTestHelpers.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <cstddef>
7 #include <limits>
8 
9 #include "DataStructures/DataBox/PrefixHelpers.hpp"
10 #include "DataStructures/DataBox/Tag.hpp"
11 #include "DataStructures/DataBox/TagName.hpp"
12 #include "DataStructures/DataVector.hpp"
18 #include "NumericalAlgorithms/LinearOperators/Divergence.tpp"
22 #include "Utilities/Numeric.hpp"
23 #include "Utilities/Tuple.hpp"
24 
25 namespace FirstOrderEllipticSolutionsTestHelpers {
26 namespace detail {
27 namespace Tags {
28 
29 template <typename Tag>
30 struct OperatorAppliedTo : db::SimpleTag, db::PrefixTag {
31  using type = typename Tag::type;
32  using tag = Tag;
33 };
34 
35 } // namespace Tags
36 
37 template <typename System, typename SolutionType, typename... Maps,
38  typename... FluxesArgs, typename... SourcesArgs,
39  typename... PrimalFields, typename... AuxiliaryFields,
40  typename... PrimalFluxes, typename... AuxiliaryFluxes>
41 void verify_solution_impl(
42  const SolutionType& solution, const Mesh<System::volume_dim>& mesh,
44  coord_map,
45  const double tolerance, const std::tuple<FluxesArgs...>& fluxes_args,
46  const std::tuple<SourcesArgs...>& sources_args,
47  tmpl::list<PrimalFields...> /*meta*/,
48  tmpl::list<AuxiliaryFields...> /*meta*/,
49  tmpl::list<PrimalFluxes...> /*meta*/,
50  tmpl::list<AuxiliaryFluxes...> /*meta*/) {
51  using all_fields = tmpl::list<PrimalFields..., AuxiliaryFields...>;
52  using all_fluxes = tmpl::list<PrimalFluxes..., AuxiliaryFluxes...>;
53  using FluxesComputer = typename System::fluxes_computer;
54  using SourcesComputer = typename System::sources_computer;
55  CAPTURE(mesh);
56 
57  const size_t num_points = mesh.number_of_grid_points();
58  const auto logical_coords = logical_coordinates(mesh);
59  const auto inertial_coords = coord_map(logical_coords);
60  const auto solution_fields = variables_from_tagged_tuple(
61  solution.variables(inertial_coords, all_fields{}));
62 
63  // Apply operator to solution fields: -div(F) + S = f
64  Variables<all_fluxes> fluxes{num_points};
65  std::apply(
66  [&fluxes, &solution_fields](const auto&... expanded_fluxes_args) {
67  FluxesComputer::apply(make_not_null(&get<PrimalFluxes>(fluxes))...,
68  expanded_fluxes_args...,
69  get<AuxiliaryFields>(solution_fields)...);
70  FluxesComputer::apply(make_not_null(&get<AuxiliaryFluxes>(fluxes))...,
71  expanded_fluxes_args...,
72  get<PrimalFields>(solution_fields)...);
73  },
74  fluxes_args);
75  Variables<db::wrap_tags_in<Tags::OperatorAppliedTo, all_fields>>
76  operator_applied_to_fields{num_points};
77  divergence(make_not_null(&operator_applied_to_fields), fluxes, mesh,
78  coord_map.inv_jacobian(logical_coords));
79  operator_applied_to_fields *= -1.;
80  std::apply(
81  [&operator_applied_to_fields, &solution_fields,
82  &fluxes](const auto&... expanded_sources_args) {
83  SourcesComputer::apply(
84  make_not_null(&get<Tags::OperatorAppliedTo<AuxiliaryFields>>(
85  operator_applied_to_fields))...,
86  expanded_sources_args..., get<PrimalFields>(solution_fields)...);
87  SourcesComputer::apply(
88  make_not_null(&get<Tags::OperatorAppliedTo<PrimalFields>>(
89  operator_applied_to_fields))...,
90  expanded_sources_args..., get<PrimalFields>(solution_fields)...,
91  get<PrimalFluxes>(fluxes)...);
92  },
93  sources_args);
94 
95  // Set RHS for auxiliary equations to minus auxiliary fields, and for primal
96  // equations to the solution's fixed sources f(x)
97  Variables<db::wrap_tags_in<::Tags::FixedSource, all_fields>> fixed_sources{
98  num_points, 0.};
100  get<AuxiliaryFields>(solution_fields))...);
101  fixed_sources *= -1.;
102  fixed_sources.assign_subset(solution.variables(
103  inertial_coords, tmpl::list<::Tags::FixedSource<PrimalFields>...>{}));
104 
105  // Check error norms against the given tolerance
106  tmpl::for_each<all_fields>([&operator_applied_to_fields, &fixed_sources,
107  &tolerance](auto field_tag_v) {
108  using field_tag = tmpl::type_from<decltype(field_tag_v)>;
109  const auto& operator_applied_to_field =
110  get<Tags::OperatorAppliedTo<field_tag>>(operator_applied_to_fields);
111  const auto& fixed_source =
112  get<::Tags::FixedSource<field_tag>>(fixed_sources);
113  double l2_error_square = 0.;
114  double linf_error = 0.;
115  for (size_t i = 0; i < operator_applied_to_field.size(); ++i) {
116  const auto error = abs(operator_applied_to_field[i] - fixed_source[i]);
117  l2_error_square += alg::accumulate(square(error), 0.) / error.size();
118  const double component_linf_error = *alg::max_element(error);
119  if (component_linf_error > linf_error) {
120  linf_error = component_linf_error;
121  }
122  }
123  const double l2_error =
124  sqrt(l2_error_square / operator_applied_to_field.size());
125  CAPTURE(db::tag_name<field_tag>());
126  CAPTURE(l2_error);
127  CAPTURE(linf_error);
128  CHECK(l2_error < tolerance);
129  CHECK(linf_error < tolerance);
130  });
131 }
132 
133 } // namespace detail
134 
135 /// \ingroup TestingFrameworkGroup
136 /// Test that the `solution` numerically solves the `System` on the given grid
137 /// for the given tolerance
138 /// @{
139 template <typename System, typename SolutionType, typename... Maps,
140  typename... FluxesArgs, typename... SourcesArgs>
142  const SolutionType& solution, const Mesh<System::volume_dim>& mesh,
144  coord_map,
145  const double tolerance, const std::tuple<FluxesArgs...>& fluxes_args,
146  const std::tuple<SourcesArgs...>& sources_args = std::tuple<>{}) {
147  detail::verify_solution_impl<System>(
148  solution, mesh, coord_map, tolerance, fluxes_args, sources_args,
149  typename System::primal_fields{}, typename System::auxiliary_fields{},
150  typename System::primal_fluxes{}, typename System::auxiliary_fluxes{});
151 }
152 
153 template <typename System, typename SolutionType, typename... Maps>
155  const SolutionType& solution, const Mesh<System::volume_dim>& mesh,
157  coord_map,
158  const double tolerance) {
159  using argument_tags = tmpl::remove_duplicates<
160  tmpl::append<typename System::fluxes_computer::argument_tags,
161  typename System::sources_computer::argument_tags>>;
162  const auto background_fields = [&solution, &mesh, &coord_map]() {
163  if constexpr (tmpl::size<argument_tags>::value > 0) {
164  const auto logical_coords = logical_coordinates(mesh);
165  const auto inertial_coords = coord_map(logical_coords);
166  const auto inv_jacobian = coord_map.inv_jacobian(logical_coords);
167  return solution.variables(inertial_coords, mesh, inv_jacobian,
168  argument_tags{});
169  } else {
170  (void)solution;
171  (void)mesh;
172  (void)coord_map;
173  return tuples::TaggedTuple<>{};
174  }
175  }();
176  const auto get_items = [](const auto&... args) {
177  return std::forward_as_tuple(args...);
178  };
179  verify_solution<System>(
180  solution, mesh, coord_map, tolerance,
181  tuples::apply<typename System::fluxes_computer::argument_tags>(
182  get_items, background_fields),
183  tuples::apply<typename System::sources_computer::argument_tags>(
184  get_items, background_fields));
185 }
186 /// @}
187 
188 /*!
189  * \ingroup TestingFrameworkGroup
190  * Test that the `solution` numerically solves the `System` on the given grid
191  * and that the discretization error decreases as expected for a smooth
192  * function.
193  *
194  * \details We expect exponential convergence for a smooth solution, so the
195  * tolerance is computed as
196  *
197  * \f{equation}
198  * C_1 \exp{\left(-C_2 * N_\mathrm{points}\right)}
199  * \f}
200  *
201  * where \f$C_1\f$ is the `tolerance_offset`, \f$C_2\f$ is the
202  * `tolerance_scaling` and \f$N_\mathrm{points}\f$ is the number of grid points
203  * per dimension.
204  */
205 template <typename System, typename SolutionType,
206  size_t Dim = System::volume_dim, typename... Maps,
207  typename PackageFluxesArgs>
209  const SolutionType& solution,
211  coord_map,
212  const double tolerance_offset, const double tolerance_scaling,
213  PackageFluxesArgs&& package_fluxes_args) {
214  INFO("Verify smooth solution");
215  for (size_t num_points = Spectral::minimum_number_of_points<
216  Spectral::Basis::Legendre, Spectral::Quadrature::GaussLobatto>;
217  num_points <=
218  Spectral::maximum_number_of_points<Spectral::Basis::Legendre>;
219  num_points++) {
220  CAPTURE(num_points);
221  const double tolerance =
222  tolerance_offset * exp(-tolerance_scaling * num_points);
223  CAPTURE(tolerance);
224  const Mesh<Dim> mesh{num_points, Spectral::Basis::Legendre,
225  Spectral::Quadrature::GaussLobatto};
226  FirstOrderEllipticSolutionsTestHelpers::verify_solution<System>(
227  solution, mesh, coord_map, tolerance, package_fluxes_args(mesh));
228  }
229 }
230 
231 /*!
232  * \ingroup TestingFrameworkGroup
233  * Test that the `solution` numerically solves the `System` on the given grid
234  * and that the discretization error decreases as a power law.
235  *
236  * \details The tolerance is computed as
237  *
238  * \f{equation}
239  * C \left(N_\mathrm{points}\right)^{-p}
240  * \f}
241  *
242  * where \f$C\f$ is the `tolerance_offset`, \f$p\f$ is the `tolerance_pow` and
243  * \f$N_\mathrm{points}\f$ is the number of grid points per dimension.
244  */
245 template <typename System, typename SolutionType,
246  size_t Dim = System::volume_dim, typename... Maps>
248  const SolutionType& solution,
250  coord_map,
251  const double tolerance_offset, const double tolerance_pow) {
252  INFO("Verify solution with power-law convergence");
253  for (size_t num_points = Spectral::minimum_number_of_points<
254  Spectral::Basis::Legendre, Spectral::Quadrature::GaussLobatto>;
255  num_points <=
256  Spectral::maximum_number_of_points<Spectral::Basis::Legendre>;
257  num_points++) {
258  CAPTURE(num_points);
259  const double tolerance = tolerance_offset * pow(num_points, -tolerance_pow);
260  CAPTURE(tolerance);
261  FirstOrderEllipticSolutionsTestHelpers::verify_solution<System>(
262  solution,
263  Mesh<Dim>{num_points, Spectral::Basis::Legendre,
264  Spectral::Quadrature::GaussLobatto},
265  coord_map, tolerance);
266  }
267 }
268 
269 } // namespace FirstOrderEllipticSolutionsTestHelpers
expand_pack
constexpr void expand_pack(Ts &&...) noexcept
Allows zero-cost unordered expansion of a parameter.
Definition: TMPL.hpp:585
FirstOrderEllipticSolutionsTestHelpers::verify_smooth_solution
void verify_smooth_solution(const SolutionType &solution, const domain::CoordinateMap< Frame::Logical, Frame::Inertial, Maps... > &coord_map, const double tolerance_offset, const double tolerance_scaling, PackageFluxesArgs &&package_fluxes_args)
Definition: FirstOrderEllipticSolutionsTestHelpers.hpp:208
variables_from_tagged_tuple
Variables< tmpl::list< Tags... > > variables_from_tagged_tuple(const tuples::TaggedTuple< Tags... > &tuple) noexcept
Definition: Variables.hpp:842
pow
constexpr decltype(auto) pow(const T &t) noexcept
Compute t^N where N is an integer (positive or negative)
Definition: ConstantExpressions.hpp:160
sqrt
auto sqrt(const TensorExpression< T, typename T::type, tmpl::list<>, tmpl::list<>, tmpl::list<>> &t)
Returns the tensor expression representing the square root of a tensor expression that evaluates to a...
Definition: SquareRoot.hpp:83
db::PrefixTag
Mark a struct as a prefix tag by inheriting from this.
Definition: Tag.hpp:103
square
constexpr decltype(auto) square(const T &x)
Compute the square of x
Definition: ConstantExpressions.hpp:55
MakeWithRandomValues.hpp
std::tuple
db::get
const auto & get(const DataBox< TagList > &box) noexcept
Retrieve the item with tag Tag from the DataBox.
Definition: DataBox.hpp:791
db::SimpleTag
Mark a struct as a simple tag by inheriting from this.
Definition: Tag.hpp:36
evolution::dg::subcell::fd::mesh
Mesh< Dim > mesh(const Mesh< Dim > &dg_mesh) noexcept
Computes the cell-centered finite-difference mesh from the DG mesh, using grid points per dimension,...
divergence
auto divergence(const Variables< FluxTags > &F, const Mesh< Dim > &mesh, const InverseJacobian< DataVector, Dim, Frame::Logical, DerivativeFrame > &inverse_jacobian) noexcept -> Variables< db::wrap_tags_in< Tags::div, FluxTags >>
Compute the (Euclidean) divergence of fluxes.
CoordinateMap.hpp
TestHelpers.hpp
domain::CoordinateMap::inv_jacobian
constexpr InverseJacobian< double, dim, SourceFrame, TargetFrame > inv_jacobian(tnsr::I< double, dim, SourceFrame > source_point, const double time=std::numeric_limits< double >::signaling_NaN(), const std::unordered_map< std::string, std::unique_ptr< domain::FunctionsOfTime::FunctionOfTime >> &functions_of_time=std::unordered_map< std::string, std::unique_ptr< domain::FunctionsOfTime::FunctionOfTime >>{}) const noexcept override
Compute the inverse Jacobian of the Maps... at the point(s) source_point
Definition: CoordinateMap.hpp:337
Tags::FixedSource
Prefix indicating a source term that is independent of dynamic variables.
Definition: Prefixes.hpp:75
FirstOrderEllipticSolutionsTestHelpers::verify_solution_with_power_law_convergence
void verify_solution_with_power_law_convergence(const SolutionType &solution, const domain::CoordinateMap< Frame::Logical, Frame::Inertial, Maps... > &coord_map, const double tolerance_offset, const double tolerance_pow)
Definition: FirstOrderEllipticSolutionsTestHelpers.hpp:247
Spectral::minimum_number_of_points
constexpr size_t minimum_number_of_points
Minimum number of possible collocation points for a quadrature type.
Definition: Spectral.hpp:128
cstddef
tuples::TaggedTuple<>
Definition: TaggedTuple.hpp:494
MakeWithValue.hpp
LogicalCoordinates.hpp
ConstantExpressions.hpp
Tuple.hpp
Mesh
Holds the number of grid points, basis, and quadrature in each direction of the computational grid.
Definition: Mesh.hpp:49
logical_coordinates
void logical_coordinates(gsl::not_null< tnsr::I< DataVector, VolumeDim, Frame::Logical > * > logical_coords, const Mesh< VolumeDim > &mesh) noexcept
Compute the logical coordinates in an Element.
limits
Tensor.hpp
make_not_null
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,...
Definition: Gsl.hpp:880
alg::max_element
constexpr decltype(auto) max_element(const Container &c)
Convenience wrapper around std::max_element.
Definition: Algorithm.hpp:327
FirstOrderEllipticSolutionsTestHelpers::verify_solution
void verify_solution(const SolutionType &solution, const Mesh< System::volume_dim > &mesh, const domain::CoordinateMap< Frame::Logical, Frame::Inertial, Maps... > coord_map, const double tolerance)
Definition: FirstOrderEllipticSolutionsTestHelpers.hpp:154
alg::accumulate
decltype(auto) accumulate(const Container &c, T init)
Convenience wrapper around std::accumulate, returns std::accumulate(begin(c), end(c),...
Definition: Numeric.hpp:54
Mesh.hpp
domain::CoordinateMap
A coordinate map or composition of coordinate maps.
Definition: CoordinateMap.hpp:240