Hands on example with Brill waves

This is a small tutorial to show how to change an already existing example. The prerequisites for this are

• You already compiled Chombo
• Able to run and submit code on your computer/cluster
• Know how to visualize data using visIT or other applications

Brill waves

You can find more background about Brill waves in Miguel Alcubierre's book "Introduction to 3+1 Numerical Relativity" on page 395 in section 10.4.2. We will mostly follow their guidance, though, for simplicity, instead of solving Psi, we solve for chi. The result is not exactly the Brill waves described in the book, though that is not the aim of this exercise.

Coding

First, we prepare to set some code up that we can manipulate.

• Get Fresh clone

  git clone https://github.com/GRChombo/GRChombo.git
  cd GRChombo/Examples/KerrBH
  git checkout 52c5e8f4e43506c8d482e789f6d2b2ac77ba114a
  

• Create new class for Initial data "BrillWaveData"

    vi BrillWave.hpp
  

  The class should contain

  • A member double m_dx
  • a templated function over data_t called "compute" that takes as an argument Cell<data_t> current_cell
  • A constructor that takes an argument a_dx and writes it on m_dx

  An example for the class described above:

    #ifndef BRILLWAVE_HPP_
    #define BRILLWAVE_HPP_
  
    #include "Cell.hpp"
    #include "Coordinates.hpp"
    #include "Tensor.hpp"
    #include "TensorAlgebra.hpp"
    #include "UserVariables.hpp" //This files needs NUM_VARS - total no. components
    #include "VarsTools.hpp"
    #include "simd.hpp"
  
    //! Class which calculates Axial symmetric GWs 
    class BrillWaveData
    {
      public:
  
        //! The constructor
        BrillWaveData( double a_dx)
        : m_dx(a_dx)
        {
        }
  
        //! Function to compute the value of all the initial vars on the grid
        template <class data_t> void compute(Cell<data_t> current_cell) const
        {
  
            BSSNVars::VarsWithGauge<data_t> vars;
            VarsTools::assign(vars, 0.); // Set only the non-zero components explicitly below
  
            // store the vars
            current_cell.store_vars(vars);
        }
  
  
      protected:
            double m_dx;
    };
  
    #endif /* BRILLWAVE_HPP_ */
  

  This is a simple initial data class that sets all values to zero.

• Include this file in KerrBHLevel.cpp. Go to line 20 and replace

  #include "KerrBH.hpp"

with

#include "BrillWave.hpp"

• Next, change line 41 from

    make_compute_pack(SetValue(0.), KerrBH(m_p.kerr_params, m_dx));
  

  to

    make_compute_pack(SetValue(0.), BrillWaveData(m_dx))
  

• Compile the code and check if everything compiles. (You can also run the code, but it will crash after the first step with all quantites being zero.)

---

Now we can get started setting up the exact shape of the initial data we plan to implement.

• First we need coordinates that we can call. Add in the compute function

    template <class data_t> void compute(Cell<data_t> current_cell) const {
     ... 
    }
  

  an initialisation of the coordinate class

    Coordinates<data_t> coords(current_cell, m_dx);
  

• Now we can call define all coordinates which we will need soon

    data_t x = coords.x;
    double y = coords.y;
    double z = coords.z;
    data_t rho2 = (x*x + y*y);
  

  Note: We treat x differently on purpose (This is important for code optimisation using vector operations).

• We can now define variables needed for the calculation

    data_t q; 
    data_t Psi4; // Psi to the power of four 
    Tensor<2,data_t> gamma; // spatial metric
  

• We initialise the metric

    FOR2(i,j){
        gamma[i][j] = 0;
    }
    FOR1(i){
        gamma[i][i] = 1;
    }
  

  Note: In the code we defined FOR loop macros for(int i = 0; i < 3; i++)

• We choose a form for q ( This is a free choice ), we propose this

     q =  0.01 * rho2 * exp(-rho2-z*z);
  

• We would need to calculate Psi4 later on inside the code, for the moment we set it to one

    Psi4 = 1;  
  

• We write the cartesian form of the metric

    gamma[0][0] = ( exp(2.*q)*x*x + y*y )*Psi4/(x*x+y*y);
    gamma[1][1] = ( x*x + exp(2.*q)*y*y )*Psi4/(x*x+y*y);
    gamma[0][1] = ( exp(2.*q) - 1. )*x*y*Psi4/(x*x+y*y);
    gamma[1][0] = gamma[0][1];
    gamma[2][2] = exp(2.*q) * Psi4;
  

• We have to decompose the ADM spatial metric gamma into BSSN components

     // BSSN decomposition
    data_t detgamma = TensorAlgebra::compute_determinant(gamma);
  
    vars.chi = 1/detgamma;
  
    FOR2(i,j){
        vars.h[i][j] = pow(vars.chi,3)*gamma[i][j];
    }
  

• Set the gauge variables

    vars.lapse = 1.0; 
  

Fix tagging criterion

• For this we manipulate the tagging function in KerrBHLevel.cpp

    void KerrBHLevel::computeTaggingCriterion(FArrayBox &tagging_criterion,
                                      const FArrayBox &current_state)
    {
    ...
    }
  

• We change the

        BoxLoops::loop(ChiTaggingCriterion(m_dx), current_state, tagging_criterion);
  

  to new fixed tagging criterion

        const double L = 2.5;
        const std::array<double, CH_SPACEDIM> center = {0,0,0};
        BoxLoops::loop(FixedGridsTaggingCriterion(m_dx,m_level,L,center), current_state, tagging_criterion);
  

• On top of the file include

    #include "FixedGridsTaggingCriterion.hpp"       
  

Fix the params file

• Open params.txt and change the following parameters. The box is a bit too large initially, we want to change the physical size

    L_full = 12
  

• Change the boundaries to be reflective, to use the symmetry of the problem

    lo_boundary = 2 2 2 
  

• There is also the need to set

    extraction_radii = 0.1 0.1 0.1
  

Note: This is due to issue that might be fixed at a later date

• Where we set to restrict the maximum refinement level to one

    max_level = 1
  

Set up relaxation routine

We want to modify the underlying equations of motion in KerrBHLevel.cpp. There we want to change the function

    void KerrBHLevel::specificEvalRHS(GRLevelData &a_soln, GRLevelData &a_rhs,
                              const double a_time)
    {
    ...
    }

where we change the evolution equations for CCZ4

    // Enforce the trace free A_ij condition and positive chi and alpha
    BoxLoops::loop(make_compute_pack(TraceARemoval(), PositiveChiAndAlpha()),
               a_soln, a_soln, INCLUDE_GHOST_CELLS);

    // Calculate CCZ4 right hand side
    BoxLoops::loop(CCZ4(m_p.ccz4_params, m_dx, m_p.sigma, m_p.formulation),
               a_soln, a_rhs, EXCLUDE_GHOST_CELLS);

we add the if statement

    if(a_time > 2.50 ) {
           <CCZ4 EVOLUTION>
    }else {
    double relaxation_speed = 0.01;
    BoxLoops::loop(ChiRelaxation( m_dx, relaxation_speed ),
               a_soln, a_rhs, EXCLUDE_GHOST_CELLS);

    }

that will use the ChiRelaxation until time 2.5 with relaxation_speed of 0.01. To include the matter folder in the compilation of the code, add in the GNUMakefile

    src_dirs := $(GRCHOMBO_SOURCE)/utils \
                $(GRCHOMBO_SOURCE)/simd  \
                $(GRCHOMBO_SOURCE)/CCZ4  \
                $(GRCHOMBO_SOURCE)/BoxUtils  \
    +           $(GRCHOMBO_SOURCE)/Matter  \
                $(GRCHOMBO_SOURCE)/GRChomboCore  \
                $(GRCHOMBO_SOURCE)/TaggingCriteria  \
                $(GRCHOMBO_SOURCE)/AMRInterpolator  \
                $(GRCHOMBO_SOURCE)/InitialConditions/BlackHoles

and add the appropriate include

    #include "ChiRelaxation.hpp"

These steps should give a spacetime consistent with GR, but the solver will take very long to provide a good solution (i.e., small hamiltonian constraint). It might need some further playing around to get a good quality spacetime. However, this guide's idea was to introduce how to change different parts of GRChombo, so we didn't spend much more time refining it. Feel free to add comments if you found suitable parameters for relaxation time/relaxation speed/max_level/maximum amplitude in q that will give a consistent simulation.