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.IO.H5 import list_observations, open_volfiles
40 : import glob
41 :
42 : h5files = glob.glob(
43 : # Point this to your volume data files
44 : "BbhVolume*.h5"
45 : )
46 : subfile_name = "/VolumeData"
47 : obs_ids, obs_times = list_observations(open_volfiles(h5files, subfile_name))
48 : ```
49 :
50 : We can also list all available variables in the volume data files:
51 :
52 :
53 : ```python
54 : import rich.columns
55 :
56 : for volfile in open_volfiles(h5files, subfile_name, obs_ids[-1]):
57 : all_vars = volfile.list_tensor_components(obs_ids[-1])
58 : # Only look into the first file. The other should have the same variables.
59 : break
60 :
61 : rich.print(rich.columns.Columns(all_vars))
62 : ```
63 :
64 :
65 : <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
66 : ExtrinsicCurvature_yy ExtrinsicCurvature_zx ExtrinsicCurvature_zy
67 : ExtrinsicCurvature_zz InertialCoordinates_x InertialCoordinates_y
68 : InertialCoordinates_z Lapse RadiallyCompressedCoordinates_x
69 : RadiallyCompressedCoordinates_y RadiallyCompressedCoordinates_z ShiftExcess_x
70 : ShiftExcess_y ShiftExcess_z Shift_x
71 : Shift_y Shift_z SpatialMetric_xx
72 : SpatialMetric_yx SpatialMetric_yy SpatialMetric_zx
73 : SpatialMetric_zy SpatialMetric_zz
74 : </pre>
75 :
76 :
77 :
78 : ### Interpolate to a slice
79 :
80 : One way of visualizing your 3D volume data in a 2D plot is with a slice. Use
81 : `interpolate_to_points` to interpolate your volume data to any set of
82 : coordinates:
83 :
84 :
85 : ```python
86 : from spectre.IO.Exporter import interpolate_to_points
87 :
88 : # Coordinates of a regular grid in the xy plane
89 : x, y = np.meshgrid(np.linspace(-15, 15, 300), np.linspace(-15, 15, 300))
90 : z = np.zeros(x.shape)
91 :
92 : # Interpolation can be slow, so it's useful to cache the result. Remember to
93 : # delete the cache file if you change anything.
94 : cache_file = "interpolated_data.npy"
95 : if os.path.isfile(cache_file):
96 : lapse = np.load(cache_file)
97 : else:
98 : lapse = np.array(
99 : interpolate_to_points(
100 : h5files,
101 : subfile_name=subfile_name,
102 : observation_id=obs_ids[-1],
103 : tensor_components=["Lapse"],
104 : target_points=[x.flatten(), y.flatten(), z.flatten()],
105 : )
106 : ).reshape(x.shape)
107 : np.save(cache_file, lapse)
108 : ```
109 :
110 : Now we can plot the 2D slice with `matplotlib`:
111 :
112 :
113 : ```python
114 : # Plot lapse
115 : plt.contourf(x, y, np.log(1 - lapse))
116 :
117 : # Plot circles for black holes
118 : ax = plt.gca()
119 : for bh_pos in [-8, 8]:
120 : ax.add_patch(
121 : plt.Circle(xy=(bh_pos, 0), radius=0.89, color="black", fill=True)
122 : )
123 :
124 : # Make plot square
125 : ax.set_aspect("equal")
126 : ```
127 :
128 : /var/folders/xp/5t2ny359187c4ljckf8sz3m80000gn/T/ipykernel_19858/3375439396.py:2: RuntimeWarning: invalid value encountered in subtract
129 : plt.contourf(x, y, np.log(1 - lapse))
130 :
131 :
132 :
133 :
134 : 
135 :
136 :
137 :
138 : ### Plot individual elements
139 :
140 : You can also iterate over elements in your volume data to plot things like
141 : element boundaries, collocation points, etc. Use `iter_elements`:
142 :
143 :
144 : ```python
145 : from spectre.IO.H5.IterElements import iter_elements
146 :
147 : # Create a 3D plot
148 : ax = plt.gcf().add_subplot(111, projection="3d")
149 : ax.axis("off")
150 :
151 : # Iterate over elements
152 : for element in iter_elements(
153 : open_volfiles(h5files, subfile_name),
154 : obs_ids=obs_ids[0],
155 : element_patterns=["B3,*"], # Only plot elements in block 3
156 : ):
157 : # Draw outline of the element by mapping the edges of the logical cube to
158 : # inertial coordinates using the element map
159 : for d in range(3):
160 : for edge in range(4):
161 : line = np.zeros((3, 100))
162 : line[d, :] = np.linspace(-1, 1, 100)
163 : line[(d + 1) % 3, :] = 2 * (edge % 2) - 1
164 : line[(d + 2) % 3, :] = 2 * (edge // 2) - 1
165 : x, y, z = element.map(line)
166 : ax.plot(x, y, z, color="black")
167 :
168 : # Make plot square
169 : ax.set_aspect("equal")
170 : ```
171 :
172 :
173 :
174 : 
175 :
176 :
177 :
178 : ### Transforming volume data
179 :
180 : Often you want to post-process volume data to compute derived quantities from
181 : the ones the simulation has written out. You can sometimes do this within tools
182 : like ParaView (e.g. using ParaView's "Calculator" filter) if the computation is
183 : simple enough and pointwise (i.e., needs no derivatives or other mesh
184 : information). If you can't get what you need in ParaView, you can use `spectre
185 : transform-vol`. It takes any Python function (a "kernel"), runs it over your
186 : volume data, and writes the result back into the files.
187 :
188 : For example, let's add the number of grid points in each element as a field
189 : that we can visualize:
190 :
191 : ```sh
192 : spectre transform-vol BbhVolume*.h5 -d VolumeData \
193 : -k spectre.Spectral.Mesh3D:extents
194 : ```
195 :
196 : You will be prompted to select an output dataset name (hit enter to select the
197 : default "Extents"). The result will be written back into the H5 files. If you
198 : now regenerate an XDMF file for the volume data you will be able to see the
199 : output in ParaView.
200 :
201 : The `transform-vol` tool supports arbitrary Python functions as kernels. Many
202 : Python functions from the `spectre` module are useful kernels. For example,
203 : here are some useful kernels:
204 :
205 : - p-refinement (number of grid points per dimension):
206 : `spectre.Spectral.Mesh3D:extents`
207 : - Relative/absolute truncation error (from power monitors):
208 : `spectre.NumericalAlgorithms.LinearOperators:relative_truncation_error`
209 : or `absolute_truncation_error`
210 : - Coordinates in different frames: `Element:grid_coordinates` or
211 : `Element:distorted_coordinates`
212 :
213 : You can find more useful kernels by browsing the `spectre.PointwiseFunctions`
214 : module: https://spectre-code.org/py/_autosummary/spectre.PointwiseFunctions.html
215 :
216 : You can also write your own kernels. Create a Python file, e.g. `kernel.py`, and
217 : write a function like this:
218 :
219 :
220 : ```python
221 : %%file kernel.py
222 : from spectre.DataStructures.Tensor import Scalar, DataVector
223 :
224 : def lapse_squared(lapse: Scalar[DataVector]) -> Scalar[DataVector]:
225 : return np.array(lapse)**2
226 : ```
227 :
228 : Now you can run this kernel over your volume data like this:
229 :
230 : ```sh
231 : spectre transform-vol BbhVolume*.h5 -d VolumeData \
232 : -e kernel.py -k lapse_squared
233 : ```
234 :
235 : You can find more details by running `spectre transform-vol --help`.
236 :
237 :
|
