Add shock initial data

Examples/PerfectFluid/ObliqueShockID.hpp

@@ -31,9 +31,12 @@ class InitialFluidData
     {
         std::array<double, CH_SPACEDIM>
             center; //!< Centre of perturbation in initial SF bubble
-        double rho0;
-        double delta;
-        double width; //!< Width of bump in initial SF bubble
+        double emax;
+        double emin;
+        double vx_in;
+        double vy_in;
+        double nn_in;
+        double L;
     };
 
     //! The constructor
@@ -47,35 +50,64 @@ class InitialFluidData
     {
         // 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;
+        const double gamma = 10.;
+        double one_o_gamma = 1. / gamma;
+        const double delta = 10.;
+        data_t de, dvx, dvy, dn;
+        data_t x_p_L, x_m_L;
+        double y_p_L, y_m_L;
+
         data_t chi_regularised = simd_max(metric_vars.chi, 1e-6);
 
-        vi[0] = 0.;
-        vi[1] = 0.;
+        // Transition function
+        // template <class data_t>
+        // data_t transD(data_t d, double delta)
+        //{
+        //  return 0.5*(1.+tanh(d/delta));
+        //}
+
+        // data_t rho = m_params.emax - m_params.emin*transD(de,delta);
+        data_t rho =
+            m_params.emax - m_params.emin * 0.5 * (1. + tanh(de / delta));
+        x_p_L = pow(x + m_params.L, gamma);
+        x_m_L = pow(x - m_params.L, gamma);
+        y_p_L = pow(y + m_params.L, gamma);
+        y_m_L = pow(y - m_params.L, gamma);
+
+        de = m_params.L - pow(x_m_L + y_m_L, one_o_gamma);
+        dvx = m_params.L - pow(x_p_L + y_m_L, one_o_gamma);
+        dvy = m_params.L - pow(x_m_L + y_p_L, one_o_gamma);
+        dn = m_params.L - pow(x_p_L + y_p_L, one_o_gamma);
+
+        Tensor<1, data_t> vi;
+        // vi[0] = m_params.vx_in*transD(dvx,delta);;
+        vi[0] = m_params.vx_in * 0.5 * (1. + tanh(dvx / delta));
+        // vi[1] = m_params.vy_in*transD(dvy,delta);;
+        vi[1] = m_params.vy_in * 0.5 * (1. + tanh(dvy / delta));
         vi[2] = 0.;
 
-        // calculate the field value
-        data_t rho0 =
-            m_params.amplitude * (exp(-pow(rr / m_params.width, 2.0))) +
-            m_params.delta;
+        // data_t nn = m_params.nn_in*transD(dn,delta) + 0.1;
+        data_t nn = m_params.nn_in * 0.5 * (1. + tanh(dn / delta)) + 0.1;
+
+        data_t eps = 0.;
         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 = rho0 * (1. + eps) / 3.;
+
+        data_t P = rho * (1. + eps) / 3.;
         data_t WW = 1. / (1. - v2);
-        data_t hh = 1. + eps + P / rho0;
+        data_t hh = 1. + eps + P / rho;
 
         data_t D = rho * sqrt(WW);
         data_t tau = rho * hh * WW - P - D;
+
+        Tensor<1, data_t> Sj;
         FOR(i)
         {
             Sj[i] = 0.;

Examples/PerfectFluid/ObliqueShockID_SimulationParameters.hpp

@@ -0,0 +1,85 @@
+/* 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
+        initial_params.L = L;
+        pp.load("G_Newton", G_Newton,
+                0.0); // for now the example neglects backreaction
+        pp.load("fluid_emax", initial_params.emax, 3.);
+        pp.load("fluid_emin", initial_params.emin, 2.97);
+        pp.load("fluid_vx_in", initial_params.vx_in, 0.97);
+        pp.load("fluid_vy_in", initial_params.vy_in, 0.97);
+        pp.load("fluid_nn_in", initial_params.nn_in, 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_ */