FUKA
|
Here lies the binary black hole initial data solver and diagnostic code. The initial data is constructed using maximal slicing with a flat spacial metric. Therefore, the dimensionless spins of the black holes are limited to approximately [-0.84, 0.84]
. The v2 code is a vast improvement over the v1 code as it uses super-imposed BH solutions to initialize the binary instead of building the binary from flat space. When comparing to FUKAv1, generating equal-mass non-rotating ID is roughly >4x
faster due to the reduced number of stages. When comparing to FUKAv1, generating arbitrary spin and unequal mass, the cost savings is roughly (4 + N)x
faster where N is, in the case of FUKAv1, the number of iterative solutions needed to achieve a given mass ratio and spin.
Note:
Mtot
below, we will be referring to the sum of the individual Christodoulou masses Mtot := (MCH_MINUS + MCH_PLUS)
where plus and minus simply refer to their location on the x-axis.q <= 1
is built into the FUKAv2 codes
CMakeLists.txt
is the file needed by CMake to compile the codescompile
is a symbolic link to the script stored in $HOME_KADATH/build_release
to ease compilingsrc
directory contains the relevant source files:solve.cpp
: the one and done solve codereader.cpp
: the reader can provide diagonstics from ID solutions that are computed from the IDsolve
for the first timempirun ./bin/Release/solve initial_bbh.info
solve
for the first timempirun ./bin/Release/solve initial_bbh.info
BBH_ECC_RED.10.0.0.1.q1.3.3.09.<info/dat>
We can deconstruct the name to make it understandable:
BBH_ECC_RED.
denotes a converged BBH solution after the eccentricity reduction stage is completed which uses 3.5th order PN estimates for the orbital frequency and radial infall velocity. This is meant to distinguish the solution from earlier stages which will be discussed later. This also distinquishes it from checkpoints that can be turned on which are saved to file during each iteration of the solver10
: separation distance in geometric units!0.0.
: In this case, both BHs are not spinning1.q1
: the total mass and mass ratio are equal to one0.0.09
: nshells = 0
for bh1
, nshells = 0
for bh2
where the resolution in each domain is 9
collocation points in the radial and theta direction with 8
points in the phi direction.info
: The info file contains all the steering parameters, values used in creating the domain decomposition, stored fields, stages, settings, and controls.dat
: The dat file is a binary file that contains the numerical space and the variable fieldsThe default configuration for BBH has a total mass of 1M and non-spinning. Aside from the diagnostics observed during the solver stage, we can use the reader to verify the ID. This can be done by running:
./bin/Release/reader BBH_ECC_RED.10.0.0.1.q1.0.0.09.info
Which results in the following:
The first two blocks contain information related to the component BHs. These details are covered in the BH README. The only additional parameter is the Center_COM
. This is the coordinate center of each object when shifted by the "center-of-mass" of the binary or, more specifically, the location of the axis of rotation for the binary that can approximate a quasi-stationary solution.
The third block contains information specifically related to the binary
res
: resolution in each numerical domain [r,theta,phi]
q
: mass ratio such that q <= 1
Separation
: coordinate separation d [d/Mtot]
Orbital Omega
: orbital frequency of the binaryE_b / mu
: Binding energy normalized by the reduced massCOM<x/y>
: shift in the coordinates to find a helical killing vectorA-COM<x/y/z>
: A numerical calculation of the COM based on the analytical prescription from Osokine+ (REF)Note: The Diff
noted by the ADM mass is the symmetric difference between the ADM and Komar mass.
Using your favorite text editor, you can open up the initial_bbh.info
. We will go through the file, but we'll discuss only the details relevant to the BBH case. For details on all the parameters you can see more in Configurator README.
Notes:
initial_bbh.info
. Using old initial data unless for very small changes in chi
is inefficient.solve fullto obtain the full Config file. Although useful for development, there is little advantage to using the full Config.
The above includes parameters that must be fixed by the user. The parameters for each BH are simply copied from the isolated solution which can be read in detail in the BH README - the same fixing applies also in the BBH.
The fixing parameters most relevant to the binary are
distance
: this is in geometric units! It is important to pick something reasonable. A general rule for a binary with a few orbits is distance = 8 * Mtot
, however this strongly depends on q
and the spins of component objectsouter_shells
: This allows for additional shells to be placed near the compactified domain. This can be helpful for more accurate quasi-equilibrium ID at lower resolution, but otherwise can be ignored and left to zerores
: global resolution of the binary!adot
: (optional) This is the radial infall velocity parameter when performing eccentricity reduction. This will be discussed more in the eccentricity reduction section belowecc_omega
: (optional) This is the fixed orbital velocity parameter used when performing eccentricity reduction. This will be discussed more in the eccentricity reduction section below Within the full Config or the output solutions, Fields
document the fields that are used in the solver and stored in the dat
file. Changing this has no impact.
For the BBH, only the TOTAL_BC
and ECC_RED
stages are relevant. Old stages available in the original v1 solver are since deprecated.
checkpoint
: this will result in checkpoints being saved to file during each solving iteration - mainly helpful for high resolution binary ID where walltimes or server failures are a concern prior to a converged solution being obtainedcorot_binary
: the objects are no longer fixed based on chi
and instead provide a corotating ID solutionfixed_lapse
: toggling this control enables a fixed lapse on the horizon - not recommendedsequences
: this toggle is enabled by default and essentially tells the driver routine to start from scratch. If this is enabled when attempting to use a previous solution, the previous fields and numerical space (i.e. the dat
file) is ignoreduse_pn
: toggle whether to always use 3.5PN estimates. It is important to ensure this is off if the user wants to specify their own adot
and global_omega
parameters by hand (e.g. for iterative eccentricity reduction)resolve
: force all implicit compact object solutions to be resolved regardless of an existing previous solutioncentralized_cos
: stores all implicit COs into $HOME_KADATH/COs
solver_max_iterations
: set the number of iterations not to exceedsolver_precision
: Determines what is the maximum precision allowed by the solver that determines whether or not the solution has convergedinitial_resolution
: by default all solutions are ran at at a default resolution of [9,9,8]
prior to regridding to a higher resolution. In the event one wants to increase the default initial_resolution
, it can be done hereNow that you've generated the simplest case and we have a better understanding of the config file, we can try something more interesting
res 11
mch 0.1
chi -0.5
mch 0.9
chi 0.85
mpirun ./bin/Release/solve initial_bh.info
This time around we see the iterative chi
increase being done for the primary BH, but overall the only changes observed are related to the isolated BH solvers (see the BH README for details), but the binary solver itself is consistent.
This results in the converged dataset of BBH_ECC_RED.10.-0.5.0.85.1.q0.111111.2.0.11.info/dat
We can of course verify that the ID matches our expectation using
./bin/Release/reader BBH_ECC_RED.10.-0.5.0.85.1.q0.111111.2.0.11.info
To generate the initial setup for BBH ID, we first need to make some guesses based on the input Christodoulou masses and the coordinate separation
COM
- An estimate of the center-of-mass: this is purely Newtonianglobal_omega
- An estimate of the orbital frequency: 3.5th PN estimateOnce these estimates are computed an interface code is ran which
global_omega
At this point, the binary numerical space and fields are constructed and the isolated solutions are interpolated onto the new grid using the idea of super-imposed solutions. Specifically:
w
is chosen such that w = distance / 2 =: decay_limit
decay_rate = exp(-(r_BH / w)^4)
where r_BH
is the coordinate distance to the respective BH centerFor example, if we wanted to compute the initial guess for the lapse at a given point x
, it would be
lapse(x) = 1. + decay_rate_BH1 * (lapse_BH1 - 1.) + decay_rate_BH2 * (lapse_BH2 - 1.)
This is then repeated for all fields in all numerical domains, with the compactified domain set to the asymptotic values of the fields, i.e. lapse = psi = 1
, shift = 0
.
The TOTAL_BC
stage solves the full XCTS system of equations consistently and all at once. The only thing that distinguishes TOTAL_BC
from ECC_RED
is that, by default, TOTAL_BC
fixes the variable global_omega
using the quasi-equilibrium assumption of Madm == Mkomar
where as ECC_RED
uses user-defined value for ecc_omega
and adot
or the built-in 3.5PN estimates. However, when running a new initial data sequence, the global_omega
is initially fixed in the TOTAL_BC
stage using the initial 3.5PN estimate prior to obtaining a quasi-equilibrium solution.
There is of course a very simple reason for this - otherwise the solution would diverge. Put simply, in ID generation of binary objects the shift is very sensitive in such systems and, as such, can cause wild changes in the solution before convergence is obtained, if at all. Additionally, the quasi-equilibrium equations for BBHs are evaluated at "infinity" on the surface of the compatified domain which is incredibly course. To get around these challenges, global_omega
is fixed for one solving stage to allow all the fields to converge to a reasonable initial solution.
After this initial solution is obtained, the fixed_omega
control is deactivated and the TOTAL_BC
stage is reran with the quasi-equilibrium constraint. This is also done to ensure an accurate COM
is found before obtaining 3.5PN estimates of the global_omega
and adot
within the ECC_RED
stage.
In this stage global_omega
and adot
are fixed to generate a less eccentric binary than when using a quasi-equilibrium approximation. After an initial 3.5PN estimate binary has been created, two new parameters appear called ecc_omega
and adot
.
ecc_omega
or adot
are not present in the config file, 3.5PN estimates will always be useduse_pn
is set to on
, the values for ecc_omega
and adot
will ALWAYS be overwrittenecc_omega
and adot
are both set and use_pn
is set to off
, these parameters will only be used in the ecc_red
stage regardless of whether previous stages are ran that change global_omega
adot
and ecc_omega
manually in the initial_bbh.info
file and these values would be used in the ECC_RED
stage. The global_omega
parameter will also be updated to reflect this once a converged solution is obtained.Once the above procedures are completed for the last stage and the input res
is higher than the initial_resolution
, the low resolution solution is interpolated onto the higher resolution numerical grid and is used as the initial guess before running the last stage again. No additional iterative procedures are required at this point.