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