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