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/DataBox/PrefixHelpers.hpp"
11 #include "DataStructures/DataBox/Tag.hpp"
12 #include "DataStructures/DataBox/TagName.hpp"
13 #include "DataStructures/DataVector.hpp"
17 #include "Elliptic/FirstOrderOperator.hpp"
20 #include "NumericalAlgorithms/LinearOperators/Divergence.tpp"
24 #include "Utilities/Numeric.hpp"
25 #include "Utilities/Tuple.hpp"
26 
27 namespace FirstOrderEllipticSolutionsTestHelpers {
28 
29 namespace Tags {
30 
31 template <typename Tag>
33  using type = typename Tag::type;
34  using tag = Tag;
35 };
36 
37 } // namespace Tags
38 
39 /// \ingroup TestingFrameworkGroup
40 /// Test that the `solution` numerically solves the `System` on the given grid
41 /// for the given tolerance
42 template <typename System, typename SolutionType, typename... Maps,
43  typename... FluxesArgs, typename... SourcesArgs>
45  const SolutionType& solution,
46  const typename System::fluxes& fluxes_computer,
47  const Mesh<System::volume_dim>& mesh,
49  coord_map,
50  const double tolerance,
51  const std::tuple<FluxesArgs...>& fluxes_args = std::tuple<>{},
52  const std::tuple<SourcesArgs...>& sources_args = std::tuple<>{}) {
53  using all_fields = typename System::fields_tag::tags_list;
54  using primal_fields = typename System::primal_fields;
55  using auxiliary_fields = typename System::auxiliary_fields;
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
64  auto fluxes = std::apply(
65  [&solution_fields,
66  &fluxes_computer](const auto&... expanded_fluxes_args) {
67  return ::elliptic::first_order_fluxes<System::volume_dim, primal_fields,
68  auxiliary_fields>(
69  solution_fields, fluxes_computer, expanded_fluxes_args...);
70  },
71  fluxes_args);
72  auto div_fluxes = divergence(std::move(fluxes), mesh,
73  coord_map.inv_jacobian(logical_coords));
74  auto sources = std::apply(
75  [&solution_fields](const auto&... expanded_sources_args) {
76  return ::elliptic::first_order_sources<primal_fields, auxiliary_fields,
77  typename System::sources>(
78  solution_fields, expanded_sources_args...);
79  },
80  sources_args);
81  Variables<db::wrap_tags_in<Tags::OperatorAppliedTo, all_fields>>
82  operator_applied_to_fields{num_points};
83  elliptic::first_order_operator(make_not_null(&operator_applied_to_fields),
84  std::move(div_fluxes), std::move(sources));
85 
86  // Set fixed sources to zero for auxiliary fields, and retrieve the fixed
87  // sources for primal fields from the solution
88  Variables<db::wrap_tags_in<::Tags::FixedSource, all_fields>> fixed_sources{
89  num_points, 0.};
90  fixed_sources.assign_subset(solution.variables(
92 
93  // Check error norms against the given tolerance
94  tmpl::for_each<all_fields>([&operator_applied_to_fields, &fixed_sources,
95  &tolerance](auto field_tag_v) {
96  using field_tag = tmpl::type_from<decltype(field_tag_v)>;
97  const auto& operator_applied_to_field =
98  get<Tags::OperatorAppliedTo<field_tag>>(operator_applied_to_fields);
99  const auto& fixed_source =
100  get<::Tags::FixedSource<field_tag>>(fixed_sources);
101  double l2_error_square = 0.;
102  double linf_error = 0.;
103  for (size_t i = 0; i < operator_applied_to_field.size(); ++i) {
104  const auto error = abs(operator_applied_to_field[i] - fixed_source[i]);
105  l2_error_square += alg::accumulate(square(error), 0.) / error.size();
106  const double component_linf_error = *alg::max_element(error);
107  if (component_linf_error > linf_error) {
108  linf_error = component_linf_error;
109  }
110  }
111  const double l2_error =
112  sqrt(l2_error_square / operator_applied_to_field.size());
113  CAPTURE(db::tag_name<field_tag>());
114  CAPTURE(l2_error);
115  CAPTURE(linf_error);
116  CHECK(l2_error < tolerance);
117  CHECK(linf_error < tolerance);
118  });
119 }
120 
121 /*!
122  * \ingroup TestingFrameworkGroup
123  * Test that the `solution` numerically solves the `System` on the given grid
124  * and that the discretization error decreases as expected for a smooth
125  * function.
126  *
127  * \details We expect exponential convergence for a smooth solution, so the
128  * tolerance is computed as
129  *
130  * \f{equation}
131  * C_1 \exp{\left(-C_2 * N_\mathrm{points}\right)}
132  * \f}
133  *
134  * where \f$C_1\f$ is the `tolerance_offset`, \f$C_2\f$ is the
135  * `tolerance_scaling` and \f$N_\mathrm{points}\f$ is the number of grid points
136  * per dimension.
137  */
138 template <typename System, typename SolutionType,
139  size_t Dim = System::volume_dim, typename... Maps,
140  typename PackageFluxesArgs>
142  const SolutionType& solution,
143  const typename System::fluxes& fluxes_computer,
145  coord_map,
146  const double tolerance_offset, const double tolerance_scaling,
147  PackageFluxesArgs&& package_fluxes_args) {
148  INFO("Verify smooth solution");
149  for (size_t num_points = Spectral::minimum_number_of_points<
150  Spectral::Basis::Legendre, Spectral::Quadrature::GaussLobatto>;
151  num_points <=
152  Spectral::maximum_number_of_points<Spectral::Basis::Legendre>;
153  num_points++) {
154  CAPTURE(num_points);
155  const double tolerance =
156  tolerance_offset * exp(-tolerance_scaling * num_points);
157  CAPTURE(tolerance);
158  const Mesh<Dim> mesh{num_points, Spectral::Basis::Legendre,
159  Spectral::Quadrature::GaussLobatto};
160  FirstOrderEllipticSolutionsTestHelpers::verify_solution<System>(
161  solution, fluxes_computer, mesh, coord_map, tolerance,
162  package_fluxes_args(mesh));
163  }
164 }
165 
166 /*!
167  * \ingroup TestingFrameworkGroup
168  * Test that the `solution` numerically solves the `System` on the given grid
169  * and that the discretization error decreases as a power law.
170  *
171  * \details The tolerance is computed as
172  *
173  * \f{equation}
174  * C \left(N_\mathrm{points}\right)^{-p}
175  * \f}
176  *
177  * where \f$C\f$ is the `tolerance_offset`, \f$p\f$ is the `tolerance_pow` and
178  * \f$N_\mathrm{points}\f$ is the number of grid points per dimension.
179  */
180 template <typename System, typename SolutionType,
181  size_t Dim = System::volume_dim, typename... Maps>
183  const SolutionType& solution,
184  const typename System::fluxes& fluxes_computer,
186  coord_map,
187  const double tolerance_offset, const double tolerance_pow) {
188  INFO("Verify solution with power-law convergence");
189  for (size_t num_points = Spectral::minimum_number_of_points<
190  Spectral::Basis::Legendre, Spectral::Quadrature::GaussLobatto>;
191  num_points <=
192  Spectral::maximum_number_of_points<Spectral::Basis::Legendre>;
193  num_points++) {
194  CAPTURE(num_points);
195  const double tolerance = tolerance_offset * pow(num_points, -tolerance_pow);
196  CAPTURE(tolerance);
197  FirstOrderEllipticSolutionsTestHelpers::verify_solution<System>(
198  solution, fluxes_computer,
199  Mesh<Dim>{num_points, Spectral::Basis::Legendre,
200  Spectral::Quadrature::GaussLobatto},
201  coord_map, tolerance);
202  }
203 }
204 
205 } // namespace FirstOrderEllipticSolutionsTestHelpers
elliptic::first_order_operator
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
variables_from_tagged_tuple
Variables< tmpl::list< Tags... > > variables_from_tagged_tuple(const tuples::TaggedTuple< Tags... > &tuple) noexcept
Definition: Variables.hpp:798
DataBoxTag.hpp
pow
constexpr decltype(auto) pow(const T &t) noexcept
Compute t^N where N is an integer (positive or negative)
Definition: ConstantExpressions.hpp:160
db::PrefixTag
Marks an item as being a prefix to another tag.
Definition: Tag.hpp:66
square
constexpr decltype(auto) square(const T &x)
Compute the square of x
Definition: ConstantExpressions.hpp:55
MakeWithRandomValues.hpp
std::tuple
db::SimpleTag
Tags for the DataBox inherit from this type.
Definition: Tag.hpp:23
CoordinateMap.hpp
std::sqrt
T sqrt(T... args)
divergence
Scalar< DataVector > divergence(const tnsr::I< DataVector, Dim, DerivativeFrame > &input, const Mesh< Dim > &mesh, const InverseJacobian< DataVector, Dim, Frame::Logical, DerivativeFrame > &inverse_jacobian) noexcept
Compute the divergence of the vector input
Definition: Divergence.cpp:26
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:333
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:123
db::apply
constexpr auto apply(F &&f, const DataBox< BoxTags > &box, Args &&... args) noexcept
Apply the invokable f with argument Tags TagsList from DataBox box
Definition: DataBox.hpp:1424
cstddef
MakeWithValue.hpp
LogicalCoordinates.hpp
Mesh::number_of_grid_points
size_t number_of_grid_points() const noexcept
The total number of grid points in all dimensions.
Definition: Mesh.hpp:127
ConstantExpressions.hpp
Tuple.hpp
FirstOrderEllipticSolutionsTestHelpers::verify_solution_with_power_law_convergence
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:182
alg::max_element
decltype(auto) max_element(const Container &c)
Convenience wrapper around std::max_element.
Definition: Algorithm.hpp:327
Mesh
Holds the number of grid points, basis, and quadrature in each direction of the computational grid.
Definition: Mesh.hpp:47
FirstOrderEllipticSolutionsTestHelpers::verify_solution
void verify_solution(const SolutionType &solution, const typename System::fluxes &fluxes_computer, const Mesh< System::volume_dim > &mesh, const domain::CoordinateMap< Frame::Logical, Frame::Inertial, Maps... > coord_map, const double tolerance, const std::tuple< FluxesArgs... > &fluxes_args=std::tuple<>{}, const std::tuple< SourcesArgs... > &sources_args=std::tuple<>{})
Definition: FirstOrderEllipticSolutionsTestHelpers.hpp:44
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.
FirstOrderEllipticSolutionsTestHelpers::Tags::OperatorAppliedTo
Definition: FirstOrderEllipticSolutionsTestHelpers.hpp:32
limits
FirstOrderEllipticSolutionsTestHelpers::verify_smooth_solution
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, PackageFluxesArgs &&package_fluxes_args)
Definition: FirstOrderEllipticSolutionsTestHelpers.hpp:141
Tensor.hpp
db::wrap_tags_in
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: PrefixHelpers.hpp:30
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
Mesh.hpp
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
domain::CoordinateMap
A coordinate map or composition of coordinate maps.
Definition: CoordinateMap.hpp:236