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/Tensor/TypeAliases.hpp"
      10             : #include "Domain/Tags.hpp"
      11             : #include "Evolution/Systems/NewtonianEuler/Sources/Source.hpp"
      12             : #include "Evolution/Systems/NewtonianEuler/TagsDeclarations.hpp"
      13             : #include "Options/String.hpp"
      14             : #include "PointwiseFunctions/AnalyticSolutions/Tags.hpp"
      15             : #include "Utilities/MakeArray.hpp"
      16             : #include "Utilities/TMPL.hpp"
      17             : 
      18             : // IWYU pragma: no_forward_declare Tensor
      19             : 
      20             : /// \cond
      21             : class DataVector;
      22             : 
      23             : namespace NewtonianEuler {
      24             : namespace Solutions {
      25             : template <size_t Dim>
      26             : struct IsentropicVortex;
      27             : }  // namespace Solutions
      28             : }  // namespace NewtonianEuler
      29             : 
      30             : namespace Tags {
      31             : struct Time;
      32             : }  // namespace Tags
      33             : 
      34             : namespace gsl {
      35             : template <typename T>
      36             : class not_null;
      37             : }  // namespace gsl
      38             : /// \endcond
      39             : 
      40             : namespace NewtonianEuler::Sources {
      41             : 
      42             : /*!
      43             :  * \brief Source generating a modified isentropic vortex.
      44             :  *
      45             :  * If Solutions::IsentropicVortex is modifed so that the flow velocity along
      46             :  * the \f$z-\f$axis is not a constant but a function of \f$z\f$, the new vortex
      47             :  * will be a solution to the 3-D Newtonian Euler equations with a source term,
      48             :  *
      49             :  * \f{align*}
      50             :  * \partial_t\rho + \partial_i F^i(\rho) &= S(\rho)\\
      51             :  * \partial_t S^i + \partial_j F^{j}(S^i) &= S(S^i)\\
      52             :  * \partial_t e + \partial_i F^i(e) &= S(e),
      53             :  * \f}
      54             :  *
      55             :  * where \f$F^i(u)\f$ is the volume flux of the conserved quantity \f$u\f$
      56             :  * (see ComputeFluxes), and
      57             :  *
      58             :  * \f{align*}
      59             :  * S(\rho) &= \rho \dfrac{dv_z}{dz}\\
      60             :  * S(S_x) &= S_x \dfrac{dv_z}{dz}\\
      61             :  * S(S_y) &= S_y \dfrac{dv_z}{dz}\\
      62             :  * S(S_z) &= 2S_z \dfrac{dv_z}{dz}\\
      63             :  * S(e) &= \left(e + p + v_z S_z\right)\dfrac{dv_z}{dz},
      64             :  * \f}
      65             :  *
      66             :  * where \f$\rho\f$ is the mass density of the vortex, \f$S_i\f$ is
      67             :  * its momentum density, \f$e\f$ is its energy density,
      68             :  * \f$v_z = v_z(z)\f$ is the \f$z-\f$component of its velocity,
      69             :  * and \f$p\f$ is its pressure. These quantities are readily obtained
      70             :  * from the primitive variables, whose expressions are those in
      71             :  * Solutions::IsentropicVortex
      72             :  *
      73             :  * Currently \f$v_z(z)=sin(z)\f$ is hard-coded.
      74             :  */
      75           1 : class VortexPerturbation : public Source<3> {
      76             :  public:
      77             :   /// The perturbation amplitude
      78           1 :   struct PerturbationAmplitude {
      79           0 :     using type = double;
      80           0 :     static constexpr Options::String help = {"The perturbation amplitude."};
      81             :   };
      82             : 
      83           0 :   using options = tmpl::list<PerturbationAmplitude>;
      84             : 
      85           0 :   static constexpr Options::String help = {
      86             :       "Source terms corresponding to a vortex perturbation. Should be used "
      87             :       "with the IsentropicVortex solution."};
      88             : 
      89           0 :   explicit VortexPerturbation(double perturbation_amplitude);
      90             : 
      91           0 :   VortexPerturbation() = default;
      92           0 :   VortexPerturbation(const VortexPerturbation& /*rhs*/) = default;
      93           0 :   VortexPerturbation& operator=(const VortexPerturbation& /*rhs*/) = default;
      94           0 :   VortexPerturbation(VortexPerturbation&& /*rhs*/) = default;
      95           0 :   VortexPerturbation& operator=(VortexPerturbation&& /*rhs*/) = default;
      96           0 :   ~VortexPerturbation() override = default;
      97             : 
      98             :   /// \cond
      99             :   explicit VortexPerturbation(CkMigrateMessage* msg);
     100             :   using PUP::able::register_constructor;
     101             :   WRAPPED_PUPable_decl_template(VortexPerturbation);
     102             :   /// \endcond
     103             : 
     104             :   // NOLINTNEXTLINE(google-runtime-references)
     105           0 :   void pup(PUP::er& p) override;
     106             : 
     107           0 :   auto get_clone() const -> std::unique_ptr<Source<3>> override;
     108             : 
     109           0 :   void operator()(
     110             :       gsl::not_null<Scalar<DataVector>*> source_mass_density_cons,
     111             :       gsl::not_null<tnsr::I<DataVector, 3>*> source_momentum_density,
     112             :       gsl::not_null<Scalar<DataVector>*> source_energy_density,
     113             :       const Scalar<DataVector>& mass_density_cons,
     114             :       const tnsr::I<DataVector, 3>& momentum_density,
     115             :       const Scalar<DataVector>& energy_density,
     116             :       const tnsr::I<DataVector, 3>& velocity,
     117             :       const Scalar<DataVector>& pressure,
     118             :       const Scalar<DataVector>& specific_internal_energy,
     119             :       const EquationsOfState::EquationOfState<false, 2>& eos,
     120             :       const tnsr::I<DataVector, 3>& coords, double time) const override;
     121             : 
     122             :  private:
     123           0 :   double perturbation_amplitude_ = std::numeric_limits<double>::signaling_NaN();
     124             : };
     125             : }  // namespace NewtonianEuler::Sources