Feature/complex sf

Source/Matter/ComplexScalarField.hpp

@@ -0,0 +1,133 @@
+/* GRChombo
+ * Copyright 2012 The GRChombo collaboration.
+ * Please refer to LICENSE in GRChombo's root directory.
+ */
+
+#ifndef COMPLEXSCALARFIELD_HPP_
+#define COMPLEXSCALARFIELD_HPP_
+
+#include "CCZ4Geometry.hpp"
+#include "DefaultPotential.hpp"
+#include "DimensionDefinitions.hpp"
+#include "FourthOrderDerivatives.hpp"
+#include "Tensor.hpp"
+#include "TensorAlgebra.hpp"
+#include "UserVariables.hpp" //This files needs NUM_VARS, total num of components
+#include "VarsTools.hpp"
+
+//!  Calculates the matter type specific elements such as the EMTensor and
+//   matter evolution
+/*!
+     This class is an example of a matter_t object which calculates the
+     matter type specific elements for the RHS update and the evaluation
+     of the constraints. This includes the Energy Momentum Tensor, and
+     the matter evolution terms. In this case, a scalar field,
+     the matter elements are phi_Re and phi_Im, and (minus) their conjugate
+     momenta, Pi_Re and Pi_Im.
+     It is templated over a potential function potential_t which the
+     user must specify in a class, although a default is provided which
+     sets dVdphi_Re, dVdphi_Im and V_of_phi to zero.
+     It assumes minimal coupling of the field to gravity.
+     \sa MatterCCZ4(), ConstraintsMatter()
+*/
+template <class potential_t = DefaultPotential> class ComplexScalarField
+{
+  protected:
+    //! The local copy of the potential
+    potential_t my_potential;
+
+  public:
+    //!  Constructor of class ComplexScalarField, inputs are the matter
+    //!  parameters.
+    ComplexScalarField(const potential_t a_potential)
+        : my_potential(a_potential)
+    {
+    }
+
+    //! Structure containing the rhs variables for the matter fields
+    template <class data_t> struct Vars
+    {
+        data_t phi_Re;
+        data_t Pi_Re;
+        data_t phi_Im;
+        data_t Pi_Im;
+
+        /// 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)
+        {
+            VarsTools::define_enum_mapping(mapping_function, c_phi_Re, phi_Re);
+            VarsTools::define_enum_mapping(mapping_function, c_Pi_Re, Pi_Re);
+            VarsTools::define_enum_mapping(mapping_function, c_phi_Im, phi_Im);
+            VarsTools::define_enum_mapping(mapping_function, c_Pi_Im, Pi_Im);
+        }
+    };
+
+    //! Structure containing the rhs variables for the matter fields requiring
+    //!  2nd derivs
+    template <class data_t> struct Diff2Vars
+    {
+        data_t phi_Re;
+        data_t phi_Im;
+
+        /// 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)
+        {
+            VarsTools::define_enum_mapping(mapping_function, c_phi_Re, phi_Re);
+            VarsTools::define_enum_mapping(mapping_function, c_phi_Im, phi_Im);
+        }
+    };
+
+    //! The function which calculates the EM Tensor, given the vars and
+    //! derivatives, including the potential
+    template <class data_t, template <typename> class vars_t>
+    emtensor_t<data_t> compute_emtensor(
+        const vars_t<data_t> &vars,          //!< the value of the variables
+        const vars_t<Tensor<1, data_t>> &d1, //!< the value of the 1st derivs
+        const Tensor<2, data_t> &h_UU, //!< the inverse metric (raised indices)
+        const Tensor<3, data_t> &chris_ULL)
+        const; //!< the conformal christoffel symbol
+
+    //! The function which calculates the EM Tensor, given the vars and
+    //! derivatives, excluding the potential
+    template <class data_t, template <typename> class vars_t>
+    static void emtensor_excl_potential(
+        emtensor_t<data_t> &out,             //!< the em tensor output
+        const vars_t<data_t> &vars,          //!< the value of the variables
+        const vars_t<Tensor<1, data_t>> &d1, //!< the value of the first derivs
+        const Tensor<2, data_t> &h_UU, //!< the inverse metric (raised indices).
+        const Tensor<3, data_t>
+            &chris_ULL); //!< the conformal christoffel symbol
+
+    //! The function which adds in the RHS for the matter field vars,
+    //! including the potential
+    template <class data_t, template <typename> class vars_t,
+              template <typename> class diff2_vars_t,
+              template <typename> class rhs_vars_t>
+    void add_matter_rhs(
+        rhs_vars_t<data_t> &total_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; //!< the value of the advection terms
+
+    //! The function which calculates the RHS for the matter field vars
+    //! excluding the potential
+    template <class data_t, template <typename> class vars_t,
+              template <typename> class diff2_vars_t,
+              template <typename> class rhs_vars_t>
+    static void matter_rhs_excl_potential(
+        rhs_vars_t<data_t> &rhs, //!< the value of the RHS terms for the sf vars
+        const vars_t<data_t> &vars, //!< the values of all the variables
+        const vars_t<Tensor<1, data_t>> &d1, //!< the 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);
+};
+
+#include "ComplexScalarField.impl.hpp"
+
+#endif /* COMPLEXSCALARFIELD_HPP_ */

Source/Matter/ComplexScalarField.impl.hpp

@@ -0,0 +1,149 @@
+/* GRChombo
+ * Copyright 2012 The GRChombo collaboration.
+ * Please refer to LICENSE in GRChombo's root directory.
+ */
+
+#if !defined(COMPLEXSCALARFIELD_HPP_)
+#error "This file should only be included through ComplexScalarField.hpp"
+#endif
+
+#ifndef COMPLEXSCALARFIELD_IMPL_HPP_
+#define COMPLEXSCALARFIELD_IMPL_HPP_
+
+// Calculate the stress energy tensor elements
+template <class potential_t>
+template <class data_t, template <typename> class vars_t>
+emtensor_t<data_t> ComplexScalarField<potential_t>::compute_emtensor(
+    const vars_t<data_t> &vars, const vars_t<Tensor<1, data_t>> &d1,
+    const Tensor<2, data_t> &h_UU, const Tensor<3, data_t> &chris_ULL) const
+{
+    emtensor_t<data_t> out;
+
+    // call the function which computes the em tensor excluding the potential
+    emtensor_excl_potential(out, vars, d1, h_UU, chris_ULL);
+
+    // set the potential values
+    data_t V_of_phi = 0.0;
+    data_t dVdphi_Re = 0.0;
+    data_t dVdphi_Im = 0.0;
+
+    // compute potential and add constributions to EM Tensor
+    my_potential.compute_potential(V_of_phi, dVdphi_Re, dVdphi_Im, vars);
+
+    out.rho += V_of_phi;
+    out.S += -3.0 * V_of_phi;
+    FOR(i, j) { out.Sij[i][j] += -vars.h[i][j] * V_of_phi / vars.chi; }
+
+    return out;
+}
+
+// Calculate the stress energy tensor elements
+template <class potential_t>
+template <class data_t, template <typename> class vars_t>
+void ComplexScalarField<potential_t>::emtensor_excl_potential(
+    emtensor_t<data_t> &out, const vars_t<data_t> &vars,
+    const vars_t<Tensor<1, data_t>> &d1, const Tensor<2, data_t> &h_UU,
+    const Tensor<3, data_t> &chris_ULL)
+{
+    // Useful quantity Vt
+    data_t Vt = -vars.Pi_Re * vars.Pi_Re - vars.Pi_Im * vars.Pi_Im;
+    FOR2(i, j)
+    {
+        Vt += gamma_UU[i][j] * d1.phi_Re[i] * d1.phi_Re[j] +
+              gamma_UU[i][j] * d1.phi_Im[i] * d1.phi_Im[j];
+    }
+
+    // Calculate components of EM Tensor
+    // S_ij = T_ij
+    FOR(i, j)
+    {
+        out.Sij[i][j] = -0.5 * vars.h[i][j] * Vt / vars.chi +
+                        d1.phi_Re[i] * d1.phi_Re[j] +
+                        d1.phi_Im[i] * d1.phi_Im[j];
+    }
+
+    // S = Tr_S_ij
+    out.S = vars.chi * TensorAlgebra::compute_trace(out.Sij, h_UU);
+
+    // S_i (note lower index) = - n^a T_ai
+    FOR(i)
+    {
+        out.Si[i] = -d1.phi_Re[i] * vars.Pi_Re - d1.phi_Im[i] * vars.Pi_Im;
+    }
+
+    // rho = n^a n^b T_ab
+    out.rho = vars.Pi_Re * vars.Pi_Re + vars.Pi_Im * vars.Pi_Im + 0.5 * Vt;
+}
+
+// Adds in the RHS for the matter vars
+template <class potential_t>
+template <class data_t, template <typename> class vars_t,
+          template <typename> class diff2_vars_t,
+          template <typename> class rhs_vars_t>
+void ComplexScalarField<potential_t>::add_matter_rhs(
+    rhs_vars_t<data_t> &total_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
+{
+    // first get the non potential part of the rhs
+    // this may seem a bit long winded, but it makes the function
+    // work for more multiple fields
+
+    // call the function for the rhs excluding the potential
+    matter_rhs_excl_potential(total_rhs, vars, d1, d2, advec);
+
+    // set the potential values
+    data_t V_of_phi = 0.0;
+    data_t dVdphi_Re = 0.0;
+    data_t dVdphi_Im = 0.0;
+    my_potential.compute_potential(V_of_phi, dVdphi_Re, dVdphi_Im, vars);
+
+    // adjust RHS for the potential term
+    total_rhs.Pi_Re += -vars.lapse * dVdphi_Re;
+    total_rhs.Pi_Im += -vars.lapse * dVdphi_Im;
+}
+
+// the RHS excluding the potential terms
+template <class potential_t>
+template <class data_t, template <typename> class vars_t,
+          template <typename> class diff2_vars_t,
+          template <typename> class rhs_vars_t>
+void ComplexScalarField<potential_t>::matter_rhs_excl_potential(
+    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)
+{
+    using namespace TensorAlgebra;
+
+    const auto h_UU = compute_inverse_sym(vars.h);
+    const auto chris = compute_christoffel(d1.h, h_UU);
+
+    // evolution equations for scalar field and (minus) its conjugate momentum
+    rhs.phi_Re = vars.lapse * vars.Pi_Re + advec.phi_Re;
+    rhs.Pi_Re = vars.lapse * vars.K * vars.Pi_Re + advec.Pi_Re;
+    rhs.phi_Im = vars.lapse * vars.Pi_Im + advec.phi_Im;
+    rhs.Pi_Im = vars.lapse * vars.K * vars.Pi_Im + advec.Pi_Im;
+
+    FOR(i, j)
+    {
+        // includes non conformal parts of chris not included in chris_ULL
+        rhs.Pi_Re +=
+            h_UU[i][j] * (-0.5 * d1.chi[j] * vars.lapse * d1.phi_Re[i] +
+                          vars.chi * vars.lapse * d2.phi_Re[i][j] +
+                          vars.chi * d1.lapse[i] * d1.phi_Re[j]);
+        rhs.Pi_Im +=
+            h_UU[i][j] * (-0.5 * d1.chi[j] * vars.lapse * d1.phi_Im[i] +
+                          vars.chi * vars.lapse * d2.phi_Im[i][j] +
+                          vars.chi * d1.lapse[i] * d1.phi_Im[j]);
+        FOR(k)
+        {
+            rhs.Pi_Re += -vars.chi * vars.lapse * h_UU[i][j] *
+                         chris.ULL[k][i][j] * d1.phi_Re[k];
+            rhs.Pi_Im += -vars.chi * vars.lapse * h_UU[i][j] *
+                         chris.ULL[k][i][j] * d1.phi_Im[k];
+        }
+    }
+}
+
+#endif /* COMPLEXSCALARFIELD_IMPL_HPP_ */