InterpolationTargetWedgeSectionTorus.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <cmath>
7 #include <cstddef>
8 
9 #include "DataStructures/DataVector.hpp"
11 #include "NumericalAlgorithms/Interpolation/SendPointsToInterpolator.hpp"
13 #include "Options/Options.hpp"
15 #include "Utilities/Requires.hpp"
16 #include "Utilities/TMPL.hpp"
18 
19 namespace PUP {
20 class er;
21 } // namespace PUP
22 namespace db {
23 template <typename TagsList>
24 class DataBox;
25 } // namespace db
26 namespace intrp {
27 namespace Tags {
28 template <typename Metavariables>
29 struct TemporalIds;
30 } // namespace Tags
31 } // namespace intrp
32 
33 namespace intrp {
34 
35 namespace OptionHolders {
36 /// \brief A solid torus of points, useful, e.g., when measuring data from an
37 /// accretion disc.
38 ///
39 /// The torus's cross section (e.g., a cut at \f$\phi=0\f$) is a wedge-like
40 /// shape bounded by \f$r_{\text{min}} \le r \le r_{\text{max}}\f$ and
41 /// \f$\theta_{\text{min}} \le \theta \le \theta_{\text{max}}\f$.
42 ///
43 /// The grid points are located on surfaces of constant \f$r\f$, \f$\theta\f$,
44 /// and \f$\phi\f$. There are `NumberRadialPoints` points in the radial
45 /// direction between `MinRadius` and `MaxRadius` (including these endpoints);
46 /// `NumberThetaPoints` points in the \f$\theta\f$ direction between `MinTheta`
47 /// and `MaxTheta` (including these endpoints); `NumberPhiPoints` points in the
48 /// \f$\phi\f$ direction (with one point always at \f$\phi=0\f$).
49 ///
50 /// By default, the points follow a Legendre Gauss-Lobatto distribution in the
51 /// \f$r\f$ and \f$\theta\f$ directions, and a uniform distribution in the
52 /// \f$\phi\f$ direction. The distribution in the \f$r\f$ (and/or \f$\theta\f$)
53 /// direction can be made uniform using the `UniformRadialGrid` (and/or
54 /// `UniformThetaGrid`) option.
55 ///
56 /// The `target_points` form a 3D mesh ordered with \f$r\f$ varying fastest,
57 /// then \f$\theta\f$, and finally \f$\phi\f$ varying slowest.
58 ///
59 /// \note Input coordinates (radii, angles) are interpreted in the frame given
60 /// by `Metavariables::domain_frame`
62  struct MinRadius {
63  using type = double;
64  static constexpr OptionString help = {"Inner radius of torus"};
65  static type lower_bound() noexcept { return 0.0; }
66  };
67  struct MaxRadius {
68  using type = double;
69  static constexpr OptionString help = {"Outer radius of torus"};
70  static type lower_bound() noexcept { return 0.0; }
71  };
72  struct MinTheta {
73  using type = double;
74  static constexpr OptionString help = {"Angle of top of wedge (radians)"};
75  static type lower_bound() noexcept { return 0.0; }
76  static type upper_bound() noexcept { return M_PI; }
77  };
78  struct MaxTheta {
79  using type = double;
80  static constexpr OptionString help = {"Angle of bottom of wedge (radians)"};
81  static type lower_bound() noexcept { return 0.0; }
82  static type upper_bound() noexcept { return M_PI; }
83  };
85  using type = size_t;
86  static constexpr OptionString help = {
87  "Number of radial points, including endpoints"};
88  static type lower_bound() noexcept { return 2; }
89  };
91  using type = size_t;
92  static constexpr OptionString help = {
93  "Number of theta points, including endpoints"};
94  static type lower_bound() noexcept { return 2; }
95  };
96  struct NumberPhiPoints {
97  using type = size_t;
98  static constexpr OptionString help = {"Number of phi points"};
99  static type lower_bound() noexcept { return 1; }
100  };
102  using type = bool;
103  static constexpr OptionString help = {
104  "Use uniform radial grid [default: LGL grid]"};
105  static type default_value() noexcept { return false; }
106  };
108  using type = bool;
109  static constexpr OptionString help = {
110  "Use uniform theta grid [default: LGL grid]"};
111  static type default_value() noexcept { return false; }
112  };
113 
114  using options =
118  static constexpr OptionString help = {
119  "A torus extending from MinRadius to MaxRadius in r, MinTheta to MaxTheta"
120  " in theta, and 2pi in phi."};
121 
122  WedgeSectionTorus(double min_radius_in, double max_radius_in,
123  double min_theta_in, double max_theta_in,
124  size_t number_of_radial_points_in,
125  size_t number_of_theta_points_in,
126  size_t number_of_phi_points_in,
127  bool use_uniform_radial_grid_in,
128  bool use_uniform_theta_grid_in,
129  const OptionContext& context = {});
130 
131  WedgeSectionTorus() = default;
132  WedgeSectionTorus(const WedgeSectionTorus& /*rhs*/) = default;
133  WedgeSectionTorus& operator=(const WedgeSectionTorus& /*rhs*/) = default;
134  WedgeSectionTorus(WedgeSectionTorus&& /*rhs*/) noexcept = default;
135  WedgeSectionTorus& operator=(WedgeSectionTorus&& /*rhs*/) noexcept = default;
136  ~WedgeSectionTorus() = default;
137 
138  // clang-tidy non-const reference pointer
139  void pup(PUP::er& p) noexcept; // NOLINT
140 
141  double min_radius;
142  double max_radius;
143  double min_theta;
144  double max_theta;
145  size_t number_of_radial_points;
146  size_t number_of_theta_points;
147  size_t number_of_phi_points;
148  bool use_uniform_radial_grid;
149  bool use_uniform_theta_grid;
150 };
151 
152 bool operator==(const WedgeSectionTorus& lhs,
153  const WedgeSectionTorus& rhs) noexcept;
154 bool operator!=(const WedgeSectionTorus& lhs,
155  const WedgeSectionTorus& rhs) noexcept;
156 
157 } // namespace OptionHolders
158 
159 namespace Actions {
160 /// \ingroup ActionsGroup
161 /// \brief Sends points in a wedge-sectioned torus to an `Interpolator`.
162 ///
163 /// Uses:
164 /// - DataBox:
165 /// - `::Tags::Domain<3, Frame>`
166 /// - `::Tags::Variables<typename
167 /// InterpolationTargetTag::vars_to_interpolate_to_target>`
168 ///
169 /// DataBox changes:
170 /// - Adds: nothing
171 /// - Removes: nothing
172 /// - Modifies:
173 /// - `Tags::IndicesOfFilledInterpPoints`
174 /// - `::Tags::Variables<typename
175 /// InterpolationTargetTag::vars_to_interpolate_to_target>`
176 ///
177 /// For requirements on InterpolationTargetTag, see InterpolationTarget
178 template <typename InterpolationTargetTag>
181  using const_global_cache_tags = tmpl::list<InterpolationTargetTag>;
182  template <typename DbTags, typename... InboxTags, typename Metavariables,
183  typename ArrayIndex, typename ActionList,
184  typename ParallelComponent,
185  Requires<tmpl::list_contains_v<
186  DbTags, typename Tags::TemporalIds<Metavariables>>> = nullptr>
187  static void apply(
188  db::DataBox<DbTags>& box,
189  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
191  const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
192  const ParallelComponent* const /*meta*/,
193  const typename Metavariables::temporal_id::type& temporal_id) noexcept {
194  const auto& options = Parallel::get<InterpolationTargetTag>(cache);
195 
196  // Compute locations of constant r/theta/phi surfaces
197  const size_t num_radial = options.number_of_radial_points;
198  const DataVector radii_1d = [&num_radial, &options ]() noexcept {
199  DataVector result(num_radial);
200  if (options.use_uniform_radial_grid) {
201  // uniform point distribution
202  for (size_t r = 0; r < num_radial; ++r) {
203  result[r] =
204  options.min_radius + (options.max_radius - options.min_radius) *
205  r / (num_radial - 1.0);
206  }
207  } else {
208  // Legendre Gauss-Lobatto point distribution
209  const double mean = 0.5 * (options.max_radius + options.min_radius);
210  const double diff = 0.5 * (options.max_radius - options.min_radius);
211  result =
212  mean + diff * Spectral::collocation_points<
213  Spectral::Basis::Legendre,
214  Spectral::Quadrature::GaussLobatto>(num_radial);
215  }
216  return result;
217  }
218  ();
219  const size_t num_theta = options.number_of_theta_points;
220  const DataVector thetas_1d = [&num_theta, &options ]() noexcept {
221  DataVector result(num_theta);
222  if (options.use_uniform_theta_grid) {
223  // uniform point distribution
224  for (size_t theta = 0; theta < num_theta; ++theta) {
225  result[theta] =
226  options.min_theta + (options.max_theta - options.min_theta) *
227  theta / (num_theta - 1.0);
228  }
229  } else {
230  // Legendre Gauss-Lobatto point distribution
231  const double mean = 0.5 * (options.max_theta + options.min_theta);
232  const double diff = 0.5 * (options.max_theta - options.min_theta);
233  result =
234  mean + diff * Spectral::collocation_points<
235  Spectral::Basis::Legendre,
236  Spectral::Quadrature::GaussLobatto>(num_theta);
237  }
238  return result;
239  }
240  ();
241  const size_t num_phi = options.number_of_phi_points;
242  const DataVector phis_1d = [&num_phi]() noexcept {
243  DataVector result(num_phi);
244  for (size_t phi = 0; phi < num_phi; ++phi) {
245  // We do NOT want a grid point at phi = 2pi, as this would duplicate the
246  // phi = 0 data. So, divide by num_phi rather than (n-1) as elsewhere.
247  result[phi] = 2.0 * M_PI * phi / num_phi;
248  }
249  return result;
250  }
251  ();
252 
253  // Take tensor product to get full 3D r/theta/phi points
254  const size_t num_total = num_radial * num_theta * num_phi;
255  DataVector radii(num_total), thetas(num_total), phis(num_total);
256  for (size_t phi = 0; phi < num_phi; ++phi) {
257  for (size_t theta = 0; theta < num_theta; ++theta) {
258  for (size_t r = 0; r < num_radial; ++r) {
259  const size_t i =
260  r + theta * num_radial + phi * num_theta * num_radial;
261  radii[i] = radii_1d[r];
262  thetas[i] = thetas_1d[theta];
263  phis[i] = phis_1d[phi];
264  }
265  }
266  }
267 
268  // Compute x/y/z coordinates
269  // Note: theta measured from +z axis, phi measured from +x axis
270  tnsr::I<DataVector, 3, typename Metavariables::domain_frame> target_points(
271  num_total);
272  get<0>(target_points) = radii * sin(thetas) * cos(phis);
273  get<1>(target_points) = radii * sin(thetas) * sin(phis);
274  get<2>(target_points) = radii * cos(thetas);
275 
276  send_points_to_interpolator<InterpolationTargetTag>(
277  box, cache, target_points, temporal_id);
278  }
279 };
280 
281 } // namespace Actions
282 } // namespace intrp
Definition: InterpolationTargetWedgeSectionTorus.hpp:107
Definition: InterpolationTargetWedgeSectionTorus.hpp:90
Definition: Strahlkorper.hpp:14
Defines class tuples::TaggedTuple.
Sends points in a wedge-sectioned torus to an Interpolator.
Definition: InterpolationTargetWedgeSectionTorus.hpp:179
const DataVector & collocation_points(const size_t num_points) noexcept
Collocation points.
Definition: Spectral.cpp:232
Definition: AddTemporalIdsToInterpolationTarget.hpp:17
Defines classes and functions for making classes creatable from input files.
Definition: InterpolationTargetWedgeSectionTorus.hpp:101
Defines the type alias Requires.
constexpr auto apply(F &&f, const DataBox< BoxTags > &box, Args &&... args)
Apply the function f with argument Tags TagsList from DataBox box
Definition: DataBox.hpp:1595
Definition: InterpolationTargetWedgeSectionTorus.hpp:72
const char *const OptionString
The string used in option structs.
Definition: Options.hpp:26
Definition: InterpolationTargetWedgeSectionTorus.hpp:84
temporal_ids on which to interpolate.
Definition: InterpolationTargetWedgeSectionTorus.hpp:29
Definition: InterpolationTargetWedgeSectionTorus.hpp:96
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:272
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&#39;s equations...
Definition: ComputeGhQuantities.cpp:22
Definition: InterpolationTargetWedgeSectionTorus.hpp:24
Definition: DataBoxTag.hpp:29
A Charm++ chare that caches constant data once per Charm++ node.
Definition: ConstGlobalCache.hpp:76
Definition: InterpolationTargetWedgeSectionTorus.hpp:62
Defines a list of useful type aliases for tensors.
Definition: InterpolationTargetWedgeSectionTorus.hpp:67
Namespace for DataBox related things.
Definition: DataBox.hpp:33
Stores a collection of function values.
Definition: DataVector.hpp:46
Information about the nested operations being performed by the parser, for use in printing errors...
Definition: Options.hpp:35
Wraps the template metaprogramming library used (brigand)
A solid torus of points, useful, e.g., when measuring data from an accretion disc.
Definition: InterpolationTargetWedgeSectionTorus.hpp:61
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
Identifies a step in the linear solver algorithm.
Definition: IterationId.hpp:25
Definition: SolvePoissonProblem.hpp:38
Defines class template ConstGlobalCache.
Definition: InterpolationTargetWedgeSectionTorus.hpp:78
Definition: ComputeTimeDerivative.hpp:28