Add VelocityPerturbation option for TovStar

src/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedTovStar.cpp

@@ -32,15 +32,15 @@ MagnetizedTovStar& MagnetizedTovStar::operator=(MagnetizedTovStar&& /*rhs*/) =
 MagnetizedTovStar::~MagnetizedTovStar() = default;
 
 MagnetizedTovStar::MagnetizedTovStar(
-    const double central_rest_mass_density,
+    const double central_rest_mass_density, const double velocity_perturbation,
     std::unique_ptr<MagnetizedTovStar::equation_of_state_type>
         equation_of_state,
     const RelativisticEuler::Solutions::TovCoordinates coordinate_system,
     std::vector<std::unique_ptr<
         grmhd::AnalyticData::InitialMagneticFields::InitialMagneticField>>
         magnetic_fields)
-    : tov_star(central_rest_mass_density, std::move(equation_of_state),
-               coordinate_system),
+    : tov_star(central_rest_mass_density, velocity_perturbation,
+               std::move(equation_of_state), coordinate_system),
       magnetic_fields_(std::move(magnetic_fields)) {}
 
 MagnetizedTovStar::MagnetizedTovStar(const MagnetizedTovStar& rhs)

src/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedTovStar.hpp

@@ -41,6 +41,7 @@ struct MagnetizedTovVariables
   using Base::operator();
   using Base::coords;
   using Base::eos;
+  using Base::perturbation;
   using Base::radial_solution;
   using Base::radius;
 
@@ -52,10 +53,11 @@ struct MagnetizedTovVariables
       const tnsr::I<DataType, 3>& local_x, const DataType& local_radius,
       const RelativisticEuler::Solutions::TovSolution& local_radial_solution,
       const EquationsOfState::EquationOfState<true, 1>& local_eos,
+      const double perturb,
       const std::vector<std::unique_ptr<
           grmhd::AnalyticData::InitialMagneticFields::InitialMagneticField>>&
           mag_fields)
-      : Base(local_x, local_radius, local_radial_solution, local_eos),
+      : Base(local_x, local_radius, local_radial_solution, local_eos, perturb),
         magnetic_fields(mag_fields) {}
 
   void operator()(
@@ -145,7 +147,7 @@ class MagnetizedTovStar : public virtual evolution::initial_data::InitialData,
   ~MagnetizedTovStar() override;
 
   MagnetizedTovStar(
-      double central_rest_mass_density,
+      double central_rest_mass_density, double velocity_perturbation,
       std::unique_ptr<EquationsOfState::EquationOfState<true, 1>>
           equation_of_state,
       RelativisticEuler::Solutions::TovCoordinates coordinate_system,

src/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/TovStar.cpp

@@ -20,11 +20,12 @@ namespace RelativisticEuler::Solutions {
 TovStar::TovStar(CkMigrateMessage* msg) : InitialData(msg) {}
 
 TovStar::TovStar(
-    const double central_rest_mass_density,
+    const double central_rest_mass_density, const double velocity_perturbation,
     std::unique_ptr<EquationsOfState::EquationOfState<true, 1>>
         equation_of_state,
     const RelativisticEuler::Solutions::TovCoordinates coordinate_system)
     : central_rest_mass_density_(central_rest_mass_density),
+      velocity_perturbation_(velocity_perturbation),
       equation_of_state_(std::move(equation_of_state)),
       coordinate_system_(coordinate_system),
       radial_solution_(*equation_of_state_, central_rest_mass_density_,
@@ -33,13 +34,15 @@ TovStar::TovStar(
 TovStar::TovStar(const TovStar& rhs)
     : evolution::initial_data::InitialData(rhs),
       central_rest_mass_density_(rhs.central_rest_mass_density_),
+      velocity_perturbation_(rhs.velocity_perturbation_),
       equation_of_state_(rhs.equation_of_state_->get_clone()),
       coordinate_system_(rhs.coordinate_system_),
       radial_solution_(*equation_of_state_, central_rest_mass_density_,
                        coordinate_system_) {}
 
 TovStar& TovStar::operator=(const TovStar& rhs) {
   central_rest_mass_density_ = rhs.central_rest_mass_density_;
+  velocity_perturbation_ = rhs.velocity_perturbation_;
   equation_of_state_ = rhs.equation_of_state_->get_clone();
   coordinate_system_ = rhs.coordinate_system_;
   radial_solution_ = RelativisticEuler::Solutions::TovSolution(
@@ -55,6 +58,7 @@ std::unique_ptr<evolution::initial_data::InitialData> TovStar::get_clone()
 void TovStar::pup(PUP::er& p) {
   InitialData::pup(p);
   p | central_rest_mass_density_;
+  p | velocity_perturbation_;
   p | equation_of_state_;
   p | coordinate_system_;
   p | radial_solution_;
@@ -514,11 +518,23 @@ void TovVariables<DataType, Region>::operator()(
 template <typename DataType, StarRegion Region>
 void TovVariables<DataType, Region>::operator()(
     const gsl::not_null<tnsr::I<DataType, 3>*> spatial_velocity,
-    const gsl::not_null<Cache*> /*cache*/,
+    [[maybe_unused]] const gsl::not_null<Cache*> cache,
     hydro::Tags::SpatialVelocity<DataType, 3> /*meta*/) const {
-  get<0>(*spatial_velocity) = 0.;
-  get<1>(*spatial_velocity) = 0.;
-  get<2>(*spatial_velocity) = 0.;
+  if constexpr (Region == StarRegion::Center or
+                Region == StarRegion::Exterior) {
+    get<0>(*spatial_velocity) = 0.;
+    get<1>(*spatial_velocity) = 0.;
+    get<2>(*spatial_velocity) = 0.;
+  } else {
+    const auto& areal_radius =
+        radial_solution.coordinate_system() ==
+                RelativisticEuler::Solutions::TovCoordinates::Isotropic
+            ? get(cache->get_var(*this, Tags::ArealRadius<DataType>{}))
+            : radius;
+    get<0>(*spatial_velocity) = perturbation * coords.get(0) / areal_radius;
+    get<1>(*spatial_velocity) = perturbation * coords.get(1) / areal_radius;
+    get<2>(*spatial_velocity) = perturbation * coords.get(2) / areal_radius;
+  }
 }
 
 template <typename DataType, StarRegion Region>

src/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/TovStar.hpp

@@ -139,16 +139,19 @@ struct TovVariables {
   TovVariables(
       const tnsr::I<DataType, 3>& local_coords, const DataType& local_radius,
       const RelativisticEuler::Solutions::TovSolution& local_radial_solution,
-      const EquationsOfState::EquationOfState<true, 1>& local_eos)
+      const EquationsOfState::EquationOfState<true, 1>& local_eos,
+      const double perturb)
       : coords(local_coords),
         radius(local_radius),
         radial_solution(local_radial_solution),
-        eos(local_eos) {}
+        eos(local_eos),
+        perturbation(perturb) {}
 
   const tnsr::I<DataType, 3>& coords;
   const DataType& radius;
   const RelativisticEuler::Solutions::TovSolution& radial_solution;
   const EquationsOfState::EquationOfState<true, 1>& eos;
+  const double perturbation;
 
   void operator()(gsl::not_null<Scalar<DataType>*> mass_over_radius,
                   gsl::not_null<Cache*> cache,
@@ -336,6 +339,13 @@ class TovStar : public virtual evolution::initial_data::InitialData,
     static type lower_bound() { return 0.; }
   };
 
+  struct VelocityPerturbation {
+    using type = double;
+    static constexpr Options::String help = {
+        "Magnitude of radially symmetric velocity perturbation for driving "
+        "oscillations."};
+  };
+
   /// Areal (Schwarzschild) or isotropic coordinates
   struct Coordinates {
     using type = RelativisticEuler::Solutions::TovCoordinates;
@@ -346,7 +356,7 @@ class TovStar : public virtual evolution::initial_data::InitialData,
   static constexpr size_t volume_dim = 3_st;
 
   using options =
-      tmpl::list<CentralDensity,
+      tmpl::list<CentralDensity, VelocityPerturbation,
                  hydro::OptionTags::InitialDataEquationOfState<true, 1>,
                  Coordinates>;
 
@@ -361,7 +371,7 @@ class TovStar : public virtual evolution::initial_data::InitialData,
   TovStar(TovStar&& /*rhs*/) = default;
   TovStar& operator=(TovStar&& /*rhs*/) = default;
   ~TovStar() override = default;
-  TovStar(double central_rest_mass_density,
+  TovStar(double central_rest_mass_density, double velocity_perturbation,
           std::unique_ptr<EquationsOfState::EquationOfState<true, 1>>
               equation_of_state,
           const RelativisticEuler::Solutions::TovCoordinates coordinate_system =
@@ -420,7 +430,11 @@ class TovStar : public virtual evolution::initial_data::InitialData,
           VarsComputer<DataType, tov_detail::StarRegion::Exterior>;
       typename ExteriorVarsComputer::Cache cache{get_size(radius)};
       ExteriorVarsComputer computer{
-          x, radius, radial_solution_, *equation_of_state_,
+          x,
+          radius,
+          radial_solution_,
+          *equation_of_state_,
+          velocity_perturbation_,
           std::forward<VarsComputerArgs>(vars_computer_args)...};
       return {cache.get_var(computer, Tags{})...};
     } else if (max(radius) <= outer_radius and
@@ -430,7 +444,11 @@ class TovStar : public virtual evolution::initial_data::InitialData,
           VarsComputer<DataType, tov_detail::StarRegion::Interior>;
       typename InteriorVarsComputer::Cache cache{get_size(radius)};
       InteriorVarsComputer computer{
-          x, radius, radial_solution_, *equation_of_state_,
+          x,
+          radius,
+          radial_solution_,
+          *equation_of_state_,
+          velocity_perturbation_,
           std::forward<VarsComputerArgs>(vars_computer_args)...};
       return {cache.get_var(computer, Tags{})...};
     } else {
@@ -472,22 +490,31 @@ class TovStar : public virtual evolution::initial_data::InitialData,
         if (get_element(radius, i) > outer_radius) {
           typename ExteriorVarsComputer::Cache cache{1};
           ExteriorVarsComputer computer{
-              x_i, get_element(radius, i), radial_solution_,
+              x_i,
+              get_element(radius, i),
+              radial_solution_,
               *equation_of_state_,
+              velocity_perturbation_,
               std::forward<VarsComputerArgs>(vars_computer_args)...};
           expand_pack(get_var(i, cache, computer, Tags{})...);
         } else if (get_element(radius, i) > center_radius_cutoff) {
           typename InteriorVarsComputer::Cache cache{1};
           InteriorVarsComputer computer{
-              x_i, get_element(radius, i), radial_solution_,
+              x_i,
+              get_element(radius, i),
+              radial_solution_,
               *equation_of_state_,
+              velocity_perturbation_,
               std::forward<VarsComputerArgs>(vars_computer_args)...};
           expand_pack(get_var(i, cache, computer, Tags{})...);
         } else {
           typename CenterVarsComputer::Cache cache{1};
           CenterVarsComputer computer{
-              x_i, get_element(radius, i), radial_solution_,
+              x_i,
+              get_element(radius, i),
+              radial_solution_,
               *equation_of_state_,
+              velocity_perturbation_,
               std::forward<VarsComputerArgs>(vars_computer_args)...};
           expand_pack(get_var(i, cache, computer, Tags{})...);
         }
@@ -521,6 +548,7 @@ class TovStar : public virtual evolution::initial_data::InitialData,
 
   double central_rest_mass_density_ =
       std::numeric_limits<double>::signaling_NaN();
+  double velocity_perturbation_ = std::numeric_limits<double>::signaling_NaN();
   std::unique_ptr<equation_of_state_type> equation_of_state_;
   RelativisticEuler::Solutions::TovCoordinates coordinate_system_{};
   RelativisticEuler::Solutions::TovSolution radial_solution_{};

src/PointwiseFunctions/AnalyticSolutions/Xcts/TovStar.hpp

@@ -195,12 +195,12 @@ class TovStar : public elliptic::analytic_data::AnalyticSolution {
   TovStar& operator=(TovStar&&) = default;
   ~TovStar() = default;
 
-  TovStar(double central_rest_mass_density,
+  TovStar(double central_rest_mass_density, double velocity_perturbation,
           std::unique_ptr<EquationsOfState::EquationOfState<true, 1>>
               equation_of_state,
           const RelativisticEuler::Solutions::TovCoordinates coordinate_system)
-      : tov_star(central_rest_mass_density, std::move(equation_of_state),
-                 coordinate_system) {}
+      : tov_star(central_rest_mass_density, velocity_perturbation,
+                 std::move(equation_of_state), coordinate_system) {}
 
   const EquationsOfState::EquationOfState<true, 1>& equation_of_state() const {
     return tov_star.equation_of_state();

tests/Unit/Elliptic/Systems/Xcts/Test_HydroQuantities.cpp

@@ -17,8 +17,8 @@
 #include "PointwiseFunctions/AnalyticSolutions/Xcts/Factory.hpp"
 #include "PointwiseFunctions/AnalyticSolutions/Xcts/TovStar.hpp"
 #include "PointwiseFunctions/Hydro/EquationsOfState/PolytropicFluid.hpp"
-#include "PointwiseFunctions/Hydro/Tags.hpp"
 #include "PointwiseFunctions/Hydro/LowerSpatialFourVelocity.hpp"
+#include "PointwiseFunctions/Hydro/Tags.hpp"
 #include "PointwiseFunctions/InitialDataUtilities/AnalyticSolution.hpp"
 #include "PointwiseFunctions/InitialDataUtilities/Background.hpp"
 #include "Utilities/ProtocolHelpers.hpp"
@@ -47,7 +47,8 @@ SPECTRE_TEST_CASE("Unit.Elliptic.Systems.Xcts.HydroQuantities",
   using hydro_tags = AnalyticData::hydro_tags<DataVector>;
   using spatial_metric_tag = gr::Tags::SpatialMetric<DataVector, 3>;
   const Solutions::TovStar tov_star{
-      1.e-3, std::make_unique<EquationsOfState::PolytropicFluid<true>>(1., 2.),
+      1.e-3, 0.0,
+      std::make_unique<EquationsOfState::PolytropicFluid<true>>(1., 2.),
       RelativisticEuler::Solutions::TovCoordinates::Schwarzschild};
   const double outer_radius = tov_star.radial_solution().outer_radius();
   const tnsr::I<DataVector, 3> x{

tests/Unit/Evolution/DgSubcell/Test_BackgroundGrVars.cpp

@@ -226,7 +226,7 @@ void test(const gsl::not_null<std::mt19937*> gen, const bool did_rollback) {
   const auto solution = []() {
     if constexpr (TestRuntimeInitialData) {
       return RelativisticEuler::Solutions::TovStar{
-          1.0e-3,
+          1.0e-3, 0.0,
           EquationsOfState::PolytropicFluid<true>{100.0, 2.0}.get_clone(),
           RelativisticEuler::Solutions::TovCoordinates::Schwarzschild};
     } else {

tests/Unit/Evolution/DiscontinuousGalerkin/Test_BackgroundGrVars.cpp

@@ -170,7 +170,7 @@ void test(const gsl::not_null<std::mt19937*> gen) {
   const auto solution = []() {
     if constexpr (TestRuntimeInitialData) {
       return RelativisticEuler::Solutions::TovStar{
-          1.0e-3,
+          1.0e-3, 0.0,
           EquationsOfState::PolytropicFluid<true>{100.0, 2.0}.get_clone(),
           RelativisticEuler::Solutions::TovCoordinates::Schwarzschild};
     } else {

tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/Actions/Test_SetInitialData.cpp

@@ -161,7 +161,7 @@ void test_set_initial_data(
 
   // We get test data from a TOV solution
   TovStar tov_star{
-      1.e-3,
+      1.e-3, 0.0,
       std::make_unique<EquationsOfState::PolytropicFluid<true>>(100., 2.)};
   const double star_radius = tov_star.radial_solution().outer_radius();
   const auto eos = tov_star.equation_of_state().promote_to_3d_eos();
@@ -381,11 +381,12 @@ SPECTRE_TEST_CASE("Unit.Evolution.Systems.GhValenciaDivClean.SetInitialData",
         "  DensityCutoff: 1.e-14",
         true, active_grid);
     test_set_initial_data(
-        TovStar{1.e-3,
+        TovStar{1.e-3, 0.0,
                 std::make_unique<EquationsOfState::PolytropicFluid<true>>(100.,
                                                                           2.)},
         "GeneralizedHarmonic(TovStar):\n"
         "  CentralDensity: 1.e-3\n"
+        "  VelocityPerturbation: 0.0\n"
         "  EquationOfState:\n"
         "    PolytropicFluid:\n"
         "      PolytropicConstant: 100.\n"

tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/BoundaryConditions/Test_DirichletAnalytic.cpp

@@ -351,6 +351,7 @@ SPECTRE_TEST_CASE(
             "  AnalyticPrescription:\n"
             "      GeneralizedHarmonic(MagnetizedTovStar):\n"
             "        CentralDensity: 1.28e-3\n"
+            "        VelocityPerturbation: 0.0\n"
             "        EquationOfState:\n"
             "          PolytropicFluid:\n"
             "            PolytropicConstant: 100.0\n"
@@ -367,7 +368,7 @@ SPECTRE_TEST_CASE(
 
     const gh::Solutions::WrappedGr<grmhd::AnalyticData::MagnetizedTovStar>
         analytic_solution_or_data{
-            1.28e-3,
+            1.28e-3, 0.0,
             std::make_unique<EquationsOfState::PolytropicFluid<true>>(100.0,
                                                                       2.0),
             RelativisticEuler::Solutions::TovCoordinates::Schwarzschild,

tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/BoundaryConditions/Test_DirichletFreeOutflow.cpp

@@ -551,6 +551,7 @@ SPECTRE_TEST_CASE(
             "  AnalyticPrescription:\n"
             "      GeneralizedHarmonic(MagnetizedTovStar):\n"
             "        CentralDensity: 1.28e-3\n"
+            "        VelocityPerturbation: 0.0\n"
             "        EquationOfState:\n"
             "          PolytropicFluid:\n"
             "            PolytropicConstant: 100.0\n"
@@ -566,7 +567,7 @@ SPECTRE_TEST_CASE(
             ->get_clone();
     const gh::Solutions::WrappedGr<grmhd::AnalyticData::MagnetizedTovStar>
         analytic_solution_or_data{
-            1.28e-3,
+            1.28e-3, 0.0,
             std::make_unique<EquationsOfState::PolytropicFluid<true>>(100.0,
                                                                       2.0),
             RelativisticEuler::Solutions::TovCoordinates::Schwarzschild,

tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Test_BoundaryConditionGhostData.cpp

@@ -339,7 +339,8 @@ SPECTRE_TEST_CASE(
   pypp::SetupLocalPythonEnvironment local_python_env{
       "PointwiseFunctions/AnalyticSolutions/"};
   const SolutionForTest solution{RelativisticEuler::Solutions::TovStar{
-      1.28e-3, EquationsOfState::PolytropicFluid<true>{100.0, 2.0}.get_clone(),
+      1.28e-3, 0.0,
+      EquationsOfState::PolytropicFluid<true>{100.0, 2.0}.get_clone(),
       RelativisticEuler::Solutions::TovCoordinates::Schwarzschild}};
   test(
       grmhd::GhValenciaDivClean::BoundaryConditions::DirichletAnalytic{

tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/Test_NeighborPackagedData.cpp

@@ -112,7 +112,7 @@ double test(const size_t num_dg_pts) {
           domain::CoordinateMaps::Identity<3>{});
 
   const gh::Solutions::WrappedGr<::RelativisticEuler::Solutions::TovStar> soln{
-      1.28e-3,
+      1.28e-3, 0.0,
       std::make_unique<EquationsOfState::PolytropicFluid<true>>(100.0, 2.0)
           ->get_clone(),
       RelativisticEuler::Solutions::TovCoordinates::Schwarzschild};

tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/Test_TimeDerivative.cpp

@@ -114,7 +114,7 @@ double test(const size_t num_dg_pts, std::optional<double> expansion_velocity,
       domain::CoordinateMaps::ProductOf3Maps<Affine, Affine, Affine>;
 
   const gh::Solutions::WrappedGr<::RelativisticEuler::Solutions::TovStar> soln{
-      1.28e-3,
+      1.28e-3, 0.0,
       std::make_unique<EquationsOfState::PolytropicFluid<true>>(100.0, 2.0)
           ->get_clone(),
       RelativisticEuler::Solutions::TovCoordinates::Schwarzschild};
@@ -267,8 +267,7 @@ double test(const size_t num_dg_pts, std::optional<double> expansion_velocity,
       neighbor_data{};
   using prims_to_reconstruct_tags = grmhd::GhValenciaDivClean::Tags::
       primitive_grmhd_and_spacetime_reconstruction_tags;
-  for (const auto & [ direction, neighbors_in_direction ] :
-       element.neighbors()) {
+  for (const auto& [direction, neighbors_in_direction] : element.neighbors()) {
     auto neighbor_logical_coords = logical_coordinates(subcell_mesh);
     neighbor_logical_coords.get(direction.dimension()) +=
         2.0 * direction.sign();

tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/Actions/Test_NumericInitialData.cpp

@@ -133,7 +133,7 @@ void test_numeric_initial_data(const NumericInitialData& initial_data,
 
   // We get test data from a TOV solution
   RelativisticEuler::Solutions::TovStar tov_star{
-      1.e-3,
+      1.e-3, 0.0,
       std::make_unique<EquationsOfState::PolytropicFluid<true>>(100., 2.)};
   const double star_radius = tov_star.radial_solution().outer_radius();
   const auto& eos = tov_star.equation_of_state();