Add fluid rhs

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 data_t, template <typename> class vars_t>
+vars_t<data_t> compute_flux(const vars_t<data_t> &vars, const int idir, const double lambda)
+{
+    const auto h_UU = TensorAlgebra::compute_inverse_sym(vars.h);
+    vars_t<data_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_ */

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 eos_t = DefaultEoS> 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 <class data_t> struct Vars

Examples/PerfectFluid/PerfectFluid.impl.hpp

@@ -86,17 +86,56 @@ void PerfectFluid<eos_t>::add_matter_rhs(
     data_t dPdrho = 0.0;
     my_eos.compute_eos(P_of_rho, dPdrho, vars);
 
+    vars_t<data_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<data_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<data_t> flux_right_p = Fluxes::compute_flux(vars_right_p, idir, m_lambda);
+
+        vars_t<data_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<data_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<data_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<data_t> flux_left_p = Fluxes::compute_flux(vars_left_p, idir, m_lambda);
+
+        vars_t<data_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<data_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_ */

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<PerfectFluidEoS>(
             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<PerfectFluidEoS, MovingPunctureGauge,

Examples/PerfectFluid/PrimitiveRecovery.hpp

@@ -15,7 +15,7 @@ namespace PrimitiveRecovery
 {
 // Primitive to conservative variables
 template <class data_t, template <typename> class vars_t>
-void PtoC(const vars_t<data_t> &vars)
+void PtoC(vars_t<data_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<data_t> &vars)
 
 // Conservative to primitive variables
 template <class data_t, template <typename> class vars_t>
-void CtoP(const vars_t<data_t> &vars)
+void CtoP(vars_t<data_t> &vars)
 {
     const auto h_UU = TensorAlgebra::compute_inverse_sym(vars.h);
     data_t E = vars.tau + vars.D;

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;

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 data_t, template <typename> class vars_t>
+vars_t<data_t> compute_source(const vars_t<data_t> &vars,
+	const vars_t<Tensor<1, data_t>> &d1)
+{
+    data_t chi_regularised = simd_max(1e-6, vars.chi);
+    const auto h_UU = TensorAlgebra::compute_inverse_sym(vars.h);
+    vars_t<data_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_ */