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 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.
Note that package vtk is required to use full visualization capabilities of heat solver data. To install a compatible version of vtk one can use pip install tidy3d[vtk] command
[2]:
wvl_um = 1.55
freq0 = td.C_0 / wvl_um
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.
[3]:
Si_n = 3.4777 # Si refraction index
Si_n_slope = 1.86e-4 # Si thermo-optic coefficient dn/dT, 1/K
Si_k = 148e-6 # 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-6 # SiO2 thermal conductivity, W/(um*K)
SiO2_s = 0.709 * 2.203e-12 # SiO2 volumetric heat capacity, J / (um^3 * K)
TiN_k = 28e-6 # TiN thermal conductivity W/(um*K)
TiN_s = 0.598 * 5.240e-12 # TiN volumetric heat capacityJ / (um^3 * K)
TiN_sigma = 2.3 # Electric conductivity of TiN, S/um
Pack material properties into Tidy3D medium classes. Since we defined optical material properties in terms of refraction index, it is convenient to create material specifications using class method Medium.from_nk() and specify perturbation models using IndexPerturbation class.
[4]:
Si = td.PerturbationMedium.from_unperturbed(
medium=td.Medium.from_nk(n=Si_n, k=0, freq=freq0),
perturbation_spec=td.IndexPerturbation(
delta_n=td.ParameterPerturbation(
heat=td.LinearHeatPerturbation(coeff=Si_n_slope, temperature_ref=300)
),
freq=freq0,
),
heat_spec=td.SolidSpec(
conductivity=Si_k,
capacity=Si_s,
),
name="Si",
)
# an alternative but equivalent way of defining a medium specification
SiO2 = td.PerturbationMedium(
permittivity=SiO2_n ** 2,
perturbation_spec=td.IndexPerturbation(
delta_n=td.ParameterPerturbation(
heat=td.LinearHeatPerturbation(coeff=SiO2_n_slope, temperature_ref=300)
),
freq=freq0,
),
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 using ClipOperation class, that allows to combine different geometries using boolean operators.
[9]:
# ring waveguide
wg_ring = td.Structure(
geometry=td.ClipOperation(
geometry_a=td.Cylinder(
center=[0, 0, z_wg],
axis=2,
radius=ring_radius + w_wg / 2.0,
length=h_wg,
),
geometry_b=td.Cylinder(
center=[0, 0, z_wg],
axis=2,
radius=ring_radius - w_wg / 2.0,
length=w_wg * 1.5,
),
operation="difference",
),
medium=Si,
name="wg_ring",
)
# heater outer radius
heater = td.Structure(
geometry=td.ClipOperation(
geometry_a=td.Cylinder(
center=[0, 0, z_heater],
axis=2,
radius=ring_radius + w_heater / 2.0,
length=h_heater,
),
geometry_b=td.GeometryGroup(
geometries=[
td.Cylinder(
center=[0, 0, z_heater],
axis=2,
radius=ring_radius - w_heater / 2.0,
length=h_heater * 1.5,
),
td.Box(
center=[0, -ring_radius, z_heater],
size=(0.5 * w_heater, 1.5 * w_heater, 1.5 * h_heater),
),
]
),
operation="difference",
),
medium=TiN,
name="heater",
)
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,
waveguide_top,
waveguide_bottom,
heater,
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])
ax[0].set_title("Pseudo-colors")
ax[1].set_title("Relative permittivity")
ax[2].set_title("Thermal conductivity")
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 15 mA we get:
[17]:
heat_current = 15e-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"], rate=volumetric_heat)
Volumetric heat rate: 0.0012477817213842055 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 dl_interface we are enforcing approximately 1 layer 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 / 2
dl_bulk = (h_box + h_clad) / 15
heat_grid_spec=td.DistanceUnstructuredGrid(
dl_interface=dl_interface,
dl_bulk=dl_bulk,
distance_interface=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 in the form of SpatialDataArray’s. 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=(0, 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")
19:17:07 CDT Created task 'thermal-ring' with task_id 'he-c2a26381-80f7-44da-8a41-691ba3c97fe2' and task_type 'HEAT'.
Tidy3D's Heat solver is currently in the beta stage. Cost of Heat simulations is subject to change in the future.
View task using web UI at 'https://tidy3d.simulation.cloud/workbench?taskId=he-c2a26381-80f7- 44da-8a41-691ba3c97fe2'.
19:17:09 CDT status = queued
To cancel the simulation, use 'web.abort(task_id)' or '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.
19:17:11 CDT status = preprocess
19:17:35 CDT Maximum FlexCredit cost: 0.025. Use 'web.real_cost(task_id)' to get the billed FlexCredit cost after a simulation run.
starting up solver
running solver
19:17:38 CDT status = postprocess
19:18:28 CDT status = success
19:18:29 CDT View simulation result at 'https://tidy3d.simulation.cloud/workbench?taskId=he-c2a26381-80f7- 44da-8a41-691ba3c97fe2'.
19:18:31 CDT 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. For more information on it please refer to the tutorial on unstructured data classes. 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 study the simulation domain will be fully covered by available temperature data.
[26]:
perturbed_scene = scene.perturbed_mediums_copy(
temperature=heat_sim_data["all"].temperature
)
19:18:37 CDT WARNING: Provided 'temperature' does not fully cover structures[0].
WARNING: Provided 'temperature' does not fully cover structures[2].
WARNING: Provided 'temperature' does not fully cover structures[3].
WARNING: Provided 'temperature' does not fully cover structures[5].
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.2]
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.9366e14], 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
)
19:18:38 CDT WARNING: 'simulation.structures[4]' is outside of the simulation domain.
WARNING: 'simulation.structures[5]' is outside of the simulation domain.
WARNING: 'simulation.structures[4]' is outside of the simulation domain.
WARNING: 'simulation.structures[5]' is outside of the simulation domain.
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)
Note that the simulation object passed to the mode solver for the perturbed case contain custom medium data. Sometimes this can lead to large amounts of uploaded/downloaded data. To avoid that, the mode solver web API by default restricts passed simulation object containing custom mediums to mode solver plane only. When that happens, a warning is displayed. Setting reduce_simulation=True forces such simulation reductions for all cases (even if no custom mediums present) and effectively
silences the warning.
[38]:
# mode_data = mode_solver.solve()
# mode_data_perturbed = mode_solver_perturbed.solve()
mode_data = run_mode_solver(mode_solver)
mode_data_perturbed = run_mode_solver(mode_solver_perturbed)
19:18:42 CDT WARNING: Simulation has 1.54e+06 time steps. The 'run_time' may be unnecessarily large, unless there are very long-lived resonances.
19:18:43 CDT Mode solver created with task_id='fdve-f9b7f9b1-90af-4fb7-bf74-2c2b75480817', solver_id='mo-da4cf69f-ce9a-48e8-8dc0-e29d675f6af9'.
19:18:45 CDT Mode solver status: queued
19:19:03 CDT Mode solver status: running
19:19:09 CDT Mode solver status: success
19:19:11 CDT WARNING: The associated 'Simulation' object contains custom mediums. It will be automatically restricted to the mode solver plane to reduce data for uploading. To force uploading the original 'Simulation' object use 'reduce_simulation=False'. Setting 'reduce_simulation=True' will force simulation reduction in all cases and silence this warning.
19:19:12 CDT 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-10cdc768-529a-4ddd-a7d4-60ebc781d5b8', solver_id='mo-1496ad99-f4c2-4f9b-809d-762e67371ef0'.
19:19:14 CDT Mode solver status: queued
19:19:16 CDT Mode solver status: running
19:19:22 CDT 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()
19:19:23 CDT WARNING: Simulation has 1.54e+06 time steps. The 'run_time' may be unnecessarily large, unless there are very long-lived resonances.
19:19:24 CDT WARNING: Simulation has 1.55e+06 time steps. The 'run_time' may be unnecessarily large, unless there are very long-lived resonances.
19:19:33 CDT Started working on Batch containing 2 tasks.
19:19:41 CDT Maximum FlexCredit cost: 15.401 for the whole batch.
Use 'Batch.real_cost()' to get the billed FlexCredit cost after the Batch has completed.
19:54:26 CDT 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"]
[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()
[ ]: