Data plotting
Contents
Data plotting#
This notebook is a tutorial on working with Tidy3D output data.
We will cover:
Accessing data.
Manipulating data.
Visualizing data.
First we import the packages we’ll need.
[1]:
%matplotlib inline
import numpy as np
import matplotlib.pylab as plt
import tidy3d as td
import tidy3d.web as web
[22:13:18] WARNING This version of Tidy3D was pip installed from the 'tidy3d-beta' repository on __init__.py:103 PyPI. Future releases will be uploaded to the 'tidy3d' repository. From now on, please use 'pip install tidy3d' instead.
INFO Using client version: 1.9.0rc1 __init__.py:121
Setup#
Creating Simulation#
First, let’s make a Simulation so we have data to plot.
We will add each possible type of monitor into the simultion to explore their output data separately.
[2]:
# simulation parameters
Lx, Ly, Lz = 5, 5, 5
min_steps_per_wvl = 32
# monitor parameters
freq0 = 2e14
freqs = np.linspace(1e14, 3e14, 11)
num_modes = 3
simulation = td.Simulation(
size=(Lx, Ly, Lz),
grid_spec=td.GridSpec.auto(min_steps_per_wvl=min_steps_per_wvl),
run_time=4e-13,
boundary_spec=td.BoundarySpec.all_sides(boundary=td.PML()),
structures=[
td.Structure(
geometry=td.Box(center=(0, 0, 0), size=(10001, 1.4, 1.5)),
medium=td.Medium(permittivity=2),
name="waveguide",
),
td.Structure(
geometry=td.Box(center=(0, 0.5, 0.5), size=(1.5, 1.4, 1.5)),
medium=td.Medium(permittivity=2),
name="scatterer",
),
],
sources=[
td.ModeSource(
source_time=td.GaussianPulse(freq0=freq0, fwidth=6e13),
center=(-2.0, 0.0, 0.0),
size=(0.0, 3, 3),
direction="+",
mode_spec=td.ModeSpec(),
mode_index=0,
)
],
monitors=[
td.FieldMonitor(
fields=["Ex", "Ey", "Ez"],
size=(td.inf, 0, td.inf),
center=(0, 0, 0),
freqs=freqs,
name="field",
),
td.FieldTimeMonitor(
fields=["Ex", "Ey", "Ez"],
size=(td.inf, 0, td.inf),
center=(0, 0, 0),
interval=200,
name="field_time",
),
td.FluxMonitor(size=(0, 3, 3), center=(2, 0, 0), freqs=freqs, name="flux"),
td.FluxTimeMonitor(
size=(0, 3, 3), center=(2, 0, 0), interval=10, name="flux_time"
),
td.ModeMonitor(
size=(0, 3, 3),
center=(2, 0, 0),
freqs=freqs,
mode_spec=td.ModeSpec(num_modes=num_modes),
name="mode",
),
],
)
[3]:
tmesh = simulation.tmesh
total_size_bytes = 0
for monitor in simulation.monitors:
monitor_grid = simulation.discretize(monitor)
num_cells = np.prod(monitor_grid.num_cells)
monitor_size = monitor.storage_size(num_cells=num_cells, tmesh=tmesh)
print(f"monitor {monitor.name} requires {monitor_size:.2e} bytes of storage.")
INFO Auto meshing using wavelength 1.4990 defined from sources. grid_spec.py:510
monitor field requires 7.06e+06 bytes of storage.
monitor field_time requires 1.06e+07 bytes of storage.
monitor flux requires 4.40e+01 bytes of storage.
monitor flux_time requires 2.65e+03 bytes of storage.
monitor mode requires 7.92e+02 bytes of storage.
Visualize Geometry#
We’ve created a simple waveguide with a defect / scattering region defined using a Box geometry.
A modal source is injected from the -x side of the simulation and we measure the mode amplitudes and flux at the +x side.
We’ve also placed a couple field monitors to visualize the field patterns.
Let’s take a look at the geometry from a few cross sections.
[4]:
_, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(14, 4))
simulation.plot(x=0.0, ax=ax1)
simulation.plot(y=0.01, ax=ax2)
simulation.plot(z=0.01, ax=ax3)
plt.show()
Make normalization Simulation#
For purposes of demonstration, let’s create another simulation without the scatterer, so we can compare what the data should look like for just the straight waveguide case.
This is as simple as making a copy of the original simulation and removing the scatterer from the list of structures.
[5]:
# get rid of scatterer for normalization
simulation0 = simulation.copy(update=dict(structures=[simulation.structures[0]]))
Running simulations#
Now we will run both simulations and load them into SimulationData objects.
Since these will be a short and simple simulations, we can use the web.run function to do them all in one line.
[6]:
sim0_data = web.run(
simulation0, task_name="straight waveguide", path="data/simulation.hdf5"
)
sim_data = web.run(
simulation, task_name="scattered waveguide", path="data/simulation.hdf5"
)
[22:13:19] INFO Auto meshing using wavelength 1.4990 defined from sources. grid_spec.py:510
Inspecting Log#
Now let’s take a look at the log to see if everything looks ok (fields decayed, etc).
[7]:
print(sim_data.log)
Simulation domain Nx, Ny, Nz: [176, 152, 152]
Applied symmetries: (0, 0, 0)
Number of computational grid points: 4.2214e+06.
Using subpixel averaging: True
Number of time steps: 6.6230e+03
Automatic shutoff factor: 1.00e-05
Time step (s): 6.0408e-17
Compute source modes time (s): 0.5539
Compute monitor modes time (s): 2.2783
Rest of setup time (s): 2.6144
Running solver for 6623 time steps...
- Time step 220 / time 1.33e-14s ( 3 % done), field decay: 1.00e+00
- Time step 264 / time 1.59e-14s ( 4 % done), field decay: 1.00e+00
- Time step 529 / time 3.20e-14s ( 8 % done), field decay: 1.00e+00
- Time step 794 / time 4.80e-14s ( 12 % done), field decay: 8.03e-04
- Time step 1059 / time 6.40e-14s ( 16 % done), field decay: 1.36e-06
Field decay smaller than shutoff factor, exiting solver.
Solver time (s): 2.7127
Post-Processing#
Now that the simulations have run and their data are loaded into the sim_data and sim0_data variables, we will explore how to access, manipulate, and visualize the data.
Accessing Data#
Original Simulation#
The SimulationData objects store a copy of the original Simulation, so it can be recovered if the SimulationData is loaded in a new session and the Simulation is no longer in memory.
[8]:
print(sim_data.simulation.size)
(5.0, 5.0, 5.0)
Monitor Data#
More importantly, the SimulationData contains a reference to the data for each of the monitors within the original Simulation. This data can be accessed directly using the name given to the monitors initially.
For example, our FluxMonitor was named 'flux' while our FluxTimeMonitor was named 'flux_time', each with a .flux attribute storing the raw data. Therefore, we can access the data as follows.
[9]:
# get the flux data from the monitor name
flux_data = sim_data["flux"].flux
flux_time_data = sim_data["flux_time"].flux
Structure of Data#
Now that we have the data loaded, let’s inspect it.
[10]:
flux_data
[10]:
<xarray.FluxDataArray (f: 11)>
array([0.591912 , 0.7745048 , 0.8931 , 0.95288116, 0.97218627,
0.9766828 , 0.9795849 , 0.9839063 , 0.9879946 , 0.9901106 ,
0.99019235], dtype=float32)
Coordinates:
* f (f) float64 1e+14 1.2e+14 1.4e+14 1.6e+14 ... 2.6e+14 2.8e+14 3e+14
Attributes:
units: W
long_name: fluxThe data is stored as a DataArray object using the xarray package. You can think of it as a dataset where the flux data is stored as a large, multi-dimensional array (like a numpy array) and the coordinates along each of the dimensions are specified so it is easy to work with.
For example:
[11]:
# flux values (W) as numpy array
print(f"shape of flux dataset = {flux_data.shape}\n.")
print(f"frequencies in dataset = {flux_data.coords.values} \n.")
print(f"flux values in dataset = {flux_data.values}\n.")
shape of flux dataset = (11,)
.
frequencies in dataset = <bound method Mapping.values of Coordinates:
* f (f) float64 1e+14 1.2e+14 1.4e+14 1.6e+14 ... 2.6e+14 2.8e+14 3e+14>
.
flux values in dataset = [0.591912 0.7745048 0.8931 0.95288116 0.97218627 0.9766828
0.9795849 0.9839063 0.9879946 0.9901106 0.99019235]
.
So we can see the datset stores a 1D array of flux values (Watts) and a 1D array of frequency values (Hz).
Manipulating Data#
We can now do some useful operations with this data, including:
[12]:
# Selecting data at a specific coordinate
flux_central_freq = flux_data.sel(f=freq0)
print(f"at central frequency, flux is {float(flux_central_freq):.2e} W")
at central frequency, flux is 9.77e-01 W
[13]:
# Selecting data at a specific index
flux_4th_freq = flux_data.isel(f=4)
print(f"at 4th frequency, flux is {float(flux_4th_freq):.2e} W")
at 4th frequency, flux is 9.72e-01 W
[14]:
# Interpolating data at a specific coordinate
flux_interp_freq = flux_data.interp(f=213.243523142e12)
print(
f"at an intermediate (not stored) frequency, flux is {float(flux_interp_freq):.2e} W"
)
at an intermediate (not stored) frequency, flux is 9.79e-01 W
[15]:
# arithmetic opertions with numbers
flux_times_2_plus_1 = 2 * flux_data + 1
[16]:
# operations with dataset (in this case, dividing by the normlization flux)
flux_transmission = flux_data / sim0_data["flux"].flux
There are mny more things you can do with the xarray package that are not covered here, so it is best to take a look at their documention for more details.
Plotting Data#
These datasets have built in plotting methods, which can be handy for quickly producing customizable plots.
Let’s plot the flux data to get a feeling.
[17]:
f, (ax1, ax2, ax3) = plt.subplots(1, 3, tight_layout=True, figsize=(10, 3))
flux_data.plot(ax=ax1)
ax1.set_title("flux vs. frequency")
ax1.set_xlim(150e12, 250e12)
# plot the ratio of flux with scatterer to without scatterer with dashed line signifying unity transmission.
flux_transmission.plot(ax=ax2, color="k")
ax2.plot(flux_data.f, np.ones_like(flux_data.f), "--")
ax2.set_ylabel("flux(f) / flux_normalize(f)")
ax2.set_xlabel("frequency (Hz)")
ax2.set_title("transmission spectrum")
ax2.set_xlim(150e12, 250e12)
ax2.set_ylim(-0.1, 1.1)
# plot the time dependence of the flux
flux_time_data.plot(color="crimson", ax=ax3)
ax3.set_title("flux vs. time")
plt.show()
Conveniently, xarray adds a lot of the plot labels and formatting automatically, which can save some time.
Complex, Multi-dimensional Data#
While the flux data is 1D and quite simple, the same approaches apply to more complex data, such as mode amplitude and field data, which can be complex-valued and muli-dimensional.
Let’s use the ModeMonitor data as an example. We’ll access the data using it’s name 'mode' in the original simulation and print out some of the metadata to examine.
[18]:
mode_data = sim_data["mode"]
mode_data.amps.shape
[18]:
(2, 11, 3)
As we can see, the data is complex-valued and contains three dimensions. - a propagation direction for the mode, which is just positive or negative (+/-) - an index into the list of modes returned by the solver given our initial ModeSpec specifications. - the frequency of the mode
Let’s plot the data just to get a feeling for it, using .sel() to split the forward and backward propagating modes onto separate plot axes.
[19]:
f, (ax1, ax2) = plt.subplots(1, 2, tight_layout=True, figsize=(8, 3))
mode_data.amps.sel(direction="+").abs.plot.line(x="f", ax=ax1)
mode_data.amps.sel(direction="-").abs.plot.line(x="f", ax=ax2)
ax1.set_xlim(150e12, 250e12)
ax2.set_xlim(150e12, 250e12)
plt.show()
As before, we can manipulate the data using basic algebra or xarray built ins, for example.
[20]:
# sum abolute value squared of the mode amplitudes to get powers
mode_powers = abs(mode_data.amps) ** 2
# select the powers at the central frequency only
powers_central_freq = mode_powers.sel(f=freq0)
# sum the powers over all of the mode indices
powers_sum_modes = powers_central_freq.sum("mode_index")
Field Data#
Data from a FieldMonitor or FieldTimeMonitor is more complicated as it contains data from several field components.
As such, when we just grab the data by name, we dont get an xarray object, but rather a FieldData container holding all of the field components.
[21]:
field_data = sim_data["field"]
field0_data = sim0_data["field"]
print(type(field_data))
<class 'tidy3d.components.data.monitor_data.FieldData'>
The field_data object contains data for each of the field components in it’s data_dict dictionary.
Let’s look at what field components are contined
[22]:
print(field_data.dict().keys())
dict_keys(['type', 'Ex', 'Ey', 'Ez', 'Hx', 'Hy', 'Hz', 'monitor', 'symmetry', 'symmetry_center', 'grid_expanded', 'grid_primal_correction', 'grid_dual_correction'])
The individual field components, themselves, can be accessed conveniently using a “dot” syntax, and are stored similarly to the flux and mode data as xarray data arrays.
[23]:
field_data.Ex
[23]:
<xarray.ScalarFieldDataArray (x: 178, y: 1, z: 154, f: 11)>
array([[[[ 0.00000000e+00+0.0000000e+00j,
0.00000000e+00+0.0000000e+00j,
-0.00000000e+00-0.0000000e+00j, ...,
0.00000000e+00+0.0000000e+00j,
0.00000000e+00-0.0000000e+00j,
-0.00000000e+00-0.0000000e+00j],
[ 0.00000000e+00+0.0000000e+00j,
0.00000000e+00+0.0000000e+00j,
-0.00000000e+00-0.0000000e+00j, ...,
0.00000000e+00+0.0000000e+00j,
0.00000000e+00-0.0000000e+00j,
-0.00000000e+00-0.0000000e+00j],
[ 6.43478160e-08+2.2920256e-08j,
-7.46410151e-08+4.1726917e-08j,
-3.07212815e-08-2.5303006e-09j, ...,
2.56144510e-08+2.5176810e-08j,
-3.49830742e-08+3.8806519e-08j,
-1.45186245e-08-1.9318170e-08j],
...,
[-5.54970967e-08-1.4926345e-08j,
...
-2.41072717e-09+5.2363074e-09j],
...,
[-6.27467078e-10-7.4389943e-09j,
1.45480961e-09+9.0445926e-09j,
2.89787055e-10-5.5823879e-09j, ...,
1.90179317e-09+1.7743294e-09j,
-9.13492282e-10-2.1816593e-09j,
4.50721593e-09-6.2089284e-10j],
[-6.59422794e-10+8.6273655e-10j,
7.13056489e-11+9.2542590e-10j,
-2.19914156e-11-5.1966337e-11j, ...,
8.63919630e-11+4.0660250e-10j,
-1.42005477e-10+1.8165988e-10j,
4.32076513e-10+1.3890591e-10j],
[ 0.00000000e+00+0.0000000e+00j,
0.00000000e+00+0.0000000e+00j,
-0.00000000e+00-0.0000000e+00j, ...,
0.00000000e+00+0.0000000e+00j,
0.00000000e+00-0.0000000e+00j,
-0.00000000e+00-0.0000000e+00j]]]], dtype=complex64)
Coordinates:
* x (x) float64 -2.913 -2.88 -2.847 -2.814 ... 2.814 2.847 2.88 2.913
* y (y) float64 0.0
* z (z) float64 -3.108 -3.061 -3.015 -2.968 ... 2.922 2.968 3.015 3.062
* f (f) float64 1e+14 1.2e+14 1.4e+14 1.6e+14 ... 2.6e+14 2.8e+14 3e+14
Attributes:
long_name: field valueIn this case, the dimensions of the data are the spatial locations on the yee lattice (x,y,z) and the frequency.
Centering Field Data#
For many advanced plots and other manipulations, it is convenient to have the field data co-located at the same positions, since they are naturally defined on separate locations on the yee lattice. For this, we provide an .at_centers() method of SimulationData to return the fields co-located at the yee grid centers.
[24]:
field0_data_centered = sim0_data.at_centers("field").interp(f=200e12)
field_data_centered = sim_data.at_centers("field").interp(f=200e12)
[22:14:20] INFO Auto meshing using wavelength 1.4990 defined from sources. grid_spec.py:510
INFO Auto meshing using wavelength 1.4990 defined from sources. grid_spec.py:510
The centered field data is stored as an xarray Dataset, which provies similar functionality as the DataArray objects, but is aware of all of the field data and provides more convenience methods, as we will explore in the next section.
Plotting Fields#
Simple#
Plotting fields requires a little bit of manipulation to get them in the right shape, but then it is straightforward.
[25]:
# get the field data on the y=0 plane at frequency 200THz
Ez_data = field_data.Ez.isel(y=0).interp(f=2e14)
Ez0_data = field0_data.Ez.isel(y=0).interp(f=2e14)
[26]:
# amplitude plot rz Ex(x,y) field on plane
_, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, tight_layout=True, figsize=(10, 6))
Ez0_data.real.plot(x="x", y="z", ax=ax1)
Ez_data.real.plot(x="x", y="z", ax=ax2)
abs(Ez0_data).plot(x="x", y="z", cmap="magma", ax=ax3)
abs(Ez_data).plot(x="x", y="z", cmap="magma", ax=ax4)
ax1.set_title("without scatterer (real Ez)")
ax2.set_title("with scatterer (real Ez)")
ax3.set_title("without scatterer (abs Ez)")
ax4.set_title("with scatterer (abs Ez)")
plt.show()
Advanced#
More advanced plotting is available and is best done with the centered field data.
Structure Overlay#
One can overlay the structure permittivity by calling plot_fields from the SimulationData object as follows:
[27]:
ax = sim_data.plot_field("field", "Ex", y=0.0, f=200e12, eps_alpha=0.2)
Quiver Plots#
Quiver plotting can be done through the centered field Dataset. We first should spatially downsample the data to be quivered to reduce clutter on the plot, which can be done by .sel() with slice(None, None, skip) range to select every skip data point along that axis.
[28]:
# downsample the field data for more clear quiver plotting
field0_data_resampled = field0_data_centered.sel(
x=slice(None, None, 7), z=slice(None, None, 7)
)
# quiver plot of \vec{E}_{x,y}(x,y) on plane with Ez(x,y) underlying.
f, ax = plt.subplots(figsize=(8, 5))
field0_data_centered.isel(y=0).Ez.real.plot.imshow(x="x", y="z", ax=ax, robust=True)
field0_data_resampled.isel(y=0).real.plot.quiver("x", "z", "Ex", "Ez", ax=ax)
plt.show()
Stream Plots#
Stream plots can be another useful feature, although they do not take the amplitude into account, so use with caution.
Time Lapse Fields#
For time-domain data, it can be convenient to see the fields at various snapshots of the simulation. For this, it is convenient to supply the row and col keyword arguments to plot, which expand all plots on a row or column.
[29]:
field_time_data = sim_data["field_time"]
# select the data and downsample and window the time coordinate to reduce number of plots.
time_data = field_time_data.Ez.isel(y=0).sel(t=slice(1e-14, 5e-14, 1))
time_data.plot(x="x", y="z", col="t")
plt.show()
Export data to csv#
You can also export xarray.DataArray to csv file using its to_dataframe() method.
[30]:
df = time_data.to_dataframe("value")
df.to_csv("./data/time_data.csv")
It will stack the array coordinates to generate a 2D table. Import the csv file to check the content:
[31]:
import pandas as pd
df = pd.read_csv("./data/time_data.csv")
display(df)
| x | z | t | y | value | |
|---|---|---|---|---|---|
| 0 | -2.929245 | -3.084856 | 1.208153e-14 | 0.0 | 0.0 |
| 1 | -2.929245 | -3.084856 | 2.416306e-14 | 0.0 | 0.0 |
| 2 | -2.929245 | -3.084856 | 3.624459e-14 | 0.0 | 0.0 |
| 3 | -2.929245 | -3.084856 | 4.832612e-14 | 0.0 | 0.0 |
| 4 | -2.929245 | -3.038042 | 1.208153e-14 | 0.0 | 0.0 |
| ... | ... | ... | ... | ... | ... |
| 109643 | 2.896226 | 3.038690 | 4.832612e-14 | 0.0 | 0.0 |
| 109644 | 2.896226 | 3.085504 | 1.208153e-14 | 0.0 | 0.0 |
| 109645 | 2.896226 | 3.085504 | 2.416306e-14 | 0.0 | 0.0 |
| 109646 | 2.896226 | 3.085504 | 3.624459e-14 | 0.0 | 0.0 |
| 109647 | 2.896226 | 3.085504 | 4.832612e-14 | 0.0 | 0.0 |
109648 rows × 5 columns
Conclusion#
This hopefully gives some sense of what kind of data post processing and visualization can be done in Tidy3D.
For more detailed information and reference, the best places to check out are
API reference for SimulationData.