Line data Source code
1 0 : # Finding apparent horizons in volume data {#examples_find_horizons}
2 :
3 : In this example we will generate some numeric volume data representing a Kerr
4 : black hole and then find its apparent horizon.
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 os
20 : import pandas as pd
21 : from copy import deepcopy
22 : from ruamel.yaml import YAML
23 : ```
24 :
25 :
26 : ```bash
27 : %%bash
28 : # Clean up output files from previous runs
29 : rm -f Kerr*.h5
30 : rm -f FindHorizons*.h5
31 : ```
32 :
33 : First, make sure you have compiled the `SolveXcts` and `FindHorizons3D`
34 : executables in `Release` mode. Put the path to your build directory below.
35 :
36 :
37 : ```python
38 : SPECTRE_BUILD_DIR = "/Users/nlf/Work/spectre/build-Default-Release"
39 : SPECTRE_HOME = "/Users/nlf/Projects/spectre/develop"
40 : ```
41 :
42 : ## Generate Kerr volume data
43 :
44 : We run the `SolveXcts` executable to generate Kerr volume data. We could also
45 : just invoke the Kerr analytic solution and write the data to disk, but the
46 : `SolveXcts` executable already writes volume data in the correct format.
47 :
48 :
49 : ```python
50 : # Load example input file
51 : kerrschild_input_file_path = os.path.join(
52 : SPECTRE_HOME, "tests/InputFiles/Xcts/KerrSchild.yaml"
53 : )
54 : yaml = YAML()
55 : with open(kerrschild_input_file_path, "r") as open_input_file:
56 : kerr_input_file = yaml.load(open_input_file)
57 :
58 : # Modify Kerr-Schild example input file
59 : # - Set some interesting Kerr parameters
60 : kerr_input_file["Background"]["KerrSchild"] = dict(
61 : Mass=1.0, Spin=[0.0, 0.0, 0.9], Center=[0.0, 0.0, 0.0]
62 : )
63 : # - Set the initial guess to the solution to converge quickly
64 : kerr_input_file["InitialGuess"] = kerr_input_file["Background"]
65 : # - Choose domain parameters
66 : domain_params = kerr_input_file["DomainCreator"]["Sphere"]
67 : domain_params["InnerRadius"] = 1.0
68 : domain_params["InitialRefinement"] = 0
69 : domain_params["InitialGridPoints"] = [8, 8, 8]
70 : # - Allow the elliptic solver to converge
71 : kerr_input_file["NonlinearSolver"]["NewtonRaphson"]["ConvergenceCriteria"][
72 : "MaxIterations"
73 : ] = 10
74 : # - Set output file names
75 : kerr_input_file["Observers"]["VolumeFileName"] = "KerrVolume"
76 : kerr_input_file["Observers"]["ReductionFileName"] = "KerrReductions"
77 :
78 : # Write modified input file
79 : with open("Kerr.yaml", "w") as open_input_file:
80 : yaml.dump(kerr_input_file, open_input_file)
81 : ```
82 :
83 :
84 : ```python
85 : SOLVE_XCTS = os.path.join(SPECTRE_BUILD_DIR, "bin/SolveXcts")
86 : !{SOLVE_XCTS} --input-file Kerr.yaml +auto-provision
87 : ```
88 :
89 : Charm++: standalone mode (not using charmrun)
90 : Charm++> Running in Multicore mode: 12 threads (PEs)
91 : Converse/Charm++ Commit ID: v6.10.2-0-g7bf00fa
92 : 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'.
93 : CharmLB> Load balancer assumes all CPUs are same.
94 : HWLOC> [0] Thread 0x7fa168aba180 bound to cpuset: 0x00000001
95 : HWLOC> [10] Thread 0x7fa148a87700 bound to cpuset: 0x00000400
96 : HWLOC> [1] Thread 0x7fa14d290700 bound to cpuset: 0x00000002
97 : HWLOC> [5] Thread 0x7fa14b28c700 bound to cpuset: 0x00000020
98 : HWLOC> [8] Thread 0x7fa149a89700 bound to cpuset: 0x00000100
99 : HWLOC> [2] Thread 0x7fa14ca8f700 bound to cpuset: 0x00000004
100 : HWLOC> [11] Thread 0x7fa148286700 bound to cpuset: 0x00000800
101 : HWLOC> [3] Thread 0x7fa14c28e700 bound to cpuset: 0x00000008
102 : HWLOC> [6] Thread 0x7fa14aa8b700 bound to cpuset: 0x00000040
103 : HWLOC> [4] Thread 0x7fa14ba8d700 bound to cpuset: 0x00000010
104 : HWLOC> [9] Thread 0x7fa149288700 bound to cpuset: 0x00000200
105 : HWLOC> [7] Thread 0x7fa14a28a700 bound to cpuset: 0x00000080
106 : Charm++> Running on 1 hosts (12 sockets x 1 cores x 1 PUs = 12-way SMP)
107 : Charm++> cpu topology info is gathered in 0.003 seconds.
108 :
109 : Executing '/Users/nlf/Work/spectre/build-Default-Release/bin/SolveXcts' using 12 processors.
110 : Charm++ startup time in seconds: 5.040069
111 : Date and time at startup: Tue Mar 1 13:38:03 2022
112 :
113 : SpECTRE Build Information:
114 : Version: 2022.02.17
115 : Compiled on host: 2a7fcfdc5d9f
116 : Compiled in directory: /Users/nlf/Work/spectre/build-Default-Release
117 : Source directory is: /Users/nlf/Projects/spectre/develop
118 : Compiled on git branch: find_horizons_exec
119 : Compiled on git revision: a702fcd88e2
120 : Linked on: Wed Feb 23 22:29:46 2022
121 :
122 : The following options differ from their suggested values:
123 :
124 : Option parsing completed.
125 : Multigrid level 0 has 6 elements in 6 blocks distributed on 12 procs.
126 : NewtonRaphson initialized with residual: 2.692374e+00
127 : Gmres initialized with residual: 2.692374e+00
128 : Gmres(1) iteration complete. Remaining residual: 1.585652e-02
129 : Gmres(2) iteration complete. Remaining residual: 3.250109e-03
130 : Gmres(3) iteration complete. Remaining residual: 6.606992e-04
131 : Gmres(4) iteration complete. Remaining residual: 1.151248e-04
132 : Gmres has converged in 4 iterations: RelativeResidual - The residual magnitude has decreased to a fraction of 0.0001 of its initial value or below (4.27596e-05).
133 : NewtonRaphson(1) iteration complete (0 globalization steps, step length 1). Remaining residual: 1.430971e-01
134 : Gmres initialized with residual: 1.430971e-01
135 : Gmres(1) iteration complete. Remaining residual: 2.573815e-03
136 : Gmres(2) iteration complete. Remaining residual: 5.880012e-04
137 : Gmres(3) iteration complete. Remaining residual: 1.058190e-04
138 : Gmres(4) iteration complete. Remaining residual: 2.932433e-05
139 : Gmres(5) iteration complete. Remaining residual: 4.199566e-06
140 : Gmres has converged in 5 iterations: RelativeResidual - The residual magnitude has decreased to a fraction of 0.0001 of its initial value or below (2.93477e-05).
141 : NewtonRaphson(2) iteration complete (0 globalization steps, step length 1). Remaining residual: 1.914763e-05
142 : Gmres initialized with residual: 1.914763e-05
143 : Gmres(1) iteration complete. Remaining residual: 2.927570e-06
144 : Gmres(2) iteration complete. Remaining residual: 7.516803e-07
145 : Gmres(3) iteration complete. Remaining residual: 1.695121e-07
146 : Gmres(4) iteration complete. Remaining residual: 3.852025e-08
147 : Gmres(5) iteration complete. Remaining residual: 6.239554e-09
148 : Gmres(6) iteration complete. Remaining residual: 3.306899e-10
149 : Gmres has converged in 6 iterations: RelativeResidual - The residual magnitude has decreased to a fraction of 0.0001 of its initial value or below (1.72705e-05).
150 : NewtonRaphson(3) iteration complete (0 globalization steps, step length 1). Remaining residual: 3.352486e-10
151 : Gmres initialized with residual: 3.352486e-10
152 : Gmres(1) iteration complete. Remaining residual: 2.969095e-11
153 : Gmres(2) iteration complete. Remaining residual: 8.169831e-12
154 : Gmres(3) iteration complete. Remaining residual: 1.702308e-12
155 : Gmres(4) iteration complete. Remaining residual: 4.723995e-13
156 : Gmres has converged in 4 iterations: AbsoluteResidual - The residual magnitude has decreased to 1e-12 or below (4.72399e-13).
157 : NewtonRaphson(4) iteration complete (0 globalization steps, step length 1). Remaining residual: 4.951316e-13
158 : NewtonRaphson has converged in 4 iterations: AbsoluteResidual - The residual magnitude has decreased to 1e-10 or below (4.95132e-13).
159 :
160 : Done!
161 : Wall time in seconds: 25.543363
162 : Date and time at completion: Tue Mar 1 13:38:24 2022
163 :
164 : [Partition 0][Node 0] End of program
165 :
166 :
167 : ## Find horizons in the volume data
168 :
169 : We run the `FindHorizons3D` executable over the generated volume data:
170 :
171 :
172 : ```python
173 : # Load example input file
174 : horizons_input_file_path = os.path.join(
175 : SPECTRE_HOME, "tests/InputFiles/FindHorizons/FindHorizons3D.yaml"
176 : )
177 : with open(horizons_input_file_path, "r") as open_input_file:
178 : horizons_input_file = yaml.load(open_input_file)
179 :
180 : # Modify Kerr-Schild example input file
181 : # - Select the same domain that we solved above
182 : horizons_input_file["DomainCreator"] = deepcopy(
183 : kerr_input_file["DomainCreator"]
184 : )
185 : del horizons_input_file["DomainCreator"]["Sphere"]["OuterBoundaryCondition"]
186 : horizons_input_file["DomainCreator"]["Sphere"]["Interior"] = {"Excise": True}
187 : # - Set importer file names
188 : importer_params = horizons_input_file["Importers"]["VolumeData"]
189 : importer_params["FileGlob"] = "KerrVolume*.h5"
190 : importer_params["ObservationValue"] = "Last"
191 : # - AH finder parameters
192 : ah_params = horizons_input_file["ApparentHorizons"]["AhA"]
193 : ah_params["InitialGuess"]["Radius"] = 2.0
194 : # - Set output file names
195 : horizons_input_file["Observers"]["VolumeFileName"] = "FindHorizonsVolume"
196 : horizons_input_file["Observers"]["ReductionFileName"] = "FindHorizonsReductions"
197 : horizons_input_file["Observers"]["SurfaceFileName"] = "FindHorizonsSurfaces"
198 :
199 : # Write modified input file
200 : with open("FindHorizons.yaml", "w") as open_input_file:
201 : yaml.dump(horizons_input_file, open_input_file)
202 : ```
203 :
204 :
205 : ```python
206 : FIND_HORIZONS = os.path.join(SPECTRE_BUILD_DIR, "bin/FindHorizons3D")
207 : !{FIND_HORIZONS} --input-file FindHorizons.yaml +auto-provision
208 : ```
209 :
210 : Charm++: standalone mode (not using charmrun)
211 : Charm++> Running in Multicore mode: 12 threads (PEs)
212 : Converse/Charm++ Commit ID: v6.10.2-0-g7bf00fa
213 : 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'.
214 : CharmLB> Load balancer assumes all CPUs are same.
215 : HWLOC> [0] Thread 0x7eff0efdb180 bound to cpuset: 0x00000001
216 : HWLOC> [2] Thread 0x7efef2fb0700 bound to cpuset: 0x00000004
217 : HWLOC> [1] Thread 0x7efef37b1700 bound to cpuset: 0x00000002
218 : HWLOC> [8] Thread 0x7efeeffaa700 bound to cpuset: 0x00000100
219 : HWLOC> [9] Thread 0x7efeef7a9700 bound to cpuset: 0x00000200
220 : HWLOC> [6] Thread 0x7efef0fac700 bound to cpuset: 0x00000040
221 : HWLOC> [7] Thread 0x7efef07ab700 bound to cpuset: 0x00000080
222 : HWLOC> [4] Thread 0x7efef1fae700 bound to cpuset: 0x00000010
223 : HWLOC> [10] Thread 0x7efeeefa8700 bound to cpuset: 0x00000400
224 : HWLOC> [5] Thread 0x7efef17ad700 bound to cpuset: 0x00000020
225 : HWLOC> [3] Thread 0x7efef27af700 bound to cpuset: 0x00000008
226 : HWLOC> [11] Thread 0x7efeee7a7700 bound to cpuset: 0x00000800
227 : Charm++> Running on 1 hosts (12 sockets x 1 cores x 1 PUs = 12-way SMP)
228 : Charm++> cpu topology info is gathered in 0.009 seconds.
229 :
230 : Executing '/Users/nlf/Work/spectre/build-Default-Release/bin/FindHorizons3D' using 12 processors.
231 : Charm++ startup time in seconds: 0.097010
232 : Date and time at startup: Tue Mar 1 13:38:25 2022
233 :
234 : SpECTRE Build Information:
235 : Version: 2022.02.17
236 : Compiled on host: 2a7fcfdc5d9f
237 : Compiled in directory: /Users/nlf/Work/spectre/build-Default-Release
238 : Source directory is: /Users/nlf/Projects/spectre/develop
239 : Compiled on git branch: find_horizons_exec
240 : Compiled on git revision: a702fcd88e2
241 : Linked on: Wed Feb 23 22:29:46 2022
242 :
243 : The following options differ from their suggested values:
244 :
245 : Option parsing completed.
246 : AhA: t=0: its=0: -5.2e-01<R<1e+00, |R|=0.5, |R_grid|=0.5, 2<r<2
247 : AhA: t=0: its=1: -2.8e-03<R<4e-01, |R|=0.2, |R_grid|=0.2, 1.757<r<1.93
248 : AhA: t=0: its=2: -1.2e-02<R<2e-01, |R|=0.1, |R_grid|=0.1, 1.653<r<1.86
249 : AhA: t=0: its=3: -1.6e-02<R<2e-01, |R|=0.07, |R_grid|=0.08, 1.591<r<1.813
250 : AhA: t=0: its=4: -2.0e-02<R<1e-01, |R|=0.05, |R_grid|=0.06, 1.55<r<1.78
251 : AhA: t=0: its=5: -2.3e-02<R<1e-01, |R|=0.04, |R_grid|=0.05, 1.521<r<1.758
252 : AhA: t=0: its=6: -4.1e-02<R<8e-02, |R|=0.03, |R_grid|=0.04, 1.501<r<1.742
253 : AhA: t=0: its=7: -5.7e-02<R<7e-02, |R|=0.02, |R_grid|=0.03, 1.486<r<1.73
254 : AhA: t=0: its=8: -6.8e-02<R<7e-02, |R|=0.01, |R_grid|=0.03, 1.475<r<1.721
255 : AhA: t=0: its=9: -7.7e-02<R<6e-02, |R|=0.01, |R_grid|=0.03, 1.467<r<1.715
256 : AhA: t=0: its=10: -8.3e-02<R<6e-02, |R|=0.008, |R_grid|=0.03, 1.461<r<1.71
257 : AhA: t=0: its=11: -8.7e-02<R<6e-02, |R|=0.006, |R_grid|=0.03, 1.456<r<1.707
258 : AhA: t=0: its=12: -9.1e-02<R<5e-02, |R|=0.005, |R_grid|=0.02, 1.453<r<1.704
259 : AhA: t=0: its=13: -9.4e-02<R<5e-02, |R|=0.004, |R_grid|=0.02, 1.45<r<1.702
260 : AhA: t=0: its=14: -9.6e-02<R<5e-02, |R|=0.003, |R_grid|=0.02, 1.448<r<1.701
261 : AhA: t=0: its=15: -9.7e-02<R<5e-02, |R|=0.002, |R_grid|=0.02, 1.447<r<1.699
262 : AhA: t=0: its=16: -9.8e-02<R<5e-02, |R|=0.002, |R_grid|=0.02, 1.446<r<1.699
263 : AhA: t=0: its=17: -9.9e-02<R<5e-02, |R|=0.001, |R_grid|=0.02, 1.445<r<1.698
264 : AhA: t=0: its=18: -1.0e-01<R<5e-02, |R|=0.0009, |R_grid|=0.02, 1.444<r<1.697
265 : AhA: t=0: its=19: -1.0e-01<R<5e-02, |R|=0.0007, |R_grid|=0.02, 1.444<r<1.697
266 : AhA: t=0: its=20: -1.0e-01<R<5e-02, |R|=0.0005, |R_grid|=0.02, 1.443<r<1.697
267 : AhA: t=0: its=21: -1.0e-01<R<5e-02, |R|=0.0004, |R_grid|=0.02, 1.443<r<1.697
268 : AhA: t=0: its=22: -1.0e-01<R<5e-02, |R|=0.0003, |R_grid|=0.02, 1.443<r<1.696
269 : AhA: t=0: its=23: -1.0e-01<R<5e-02, |R|=0.0002, |R_grid|=0.02, 1.443<r<1.696
270 :
271 : Done!
272 : Wall time in seconds: 5.510678
273 : Date and time at completion: Tue Mar 1 13:38:31 2022
274 :
275 : [Partition 0][Node 0] End of program
276 :
277 :
278 : ## Horizon quantities
279 :
280 : The horizon finder measures a few surface quantities:
281 :
282 :
283 : ```python
284 : # These routines read the data and process them a bit. You can skip them to see
285 : # the results below.
286 :
287 :
288 : def load_dataset(subfile):
289 : legend = subfile.attrs["Legend"]
290 : return pd.DataFrame(data=subfile, columns=legend).set_index(legend[0])
291 : ```
292 :
293 :
294 : ```python
295 : with h5py.File("FindHorizonsReductions.h5", "r") as open_h5_file:
296 : ah = load_dataset(open_h5_file["AhA.dat"])
297 : ```
298 :
299 :
300 : ```python
301 : print(ah.to_string())
302 : ```
303 :
304 : Area IrreducibleMass MaxRicciScalar MinRicciScalar ChristodoulouMass DimensionlessSpinMagnitude
305 : Time
306 : 0.0 36.03551 0.846702 1.734013 -0.082303 0.96325 0.83824
307 :
|
