Add different initial data options and params

Examples/PerfectFluid/GaussianID.hpp

@@ -0,0 +1,99 @@
+/* GRChombo
+ * Copyright 2012 The GRChombo collaboration.
+ * Please refer to LICENSE in GRChombo's root directory.
+ */
+
+#ifndef INITIALFLUIDDATA_HPP_
+#define INITIALFLUIDDATA_HPP_
+
+#include "ADMConformalVars.hpp"
+#include "Cell.hpp"
+#include "Coordinates.hpp"
+#include "FluidCCZ4RHS.hpp"
+#include "PerfectFluid.hpp"
+#include "Tensor.hpp"
+#include "UserVariables.hpp" //This files needs NUM_VARS - total no. components
+#include "VarsTools.hpp"
+#include "simd.hpp"
+
+//! Class which sets the initial fluid matter config
+// aka Kevin data
+class InitialFluidData
+{
+    // Now the non grid ADM vars
+    template <class data_t>
+    using MetricVars = ADMConformalVars::VarsWithGauge<data_t>;
+
+  public:
+    //! A structure for the input params for fluid properties and initial
+    //! conditions
+    struct params_t
+    {
+        std::array<double, CH_SPACEDIM>
+            center; //!< Centre of perturbation in initial SF bubble
+        double rho0;
+        double awidth;
+        double delta;
+    };
+
+    //! The constructor
+    InitialFluidData(params_t a_params, double a_dx)
+        : m_dx(a_dx), m_params(a_params)
+    {
+    }
+
+    //! Function to compute the value of all the initial vars on the grid
+    template <class data_t> void compute(Cell<data_t> current_cell) const
+    {
+        // where am i?
+        Coordinates<data_t> coords(current_cell, m_dx, m_params.center);
+        data_t rr = coords.get_radius();
+        data_t rr2 = rr * rr;
+
+        const auto metric_vars = current_cell.template load_vars<MetricVars>();
+
+        data_t x = coords.x;
+        double y = coords.y;
+        double z = coords.z;
+
+        Tensor<1, data_t> vi, Sj;
+        data_t chi_regularised = simd_max(metric_vars.chi, 1e-6);
+
+        vi[0] = 0.;
+        vi[1] = 0.;
+        vi[2] = 0.;
+
+        // calculate the field value
+        data_t rho = m_params.rho0 * (exp(-pow(rr / m_params.awidth, 2.0))) +
+                     m_params.delta;
+        data_t v2 = 0.;
+        FOR(i, j) v2 += metric_vars.h[i][j] * vi[i] * vi[j] / chi_regularised;
+        data_t eps = 0.;
+        data_t P = rho * (1. + eps) / 3.;
+        data_t WW = 1. / (1. - v2);
+        data_t hh = 1. + eps + P / rho;
+
+        data_t D = rho * sqrt(WW);
+        data_t tau = rho * hh * WW - P - D;
+        FOR(i)
+        {
+            Sj[i] = 0.;
+            FOR(j)
+            Sj[i] +=
+                rho * hh * WW * metric_vars.h[i][j] * vi[j] / chi_regularised;
+        }
+
+        // store the vars
+        current_cell.store_vars(rho, c_rho);
+        current_cell.store_vars(vi, GRInterval<c_vi1, c_vi3>());
+        current_cell.store_vars(D, c_D);
+        current_cell.store_vars(Sj, GRInterval<c_Sj1, c_Sj3>());
+        current_cell.store_vars(tau, c_tau);
+    }
+
+  protected:
+    double m_dx;
+    const params_t m_params;
+};
+
+#endif /* INITIALFLUIDDATA_HPP_ */

Examples/PerfectFluid/GaussianID_SimulationParameters.hpp

@@ -0,0 +1,82 @@
+/* GRChombo
+ * Copyright 2012 The GRChombo collaboration.
+ * Please refer to LICENSE in GRChombo's root directory.
+ */
+
+#ifndef SIMULATIONPARAMETERS_HPP_
+#define SIMULATIONPARAMETERS_HPP_
+
+// General includes
+#include "GRParmParse.hpp"
+#include "SimulationParametersBase.hpp"
+
+// Problem specific includes:
+#include "EoS.hpp"
+#include "InitialFluidData.hpp"
+#include "Minkowski.hpp"
+
+class SimulationParameters : public SimulationParametersBase
+{
+  public:
+    SimulationParameters(GRParmParse &pp) : SimulationParametersBase(pp)
+    {
+        // read the problem specific params
+        read_params(pp);
+        check_params();
+    }
+
+    void read_params(GRParmParse &pp)
+    {
+        // Initial scalar field data
+        initial_params.center =
+            center; // already read in SimulationParametersBase
+        pp.load("G_Newton", G_Newton,
+                0.0); // for now the example neglects backreaction
+        pp.load("fluid_rho0", initial_params.rho0, 1.0);
+        pp.load("fluid_width", initial_params.awidth, 0.1);
+        pp.load("fluid_delta", initial_params.delta, 0.2);
+        pp.load("lambda", lambda, 1.); // eigenvalue for numerical flux
+        pp.load("eos_w", eos_params.eos_w, 1. / 3.);
+
+        // Initial Kerr data
+        pp.load("kerr_mass", kerr_params.mass, 1.0);
+        pp.load("kerr_spin", kerr_params.spin, 0.0);
+        pp.load("kerr_center", kerr_params.center, center);
+    }
+
+    void check_params()
+    {
+        //        warn_parameter("scalar_mass", potential_params.scalar_mass,
+        //                potential_params.scalar_mass <
+        //                   0.2 / coarsest_dx / dt_multiplier,
+        //               "oscillations of scalar field do not appear to be "
+        //               "resolved on coarsest level");
+        // warn_parameter("scalar_width", initial_params.width,
+        //               initial_params.width < 0.5 * L,
+        //               "is greater than half the domain size");
+        warn_parameter("kerr_mass", kerr_params.mass, kerr_params.mass >= 0.0,
+                       "should be >= 0.0");
+        check_parameter("kerr_spin", kerr_params.spin,
+                        std::abs(kerr_params.spin) <= kerr_params.mass,
+                        "must satisfy |a| <= M = " +
+                            std::to_string(kerr_params.mass));
+        FOR(idir)
+        {
+            std::string name = "kerr_center[" + std::to_string(idir) + "]";
+            warn_parameter(
+                name, kerr_params.center[idir],
+                (kerr_params.center[idir] >= 0) &&
+                    (kerr_params.center[idir] <= (ivN[idir] + 1) * coarsest_dx),
+                "should be within the computational domain");
+        }
+    }
+
+    // Initial data for matter and potential and BH
+    double G_Newton;
+    double lambda;
+    InitialFluidData::params_t initial_params;
+    EoS::params_t eos_params;
+    Minkowski::params_t kerr_params;
+};
+
+#endif /* SIMULATIONPARAMETERS_HPP_ */

Examples/PerfectFluid/InitialFluidData.hpp

@@ -31,11 +31,12 @@ class InitialFluidData
     {
         std::array<double, CH_SPACEDIM>
             center; //!< Centre of perturbation in initial SF bubble
+        double L;
         double rho0;
         double uflow;
         double amplitude;
         double awidth;
-        double sigma;
+        double asigma;
         std::array<double, 2> ycenter;
     };
 
@@ -63,9 +64,9 @@ class InitialFluidData
         vi[0] = m_params.uflow *
                 (tanh((y - m_params.ycenter[0]) / m_params.awidth) -
                  tanh((y - m_params.ycenter[1]) / m_params.awidth) - 1.);
-        vi[1] = m_params.amplitude * sin(2. * M_PI * x) *
-                (exp(-pow((y - m_params.ycenter[0]) / m_params.sigma, 2)) +
-                 exp(-pow((y - m_params.ycenter[1]) / m_params.sigma, 2)));
+        vi[1] = m_params.amplitude * sin(2. * M_PI * x / m_params.L) *
+                (exp(-pow((y - m_params.ycenter[0]) / m_params.asigma, 2)) +
+                 exp(-pow((y - m_params.ycenter[1]) / m_params.asigma, 2)));
         vi[2] = 0.;
 
         data_t nn =

Examples/PerfectFluid/KelvinHelmhotzID.hpp

@@ -0,0 +1,107 @@
+/* GRChombo
+ * Copyright 2012 The GRChombo collaboration.
+ * Please refer to LICENSE in GRChombo's root directory.
+ */
+
+#ifndef INITIALFLUIDDATA_HPP_
+#define INITIALFLUIDDATA_HPP_
+
+#include "ADMConformalVars.hpp"
+#include "Cell.hpp"
+#include "Coordinates.hpp"
+#include "FluidCCZ4RHS.hpp"
+#include "PerfectFluid.hpp"
+#include "Tensor.hpp"
+#include "UserVariables.hpp" //This files needs NUM_VARS - total no. components
+#include "VarsTools.hpp"
+#include "simd.hpp"
+
+//! Class which sets the initial fluid matter config
+// aka Kevin data
+class InitialFluidData
+{
+    // Now the non grid ADM vars
+    template <class data_t>
+    using MetricVars = ADMConformalVars::VarsWithGauge<data_t>;
+
+  public:
+    //! A structure for the input params for fluid properties and initial
+    //! conditions
+    struct params_t
+    {
+        std::array<double, CH_SPACEDIM>
+            center; //!< Centre of perturbation in initial SF bubble
+        double L;
+        double rho0;
+        double uflow;
+        double amplitude;
+        double awidth;
+        double asigma;
+        std::array<double, 2> ycenter;
+    };
+
+    //! The constructor
+    InitialFluidData(params_t a_params, double a_dx)
+        : m_dx(a_dx), m_params(a_params)
+    {
+    }
+
+    //! Function to compute the value of all the initial vars on the grid
+    template <class data_t> void compute(Cell<data_t> current_cell) const
+    {
+        // where am i?
+        Coordinates<data_t> coords(current_cell, m_dx, m_params.center);
+
+        const auto metric_vars = current_cell.template load_vars<MetricVars>();
+
+        data_t x = coords.x;
+        double y = coords.y;
+        double z = coords.z;
+
+        Tensor<1, data_t> vi, Sj;
+        data_t chi_regularised = simd_max(metric_vars.chi, 1e-6);
+
+        vi[0] = m_params.uflow *
+                (tanh((y - m_params.ycenter[0]) / m_params.awidth) -
+                 tanh((y - m_params.ycenter[1]) / m_params.awidth) - 1.);
+        vi[1] = m_params.amplitude * sin(2. * M_PI * x / m_params.L) *
+                (exp(-pow((y - m_params.ycenter[0]) / m_params.asigma, 2)) +
+                 exp(-pow((y - m_params.ycenter[1]) / m_params.asigma, 2)));
+        vi[2] = 0.;
+
+        data_t nn =
+            1. + 0.5 * (tanh((y - m_params.ycenter[0]) / m_params.awidth) -
+                        tanh((y - m_params.ycenter[1]) / m_params.awidth));
+        data_t eps = 0.;
+
+        data_t rho = m_params.rho0;
+        data_t v2 = 0.;
+        FOR(i, j) v2 += metric_vars.h[i][j] * vi[i] * vi[j] / chi_regularised;
+        data_t P = rho * (1. + eps) / 3.;
+        data_t WW = 1. / (1. - v2);
+        data_t hh = 1. + eps + P / rho;
+
+        data_t D = rho * sqrt(WW);
+        data_t tau = rho * hh * WW - P - D;
+        FOR(i)
+        {
+            Sj[i] = 0.;
+            FOR(j)
+            Sj[i] +=
+                rho * hh * WW * metric_vars.h[i][j] * vi[j] / chi_regularised;
+        }
+
+        // store the vars
+        current_cell.store_vars(rho, c_rho);
+        current_cell.store_vars(vi, GRInterval<c_vi1, c_vi3>());
+        current_cell.store_vars(D, c_D);
+        current_cell.store_vars(Sj, GRInterval<c_Sj1, c_Sj3>());
+        current_cell.store_vars(tau, c_tau);
+    }
+
+  protected:
+    double m_dx;
+    const params_t m_params;
+};
+
+#endif /* INITIALFLUIDDATA_HPP_ */