From 706134bb0d5e7b0270e964fceab9f89adf7ff305 Mon Sep 17 00:00:00 2001 From: dinatraykova Date: Thu, 31 Aug 2023 18:39:23 +0200 Subject: [PATCH] Add shock initial data --- Examples/PerfectFluid/ObliqueShockID.hpp | 62 ++++++++++---- .../ObliqueShockID_SimulationParameters.hpp | 85 +++++++++++++++++++ 2 files changed, 132 insertions(+), 15 deletions(-) create mode 100644 Examples/PerfectFluid/ObliqueShockID_SimulationParameters.hpp diff --git a/Examples/PerfectFluid/ObliqueShockID.hpp b/Examples/PerfectFluid/ObliqueShockID.hpp index da06ff695..c06346345 100644 --- a/Examples/PerfectFluid/ObliqueShockID.hpp +++ b/Examples/PerfectFluid/ObliqueShockID.hpp @@ -31,9 +31,12 @@ class InitialFluidData { std::array 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,8 +50,6 @@ class InitialFluidData { // where am i? Coordinates 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(); @@ -56,26 +57,57 @@ class InitialFluidData 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 + // 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.; diff --git a/Examples/PerfectFluid/ObliqueShockID_SimulationParameters.hpp b/Examples/PerfectFluid/ObliqueShockID_SimulationParameters.hpp new file mode 100644 index 000000000..5dd55a713 --- /dev/null +++ b/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_ */