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 
34 /// Holds code that is shared between multiple tests. Currently used by
35 /// - Test_InterpolateToTarget
36 /// - Test_InterpolateWithoutInterpolatorComponent.
38 
39 namespace Tags {
40 // Simple Variables tag for test.
42  using type = Scalar<DataVector>;
43 };
44 // Tag holding inertial-frame target points for test.
46  using type = tnsr::I<DataVector, 3, Frame::Inertial>;
47 };
48 // ComputeItem for test.
50  using type = Scalar<DataVector>;
51 };
53  static Scalar<DataVector> function(const Scalar<DataVector>& x) noexcept {
54  auto result = make_with_value<Scalar<DataVector>>(x, 0.0);
55  get<>(result) = get<>(x) * 2.0;
56  return result;
57  }
58  using argument_tags = tmpl::list<TestSolution>;
59  using base = MultiplyByTwo;
60 };
61 } // namespace Tags
62 
63 template <typename TagName, typename VarsTagList>
64 void fill_variables(
65  const gsl::not_null<Variables<VarsTagList>*> vars,
66  const tnsr::I<DataVector, 3, Frame::Inertial>& coords) noexcept {
67  // Some analytic solution used to fill the volume data and
68  // to test the interpolated data.
69  auto& solution = get<TagName>(*vars);
70  get(solution) =
71  (2.0 * get<0>(coords) + 3.0 * get<1>(coords) + 5.0 * get<2>(coords));
72 
73  if constexpr (std::is_same_v<TagName, Tags::MultiplyByTwo>) {
74  get(solution) *= 2.0;
75  } else if constexpr (not std::is_same_v<TagName, Tags::TestSolution>) {
76  ERROR("Do not understand the given TagName");
77  }
78 }
79 
80 template <typename InterpolationTargetTag>
82  template <
83  typename ParallelComponent, typename DbTags, typename Metavariables,
84  typename ArrayIndex, typename TemporalId,
86  static void apply(
89  const ArrayIndex& /*array_index*/,
90  const std::vector<Variables<
91  typename InterpolationTargetTag::vars_to_interpolate_to_target>>&
92  vars_src,
93  const std::vector<std::vector<size_t>>& global_offsets,
94  const TemporalId& /*temporal_id*/) noexcept {
95  CHECK(global_offsets.size() == vars_src.size());
96  // global_offsets and vars_src always have a size of 1 for calls
97  // directly from the elements; the outer vector is used only by
98  // the Interpolator parallel component.
99  CHECK(global_offsets.size() == 1);
100 
101  // Here we have received only some of the points.
102  const size_t num_pts_received = global_offsets[0].size();
103 
104  // Create a new target_points containing only the ones we have received.
105  const auto& all_target_points = db::get<Tags::TestTargetPoints>(box);
106  tnsr::I<DataVector, 3, Frame::Inertial> target_points(num_pts_received);
107  for (size_t i = 0; i < num_pts_received; ++i) {
108  for (size_t d = 0; d < 3; ++d) {
109  target_points.get(d)[i] =
110  all_target_points.get(d)[global_offsets[0][i]];
111  }
112  }
113 
114  static_assert(
115  tmpl::count<typename InterpolationTargetTag::
116  vars_to_interpolate_to_target>::value == 1,
117  "For this test, we assume only a single interpolated variable");
118  using solution_tag = tmpl::front<
119  typename InterpolationTargetTag::vars_to_interpolate_to_target>;
120  const auto& test_solution = get<solution_tag>(vars_src[0]);
121 
122  // Expected solution
123  Variables<tmpl::list<solution_tag>> expected_vars(vars_src[0].size());
124  fill_variables<solution_tag>(make_not_null(&expected_vars), target_points);
125  const auto& expected_solution = get<solution_tag>(expected_vars);
126 
127  // We don't have that many points, so interpolation is good for
128  // only a few digits.
129  Approx custom_approx = Approx::custom().epsilon(1.e-5).scale(1.0);
130  CHECK_ITERABLE_CUSTOM_APPROX(test_solution, expected_solution,
131  custom_approx);
132  }
133 };
134 
135 template <typename Metavariables, typename InterpolationTargetTag>
137  using metavariables = Metavariables;
139  using array_index = size_t;
140  using component_being_mocked =
142  using simple_tags = tmpl::list<Tags::TestTargetPoints>;
143  using phase_dependent_action_list = tmpl::list<Parallel::PhaseActions<
144  typename Metavariables::Phase, Metavariables::Phase::Initialization,
145  tmpl::list<ActionTesting::InitializeDataBox<simple_tags>>>>;
146  using replace_these_simple_actions =
148  InterpolationTargetTag>>;
149  using with_these_simple_actions = tmpl::list<
151 };
152 
153 template <typename DomainCreator>
155 make_volume_data_and_mesh(const DomainCreator& domain_creator,
156  const Domain<3>& domain,
157  const ElementId<3>& element_id) noexcept {
158  const auto& block = domain.blocks()[element_id.block_id()];
159  Mesh<3> mesh{domain_creator.initial_extents()[element_id.block_id()],
160  Spectral::Basis::Legendre, Spectral::Quadrature::GaussLobatto};
161  if (block.is_time_dependent()) {
162  ERROR("The block must be time-independent");
163  }
164 
165  // create volume data without the compute_items_on_source
166  ElementMap<3, Frame::Inertial> map{element_id,
167  block.stationary_map().get_clone()};
168  const auto inertial_coords = map(logical_coordinates(mesh));
169  Variables<tmpl::list<Tags::TestSolution>> vars(mesh.number_of_grid_points());
170  fill_variables<Tags::TestSolution>(make_not_null(&vars), inertial_coords);
171  return std::make_tuple(std::move(vars), mesh);
172 }
173 
174 // This is the main test; Takes a functor as an argument so that
175 // different tests can call it to do slightly different things.
176 template <typename Metavariables, typename elem_component, typename Functor>
177 void test_interpolate_on_element(
178  Functor initialize_elements_and_queue_simple_actions) noexcept {
179  using metavars = Metavariables;
180  using target_component =
181  mock_interpolation_target<metavars,
182  typename metavars::InterpolationTargetA>;
183 
184  const auto domain_creator =
185  domain::creators::Shell(0.9, 2.9, 2, {{7, 7}}, false);
186  const auto domain = domain_creator.create_domain();
187  Slab slab(0.0, 1.0);
188  TimeStepId temporal_id(true, 0, Time(slab, Rational(11, 15)));
189 
190  // Create Element_ids.
191  std::vector<ElementId<3>> element_ids{};
192  for (const auto& block : domain.blocks()) {
193  const auto initial_ref_levs =
194  domain_creator.initial_refinement_levels()[block.id()];
195  auto elem_ids = initial_element_ids(block.id(), initial_ref_levs);
196  element_ids.insert(element_ids.end(), elem_ids.begin(), elem_ids.end());
197  }
198 
199  // Create target points and interp_point_info
200  const size_t num_points = 10;
201  tnsr::I<DataVector, 3, Frame::Inertial> target_points(num_points);
202  const typename intrp::Tags::InterpPointInfo<
203  metavars>::type interp_point_info = [&domain, &target_points]() {
204  MAKE_GENERATOR(gen);
205  std::uniform_real_distribution<> r_dist(0.9001, 2.8999);
206  std::uniform_real_distribution<> theta_dist(0.0, M_PI);
207  std::uniform_real_distribution<> phi_dist(0.0, 2 * M_PI);
208  for (size_t i = 0; i < num_points; ++i) {
209  const double r = r_dist(gen);
210  const double theta = theta_dist(gen);
211  const double phi = phi_dist(gen);
212  get<0>(target_points)[i] = r * sin(theta) * cos(phi);
213  get<1>(target_points)[i] = r * sin(theta) * sin(phi);
214  get<2>(target_points)[i] = r * cos(theta);
215  }
216  typename intrp::Tags::InterpPointInfo<metavars>::type interp_point_info_l{};
217  get<intrp::Vars::PointInfoTag<typename metavars::InterpolationTargetA, 3>>(
218  interp_point_info_l) = block_logical_coordinates(domain, target_points);
219  return interp_point_info_l;
220  }();
221 
222  // Emplace target component.
225  metavars::Phase::Initialization);
226  ActionTesting::emplace_component_and_initialize<target_component>(
227  &runner, 0, {target_points});
228 
229  initialize_elements_and_queue_simple_actions(domain_creator, domain,
230  element_ids, interp_point_info,
231  runner, temporal_id);
232 
233  // Only some of the actions/events just invoked on elements (those elements
234  // which contain target points) will queue a simple action on the
235  // InterpolationTarget. Invoke those simple actions now.
236  while (not ActionTesting::is_simple_action_queue_empty<target_component>(
237  runner, 0)) {
238  runner.template invoke_queued_simple_action<target_component>(0);
239  }
240 }
241 } // namespace InterpolateOnElementTestHelpers
Parallel::ConstGlobalCache
Definition: ElementReceiveInterpPoints.hpp:16
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:1920
DataBoxTag.hpp
ActionTesting::MockRuntimeSystem
Definition: ActionTesting.hpp:1613
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:689
std::cos
T cos(T... args)
std::vector
TestingFramework.hpp
Domain.hpp
random
InterpolateOnElementTestHelpers::Tags::TestSolution
Definition: InterpolateOnElementTestHelpers.hpp:41
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
Spectral.hpp
InterpolateOnElementTestHelpers::Tags::TestTargetPoints
Definition: InterpolateOnElementTestHelpers.hpp:45
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:49
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
block_logical_coordinates
std::vector< block_logical_coord_holder< Dim > > block_logical_coordinates(const Domain< Dim > &domain, const tnsr::I< DataVector, Dim, Frame::Inertial > &x) noexcept
Definition: BlockLogicalCoordinates.cpp:25
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:81
InterpolateOnElementTestHelpers
Holds code that is shared between multiple tests. Currently used by.
Definition: InterpolateOnElementTestHelpers.hpp:37
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:136
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
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:1588
Mesh
Holds the number of grid points, basis, and quadrature in each direction of the computational grid.
Definition: Mesh.hpp:50
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
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
InterpolateOnElementTestHelpers::Tags::MultiplyByTwoComputeItem
Definition: InterpolateOnElementTestHelpers.hpp:52
Rational
Definition: Rational.hpp:24
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:111