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 
10 #include "DataStructures/DataVector.hpp"
14 #include "Domain/Mesh.hpp"
15 #include "Elliptic/FirstOrderOperator.hpp"
16 #include "NumericalAlgorithms/LinearOperators/Divergence.tpp"
19 #include "Utilities/Numeric.hpp"
22 
24 
25 namespace Tags {
26 
27 template <typename Tag>
29  using type = typename Tag::type;
30  using tag = Tag;
31 };
32 
33 } // namespace Tags
34 
35 /// \ingroup TestingFrameworkGroup
36 /// Test that the `solution` numerically solves the `System` on the given grid
37 /// for the given tolerance
38 template <typename System, typename SolutionType, typename... FluxesArgs,
39  typename... SourcesArgs, typename... Maps>
41  const SolutionType& solution,
42  const typename System::fluxes& fluxes_computer,
43  const FluxesArgs&... fluxes_args, const SourcesArgs&... sources_args,
44  const Mesh<System::volume_dim>& mesh,
46  coord_map,
47  const double tolerance) {
48  constexpr size_t volume_dim = System::volume_dim;
49  using all_fields = db::get_variables_tags_list<typename System::fields_tag>;
50  using primal_fields = typename System::primal_fields;
51  using auxiliary_fields = typename System::auxiliary_fields;
52 
53  const size_t num_points = mesh.number_of_grid_points();
54  const auto logical_coords = logical_coordinates(mesh);
55  const auto inertial_coords = coord_map(logical_coords);
56  const auto solution_fields = variables_from_tagged_tuple(
57  solution.variables(inertial_coords, all_fields{}));
58 
59  // Apply operator to solution fields
60  auto fluxes =
61  elliptic::first_order_fluxes<volume_dim, primal_fields, auxiliary_fields>(
62  solution_fields, fluxes_computer, fluxes_args...);
63  auto div_fluxes = divergence(std::move(fluxes), mesh,
64  coord_map.inv_jacobian(logical_coords));
65  auto sources = elliptic::first_order_sources<primal_fields, auxiliary_fields,
66  typename System::sources>(
67  solution_fields, sources_args...);
68  Variables<db::wrap_tags_in<Tags::OperatorAppliedTo, all_fields>>
69  operator_applied_to_fields{num_points};
70  elliptic::first_order_operator(make_not_null(&operator_applied_to_fields),
71  std::move(div_fluxes), std::move(sources));
72 
73  // Set fixed sources to zero for auxiliary fields, and retrieve the fixed
74  // sources for primal fields from the solution
75  Variables<db::wrap_tags_in<::Tags::FixedSource, all_fields>> fixed_sources{
76  num_points, 0.};
77  fixed_sources.assign_subset(solution.variables(
78  inertial_coords,
80 
81  // Check error norms against the given tolerance
82  tmpl::for_each<all_fields>([&operator_applied_to_fields, &fixed_sources,
83  &tolerance](auto field_tag_v) {
84  using field_tag = tmpl::type_from<decltype(field_tag_v)>;
85  const auto& operator_applied_to_field =
86  get<Tags::OperatorAppliedTo<field_tag>>(operator_applied_to_fields);
87  const auto& fixed_source =
88  get<::Tags::FixedSource<field_tag>>(fixed_sources);
89  double l2_error_square = 0.;
90  double linf_error = 0.;
91  for (size_t i = 0; i < operator_applied_to_field.size(); ++i) {
92  const auto error = abs(operator_applied_to_field[i] - fixed_source[i]);
93  l2_error_square += alg::accumulate(square(error), 0.) / error.size();
94  const double component_linf_error = *alg::max_element(error);
95  if (component_linf_error > linf_error) {
96  linf_error = component_linf_error;
97  }
98  }
99  const double l2_error =
100  sqrt(l2_error_square / operator_applied_to_field.size());
101  CAPTURE(db::tag_name<field_tag>());
102  CAPTURE(l2_error);
103  CAPTURE(linf_error);
104  CHECK(l2_error < tolerance);
105  CHECK(linf_error < tolerance);
106  });
107 }
108 
109 /*!
110  * \ingroup TestingFrameworkGroup
111  * Test that the `solution` numerically solves the `System` on the given grid
112  * and that the discretization error decreases as expected for a smooth
113  * function.
114  *
115  * \details We expect exponential convergence for a smooth solution, so the
116  * tolerance is computed as
117  *
118  * \f{equation}
119  * C_1 \exp{\left(-C_2 * N_\mathrm{points}\right)}
120  * \f}
121  *
122  * where \f$C_1\f$ is the `tolerance_offset`, \f$C_2\f$ is the
123  * `tolerance_scaling` and \f$N_\mathrm{points}\f$ is the number of grid points
124  * per dimension.
125  */
126 template <typename System, typename SolutionType,
127  size_t Dim = System::volume_dim, typename... Maps>
129  const SolutionType& solution,
130  const typename System::fluxes& fluxes_computer,
132  coord_map,
133  const double tolerance_offset, const double tolerance_scaling) {
134  INFO("Verify smooth solution");
135  for (size_t num_points = Spectral::minimum_number_of_points<
136  Spectral::Basis::Legendre, Spectral::Quadrature::GaussLobatto>;
137  num_points <=
138  Spectral::maximum_number_of_points<Spectral::Basis::Legendre>;
139  num_points++) {
140  CAPTURE(num_points);
141  const double tolerance =
142  tolerance_offset * exp(-tolerance_scaling * num_points);
143  CAPTURE(tolerance);
144  FirstOrderEllipticSolutionsTestHelpers::verify_solution<System>(
145  solution, fluxes_computer,
146  Mesh<Dim>{num_points, Spectral::Basis::Legendre,
147  Spectral::Quadrature::GaussLobatto},
148  coord_map, tolerance);
149  }
150 }
151 
152 /*!
153  * \ingroup TestingFrameworkGroup
154  * Test that the `solution` numerically solves the `System` on the given grid
155  * and that the discretization error decreases as a power law.
156  *
157  * \details The tolerance is computed as
158  *
159  * \f{equation}
160  * C \left(N_\mathrm{points}\right)^{-p}
161  * \f}
162  *
163  * where \f$C\f$ is the `tolerance_offset`, \f$p\f$ is the `tolerance_pow` and
164  * \f$N_\mathrm{points}\f$ is the number of grid points per dimension.
165  */
166 template <typename System, typename SolutionType,
167  size_t Dim = System::volume_dim, typename... Maps>
169  const SolutionType& solution,
170  const typename System::fluxes& fluxes_computer,
172  coord_map,
173  const double tolerance_offset, const double tolerance_pow) {
174  INFO("Verify solution with power-law convergence");
175  for (size_t num_points = Spectral::minimum_number_of_points<
176  Spectral::Basis::Legendre, Spectral::Quadrature::GaussLobatto>;
177  num_points <=
178  Spectral::maximum_number_of_points<Spectral::Basis::Legendre>;
179  num_points++) {
180  CAPTURE(num_points);
181  const double tolerance = tolerance_offset * pow(num_points, -tolerance_pow);
182  CAPTURE(tolerance);
183  FirstOrderEllipticSolutionsTestHelpers::verify_solution<System>(
184  solution, fluxes_computer,
185  Mesh<Dim>{num_points, Spectral::Basis::Legendre,
186  Spectral::Quadrature::GaussLobatto},
187  coord_map, tolerance);
188  }
189 }
190 
191 } // namespace FirstOrderEllipticSolutionsTestHelpers
decltype(auto) max_element(const Container &c)
Convenience wrapper around std::max_element.
Definition: Algorithm.hpp:277
void verify_solution(const SolutionType &solution, const typename System::fluxes &fluxes_computer, const FluxesArgs &... fluxes_args, const SourcesArgs &... sources_args, const Mesh< System::volume_dim > &mesh, const domain::CoordinateMap< Frame::Logical, Frame::Inertial, Maps... > coord_map, const double tolerance)
Test that the solution numerically solves the System on the given grid for the given tolerance...
Definition: FirstOrderEllipticSolutionsTestHelpers.hpp:40
void verify_solution_with_power_law_convergence(const SolutionType &solution, const typename System::fluxes &fluxes_computer, const domain::CoordinateMap< Frame::Logical, Frame::Inertial, Maps... > &coord_map, const double tolerance_offset, const double tolerance_pow)
Definition: FirstOrderEllipticSolutionsTestHelpers.hpp:168
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:266
void first_order_operator(const gsl::not_null< Variables< OperatorTags > *> operator_applied_to_vars, const Variables< DivFluxesTags > &div_fluxes, const Variables< SourcesTags > &sources) noexcept
Compute the bulk contribution to the operator represented by the OperatorTags.
Definition: FirstOrderOperator.hpp:161
Holds the number of grid points, basis, and quadrature in each direction of the computational grid...
Definition: Mesh.hpp:49
Helper functions for data structures used in unit tests.
Tags for the DataBox inherit from this type.
Definition: DataBoxTag.hpp:64
tmpl::transform< TagList, tmpl::bind< Wrapper, tmpl::_1, tmpl::pin< Args >... > > wrap_tags_in
Create a new list of Tags by wrapping each tag in TagList using the Wrapper.
Definition: DataBoxTag.hpp:515
A coordinate map or composition of coordinate maps.
Definition: CoordinateMap.hpp:183
void verify_smooth_solution(const SolutionType &solution, const typename System::fluxes &fluxes_computer, const domain::CoordinateMap< Frame::Logical, Frame::Inertial, Maps... > &coord_map, const double tolerance_offset, const double tolerance_scaling)
Definition: FirstOrderEllipticSolutionsTestHelpers.hpp:128
Definition: FirstOrderEllipticSolutionsTestHelpers.hpp:23
void first_order_sources(gsl::not_null< Variables< db::wrap_tags_in<::Tags::Source, VarsTags >> *> sources, const Variables< VarsTags > &vars, const SourcesArgs &... sources_args) noexcept
Compute the sources for the first-order formulation of elliptic systems.
Definition: FirstOrderOperator.hpp:125
decltype(auto) accumulate(const Container &c, T init)
Convenience wrapper around std::accumulate, returns std::accumulate(begin(c), end(c), init).
Definition: Numeric.hpp:54
size_t number_of_grid_points() const noexcept
The total number of grid points in all dimensions.
Definition: Mesh.hpp:124
Defines functions logical_coordinates and interface_logical_coordinates.
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.
Define simple functions for constant expressions.
Definition: DataBoxTag.hpp:29
Defines the class template Mesh.
Defines class CoordinateMap.
Defines classes for Tensor.
tnsr::I< DataVector, VolumeDim, Frame::Logical > logical_coordinates(const Mesh< VolumeDim > &mesh) noexcept
Compute the logical coordinates in an Element.
Definition: LogicalCoordinates.cpp:20
constexpr size_t minimum_number_of_points
Minimum number of possible collocation points for a quadrature type.
Definition: Spectral.hpp:95
Definition: FirstOrderEllipticSolutionsTestHelpers.hpp:28
Commonly used routines, functions and definitions shared amongst unit tests.
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:879
Marks an item as being a prefix to another tag.
Definition: DataBoxTag.hpp:111
Defines classes SimpleTag, PrefixTag, ComputeTag and several functions for retrieving tag info...
Variables< tmpl::list< Tags... > > variables_from_tagged_tuple(const tuples::TaggedTuple< Tags... > &tuple) noexcept
Construct a variables from the Tensors in a TaggedTuple.
Definition: Variables.hpp:831
decltype(auto) constexpr square(const T &x)
Compute the square of x
Definition: ConstantExpressions.hpp:52
decltype(auto) constexpr pow(const T &t) noexcept
Compute t^N where N is an integer (positive or negative)
Definition: ConstantExpressions.hpp:157
Defines make_with_value.