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  using Cache =
30  Tags::Strain<3>, Tags::MinusStress<3>,
31  Tags::PotentialEnergyDensity<3>,
33 
34  const tnsr::I<DataType, 3>& x;
35  const double beam_width;
36  const ConstitutiveRelations::IsotropicHomogeneous<3>& constitutive_relation;
37  const size_t integration_intervals;
38  const double absolute_tolerance;
39  const double relative_tolerance;
40 
41  void operator()(gsl::not_null<tnsr::I<DataType, 3>*> displacement,
43  Tags::Displacement<3> /*meta*/) const noexcept;
44  void operator()(gsl::not_null<tnsr::ii<DataType, 3>*> strain,
46  Tags::Strain<3> /*meta*/) const noexcept;
47  void operator()(gsl::not_null<tnsr::II<DataType, 3>*> minus_stress,
49  Tags::MinusStress<3> /*meta*/) const noexcept;
52  Tags::PotentialEnergyDensity<3> /*meta*/) const noexcept;
53  void operator()(
54  gsl::not_null<tnsr::I<DataType, 3>*> fixed_source_for_displacement,
56  ::Tags::FixedSource<Tags::Displacement<3>> /*meta*/) const noexcept;
57 };
58 } // namespace detail
59 
60 /// \cond
61 template <typename Registrars>
62 struct HalfSpaceMirror;
63 
64 namespace Registrars {
66 } // namespace Registrars
67 /// \endcond
68 
69 /*!
70  * \brief The solution for a half-space mirror deformed by a laser beam.
71  *
72  * \details This solution is mapping (via the fluctuation dissipation theorem)
73  * thermal noise to an elasticity problem where a normally incident and
74  * axisymmetric laser beam with a Gaussian beam profile acts on the face of a
75  * semi-infinite mirror. Here we assume the face to be at \f$z = 0\f$ and the
76  * material to extend to \f$+\infty\f$ in the z-direction as well as for the
77  * mirror diameter to be comparatively large to the `beam width`. The mirror
78  * material is characterized by an isotropic homogeneous constitutive relation
79  * \f$Y^{ijkl}\f$ (see
80  * `Elasticity::ConstitutiveRelations::IsotropicHomogeneous`). In this scenario,
81  * the auxiliary elastic problem has an applied pressure distribution equal to
82  * the laser beam intensity profile \f$p(r)\f$ (see Eq. (11.94) and Eq. (11.95)
83  * in \cite ThorneBlandford2017 with F = 1 and the time dependency dropped)
84  *
85  * \f{align}
86  * T^{zr} &= T^{rz} = 0 \\
87  * T^{zz} &= p(r) = \frac{e^{-\frac{r^2}{r_0^2}}}{\pi r_0^2}\text{.}
88  * \f}
89  *
90  * in the form of a Neumann boundary condition to the face of the mirror. We
91  * find that this stress in cylinder coordinates is produced by the displacement
92  * field
93  *
94  * \f{align}
95  * \xi_{r} &= \frac{1}{2 \mu} \int_0^{\infty} dk J_1(kr)e^{(-kz)}\left(1 -
96  * \frac{\lambda + 2\mu}{\lambda + \mu} + kz \right) \tilde{p}(k) \\
97  * \xi_{\phi} &= 0 \\
98  * \xi_{z} &= \frac{1}{2 \mu} \int_0^{\infty} dk J_0(kr)e^{(-kz)}\left(1 +
99  * \frac{\mu}{\lambda + \mu} + kz \right) \tilde{p}(k)
100  * \f}
101  *
102  * and the strain
103  *
104  * \f{align}
105  * \Theta &= \frac{1}{2 \mu} \int_0^{\infty} dk
106  * J_0(kr) k e^{(-kz)}\left(\frac{-2\mu}{\lambda + \mu}\right) \tilde{p}(k) \\
107  * S_{rr} &= \Theta - S_{\phi\phi} - S_{zz} \\
108  * S_{\phi\phi} &= \frac{\xi_{r}}{r} \\
109  * S_{(rz)} &= -\frac{1}{2 \mu} \int_0^{\infty} dk J_1(kr) k e^{(-kz)}\left(kz
110  * \right) \tilde{p}(k) \\
111  * S_{zz} &= \frac{1}{2 \mu} \int_0^{\infty} dk
112  * J_0(kr) k e^{(-kz)}\left(-\frac{\mu}{\lambda + \mu} - kz \right) \tilde{p}(k)
113  * \f}
114  *
115  * (see Eqs. (11 a) - (11 c) and (13 a) - (13 e), with (13 c) swapped in favor
116  * of (12 c) in \cite Lovelace2007tn), where \f$\tilde{p}(k)= \frac{1}{2\pi}
117  * e^{-(\frac{kr_0}{2})^2}\f$ is the Hankel-Transform of the lasers intensity
118  * profile and \f$ \Theta = \mathrm{Tr}(S)\f$ the materials expansion.
119  *
120  */
121 template <typename Registrars =
122  tmpl::list<Solutions::Registrars::HalfSpaceMirror>>
123 class HalfSpaceMirror : public AnalyticSolution<3, Registrars> {
124  public:
127 
128  struct BeamWidth {
129  using type = double;
130  static constexpr Options::String help{
131  "The lasers beam width r_0 with FWHM = 2*sqrt(ln 2)*r_0"};
132  static type lower_bound() noexcept { return 0.0; }
133  };
134 
135  struct Material {
137  static constexpr Options::String help{
138  "The material properties of the beam"};
139  };
140 
142  using type = size_t;
143  static constexpr Options::String help{
144  "Workspace size for numerical integrals. Increase if integrals fail to "
145  "reach the prescribed tolerance at large distances relative to the "
146  "beam width. The suggested values for workspace size and tolerances "
147  "should accommodate distances of up to ~100 beam widths."};
148  static type lower_bound() noexcept { return 1; }
149  static type suggested_value() noexcept { return 350; }
150  };
151 
153  using type = double;
154  static constexpr Options::String help{
155  "Absolute tolerance for numerical integrals"};
156  static type lower_bound() noexcept { return 0.; }
157  static type suggested_value() noexcept { return 1e-12; }
158  };
159 
161  using type = double;
162  static constexpr Options::String help{
163  "Relative tolerance for numerical integrals"};
164  static type lower_bound() noexcept { return 0.; }
165  static type upper_bound() noexcept { return 1.; }
166  static type suggested_value() noexcept { return 1e-10; }
167  };
168 
169  using options = tmpl::list<BeamWidth, Material, IntegrationIntervals,
171  static constexpr Options::String help{
172  "A semi-infinite mirror on which a laser introduces stress perpendicular "
173  "to the mirrors surface."};
174 
175  HalfSpaceMirror() = default;
176  HalfSpaceMirror(const HalfSpaceMirror&) noexcept = default;
177  HalfSpaceMirror& operator=(const HalfSpaceMirror&) noexcept = default;
178  HalfSpaceMirror(HalfSpaceMirror&&) noexcept = default;
179  HalfSpaceMirror& operator=(HalfSpaceMirror&&) noexcept = default;
180  ~HalfSpaceMirror() noexcept override = default;
181 
182  /// \cond
183  explicit HalfSpaceMirror(CkMigrateMessage* /*unused*/) noexcept {}
184  using PUP::able::register_constructor;
185  WRAPPED_PUPable_decl_template(HalfSpaceMirror); // NOLINT
186  /// \endcond
187 
188  HalfSpaceMirror(double beam_width,
189  constitutive_relation_type constitutive_relation,
190  size_t integration_intervals = 350,
191  double absolute_tolerance = 1e-12,
192  double relative_tolerance = 1e-10) noexcept
193  : beam_width_(beam_width),
194  constitutive_relation_(std::move(constitutive_relation)),
195  integration_intervals_(integration_intervals),
196  absolute_tolerance_(absolute_tolerance),
197  relative_tolerance_(relative_tolerance) {}
198 
199  double beam_width() const noexcept { return beam_width_; }
200  size_t integration_intervals() const noexcept {
201  return integration_intervals_;
202  }
203  double absolute_tolerance() const noexcept { return absolute_tolerance_; }
204  double relative_tolerance() const noexcept { return relative_tolerance_; }
205 
206  const constitutive_relation_type& constitutive_relation() const
207  noexcept override {
208  return constitutive_relation_;
209  }
210 
211  template <typename DataType, typename... RequestedTags>
212  tuples::TaggedTuple<RequestedTags...> variables(
213  const tnsr::I<DataType, 3>& x,
214  tmpl::list<RequestedTags...> /*meta*/) const noexcept {
215  for (size_t i = 0; i < get_size(get<2>(x)); i++) {
216  if (UNLIKELY(get_element(get<2>(x), i) < 0)) {
217  ERROR(
218  "The HalfSpaceMirror solution is not defined for negative values "
219  "of z.");
220  }
221  }
222  using VarsComputer = detail::HalfSpaceMirrorVariables<DataType>;
223  typename VarsComputer::Cache cache{
224  get_size(*x.begin()),
225  VarsComputer{x, beam_width_, constitutive_relation_,
226  integration_intervals_, absolute_tolerance_,
227  relative_tolerance_}};
228  return {cache.get_var(RequestedTags{})...};
229  }
230 
231  // clang-tidy: no pass by reference
232  void pup(PUP::er& p) noexcept override { // NOLINT
233  p | beam_width_;
234  p | constitutive_relation_;
235  p | integration_intervals_;
236  p | absolute_tolerance_;
237  p | relative_tolerance_;
238  }
239 
240  private:
241  double beam_width_{std::numeric_limits<double>::signaling_NaN()};
242  constitutive_relation_type constitutive_relation_{};
243  size_t integration_intervals_{std::numeric_limits<size_t>::max()};
244  double absolute_tolerance_{std::numeric_limits<double>::signaling_NaN()};
245  double relative_tolerance_{std::numeric_limits<double>::signaling_NaN()};
246 };
247 
248 /// \cond
249 template <typename Registrars>
250 PUP::able::PUP_ID HalfSpaceMirror<Registrars>::my_PUP_ID = 0; // NOLINT
251 /// \endcond
252 
253 template <typename Registrars>
254 bool operator==(const HalfSpaceMirror<Registrars>& lhs,
255  const HalfSpaceMirror<Registrars>& rhs) noexcept {
256  return lhs.beam_width() == rhs.beam_width() and
257  lhs.constitutive_relation() == rhs.constitutive_relation() and
258  lhs.integration_intervals() == rhs.integration_intervals() and
259  lhs.absolute_tolerance() == rhs.absolute_tolerance() and
260  lhs.relative_tolerance() == rhs.relative_tolerance();
261 }
262 
263 template <typename Registrars>
264 bool operator!=(const HalfSpaceMirror<Registrars>& lhs,
265  const HalfSpaceMirror<Registrars>& rhs) noexcept {
266  return not(lhs == rhs);
267 }
268 
269 } // 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:128
Elasticity::Solutions::HalfSpaceMirror::AbsoluteTolerance
Definition: HalfSpaceMirror.hpp:152
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:160
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:123
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:380
Scalar
Tensor< T, Symmetry<>, index_list<> > Scalar
Definition: TypeAliases.hpp:21
Elasticity::Solutions::HalfSpaceMirror::Material
Definition: HalfSpaceMirror.hpp:135
Elasticity::ConstitutiveRelations::IsotropicHomogeneous< 3 >
limits
Elasticity::Solutions::HalfSpaceMirror::IntegrationIntervals
Definition: HalfSpaceMirror.hpp:141
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: ReadSpecThirdOrderPiecewisePolynomial.hpp:13