HalfSpaceMirror.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 
9 #include "DataStructures/CachedTempBuffer.hpp"
12 #include "Elliptic/Systems/Elasticity/Tags.hpp"
13 #include "Options/Options.hpp"
15 #include "PointwiseFunctions/AnalyticSolutions/Elasticity/AnalyticSolution.hpp"
16 #include "PointwiseFunctions/Elasticity/ConstitutiveRelations/IsotropicHomogeneous.hpp"
17 #include "Utilities/ContainerHelpers.hpp"
19 #include "Utilities/Registration.hpp"
20 #include "Utilities/TMPL.hpp"
21 #include "Utilities/TaggedTuple.hpp"
22 
23 namespace Elasticity::Solutions {
24 
25 namespace detail {
26 template <typename DataType>
27 struct HalfSpaceMirrorVariables {
28  struct DisplacementR {
29  using type = Scalar<DataType>;
30  };
31  using Cache =
32  CachedTempBuffer<HalfSpaceMirrorVariables, DisplacementR,
33  Tags::Displacement<3>, Tags::Strain<3>,
34  Tags::MinusStress<3>, Tags::PotentialEnergyDensity<3>,
36 
37  const tnsr::I<DataType, 3>& x;
38  const double beam_width;
39  const ConstitutiveRelations::IsotropicHomogeneous<3>& constitutive_relation;
40  const size_t integration_intervals;
41  const double absolute_tolerance;
42  const double relative_tolerance;
43 
44  void operator()(gsl::not_null<Scalar<DataType>*> displacement_r,
46  DisplacementR /*meta*/) const noexcept;
47  void operator()(gsl::not_null<tnsr::I<DataType, 3>*> displacement,
49  Tags::Displacement<3> /*meta*/) const noexcept;
50  void operator()(gsl::not_null<tnsr::ii<DataType, 3>*> strain,
52  Tags::Strain<3> /*meta*/) const noexcept;
53  void operator()(gsl::not_null<tnsr::II<DataType, 3>*> minus_stress,
55  Tags::MinusStress<3> /*meta*/) const noexcept;
58  Tags::PotentialEnergyDensity<3> /*meta*/) const noexcept;
59  void operator()(
60  gsl::not_null<tnsr::I<DataType, 3>*> fixed_source_for_displacement,
62  ::Tags::FixedSource<Tags::Displacement<3>> /*meta*/) const noexcept;
63 };
64 } // namespace detail
65 
66 /// \cond
67 template <typename Registrars>
68 struct HalfSpaceMirror;
69 
70 namespace Registrars {
72 } // namespace Registrars
73 /// \endcond
74 
75 /*!
76  * \brief The solution for a half-space mirror deformed by a laser beam.
77  *
78  * \details This solution is mapping (via the fluctuation dissipation theorem)
79  * thermal noise to an elasticity problem where a normally incident and
80  * axisymmetric laser beam with a Gaussian beam profile acts on the face of a
81  * semi-infinite mirror. Here we assume the face to be at \f$z = 0\f$ and the
82  * material to extend to \f$+\infty\f$ in the z-direction as well as for the
83  * mirror diameter to be comparatively large to the `beam width`. The mirror
84  * material is characterized by an isotropic homogeneous constitutive relation
85  * \f$Y^{ijkl}\f$ (see
86  * `Elasticity::ConstitutiveRelations::IsotropicHomogeneous`). In this scenario,
87  * the auxiliary elastic problem has an applied pressure distribution equal to
88  * the laser beam intensity profile \f$p(r)\f$ (see Eq. (11.94) and Eq. (11.95)
89  * in \cite ThorneBlandford2017 with F = 1 and the time dependency dropped)
90  *
91  * \f{align}
92  * T^{zr} &= T^{rz} = 0 \\
93  * T^{zz} &= p(r) = \frac{e^{-\frac{r^2}{r_0^2}}}{\pi r_0^2}\text{.}
94  * \f}
95  *
96  * in the form of a Neumann boundary condition to the face of the mirror. We
97  * find that this stress in cylinder coordinates is produced by the displacement
98  * field
99  *
100  * \f{align}
101  * \xi_{r} &= \frac{1}{2 \mu} \int_0^{\infty} dk J_1(kr)e^{(-kz)}\left(1 -
102  * \frac{\lambda + 2\mu}{\lambda + \mu} + kz \right) \tilde{p}(k) \\
103  * \xi_{\phi} &= 0 \\
104  * \xi_{z} &= \frac{1}{2 \mu} \int_0^{\infty} dk J_0(kr)e^{(-kz)}\left(1 +
105  * \frac{\mu}{\lambda + \mu} + kz \right) \tilde{p}(k)
106  * \f}
107  *
108  * and the strain
109  *
110  * \f{align}
111  * \Theta &= \frac{1}{2 \mu} \int_0^{\infty} dk
112  * J_0(kr) k e^{(-kz)}\left(\frac{-2\mu}{\lambda + \mu}\right) \tilde{p}(k) \\
113  * S_{rr} &= \Theta - S_{\phi\phi} - S_{zz} \\
114  * S_{\phi\phi} &= \frac{\xi_{r}}{r} \\
115  * S_{(rz)} &= -\frac{1}{2 \mu} \int_0^{\infty} dk J_1(kr) k e^{(-kz)}\left(kz
116  * \right) \tilde{p}(k) \\
117  * S_{zz} &= \frac{1}{2 \mu} \int_0^{\infty} dk
118  * J_0(kr) k e^{(-kz)}\left(-\frac{\mu}{\lambda + \mu} - kz \right) \tilde{p}(k)
119  * \f}
120  *
121  * (see Eqs. (11 a) - (11 c) and (13 a) - (13 e), with (13 c) swapped in favor
122  * of (12 c) in \cite Lovelace2007tn), where \f$\tilde{p}(k)= \frac{1}{2\pi}
123  * e^{-(\frac{kr_0}{2})^2}\f$ is the Hankel-Transform of the lasers intensity
124  * profile and \f$ \Theta = \mathrm{Tr}(S)\f$ the materials expansion.
125  *
126  */
127 template <typename Registrars =
128  tmpl::list<Solutions::Registrars::HalfSpaceMirror>>
129 class HalfSpaceMirror : public AnalyticSolution<3, Registrars> {
130  public:
133 
134  struct BeamWidth {
135  using type = double;
136  static constexpr Options::String help{
137  "The lasers beam width r_0 with FWHM = 2*sqrt(ln 2)*r_0"};
138  static type lower_bound() noexcept { return 0.0; }
139  };
140 
141  struct Material {
143  static constexpr Options::String help{
144  "The material properties of the beam"};
145  };
146 
148  using type = size_t;
149  static constexpr Options::String help{
150  "Workspace size for numerical integrals. Increase if integrals fail to "
151  "reach the prescribed tolerance at large distances relative to the "
152  "beam width. The suggested values for workspace size and tolerances "
153  "should accommodate distances of up to ~100 beam widths."};
154  static type lower_bound() noexcept { return 1; }
155  static type suggested_value() noexcept { return 350; }
156  };
157 
159  using type = double;
160  static constexpr Options::String help{
161  "Absolute tolerance for numerical integrals"};
162  static type lower_bound() noexcept { return 0.; }
163  static type suggested_value() noexcept { return 1e-12; }
164  };
165 
167  using type = double;
168  static constexpr Options::String help{
169  "Relative tolerance for numerical integrals"};
170  static type lower_bound() noexcept { return 0.; }
171  static type upper_bound() noexcept { return 1.; }
172  static type suggested_value() noexcept { return 1e-10; }
173  };
174 
175  using options = tmpl::list<BeamWidth, Material, IntegrationIntervals,
177  static constexpr Options::String help{
178  "A semi-infinite mirror on which a laser introduces stress perpendicular "
179  "to the mirrors surface."};
180 
181  HalfSpaceMirror() = default;
182  HalfSpaceMirror(const HalfSpaceMirror&) noexcept = default;
183  HalfSpaceMirror& operator=(const HalfSpaceMirror&) noexcept = default;
184  HalfSpaceMirror(HalfSpaceMirror&&) noexcept = default;
185  HalfSpaceMirror& operator=(HalfSpaceMirror&&) noexcept = default;
186  ~HalfSpaceMirror() noexcept override = default;
187 
188  /// \cond
189  explicit HalfSpaceMirror(CkMigrateMessage* /*unused*/) noexcept {}
190  using PUP::able::register_constructor;
191  WRAPPED_PUPable_decl_template(HalfSpaceMirror); // NOLINT
192  /// \endcond
193 
194  HalfSpaceMirror(double beam_width,
195  constitutive_relation_type constitutive_relation,
196  size_t integration_intervals = 350,
197  double absolute_tolerance = 1e-12,
198  double relative_tolerance = 1e-10) noexcept
199  : beam_width_(beam_width),
200  constitutive_relation_(std::move(constitutive_relation)),
201  integration_intervals_(integration_intervals),
202  absolute_tolerance_(absolute_tolerance),
203  relative_tolerance_(relative_tolerance) {}
204 
205  double beam_width() const noexcept { return beam_width_; }
206  size_t integration_intervals() const noexcept {
207  return integration_intervals_;
208  }
209  double absolute_tolerance() const noexcept { return absolute_tolerance_; }
210  double relative_tolerance() const noexcept { return relative_tolerance_; }
211 
212  const constitutive_relation_type& constitutive_relation() const
213  noexcept override {
214  return constitutive_relation_;
215  }
216 
217  template <typename DataType, typename... RequestedTags>
218  tuples::TaggedTuple<RequestedTags...> variables(
219  const tnsr::I<DataType, 3>& x,
220  tmpl::list<RequestedTags...> /*meta*/) const noexcept {
221  for (size_t i = 0; i < get_size(get<2>(x)); i++) {
222  if (UNLIKELY(get_element(get<2>(x), i) < 0)) {
223  ERROR(
224  "The HalfSpaceMirror solution is not defined for negative values "
225  "of z.");
226  }
227  }
228  using VarsComputer = detail::HalfSpaceMirrorVariables<DataType>;
229  typename VarsComputer::Cache cache{
230  get_size(*x.begin()),
231  VarsComputer{x, beam_width_, constitutive_relation_,
232  integration_intervals_, absolute_tolerance_,
233  relative_tolerance_}};
234  return {cache.get_var(RequestedTags{})...};
235  }
236 
237  // clang-tidy: no pass by reference
238  void pup(PUP::er& p) noexcept override { // NOLINT
239  p | beam_width_;
240  p | constitutive_relation_;
241  p | integration_intervals_;
242  p | absolute_tolerance_;
243  p | relative_tolerance_;
244  }
245 
246  private:
247  double beam_width_{std::numeric_limits<double>::signaling_NaN()};
248  constitutive_relation_type constitutive_relation_{};
249  size_t integration_intervals_{std::numeric_limits<size_t>::max()};
250  double absolute_tolerance_{std::numeric_limits<double>::signaling_NaN()};
251  double relative_tolerance_{std::numeric_limits<double>::signaling_NaN()};
252 };
253 
254 /// \cond
255 template <typename Registrars>
256 PUP::able::PUP_ID HalfSpaceMirror<Registrars>::my_PUP_ID = 0; // NOLINT
257 /// \endcond
258 
259 template <typename Registrars>
260 bool operator==(const HalfSpaceMirror<Registrars>& lhs,
261  const HalfSpaceMirror<Registrars>& rhs) noexcept {
262  return lhs.beam_width() == rhs.beam_width() and
263  lhs.constitutive_relation() == rhs.constitutive_relation() and
264  lhs.integration_intervals() == rhs.integration_intervals() and
265  lhs.absolute_tolerance() == rhs.absolute_tolerance() and
266  lhs.relative_tolerance() == rhs.relative_tolerance();
267 }
268 
269 template <typename Registrars>
270 bool operator!=(const HalfSpaceMirror<Registrars>& lhs,
271  const HalfSpaceMirror<Registrars>& rhs) noexcept {
272  return not(lhs == rhs);
273 }
274 
275 } // namespace Elasticity::Solutions
Elasticity::Solutions::AnalyticSolution
Base class for analytic solutions of the linear Elasticity equations.
Definition: AnalyticSolution.hpp:34
CharmPupable.hpp
UNLIKELY
#define UNLIKELY(x)
Definition: Gsl.hpp:73
std::rel_ops::operator!=
T operator!=(T... args)
Options.hpp
Elasticity::strain
void strain(gsl::not_null< tnsr::ii< DataType, Dim > * > strain, const tnsr::iJ< DataType, Dim > &deriv_displacement) noexcept
The symmetric strain on a flat background in Cartesian coordinates.
Error.hpp
Elasticity::Solutions::HalfSpaceMirror::BeamWidth
Definition: HalfSpaceMirror.hpp:134
Elasticity::Solutions::HalfSpaceMirror::AbsoluteTolerance
Definition: HalfSpaceMirror.hpp:158
get_size
decltype(auto) get_size(const T &t, SizeFunction size=GetContainerSize{}) noexcept
Retrieve the size of t if t.size() is a valid expression, otherwise if T is fundamental or a std::com...
Definition: ContainerHelpers.hpp:145
Registration::Registrar
A template for defining a registrar.
Definition: Registration.hpp:42
Tags::FixedSource
Prefix indicating a source term that is independent of dynamic variables.
Definition: Prefixes.hpp:75
ERROR
#define ERROR(m)
prints an error message to the standard error stream and aborts the program.
Definition: Error.hpp:37
WRAPPED_PUPable_decl_template
#define WRAPPED_PUPable_decl_template(className)
Mark derived classes as serializable.
Definition: CharmPupable.hpp:22
Elasticity::Solutions::HalfSpaceMirror::RelativeTolerance
Definition: HalfSpaceMirror.hpp:166
cstddef
CachedTempBuffer
Definition: CachedTempBuffer.hpp:29
tuples::TaggedTuple
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:271
Elasticity::Solutions::HalfSpaceMirror
The solution for a half-space mirror deformed by a laser beam.
Definition: HalfSpaceMirror.hpp:129
std::numeric_limits::signaling_NaN
T signaling_NaN(T... args)
Elasticity::Solutions
Analytic solutions of the linear Elasticity equations.
Definition: AnalyticSolution.hpp:18
Elasticity::potential_energy_density
void potential_energy_density(gsl::not_null< Scalar< DataVector > * > potential_energy_density, const tnsr::ii< DataVector, Dim > &strain, const tnsr::I< DataVector, Dim > &coordinates, const ConstitutiveRelations::ConstitutiveRelation< Dim > &constitutive_relation) noexcept
The potential energy density stored in the deformation of the elastic material (see Eq....
ActionTesting::cache
Parallel::GlobalCache< Metavariables > & cache(MockRuntimeSystem< Metavariables > &runner, const ArrayIndex &array_index) noexcept
Returns the GlobalCache of Component with index array_index.
Definition: MockRuntimeSystemFreeFunctions.hpp:382
Scalar
Tensor< T, Symmetry<>, index_list<> > Scalar
Definition: TypeAliases.hpp:21
Elasticity::Solutions::HalfSpaceMirror::Material
Definition: HalfSpaceMirror.hpp:141
Elasticity::ConstitutiveRelations::IsotropicHomogeneous< 3 >
limits
Elasticity::Solutions::HalfSpaceMirror::IntegrationIntervals
Definition: HalfSpaceMirror.hpp:147
Options::String
const char *const String
The string used in option structs.
Definition: Options.hpp:32
Tensor.hpp
std::numeric_limits::max
T max(T... args)
Prefixes.hpp
get_element
decltype(auto) get_element(T &t, const size_t i, SubscriptFunction at=GetContainerElement{}) noexcept
Returns the ith element if T has a subscript operator, otherwise if T is fundamental or a std::comple...
Definition: ContainerHelpers.hpp:117
TMPL.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecPiecewisePolynomial.hpp:13