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