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"
25 namespace FirstOrderEllipticSolutionsTestHelpers {
29 template <
typename Tag>
31 using type =
typename Tag::type;
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(
47 tmpl::list<PrimalFields...> ,
48 tmpl::list<AuxiliaryFields...> ,
49 tmpl::list<PrimalFluxes...> ,
50 tmpl::list<AuxiliaryFluxes...> ) {
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;
59 const auto inertial_coords = coord_map(logical_coords);
61 solution.variables(inertial_coords, all_fields{}));
64 Variables<all_fluxes> fluxes{num_points};
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)...);
75 Variables<db::wrap_tags_in<Tags::OperatorAppliedTo, all_fields>>
76 operator_applied_to_fields{num_points};
79 operator_applied_to_fields *= -1.;
81 [&operator_applied_to_fields, &solution_fields,
82 &fluxes](
const auto&... expanded_sources_args) {
83 SourcesComputer::apply(
85 operator_applied_to_fields))...,
86 expanded_sources_args..., get<PrimalFields>(solution_fields)...);
87 SourcesComputer::apply(
89 operator_applied_to_fields))...,
90 expanded_sources_args..., get<PrimalFields>(solution_fields)...,
91 get<PrimalFluxes>(fluxes)...);
97 Variables<db::wrap_tags_in<::Tags::FixedSource, all_fields>> fixed_sources{
100 get<AuxiliaryFields>(solution_fields))...);
101 fixed_sources *= -1.;
102 fixed_sources.assign_subset(solution.variables(
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]);
119 if (component_linf_error > linf_error) {
120 linf_error = component_linf_error;
123 const double l2_error =
124 sqrt(l2_error_square / operator_applied_to_field.size());
125 CAPTURE(db::tag_name<field_tag>());
128 CHECK(l2_error < tolerance);
129 CHECK(linf_error < tolerance);
138 template <
typename System,
typename SolutionType,
typename... Maps,
139 typename... FluxesArgs,
typename... SourcesArgs>
144 const double tolerance,
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{});
170 template <
typename System,
typename SolutionType,
171 size_t Dim = System::volume_dim,
typename... Maps,
172 typename PackageFluxesArgs>
174 const SolutionType& solution,
177 const double tolerance_offset,
const double tolerance_scaling,
178 PackageFluxesArgs&& package_fluxes_args) {
179 INFO(
"Verify smooth solution");
181 Spectral::Basis::Legendre, Spectral::Quadrature::GaussLobatto>;
183 Spectral::maximum_number_of_points<Spectral::Basis::Legendre>;
186 const double tolerance =
187 tolerance_offset * exp(-tolerance_scaling * num_points);
189 const Mesh<Dim> mesh{num_points, Spectral::Basis::Legendre,
190 Spectral::Quadrature::GaussLobatto};
191 FirstOrderEllipticSolutionsTestHelpers::verify_solution<System>(
192 solution, mesh, coord_map, tolerance, package_fluxes_args(mesh));
210 template <
typename System,
typename SolutionType,
211 size_t Dim = System::volume_dim,
typename... Maps>
213 const SolutionType& solution,
216 const double tolerance_offset,
const double tolerance_pow) {
217 INFO(
"Verify solution with power-law convergence");
219 Spectral::Basis::Legendre, Spectral::Quadrature::GaussLobatto>;
221 Spectral::maximum_number_of_points<Spectral::Basis::Legendre>;
224 const double tolerance = tolerance_offset *
pow(num_points, -tolerance_pow);
226 FirstOrderEllipticSolutionsTestHelpers::verify_solution<System>(
228 Mesh<Dim>{num_points, Spectral::Basis::Legendre,
229 Spectral::Quadrature::GaussLobatto},
230 coord_map, tolerance);