From 4a389aab89984467e50b3e78d45612d375943e8b Mon Sep 17 00:00:00 2001 From: llibert94 Date: Fri, 28 Jul 2023 18:28:07 +0100 Subject: [PATCH] Add fluid rhs --- Examples/PerfectFluid/Fluxes.hpp | 61 +++++++++++++++++++ Examples/PerfectFluid/PerfectFluid.hpp | 10 ++- Examples/PerfectFluid/PerfectFluid.impl.hpp | 55 ++++++++++++++--- Examples/PerfectFluid/PerfectFluidLevel.cpp | 4 +- Examples/PerfectFluid/PrimitiveRecovery.hpp | 4 +- .../PerfectFluid/SimulationParameters.hpp | 2 + Examples/PerfectFluid/Sources.hpp | 55 +++++++++++++++++ 7 files changed, 177 insertions(+), 14 deletions(-) create mode 100644 Examples/PerfectFluid/Fluxes.hpp create mode 100644 Examples/PerfectFluid/Sources.hpp diff --git a/Examples/PerfectFluid/Fluxes.hpp b/Examples/PerfectFluid/Fluxes.hpp new file mode 100644 index 000000000..5e021648a --- /dev/null +++ b/Examples/PerfectFluid/Fluxes.hpp @@ -0,0 +1,61 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef FLUXES_HPP_ +#define FLUXES_HPP_ + +#include "DimensionDefinitions.hpp" +#include "Tensor.hpp" +#include "TensorAlgebra.hpp" +#include "simd.hpp" + +namespace Fluxes +{ +// Primitive to conservative variables +template class vars_t> +vars_t compute_flux(const vars_t &vars, const int idir, const double lambda) +{ + const auto h_UU = TensorAlgebra::compute_inverse_sym(vars.h); + vars_t out; + out.D = (vars.lapse * vars.vi[idir] - vars.shift[idir]) * vars.D - lambda * vars.D; + Tensor<1, data_t> Sj_U; + FOR(i) { + Sj_U[i] = 0; + FOR(j) Sj_U[i] += vars.chi * h_UU[i][j] * vars.Sj[j]; + } + + data_t chi_regularised = simd_max(1e-6, vars.chi); + Tensor<1, data_t> vi_D; + FOR(i) + { + vi_D[i] = 0.; + FOR(j) + { + vi_D[i] += vars.h[i][j] * vars.vi[j] / chi_regularised; + } + } + + data_t v2 = 0.; + FOR(i) v2 += vars.vi[i] * vi_D[i]; + + data_t P = vars.rho * (1. + vars.eps) / 3.; //for now we assume a conformal fluid + data_t WW = 1. / (1. - v2); + data_t hh = 1. + vars.eps + P / vars.rho; + + FOR(j) + { + out.Sj[j] = vars.lapse * vars.rho * hh * WW * vars.vi[idir] * vi_D[j] - + vars.shift[idir] * vars.Sj[j] - lambda * vars.Sj[j]; + FOR(k) out.Sj[j] += vars.lapse * P * h_UU[idir][k] * vars.h[j][k]; + } + + out.tau = vars.lapse * (Sj_U[idir] - vars.D * vars.vi[idir]) - + vars.shift[idir] * vars.tau - lambda * vars.tau; + + return out; +} + +} +#endif /* FLUXES_HPP_ */ diff --git a/Examples/PerfectFluid/PerfectFluid.hpp b/Examples/PerfectFluid/PerfectFluid.hpp index 0a7e8641a..94709fb21 100644 --- a/Examples/PerfectFluid/PerfectFluid.hpp +++ b/Examples/PerfectFluid/PerfectFluid.hpp @@ -9,12 +9,15 @@ #include "CCZ4Geometry.hpp" #include "DefaultEoS.hpp" #include "DimensionDefinitions.hpp" +#include "Fluxes.hpp" #include "FourthOrderDerivatives.hpp" -#include "WENODerivatives.hpp" +#include "PrimitiveRecovery.hpp" +#include "Sources.hpp" #include "Tensor.hpp" #include "TensorAlgebra.hpp" #include "UserVariables.hpp" //This files needs NUM_VARS, total num of components #include "VarsTools.hpp" +#include "WENODerivatives.hpp" //! Calculates the matter type specific elements such as the EMTensor and // matter evolution @@ -34,11 +37,14 @@ template class PerfectFluid { protected: //! The local copy of the potential + double m_dx; + double m_lambda; eos_t my_eos; public: //! Constructor of class PerfectFluid, inputs are the matter parameters. - PerfectFluid(const eos_t a_eos) : my_eos(a_eos) {} + PerfectFluid(double a_dx, const double a_lambda,const eos_t a_eos) + : m_dx(a_dx), m_lambda(a_lambda), my_eos(a_eos) {} //! Structure containing the rhs variables for the matter fields template struct Vars diff --git a/Examples/PerfectFluid/PerfectFluid.impl.hpp b/Examples/PerfectFluid/PerfectFluid.impl.hpp index 2a5561460..400bda7cf 100644 --- a/Examples/PerfectFluid/PerfectFluid.impl.hpp +++ b/Examples/PerfectFluid/PerfectFluid.impl.hpp @@ -86,17 +86,56 @@ void PerfectFluid::add_matter_rhs( data_t dPdrho = 0.0; my_eos.compute_eos(P_of_rho, dPdrho, vars); + vars_t source = Sources::compute_source(vars, lp); // evolution equations for the fluid conservative variables - rhs.D = 0.0; - rhs.tau = 0.0; - rhs.rho = 0.0; - rhs.eps = 0.0; + rhs.D = source.D; + FOR(i) rhs.Sj[i] = source.Sj[i]; + rhs.tau = source.tau; - FOR(i) - { - rhs.Sj[i] = 0.0; - rhs.vi[i] = 0.0; + FOR(idir) { + vars_t vars_right_p; + vars_right_p.rho = rp.rho[idir]; + vars_right_p.eps = rp.eps[idir]; + FOR(j) vars_right_p.vi[j] = rp.vi[j][idir]; + PrimitiveRecovery::PtoC(vars_right_p); + vars_t flux_right_p = Fluxes::compute_flux(vars_right_p, idir, m_lambda); + + vars_t vars_right_m; + vars_right_m.rho = rm.rho[idir]; + vars_right_m.eps = rm.eps[idir]; + FOR(j) vars_right_m.vi[j] = rm.vi[j][idir]; + PrimitiveRecovery::PtoC(vars_right_m); + vars_t flux_right_m = Fluxes::compute_flux(vars_right_m, idir, m_lambda); + + rhs.D -= 1. / (2. * m_dx) * (flux_right_p.D + flux_right_m.D); + FOR(j) rhs.Sj[j] -= 1. / (2. * m_dx) * (flux_right_p.Sj[j] + flux_right_m.Sj[j]); + rhs.tau -= 1. / (2. * m_dx) * (flux_right_p.tau + flux_right_m.tau); + } + + FOR(idir) { + vars_t vars_left_p; + vars_left_p.rho = lp.rho[idir]; + vars_left_p.eps = lp.eps[idir]; + FOR(j) vars_left_p.vi[j] = lp.vi[j][idir]; + PrimitiveRecovery::PtoC(vars_left_p); + vars_t flux_left_p = Fluxes::compute_flux(vars_left_p, idir, m_lambda); + + vars_t vars_left_m; + vars_left_m.rho = lm.rho[idir]; + vars_left_m.eps = lm.eps[idir]; + FOR(j) vars_left_m.vi[j] = lm.vi[j][idir]; + PrimitiveRecovery::PtoC(vars_left_m); + vars_t flux_left_m = Fluxes::compute_flux(vars_left_m, idir, m_lambda); + + rhs.D += 1. / (2. * m_dx) * (flux_left_p.D + flux_left_m.D); + FOR(j) rhs.Sj[j] += 1. / (2. * m_dx) * (flux_left_p.Sj[j] + flux_left_m.Sj[j]); + rhs.tau += 1. / (2. * m_dx) * (flux_left_p.tau + flux_left_m.tau); } + + rhs.rho = 0.; + rhs.eps = 0.; + FOR(i) rhs.vi[i] = 0.; + } #endif /* SCALARFIELD_IMPL_HPP_ */ diff --git a/Examples/PerfectFluid/PerfectFluidLevel.cpp b/Examples/PerfectFluid/PerfectFluidLevel.cpp index 5ca0a8045..acf01aeca 100644 --- a/Examples/PerfectFluid/PerfectFluidLevel.cpp +++ b/Examples/PerfectFluid/PerfectFluidLevel.cpp @@ -72,7 +72,7 @@ void PerfectFluidLevel::prePlotLevel() { fillAllGhosts(); EoS eos(m_p.eos_params); - PerfectFluidEoS perfect_fluid(eos); + PerfectFluidEoS perfect_fluid(m_dx, m_p.lambda, eos); BoxLoops::loop( MatterConstraints( perfect_fluid, m_dx, m_p.G_Newton, c_Ham, Interval(c_Mom, c_Mom)), @@ -92,7 +92,7 @@ void PerfectFluidLevel::specificEvalRHS(GRLevelData &a_soln, GRLevelData &a_rhs, // Calculate MatterCCZ4 right hand side with matter_t = ScalarField EoS eos(m_p.eos_params); - PerfectFluidEoS perfect_fluid(eos); + PerfectFluidEoS perfect_fluid(m_dx, m_p.lambda, eos); if (m_p.max_spatial_derivative_order == 4) { FluidCCZ4RHS class vars_t> -void PtoC(const vars_t &vars) +void PtoC(vars_t &vars) { data_t chi_regularised = simd_max(1e-6, vars.chi); Tensor<1, data_t> vi_D; @@ -44,7 +44,7 @@ void PtoC(const vars_t &vars) // Conservative to primitive variables template class vars_t> -void CtoP(const vars_t &vars) +void CtoP(vars_t &vars) { const auto h_UU = TensorAlgebra::compute_inverse_sym(vars.h); data_t E = vars.tau + vars.D; diff --git a/Examples/PerfectFluid/SimulationParameters.hpp b/Examples/PerfectFluid/SimulationParameters.hpp index 0a8d921ab..c421c76f1 100644 --- a/Examples/PerfectFluid/SimulationParameters.hpp +++ b/Examples/PerfectFluid/SimulationParameters.hpp @@ -38,6 +38,7 @@ class SimulationParameters : public SimulationParametersBase pp.load("fluid_awidth", initial_params.awidth, 0.05); pp.load("fluid_sigma", initial_params.sigma, 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 @@ -75,6 +76,7 @@ class SimulationParameters : public SimulationParametersBase // 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; diff --git a/Examples/PerfectFluid/Sources.hpp b/Examples/PerfectFluid/Sources.hpp new file mode 100644 index 000000000..114aaf091 --- /dev/null +++ b/Examples/PerfectFluid/Sources.hpp @@ -0,0 +1,55 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef SOURCES_HPP_ +#define SOURCES_HPP_ + +#include "DimensionDefinitions.hpp" +#include "Tensor.hpp" +#include "TensorAlgebra.hpp" +#include "simd.hpp" + +namespace Sources +{ +// Primitive to conservative variables +template class vars_t> +vars_t compute_source(const vars_t &vars, + const vars_t> &d1) +{ + data_t chi_regularised = simd_max(1e-6, vars.chi); + const auto h_UU = TensorAlgebra::compute_inverse_sym(vars.h); + vars_t out; + + data_t v2 = 0.; + FOR(i,j) v2 += vars.h[i][j] * vars.vi[j] * vars.vi[i] / chi_regularised; + + data_t P = vars.rho * (1. + vars.eps) / 3.; //for now we assume a conformal fluid + data_t WW = 1. / (1. - v2); + data_t hh = 1. + vars.eps + P / vars.rho; + + out.D = 0.; + FOR(j) { + out.Sj[j] = - (vars.tau + vars.D) * d1.lapse[j]; + FOR(i) { + out.Sj[j] += vars.Sj[i] * d1.shift[i][j]; + FOR(k) { + out.Sj[j] += vars.lapse / 2. * (d1.h[i][k][j] - vars.h[i][k] * + d1.chi[j] / chi_regularised) * (vars.rho * hh * WW * + vars.vi[i] * vars.vi[k] / chi_regularised + P * h_UU[i][k]); + } + } + } + out.tau = 0.; + FOR(i,j) { + out.tau += vars.lapse * (vars.A[i][j] + vars.h[i][j] / 3. * vars.K) * + (vars.rho * hh * WW * vars.vi[i] * vars.vi[j] / chi_regularised + + P * h_UU[i][j]) - vars.chi * h_UU[i][j] * vars.Sj[i] * d1.lapse[j]; + } + + return out; +} + +} +#endif /* FLUXES_HPP_ */