Aerodynamics & Electromagnetic Simulation | Flexcompute
  • Sign up for Free
  • Sign up for Free
  • Sign in
  • Talk to an Expert
Aerodynamics & Electromagnetic Simulation | Flexcompute
  • Flow360 (CFD)
    • Flow360 Solver
    • Flow360 Documentation
    • Publications
    • Technical Support
    • Contact Us
  • Tidy3D (EM)
    • Tidy3D Solver
    • Tidy3D Documentation
    • Tidy3D Knowledge Base
    • Learning FDTD through Examples
    • White Paper
    • Publications
    • Technical Support
    • Contact Us
  • Solution
    • Flow360 Solution
  • Resources
    • FDTD 101
    • Aircraft CFD 101
    • CFD Essentials
    • Newsroom
    • Blog
    • Python FDTD
    • FDTD
    • FDTD Acceleration
  • About
    • Vision
    • Story
    • Team
    • Careers
    • Contact Us
  • Login
    • Flow360
    • Tidy3D
Learning FDTD through Examples  / Uniform grating coupler...
Sign up now to simulate .ipynb

Uniform grating coupler¶

[1]:
# basic imports
import numpy as np
import matplotlib.pylab as plt

# tidy3d imports
import tidy3d as td
import tidy3d.web as web

Problem Setup¶

In this example, we model a 3D grating coupler in a Silicon on Insulator (SOI) platform.

A basic schematic of the design is shown below. The simulation is about 19um x 4um x 5um with a wavelength of 1.55um and takes about 1 minute to simulate 10,000 time steps.

In the simulation, we inject a modal source into the waveguide and propagate it towards the grating structure. The radiation from the grating coupler is then measured with a near field monitor and we use a far field projection to inspect the angular dependence of the radiation.

[2]:
# basic parameters (note, all length units are microns)
nm = 1e-3
wavelength = 1550 * nm

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

# waveguide
wg_width = 400 * nm
wg_height = 220 * nm
wg_length = 2 * wavelength

# surrounding
sub_height = 2.0
air_height = 2.0
buffer = 0.5 * wavelength

# coupler
cp_width = 4 * wavelength
cp_length = 8 * wavelength
taper_length = 6 * wavelength
[3]:
# sizes
Lx = buffer + wg_length + taper_length + cp_length
Ly = buffer + cp_width + buffer
Lz = sub_height + wg_height + air_height
sim_size = [Lx, Ly, Lz]

# convenience variables to store center of coupler and waveguide
wg_center_x = +Lx / 2 - buffer - (wg_length + taper_length) / 2
cp_center_x = -Lx / 2 + buffer + cp_length / 2
wg_center_z = -Lz / 2 + sub_height + wg_height / 2
cp_center_z = -Lz / 2 + sub_height + wg_height / 2

# materials
Air = td.Medium(permittivity=1.0)
Si = td.Medium(permittivity=3.47**2)
SiO2 = td.Medium(permittivity=1.44**2)

# source parameters
freq0 = td.C_0 / wavelength
fwidth = freq0 / 10
run_time = 100 / fwidth

Mode Solve¶

To determine the pitch of the waveguide for a given design angle, we need to compute the effective index of the waveguide mode being coupled into. For this, we set up a simple simulation of the coupler region and use the mode solver to get the effective index. We will not run this simulation, we just add a ModeMonitor object in order to construct and run the mode solver below, and get the effective index of the wide-waveguide region.

Note: because the simulation is just being used for the mode solver, we can safely ignore the warning about lack of sources.

[4]:
# grating parameters
design_theta_deg = -30
design_theta_rad = np.pi * design_theta_deg / 180
grating_height = 70 * nm

# do a mode solve to get neff of the coupler

sub = td.Structure(
    geometry=td.Box(center=[0, 0, -Lz / 2], size=[td.inf, td.inf, 2 * sub_height]),
    medium=SiO2,
    name="substrate",
)

cp = td.Structure(
    geometry=td.Box(
        center=[0, 0, cp_center_z - grating_height / 4],
        size=[td.inf, cp_width, wg_height - grating_height / 2],
    ),
    medium=Si,
    name="coupler",
)

mode_plane = td.Box(center=(0, 0, 0), size=(0, 8 * cp_width, 8 * wg_height))


sim = td.Simulation(
    size=sim_size,
    grid_spec=grid_spec,
    structures=[sub, cp],
    sources=[],
    monitors=[],
    boundary_spec=td.BoundarySpec.all_sides(boundary=td.PML()),
    run_time=1e-12,
)

f, (ax1, ax2, ax3) = plt.subplots(1, 3, tight_layout=True, figsize=(15, 3))
sim.plot(x=0, ax=ax1)
sim.plot(y=0, ax=ax2)
sim.plot(z=0, ax=ax3)
plt.show()
[07:52:00] WARNING: No sources in simulation.                                                     simulation.py:584

Compute Effective index for Grating Pitch Design¶

[5]:
from tidy3d.plugins.mode import ModeSolver

ms = ModeSolver(
    simulation=sim, plane=mode_plane, mode_spec=td.ModeSpec(), freqs=[freq0]
)
mode_output = ms.solve()
           WARNING: Use the remote mode solver with subpixel averaging for better accuracy       mode_solver.py:118
           through 'tidy3d.plugins.mode.web.run(...)'.                                                             
[6]:
f, (ax1, ax2, ax3) = plt.subplots(1, 3, tight_layout=True, figsize=(15, 3))
ms.plot_field("Ex", val="real", robust=False, ax=ax1)
ms.plot_field("Ey", val="real", robust=False, ax=ax2)
ms.plot_field("Ez", val="real", robust=False, ax=ax3)
plt.show()
[07:52:01] WARNING: No sources in simulation.                                                     simulation.py:584
[7]:
neff = float(mode_output.n_eff)
print(f"effective index = {neff:.4f}")
effective index = 2.6885

Create Simulation¶

Now we set up the grating coupler to simulate in Tidy3D.

[8]:
# gratings
pitch = wavelength / (neff - np.sin(abs(design_theta_rad)))
grating_length = pitch / 2.0
num_gratings = int(cp_length / pitch)

sub = td.Structure(
    geometry=td.Box(
        center=[0, 0, -Lz / 2],
        size=[td.inf, td.inf, 2 * sub_height],
    ),
    medium=SiO2,
    name="substrate",
)

wg = td.Structure(
    geometry=td.Box(
        center=[wg_center_x, 0, wg_center_z],
        size=[buffer + wg_length + taper_length + cp_length / 2, wg_width, wg_height],
    ),
    medium=Si,
    name="waveguide",
)

cp = td.Structure(
    geometry=td.Box(
        center=[cp_center_x, 0, cp_center_z],
        size=[cp_length, cp_width, wg_height],
    ),
    medium=Si,
    name="coupler",
)

tp = td.Structure(
    geometry=td.PolySlab(
        vertices=[
            [cp_center_x + cp_length / 2 + taper_length, +wg_width / 2],
            [cp_center_x + cp_length / 2 + taper_length, -wg_width / 2],
            [cp_center_x + cp_length / 2, -cp_width / 2],
            [cp_center_x + cp_length / 2, +cp_width / 2],
        ],
        slab_bounds=(wg_center_z - wg_height / 2, wg_center_z + wg_height / 2),
        axis=2,
    ),
    medium=Si,
    name="taper",
)

grating_left_x = cp_center_x - cp_length / 2
gratings = [
    td.Structure(
        geometry=td.Box(
            center=[
                grating_left_x + (i + 0.5) * pitch,
                0,
                cp_center_z + wg_height / 2 - grating_height / 2,
            ],
            size=[grating_length, cp_width, grating_height],
        ),
        medium=Air,
        name=f"{i}th_grating",
    )
    for i in range(num_gratings)
]
[9]:
# distance to near field monitor
nf_offset = 50 * nm

plane_monitor = td.FieldMonitor(
    center=[0, 0, cp_center_z],
    size=[td.inf, td.inf, 0],
    freqs=[freq0],
    name="full_domain_fields",
)

rad_monitor = td.FieldMonitor(
    center=[0, 0, 0], size=[td.inf, 0, td.inf], freqs=[freq0], name="radiated_fields"
)

near_field_monitor = td.FieldMonitor(
    center=[cp_center_x, 0, cp_center_z + wg_height / 2 + nf_offset],
    size=[td.inf, td.inf, 0],
    freqs=[freq0],
    name="radiated_near_fields",
)
[10]:
sim = td.Simulation(
    size=sim_size,
    grid_spec=grid_spec,
    structures=[sub, wg, cp, tp] + gratings,
    sources=[],
    monitors=[plane_monitor, rad_monitor, near_field_monitor],
    run_time=run_time,
    boundary_spec=td.BoundarySpec.all_sides(boundary=td.PML()),
)
[07:52:03] WARNING: No sources in simulation.                                                     simulation.py:584

Make Modal Source¶

[11]:
source_plane = td.Box(
    center=[Lx / 2 - buffer, 0, cp_center_z],
    size=[0, 8 * wg_width, 8 * wg_height],
)

ms = ModeSolver(
    simulation=sim, plane=source_plane, mode_spec=td.ModeSpec(), freqs=[freq0]
)
mode_output = ms.solve()

f, (ax1, ax2, ax3) = plt.subplots(1, 3, tight_layout=True, figsize=(15, 3))
ms.plot_field("Ex", val="real", robust=False, ax=ax1)
ms.plot_field("Ey", val="real", robust=False, ax=ax2)
ms.plot_field("Ez", val="real", robust=False, ax=ax3)
plt.show()
           WARNING: Use the remote mode solver with subpixel averaging for better accuracy       mode_solver.py:118
           through 'tidy3d.plugins.mode.web.run(...)'.                                                             
           WARNING: No sources in simulation.                                                     simulation.py:584
[12]:
source_time = td.GaussianPulse(freq0=freq0, fwidth=fwidth)
mode_src = ms.to_source(mode_index=0, direction="-", source_time=source_time)
sim = sim.copy(update={"sources": [mode_src]})
[13]:
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, tight_layout=True, figsize=(14, 3))
sim.plot(x=0, ax=ax1)
sim.plot(y=0.01, ax=ax2)
sim.plot(z=0.1, ax=ax3)
plt.show()
[14]:
mode_src.help()
╭───────────────────────────────── <class 'tidy3d.components.source.ModeSource'> ─────────────────────────────────╮
│ Injects current source to excite modal profile on finite extent plane.                                          │
│                                                                                                                 │
│ ╭─────────────────────────────────────────────────────────────────────────────────────────────────────────────╮ │
│ │ ModeSource(                                                                                                 │ │
│ │ │   type='ModeSource',                                                                                      │ │
│ │ │   center=(12.012500000000001, 0.0, -3.191891195797325e-16),                                               │ │
│ │ │   size=(0.0, 3.2, 1.76),                                                                                  │ │
│ │ │   source_time=GaussianPulse(                                                                              │ │
│ │ │   │   amplitude=1.0,                                                                                      │ │
│ │ │   │   phase=0.0,                                                                                          │ │
│ │ │   │   type='GaussianPulse',                                                                               │ │
│ │ │   │   freq0=193414489032258.06,                                                                           │ │
│ │ │   │   fwidth=19341448903225.805,                                                                          │ │
│ │ │   │   offset=5.0                                                                                          │ │
│ │ │   ),                                                                                                      │ │
│ │ │   name=None,                                                                                              │ │
│ │ │   num_freqs=1,                                                                                            │ │
│ │ │   direction='-',                                                                                          │ │
│ │ │   mode_spec=ModeSpec(                                                                                     │ │
│ │ │   │   num_modes=1,                                                                                        │ │
│ │ │   │   target_neff=None,                                                                                   │ │
│ │ │   │   num_pml=(0, 0),                                                                                     │ │
│ │ │   │   filter_pol=None,                                                                                    │ │
│ │ │   │   angle_theta=0.0,                                                                                    │ │
│ │ │   │   angle_phi=0.0,                                                                                      │ │
│ │ │   │   precision='single',                                                                                 │ │
│ │ │   │   bend_radius=None,                                                                                   │ │
│ │ │   │   bend_axis=None,                                                                                     │ │
│ │ │   │   track_freq='central',                                                                               │ │
│ │ │   │   group_index_step=False,                                                                             │ │
│ │ │   │   type='ModeSpec'                                                                                     │ │
│ │ │   ),                                                                                                      │ │
│ │ │   mode_index=0                                                                                            │ │
│ │ )                                                                                                           │ │
│ ╰─────────────────────────────────────────────────────────────────────────────────────────────────────────────╯ │
│                                                                                                                 │
│      angle_phi = 0.0                                                                                            │
│    angle_theta = 0.0                                                                                            │
│   bounding_box = Box(                                                                                           │
│                      type='Box',                                                                                │
│                      center=(                                                                                   │
│                          12.012500000000001,                                                                    │
│                          0.0,                                                                                   │
│                          -3.3306690738754696e-16                                                                │
│                      ),                                                                                         │
│                      size=(0.0, 3.2, 1.76)                                                                      │
│                  )                                                                                              │
│         bounds = (                                                                                              │
│                      (12.012500000000001, -1.6, -0.8800000000000003),                                           │
│                      (12.012500000000001, 1.6, 0.8799999999999997)                                              │
│                  )                                                                                              │
│         center = (12.012500000000001, 0.0, -3.191891195797325e-16)                                              │
│      direction = '-'                                                                                            │
│ frequency_grid = array([1.93414489e+14])                                                                        │
│       geometry = Box(                                                                                           │
│                      type='Box',                                                                                │
│                      center=(12.012500000000001, 0.0, -3.191891195797325e-16),                                  │
│                      size=(0.0, 3.2, 1.76)                                                                      │
│                  )                                                                                              │
│ injection_axis = 0                                                                                              │
│     mode_index = 0                                                                                              │
│      mode_spec = ModeSpec(                                                                                      │
│                      num_modes=1,                                                                               │
│                      target_neff=None,                                                                          │
│                      num_pml=(0, 0),                                                                            │
│                      filter_pol=None,                                                                           │
│                      angle_theta=0.0,                                                                           │
│                      angle_phi=0.0,                                                                             │
│                      precision='single',                                                                        │
│                      bend_radius=None,                                                                          │
│                      bend_axis=None,                                                                            │
│                      track_freq='central',                                                                      │
│                      group_index_step=False,                                                                    │
│                      type='ModeSpec'                                                                            │
│                  )                                                                                              │
│           name = None                                                                                           │
│      num_freqs = 1                                                                                              │
│    plot_params = PlotParams(                                                                                    │
│                      alpha=0.4,                                                                                 │
│                      edgecolor='limegreen',                                                                     │
│                      facecolor='limegreen',                                                                     │
│                      fill=True,                                                                                 │
│                      hatch=None,                                                                                │
│                      linewidth=3.0,                                                                             │
│                      type='PlotParams'                                                                          │
│                  )                                                                                              │
│           size = (0.0, 3.2, 1.76)                                                                               │
│    source_time = GaussianPulse(                                                                                 │
│                      amplitude=1.0,                                                                             │
│                      phase=0.0,                                                                                 │
│                      type='GaussianPulse',                                                                      │
│                      freq0=193414489032258.06,                                                                  │
│                      fwidth=19341448903225.805,                                                                 │
│                      offset=5.0                                                                                 │
│                  )                                                                                              │
│           type = 'ModeSource'                                                                                   │
│      zero_dims = [0]                                                                                            │
╰─────────────────────────────────────────────────────────────────────────────────────────────────────────────────╯

Run Simulation¶

Run the simulation and plot the field patterns

[15]:
# create a project, upload to our server to run
job = web.Job(simulation=sim, task_name="grating_coupler", verbose=True)
sim_data = job.run(path="data/grating_coupler.hdf5")
print(sim_data.log)
[07:52:06] Created task 'grating_coupler' with task_id 'fdve-daa5cf67-39bd-4443-adb1-f1e27b33411bv1'. webapi.py:148
           View task using web UI at                                                                  webapi.py:150
           'https://tidy3d.simulation.cloud/workbench?taskId=fdve-daa5cf67-39bd-4443-adb1-f1e27b33411              
           bv1'.                                                                                                   
Output()


[07:52:11] status = queued                                                                            webapi.py:280
Output()
[07:52:14] status = preprocess                                                                        webapi.py:274

[07:52:22] Maximum FlexCredit cost: 1.623. Use 'web.real_cost(task_id)' to get the billed FlexCredit  webapi.py:297
           cost after a simulation run.                                                                            
           starting up solver                                                                         webapi.py:301
           running solver                                                                             webapi.py:311
Output()
[07:54:45] early shutoff detected, exiting.                                                           webapi.py:325


           status = postprocess                                                                       webapi.py:342
Output()
[07:55:06] status = success                                                                           webapi.py:349

Output()


[07:55:14] loading SimulationData from data/grating_coupler.hdf5                                      webapi.py:521
Simulation domain Nx, Ny, Nz: [1157, 334, 102]
Applied symmetries: (0, 0, 0)
Number of computational grid points: 4.0500e+07.
Using subpixel averaging: True
Number of time steps: 1.3521e+05
Automatic shutoff factor: 1.00e-05
Time step (s): 3.8240e-17


Compute source modes time (s):     0.4786
Compute monitor modes time (s):    0.0030
Rest of setup time (s):            6.3535

Running solver for 135206 time steps...
- Time step   1076 / time 4.11e-14s (  0 % done), field decay: 1.00e+00
- Time step   5408 / time 2.07e-13s (  4 % done), field decay: 1.00e+00
- Time step  10816 / time 4.14e-13s (  8 % done), field decay: 1.75e-02
- Time step  16224 / time 6.20e-13s ( 12 % done), field decay: 3.27e-03
- Time step  21632 / time 8.27e-13s ( 16 % done), field decay: 2.48e-04
- Time step  27041 / time 1.03e-12s ( 20 % done), field decay: 6.92e-05
- Time step  32449 / time 1.24e-12s ( 24 % done), field decay: 2.11e-05
- Time step  37857 / time 1.45e-12s ( 28 % done), field decay: 6.44e-06
Field decay smaller than shutoff factor, exiting solver.

Solver time (s):                   133.6126
Data write time (s):               0.1986

[16]:
fig, (ax1, ax2) = plt.subplots(2, 1, tight_layout=True, figsize=(14, 8))
sim_data.plot_field("full_domain_fields", "Ey", f=freq0, ax=ax1)
sim_data.plot_field("radiated_fields", "Ey", f=freq0, ax=ax2)
plt.show()

Far Field Projection¶

Now we use the Tidy3D's FieldProjector to compute the anglular dependence of the far field scattering based on the near field monitor.

[17]:
# create range of angles to probe (note: polar coordinates, theta = 0 corresponds to vertical (z axis))
num_angles = 1101
thetas = np.linspace(-np.pi / 2, np.pi / 2, num_angles)

# make a near-to-far monitor specifying the observation angles and frequencies of interest
monitor_n2f = td.FieldProjectionAngleMonitor(
    center=near_field_monitor.center,
    size=near_field_monitor.size,
    normal_dir="+",
    freqs=[freq0],
    theta=thetas,
    phi=[0.0],
    name="n2f",
)

# make a near field to far field projector with the near field monitor data
near_field_surface = td.FieldProjectionSurface(
    monitor=near_field_monitor, normal_dir="+"
)
n2f = td.FieldProjector(sim_data=sim_data, surfaces=[near_field_surface])

# compute the far_fields
far_fields = n2f.project_fields(monitor_n2f)

# Compute the scattered cross section
Ps = np.abs(far_fields.radar_cross_section.sel(f=freq0).values[0, ...])
Output()


[18]:
# plot the angle dependence
fig, ax = plt.subplots(subplot_kw={"projection": "polar"}, figsize=(5, 5))
ax.plot(thetas, Ps, label="measured")
ax.plot(
    [design_theta_rad, design_theta_rad],
    [0, np.max(Ps) * 0.7],
    "r--",
    alpha=0.8,
    label="design angle",
)
ax.set_theta_zero_location("N")
ax.set_theta_direction(-1)
ax.set_yticklabels([])
ax.set_title("Scattered Cross-section (arb. units)", va="bottom")
plt.legend()
plt.show()