Merge changes from the two branches

Examples/Fluid_Kerr_merged/ConservedQuantities.hpp

@@ -0,0 +1,42 @@
+/* GRChombo
+ * Copyright 2012 The GRChombo collaboration.
+ * Please refer to LICENSE in GRChombo's root directory.
+ */
+
+#ifndef CONSERVEDQUANTITIES_HPP_
+#define CONSERVEDQUANTITIES_HPP_
+
+#include "DimensionDefinitions.hpp"
+#include "Tensor.hpp"
+#include "TensorAlgebra.hpp"
+#include "simd.hpp"
+
+namespace ConservedQuantities
+{
+// Primitive to conservative variables
+template <class data_t, template <typename> class vars_t>
+void PtoC(const data_t P_of_rho, vars_t<data_t> &vars)
+{
+    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 WW = 1. / (1. - v2);
+    data_t hh = 1. + vars.eps + P_of_rho / vars.rho;
+
+    vars.D = vars.rho * sqrt(WW);
+    vars.tau = vars.rho * hh * WW - P_of_rho - vars.D;
+
+    // S_j (note lower index) = - n^a T_ai
+    FOR(i) { vars.Sj[i] = vars.rho * hh * WW * vi_D[i]; }
+}
+
+} // namespace ConservedQuantities
+#endif /* CONSERVEDQUANTITIES_HPP_ */

Examples/Fluid_Kerr_merged/DefaultEoS.hpp

@@ -0,0 +1,26 @@
+/* GRChombo
+ * Copyright 2012 The GRChombo collaboration.
+ * Please refer to LICENSE in GRChombo's root directory.
+ */
+
+#ifndef DEFAULTEOS_HPP_
+#define DEFAULTEOS_HPP_
+
+#include "Tensor.hpp"
+#include "simd.hpp"
+
+class DefaultEoS
+{
+  public:
+    //! The constructor
+    DefaultEoS() {}
+
+    //! Set the pressure of the perfect fluid here to zero
+    template <class data_t, template <typename> class vars_t>
+    void compute_eos(data_t &P_of_rho, const vars_t<data_t> &vars) const
+    {
+        P_of_rho = vars.rho * (1. + vars.eps) / 3.;
+    }
+};
+
+#endif /* DEFAULTEOS_HPP_ */

Examples/Fluid_Kerr_merged/DiagnosticVariables.hpp

@@ -0,0 +1,29 @@
+/* GRChombo
+ * Copyright 2012 The GRChombo collaboration.
+ * Please refer to LICENSE in GRChombo's root directory.
+ */
+
+#ifndef DIAGNOSTICVARIABLES_HPP
+#define DIAGNOSTICVARIABLES_HPP
+
+// assign an enum to each variable
+// todo: here add pressure, density, etc.
+enum
+{
+    c_Ham,
+    c_Mom1,
+    c_Mom2,
+    c_Mom3,
+
+    NUM_DIAGNOSTIC_VARS
+};
+
+namespace DiagnosticVariables
+{
+static const std::array<std::string, NUM_DIAGNOSTIC_VARS> variable_names = {
+    "Ham",
+
+    "Mom1", "Mom2", "Mom3"};
+}
+
+#endif /* DIAGNOSTICVARIABLES_HPP */

Examples/Fluid_Kerr_merged/EoS.hpp

@@ -0,0 +1,29 @@
+/* GRChombo
+ * Copyright 2012 The GRChombo collaboration.
+ * Please refer to LICENSE in GRChombo's root directory.
+ */
+
+#ifndef EOS_HPP_
+#define EOS_HPP_
+
+#include "simd.hpp"
+
+class EoS
+{
+  public:
+    //! The constructor
+    EoS() {}
+
+    //! Set the pressure of Gamma-2 polytrope
+    template <class data_t, template <typename> class vars_t>
+    void compute_eos(data_t &P_of_rho, const vars_t<data_t> &vars) const
+    {
+        // The pressure value as a function of rho
+        const double K = 1./3.;//100.;
+        const double n = 1.;
+        const double Gamma = 1.;// + 1. / n;
+        P_of_rho = K * pow(vars.rho, Gamma);
+    }
+};
+
+#endif /* EOS_HPP_ */

Examples/Fluid_Kerr_merged/FluidCCZ4RHS.hpp

@@ -0,0 +1,119 @@
+/* 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 "MovingPunctureGauge.hpp"
+#include "Tensor.hpp"
+#include "TensorAlgebra.hpp"
+#include "UserVariables.hpp" //This files needs NUM_VARS - total number of components
+#include "VarsTools.hpp"
+#include "WENODerivatives.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 weno_t = WENODerivatives>
+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
+    weno_t m_weno;
+    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/Fluid_Kerr_merged/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, class weno_t>
+FluidCCZ4RHS<matter_t, gauge_t, deriv_t, weno_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, class weno_t>
+template <class data_t>
+void FluidCCZ4RHS<matter_t, gauge_t, deriv_t, weno_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, lm, lp, rm, rp);
+
+    // 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, class weno_t>
+template <class data_t>
+void FluidCCZ4RHS<matter_t, gauge_t, deriv_t, weno_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_ */