Line data    Source code

       1           0 : # Binary black hole initial data {#example_bbh_id}
       2             : 
       3             : In this example we run the elliptic solver to compute initial data for binary
       4             : black holes.
       5             : 
       6             : 
       7             : ```python
       8             : # Distributed under the MIT License.
       9             : # See LICENSE.txt for details.
      10             : 
      11             : # Dependencies:
      12             : %pip install numpy matplotlib pandas 'h5py>=3.0.0' ruamel.yaml
      13             : ```
      14             : 
      15             : 
      16             : ```python
      17             : import h5py
      18             : import matplotlib.pyplot as plt
      19             : import multiprocessing
      20             : import numpy as np
      21             : import os
      22             : import pandas as pd
      23             : from ruamel.yaml import YAML
      24             : 
      25             : plt.style.use("../plots.mplstyle")
      26             : ```
      27             : 
      28             : 
      29             : ```bash
      30             : %%bash
      31             : # Clean up output files from previous runs
      32             : rm -f Bbh*.h5
      33             : ```
      34             : 
      35             : First, make sure you have compiled the `SolveXcts` executable. Put in the path
      36             : to your build directory below. Make sure you have compiled in `Release` mode.
      37             : 
      38             : 
      39             : ```python
      40             : SPECTRE_BUILD_DIR = "/Users/nlf/Work/spectre/build-Default-Release"
      41             : SPECTRE_HOME = "/Users/nlf/Projects/spectre/develop"
      42             : ```
      43             : 
      44             : 
      45             : ## Setup input file
      46             : 
      47             : We set up an input file based on the
      48             : `tests/InputFiles/Xcts/BinaryBlackHole.yaml` example:
      49             : 
      50             : 
      51             : ```python
      52             : # Load example input file
      53             : load_input_file_path = os.path.join(
      54             :     SPECTRE_HOME, "tests/InputFiles/Xcts/BinaryBlackHole.yaml"
      55             : )
      56             : yaml = YAML()
      57             : with open(load_input_file_path, "r") as open_input_file:
      58             :     input_file = yaml.load(open_input_file)
      59             : 
      60             : # Modify example input file
      61             : # - Set output file names
      62             : input_file["Observers"]["VolumeFileName"] = "BbhVolume"
      63             : input_file["Observers"]["ReductionFileName"] = "BbhReductions"
      64             : 
      65             : # Write modified input file
      66             : with open("Bbh.yaml", "w") as open_input_file:
      67             :     yaml.dump(input_file, open_input_file)
      68             : ```
      69             : 
      70             : ## Run executable
      71             : 
      72             : We pass the input file to the `SolveXcts` executable to solve the elliptic
      73             : problem. It will take a few minutes to complete on ~10 cores. Adapt the command
      74             : below to your system, or run the `SolveXcts` executable with the `Bbh.yaml`
      75             : input file manually.
      76             : 
      77             : 
      78             : ```python
      79             : NUM_CORES = multiprocessing.cpu_count()
      80             : SOLVE_XCTS = os.path.join(SPECTRE_BUILD_DIR, "bin/SolveXcts")
      81             : !{SOLVE_XCTS} --input-file Bbh.yaml +p {NUM_CORES}
      82             : ```
      83             : 
      84             :     Charm++: standalone mode (not using charmrun)
      85             :     Charm++> Running in Multicore mode: 12 threads (PEs)
      86             :     Converse/Charm++ Commit ID: v6.10.2-0-g7bf00fa
      87             :     Warning> Randomization of virtual memory (ASLR) is turned on in the kernel, thread migration may not work! Run 'echo 0 > /proc/sys/kernel/randomize_va_space' as root to disable it, or try running with '+isomalloc_sync'.
      88             :     CharmLB> Load balancer assumes all CPUs are same.
      89             :     Charm++> Running on 1 hosts (12 sockets x 1 cores x 1 PUs = 12-way SMP)
      90             :     Charm++> cpu topology info is gathered in 0.002 seconds.
      91             :     
      92             :     Executing '/Users/nlf/Work/spectre/build-Default-Release/bin/SolveXcts' using 12 processors.
      93             :     Charm++ startup time in seconds: 3.621687
      94             :     Date and time at startup: Tue Feb 22 15:08:18 2022
      95             :     
      96             :     SpECTRE Build Information:
      97             :     Version:                      2022.02.17
      98             :     Compiled on host:             98c21919c2ac
      99             :     Compiled in directory:        /Users/nlf/Work/spectre/build-Default-Release
     100             :     Source directory is:          /Users/nlf/Projects/spectre/develop
     101             :     Compiled on git branch:       xcts_input_files
     102             :     Compiled on git revision:     db7100f700b
     103             :     Linked on:                    Thu Jan 27 17:45:20 2022
     104             :     
     105             :     The following options differ from their suggested values:
     106             :     
     107             :     Option parsing completed.
     108             :     Multigrid level 0 has 232 elements in 54 blocks distributed on 12 procs.
     109             :     Multigrid level 1 has 54 elements in 54 blocks distributed on 12 procs.
     110             :     NewtonRaphson initialized with residual: 1.058476e+01
     111             :     Gmres initialized with residual: 1.058476e+01
     112             :     Gmres(1) iteration complete. Remaining residual: 2.003738e-01
     113             :     Gmres(2) iteration complete. Remaining residual: 1.206809e-01
     114             :     Gmres(3) iteration complete. Remaining residual: 8.694367e-02
     115             :     Gmres(4) iteration complete. Remaining residual: 7.275577e-02
     116             :     Gmres(5) iteration complete. Remaining residual: 4.520779e-02
     117             :     Gmres(6) iteration complete. Remaining residual: 1.464362e-02
     118             :     Gmres(7) iteration complete. Remaining residual: 6.754772e-03
     119             :     Gmres has converged in 7 iterations: RelativeResidual - The residual magnitude has decreased to a fraction of 0.001 of its initial value or below (0.00063816).
     120             :     NewtonRaphson(1) iteration complete (0 globalization steps, step length 1). Remaining residual: 2.865222e+00
     121             :     Gmres initialized with residual: 2.865222e+00
     122             :     Gmres(1) iteration complete. Remaining residual: 4.062729e-02
     123             :     Gmres(2) iteration complete. Remaining residual: 2.455538e-02
     124             :     Gmres(3) iteration complete. Remaining residual: 1.602018e-02
     125             :     Gmres(4) iteration complete. Remaining residual: 1.181106e-02
     126             :     Gmres(5) iteration complete. Remaining residual: 5.591212e-03
     127             :     Gmres(6) iteration complete. Remaining residual: 2.012421e-03
     128             :     Gmres has converged in 6 iterations: RelativeResidual - The residual magnitude has decreased to a fraction of 0.001 of its initial value or below (0.000702361).
     129             :     NewtonRaphson(2) iteration complete (0 globalization steps, step length 1). Remaining residual: 8.164886e-03
     130             :     Gmres initialized with residual: 8.164886e-03
     131             :     Gmres(1) iteration complete. Remaining residual: 1.279278e-03
     132             :     Gmres(2) iteration complete. Remaining residual: 7.680224e-04
     133             :     Gmres(3) iteration complete. Remaining residual: 6.026001e-04
     134             :     Gmres(4) iteration complete. Remaining residual: 5.116235e-04
     135             :     Gmres(5) iteration complete. Remaining residual: 3.771227e-04
     136             :     Gmres(6) iteration complete. Remaining residual: 1.521095e-04
     137             :     Gmres(7) iteration complete. Remaining residual: 4.029530e-05
     138             :     Gmres(8) iteration complete. Remaining residual: 1.632393e-05
     139             :     Gmres(9) iteration complete. Remaining residual: 5.028700e-06
     140             :     Gmres has converged in 9 iterations: RelativeResidual - The residual magnitude has decreased to a fraction of 0.001 of its initial value or below (0.000615893).
     141             :     NewtonRaphson(3) iteration complete (0 globalization steps, step length 1). Remaining residual: 1.690461e-03
     142             :     Gmres initialized with residual: 1.690461e-03
     143             :     Gmres(1) iteration complete. Remaining residual: 6.584416e-06
     144             :     Gmres(2) iteration complete. Remaining residual: 5.651935e-06
     145             :     Gmres(3) iteration complete. Remaining residual: 5.288583e-06
     146             :     Gmres(4) iteration complete. Remaining residual: 4.922907e-06
     147             :     Gmres(5) iteration complete. Remaining residual: 4.250901e-06
     148             :     Gmres(6) iteration complete. Remaining residual: 1.596286e-06
     149             :     Gmres has converged in 6 iterations: RelativeResidual - The residual magnitude has decreased to a fraction of 0.001 of its initial value or below (0.00094429).
     150             :     NewtonRaphson(4) iteration complete (0 globalization steps, step length 1). Remaining residual: 1.606242e-06
     151             :     Gmres initialized with residual: 1.606242e-06
     152             :     Gmres(1) iteration complete. Remaining residual: 7.008839e-07
     153             :     Gmres(2) iteration complete. Remaining residual: 5.285042e-07
     154             :     Gmres(3) iteration complete. Remaining residual: 4.919028e-07
     155             :     Gmres(4) iteration complete. Remaining residual: 3.871318e-07
     156             :     Gmres(5) iteration complete. Remaining residual: 1.603757e-07
     157             :     Gmres(6) iteration complete. Remaining residual: 2.922369e-08
     158             :     Gmres(7) iteration complete. Remaining residual: 1.438336e-08
     159             :     Gmres(8) iteration complete. Remaining residual: 5.271254e-09
     160             :     Gmres(9) iteration complete. Remaining residual: 2.125873e-09
     161             :     Gmres(10) iteration complete. Remaining residual: 1.044162e-09
     162             :     Gmres has converged in 10 iterations: RelativeResidual - The residual magnitude has decreased to a fraction of 0.001 of its initial value or below (0.000650065).
     163             :     NewtonRaphson(5) iteration complete (0 globalization steps, step length 1). Remaining residual: 1.066508e-09
     164             :     Gmres initialized with residual: 1.066508e-09
     165             :     Gmres(1) iteration complete. Remaining residual: 8.000055e-10
     166             :     Gmres(2) iteration complete. Remaining residual: 5.991477e-10
     167             :     Gmres(3) iteration complete. Remaining residual: 4.073811e-10
     168             :     Gmres(4) iteration complete. Remaining residual: 2.559880e-10
     169             :     Gmres(5) iteration complete. Remaining residual: 1.242268e-10
     170             :     Gmres(6) iteration complete. Remaining residual: 3.521815e-11
     171             :     Gmres has converged in 6 iterations: AbsoluteResidual - The residual magnitude has decreased to 1e-10 or below (3.52181e-11).
     172             :     NewtonRaphson(6) iteration complete (0 globalization steps, step length 1). Remaining residual: 3.589802e-11
     173             :     NewtonRaphson has converged in 6 iterations: AbsoluteResidual - The residual magnitude has decreased to 1e-10 or below (3.5898e-11).
     174             :     
     175             :     Done!
     176             :     Wall time in seconds: 769.172052
     177             :     Date and time at completion: Tue Feb 22 15:21:03 2022
     178             :     
     179             :     [Partition 0][Node 0] End of program
     180             : 
     181             : 
     182             : ## Load initial data into evolutions
     183             : 
     184             : The executable has created H5 output files with volume data and diagnostics:
     185             : 
     186             : 
     187             : ```python
     188             : !ls *.h5
     189             : ```
     190             : 
     191             :     BbhReductions.h5  BbhVolume0.h5
     192             : 
     193             : 
     194             : The volume data in the `BbhVolume*.h5` files (one per node) can be imported into
     195             : evolution executables. It contains the following data:
     196             : 
     197             : 
     198             : ```python
     199             : !h5ls -r BbhVolume*.h5
     200             : ```
     201             : 
     202             :     /                        Group
     203             :     /VolumeData.vol          Group
     204             :     /VolumeData.vol/ObservationId1423324405707320626 Group
     205             :     /VolumeData.vol/ObservationId1423324405707320626/ConformalFactor Dataset {61344}
     206             :     /VolumeData.vol/ObservationId1423324405707320626/ExtrinsicCurvature_xx Dataset {61344}
     207             :     /VolumeData.vol/ObservationId1423324405707320626/ExtrinsicCurvature_yx Dataset {61344}
     208             :     /VolumeData.vol/ObservationId1423324405707320626/ExtrinsicCurvature_yy Dataset {61344}
     209             :     /VolumeData.vol/ObservationId1423324405707320626/ExtrinsicCurvature_zx Dataset {61344}
     210             :     /VolumeData.vol/ObservationId1423324405707320626/ExtrinsicCurvature_zy Dataset {61344}
     211             :     /VolumeData.vol/ObservationId1423324405707320626/ExtrinsicCurvature_zz Dataset {61344}
     212             :     /VolumeData.vol/ObservationId1423324405707320626/HamiltonianConstraint Dataset {61344}
     213             :     /VolumeData.vol/ObservationId1423324405707320626/InertialCoordinates_x Dataset {61344}
     214             :     /VolumeData.vol/ObservationId1423324405707320626/InertialCoordinates_y Dataset {61344}
     215             :     /VolumeData.vol/ObservationId1423324405707320626/InertialCoordinates_z Dataset {61344}
     216             :     /VolumeData.vol/ObservationId1423324405707320626/Lapse Dataset {61344}
     217             :     /VolumeData.vol/ObservationId1423324405707320626/MomentumConstraint_x Dataset {61344}
     218             :     /VolumeData.vol/ObservationId1423324405707320626/MomentumConstraint_y Dataset {61344}
     219             :     /VolumeData.vol/ObservationId1423324405707320626/MomentumConstraint_z Dataset {61344}
     220             :     /VolumeData.vol/ObservationId1423324405707320626/ShiftExcess_x Dataset {61344}
     221             :     /VolumeData.vol/ObservationId1423324405707320626/ShiftExcess_y Dataset {61344}
     222             :     /VolumeData.vol/ObservationId1423324405707320626/ShiftExcess_z Dataset {61344}
     223             :     /VolumeData.vol/ObservationId1423324405707320626/Shift_x Dataset {61344}
     224             :     /VolumeData.vol/ObservationId1423324405707320626/Shift_y Dataset {61344}
     225             :     /VolumeData.vol/ObservationId1423324405707320626/Shift_z Dataset {61344}
     226             :     /VolumeData.vol/ObservationId1423324405707320626/SpatialMetric_xx Dataset {61344}
     227             :     /VolumeData.vol/ObservationId1423324405707320626/SpatialMetric_yx Dataset {61344}
     228             :     /VolumeData.vol/ObservationId1423324405707320626/SpatialMetric_yy Dataset {61344}
     229             :     /VolumeData.vol/ObservationId1423324405707320626/SpatialMetric_zx Dataset {61344}
     230             :     /VolumeData.vol/ObservationId1423324405707320626/SpatialMetric_zy Dataset {61344}
     231             :     /VolumeData.vol/ObservationId1423324405707320626/SpatialMetric_zz Dataset {61344}
     232             :     /VolumeData.vol/ObservationId1423324405707320626/bases Dataset {696}
     233             :     /VolumeData.vol/ObservationId1423324405707320626/connectivity Dataset {294400}
     234             :     /VolumeData.vol/ObservationId1423324405707320626/grid_names Dataset {5272}
     235             :     /VolumeData.vol/ObservationId1423324405707320626/quadratures Dataset {696}
     236             :     /VolumeData.vol/ObservationId1423324405707320626/total_extents Dataset {696}
     237             :     /src.tar.gz              Dataset {5307433}
     238             : 
     239             : 
     240             : - `Lapse` ($\alpha$), `Shift` ($\beta^i$), `SpatialMetric` ($\gamma_{ij}$) and
     241             :   `ExtrinsicCurvature` ($K_{ij}$): These quantities solve the Einstein
     242             :   constraint equations in "corotating" coordinates, such that, e.g.,
     243             :   $\partial_t K = 0$. The black holes remain approximately at the same
     244             :   coordinate positions when starting to evolve the initial data in these
     245             :   coordinates. Note that $\beta^i \propto r$ in these coordinates, which can be
     246             :   large and hence numerically undesireable.
     247             : - `ShiftExcess` ($\beta_\text{excess}^i$): This shift vector excludes the
     248             :   rotational part, so it is asymptotically small. Import `ShiftExcess` instead
     249             :   of `Shift` in an evolution to obtain coordinates in which the black holes are
     250             :   orbiting. Note that only
     251             :   \begin{equation}
     252             :     \beta^i = \beta_\text{background}^i + \beta_\text{excess}^i
     253             :   \end{equation}
     254             :   fulfills all coordinate conditions imposed on the initial data, such as
     255             :   $\partial_t K = 0$. The background shift is typically
     256             :   \begin{equation}
     257             :     \beta_\text{background}^i = \epsilon^{ijk} \Omega^{j} x^k + \dot{a}_0 x^i
     258             :     \text{,}
     259             :   \end{equation}
     260             :   where $\Omega$ is the orbital angular velocity and $\dot{a}_0$ is the
     261             :   expansion parameter.
     262             :   See the [Xcts::AnalyticData::Binary](https://spectre-code.org/classXcts_1_1AnalyticData_1_1Binary.html)
     263             :   class for details.
     264             : 
     265             : ## Plot diagnostics
     266             : 
     267             : We plot the diagnostics in `BbhReductions.h5` to see what happened during the
     268             : elliptic solve:
     269             : 
     270             : 
     271             : ```python
     272             : # These routines read the data and process them a bit. You can skip to the plot
     273             : # below to see the results.
     274             : 
     275             : 
     276             : def split_iteration_sequence(data):
     277             :     left_bounds = np.where(data.index == 0)[0]
     278             :     right_bounds = list(left_bounds[1:]) + [-1]
     279             :     return [data.iloc[i:j] for i, j in zip(left_bounds, right_bounds)]
     280             : 
     281             : 
     282             : def load_dataset(subfile):
     283             :     legend = subfile.attrs["Legend"]
     284             :     return pd.DataFrame(data=subfile, columns=legend).set_index(legend[0])
     285             : 
     286             : 
     287             : with h5py.File("BbhReductions.h5", "r") as h5_file:
     288             :     nonlinear_residuals = load_dataset(h5_file["NewtonRaphsonResiduals.dat"])
     289             :     all_linear_residuals = split_iteration_sequence(
     290             :         load_dataset(h5_file["GmresResiduals.dat"])
     291             :     )
     292             : ```
     293             : 
     294             : 
     295             : ```python
     296             : # Plot linear residuals
     297             : n = 0
     298             : lin_label = "Linear residual"
     299             : for linear_residuals in all_linear_residuals:
     300             :     plt.semilogy(
     301             :         linear_residuals.index + n,
     302             :         linear_residuals["Residual"]
     303             :         / all_linear_residuals[0]["Residual"].iloc[0],
     304             :         color="black",
     305             :         label=lin_label,
     306             :         marker=".",
     307             :         markersize=2.5,
     308             :         zorder=30,
     309             :     )
     310             :     n += len(linear_residuals) - 1
     311             :     lin_label = None
     312             : 
     313             : # Plot nonlinear residuals
     314             : nonlin_xticks = [0] + list(
     315             :     np.cumsum([len(l) - 1 for l in all_linear_residuals])
     316             : )
     317             : plt.semilogy(
     318             :     nonlin_xticks,
     319             :     nonlinear_residuals["Residual"] / nonlinear_residuals["Residual"].iloc[0],
     320             :     color="black",
     321             :     ls="dotted",
     322             :     marker=".",
     323             :     label="Nonlinear residual",
     324             :     zorder=40,
     325             : )
     326             : plt.xticks(nonlin_xticks)
     327             : xticks = (
     328             :     range(len(all_linear_residuals))
     329             :     if nonlinear_residuals is None
     330             :     else nonlin_xticks
     331             : )
     332             : 
     333             : plt.xlabel("Cumulative linear solver iteration")
     334             : plt.legend();
     335             : ```
     336             : 
     337             : 
     338             :     
     339             : ![png](BbhInitialData_files/BbhInitialData_17_0.png)
     340             :     
     341             : 
     342             : 
     343             : The plot shows the convergence of the nonlinear solver (dotted line), along with
     344             : the convergence of the linear solver that runs in each nonlinear iteration
     345             : (solid line). A few things to note:
     346             : 
     347             : - The residual converges down to almost machine precision, but that doesn't
     348             :   reflect the discretization error of the solution. It only shows that we have
     349             :   solved the discretized problem very accurately. To get an idea of the
     350             :   discretization error we would have to look at quantities such as constraint
     351             :   norms.
     352             : - The Newton-Raphson nonlinear solver converges slowly at first, and then begins
     353             :   to converge quadratically once we are closer to the solution and hence the
     354             :   linearization is more accurate.
     355             : - The linear solver needs only a few steps to converge, which is a feature of
     356             :   our multigrid-Schwarz preconditioner.
     357             : 
     358             : You can read more about the elliptic solver in this paper:
     359             : 
     360             : - N. L. Vu _et al_., A scalable elliptic solver with task-based parallelism for
     361             :   the SpECTRE numerical relativity code (2022),
     362             :   [arXiv:2111.06767](https://arxiv.org/abs/2111.06767)
     363             : 
     364             : 
     365             : ```python
     366             : 
     367             : ```