InterpolateOnElementTestHelpers.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
7 
8 #include <cmath>
9 #include <cstddef>
10 #include <random>
11 
14 #include "Domain/BlockLogicalCoordinates.hpp"
15 #include "Domain/Creators/Shell.hpp"
16 #include "Domain/Domain.hpp"
17 #include "Domain/ElementLogicalCoordinates.hpp"
18 #include "Domain/ElementMap.hpp"
20 #include "Domain/Structure/InitialElementIds.hpp"
21 #include "Framework/ActionTesting.hpp"
23 #include "NumericalAlgorithms/Interpolation/Actions/InterpolationTargetVarsFromElement.hpp"
24 #include "NumericalAlgorithms/Interpolation/InterpolationTarget.hpp"
25 #include "NumericalAlgorithms/Interpolation/PointInfoTag.hpp"
26 #include "NumericalAlgorithms/Interpolation/Tags.hpp"
28 #include "Parallel/PhaseDependentActionList.hpp"
29 #include "Time/Slab.hpp"
30 #include "Time/Tags.hpp"
31 #include "Time/Time.hpp"
32 #include "Time/TimeStepId.hpp"
33 #include "Utilities/Gsl.hpp"
34 #include "Utilities/TMPL.hpp"
35 
36 /// Holds code that is shared between multiple tests. Currently used by
37 /// - Test_InterpolateToTarget
38 /// - Test_InterpolateWithoutInterpolatorComponent.
40 
41 namespace Tags {
42 // Simple Variables tag for test.
44  using type = Scalar<DataVector>;
45 };
46 // Tag holding inertial-frame target points for test.
48  using type = tnsr::I<DataVector, 3, Frame::Inertial>;
49 };
51  using type = Scalar<DataVector>;
52 };
53 // Compute tag for test.
55  static void function(const gsl::not_null<Scalar<DataVector>*> result,
56  const Scalar<DataVector>& x) noexcept {
57  get<>(*result) = get<>(x) * 2.0;
58  }
59  using return_type = Scalar<DataVector>;
60  using argument_tags = tmpl::list<TestSolution>;
61  using base = MultiplyByTwo;
62 };
63 } // namespace Tags
64 
65 template <typename TagName, typename VarsTagList>
66 void fill_variables(
67  const gsl::not_null<Variables<VarsTagList>*> vars,
68  const tnsr::I<DataVector, 3, Frame::Inertial>& coords) noexcept {
69  // Some analytic solution used to fill the volume data and
70  // to test the interpolated data.
71  auto& solution = get<TagName>(*vars);
72  get(solution) =
73  (2.0 * get<0>(coords) + 3.0 * get<1>(coords) + 5.0 * get<2>(coords));
74 
75  if constexpr (std::is_same_v<TagName, Tags::MultiplyByTwo>) {
76  get(solution) *= 2.0;
77  } else if constexpr (not std::is_same_v<TagName, Tags::TestSolution>) {
78  ERROR("Do not understand the given TagName");
79  }
80 }
81 
82 template <typename InterpolationTargetTag>
84  template <
85  typename ParallelComponent, typename DbTags, typename Metavariables,
86  typename ArrayIndex, typename TemporalId,
88  static void apply(
91  const ArrayIndex& /*array_index*/,
92  const std::vector<Variables<
93  typename InterpolationTargetTag::vars_to_interpolate_to_target>>&
94  vars_src,
95  const std::vector<std::vector<size_t>>& global_offsets,
96  const TemporalId& /*temporal_id*/) noexcept {
97  CHECK(global_offsets.size() == vars_src.size());
98  // global_offsets and vars_src always have a size of 1 for calls
99  // directly from the elements; the outer vector is used only by
100  // the Interpolator parallel component.
101  CHECK(global_offsets.size() == 1);
102 
103  // Here we have received only some of the points.
104  const size_t num_pts_received = global_offsets[0].size();
105 
106  // Create a new target_points containing only the ones we have received.
107  const auto& all_target_points = db::get<Tags::TestTargetPoints>(box);
108  tnsr::I<DataVector, 3, Frame::Inertial> target_points(num_pts_received);
109  for (size_t i = 0; i < num_pts_received; ++i) {
110  for (size_t d = 0; d < 3; ++d) {
111  target_points.get(d)[i] =
112  all_target_points.get(d)[global_offsets[0][i]];
113  }
114  }
115 
116  static_assert(
117  tmpl::count<typename InterpolationTargetTag::
118  vars_to_interpolate_to_target>::value == 1,
119  "For this test, we assume only a single interpolated variable");
120  using solution_tag = tmpl::front<
121  typename InterpolationTargetTag::vars_to_interpolate_to_target>;
122  const auto& test_solution = get<solution_tag>(vars_src[0]);
123 
124  // Expected solution
125  Variables<tmpl::list<solution_tag>> expected_vars(vars_src[0].size());
126  fill_variables<solution_tag>(make_not_null(&expected_vars), target_points);
127  const auto& expected_solution = get<solution_tag>(expected_vars);
128 
129  // We don't have that many points, so interpolation is good for
130  // only a few digits.
131  Approx custom_approx = Approx::custom().epsilon(1.e-5).scale(1.0);
132  CHECK_ITERABLE_CUSTOM_APPROX(test_solution, expected_solution,
133  custom_approx);
134  }
135 };
136 
137 template <typename Metavariables, typename InterpolationTargetTag>
139  using metavariables = Metavariables;
141  using array_index = size_t;
142  using component_being_mocked =
144  using simple_tags = tmpl::list<Tags::TestTargetPoints>;
145  using phase_dependent_action_list = tmpl::list<Parallel::PhaseActions<
146  typename Metavariables::Phase, Metavariables::Phase::Initialization,
147  tmpl::list<ActionTesting::InitializeDataBox<simple_tags>>>>;
148  using replace_these_simple_actions =
150  InterpolationTargetTag>>;
151  using with_these_simple_actions = tmpl::list<
153 };
154 
155 template <typename DomainCreator>
157 make_volume_data_and_mesh(const DomainCreator& domain_creator,
158  const Domain<3>& domain,
159  const ElementId<3>& element_id) noexcept {
160  const auto& block = domain.blocks()[element_id.block_id()];
161  Mesh<3> mesh{domain_creator.initial_extents()[element_id.block_id()],
162  Spectral::Basis::Legendre, Spectral::Quadrature::GaussLobatto};
163  if (block.is_time_dependent()) {
164  ERROR("The block must be time-independent");
165  }
166 
167  // create volume data without the compute_items_on_source
168  ElementMap<3, Frame::Inertial> map{element_id,
169  block.stationary_map().get_clone()};
170  const auto inertial_coords = map(logical_coordinates(mesh));
171  Variables<tmpl::list<Tags::TestSolution>> vars(mesh.number_of_grid_points());
172  fill_variables<Tags::TestSolution>(make_not_null(&vars), inertial_coords);
173  return std::make_tuple(std::move(vars), mesh);
174 }
175 
176 // This is the main test; Takes a functor as an argument so that
177 // different tests can call it to do slightly different things.
178 template <typename Metavariables, typename elem_component, typename Functor,
179  typename... GlobalCacheTypes>
180 void test_interpolate_on_element(
181  Functor initialize_elements_and_queue_simple_actions,
182  GlobalCacheTypes... global_cache_items) noexcept {
183  using metavars = Metavariables;
184  using target_component =
185  mock_interpolation_target<metavars,
186  typename metavars::InterpolationTargetA>;
187 
188  const auto domain_creator =
189  domain::creators::Shell(0.9, 2.9, 2, {{7, 7}}, false);
190  const auto domain = domain_creator.create_domain();
191  Slab slab(0.0, 1.0);
192  TimeStepId temporal_id(true, 0, Time(slab, Rational(11, 15)));
193 
194  // Create Element_ids.
195  std::vector<ElementId<3>> element_ids{};
196  for (const auto& block : domain.blocks()) {
197  const auto initial_ref_levs =
198  domain_creator.initial_refinement_levels()[block.id()];
199  auto elem_ids = initial_element_ids(block.id(), initial_ref_levs);
200  element_ids.insert(element_ids.end(), elem_ids.begin(), elem_ids.end());
201  }
202 
203  // Create target points and interp_point_info
204  const size_t num_points = 10;
205  tnsr::I<DataVector, 3, Frame::Inertial> target_points(num_points);
206  const typename intrp::Tags::InterpPointInfo<
207  metavars>::type interp_point_info = [&domain, &target_points]() {
208  MAKE_GENERATOR(gen);
209  std::uniform_real_distribution<> r_dist(0.9001, 2.8999);
210  std::uniform_real_distribution<> theta_dist(0.0, M_PI);
211  std::uniform_real_distribution<> phi_dist(0.0, 2 * M_PI);
212  for (size_t i = 0; i < num_points; ++i) {
213  const double r = r_dist(gen);
214  const double theta = theta_dist(gen);
215  const double phi = phi_dist(gen);
216  get<0>(target_points)[i] = r * sin(theta) * cos(phi);
217  get<1>(target_points)[i] = r * sin(theta) * sin(phi);
218  get<2>(target_points)[i] = r * cos(theta);
219  }
220  typename intrp::Tags::InterpPointInfo<metavars>::type interp_point_info_l{};
221  get<intrp::Vars::PointInfoTag<typename metavars::InterpolationTargetA, 3>>(
222  interp_point_info_l) = block_logical_coordinates(domain, target_points);
223  return interp_point_info_l;
224  }();
225 
226  // Emplace target component.
228  tuples::tagged_tuple_from_typelist<
230  std::move(global_cache_items)...}};
232  metavars::Phase::Initialization);
233  ActionTesting::emplace_component_and_initialize<target_component>(
234  &runner, 0, {target_points});
235 
236  initialize_elements_and_queue_simple_actions(domain_creator, domain,
237  element_ids, interp_point_info,
238  runner, temporal_id);
239 
240  // Only some of the actions/events just invoked on elements (those elements
241  // which contain target points) will queue a simple action on the
242  // InterpolationTarget. Invoke those simple actions now.
243  while (not ActionTesting::is_simple_action_queue_empty<target_component>(
244  runner, 0)) {
245  runner.template invoke_queued_simple_action<target_component>(0);
246  }
247 }
248 } // namespace InterpolateOnElementTestHelpers
intrp::Tags::InterpPointInfo
The following tag is for the case in which interpolation bypasses the Interpolator ParallelComponent....
Definition: PointInfoTag.hpp:44
db::ComputeTag
Marks a DataBoxTag as being a compute item that executes a function.
Definition: Tag.hpp:109
ActionTesting::set_phase
void set_phase(const gsl::not_null< MockRuntimeSystem< Metavariables > * > runner, const typename Metavariables::Phase &phase) noexcept
Set the phase of all parallel components to phase
Definition: ActionTesting.hpp:1921
DataBoxTag.hpp
ActionTesting::MockRuntimeSystem
Definition: ActionTesting.hpp:1614
Slab
Definition: Slab.hpp:27
get
constexpr Tag::type & get(Variables< TagList > &v) noexcept
Return Tag::type pointing into the contiguous array.
Definition: Variables.hpp:639
Parallel::GlobalCache
Definition: ElementReceiveInterpPoints.hpp:16
std::cos
T cos(T... args)
std::vector
TestingFramework.hpp
Domain.hpp
random
InterpolateOnElementTestHelpers::Tags::TestSolution
Definition: InterpolateOnElementTestHelpers.hpp:43
std::tuple
db::SimpleTag
Tags for the DataBox inherit from this type.
Definition: Tag.hpp:23
CHECK_ITERABLE_CUSTOM_APPROX
#define CHECK_ITERABLE_CUSTOM_APPROX(a, b, appx)
Same as CHECK_ITERABLE_APPROX with user-defined Approx. The third argument should be of type Approx.
Definition: TestingFramework.hpp:151
cmath
InterpolateOnElementTestHelpers::Tags::MultiplyByTwoCompute
Definition: InterpolateOnElementTestHelpers.hpp:54
Spectral.hpp
InterpolateOnElementTestHelpers::Tags::TestTargetPoints
Definition: InterpolateOnElementTestHelpers.hpp:47
std::uniform_real_distribution
TestHelpers.hpp
ElementId
An ElementId uniquely labels an Element.
Definition: ElementId.hpp:49
ERROR
#define ERROR(m)
prints an error message to the standard error stream and aborts the program.
Definition: Error.hpp:36
InterpolateOnElementTestHelpers::Tags::MultiplyByTwo
Definition: InterpolateOnElementTestHelpers.hpp:50
Parallel::PhaseActions
List of all the actions to be executed in the specified phase.
Definition: PhaseDependentActionList.hpp:16
intrp::Actions::InterpolationTargetVarsFromElement
Receives interpolated variables from an Element on a subset of the target points.
Definition: InterpolationTargetVarsFromElement.hpp:61
domain::creators::Shell
Creates a 3D Domain in the shape of a hollow spherical shell consisting of six wedges.
Definition: Shell.hpp:41
DataBox.hpp
InterpolateOnElementTestHelpers::MockInterpolationTargetVarsFromElement
Definition: InterpolateOnElementTestHelpers.hpp:83
InterpolateOnElementTestHelpers
Holds code that is shared between multiple tests. Currently used by.
Definition: InterpolateOnElementTestHelpers.hpp:39
cstddef
LogicalCoordinates.hpp
ElementMap
The CoordinateMap for the Element from the Logical frame to the TargetFrame
Definition: ElementMap.hpp:33
InterpolateOnElementTestHelpers::mock_interpolation_target
Definition: InterpolateOnElementTestHelpers.hpp:138
Domain
A wrapper around a vector of Blocks that represent the computational domain.
Definition: Domain.hpp:38
GeneralizedHarmonic::phi
void phi(gsl::not_null< tnsr::iaa< DataType, SpatialDim, Frame > * > phi, const Scalar< DataType > &lapse, const tnsr::i< DataType, SpatialDim, Frame > &deriv_lapse, const tnsr::I< DataType, SpatialDim, Frame > &shift, const tnsr::iJ< DataType, SpatialDim, Frame > &deriv_shift, const tnsr::ii< DataType, SpatialDim, Frame > &spatial_metric, const tnsr::ijj< DataType, SpatialDim, Frame > &deriv_spatial_metric) noexcept
Computes the auxiliary variable used by the generalized harmonic formulation of Einstein's equations...
Time.hpp
block_logical_coordinates
auto block_logical_coordinates(const Domain< Dim > &domain, const tnsr::I< DataVector, Dim, Frame > &x, 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 >>{}) noexcept -> std::vector< boost::optional< IdPair< domain::BlockId, tnsr::I< double, Dim, ::Frame::Logical >>>>
TimeStepId
Definition: TimeStepId.hpp:25
std::sin
T sin(T... args)
DomainCreator
Base class for creating Domains from an option string.
Definition: DomainCreator.hpp:88
MAKE_GENERATOR
#define MAKE_GENERATOR(...)
MAKE_GENERATOR(NAME [, SEED]) declares a variable of name NAME containing a generator of type std::mt...
Definition: TestHelpers.hpp:417
ActionTesting::MockArrayChare
A mock class for the CMake-generated Parallel::Algorithms::Array
Definition: ActionTesting.hpp:1589
Mesh
Holds the number of grid points, basis, and quadrature in each direction of the computational grid.
Definition: Mesh.hpp:47
Parallel::get_const_global_cache_tags
tmpl::remove_duplicates< tmpl::flatten< tmpl::list< typename detail::get_const_global_cache_tags_from_parallel_struct< Metavariables >::type, tmpl::transform< typename Metavariables::component_list, detail::get_const_global_cache_tags_from_parallel_struct< tmpl::_1 > >, tmpl::transform< typename Metavariables::component_list, detail::get_const_global_cache_tags_from_pdal< tmpl::_1 > >> >> get_const_global_cache_tags
Given the metavariables, get a list of the unique tags that will specify the items in the GlobalCache...
Definition: ParallelComponentHelpers.hpp:83
TimeStepId.hpp
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.
Scalar
Tensor< T, Symmetry<>, index_list<> > Scalar
Definition: TypeAliases.hpp:21
Gsl.hpp
initial_element_ids
std::vector< ElementId< VolumeDim > > initial_element_ids(const std::vector< std::array< size_t, VolumeDim >> &initial_refinement_levels) noexcept
Create the ElementIds of the initial computational domain.
Definition: InitialElementIds.cpp:67
Time
Definition: Time.hpp:29
Slab.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
Requires
typename Requires_detail::requires_impl< B >::template_error_type_failed_to_meet_requirements_on_template_parameters Requires
Express requirements on the template parameters of a function or class, replaces std::enable_if_t
Definition: Requires.hpp:67
db::DataBox
Definition: InterpolationTargetWedgeSectionTorus.hpp:24
TMPL.hpp
Rational
Definition: Rational.hpp:27
gsl::not_null
Require a pointer to not be a nullptr
Definition: Gsl.hpp:183
intrp::InterpolationTarget
ParallelComponent representing a set of points to be interpolated to and a function to call upon inte...
Definition: InterpolationTarget.hpp:121