# Mode Solver#

Run this notebook in your browser using Binder.

This tutorial shows how to use the mode solver plugin in tidy3d.

[1]:

import numpy as np
import matplotlib.pylab as plt

import tidy3d as td
from tidy3d.constants import C_0
import tidy3d.web as web

from tidy3d.plugins import ModeSolver


## Setup#

We first set up the mode solver with information about our system. We start by setting parameters

[2]:

# size of simulation domain
Lx, Ly, Lz = 6, 6, 6
dl = .05

# waveguide information
wg_width = 1.5
wg_height = 1.0
wg_permittivity = 4.0

# central frequency
wvl_um = 2.0
freq0 = C_0 / wvl_um
fwidth = freq0 / 3

# run_time in ps
run_time = 1e-12

# automatic grid specification
grid_spec = td.GridSpec.auto(min_steps_per_wvl=20, wavelength=wvl_um)


Then we set up a simulation, in this case including a straight waveguide

[3]:

waveguide = td.Structure(
geometry = td.Box(size=(wg_width, td.inf, wg_height)),
medium = td.Medium(permittivity=wg_permittivity)
)

sim = td.Simulation(
size=(Lx, Ly, Lz),
grid_spec=grid_spec,
pml_layers=(None, None, None),
structures=[waveguide],
run_time = run_time,
)

ax = sim.plot(z=0)

[16:11:47] WARNING  No sources in simulation.                               simulation.py:435

<Figure size 432x288 with 1 Axes>


## Initialize Mode Solver#

With our system defined, we can now create our mode solver. We first need to specify on what plane we want to solve the modes using a td.Box() object.

[4]:

plane = td.Box(
center=(0,0,0),
size=(4, 0, 3.5)
)


The mode solver can now compute the modes given a ModeSpec object that specifies everything about the modes we’re looking for, for example:

• num_modes: how many modes to compute.

• target_neff: float, default=None, initial guess for the effective index of the mode; if not specified, the modes with the largest real part of the effective index are computed.

The full list of specification parameters can be found here.

[5]:

mode_spec = td.ModeSpec(
num_modes=3,
target_neff=2.0,
)


We can also specify a list of frequencies at which to solvefor the modes.

[6]:

num_freqs = 11
f0_ind = num_freqs // 2
freqs = np.linspace(freq0 - fwidth / 2, freq0 + fwidth / 2, num_freqs)


Finally, we can initialize the ModeSolver, and call the solve method.

[7]:

mode_solver = ModeSolver(
simulation=sim,
plane=plane,
mode_spec=mode_spec,
freqs=freqs,
)
mode_data = mode_solver.solve()


## Visualizing Mode Data#

The mode_info object contains information about the effective index of the mode and the field profiles, as well as the mode_spec that was used in the solver. The effective index data and the field profile data is in the form of xarray DataArrays.

We can for example plot the real part of the effective index for all three modes as follows.

[8]:

fig, ax = plt.subplots(1)
n_eff = mode_data.n_eff # real part of the effective mode index
n_eff.plot.line(x='f');

<Figure size 432x288 with 1 Axes>


The raw data can also be accessed.

[9]:

n_complex = mode_data.n_complex # complex effective index as a DataArray
n_eff = mode_data.n_eff.values  # real part of the effective index as numpy array
k_eff = mode_data.k_eff.values  # imag part of the effective index as numpy array

print(f'first mode effective index at freq0: n_eff = {n_eff[f0_ind, 0]:.2f}, k_eff = {k_eff[f0_ind, 0]:.2e}')

first mode effective index at freq0: n_eff = 1.77, k_eff = 4.05e-19


The fields stored in mode_data.fields can be visualized using in-built xarray methods.

[10]:

field_data = mode_data.fields

f, (ax1, ax2) = plt.subplots(1, 2, tight_layout=True, figsize=(10, 3))
abs(field_data.Ex.isel(mode_index=0, f=f0_ind)).plot(x='x', y='z', ax=ax1, cmap="magma")
abs(field_data.Ez.isel(mode_index=0, f=f0_ind)).plot(x='x', y='z', ax=ax2, cmap='magma')

ax1.set_title('|Ex(x, y)|')
ax1.set_aspect("equal")
ax2.set_title('|Ez(x, y)|')
ax2.set_aspect("equal")
plt.show()

<Figure size 720x216 with 4 Axes>


Alternatively, we can use the in-built plot_field method of mode_data, which also allows us to overlay the structures in the simulation. The image also looks slightly different because we have set robust=True by default, which scales the colorbar to between the 2nd and 98th percentile of the data.

[11]:

f, (ax1, ax2) = plt.subplots(1, 2, tight_layout=True, figsize=(10, 3))
mode_data.plot_field("Ex", "abs", mode_index=0, freq=freq0, ax=ax1)
mode_data.plot_field("Ez", "abs", mode_index=0, freq=freq0, ax=ax2)
plt.show()

<Figure size 720x216 with 4 Axes>


## Choosing the mode of interest#

We can also look at the other modes that were computed.

[12]:

mode_index = 1
f, (ax1, ax2) = plt.subplots(1, 2, tight_layout=True, figsize=(10, 3))
mode_data.plot_field("Ex", "abs", mode_index=mode_index, freq=freq0, ax=ax1)
mode_data.plot_field("Ez", "abs", mode_index=mode_index, freq=freq0, ax=ax2);
plt.show()

<Figure size 720x216 with 4 Axes>


This looks like an Ez-dominant mode. Finally, next-order mode has mixed polarization.

[13]:

mode_index = 2
f, (ax1, ax2) = plt.subplots(1, 2, tight_layout=True, figsize=(10, 3))
mode_data.plot_field("Ex", "abs", mode_index=mode_index, freq=freq0, ax=ax1)
mode_data.plot_field("Ez", "abs", mode_index=mode_index, freq=freq0, ax=ax2);
plt.show()

<Figure size 720x216 with 4 Axes>


## Exporting Results#

This looks promising!

Now we can choose the mode specifications to use in our mode source and mode monitors. These can be created separately, can be exported directly from the mode solver, for example:

[14]:

# Makes a modal source with geometry of plane with modes specified by mode_spec and a selected mode_index
source_time = td.GaussianPulse(freq0=freq0, fwidth=fwidth)
mode_src = mode_solver.to_source(mode_index=2, source_time=source_time, direction='-')

# Makes a mode monitor with geometry of plane.
mode_mon = mode_solver.to_monitor(name='mode', freqs=freqs)
# Offset the monitor along the propagation direction
mode_mon.center = (0, -2, 0)

[15]:

# In-plane field monitor, slightly offset along x
monitor = td.FieldMonitor(
center=(0, 0, 0.1),
size=(td.inf, td.inf, 0),
freqs=[freq0],
name='field'
)

sim = td.Simulation(
size=(Lx, Ly, Lz),
grid_spec=grid_spec,
run_time=run_time,
pml_layers=(td.PML(), td.PML(), td.PML()),
structures=[waveguide],
sources=[mode_src],
monitors=[monitor, mode_mon]
)

sim.plot(z=0);

<Figure size 432x288 with 1 Axes>

[16]:

job = web.Job(simulation=sim, task_name='mode_simulation')
sim_data = job.run(path='data/simulation_data.hdf5')

[16:11:54] INFO     Using Tidy3D credentials from stored file                      auth.py:74

[16:11:57] INFO     Uploaded task 'mode_simulation' with task_id                webapi.py:120
'e2e531f7-2096-42ec-936d-4492860db2e4'.

[16:12:01] INFO     Maximum flex unit cost: 0.20                                webapi.py:253

           INFO     status = queued                                             webapi.py:262

[16:12:12] INFO     status = preprocess                                         webapi.py:274

[16:12:26] INFO     starting up solver                                          webapi.py:278

[16:12:51] INFO     running solver                                              webapi.py:284

           INFO     status = postprocess                                        webapi.py:301

[16:13:05] INFO     status = success                                            webapi.py:307

[16:13:06] INFO     downloading file "monitor_data.hdf5" to                     webapi.py:537
"data/simulation_data.hdf5"

[16:13:07] INFO     loading SimulationData from data/simulation_data.hdf5       webapi.py:369


We can now plot the in-plane field and the modal amplitudes. Since we injected mode 2 and we just have a straight waveguide, all the power recorded by the modal monitor is in mode 2, going backwards.

[17]:

fig, ax = plt.subplots(1, 2, figsize=(10, 4))
sim_data.plot_field("field", "Ez", freq=freq0, ax=ax[0])
sim_data['mode']['amps'].sel(direction='-').abs.plot.line(x='f', ax=ax[1]);

<Figure size 720x288 with 3 Axes>


## Storing server-side computed modes#

We can also use a ModeFieldMonitor to store the modes as they are computed server-side. This is illustrated below. We will also request in the mode specification that the modes are ordered by their tm_fraction. In this particular simulation, this means the modes with largest integrated intensity of their Ez component. Similarly we could select te_fraction to sort by the Ex component in the current example.

[18]:

mode_spec.sort_by = "tm_fraction"

# Update mode source to use the highest-tm-fraction mode
mode_src.mode_spec = mode_spec
mode_src.mode_index = 0

# Update mode monitor to use the tm_fraction ordered mode_spec
mode_mon.mode_spec = mode_spec

# New monitor to record the modes computed at the mode decomposition monitor location
mode_solver_mon = td.ModeFieldMonitor(
center=mode_mon.center,
size=mode_mon.size,
freqs=mode_mon.freqs,
mode_spec=mode_spec,
name="mode_fields"
)

sim = td.Simulation(
size=(Lx, Ly, Lz),
grid_spec=grid_spec,
run_time=run_time,
pml_layers=(td.PML(), td.PML(), td.PML()),
structures=[waveguide],
sources=[mode_src],
monitors=[monitor, mode_mon, mode_solver_mon]
)

[19]:

job = web.Job(simulation=sim, task_name='mode_simulation')
sim_data = job.run(path='data/simulation_data.hdf5')

[16:13:11] INFO     Uploaded task 'mode_simulation' with task_id                webapi.py:120

[16:13:15] INFO     Maximum flex unit cost: 0.20                                webapi.py:253

           INFO     status = queued                                             webapi.py:262

[16:13:24] INFO     status = preprocess                                         webapi.py:274

[16:13:40] INFO     starting up solver                                          webapi.py:278

[16:14:05] INFO     running solver                                              webapi.py:284

           INFO     status = postprocess                                        webapi.py:301

[16:14:19] INFO     status = success                                            webapi.py:307

[16:14:20] INFO     downloading file "monitor_data.hdf5" to                     webapi.py:537
"data/simulation_data.hdf5"

[16:14:21] INFO     loading SimulationData from data/simulation_data.hdf5       webapi.py:369


Note the different ordering of the recorded modes compared to what we saw above, which used the default ordering "largest_neff".

[20]:

fig, ax = plt.subplots(1)
n_eff = sim_data["mode"].n_eff # real part of the effective mode index
n_eff.plot.line(x='f');

<Figure size 432x288 with 1 Axes>


Now the fundamental Ez-polarized mode is injected, and as before it is the only one that the mode monitor records any intensity in.

[21]:

fig, ax = plt.subplots(1, 2, figsize=(10, 4))
sim_data.plot_field("field", "Ez", freq=freq0, ax=ax[0])
sim_data['mode']['amps'].sel(direction='-').abs.plot.line(x='f', ax=ax[1]);

<Figure size 720x288 with 3 Axes>


We can also have a look at the mode fields stored in the ModeFieldMonitor either directly using xarray methods as above, or using the Tidy3D SimulationData in-built field plotting. Note that here we have to specificall

[22]:

fig, ax = plt.subplots(1, 2, figsize=(12, 4))
sim_data.plot_field("mode_fields", "Ex", freq=freq0, val="abs", mode_index=0, ax=ax[0])
sim_data.plot_field("mode_fields", "Ez", freq=freq0, val="abs", mode_index=0, ax=ax[1]);

<Figure size 864x288 with 4 Axes>


## Notes / Considerations#

• This mode solver runs locally, which means it does not require credits to run.

• It also means that the mode solver does not use subpixel-smoothening, even if this is specified in the simulation. On the other hand, when the modes are stored in a ModeFieldMonitor during a simulation run, the subpixel averaging is applied. Therefore, the latter results may not perfectly match those from a local run, and are more accurate.

• Symmetries are applied if they are defined in the simulation and the mode plane center sits on the simulation center.