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