Thermally tuned ring resonator#
Note: the cost of running the entire notebook is larger than 10 FlexCredits.
This example demonstrates the usage of the Tidy3D’s heat solver feature to study the influence of heat propagation on optical properties of a simple ring resonator element.
Our simulation setup will consist of two straight rectangular waveguides coupled by a ring waveguide made of Si. These elements are embedded into a SiO2 cladding on a Si wafer. The bottom of of the wafer is maintained at \(T=300\) K and the top of the cladding is exposed to air, also maintained at \(T=300\) K. The temperature field is controlled by a TiN ring heater above the ring waveguide, which is also embedded in the silicon dioxide cladding.
[1]:
import numpy as np
from matplotlib import pyplot as plt
import tidy3d as td
from tidy3d import web
This study will involve using two different solvers, thermal and optical ones. The simulation setups for these solvers are defined through two different classes HeatSimulation and Simulation, respectively. However, thanks to a general simulation description class Scene, the geometric setup of the problem can be easily transferred from one simulation to another.
Material Specification#
To create simulation setups that can be used in different solvers we need to create material specifications that contain all relevant information for each of the solvers. Specifically, when performing coupled thermal and optic simulations, each material definition will contain up to three different characteristic: 1. optic properties such as permittivity and conductivity 2. thermal properties such as thermal conductivity 3. response of optic properties to changes in temperature
Let us, for example, create a specification for silicon and silicon dioxide materials that contains all three categories of information. First, we point out that optic mediums containing response models to thermal fields (that is, information from categories 1. and 3.) can be created by using classes PerturbationMedium and PerturbationPoleResidue for non-dispersive and
dispersive material models, respectively. Second, thermal properties (that is, information from category 2.) can be specified in any material class by using field heat_spec.
Here, we will use a non-dispersive material model with a linear response of permittivity to temperature. Specification SolidSpec is used to specify thermal properties of the material.
Let us first define values of relevant material properties.
[2]:
Si_n = 3.4777 # Si refraction index
Si_n_slope = 1.86e-4 # Si thermo-optic coefficient dn/dT, 1/K
Si_k = 148e-12 # Si thermal conductivity, W / (um * K)
Si_s = 0.71 * 2.33e-12 # Si volumetric heat capacity, J / (um^3 * K)
SiO2_n = 1.444 # SiO2 refraction index
SiO2_n_slope = 1e-5 # SiO2 thermo-optic coefficient dn/dT, 1/K
SiO2_k = 1.38e-12 # SiO2 thermal conductivity, W/(um*K)
SiO2_s = 0.709 * 2.203e-12 # SiO2 volumetric heat capacity, J / (um^3 * K)
TiN_k = 28e-12 # TiN thermal conductivity W/(um*K)
TiN_s = 0.598 * 5240e-12 # TiN volumetric heat capacityJ / (um^3 * K)
TiN_sigma = 2.3e6 # Electric conductivity of TiN, S/m
Convert refraction indices and their thermal sensitivities into permittivity values.
[3]:
Si_eps = Si_n**2
Si_eps_slope = 2 * Si_n * Si_n_slope
SiO2_eps = SiO2_n**2
SiO2_eps_slope = 2 * SiO2_n * SiO2_n_slope
Pack material properties into Tidy3D medium classes.
[4]:
Si = td.PerturbationMedium(
permittivity=Si_eps,
permittivity_perturbation=td.ParameterPerturbation(
heat=td.LinearHeatPerturbation(coeff=Si_eps_slope, temperature_ref=300)
),
heat_spec=td.SolidSpec(
conductivity=Si_k,
capacity=Si_s,
),
name="Si",
)
SiO2 = td.PerturbationMedium(
permittivity=SiO2_eps,
permittivity_perturbation=td.ParameterPerturbation(
heat=td.LinearHeatPerturbation(coeff=SiO2_eps_slope, temperature_ref=300)
),
heat_spec=td.SolidSpec(
conductivity=SiO2_k,
capacity=SiO2_s,
),
name="SiO2",
)
TiN = td.PECMedium(
heat_spec=td.SolidSpec(
conductivity=TiN_k,
capacity=TiN_s,
),
name="TiN",
)
air = td.Medium(heat_spec=td.FluidSpec(), name="air")
Note that we model TiN as a perfect electric conductor and assume its optic properties are independent of temperature. Thus, a simple non-Perturbation medium class is used.
Structures#
Now we create structures that constitute our simulation setup. We start with defining basic geometric parameters.
[5]:
# in um
h_clad = 2.8 # thickness of cladding
h_box = 2 # thickness of buried oxide
h_wafer = 0.5 # wafer thickness
# waveguide parameters
w_wg = 0.5 # width
h_wg = 0.22 # thickness
couple_gap = 0.1 # between ring and straight waveguides
ring_radius = 5
# heater parameters
h_heater = 0.14 # thickness
w_heater = 2 # width
d_heater = 2 # distance between heater and waveguide
We will place the bottom surfaces of waveguides at \(z=0\) and the center of the ring resonator at \(x=0\), \(y=0\). Provided that, let us compute auxiliary variable specifying vertical positions of different elements.
[6]:
# z positions of wafer, heater, and waveguide
z_wafer = -h_box - h_wafer / 2
z_heater = h_wg + d_heater + h_heater / 2
z_wg = h_wg / 2
# y displacement of straight waveguides
y_wg = ring_radius + w_wg / 2 + couple_gap + w_wg / 2
# thickness and z position of box and cladding combined
# make structures overlap to ensure no gaps due to numerical roundoff errors
box_clad_zmax = h_clad
box_clad_zmin = -h_box - h_wafer / 2
h_box_clad = box_clad_zmax - box_clad_zmin
z_box_clad = (box_clad_zmax + box_clad_zmin) / 2
First, we create structures representing the wafer and the SiO2 layer.
[7]:
wafer = td.Structure(
geometry=td.Box(
center=[0, 0, z_wafer],
size=[td.inf, td.inf, h_wafer],
),
medium=Si,
name="wafer",
)
box_clad = td.Structure(
geometry=td.Box(
center=[0, 0, z_box_clad],
size=[td.inf, td.inf, h_box_clad],
),
medium=SiO2,
name="clad",
)
Straight waveguide are simply created using Box class.
[8]:
waveguide_top = td.Structure(
geometry=td.Box(
center=[0, y_wg, z_wg],
size=[td.inf, w_wg, h_wg],
),
medium=Si,
name="waveguide_top",
)
waveguide_bottom = td.Structure(
geometry=td.Box(
center=[0, -y_wg, z_wg],
size=[td.inf, w_wg, h_wg],
),
medium=Si,
name="waveguide_bottom",
)
Ring resonator and ring heater are created by, first, creating cylinders made of required materials and, second, creating smaller cylinders made of SiO2 that will define inner radius of these structures. The gap in the heater ring is also created using an auxiliary structure made of SiO2
[9]:
# outer ring radius
wg_ring_outer = td.Structure(
geometry=td.Cylinder(
center=[0, 0, z_wg],
axis=2,
radius=ring_radius + w_wg / 2.0,
length=h_wg,
),
medium=Si,
name="wg_ring_outer",
)
# inner ring radius
wg_ring_inner = td.Structure(
geometry=td.Cylinder(
center=[0, 0, z_wg],
axis=2,
radius=ring_radius - w_wg / 2.0,
length=w_wg * 1.5,
),
medium=SiO2,
name="wg_ring_inner",
)
# heater outer radius
heater_outer = td.Structure(
geometry=td.Cylinder(
center=[0, 0, z_heater],
axis=2,
radius=ring_radius + w_heater / 2.0,
length=h_heater,
),
medium=TiN,
name="heater_outer",
)
# heater inner radius
heater_inner = td.Structure(
geometry=td.Cylinder(
center=[0, 0, z_heater],
axis=2,
radius=ring_radius - w_heater / 2.0,
length=h_heater * 1.5,
),
medium=SiO2,
name="heater_inner",
)
#
heater_gap = td.Structure(
geometry=td.Box(
center=[0, -ring_radius, z_heater],
size=(0.5 * w_heater, 1.5 * w_heater, 1.5 * h_heater),
),
medium=SiO2,
name="heater_gap",
)
Scene#
Now the structure can be combined together into a Scene along with simulation size information.
[10]:
scene = td.Scene(
structures=[
box_clad,
wg_ring_outer,
wg_ring_inner,
waveguide_top,
waveguide_bottom,
heater_outer,
heater_inner,
heater_gap,
wafer
],
medium=air,
)
Visualization to confirm problem setup.
[11]:
fig, ax = plt.subplots(1, 3, figsize=(10, 3))
scene.plot(x=0, ax=ax[0])
scene.plot(z=z_wg, ax=ax[1])
scene.plot(z=z_heater, ax=ax[2])
plt.tight_layout()
plt.show()
We can also visualize structures according to their physical optic and thermal parameters.
[12]:
fig, ax = plt.subplots(1, 3, figsize=(20,5))
scene.plot(x=0, ax=ax[0])
scene.plot_eps(x=0, ax=ax[1])
scene.plot_heat_conductivity(x=0, ax=ax[2])
plt.tight_layout()
plt.show()
Heat Simulation#
To define a heat simulation setup (HeatSimulation) based on the created Scene object we need to define heat related specifications.
We start with specifying the simulation box dimensions. To do that we first define values for buffer space around created structures. We also define variable heat_sim_buffer that define how much larger a heat simulation is compared to an optic one. This is done to make sure that when an optic simulation is constructed we have valid temperature values in the PML regions as well. This will be applied only in x and y directions because the heat simulation is restricted to the region from the
bottom of the wafer to the top of the cladding. This is because heat equation is not solved in mediums with heat_spec=FluidSpec(), which is true for medium air defined above.
[13]:
x_buffer = 2 # buffer in x direction (from ring resonator)
y_buffer = 1 # buffer in y direction (from straight waveguides)
z_buffer = 1 # buffer in z direction
# difference between thermal and optic simulations
heat_sim_buffer = 3
Based on that compute simulation box size and its center along z.
[14]:
z_sim = (h_clad - h_box - h_wafer) / 2
h_sim = h_box + h_clad + h_wafer + 2 * z_buffer
l_heat_sim = 2 * ring_radius + 2 * x_buffer + 2 * heat_sim_buffer
w_heat_sim = 2 * ring_radius + w_heater + 2 * y_buffer + 2 * heat_sim_buffer
Before proceeding we can inspect the chosen simulation box dimensions using Scene plotting functions with custom vertical and horizontal extents.
[15]:
fig, ax = plt.subplots(1, 3, figsize=(10, 3))
scene.plot(
x=0, ax=ax[0], hlim=[-w_heat_sim / 2, w_heat_sim / 2], vlim=[z_sim - h_sim / 2, z_sim + h_sim / 2]
)
scene.plot(
z=z_wg, ax=ax[1], hlim=[-l_heat_sim / 2, l_heat_sim / 2], vlim=[-w_heat_sim / 2, w_heat_sim / 2]
)
scene.plot(
z=z_heater, ax=ax[2], hlim=[-l_heat_sim / 2, l_heat_sim / 2], vlim=[-w_heat_sim / 2, w_heat_sim / 2]
)
plt.tight_layout()
plt.show()
Having inspected that, let us proceed to other heat simulation specifications.
Create boundary conditions (HeatBoundarySpec) corresponding to temperature (TemperatureBC) boundary conditions on the bottom of the wafer, and convective (ConvectionBC) boundary condition at the cladding/air interface. Note that the convection heat transfer coefficient must be provided in units of W / (um^2 K). To specify the placement of these boundary conditions we use MediumMediumInterface specification.
[16]:
bc_air = td.HeatBoundarySpec(
placement=td.MediumMediumInterface(mediums=[air.name, SiO2.name]),
condition=td.ConvectionBC(ambient_temperature=300, transfer_coeff=10 * 1e-12),
)
bc_wafer = td.HeatBoundarySpec(
placement=td.MediumMediumInterface(mediums=[air.name, Si.name]),
condition=td.TemperatureBC(temperature=300),
)
Create a uniform heat source (UniformHeatSource) inside the ring heater to represent the Joule heating. For a given electric current in the heater \(I\) we can approximately calculate the volumetric Joule heat generation using formula \(\frac{dP}{dV} = \frac{1}{\sigma_{TiN}} \left( \frac{I}{w_{heater} h_{heater}} \right)^2\). For a current of 50 mA we get:
[17]:
heat_current = 30e-3 # A
volumetric_heat = (heat_current / w_heater / h_heater) ** 2 / TiN_sigma # W / um^3
print("Volumetric heat rate: ", volumetric_heat, "W / um^3")
heater_source = td.UniformHeatSource(structures=["heater_outer"], rate=volumetric_heat)
Volumetric heat rate: 4.991126885536821e-09 W / um^3
For spatial discretization we will use class DistanceUnstructuredGrid that specifies desired mesh size based on the distance from material interfaces. Given that h_heater is the smallest feature size of structures, we will select the mesh size near interfaces based on that. Mesh size away from interfaces is selected such that the coarse mesh would have approximately 15 points in the z-direction across the BOX + cladding region.
Selecting distance_interface as 2 * dl_interface we are enforcing approximately 2 layers of finest cells around interfaces. And distance_bulk=h_box makes the transition from the finest grid size to the coarsest grid size happen on the order of h_box. We also exclude structure "wafer" from enforcing dl_interafce since we are primarily interested in the temperature distribution in the waveguide region.
[18]:
dl_interface = h_heater / 4
dl_bulk = (h_box + h_clad) / 15
heat_grid_spec=td.DistanceUnstructuredGrid(
dl_interface=dl_interface,
dl_bulk=dl_bulk,
distance_interface=2 * dl_interface,
distance_bulk=h_box,
non_refined_structures=["wafer"],
)
We will use the results of heat simulation in a subsequent optical simulation. Thus we will record a three-dimensional temperature distribution throughout the simulation domain with a TemperatureMonitor.
[19]:
temp_mnt = td.TemperatureMonitor(size=(td.inf, td.inf, td.inf), name="all")
By default for convenience purposes, heat monitors return data on equivalent Cartesian grids. However, if the unprocessed data and/or the original unstructured grid are of interest, one can request that by setting unstructured=True. Additionally, the option conformal=True allows to request the monitor to be meshed conformally on the unstructured grid. This is useful for visualization purposes, because plots of slices of tetrahedral grids maybe be confusing for visual inspection. Here, we
will create monitors with and without conformal meshing for demonstration purposes. In total, we create 3 planar monitors: horizontal through the waveguie structure, horizontal through the heater, and a vertical one.
[20]:
temp_mnt_ugrid_wg = td.TemperatureMonitor(size=(td.inf, td.inf, 0), center=(0, 0, z_wg), name="ugrid_wg", unstructured=True, conformal=True)
temp_mnt_ugrid_heater = td.TemperatureMonitor(size=(td.inf, td.inf, 0), center=(0, 0, z_heater), name="ugrid_heater", unstructured=True)
temp_mnt_ugrid_vert = td.TemperatureMonitor(size=(0, td.inf, td.inf), center=(1e-3, 0, 0), name="ugrid_vert", unstructured=True)
Combining everything into a HeatSimulation instance. Instead of specifying structures and background medium by hand, we use a convenience class method from_scene(). Note that from the point of view of heat distribution this problem is symmetric with respect to \(x = 0\). Thus, we can use field symmetry to reduce the size, time, and cost for the heat simulation.
[21]:
# create simulation
sim = td.HeatSimulation.from_scene(
scene=scene,
size=[l_heat_sim, w_heat_sim, h_sim],
center=(0, 0, z_sim),
grid_spec=heat_grid_spec,
sources=[heater_source],
boundary_spec=[bc_air, bc_wafer],
symmetry=(1, 0, 0),
monitors=[temp_mnt, temp_mnt_ugrid_wg, temp_mnt_ugrid_heater, temp_mnt_ugrid_vert],
)
Visualization of HeatSimulation object displays: - thermal conductivity of structures (in grayscale colors), - specified boundary conditions (colored thick lines: yellow for TemperatureBC, green for HeatFluxBC, and red for ConvectionBC), - heat sources (in colored dotted hatching), - and monitors (in transparent yellow color).
[22]:
# plot the two simulations
fig, ax = plt.subplots(1, 2, figsize=(10, 3))
sim.plot_heat_conductivity(x=0.0, ax=ax[0])
sim.plot_heat_conductivity(z=z_wg, ax=ax[1], monitor_alpha=0)
plt.show()
Now we can submit our heat simulation setup for solving on our servers.
[23]:
heat_sim_data = web.run(sim, "thermal-ring")
07:58:04 CST Created task 'thermal-ring' with task_id 'he-adc79df3-0fa7-48cc-98b8-d6dd83558f8d' and task_type 'HEAT'.
Tidy3D's heat solver is currently in the beta stage. All heat simulations are charged a flat fee of 0.025 FlexCredit.
07:58:05 CST Heat solver status: queued
07:58:14 CST Heat solver status: preprocess
07:59:35 CST Heat solver status: running
07:59:40 CST Heat solver status: postprocess
08:00:23 CST Heat solver status: success
08:00:48 CST loading simulation from simulation_data.hdf5
The result of this calculation are returned as an HeatSimulationData instance containing the calculated temperature field. Both unstructured and structured monitor can be used for visualizing resulting temperature fields.
[24]:
fig, ax = plt.subplots(2, 3, figsize=(20, 10))
heat_sim_data.plot_field("all", x=0, ax=ax[0, 0])
heat_sim_data.plot_field("all", z=z_wg, ax=ax[0, 1])
heat_sim_data.plot_field("all", z=z_heater, ax=ax[0, 2])
heat_sim_data.plot_field("ugrid_vert", ax=ax[1, 0])
heat_sim_data.plot_field("ugrid_wg", ax=ax[1, 1])
heat_sim_data.plot_field("ugrid_heater", ax=ax[1, 2])
plt.tight_layout()
plt.show()
The unstructured monitors can also be used to visualize the unstructured grid used by heat solver. Note the difference between resulting grids for conformal and non-conformal monitors. In the former case the monitor plane conforms to finite element cells, while in the latter case the resulting dataset is the result of slicing unstructured tetrahedral grid.
[25]:
fig, ax = plt.subplots(1, 3, figsize=(15, 5))
sim.plot_structures_heat_conductivity(z=z_wg, ax=ax[0], alpha=0.5)
heat_sim_data["ugrid_wg"].temperature.plot(ax=ax[0], field=False)
ax[0].set_xlim([-2, 2])
ax[0].set_ylim([4, 7])
sim.plot_structures_heat_conductivity(z=z_heater, ax=ax[1], alpha=0.5)
heat_sim_data["ugrid_heater"].temperature.plot(ax=ax[1], field=False)
ax[1].set_xlim([-2, 2])
ax[1].set_ylim([3, 7])
heat_sim_data["ugrid_vert"].temperature.plot(ax=ax[2], cmap="coolwarm")
ax[2].set_xlim([3, 7])
plt.tight_layout()
plt.show()
Tidy3D provides a basic functionality for visualization of unstructured data. If a more thorough inspection is required we recommend to export the results into a vtu file using .to_vtu()method, for example:
heat_sim_data["ugrid_wg"].temperature.to_vtu("heat_mnt.vtu")
and processing it using dedicated visualization software like Paraview.
Perturbed Scene#
The obtained temperature field information can be applied to mediums present in the scene so that their optical properties are appropriately modified. This can be done by using function scene.perturbed_mediums_copy(temperature: SpatialDataArray). Note that in this case it triggers several warnings because some of the structures are extending well beyond temperature field. This is expected because the heat equation was solved on a bound region. Also, this will not have an adverse effect on
optic simulations, because in the optic simulations we are going to run the simulation domain will be fully covered y available temperature data.
[26]:
perturbed_scene = scene.perturbed_mediums_copy(
temperature=heat_sim_data["all"].temperature
)
08:00:57 CST WARNING: Provided 'temperature' does not fully cover structures[0].
08:00:59 CST WARNING: Provided 'temperature' does not fully cover structures[3].
WARNING: Provided 'temperature' does not fully cover structures[4].
WARNING: Provided 'temperature' does not fully cover structures[8].
We can compare the permittivity fields of the original and perturbed scenes. Note that both ring resonator and waveguides are affected by the heat.
[27]:
fig, ax = plt.subplots(1, 2, figsize=(10, 5))
eps_lim = [12.08, 12.18]
scene.plot_structures_eps(z=z_wg, eps_lim=eps_lim, ax=ax[0])
perturbed_scene.plot_structures_eps(z=z_wg, eps_lim=eps_lim, ax=ax[1])
plt.tight_layout()
plt.show()
Optic simulations#
Now one can run optical simulation on both scenes to investigate the influence of heat on electromagnetic response.
As mentioned above, compared to the heat simulation, the optic simulation will have smaller dimensions in \(x\) and \(y\) directions. Also, since heater and wafer structures are not expected to have significant influence on propagation in the waveguides we will set the size of simulation in z direction to 3 \(\mu\)m.
[28]:
l_optic_sim = l_heat_sim - 2 * heat_sim_buffer
w_optic_sim = w_heat_sim - 2 * heat_sim_buffer
h_optic_sim = 3
Define the wavelength range of interest and corresponding frequency values.
[29]:
# wavelength range of interest
lambda_beg = 1.5
lambda_end = 1.6
# define pulse parameters
freq_beg = td.C_0 / lambda_end
freq_end = td.C_0 / lambda_beg
freq0 = (freq_beg + freq_end) / 2
fwidth = (freq_end - freq0) / 1.5
We will inject a waveguide mode into the upper left arm of the structure.
[30]:
# calculate x-distance for mode source and mode monitors
wg_insert_x = ring_radius + 0.9 * x_buffer
mode_plane = td.Box(
center=[-wg_insert_x, y_wg, z_wg],
size=[0, 2, 2],
)
mode_source = td.ModeSource(
size=mode_plane.size,
center=mode_plane.center,
source_time=td.GaussianPulse(freq0=freq0, fwidth=fwidth),
mode_index=0,
direction="+",
num_freqs=7, # using 7 (Chebyshev) points to approximate frequency dependence
)
We will monitor the field distribution at the central wavelength in the xy plane passing through waveguide structures and mode amplitudes in the through and drop ports.
[31]:
# monitor steady state fields at central frequency over whole domain
# frequency is chosen based on simulation results for clearer demonstration
field_monitor = td.FieldMonitor(
center=[0, 0, z_wg], size=[td.inf, td.inf, 0], freqs=[1.93831331e+14], name="field"
)
# monitor the mode amps on the output waveguide
lambdas_measure = np.linspace(lambda_beg, lambda_end, 301)
freqs_measure = td.C_0 / lambdas_measure[::-1]
through_mode_mnt = td.ModeMonitor(
size=mode_plane.size,
center=[wg_insert_x, y_wg, z_wg],
freqs=freqs_measure,
mode_spec=td.ModeSpec(num_modes=2),
name="mode through",
)
drop_mode_mnt = td.ModeMonitor(
size=mode_plane.size,
center=[-wg_insert_x, -y_wg, z_wg],
freqs=freqs_measure,
mode_spec=td.ModeSpec(num_modes=2),
name="mode drop",
)
Define the remaining optic simulation parameters such as grid resolution and total simulation time, and create Simulation instances for both the original and perturbed simulation scenes.
[32]:
min_steps_per_wvl = 15
run_time = 8e-11
optic_sim = td.Simulation.from_scene(
scene=scene,
size=[l_optic_sim, w_optic_sim, h_optic_sim],
center=(0, 0, z_sim),
grid_spec=td.GridSpec.auto(
min_steps_per_wvl=min_steps_per_wvl, wavelength=td.C_0 / freq0
),
sources=[mode_source],
monitors=[field_monitor, through_mode_mnt, drop_mode_mnt],
shutoff=1e-7,
run_time=run_time,
)
perturbed_optic_sim = optic_sim.perturbed_mediums_copy(
temperature=heat_sim_data["all"].temperature
)
08:01:00 CST WARNING: 'simulation.structures[5]'is completely outside of simulation domain.
WARNING: 'simulation.structures[6]'is completely outside of simulation domain.
WARNING: 'simulation.structures[7]'is completely outside of simulation domain.
WARNING: 'simulation.structures[8]'is completely outside of simulation domain.
WARNING: Structure at structures[7] was detected as being less than half of a central wavelength from a PML on side y-min. To avoid inaccurate results or divergence, please increase gap between any structures and PML or fully extend structure through the pml.
08:01:01 CST WARNING: 'simulation.structures[5]'is completely outside of simulation domain.
WARNING: 'simulation.structures[6]'is completely outside of simulation domain.
WARNING: 'simulation.structures[7]'is completely outside of simulation domain.
WARNING: 'simulation.structures[8]'is completely outside of simulation domain.
WARNING: Structure at structures[7] was detected as being less than half of a central wavelength from a PML on side y-min. To avoid inaccurate results or divergence, please increase gap between any structures and PML or fully extend structure through the pml.
As one can see, this triggered several warnings about structures located outside of simulation domain. This is because assembled scene contains objects not relevant for optic simulations such as the heater and the substrate. By plotting simulation domain one can confirm that those structures are still present and, in fact, located in the PML layers what could be detrimental for simulation stability. Recall that for best performance absorbing PML layers must be translationally invariant along normal directions.
[33]:
fig, ax = plt.subplots(2, 2, figsize=(10, 6))
optic_sim.plot_eps(z=z_wg, ax=ax[0, 0])
optic_sim.plot_eps(x=0, ax=ax[1, 0])
perturbed_optic_sim.plot_eps(z=z_wg, ax=ax[0, 1])
perturbed_optic_sim.plot_eps(x=0, ax=ax[1, 1])
plt.tight_layout()
plt.show()
To eliminate unnecessary structures one can use method Simulation.subsection with flag remove_outside_structures=True (default: True). Additionally, one can also pass flag remove_outside_custom_mediums=True (default: False) to remove unnecessary custom medium data outside of simulation region (which includes PML layers), thus, reducing the amount of data needed to be uploaded.
[34]:
optic_sim = optic_sim.subsection(
region=optic_sim.bounding_box, remove_outside_structures=True, remove_outside_custom_mediums=True
)
perturbed_optic_sim = perturbed_optic_sim.subsection(
region=optic_sim.bounding_box, remove_outside_structures=True, remove_outside_custom_mediums=True
)
fig, ax = plt.subplots(2, 2, figsize=(10, 6))
optic_sim.plot_eps(z=z_wg, ax=ax[0, 0])
optic_sim.plot_eps(x=0, ax=ax[1, 0])
perturbed_optic_sim.plot_eps(z=z_wg, ax=ax[0, 1])
perturbed_optic_sim.plot_eps(x=0, ax=ax[1, 1])
plt.tight_layout()
plt.show()
Before proceeding to running actual FDTD simulation, we can take a look at how waveguide modes are affected in the ring resonator.
[35]:
from tidy3d.plugins.mode import ModeSolver
from tidy3d.plugins.mode.web import run as run_mode_solver
[36]:
ring_mode_plane = td.Box(size=(2, 0, 1.9), center=(-ring_radius, 0, z_wg))
mode_solver_freqs = np.linspace(freq_beg, freq_end, 11)
[37]:
mode_solver = ModeSolver(
simulation=optic_sim,
plane=ring_mode_plane,
mode_spec=td.ModeSpec(num_modes=2, bend_radius=ring_radius, bend_axis=1),
freqs=mode_solver_freqs,
direction="+",
)
mode_solver_perturbed = mode_solver.updated_copy(simulation=perturbed_optic_sim)
[38]:
mode_data = run_mode_solver(mode_solver)
mode_data_perturbed = run_mode_solver(mode_solver_perturbed, reduce_simulation=True)
08:01:13 CST WARNING: Simulation has 1.46e+06 time steps. The 'run_time' may be unnecessarily large, unless there are very long-lived resonances.
Mode solver created with task_id='fdve-e661a41a-663b-45d1-9492-9ca79a8cb83b', solver_id='mo-458b5ab3-2b75-4db8-a25f-b8e89490a55f'.
08:01:16 CST Mode solver status: queued
08:01:34 CST Mode solver status: running
08:01:44 CST Mode solver status: success
08:01:46 CST WARNING: Simulation has 1.46e+06 time steps. The 'run_time' may be unnecessarily large, unless there are very long-lived resonances.
08:01:47 CST Mode solver created with task_id='fdve-3083f771-7900-47d4-814d-991171c6b1b2', solver_id='mo-bdb1dd9c-695d-4def-9d00-58fec7fc64dd'.
08:01:50 CST Mode solver status: queued
08:02:05 CST Mode solver status: running
08:02:16 CST Mode solver status: success
In this case additional heat increased the propagation index of the waveguide mode of interest.
[39]:
mode_data.n_complex.sel(mode_index=0).real.plot()
mode_data_perturbed.n_complex.sel(mode_index=0).real.plot()
plt.legend(["no heat", "with heat"])
plt.show()
Now let us submit both optic simulations for calculation to our servers.
[40]:
batch = web.Batch(
simulations={
"ring_resonator_no_heat": optic_sim,
"ring_resonator_with_heat": perturbed_optic_sim,
}
)
batch_data = batch.run()
08:03:01 CST WARNING: Simulation has 1.55e+06 time steps. The 'run_time' may be unnecessarily large, unless there are very long-lived resonances.
08:03:02 CST Created task 'ring_resonator_no_heat' with task_id 'fdve-fdb895a0-053d-41a4-a472-bd6a4920aa3d' and task_type 'FDTD'.
View task using web UI at 'https://tidy3d.simulation.cloud/workbench?taskId=fdve-fdb895a0-053 d-41a4-a472-bd6a4920aa3d'.
WARNING: Simulation has 1.55e+06 time steps. The 'run_time' may be unnecessarily large, unless there are very long-lived resonances.
08:03:03 CST Created task 'ring_resonator_with_heat' with task_id 'fdve-01fb4ab9-11c3-4d50-a41c-d6a1b50b207f' and task_type 'FDTD'.
View task using web UI at 'https://tidy3d.simulation.cloud/workbench?taskId=fdve-01fb4ab9-11c 3-4d50-a41c-d6a1b50b207f'.
08:04:52 CST Started working on Batch.
08:05:15 CST Maximum FlexCredit cost: 15.492 for the whole batch.
Use 'Batch.real_cost()' to get the billed FlexCredit cost after the Batch has completed.
09:21:58 CST Batch complete.
Visualization of the change in steady-state field distribution at the central frequency.
[41]:
no_heat_data = batch_data["ring_resonator_no_heat"]
with_heat_data = batch_data["ring_resonator_with_heat"]
09:38:33 CST loading simulation from ./fdve-fdb895a0-053d-41a4-a472-bd6a4920aa3d.hdf5
09:39:53 CST loading simulation from ./fdve-01fb4ab9-11c3-4d50-a41c-d6a1b50b207f.hdf5
[42]:
fig, ax = plt.subplots(1, 2, figsize=(10, 4))
no_heat_data.plot_field("field", "E", "abs", z=z_wg, ax=ax[0])
with_heat_data.plot_field("field", "E", "abs", z=z_wg, ax=ax[1])
plt.tight_layout()
plt.show()
Visualization of the change in mode amplitudes in through and drop ports due to thermal effects. As expected, heating the ring resonator results in shifting of the resonance frequencies.
[43]:
fig, ax = plt.subplots(1, 2, figsize=(10, 4))
no_heat_data["mode through"].amps.sel(mode_index=0, direction="+").abs.plot(x="f", ax=ax[0])
with_heat_data["mode through"].amps.sel(mode_index=0, direction="+").abs.plot(x="f", ax=ax[0])
ax[0].legend(["no heat", "with heat"])
ax[0].set_title("Through port")
ax[0].set_ylim([0, 1])
no_heat_data["mode drop"].amps.sel(mode_index=0, direction="-").abs.plot(x="f", ax=ax[1])
with_heat_data["mode drop"].amps.sel(mode_index=0, direction="-").abs.plot(x="f", ax=ax[1])
ax[1].legend(["no heat", "with heat"])
ax[1].set_title("Drop port")
ax[1].set_ylim([0, 1])
plt.tight_layout()
plt.show()
[ ]: