Adding Cosmology Example

.gitignore

@@ -10,6 +10,7 @@ ipo_out.optrpt
 # Scheduler files
 jobscript*
 *.pbs
+*.sh
 
 # Example output
 Weyl_integral*
@@ -36,4 +37,5 @@ log*
 .*.sw[n-p]
 *.dSYM
 *.patch
-*.log
+*.log
+*.err

Examples/ScalarFieldCosmo/ConstraintsExtraction.hpp

@@ -0,0 +1,116 @@
+/* GRChombo
+ * Copyright 2012 The GRChombo collaboration.
+ * Please refer to LICENSE in GRChombo's root directory.
+ */
+
+#ifndef CONSTRAINTSEXTRACTION_HPP_
+#define CONSTRAINTSEXTRACTION_HPP_
+
+#include "AMRInterpolator.hpp"
+#include "InterpolationQuery.hpp"
+#include "Lagrange.hpp"
+#include "SimulationParametersBase.hpp"
+#include "SmallDataIO.hpp"
+#include "SphericalHarmonics.hpp"
+#include "UserVariables.hpp"
+#include <fstream>
+#include <iomanip>
+#include <iostream>
+#include <sstream>
+#include <string>
+
+//!  The class allows extraction of constraints (Ham and Mom) along x-axis
+
+class ConstraintsExtraction
+{
+  private:
+    //! Params for extraction
+    const int m_comp1;
+    const int m_comp2;
+    const int m_num_points;
+    const double m_L;
+    const std::array<double, CH_SPACEDIM>
+        m_origin; // Origin of an extraction line
+    const double m_dt;
+    const double m_time;
+
+  public:
+    //! The constructor
+    ConstraintsExtraction(int a_comp1, int a_comp2, int a_num_points,
+                          double a_L, std::array<double, CH_SPACEDIM> a_origin,
+                          double a_dt, double a_time)
+        : m_comp1(a_comp1), m_comp2(a_comp2), m_num_points(a_num_points),
+          m_origin(a_origin), m_L(a_L), m_dt(a_dt), m_time(a_time)
+    {
+    }
+
+    //! Destructor
+    ~ConstraintsExtraction() {}
+
+    //! Execute the query
+    void execute_query(AMRInterpolator<Lagrange<2>> *a_interpolator,
+                       std::string a_file_prefix) const
+    {
+        CH_TIME("CustomExtraction::execute_query");
+        if (a_interpolator == nullptr)
+        {
+            MayDay::Error("Interpolator has not been initialised.");
+        }
+        std::vector<double> interp_ham_data(m_num_points);
+        std::vector<double> interp_mom_data(m_num_points);
+        std::vector<double> interp_x(m_num_points);
+        std::vector<double> interp_y(m_num_points);
+        std::vector<double> interp_z(m_num_points);
+
+        // Work out the coordinates
+        // go out along x-axis from 0 to L
+        for (int idx = 0; idx < m_num_points; ++idx)
+        {
+            interp_x[idx] = (double(idx) / double(m_num_points) * m_L);
+            interp_y[idx] = m_origin[1];
+            interp_z[idx] = m_origin[2];
+        }
+
+        // set up the query
+        InterpolationQuery query_ham(m_num_points);
+        query_ham.setCoords(0, interp_x.data())
+            .setCoords(1, interp_y.data())
+            .setCoords(2, interp_z.data())
+            .addComp(m_comp1, interp_ham_data.data(), Derivative::LOCAL,
+                     VariableType::diagnostic); // evolution or diagnostic
+        InterpolationQuery query_mom(m_num_points);
+        query_mom.setCoords(0, interp_x.data())
+            .setCoords(1, interp_y.data())
+            .setCoords(2, interp_z.data())
+            .addComp(m_comp2, interp_mom_data.data(), Derivative::LOCAL,
+                     VariableType::diagnostic); // evolution or diagnostic
+
+        // submit the query
+        a_interpolator->interp(query_ham);
+        a_interpolator->interp(query_mom);
+
+        // now write out
+        bool first_step = (m_time == 0.0);
+        double restart_time = 0.0;
+        SmallDataIO output_file(a_file_prefix, m_dt, m_time, restart_time,
+                                SmallDataIO::APPEND, first_step);
+
+        std::vector<std::string> header_line(2);
+        if (first_step)
+        {
+            header_line[0] = "Ham";
+            header_line[1] = "Mom";
+            output_file.write_header_line(header_line, "x");
+        }
+
+        for (int idx = 0; idx < m_num_points; ++idx)
+        {
+            std::vector<double> data(2);
+            data[0] = interp_ham_data[idx];
+            data[1] = interp_mom_data[idx];
+            output_file.write_data_line(data, interp_x[idx]);
+        }
+    }
+};
+
+#endif /* CONSTRAINTSEXTRACTION_HPP_ */

Examples/ScalarFieldCosmo/CosmoDiagnostics.hpp

@@ -0,0 +1,65 @@
+/* GRChombo
+ * Copyright 2012 The GRChombo collaboration.
+ * Please refer to LICENSE in GRChombo's root directory.
+ */
+
+#ifndef COSMODIAGNOSTICS_HPP_
+#define COSMODIAGNOSTICS_HPP_
+
+#include "CCZ4Geometry.hpp"
+#include "CCZ4Vars.hpp"
+#include "Cell.hpp"
+#include "Coordinates.hpp"
+#include "DiagnosticVariables.hpp"
+#include "FourthOrderDerivatives.hpp"
+#include "GRInterval.hpp"
+#include "NewMatterConstraints.hpp"
+#include "Tensor.hpp"
+#include "simd.hpp"
+#include <array>
+
+//! Calculates all relevant variables, which
+//! are stored as diagnostics
+
+template <class matter_t>
+class CosmoDiagnostics // public MatterConstraints<matter_t>
+{
+  public:
+    // Inherit the variable definitions from CCZ4 + matter_t
+    template <class data_t>
+    using CCZ4Vars = typename CCZ4::template Vars<data_t>;
+
+    template <class data_t>
+    using MatterVars = typename matter_t::template Vars<data_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);
+        }
+    };
+
+    //! Constructor of class CosmoDiagnostics
+    CosmoDiagnostics(const matter_t a_matter, double dx, double G_Newton)
+        : m_matter(a_matter), m_deriv(dx)
+    {
+    }
+
+    //! The compute member which calculates the diagnostics at each point in the
+    //! box
+    template <class data_t> void compute(Cell<data_t> current_cell) const;
+
+  protected:
+    matter_t m_matter; //!< The matter object, e.g. a scalar field
+    const FourthOrderDerivatives m_deriv;
+};
+
+#include "CosmoDiagnostics.impl.hpp"
+
+#endif /* COSMODIAGNOSTICS_HPP_ */

Examples/ScalarFieldCosmo/CosmoDiagnostics.impl.hpp

@@ -0,0 +1,54 @@
+/* GRChombo
+ * Copyright 2012 The GRChombo collaboration.
+ * Please refer to LICENSE in GRChombo's root directory.
+ */
+
+#if !defined(COSMODIAGNOSTICS_HPP_)
+#error "This file should only be included through CosmoDiagnostics.hpp"
+#endif
+
+#ifndef COSMODIAGNOSTICS_IMPL_HPP_
+#define COSMODIAGNOSTICS_IMPL_HPP_
+#include "DimensionDefinitions.hpp"
+
+template <class matter_t>
+template <class data_t>
+void CosmoDiagnostics<matter_t>::compute(Cell<data_t> current_cell) const
+{
+    // Load local vars and calculate derivs
+    const auto vars = current_cell.template load_vars<Vars>();
+    const auto d1 = m_deriv.template diff1<Vars>(current_cell);
+    const auto d2 = m_deriv.template diff2<Vars>(current_cell);
+
+    // Inverse metric and Christoffel symbol
+    const auto h_UU = TensorAlgebra::compute_inverse_sym(vars.h);
+    const auto chris = TensorAlgebra::compute_christoffel(d1.h, h_UU);
+
+    // Define quantities
+    data_t rho;
+    data_t sqrt_gamma;
+    data_t S;
+    data_t rho_scaled;
+    data_t S_scaled;
+    data_t K_scaled;
+    data_t a;
+
+    // Energy Momentum Tensor
+    const auto emtensor = m_matter.compute_emtensor(vars, d1, h_UU, chris.ULL);
+
+    sqrt_gamma = pow(vars.chi, -3. / 2.);
+    rho = emtensor.rho;
+    S = emtensor.S;
+    K_scaled = vars.K / pow(vars.chi, 3. / 2.);
+    rho_scaled = emtensor.rho / pow(vars.chi, 3. / 2.);
+    S_scaled = emtensor.S / pow(vars.chi, 3. / 2.);
+
+    // Write the diagnostics into the output vector of cells
+    current_cell.store_vars(sqrt_gamma, c_sqrt_gamma);
+    current_cell.store_vars(rho, c_rho);
+    current_cell.store_vars(rho_scaled, c_rho_scaled);
+    current_cell.store_vars(S_scaled, c_S_scaled);
+    current_cell.store_vars(K_scaled, c_K_scaled);
+}
+
+#endif /* COSMODIAGNOSTICS_IMPL_HPP_ */

Examples/ScalarFieldCosmo/CustomExtraction.hpp

@@ -0,0 +1,121 @@
+/* GRChombo
+ * Copyright 2012 The GRChombo collaboration.
+ * Please refer to LICENSE in GRChombo's root directory.
+ */
+
+#ifndef CUSTOMEXTRACTION_HPP_
+#define CUSTOMEXTRACTION_HPP_
+
+#include "AMRInterpolator.hpp"
+#include "InterpolationQuery.hpp"
+#include "Lagrange.hpp"
+#include "SimulationParametersBase.hpp"
+#include "SmallDataIO.hpp"
+#include "SphericalHarmonics.hpp"
+#include "UserVariables.hpp"
+#include <fstream>
+#include <iomanip>
+#include <iostream>
+#include <sstream>
+#include <string>
+
+// Function to convert double to string with precision
+std::string to_string_with_precision(double value, int precision)
+{
+    std::ostringstream out;
+    out << std::fixed << std::setprecision(precision) << value;
+    return out.str();
+}
+
+//!  The class allows extraction of the values of components at
+//!  specified points
+
+//!  In this example, we do a lineout data extraction along x-axis
+
+class CustomExtraction
+{
+  private:
+    //! Params for extraction
+    const int m_comp;
+    const int m_num_points;
+    const double m_L;
+    const std::array<double, CH_SPACEDIM>
+        m_origin; // Origin of an extraction line
+    const double m_dt;
+    const double m_time;
+
+  public:
+    //! The constructor
+    CustomExtraction(int a_comp, int a_num_points, double a_L,
+                     std::array<double, CH_SPACEDIM> a_origin, double a_dt,
+                     double a_time)
+        : m_comp(a_comp), m_num_points(a_num_points), m_origin(a_origin),
+          m_L(a_L), m_dt(a_dt), m_time(a_time)
+    {
+    }
+
+    //! Destructor
+    ~CustomExtraction() {}
+
+    //! Execute the query
+    void execute_query(AMRInterpolator<Lagrange<4>> *a_interpolator,
+                       std::string a_file_prefix) const
+    {
+        CH_TIME("CustomExtraction::execute_query");
+        if (a_interpolator == nullptr)
+        {
+            MayDay::Error("Interpolator has not been initialised.");
+        }
+        std::vector<double> interp_var_data(m_num_points);
+        std::vector<double> interp_x(m_num_points);
+        std::vector<double> interp_y(m_num_points);
+        std::vector<double> interp_z(m_num_points);
+
+        // Work out the coordinates
+        // go out along x-axis from m_origin to L
+        for (int idx = 0; idx < m_num_points; ++idx)
+        {
+            interp_x[idx] =
+                m_origin[0] + (double(idx) / double(m_num_points) * m_L);
+            interp_y[idx] = m_origin[1];
+            interp_z[idx] = m_origin[2];
+        }
+
+        // set up the query
+        InterpolationQuery query(m_num_points);
+        query.setCoords(0, interp_x.data())
+            .setCoords(1, interp_y.data())
+            .setCoords(2, interp_z.data())
+            .addComp(m_comp, interp_var_data.data(), Derivative::LOCAL,
+                     VariableType::diagnostic); // evolution or diagnostic
+
+        // submit the query
+        a_interpolator->interp(query);
+
+        // now write out
+        bool first_step = (m_time == 0.0);
+        double restart_time = 0.0;
+        SmallDataIO output_file(a_file_prefix, m_dt, m_time, restart_time,
+                                SmallDataIO::APPEND, first_step);
+
+        std::vector<std::string> header_line(m_num_points);
+
+        if (first_step)
+        {
+            for (int i = 0; i < m_num_points; ++i)
+            {
+                // header_line[i] =
+                //     "p" + std::to_string(i + 1) + "(" +
+                //     to_string_with_precision(interp_x[i], 2) + "," +
+                //     to_string_with_precision(m_origin[1], 2) + "," +
+                //     to_string_with_precision(m_origin[2], 2) + ")";
+                header_line[i] =
+                    "x = " + to_string_with_precision(interp_x[i], 2);
+            }
+            output_file.write_header_line(header_line);
+        }
+        output_file.write_time_data_line(interp_var_data);
+    }
+};
+
+#endif /* CUSTOMEXTRACTION_HPP_ */