Bimetric 3+1 Toolkit implements the standard 3+1 evolution for spherically symmetric bimetric and GR spacetimes.

The project goal is the construction of a computational toolkit for the evolution of 3+1 equations of modified theories of gravity in spherical symmetry, at the same time resembling the functionality of the EinsteinToolkit to ease transition to the full scale projects. The current implementation handles bimetric and GR spacetimes.

Note that the toolkit can handle any system of differential equations. For a pedagogical example how to use the toolkit, see Gowdy spacetimes solver.

Version:: 0.28, inception 2018-04-23, last modified 2018-09-18 09:05

Author:: Mikica Kocic

Source code:: bimetric-ss.zip
A3 overview: Bimetric equations in standard 3+1 form in spherical symmetry

Main file:: bim-solver.cpp, see also Module and File Overview

Note:: The 3+1 decomposition is based on [7]. The numerical algorithms are based on [9], [3], and [12].

Example simulation (dust collapse testbed)

The simulation is compliant with the results of [8].

Main Features

Formalism

Bimetric equations in standard 3+1 form (with the evolution of p and the explicit bimetric lapse ratio).

Matter equations for the perfect fluid in conservative form.

Equations regularized for spherical symmetry.

Gauge setup

Lapse:
- Maximal slicing (boundary value problem at each slice, see maximalSlice.h).
- Bona-Massó slicing condition and the K-driver (evolution).
- Algebraic slicing (normal coordinates, and (1+erf)^{-2}).

Shift:
- Planned: Minimal distortion and Gamma-driver.

Boundary conditions

Imposed parity conditions for local flattness at r = 0 (on the inner boundary).

Extrapolation (linear 2nd, or 4th order) or Sommerfeld outgoing wave radiative condition on the outer boundary.

Spatial discretization

2nd, 4th, or 6th order centered differences on a staggered grid (see finiteDifferences.h).

Planned: Upwind differences on shift advection terms.

Temporal discretization

Method of Lines:
- Runge-Kutta: 2, 3, and 4 steps
- Iterated Crank-Nicolson (ICN): 2 and 3 steps
- Averaged ICN: 2 and 3 steps
- Generic N-step algorithm with arbitrary coefficients (see methodOfLines.h)

Kreiss-Oliger dissipation (2nd and 4th order) in all of the evolution equations.

Courant-Friedrichs-Lewy factor (CFL) 0.5 as default.

Grid-driver code

Uniform spatial grid (see gridDriver.h).

Numerics

Implemented classes: Matrix, Vector, BandLUDecomposition (band-diagonal matrix LU decomposition with Crout factorization), CubicSpline (normal cubic spline interpolation), arbirary FD extrapolation, and arbitrary FD derivatives.

64-bit and 128-bit floating point

C++ code with OpenMP and MPI support adapted for high-performance computing.

Planned: Transition to Cactus

Horizons

Apparent horizon finder.

Initial Data

Minkowski GR (opt. with a gauge wave)

Bimetric Minkowski (opt. with a gauge wave)

GR collisionless matter (dust) with Gaussian profile (as the main testbed).

Bimetric Gaussian dust with a "polytropic" conformally flat g and f.

Implementation details

Main Class

BimetricEvolve: Encapsulates a 3+1 evolution solver for bimetric spacetimes.

Modules across the source files

     +-------------------------------------------------------+    +-----------------+
     |                  Bimetric Evolution                   |    |  Bimetric Model |
     +-------------------------------------------------------+    +-----------------+
     | bim-solver.cpp   eomEvolution.h    eomGauge.h         |    | bimetricModel.h |
     | maximalSlice.h   eomMatter.h       eomLapseRatios.h   |    +-----------------+
     |                  eomSources.h      eomMiscEquations.h |
     |                  eomConstraints.h                     |
     +-------------------------------------------------------+
   
  +-----------------+   +-------------------+   +----------------+   +---------------+
  |   Grid Driver   |   |   Initial Data    |   | MoL Integrator |   |  Output Data  |
  +-----------------+   +-------------------+   +----------------+   +---------------+
  | gridDriver.h    |   | gridInitialData.h |   |  integrator.h  |   | gridOutput.h  |
  | gridFunctions.h |   +-------------------+   +----------------+   +---------------+
  +-----------------+
         +---------------------+      +--------------+      +-----------------+
         |  Numerical Methods  |      |  Data Types  |      |    Utilities    |
         +---------------------+      +--------------+      +-----------------+
         | finiteDifferences.h |      | matrix.h     |      | trackUsedTime.h |
         | methodOfLines.h     |      | dataTypes.h  |      | mpiWorld.h      |
         | cubicSpline.h       |      +--------------+      | mpiDummyWorld.h |
         | bandSol.h           |                            | paramsHolder.h  |
         +---------------------+                            | slog.h          |
                                                            +-----------------+

Data Flow

  +---------------------------------+              +----------------------------+
  | Mathematica Notebook            |------------>>| Initial Data & Parameters  |
  | bimetric-ss-in-3+1 (cpp solver) |  Solve the   |  input.dat  (initial data) |
  +---------------------------------+  constraint  |  config.ini (parameters)   |
            |   Export the             equations   +----------------------------+
            |   C++ code                                         |
            V                                                    V
  +---------------------------------+  eom-std.h   +----------------------------+
  | eomEvolution.h eomMatter.h      |------------>>|       bim-solver.cpp       |
  | eomSources.h   eomConstraints.h |              +----------------------------+
  | eomGauge.h     eomLapseRatios.h |                            |
  | eomMiscEquations.h              |                            V
  +---------------------------------+              +----------------------------+
                                                   |         output.dat         |
     See also:  maximalSlice.h                     +----------------------------+
      
     Numerics:  methodOfLines.h, finiteDifferences.h,
                dataTypes.h, matrix.h, bandSol.h, cubicSpline.h

Regularization of spherically symmetric evolution codes

The equations are regularized according to: [1], [2], [4], [5], [6], [10], and [11].

A nice overview of the parity conditions at r = 0 can be found in [11].

Also, parity conditions for components of vectors and tensors can be found in Table 1 in [4].

Finite Difference Approximation

At the core of finite difference approximation is a discretization of the spacetime, or a numerical grid (shown below). Then, a function f(t,r) is represented by values at a discrete set of points (m,n).

                p in M          <-->     (t,r)      <-->        (m,n)
         (point in a manifold)        (in a chart)        (discretized point)
                                           |                      |
                                           V                      V
                                         f(t,r)                 f(m,n)

The points (m,n) are called the grid points or nodes. The grid points are chosen to reside at the center of grid cells (the other approach would be to associate the grid points with the vertices of the grid cells). The data structure GridPoint is a wrapper for the set of fields {f,...} that arekept in the memory for the point (m,n).The grid consists of the interior points, and the virtual points at the boundary. The virtual points are used to impose the boundary conditions. A cell-centered uniform grid with nLen grid points in each grid row is shown below (see also [3], p.192).

         left virtual                grid cell        grid point     right virtual
            point(s)                 .---^---.            |             point(s)
      .        |               n-1  '    n    '  n+1      |                |
      :   +----|----+-------+-------+---------+-------+---|---+-------+----|----+
     m+1  |    O    |       |       | (m+1,n) |       |   X   |       |    O    |
          +---------+-------+-------.=========.-------+-------+-------+---------+
      m   |         |       |(m,n-1)|  (m,n)  |(m,n+1)|       |       |         |
          +---------+-------+-------'========='-------+-------+-------+---------+
     m-1  |         |       |       |         |       |       |       |         |
      .   +---------'-------+-------+---------+-------+-------+-------'---------+
      :        0    |   1       2      . . .                    nLen  |  nLen+1
                    ^                                                 ^
                   r_min                                            r_max

The cell-centered uniform grid implies:

     r = r_min + ( n - 1/2 ) delta_r
     delta_r = ( r_max - r_min ) / nLen

if nGhost = 1, n ranges between 1 and nLen; the boundary virtual points are grid[m][0] and grid[m][nLen+1]

m ranges between 1 and mLen, and wraps around; the slices grid[mLen+1]... are used for intermediate integration steps.

In parallel simulations, the grid is split across processes using MPI, allowing storage of large grid functions and parallel work on each part.

Compiling the code

The code can use OpenMP for almost all parallel loops on the spatial grid. These loops are indicated by OMP_parallel_for, which expands to:

     _Pragma("omp parallel for")   i.e., into: #pragma omp parallel for

To compile with enabled parallelism, use the -fopenmp switch:

     g++ -Wall -m64 -std=c++14 -O3 -fopenmp bim-solver.cpp -o bim-solver

The code can also be compiled to be executed on a cluster with the MPI support. To compile, use the mpic++ compiler with the -D_USEMPI switch:

     mpic++  -D_USEMPI  -Wall -m64 ...

Then use mpiexec to run the resulting executable on a cluster.

Warning:: The maximal slicing is not compliant with MPI since the boundary value problem for the slicing differential equation requires access to the whole grid!

References

[1]

Miguel Alcubierre and Jose A. Gonzalez. Regularization of spherically symmetric evolution codes in numerical relativity. Comput. Phys. Commun., 167:76–84, 2005.

[2]

Miguel Alcubierre and Martha D. Mendez. Formulations of the 3+1 evolution equations in curvilinear coordinates. Gen. Rel. Grav., 43:2769–2806, 2011.

[3]

Baumgarte and Shapiro. Numerical Relativity: Solving Einstein's Equations on the Computer. Cambridge University Press, 2010.

[4]

Baumgarte, Montero, Cordero-Carrión, and Müller. Numerical relativity in spherical polar coordinates: Evolution calculations with the BSSN formulation. Phys. Rev. D, 87:044026, Feb 2013.

[5]

Isabel Cordero-Carrión and Pablo Cerda-Duran. Partially implicit Runge-Kutta methods for wave-like equations in spherical-type coordinates. 2012.

[6]

Mewes et al. Numerical relativity in spherical coordinates with the Einstein Toolkit. Phys. Rev., D97(8):084059, 2018.

[7]

Mikica Kocic. Geometric mean of bimetric spacetimes. 2018.

[8]

Nakamura, T. and Maeda, M and Miyama, S. and Sasaki, M. General Relativistic Collapse of an Axially Symmetric Star. I The Formulation and the Initial Value Equations. Progress of Theoretical Physics, 63(4):1229–1244, 1980.

[9]

W.H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery. Numerical Recipes 3rd Edition: The Art of Scientific Computing. Cambridge University Press, 2007.

[10]

Ian Ruchlin, Zachariah B. Etienne, and Thomas W. Baumgarte. SENR/NRPy+: Numerical Relativity in Singular Curvilinear Coordinate Systems. Phys. Rev., D97(6):064036, 2018.

[11]

Milton Ruiz, Miguel Alcubierre, and Dario Nunez. Regularization of spherical and axisymmetric evolution codes in numerical relativity. Gen. Rel. Grav., 40:159–182, 2008.

[12]

Masaru Shibata. Numerical Relativity. 100 Years Of General Relativity. World Scientific Publishing Company, 2015.

Copyright:: GNU General Public License (GPLv3).