Line data Source code
1 0 : ## Visualization with Python {#tutorial_vis_python}
2 :
3 : This tutorial shows you how to load 3D volume data in Python and visualize it
4 : using `matplotlib`.
5 :
6 : First, you need access to the `spectre` Python bindings. A simple way to do this
7 : is to run a Jupyter server from the build directory:
8 :
9 : ```sh
10 : ./bin/python-spectre -m jupyterlab
11 : ```
12 :
13 : You can find detailed instruction in the tutorial on \ref spectre_using_python.
14 :
15 : > Note for VSCode users: You can also select this Jupyter server as kernel for
16 : > notebooks running in VSCode (see [docs](https://code.visualstudio.com/docs/datascience/jupyter-notebooks#_connect-to-a-remote-jupyter-server)).
17 :
18 : Now we can import modules from the `spectre` Python package.
19 :
20 :
21 : ```python
22 : # Distributed under the MIT License.
23 : # See LICENSE.txt for details.
24 :
25 : import matplotlib.pyplot as plt
26 : import numpy as np
27 : import os
28 :
29 : plt.style.use("../Examples/plots.mplstyle")
30 : ```
31 :
32 : ## Select an observation
33 :
34 : First, we open the H5 volume data files from a simulation and get a list of all
35 : available observations and their times:
36 :
37 :
38 : ```python
39 : from spectre.Visualization.OpenVolfiles import open_volfiles
40 : from spectre.Visualization.ReadH5 import list_observations
41 : import glob
42 :
43 : h5files = glob.glob(
44 : # Point this to your volume data files
45 : "BbhVolume*.h5"
46 : )
47 : subfile_name = "/VolumeData"
48 : obs_ids, obs_times = list_observations(open_volfiles(h5files, subfile_name))
49 : ```
50 :
51 : We can also list all available variables in the volume data files:
52 :
53 :
54 : ```python
55 : import rich.columns
56 :
57 : for volfile in open_volfiles(h5files, subfile_name, obs_ids[-1]):
58 : all_vars = volfile.list_tensor_components(obs_ids[-1])
59 : # Only look into the first file. The other should have the same variables.
60 : break
61 :
62 : rich.print(rich.columns.Columns(all_vars))
63 : ```
64 :
65 :
66 : <pre style="white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace">ConformalFactor ExtrinsicCurvature_xx ExtrinsicCurvature_yx
67 : ExtrinsicCurvature_yy ExtrinsicCurvature_zx ExtrinsicCurvature_zy
68 : ExtrinsicCurvature_zz InertialCoordinates_x InertialCoordinates_y
69 : InertialCoordinates_z Lapse RadiallyCompressedCoordinates_x
70 : RadiallyCompressedCoordinates_y RadiallyCompressedCoordinates_z ShiftExcess_x
71 : ShiftExcess_y ShiftExcess_z Shift_x
72 : Shift_y Shift_z SpatialMetric_xx
73 : SpatialMetric_yx SpatialMetric_yy SpatialMetric_zx
74 : SpatialMetric_zy SpatialMetric_zz
75 : </pre>
76 :
77 :
78 :
79 : ### Interpolate to a slice
80 :
81 : One way of visualizing your 3D volume data in a 2D plot is with a slice. Use
82 : `interpolate_to_points` to interpolate your volume data to any set of
83 : coordinates:
84 :
85 :
86 : ```python
87 : from spectre.IO.Exporter import interpolate_to_points
88 :
89 : # Coordinates of a regular grid in the xy plane
90 : x, y = np.meshgrid(np.linspace(-15, 15, 300), np.linspace(-15, 15, 300))
91 : z = np.zeros(x.shape)
92 :
93 : # Interpolation can be slow, so it's useful to cache the result. Remember to
94 : # delete the cache file if you change anything.
95 : cache_file = "interpolated_data.npy"
96 : if os.path.isfile(cache_file):
97 : lapse = np.load(cache_file)
98 : else:
99 : lapse = np.array(
100 : interpolate_to_points(
101 : h5files,
102 : subfile_name=subfile_name,
103 : observation_id=obs_ids[-1],
104 : tensor_components=["Lapse"],
105 : target_points=[x.flatten(), y.flatten(), z.flatten()],
106 : )
107 : ).reshape(x.shape)
108 : np.save(cache_file, lapse)
109 : ```
110 :
111 : Now we can plot the 2D slice with `matplotlib`:
112 :
113 :
114 : ```python
115 : # Plot lapse
116 : plt.contourf(x, y, np.log(1 - lapse))
117 :
118 : # Plot circles for black holes
119 : ax = plt.gca()
120 : for bh_pos in [-8, 8]:
121 : ax.add_patch(
122 : plt.Circle(xy=(bh_pos, 0), radius=0.89, color="black", fill=True)
123 : )
124 :
125 : # Make plot square
126 : ax.set_aspect("equal")
127 : ```
128 :
129 : /var/folders/xp/5t2ny359187c4ljckf8sz3m80000gn/T/ipykernel_19858/3375439396.py:2: RuntimeWarning: invalid value encountered in subtract
130 : plt.contourf(x, y, np.log(1 - lapse))
131 :
132 :
133 :
134 :
135 : 
136 :
137 :
138 :
139 : ### Plot individual elements
140 :
141 : You can also iterate over elements in your volume data to plot things like
142 : element boundaries, collocation points, etc. Use `iter_elements`:
143 :
144 :
145 : ```python
146 : from spectre.IO.H5.IterElements import iter_elements
147 :
148 : # Create a 3D plot
149 : ax = plt.gcf().add_subplot(111, projection="3d")
150 : ax.axis("off")
151 :
152 : # Iterate over elements
153 : for element in iter_elements(
154 : open_volfiles(h5files, subfile_name),
155 : obs_ids=obs_ids[0],
156 : element_patterns=["B3,*"], # Only plot elements in block 3
157 : ):
158 : # Draw outline of the element by mapping the edges of the logical cube to
159 : # inertial coordinates using the element map
160 : for d in range(3):
161 : for edge in range(4):
162 : line = np.zeros((3, 100))
163 : line[d, :] = np.linspace(-1, 1, 100)
164 : line[(d + 1) % 3, :] = 2 * (edge % 2) - 1
165 : line[(d + 2) % 3, :] = 2 * (edge // 2) - 1
166 : x, y, z = element.map(line)
167 : ax.plot(x, y, z, color="black")
168 :
169 : # Make plot square
170 : ax.set_aspect("equal")
171 : ```
172 :
173 :
174 :
175 : 
176 :
177 :
178 :
179 :
180 : ```python
181 :
182 : ```
|
