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