Advanced monitor data (xarray) manipulation and visualization#

In the previous tutorial, we have demonstrated the basics of data manipulation and visualization. As discussed, in Tidy3D, we represent multi-dimensional data (such as arrays) using xarray.DataArray objects. It’s important to note that the xarray package provides a lot of functionalities that can greatly simplify our data postprocessing and visualization when used appropriately.

In this tutorial, we will first briefly review some basic xarray features and then delve into some advanced features such as performing interpolation and integration. This tutorial is designed for users who are already familiar with Tidy3D’s basic workflow and would like to learn more about data manipulation and visualization. If you are very new to Tidy3D, we suggest checking out the quick start, the first walkthrough, and the working principles tutorials first.

It’s important to note that the xarray library contains many more functionalities that are not demonstrated here. We encourage users to reference the documentation to see what is possible. We also recommend users to watch the introduction video on xarray if you would like to learn more about it. After overcoming the initial learning curve, our users all find data manipulation with xarray to be extremely flexible and efficient.

[1]:

import numpy as np
import matplotlib.pyplot as plt
import os

import tidy3d as td
from tidy3d import web

Load Simulation Data#

Since this tutorial is focused on data analysis, we will skip the simulation setup part and only load the simulation result from a previously obtained hdf5 file. We load the hdf5 file to a variable sim_data using the from_hdf5 method of SimulationData. This is completely analogous to the situation where you just finished running a simulation.

[2]:

sim_data = td.SimulationData.from_hdf5("misc/mmi_simulation_data.h5")

09:48:35 Eastern Daylight Time WARNING: updating Simulation from 2.4 to 2.5

The loaded data contains two parts: simulation and monitor_data. simulation is simply the original Simulation. We can use it to visualize the simulation in 3D for example.

As we can see, this simulation contains a 1 by 2 MMI power splitter.

[3]:

sim_data.simulation.plot_3d()

monitor_data, as the name suggests, contains data from all monitors defined in the simulation. monitor_data is a dictionary with the keys being the names of the monitors. Let’s print out the keys to see what they are.

Here we have three monitors named Flux monitor, Mode monitor, and Field monitor. In the following sections, we will examine them individually and use them to demonstrate various functionalities of xarray.

[4]:

sim_data.monitor_data.keys()

[4]:

dict_keys(['Flux monitor', 'Mode monitor', 'Field monitor'])

FluxMonitor Data#

We can access the data for each monitor by indexing into the SimulationData with the monitor .name. The monitor dataset contains one or more xarray.DataArray objects storing the actual raw data. For the flux monitor data, we can access the raw flux data as a function of frequency with .flux. As most data are multidimensional, it’s often very helpful to print out the data and directly inspect its structure.

Here we can see that the flux is an array with only one coordinate f, which is frequency in Hz. The unit of flux is W.

[5]:

sim_data["Flux monitor"].flux

[5]:

<xarray.FluxDataArray (f: 50)>
array([0.49619535, 0.4959998 , 0.4958456 , 0.49578664, 0.4958678 ,
       0.49604478, 0.49626908, 0.49647027, 0.4966077 , 0.49661922,
       0.49642277, 0.49594647, 0.4951819 , 0.4942104 , 0.4932421 ,
       0.4924945 , 0.4921227 , 0.49211392, 0.49230942, 0.4924388 ,
       0.49224442, 0.49157965, 0.4904679 , 0.48906133, 0.4876103 ,
       0.48636082, 0.4854779 , 0.4849979 , 0.48482403, 0.48475057,
       0.48453504, 0.48398182, 0.48302355, 0.48175392, 0.48037806,
       0.47914445, 0.47820777, 0.47757253, 0.47708204, 0.47651216,
       0.47562975, 0.47437435, 0.4728097 , 0.47111368, 0.4695425 ,
       0.46828848, 0.4674291 , 0.4669524 , 0.46673352, 0.46655577],
      dtype=float32)
Coordinates:
  * f        (f) float64 1.999e+14 1.996e+14 1.993e+14 ... 1.876e+14 1.874e+14
Attributes:
    units:      W
    long_name:  flux

xarray.FluxDataArray

f: 50

0.4962 0.496 0.4958 0.4958 0.4959 ... 0.4674 0.467 0.4667 0.4666

array([0.49619535, 0.4959998 , 0.4958456 , 0.49578664, 0.4958678 ,
       0.49604478, 0.49626908, 0.49647027, 0.4966077 , 0.49661922,
       0.49642277, 0.49594647, 0.4951819 , 0.4942104 , 0.4932421 ,
       0.4924945 , 0.4921227 , 0.49211392, 0.49230942, 0.4924388 ,
       0.49224442, 0.49157965, 0.4904679 , 0.48906133, 0.4876103 ,
       0.48636082, 0.4854779 , 0.4849979 , 0.48482403, 0.48475057,
       0.48453504, 0.48398182, 0.48302355, 0.48175392, 0.48037806,
       0.47914445, 0.47820777, 0.47757253, 0.47708204, 0.47651216,
       0.47562975, 0.47437435, 0.4728097 , 0.47111368, 0.4695425 ,
       0.46828848, 0.4674291 , 0.4669524 , 0.46673352, 0.46655577],
      dtype=float32)

Coordinates: (1)

(f)

float64

1.999e+14 1.996e+14 ... 1.874e+14

units :: Hz
long_name :: frequency

array([1.998616e+14, 1.995901e+14, 1.993193e+14, 1.990492e+14, 1.987798e+14,
       1.985112e+14, 1.982433e+14, 1.979762e+14, 1.977097e+14, 1.974440e+14,
       1.971789e+14, 1.969146e+14, 1.966510e+14, 1.963881e+14, 1.961259e+14,
       1.958644e+14, 1.956036e+14, 1.953435e+14, 1.950841e+14, 1.948253e+14,
       1.945673e+14, 1.943099e+14, 1.940532e+14, 1.937972e+14, 1.935419e+14,
       1.932872e+14, 1.930333e+14, 1.927799e+14, 1.925273e+14, 1.922753e+14,
       1.920239e+14, 1.917732e+14, 1.915232e+14, 1.912738e+14, 1.910251e+14,
       1.907770e+14, 1.905296e+14, 1.902828e+14, 1.900366e+14, 1.897911e+14,
       1.895462e+14, 1.893019e+14, 1.890583e+14, 1.888153e+14, 1.885729e+14,
       1.883312e+14, 1.880900e+14, 1.878495e+14, 1.876096e+14, 1.873703e+14])

Indexes: (1)

PandasIndex

PandasIndex(Index([199861643525342.56, 199590092379248.34,  199319278142641.5,
        199049197819955.0, 198779848431835.97,  198511227015036.2,
       198243330622303.34, 197976156322273.28,  197709701199363.1,
       197443962353665.03,  197178936900841.3,  196914621972019.8,
       196651014713690.44, 196388112287602.62, 196125911870663.25,
        195864410654835.7, 195603605847039.62, 195343494669051.53,
       195084074357406.06,  194825342163298.1, 194567295352485.75,
       194309931205194.12, 194053247016019.47, 193797240093834.75,
        193541907761695.3, 193287247356745.72, 193033256230127.12,
       192779931746885.47, 192527271285880.38,  192275272239694.7,
       192023932014544.75, 191773248030191.56,  191523217719852.3,
        191273838530112.9,  191025107920841.0, 190777023365099.62,
       190529582349061.88, 190282782371925.78, 190036620945830.16,
       189791095595770.94, 189546203859518.34,  189301943287534.4,
        189058311442891.5, 188815305901191.12, 188572924250483.56,
       188331164091188.06,  188090023036013.7,  187849498709880.7,
       187609588749842.53, 187370290805008.84],
      dtype='float64', name='f'))

Attributes: (2)
units :
W
long_name :
flux

The data array has .values and .coords. The .values are just the raw, unlabelled data. The .coords are a dictionary mapping each dimension of the data to a set of data. For example, mapping the flux values to a specific frequency. We can inspect them by printing them out.

[6]:

sim_data["Flux monitor"].flux.values

[6]:

array([0.49619535, 0.4959998 , 0.4958456 , 0.49578664, 0.4958678 ,
       0.49604478, 0.49626908, 0.49647027, 0.4966077 , 0.49661922,
       0.49642277, 0.49594647, 0.4951819 , 0.4942104 , 0.4932421 ,
       0.4924945 , 0.4921227 , 0.49211392, 0.49230942, 0.4924388 ,
       0.49224442, 0.49157965, 0.4904679 , 0.48906133, 0.4876103 ,
       0.48636082, 0.4854779 , 0.4849979 , 0.48482403, 0.48475057,
       0.48453504, 0.48398182, 0.48302355, 0.48175392, 0.48037806,
       0.47914445, 0.47820777, 0.47757253, 0.47708204, 0.47651216,
       0.47562975, 0.47437435, 0.4728097 , 0.47111368, 0.4695425 ,
       0.46828848, 0.4674291 , 0.4669524 , 0.46673352, 0.46655577],
      dtype=float32)

[7]:

sim_data["Flux monitor"].flux.coords

[7]:

Coordinates:
  * f        (f) float64 1.999e+14 1.996e+14 1.993e+14 ... 1.876e+14 1.874e+14

xarray allows us to quickly plot the data. In this case, the plot will just be a 1D line since the data only contains one coordinate. Note that the x and y axis labels are automatically created accordingly.

[8]:

sim_data["Flux monitor"].flux.plot()
plt.show()

../_images/notebooks_XarrayTutorial_17_0.png

In many cases, we would like to customize the plot. For example, the transmission is often plotted in dB scale. Instead of frequency in Hz as the x axis, wavelength in micron is more standard. In addition, we might want to plot the curve in a certain style. This can be achieved by first extracting the data from the DataArray and manually plotting them using matplotlib or other visualization library of choice.

[9]:

# extract frequecies from the monitor data
freqs = sim_data["Flux monitor"].flux.f

# extract the transmission flux from the monitor data
T = sim_data["Flux monitor"].flux

# plot the transmission in dB
plt.plot(td.C_0 / freqs, 10 * np.log10(T), "red", linewidth=3)
plt.xlabel("Wavelength ($\mu m$)")
plt.ylabel("Transmission (dB)")
plt.ylim(-10, 0)
plt.show()

../_images/notebooks_XarrayTutorial_19_0.png

Saving Data as a CSV File#

After running a simulation, often we want to save the data for future use to avoid rerunning the same simulation. In Tidy3D, after a simulation is finished, the hdf5 data file will be automatically downloaded to the specified path. We can load it in the future as demonstrated at the beginning of this notebook.

Sometimes users might prefer to only save certain data in other formats such as csv files. Here we demonstrate how to save the transmission flux spectrum. First we define a numpy array to store the data and then use the savetxt method to save it as a csv file. This csv file will contain two columns with the first column being the frequencies and the second column being the flux.

[10]:

flux_data = np.zeros((len(freqs), 2))
flux_data[:, 0] = sim_data["Flux monitor"].flux.f
flux_data[:, 1] = sim_data["Flux monitor"].flux

# ensure the directory exists
if not os.path.exists("data"):
    os.makedirs("data")

np.savetxt("data/flux_data.csv", flux_data, delimiter=",")

Integration#

Sometimes the postprocessing of the simulation data requires the calculation of certain integrals. A common example is the calculation of the overlap integral of two modes. For instance, when designing an edge coupler, we might want to calculate the overlap integral of the fiber mode and the waveguide mode. In Tidy3D, we provide a convenient function outer_dot that does this integral.

To demonstrate this, we load data from two hdf5 files. They contain the field profiles of a Gaussian mode and a waveguide mode.

[22]:

waveguide_mode_data = td.ModeSolverData.from_hdf5("misc/waveguide_mode_data.h5")
gaussian_beam_data = td.SimulationData.from_hdf5("misc/gaussian_beam_data.h5")

09:48:37 Eastern Daylight Time WARNING: updating Simulation from 2.4 to 2.5

Before calculating the overlap, let’s plot the fields to visualize them.

[23]:

waveguide_mode_data.Ey.abs.isel(f=5).plot(x="y", y="z", vmin=0, vmax=30, cmap="hot")
plt.show()

../_images/notebooks_XarrayTutorial_53_0.png

[24]:

gaussian_beam_data["monitor_0"].Ey.abs.isel(f=5).plot(x="y", y="z", vmin=0, vmax=30, cmap="hot")
plt.show()

../_images/notebooks_XarrayTutorial_54_0.png

From the field profiles, we can see that the waveguide mode and the Gaussian mode have a similar mode area so they should have a significant overlap.

Now we calculate the overlap using outer_dot and plot the result in dB.

[25]:

overlap = waveguide_mode_data.outer_dot(gaussian_beam_data["monitor_0"])
plt.plot(td.C_0 / overlap.f, 20 * np.log10(np.abs(overlap)))
plt.xlabel("Wavelength ($\mu m$)")
plt.ylabel("Overlap integral (dB)")
plt.ylim(-2, 0)
plt.show()
plt.show()

../_images/notebooks_XarrayTutorial_56_0.png

For demonstration purposes, we can also perform the integral manually using the integrate method in xarray. The overlap integral is given by

\[\frac{1}{4} \left(\mathbf{E}_1 \times \mathbf{H}_2^* - \mathbf{H}_1 \times \mathbf{E}_2^* \right) dA,\]

where \(\mathbf{E}_1\) and \(\mathbf{H}_1\) are the electric and magnetic fields of the first mode, \(\mathbf{E}_2\) and \(\mathbf{H}_2\) are the fields of the second mode, and * denotes complex conjugation.

One complication here is that the fields of the two modes might not share the same grid as they can be calculated from different simulations. The computation between different arrays does require them to share the same coordinates. In this situation, we can use the interp_like method to easily interpolate one array to the grid of the other array.

As expected, the result calculated this way is the same as that calculated from outer_dot above.

[26]:

Ey_wg = waveguide_mode_data.Ey.squeeze(drop=True)
Ez_wg = waveguide_mode_data.Ez.squeeze(drop=True)
Hy_wg = waveguide_mode_data.Hy.squeeze(drop=True)
Hz_wg = waveguide_mode_data.Hz.squeeze(drop=True)

Ey_gauss = gaussian_beam_data["monitor_0"].Ey.interp_like(Ey_wg)
Ez_gauss = gaussian_beam_data["monitor_0"].Ez.interp_like(Ey_wg)
Hy_gauss = gaussian_beam_data["monitor_0"].Hy.interp_like(Ey_wg)
Hz_gauss = gaussian_beam_data["monitor_0"].Hz.interp_like(Ey_wg)

e_self_x_h_other = Ey_wg * Hz_gauss - Ez_wg * Hy_gauss
h_self_x_e_other = Hy_wg * Ez_gauss - Hz_wg * Ey_gauss

integrand = 0.25 * (e_self_x_h_other - h_self_x_e_other)

overlap = integrand.integrate(coord=["y", "z"])

plt.plot(td.C_0 / overlap.f, 20 * np.log10(np.abs(overlap)))
plt.xlabel("Wavelength ($\mu m$)")
plt.ylabel("Overlap integral (dB)")
plt.ylim(-2, 0)
plt.show()
plt.show()

../_images/notebooks_XarrayTutorial_58_0.png

Besides the overlap integral shown above, in other example notebooks we also demonstrated how to perform integrals to calculate the cavity mode volume and the nonlinear feedback parameter.

Closing Remarks#

In this tutorial, we demonstrated various futures commonly used in Tidy3D’s data postprocessing and visualization. We used data from FieldMonitor, FluxMonitor, ModeMonitor, and ModeSolver as examples. However, the methods are general and equally applicable to data from DiffractionMonitor, FieldProjectionCartesianMonitor, and so on.

[ ]:

Advanced monitor data (xarray) manipulation and visualization

Contents

Advanced monitor data (xarray) manipulation and visualization#

Load Simulation Data#

FluxMonitor Data#

Saving Data as a CSV File#

ModeMonitor Data#

Data Selection/Slicing#

FieldMonitor Data#

Data Interpolation#

Integration#

Closing Remarks#