Add fluid matter rhs class

Examples/PerfectFluid/FluidCCZ4RHS.hpp

@@ -0,0 +1,118 @@
+/* GRChombo
+ * Copyright 2012 The GRChombo collaboration.
+ * Please refer to LICENSE in GRChombo's root directory.
+ */
+
+#ifndef FLUIDCCZ4RHS_HPP_
+#define FLUIDCCZ4RHS_HPP_
+
+#include "CCZ4Geometry.hpp"
+#include "CCZ4RHS.hpp"
+#include "Cell.hpp"
+#include "FourthOrderDerivatives.hpp"
+#include "WENODerivatives.hpp"
+#include "MovingPunctureGauge.hpp"
+#include "Tensor.hpp"
+#include "TensorAlgebra.hpp"
+#include "UserVariables.hpp" //This files needs NUM_VARS - total number of components
+#include "VarsTools.hpp"
+#include "simd.hpp"
+
+//!  Calculates RHS using CCZ4 including matter terms, and matter variable
+//!  evolution
+/*!
+     The class calculates the RHS evolution for all the variables. It inherits
+   from the CCZ4RHS class, which it uses to do the non matter evolution of
+   variables. It then adds in the additional matter terms to the CCZ4 evolution
+   (those including the stress energy tensor), and calculates the evolution of
+   the matter variables. It does not assume a specific form of matter but is
+   templated over a matter class matter_t. Please see the class ScalarField as
+   an example of a matter_t. \sa CCZ4RHS(), ScalarField()
+*/
+
+template <class matter_t, class gauge_t = MovingPunctureGauge,
+          class deriv_t = FourthOrderDerivatives>
+class FluidCCZ4RHS : public CCZ4RHS<gauge_t, deriv_t>
+{
+  public:
+    // Use this alias for the same template instantiation as this class
+    using CCZ4 = CCZ4RHS<gauge_t, deriv_t>;
+
+    using params_t = CCZ4_params_t<typename gauge_t::params_t>;
+
+    template <class data_t>
+    using MatterVars = typename matter_t::template Vars<data_t>;
+
+    template <class data_t>
+    using MatterDiff2Vars = typename matter_t::template Diff2Vars<data_t>;
+
+    template <class data_t>
+    using CCZ4Vars = typename CCZ4::template Vars<data_t>;
+
+    template <class data_t>
+    using CCZ4Diff2Vars = typename CCZ4::template Diff2Vars<data_t>;
+
+    // Inherit the variable definitions from CCZ4RHS + matter_t
+    template <class data_t>
+    struct Vars : public CCZ4Vars<data_t>, public MatterVars<data_t>
+    {
+        /// Defines the mapping between members of Vars and Chombo grid
+        /// variables (enum in User_Variables)
+        template <typename mapping_function_t>
+        void enum_mapping(mapping_function_t mapping_function)
+        {
+            CCZ4Vars<data_t>::enum_mapping(mapping_function);
+            MatterVars<data_t>::enum_mapping(mapping_function);
+        }
+    };
+
+    template <class data_t>
+    struct Diff2Vars : public CCZ4Diff2Vars<data_t>,
+                       public MatterDiff2Vars<data_t>
+    {
+        /// Defines the mapping between members of Vars and Chombo grid
+        /// variables (enum in User_Variables)
+        template <typename mapping_function_t>
+        void enum_mapping(mapping_function_t mapping_function)
+        {
+            CCZ4Diff2Vars<data_t>::enum_mapping(mapping_function);
+            MatterDiff2Vars<data_t>::enum_mapping(mapping_function);
+        }
+    };
+
+    //!  Constructor of class MatterCCZ4
+    /*!
+       Inputs are the grid spacing, plus the CCZ4 evolution parameters and a
+       matter object. It also takes the dissipation parameter sigma, and allows
+       the formulation to be toggled between CCZ4 and BSSN. The default is CCZ4.
+       It allows the user to set the value of Newton's constant, which is set to
+       one by default.
+    */
+    FluidCCZ4RHS(matter_t a_matter, params_t a_params, double a_dx,
+                  double a_sigma, int a_formulation = CCZ4RHS<>::USE_CCZ4,
+                  double a_G_Newton = 1.0);
+
+    //!  The compute member which calculates the RHS at each point in the box
+    //!  \sa matter_rhs_equation()
+    template <class data_t> void compute(Cell<data_t> current_cell) const;
+
+  protected:
+    //! The function which adds in the EM Tensor terms to the CCZ4 rhs \sa
+    //! compute()
+    template <class data_t>
+    void add_emtensor_rhs(
+        Vars<data_t>
+            &matter_rhs, //!< the RHS data for each variable at that point.
+        const Vars<data_t> &vars, //!< the value of the variables at the point.
+        const Vars<Tensor<1, data_t>>
+            &d1 //!< the value of the first derivatives of the variables.
+    ) const;
+
+    // Class members
+    matter_t my_matter;      //!< The matter object, e.g. a scalar field.
+    const double m_G_Newton; //!< Newton's constant, set to one by default.
+};
+
+#include "FluidCCZ4RHS.impl.hpp"
+
+#endif /* FLUIDCCZ4RHS_HPP_ */

Examples/PerfectFluid/FluidCCZ4RHS.impl.hpp

@@ -0,0 +1,120 @@
+/* GRChombo
+ * Copyright 2012 The GRChombo collaboration.
+ * Please refer to LICENSE in GRChombo's root directory.
+ */
+
+#if !defined(FLUIDCCZ4RHS_HPP_)
+#error "This file should only be included through FluidCCZ4RHS.hpp"
+#endif
+
+#ifndef FLUIDCCZ4RHS_IMPL_HPP_
+#define FLUIDCCZ4RHS_IMPL_HPP_
+#include "DimensionDefinitions.hpp"
+
+template <class matter_t, class gauge_t, class deriv_t>
+FluidCCZ4RHS<matter_t, gauge_t, deriv_t>::FluidCCZ4RHS(
+    matter_t a_matter, CCZ4_params_t<typename gauge_t::params_t> a_params,
+    double a_dx, double a_sigma, int a_formulation, double a_G_Newton)
+    : CCZ4RHS<gauge_t, deriv_t>(a_params, a_dx, a_sigma, a_formulation,
+                                0.0 /*No cosmological constant*/),
+      my_matter(a_matter), m_G_Newton(a_G_Newton)
+{
+}
+
+template <class matter_t, class gauge_t, class deriv_t>
+template <class data_t>
+void FluidCCZ4RHS<matter_t, gauge_t, deriv_t>::compute(
+    Cell<data_t> current_cell) const
+{
+    // copy data from chombo gridpoint into local variables
+    const auto matter_vars = current_cell.template load_vars<Vars>();
+    const auto d1 = this->m_deriv.template diff1<Vars>(current_cell);
+    const auto d2 = this->m_deriv.template diff2<Diff2Vars>(current_cell);
+    const auto advec =
+        this->m_deriv.template advection<Vars>(current_cell, matter_vars.shift);
+    ///const auto lm = this->m_weno.template get_Pface<Vars>
+      //      (current_cell, WENODerivatives::LEFT_MINUS);
+    //  (current_cell,0);
+    //const auto lp = this->m_weno.template get_Pface<Vars>
+      // (current_cell, WENODerivatives::LEFT_PLUS);
+    // (current_cell,1);
+    //const auto rm = this->m_weno.template get_Pface<Vars>
+      // (current_cell, WENODerivatives::RIGHT_MINUS);
+      //(current_cell,2);
+    //const auto rp = this->m_weno.template get_Pface<Vars>
+      //      (current_cell, WENODerivatives::RIGHT_PLUS);
+      //(current_cell,3);
+
+    // Call CCZ4 RHS - work out RHS without matter, no dissipation
+    Vars<data_t> matter_rhs;
+    this->rhs_equation(matter_rhs, matter_vars, d1, d2, advec);
+
+    // add RHS matter terms from EM Tensor
+    add_emtensor_rhs(matter_rhs, matter_vars, d1);
+
+    // add evolution of matter fields themselves
+    my_matter.add_matter_rhs(matter_rhs, matter_vars, d1, d1, d1, d1);
+
+    // Add dissipation to all terms
+    this->m_deriv.add_dissipation(matter_rhs, current_cell, this->m_sigma);
+
+    // Write the rhs into the output FArrayBox
+    current_cell.store_vars(matter_rhs);
+}
+
+// Function to add in EM Tensor matter terms to CCZ4 rhs
+template <class matter_t, class gauge_t, class deriv_t>
+template <class data_t>
+void FluidCCZ4RHS<matter_t, gauge_t, deriv_t>::add_emtensor_rhs(
+    Vars<data_t> &matter_rhs, const Vars<data_t> &matter_vars,
+    const Vars<Tensor<1, data_t>> &d1) const
+{
+    using namespace TensorAlgebra;
+
+    const auto h_UU = compute_inverse_sym(matter_vars.h);
+    const auto chris = compute_christoffel(d1.h, h_UU);
+
+    // Calculate elements of the decomposed stress energy tensor
+    const auto emtensor =
+        my_matter.compute_emtensor(matter_vars, d1, h_UU, chris.ULL);
+
+    // Update RHS for K and Theta depending on formulation
+    if (this->m_formulation == CCZ4RHS<>::USE_BSSN)
+    {
+        matter_rhs.K += 4.0 * M_PI * m_G_Newton * matter_vars.lapse *
+                        (emtensor.S + emtensor.rho);
+        matter_rhs.Theta += 0.0;
+    }
+    else
+    {
+        matter_rhs.K += 4.0 * M_PI * m_G_Newton * matter_vars.lapse *
+                        (emtensor.S - 3 * emtensor.rho);
+        matter_rhs.Theta +=
+            -8.0 * M_PI * m_G_Newton * matter_vars.lapse * emtensor.rho;
+    }
+
+    // Update RHS for other variables
+    Tensor<2, data_t> Sij_TF = emtensor.Sij;
+    make_trace_free(Sij_TF, matter_vars.h, h_UU);
+
+    FOR(i, j)
+    {
+        matter_rhs.A[i][j] += -8.0 * M_PI * m_G_Newton * matter_vars.chi *
+                              matter_vars.lapse * Sij_TF[i][j];
+    }
+
+    FOR(i)
+    {
+        data_t matter_term_Gamma = 0.0;
+        FOR(j)
+        {
+            matter_term_Gamma += -16.0 * M_PI * m_G_Newton * matter_vars.lapse *
+                                 h_UU[i][j] * emtensor.Si[j];
+        }
+
+        matter_rhs.Gamma[i] += matter_term_Gamma;
+        matter_rhs.B[i] += matter_term_Gamma;
+    }
+}
+
+#endif /* FLUIDCCZ4RHS_IMPL_HPP_ */

Examples/PerfectFluid/InitialFluidData.hpp

@@ -9,7 +9,7 @@
 #include "ADMConformalVars.hpp"
 #include "Cell.hpp"
 #include "Coordinates.hpp"
-#include "MatterCCZ4RHS.hpp"
+#include "FluidCCZ4RHS.hpp"
 #include "PerfectFluid.hpp"
 #include "Tensor.hpp"
 #include "UserVariables.hpp" //This files needs NUM_VARS - total no. components
@@ -53,9 +53,8 @@ class InitialFluidData
         // where am i?
         Coordinates<data_t> coords(current_cell, m_dx, m_params.center);
 
-	//	MetricVars<data_t> metric_vars;
-	//m_background.compute_metric_background(metric_vars, coords);
-
+	MetricVars<data_t> metric_vars;
+
 	data_t x = coords.x;
 	double y = coords.y;
 	double z = coords.z;

Examples/PerfectFluid/PerfectFluid.hpp

@@ -102,9 +102,10 @@ template <class eos_t = DefaultEoS> class PerfectFluid
     void add_matter_rhs(
         rhs_vars_t<data_t> &rhs,             //!< value of the RHS for all vars
         const vars_t<data_t> &vars,          //!< value of the variables
-        const vars_t<Tensor<1, data_t>> &d1, //!< value of the 1st derivs
-        const diff2_vars_t<Tensor<2, data_t>> &d2, //!< value of the 2nd derivs
-        const vars_t<data_t> &advec)
+        const vars_t<Tensor<1, data_t>> &lm, //!< value of the left_minus weno
+	const vars_t<Tensor<1, data_t>> &lp, //!< value of the left_plus weno
+	const vars_t<Tensor<1, data_t>> &rm, //!< value of the right_minus weno
+	const vars_t<Tensor<1, data_t>> &rp) //!< value of the right_plus weno
         const; //!< the value of the advection terms
 
 };

Examples/PerfectFluid/PerfectFluid.impl.hpp

@@ -70,9 +70,10 @@ template <class data_t, template <typename> class vars_t,
           template <typename> class rhs_vars_t>
 void PerfectFluid<eos_t>::add_matter_rhs(
     rhs_vars_t<data_t> &rhs, const vars_t<data_t> &vars,
-    const vars_t<Tensor<1, data_t>> &d1,
-    const diff2_vars_t<Tensor<2, data_t>> &d2,
-    const vars_t<data_t> &advec) const
+    const vars_t<Tensor<1, data_t>> &lm,
+    const vars_t<Tensor<1, data_t>> &lp,
+    const vars_t<Tensor<1, data_t>> &rm,
+    const vars_t<Tensor<1, data_t>> &rp) const
 {
     using namespace TensorAlgebra;
 

Examples/PerfectFluid/PerfectFluidLevel.cpp

@@ -13,7 +13,7 @@
 #include "TraceARemoval.hpp"
 
 // For RHS update
-#include "MatterCCZ4RHS.hpp"
+#include "FluidCCZ4RHS.hpp"
 
 // For constraints calculation
 #include "NewMatterConstraints.hpp"
@@ -54,6 +54,8 @@ void PerfectFluidLevel::initialData()
 
     // First set everything to zero then initial conditions for scalar field -
     // here a Kerr BH and a scalar field profile
+    //    EoS eos(m_p.eos_params);
+    //PerfectFluidEoS perfect_fluid(eos);
     BoxLoops::loop(
         make_compute_pack(SetValue(0.), Minkowski(m_p.kerr_params, m_dx),
                           InitialFluidData(m_p.initial_params, m_dx)),
@@ -93,16 +95,18 @@ void PerfectFluidLevel::specificEvalRHS(GRLevelData &a_soln, GRLevelData &a_rhs,
     PerfectFluidEoS perfect_fluid(eos);
     if (m_p.max_spatial_derivative_order == 4)
     {
-        MatterCCZ4RHS<PerfectFluidEoS, MovingPunctureGauge,
-                      FourthOrderDerivatives>
+        FluidCCZ4RHS<PerfectFluidEoS, MovingPunctureGauge,
+		     //	     FourthOrderDerivatives, WENODerivatives>
+		     FourthOrderDerivatives>
             my_ccz4_matter(perfect_fluid, m_p.ccz4_params, m_dx, m_p.sigma,
                            m_p.formulation, m_p.G_Newton);
         BoxLoops::loop(my_ccz4_matter, a_soln, a_rhs, EXCLUDE_GHOST_CELLS);
     }
     else if (m_p.max_spatial_derivative_order == 6)
     {
-        MatterCCZ4RHS<PerfectFluidEoS, MovingPunctureGauge,
-                      SixthOrderDerivatives>
+        FluidCCZ4RHS<PerfectFluidEoS, MovingPunctureGauge,
+		     // SixthOrderDerivatives, WENODerivatives>
+		     FourthOrderDerivatives>
             my_ccz4_matter(perfect_fluid, m_p.ccz4_params, m_dx, m_p.sigma,
                            m_p.formulation, m_p.G_Newton);
         BoxLoops::loop(my_ccz4_matter, a_soln, a_rhs, EXCLUDE_GHOST_CELLS);