Add white noise in pressure in the disk and normalize the density

src/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedFmDisk.cpp

@@ -28,10 +28,10 @@ MagnetizedFmDisk::MagnetizedFmDisk(
     const double bh_mass, const double bh_dimless_spin,
     const double inner_edge_radius, const double max_pressure_radius,
     const double polytropic_constant, const double polytropic_exponent,
-    const double threshold_density, const double inverse_plasma_beta,
-    const size_t normalization_grid_res)
+    const double noise, const double threshold_density,
+    const double inverse_plasma_beta, const size_t normalization_grid_res)
     : fm_disk_(bh_mass, bh_dimless_spin, inner_edge_radius, max_pressure_radius,
-               polytropic_constant, polytropic_exponent),
+               polytropic_constant, polytropic_exponent, noise),
       threshold_density_(threshold_density),
       inverse_plasma_beta_(inverse_plasma_beta),
       normalization_grid_res_(normalization_grid_res),

src/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedFmDisk.hpp

@@ -150,7 +150,7 @@ class MagnetizedFmDisk : public virtual evolution::initial_data::InitialData,
   MagnetizedFmDisk(
       double bh_mass, double bh_dimless_spin, double inner_edge_radius,
       double max_pressure_radius, double polytropic_constant,
-      double polytropic_exponent, double threshold_density,
+      double polytropic_exponent, double noise, double threshold_density,
       double inverse_plasma_beta,
       size_t normalization_grid_res = BFieldNormGridRes::suggested_value());
 

src/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/FishboneMoncriefDisk.cpp

@@ -7,6 +7,7 @@
 #include <cmath>
 #include <cstddef>
 #include <pup.h>
+#include <random>
 
 #include "DataStructures/DataVector.hpp"
 #include "DataStructures/Tensor/EagerMath/DotProduct.hpp"
@@ -31,13 +32,15 @@ FishboneMoncriefDisk::FishboneMoncriefDisk(const double bh_mass,
                                            const double inner_edge_radius,
                                            const double max_pressure_radius,
                                            const double polytropic_constant,
-                                           const double polytropic_exponent)
+                                           const double polytropic_exponent,
+                                           const double noise)
     : bh_mass_(bh_mass),
       bh_spin_a_(bh_mass * bh_dimless_spin),
       inner_edge_radius_(bh_mass * inner_edge_radius),
       max_pressure_radius_(bh_mass * max_pressure_radius),
       polytropic_constant_(polytropic_constant),
       polytropic_exponent_(polytropic_exponent),
+      noise_(noise),
       equation_of_state_{polytropic_constant_, polytropic_exponent_},
       background_spacetime_{
           bh_mass_, {{0.0, 0.0, bh_dimless_spin}}, {{0.0, 0.0, 0.0}}} {
@@ -52,6 +55,13 @@ FishboneMoncriefDisk::FishboneMoncriefDisk(const double bh_mass,
       (square(bh_spin_a_) - 2.0 * a_sqrt_m * sqrt_rmax + rmax_squared) /
       (2.0 * a_sqrt_m * rmax_sqrt_rmax +
        (rmax - 3.0 * bh_mass_) * rmax_squared);
+
+  // compute rho_max_ here:
+  const double potential_at_rmax = potential(rmax_squared, 1.0);
+  const double potential_at_rin = potential(square(inner_edge_radius_), 1.0);
+  const double h_at_rmax = exp(potential_at_rin - potential_at_rmax);
+  rho_max_ = get(equation_of_state_.rest_mass_density_from_enthalpy(
+      Scalar<double>{h_at_rmax}));
 }
 
 std::unique_ptr<evolution::initial_data::InitialData>
@@ -68,6 +78,8 @@ void FishboneMoncriefDisk::pup(PUP::er& p) {
   p | polytropic_constant_;
   p | polytropic_exponent_;
   p | angular_momentum_;
+  p | rho_max_;
+  p | noise_;
   p | equation_of_state_;
   p | background_spacetime_;
 }
@@ -140,6 +152,7 @@ FishboneMoncriefDisk::variables(
   variables_impl(vars, [&rest_mass_density, &specific_enthalpy, this](
                            const size_t s, const double /*potential_at_s*/) {
     get_element(get(rest_mass_density), s) =
+        (1 / rho_max_) *
         get(equation_of_state_.rest_mass_density_from_enthalpy(
             Scalar<double>{get_element(get(specific_enthalpy), s)}));
   });
@@ -183,14 +196,18 @@ FishboneMoncriefDisk::variables(
     const tnsr::I<DataType, 3>& x,
     tmpl::list<hydro::Tags::Pressure<DataType>> /*meta*/,
     gsl::not_null<IntermediateVariables<DataType>*> vars) const {
+  std::random_device rd;
+  std::mt19937 gen(rd());
+  std::uniform_real_distribution<> dis(-noise_, noise_);
   const auto rest_mass_density = get<hydro::Tags::RestMassDensity<DataType>>(
       variables(x, tmpl::list<hydro::Tags::RestMassDensity<DataType>>{}, vars));
   auto pressure = make_with_value<Scalar<DataType>>(x, 0.0);
-  variables_impl(vars, [&pressure, &rest_mass_density, this](
+  variables_impl(vars, [&pressure, &rest_mass_density, &gen, &dis, this](
                            const size_t s, const double /*potential_at_s*/) {
     get_element(get(pressure), s) =
+        (1 / rho_max_) * (1 + dis(gen)) *
         get(equation_of_state_.pressure_from_density(
-            Scalar<double>{get_element(get(rest_mass_density), s)}));
+            Scalar<double>{rho_max_ * get_element(get(rest_mass_density), s)}));
   });
   return {std::move(pressure)};
 }
@@ -208,7 +225,7 @@ FishboneMoncriefDisk::variables(
                            const size_t s, const double /*potential_at_s*/) {
     get_element(get(specific_internal_energy), s) =
         get(equation_of_state_.specific_internal_energy_from_density(
-            Scalar<double>{get_element(get(rest_mass_density), s)}));
+            Scalar<double>{rho_max_ * get_element(get(rest_mass_density), s)}));
   });
   return {std::move(specific_internal_energy)};
 }
@@ -358,7 +375,8 @@ bool operator==(const FishboneMoncriefDisk& lhs,
          lhs.inner_edge_radius_ == rhs.inner_edge_radius_ and
          lhs.max_pressure_radius_ == rhs.max_pressure_radius_ and
          lhs.polytropic_constant_ == rhs.polytropic_constant_ and
-         lhs.polytropic_exponent_ == rhs.polytropic_exponent_;
+         lhs.polytropic_exponent_ == rhs.polytropic_exponent_ and
+         lhs.noise_ == rhs.noise_ and lhs.rho_max_ == rhs.rho_max_;
 }
 
 bool operator!=(const FishboneMoncriefDisk& lhs,

src/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/FishboneMoncriefDisk.hpp

@@ -126,6 +126,9 @@ namespace Solutions {
  * p = K\rho^\gamma.
  * \f}
  *
+ * Following \cite Porth2016rfi, the rest mass density is normalized
+ * to be 1. at maximum.
+ *
  * Once all the variables are known in Boyer-Lindquist (or Kerr) coordinates, it
  * is straightforward to write them in Cartesian Kerr-Schild coordinates. The
  * coordinate transformation in gr::KerrSchildCoords helps read the Jacobian
@@ -213,10 +216,20 @@ class FishboneMoncriefDisk
         "The polytropic exponent of the fluid."};
     static type lower_bound() { return 1.; }
   };
+  /// The magnitude of noise added to pressure/energy of the Disk
+  /// to drive MRI.
+  struct Noise {
+    using type = double;
+    static constexpr Options::String help = {
+        "The magnitude of the white noise perturbation added to "
+        "pressure to excite MRI in the disk."};
+    static type lower_bound() { return 0.; }
+    static type upper_bound() { return 1.; }
+  };
 
   using options =
       tmpl::list<BhMass, BhDimlessSpin, InnerEdgeRadius, MaxPressureRadius,
-                 PolytropicConstant, PolytropicExponent>;
+                 PolytropicConstant, PolytropicExponent, Noise>;
   static constexpr Options::String help = {
       "Fluid disk orbiting a Kerr black hole."};
 
@@ -230,7 +243,8 @@ class FishboneMoncriefDisk
 
   FishboneMoncriefDisk(double bh_mass, double bh_dimless_spin,
                        double inner_edge_radius, double max_pressure_radius,
-                       double polytropic_constant, double polytropic_exponent);
+                       double polytropic_constant, double polytropic_exponent,
+                       double noise);
 
   auto get_clone() const
       -> std::unique_ptr<evolution::initial_data::InitialData> override;
@@ -410,6 +424,8 @@ class FishboneMoncriefDisk
   double polytropic_constant_ = std::numeric_limits<double>::signaling_NaN();
   double polytropic_exponent_ = std::numeric_limits<double>::signaling_NaN();
   double angular_momentum_ = std::numeric_limits<double>::signaling_NaN();
+  double rho_max_ = std::numeric_limits<double>::signaling_NaN();
+  double noise_ = std::numeric_limits<double>::signaling_NaN();
   EquationsOfState::PolytropicFluid<true> equation_of_state_{};
   gr::Solutions::SphericalKerrSchild background_spacetime_{};
 };

tests/InputFiles/GrMhd/ValenciaDivClean/FishboneMoncriefDisk.yaml

@@ -55,6 +55,7 @@ InitialData: &initial_data
     MaxPressureRadius: 12.0
     PolytropicConstant: 0.001
     PolytropicExponent: 1.3333333333333333333333
+    Noise: 0.0
 
 EquationOfState:
   FromInitialData

tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedFmDisk.py

@@ -14,6 +14,7 @@ def rest_mass_density(
     dimless_r_max,
     polytropic_constant,
     polytropic_exponent,
+    noise,
     threshold_density,
     plasma_beta,
 ):
@@ -27,6 +28,7 @@ def rest_mass_density(
         dimless_r_max,
         polytropic_constant,
         polytropic_exponent,
+        noise,
     )
 
 
@@ -38,6 +40,7 @@ def spatial_velocity(
     dimless_r_max,
     polytropic_constant,
     polytropic_exponent,
+    noise,
     threshold_density,
     plasma_beta,
 ):
@@ -51,6 +54,7 @@ def spatial_velocity(
         dimless_r_max,
         polytropic_constant,
         polytropic_exponent,
+        noise,
     )
 
 
@@ -62,6 +66,7 @@ def specific_internal_energy(
     dimless_r_max,
     polytropic_constant,
     polytropic_exponent,
+    noise,
     threshold_density,
     plasma_beta,
 ):
@@ -75,6 +80,7 @@ def specific_internal_energy(
         dimless_r_max,
         polytropic_constant,
         polytropic_exponent,
+        noise,
     )
 
 
@@ -86,6 +92,7 @@ def pressure(
     dimless_r_max,
     polytropic_constant,
     polytropic_exponent,
+    noise,
     threshold_density,
     plasma_beta,
 ):
@@ -99,6 +106,7 @@ def pressure(
         dimless_r_max,
         polytropic_constant,
         polytropic_exponent,
+        noise,
     )
 
 
@@ -111,6 +119,7 @@ def magnetic_potential(
     rmax,
     polytropic_constant,
     polytropic_exponent,
+    noise,
     threshold_rest_mass_density,
 ):
     l = fm_disk.angular_momentum(m, a, rmax)
@@ -135,6 +144,7 @@ def magnetic_field(
     dimless_r_max,
     polytropic_constant,
     polytropic_exponent,
+    noise,
     threshold_density,
     plasma_beta,
 ):
@@ -152,6 +162,7 @@ def magnetic_field(
         dimless_r_max,
         polytropic_constant,
         polytropic_exponent,
+        noise,
     )
     x_max = np.array([r_max, 0.0, 0.0])
     threshold_rho = threshold_density * (
@@ -164,6 +175,7 @@ def magnetic_field(
             dimless_r_max,
             polytropic_constant,
             polytropic_exponent,
+            noise,
         )
     )
     x_ks = [0.0, 0.0, 0.0]
@@ -191,6 +203,7 @@ def magnetic_field(
                     r_max,
                     polytropic_constant,
                     polytropic_exponent,
+                    noise,
                     threshold_rho,
                 )
                 - magnetic_potential(
@@ -202,6 +215,7 @@ def magnetic_field(
                     r_max,
                     polytropic_constant,
                     polytropic_exponent,
+                    noise,
                     threshold_rho,
                 )
             )
@@ -220,6 +234,7 @@ def magnetic_field(
                 r_max,
                 polytropic_constant,
                 polytropic_exponent,
+                noise,
                 threshold_rho,
             )
             - magnetic_potential(
@@ -231,6 +246,7 @@ def magnetic_field(
                 r_max,
                 polytropic_constant,
                 polytropic_exponent,
+                noise,
                 threshold_rho,
             )
         ) * prefactor
@@ -242,7 +258,7 @@ def magnetic_field(
         # in test_variables function of Test_MagnetizedFmDisk.cpp.
         # Note: we are using lower (than default)
         # b_field_normalization for a faster testing time.
-        normalization = 2.077412435947072
+        normalization = 8.49943510498341
         result = normalization * ks_coords.cartesian_from_spherical_ks(
             result, x_ks, bh_mass, bh_dimless_spin
         )
@@ -259,6 +275,7 @@ def divergence_cleaning_field(
     dimless_r_max,
     polytropic_constant,
     polytropic_exponent,
+    noise,
     threshold_density,
     plasma_beta,
 ):
@@ -273,6 +290,7 @@ def lorentz_factor(
     dimless_r_max,
     polytropic_constant,
     polytropic_exponent,
+    noise,
     threshold_density,
     plasma_beta,
 ):
@@ -286,6 +304,7 @@ def lorentz_factor(
         dimless_r_max,
         polytropic_constant,
         polytropic_exponent,
+        noise,
     )
 
 
@@ -297,6 +316,7 @@ def specific_enthalpy(
     dimless_r_max,
     polytropic_constant,
     polytropic_exponent,
+    noise,
     threshold_density,
     plasma_beta,
 ):
@@ -310,4 +330,5 @@ def specific_enthalpy(
         dimless_r_max,
         polytropic_constant,
         polytropic_exponent,
+        noise,
     )