Apparent horizon finder

Here we give some brief instructions on using the Apparent Horizon finder feature of the code. More detailed (but potentially out of date) instructions are given in the Further Resources section below.

Building with the AH Finder

Note that using this feature requires modifications to the build system, in particular to use the PETSc library.

PETSc setup

If you are very lucky, the module file on your cluster may just work automatically, so you can simply do

module load <YourPETScModule>

If this doesn't work, check to see if the PETSC_DIR and PETSC_ARCH environment modules are set by it by either checking

module show <YourPETScModule>

or trying

echo $PETSC_DIR

and similarly for $PETSC_ARCH. Note that depending on how PETSc has been build, this last variable might need to be the empty string anyway.

If these variables aren't set by the module, you can set them appropriately in your shell initialization scripts (e.g. .bashrc), for example

export PETSC_DIR=/path/to/petsc-3.14.4
export PETSC_ARCH=linux-x86-rome

We will assume these variables are set correctly below.

If you are unlucky, you may need to build PETSc yourself using the instructions on their website. A page of miscellaneous tips on how to install/use this for machines we have used is given at Use and install PETSc.

Make.defs.local modifications

You may need to make the following three changes to your Chombo Make.defs.local file (the one mentioned in the instructions here).

1. Add the PETSc include paths so the compiler knows where to look for the PETSc header files. This can be done by adding the following to the XTRACPPFLAGS variable:

   XTRACPPFLAGS = -I${PETSC_DIR}/include -I${PETSC_DIR}/${PETSC_ARCH}/include

   Note that this might not be necessary if your PETSc module appends to the CPATH environment variable.
2. Link with the PETSc libraries by adding the following to the syslibflags variable:

    syslibflags = <other library flags> -L${PETSC_DIR}/${PETSC_ARCH}/lib -lpetsc

   Note that the -L${PETSC_DIR}/${PETSC_ARCH}/lib bit might not be necessary if your PETSc module suitably updates the LIBRARY_PATH environment variable.
3. Enable the AH Finder code path by defining the USE_AHFINDER macro. This can be done by adding the following to the XTRACPPFLAGS variable

   XTRACPPFLAGS = <other preprocessor flags> -DUSE_AHFINDER

   so after steps 1 and 3, it will look something like

   XTRACPPFLAGS = -I${PETSC_DIR}/include -I${PETSC_DIR}/${PETSC_ARCH}/include -DUSE_AHFINDER

   assuming there aren't any other flags passed to this variable.

Then make realclean before making Chombo again as usual per the instructions here. We wish you the best of luck.

Updates to the Example files

You will then need to set up the AHFinder in the Main_MyExample.cpp file. Make sure the AMR class used is BHAMR not GRAMR and that an AMRInterpolator is being set up. This should be the case if you are working from recent versions of the code, but if not, see the most recent Examples (e.g. KerrBH) to see how this is done there. You may also need to #include "BHAMR.hpp" in the ExampleLevel.hpp file and dynamically cast the m_gr_amr object to it. Again, see the KerrBH files as an example.q

You also need to add a block of code to initialise the AHFinder:

#ifdef USE_AHFINDER
    if (sim_params.AH_activate)
    {
        AHSurfaceGeometry sphere(sim_params.kerr_params.center);
        bh_amr.m_ah_finder.add_ah(sphere, sim_params.AH_initial_guess,
                                  sim_params.AH_params);
    }
#endif

For the use of $\chi$ contours to initialise, see the KerrBH example Main_KerrBH.cpp file.

To actually call the AHFinder, we need to add code where it should run, which is usually in the specificPostTimestep() function, on a specific level.

Therefore in the ExampleLevel.cpp file, in the specificPostTimestep() function, we add:

#ifdef USE_AHFINDER
    if (m_p.AH_activate && m_level == m_p.AH_params.level_to_run)
        m_bh_amr.m_ah_finder.solve(m_dt, m_time, m_restart_time);
#endif

Note that the solve function automatically manages multiple surfaces if they are specified initially, and as mergers (again see the BinaryBH example).

Updates to params.txt

Look at the KerrBH and BinaryBH example for the list of AH parameters that need to be added. The naming should be mostly self explanatory - note that many are commented out and the default value is given. More information on their usage can be found in AHParams.hpp file in the ApparentHorizonFinder Source folder.

The majority of the AHFinder params are loaded in the SimulationParametersBase class, so assuming you inherit from this you will not need to explicitly load them. However, since it is example specific, you will need to add an initial guess for the horizons in your SimulationParameters class - for example in the KerrBH case we add:

#ifdef USE_AHFINDER
        pp.load("AH_initial_guess", AH_initial_guess, 0.5 * kerr_params.mass);
#endif

and

#ifdef USE_AHFINDER
    double AH_initial_guess;
#endif

to specify the initial guess for the radius of the AH surface, as it is a single surface. In the BinaryBH example you can see how to add multiple surfaces where there are two BHs.

Processing the outputs

The AHFinder outputs some useful commentary into the pout.* files, along the lines of:

Setting Original Guess to f=0.244239
Solver converged. Horizon FOUND.
mass = 0.500228
spin = 1.55618e-05
Setting Original Guess to f=0.244239
Solver converged. Horizon FOUND.
mass = 0.500228
spin = 1.55618e-05
BHs #0 and #1 at distance = 12.2138 > minimum distance = 1.95392. Skipping solve for merger...

Data is also output in a ‘coords’ file for each AH and for each step, containing the coordinates of the AH surface, and a ‘stats’ file for each AH, containing the physical information about it (e.g. its area and spin) for all timesteps. These data files can simply be plotted using gnuplot or python as usual, e.g. for the KerrBH example, one can use

import numpy as np
from pylab import savefig
import matplotlib.pyplot as plt
import matplotlib.cm as cm

stats = np.loadtxt('data/stats_AH1.dat')
timedata = stats[:,0]
dt = timedata[1] - timedata[0]
Spin = stats[:,5]
Mass = stats[:,3]
DeltaSpin = Spin-Spin[0]
DeltaMass = Mass-Mass[0]

plt.plot(timedata,DeltaSpin,color=cm.Reds(8./10.,1.),
         label='$a - a(t=0)$')
plt.plot(timedata,DeltaMass,color=cm.Blues(7./10.,1),ls='--',
         label=r'$M - M(t=0)$')
plt.ylim(-0.00015,0.00015)
plt.legend(ncol=2,fontsize=14, bbox_to_anchor=(0., 0.95, 1., 0.102), loc='lower left')
plt.xlabel(r'$t/M$', fontsize=14)
plt.savefig('DeltaMassAndSpin.png',dpi=256, bbox_inches='tight',pad_inches = 0.1)
plt.close()

To produce a plot like this:

An example for plotting and visualising the surfaces can be found here.

An example of running the AHFinder on existing hdf5 files is here.

Further resources

The following Google drive link (more up to date) or pdf slides illustrate how to use the apparent horizon finder in more detail (thanks to Tiago França for these!) There are also some more advanced tips on using the AHFinder.

• AHFinder Google Slides
• AHFinder PDF slides