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