From 94f6b2b9873f6737e516305da8c72b7dd990b760 Mon Sep 17 00:00:00 2001 From: dinatraykova Date: Thu, 31 Aug 2023 17:43:30 +0200 Subject: [PATCH] Add different initial data options and params --- Examples/PerfectFluid/GaussianID.hpp | 99 ++++++++++++++++ .../GaussianID_SimulationParameters.hpp | 82 ++++++++++++++ Examples/PerfectFluid/InitialFluidData.hpp | 9 +- Examples/PerfectFluid/KelvinHelmhotzID.hpp | 107 ++++++++++++++++++ .../KelvinHelmhotzID_SimulationParameters.hpp | 86 ++++++++++++++ .../PerfectFluid/SimulationParameters.hpp | 5 +- Examples/PerfectFluid/WENODerivatives.hpp | 16 +-- 7 files changed, 391 insertions(+), 13 deletions(-) create mode 100644 Examples/PerfectFluid/GaussianID.hpp create mode 100644 Examples/PerfectFluid/GaussianID_SimulationParameters.hpp create mode 100644 Examples/PerfectFluid/KelvinHelmhotzID.hpp create mode 100644 Examples/PerfectFluid/KelvinHelmhotzID_SimulationParameters.hpp diff --git a/Examples/PerfectFluid/GaussianID.hpp b/Examples/PerfectFluid/GaussianID.hpp new file mode 100644 index 000000000..5b9e9d716 --- /dev/null +++ b/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 + using MetricVars = ADMConformalVars::VarsWithGauge; + + public: + //! A structure for the input params for fluid properties and initial + //! conditions + struct params_t + { + std::array + 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 void compute(Cell current_cell) const + { + // 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(); + + 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()); + current_cell.store_vars(D, c_D); + current_cell.store_vars(Sj, GRInterval()); + current_cell.store_vars(tau, c_tau); + } + + protected: + double m_dx; + const params_t m_params; +}; + +#endif /* INITIALFLUIDDATA_HPP_ */ diff --git a/Examples/PerfectFluid/GaussianID_SimulationParameters.hpp b/Examples/PerfectFluid/GaussianID_SimulationParameters.hpp new file mode 100644 index 000000000..687935eae --- /dev/null +++ b/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_ */ diff --git a/Examples/PerfectFluid/InitialFluidData.hpp b/Examples/PerfectFluid/InitialFluidData.hpp index 784dcaea3..d4d49b546 100644 --- a/Examples/PerfectFluid/InitialFluidData.hpp +++ b/Examples/PerfectFluid/InitialFluidData.hpp @@ -31,11 +31,12 @@ class InitialFluidData { std::array center; //!< Centre of perturbation in initial SF bubble + double L; double rho0; double uflow; double amplitude; double awidth; - double sigma; + double asigma; std::array 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 = diff --git a/Examples/PerfectFluid/KelvinHelmhotzID.hpp b/Examples/PerfectFluid/KelvinHelmhotzID.hpp new file mode 100644 index 000000000..d4d49b546 --- /dev/null +++ b/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 + using MetricVars = ADMConformalVars::VarsWithGauge; + + public: + //! A structure for the input params for fluid properties and initial + //! conditions + struct params_t + { + std::array + center; //!< Centre of perturbation in initial SF bubble + double L; + double rho0; + double uflow; + double amplitude; + double awidth; + double asigma; + std::array 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 void compute(Cell current_cell) const + { + // where am i? + Coordinates coords(current_cell, m_dx, m_params.center); + + const auto metric_vars = current_cell.template load_vars(); + + 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()); + current_cell.store_vars(D, c_D); + current_cell.store_vars(Sj, GRInterval()); + current_cell.store_vars(tau, c_tau); + } + + protected: + double m_dx; + const params_t m_params; +}; + +#endif /* INITIALFLUIDDATA_HPP_ */ diff --git a/Examples/PerfectFluid/KelvinHelmhotzID_SimulationParameters.hpp b/Examples/PerfectFluid/KelvinHelmhotzID_SimulationParameters.hpp new file mode 100644 index 000000000..296290621 --- /dev/null +++ b/Examples/PerfectFluid/KelvinHelmhotzID_SimulationParameters.hpp @@ -0,0 +1,86 @@ +/* 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_rho", initial_params.rho0, 1.0); + pp.load("fluid_uflow", initial_params.uflow, 1. / (4. * sqrt(3))); + pp.load("fluid_amplitude", initial_params.amplitude, 0.01); + pp.load("fluid_width", initial_params.awidth, 0.05); + pp.load("fluid_sigma", initial_params.asigma, 0.2); + pp.load("fluid_ycenter", initial_params.ycenter, {-0.5, 0.5}); + 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_ */ diff --git a/Examples/PerfectFluid/SimulationParameters.hpp b/Examples/PerfectFluid/SimulationParameters.hpp index ebb1c5c69..296290621 100644 --- a/Examples/PerfectFluid/SimulationParameters.hpp +++ b/Examples/PerfectFluid/SimulationParameters.hpp @@ -30,13 +30,14 @@ class SimulationParameters : public SimulationParametersBase // 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_rho", initial_params.rho0, 1.0); pp.load("fluid_uflow", initial_params.uflow, 1. / (4. * sqrt(3))); pp.load("fluid_amplitude", initial_params.amplitude, 0.01); - pp.load("fluid_awidth", initial_params.awidth, 0.05); - pp.load("fluid_sigma", initial_params.sigma, 0.2); + pp.load("fluid_width", initial_params.awidth, 0.05); + pp.load("fluid_sigma", initial_params.asigma, 0.2); pp.load("fluid_ycenter", initial_params.ycenter, {-0.5, 0.5}); pp.load("lambda", lambda, 1.); // eigenvalue for numerical flux pp.load("eos_w", eos_params.eos_w, 1. / 3.); diff --git a/Examples/PerfectFluid/WENODerivatives.hpp b/Examples/PerfectFluid/WENODerivatives.hpp index aa852129e..ff958975c 100644 --- a/Examples/PerfectFluid/WENODerivatives.hpp +++ b/Examples/PerfectFluid/WENODerivatives.hpp @@ -135,14 +135,16 @@ class WENODerivatives const auto in_index = current_cell.get_in_index(); const auto strides = current_cell.get_box_pointers().m_in_stride; vars_t> d1; - d1.enum_mapping([&](const int &ivar, Tensor<1, data_t> &var) { - FOR(idir) + d1.enum_mapping( + [&](const int &ivar, Tensor<1, data_t> &var) { - var[idir] = get_Pface( - current_cell.get_box_pointers().m_in_ptr[ivar], in_index, - strides[idir], dir_switch); - } - }); + FOR(idir) + { + var[idir] = get_Pface( + current_cell.get_box_pointers().m_in_ptr[ivar], + in_index, strides[idir], dir_switch); + } + }); return d1; }