Performing visualization of simulation data#

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]:

import numpy as np
import matplotlib.pylab as plt

import tidy3d as td
import tidy3d.web as web

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",
        ),
    ],
)

[15:20:57] WARNING: Default value for the field monitor           monitor.py:261
           'colocate' setting has changed to 'True' in Tidy3D                   
           2.4.0. All field components will be colocated to the                 
           grid boundaries. Set to 'False' to get the raw fields                
           on the Yee grid instead.

           WARNING: Default value for the field monitor           monitor.py:261
           'colocate' setting has changed to 'True' in Tidy3D                   
           2.4.0. All field components will be colocated to the                 
           grid boundaries. Set to 'False' to get the raw fields                
           on the Yee grid instead.

[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.")

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/simulation0.hdf5",
    verbose=True,
)

sim_data = web.run(
    simulation,
    task_name="scattered waveguide",
    path="data/simulation.hdf5",
    verbose=True,
)

[15:20:58] Created task 'straight waveguide' with task_id          webapi.py:188
           'fdve-a81e80bb-710f-4f9c-80e0-95da36d80335v1'.

           View task using web UI at                               webapi.py:190
           'https://tidy3d.simulation.cloud/workbench?taskId=fdve-              
           a81e80bb-710f-4f9c-80e0-95da36d80335v1'.

[15:20:59] status = queued                                         webapi.py:361

[15:21:03] status = preprocess                                     webapi.py:355

[15:21:08] Maximum FlexCredit cost: 0.025. Use                     webapi.py:341
           'web.real_cost(task_id)' to get the billed FlexCredit                
           cost after a simulation run.

           starting up solver                                      webapi.py:377

           running solver                                          webapi.py:386

           To cancel the simulation, use 'web.abort(task_id)' or   webapi.py:387
           'web.delete(task_id)' or abort/delete the task in the                
           web UI. Terminating the Python script will not stop the              
           job running on the cloud.

[15:21:14] early shutoff detected, exiting.                        webapi.py:404

           status = postprocess                                    webapi.py:420

[15:21:35] status = success                                        webapi.py:427

[15:21:36] loading SimulationData from data/simulation0.hdf5       webapi.py:591

           Created task 'scattered waveguide' with task_id         webapi.py:188
           'fdve-d66a2a3e-747f-4c19-a006-a44490a50d38v1'.

           View task using web UI at                               webapi.py:190
           'https://tidy3d.simulation.cloud/workbench?taskId=fdve-              
           d66a2a3e-747f-4c19-a006-a44490a50d38v1'.

[15:21:37] status = queued                                         webapi.py:361

[15:21:42] status = preprocess                                     webapi.py:355

[15:21:46] Maximum FlexCredit cost: 0.025. Use                     webapi.py:341
           'web.real_cost(task_id)' to get the billed FlexCredit                
           cost after a simulation run.

           starting up solver                                      webapi.py:377

           running solver                                          webapi.py:386

           To cancel the simulation, use 'web.abort(task_id)' or   webapi.py:387
           'web.delete(task_id)' or abort/delete the task in the                
           web UI. Terminating the Python script will not stop the              
           job running on the cloud.

[15:21:53] early shutoff detected, exiting.                        webapi.py:404

           status = postprocess                                    webapi.py:420

[15:22:13] status = success                                        webapi.py:427

[15:22:14] loading SimulationData from data/simulation.hdf5        webapi.py:591

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.4482
Compute monitor modes time (s):    2.5589
Rest of setup time (s):            2.3952

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):                   0.9497
Data write time (s):               0.0540

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.5919975 , 0.7746099 , 0.8931861 , 0.9529229 , 0.97218525,
       0.9766577 , 0.97956336, 0.98391604, 0.9880343 , 0.99017584,
       0.9902634 ], 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:  flux

The 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.5919975  0.7746099  0.8931861  0.9529229  0.97218525 0.9766577
 0.97956336 0.98391604 0.9880343  0.99017584 0.9902634 ]
.

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")

Export data to csv#

In many cases, you might want to save the monitor data as a .csv file for future use. This can be done in multiple ways. Here we demonstrate two methods as examples.

For saving a 2D array such as a field plot at a specific frequency or time instance, we can first convert the xarray.DataArray to a numpy.array using the to_numpy() method and then using numpy’s savetxt() method to save the 2D numpy array as a .csv file. The spatial coordinates can be saved as separate .csv files in the same way.

[31]:

# save Ez to a numpy array
Ez = Ez0_data.real.squeeze().to_numpy()

# save x coordinate to a numpy array
x_array = Ez0_data.real.coords['x'].to_numpy()

# save z coordinate to a numpy array
z_array =  Ez0_data.real.coords['z'].to_numpy()

# save the numpy arries to .csv files
np.savetxt('./data/Ez.csv', Ez, delimiter=',')
np.savetxt('./data/x.csv', x_array)
np.savetxt('./data/z.csv', z_array)

If the data you want to save is higher dimensional such as the field distribution at multiple time instances, you can consider export the xarray.DataArray to a csv file using its to_dataframe() method.

[32]:

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:

[33]:

import pandas as pd

df = pd.read_csv("./data/time_data.csv")
display(df)

	x	z	t	y	value
0	-2.896226	-3.061435	1.208153e-14	0.0	0.0
1	-2.896226	-3.061435	2.416306e-14	0.0	0.0
2	-2.896226	-3.061435	3.624459e-14	0.0	0.0
3	-2.896226	-3.061435	4.832612e-14	0.0	0.0
4	-2.896226	-3.014649	1.208153e-14	0.0	0.0
...	...	...	...	...	...
108319	2.896226	3.015268	4.832612e-14	0.0	0.0
108320	2.896226	3.062111	1.208153e-14	0.0	0.0
108321	2.896226	3.062111	2.416306e-14	0.0	0.0
108322	2.896226	3.062111	3.624459e-14	0.0	0.0
108323	2.896226	3.062111	4.832612e-14	0.0	0.0

108324 rows × 5 columns

Load data#

When running the simulations, we have specified the path to save the simulation data, i.e. path="data/simulation0.hdf5". It means the simulation data is saved as a .hdf5 file, which can be loaded into a SimulationData object in the future for analysis. Here we demonstrate loading the .hdf5 files using the from_file method. If you performed a mode solving simulation instead of a FDTD simulation, ModeSolverData can be loaded in a similar fashion.

[34]:

sim_data0 = td.SimulationData.from_file("data/simulation0.hdf5")
sim_data = td.SimulationData.from_file("data/simulation.hdf5")

After the simulation data is loaded, you can perform postprocessing as discussed above such as plotting the field distribution by using the plot_field method.

[35]:

fig, (ax1, ax2) = plt.subplots(1,2,tight_layout=True, figsize=(8, 4))
sim_data0.plot_field(field_monitor_name="field", field_name='E', val='abs', f=2e14, ax=ax1)
sim_data.plot_field(field_monitor_name="field", field_name='E', val='abs', f=2e14, ax=ax2)
plt.show()

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.
Xarray documentation.

If you are new to the finite-difference time-domain (FDTD) method, we highly recommend going through our FDTD101 tutorials.

[ ]:

Performing visualization of simulation data

Contents

Performing visualization of simulation data#

Setup#

Creating Simulation#

Visualize Geometry#

Make normalization Simulation#

Running simulations#

Inspecting Log#

Post-Processing#

Accessing Data#

Original Simulation#

Monitor Data#

Structure of Data#

Manipulating Data#

Plotting Data#

Complex, Multi-dimensional Data#

Field Data#

Colocating Field Data#

Plotting Fields#

Simple#

Advanced#

Structure Overlay#

Quiver Plots#

Stream Plots#

Time Lapse Fields#

Export data to csv#

Load data#

Conclusion#