Defining spatially-varying sources#
In this example we illustrate how to incorporate a source with a custom spatial dependence in a Tidy3D simulation. We will illustrate this both using data taken from a monitor from another Tidy3D simulation and using a simple array defined from scratch. This notebook will cover both custom “field” sources and custom “current” sources.
The custom field source can be used to inject a specific (E, H) field distribution on a plane, e.g. coming from another simulation. Internally, we use the equivalence principle to compute the actual source currents (all sources in FDTD have to be converted to current sources). Because of this, the custom field source will only produce reliable results if the provided fields decay by the edges of the source plane, or if they extend through the simulation boundaries, and are well-matched to those boundaries.
The custom current source is simpler, and can be used to simply provide a specific electric current and magnetic current distribution within a volume. The component labels Ex, Ey, Ez, Hx, Hy, Hz in the dataset are directly assigned as electric/magnetic current amplitudes Jx, Jy, Jz, Mx, My, Mz.
If you are new to the finite-difference time-domain (FDTD) method, we highly recommend going through our FDTD101 tutorials.
[1]:
# standard python imports
import numpy as np
import matplotlib.pyplot as plt
# tidy3D import
import tidy3d as td
from tidy3d import web
Starting simulation with an in-built source#
We will first run a simulation with a GaussianBeam source propagating in empty space. The definition of the source below is for a converging Gaussian beam (negative waist_distance parameter).
[2]:
# Free-space wavelength (in um) and frequency (in Hz)
lambda0 = 1.55
freq0 = td.C_0 / lambda0
fwidth = freq0 / 10
# Simulation size and run time
sim_size = [10, 10, 10]
run_time = 20 / fwidth
# Grid specification
grid_spec = td.GridSpec.auto(wavelength=lambda0)
# In-built GaussianBeam source
pulse = td.GaussianPulse(freq0=freq0, fwidth=fwidth)
waist_radius = 2
src_pos = -4
gaussian_source = td.GaussianBeam(
source_time=pulse,
center=(src_pos, 0, 0),
size=(0, 10, 10),
waist_radius=waist_radius,
waist_distance=-8,
direction="+",
)
We will use monitors to record the fields at the beam waist along the propagation direction, and inject those in another simulation using a CustomFieldSource. The best way to achieve this is to take into account the numerical grid used in the FDTD algorithm. Specifically, we need to record the raw fields at their Yee grid locations in order to achieve the highest accuracy in injecting the
recorded data. For this purpose, we need to set colocate=False in the monitor definition. Furthermore, it is best to set a slightly nonzero size in the propagation direction. A monitor with a zero size will interpolate the fields to the exact monitor location along the zero-size dimension, which loses a bit of information. On the other hand, a monitor with a small but nonzero size along the propagation direction will store the fields exactly as they are in the simulation. We will illustrate
the difference in the results below.
[3]:
# Monitor for propagation
mnt_xy = td.FieldMonitor(
center=(0, 0, 0), size=(10, 10, 0), freqs=[freq0], name="field_xy"
)
# Monitors for forward and backward flux
mnt_flux_pos = src_pos - 0.5
mnt_flux_f = td.FluxMonitor(
center=(-mnt_flux_pos, 0, 0), size=(0, td.inf, td.inf), freqs=[freq0], name="flux_f"
)
mnt_flux_b = td.FluxMonitor(
center=(mnt_flux_pos, 0, 0), size=(0, td.inf, td.inf), freqs=[freq0], name="flux_b"
)
# Monitor to be used as custom source, 0D along x
mnt_yz_1 = td.FieldMonitor(
center=(-src_pos, 0, 0), size=(0, 8, 8), freqs=[freq0], colocate=False, name="yz_zero_size_x"
)
# Monitor to be used as custom source, small nonzero along x
mnt_yz_2 = td.FieldMonitor(
center=(-src_pos, 0, 0), size=(1e-5, 8, 8), freqs=[freq0], colocate=False, name="yz_nonzero_size_x"
)
[10:29:40] WARNING: Default value for the field monitor monitor.py:261 'colocate' setting has changed to 'True' in Tidy3D 2.4.0. All field components will be colocated to the grid boundaries. Set to 'False' to get the raw fields on the Yee grid instead.
[4]:
sim = td.Simulation(
size=sim_size,
grid_spec=grid_spec,
sources=[gaussian_source],
monitors=[mnt_xy, mnt_flux_f, mnt_flux_b, mnt_yz_1, mnt_yz_2],
run_time=run_time,
boundary_spec=td.BoundarySpec.all_sides(boundary=td.PML()),
)
sim.plot(z=0)
plt.show()
[5]:
sim_data = web.run(sim, task_name="free space gaussian", verbose=True)
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- c96bb6c3-fc96-4bcd-91f2-47134caabfbdv1'.
[6]:
fig, ax = plt.subplots(1)
sim_data.plot_field("field_xy", field_name="Ey", val="abs", ax=ax)
ax.set_xlim([-sim.size[0] / 2, sim.size[0] / 2])
ax.set_ylim([-sim.size[1] / 2, sim.size[1] / 2])
ax.set_title(
f"Flux fwd={float(sim_data['flux_f'].flux):1.2e}, bck={float(sim_data['flux_b'].flux):1.2e}"
)
plt.show()
Custom source from simulation data#
Now we can use the recorded data as an input to a CustomFieldSource in another simulation. We will run two simulations to illustrate the difference between zero and nonzero size of the monitor along x. Note that we can create the source directly using the FieldData method to_source.
Note that we can center the source anywhere we want. The size of the source is by default taken from the FieldData used to create it. However, the source size needs to be 0 along one of the directions, so in the case of the nonzero_size_x data, we have to manually reset it.
[7]:
custom_src_1 = sim_data["yz_zero_size_x"].to_source(source_time=pulse, center=(0, 1, 0))
custom_src_2 = sim_data["yz_nonzero_size_x"].to_source(
source_time=pulse,
center=(0, 1, 0),
# size is by default taken from the data, but it must be reset to zero along x
size=(0, 8, 8),
)
sim_1 = td.Simulation(
size=sim_size,
grid_spec=grid_spec,
sources=[custom_src_1],
monitors=[mnt_xy, mnt_flux_f, mnt_flux_b],
run_time=run_time,
boundary_spec=td.BoundarySpec.all_sides(boundary=td.PML()),
)
sim_2 = td.Simulation(
size=sim_size,
grid_spec=grid_spec,
sources=[custom_src_2],
monitors=[mnt_xy, mnt_flux_f, mnt_flux_b],
run_time=run_time,
boundary_spec=td.BoundarySpec.all_sides(boundary=td.PML()),
)
[8]:
batch = web.Batch(simulations={"zero_x": sim_1, "nonzero_x": sim_2}, verbose=True)
batch_results = batch.run(path_dir="data")
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- 9e2c69b7-594e-4965-8f3c-0e3a3b999043v1'.
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- 1670879d-36a6-43b9-bc6e-73122dc47e9av1'.
[10:32:19] Started working on Batch. container.py:475
[10:32:29] Maximum FlexCredit cost: 0.050 for the whole batch. container.py:479 Use 'Batch.real_cost()' to get the billed FlexCredit cost after the Batch has completed.
[10:32:53] Batch complete. container.py:522
Below we plot the injected field in the two simulations. As can be seen, there is a bit of backwards ppower injected in the simulation in which we used a FieldMonitor with size 0 along x, while the injection is extremely clean in the case where the size of the monitor was slightly larger than zero. As mentioned above, this is because the fields from the original simulation are captured exactly as they are on the numerical grid.
[9]:
fig, ax = plt.subplots(1, 2, figsize=(10, 4))
sim_data = batch_results["zero_x"]
sim_data.plot_field("field_xy", field_name="Ey", val="abs", ax=ax[0])
ax[0].set_xlim([-sim.size[0] / 2, sim.size[0] / 2])
ax[0].set_ylim([-sim.size[1] / 2, sim.size[1] / 2])
title = f"Flux: fwd={float(sim_data['flux_f'].flux):1.2e}, bck={float(sim_data['flux_b'].flux):1.2e}"
title += "\nzero-size FieldData along x"
ax[0].set_title(title)
sim_data = batch_results["nonzero_x"]
sim_data.plot_field("field_xy", field_name="Ey", val="abs", ax=ax[1])
ax[1].set_xlim([-sim.size[0] / 2, sim.size[0] / 2])
ax[1].set_ylim([-sim.size[1] / 2, sim.size[1] / 2])
title = f"Flux: fwd={float(sim_data['flux_f'].flux):1.2e}, bck={float(sim_data['flux_b'].flux):1.2e}"
title += "\nnonzero-size FieldData along x"
ax[1].set_title(title)
plt.show()
Importing arbitrary fields#
Next we show how to import arbitrary fields. We create numpy arrays to emulate the same Gaussian beam input. Data can be imported from file in exactly the same way, once it is read into a numpy array. We create two types of datasets: the first one contains the Ey field only, while the second one contains Ey and Hz. We can provide any combination of field components to the custom source, as long as at least one of the tangential fields is defined. However, as we will see, providing
both E and H is required to make the source directional.
Note: the fields are assumed to be defined w.r.t. the local coordinates of the source, that is to say with the coordinate origin lying at the source center. This makes it easy to specify a source profile and use that anywhere in a simulation. We illustrate this below by defining the dataset to be centered around (0, 0, 0), while the source is placed at a different location in the simulation.
[10]:
# Scalar gaussian field with the same waist_radius as the beam source above
ys, zs = np.linspace(-4, 4, 101), np.linspace(-4, 4, 101)
y_grid, z_grid = np.meshgrid(ys, zs)
scalar_gaussian = np.exp(-(y_grid**2 + z_grid**2) / waist_radius**2)
# Field dataset defining only the E-field
dataset_E = td.FieldDataset(
Ey=td.ScalarFieldDataArray(
scalar_gaussian[None, ..., None],
coords={
"x": [0],
"y": ys,
"z": zs,
"f": [freq0],
},
)
)
# Field dataset defining both E and H
dataset_EH = td.FieldDataset(
Ey=td.ScalarFieldDataArray(
scalar_gaussian[None, ..., None],
coords={
"x": [0],
"y": ys,
"z": zs,
"f": [freq0],
},
),
Hz=td.ScalarFieldDataArray(
scalar_gaussian[None, ..., None] / td.ETA_0,
coords={
"x": [0],
"y": ys,
"z": zs,
"f": [freq0],
},
),
)
custom_src_3 = td.CustomFieldSource(
source_time=pulse,
center=(-1, 1, 0),
size=(0, 8, 8),
field_dataset=dataset_E,
)
custom_src_4 = td.CustomFieldSource(
source_time=pulse,
center=(-1, 1, 0),
size=(0, 8, 8),
field_dataset=dataset_EH,
)
sim_3 = td.Simulation(
size=sim_size,
grid_spec=grid_spec,
sources=[custom_src_3],
monitors=[mnt_xy, mnt_flux_f, mnt_flux_b],
run_time=run_time,
boundary_spec=td.BoundarySpec.all_sides(boundary=td.PML()),
)
sim_4 = td.Simulation(
size=sim_size,
grid_spec=grid_spec,
sources=[custom_src_4],
monitors=[mnt_xy, mnt_flux_f, mnt_flux_b],
run_time=run_time,
boundary_spec=td.BoundarySpec.all_sides(boundary=td.PML()),
)
[11]:
batch = web.Batch(simulations={"custom_E": sim_3, "custom_EH": sim_4}, verbose=True)
batch_results = batch.run(path_dir="data")
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- 824fe08e-7cdc-438b-a5ff-04609c1e5cc1v1'.
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- e3eccfef-14da-41c8-874b-ab7b1c780344v1'.
[10:32:57] Started working on Batch. container.py:475
[10:33:06] Maximum FlexCredit cost: 0.050 for the whole batch. container.py:479 Use 'Batch.real_cost()' to get the billed FlexCredit cost after the Batch has completed.
[10:33:25] Batch complete. container.py:522
As can be seen below, if we only provide the E field to the custom source (or only the H field), the source is not directional, but instead equal power is injected in both directions. We can make it directional by providing both fields, provided that they are set with the correct phase offset. In this example, multiplying one of the two fields by -1 will make the source inject in the backwards direction.
Note also that in the case of custom fields, the flux is not automatically normalized. This could however be done by hand if both the E and the H fields are known, for example using the Tidy3D FieldData property flux.
[12]:
fig, ax = plt.subplots(1, 2, figsize=(10, 4))
sim_data = batch_results["custom_E"]
sim_data.plot_field("field_xy", field_name="Ey", val="abs", ax=ax[0])
ax[0].set_xlim([-sim.size[0] / 2, sim.size[0] / 2])
ax[0].set_ylim([-sim.size[1] / 2, sim.size[1] / 2])
title = f"Flux: fwd={float(sim_data['flux_f'].flux):1.2e}, bck={float(sim_data['flux_b'].flux):1.2e}"
title += "\nOnly E field supplied"
ax[0].set_title(title)
sim_data = batch_results["custom_EH"]
sim_data.plot_field("field_xy", field_name="Ey", val="abs", ax=ax[1])
ax[1].set_xlim([-sim.size[0] / 2, sim.size[0] / 2])
ax[1].set_ylim([-sim.size[1] / 2, sim.size[1] / 2])
title = f"Flux: fwd={float(sim_data['flux_f'].flux):1.2e}, bck={float(sim_data['flux_b'].flux):1.2e}"
title += "\nE and H fields supplied"
ax[1].set_title(title)
plt.show()
Notes on CustomFieldSource#
Only the field components tangential to the custom source plane are needed and used in the simulation. Due to the equivalence principle, these fully define the currents that need to be injected. This is not to say that the normal components of the data (
Ex,Hxin our example) is lost or not injected. It is merely not needed as it can be uniquely obtained using the tangential components.Source data can be imported from file just as shown here, after the data is imported as a numpy array using standard numpy functions like loadtxt.
If the data is not coming from a Tidy3D simulation, the normalization is likely going to be arbitrary and the directionality of the source will likely not be perfect, even if both the
EandHfields are provided. An empty normalizing run may be needed to accurately normalize results.
Custom Current Sources#
Alternatively, one might want to inject a raw electric and magnetic current distribution within the simulation. For this, we introduce a CustomCurrentSource object.
The syntax is very similar to CustomFieldSource, except instead of a field_dataset, the source accepts a current_dataset. This dataset still contains E{x,y,z} and H{x,y,z} field components, which correspond to J and M components respectively.
There are also fewer constraints on the data requirements for CustomCurrentSource. It can be volumetric or planar without requiring tangential components.
Finally, note that the dataset is still defined w.r.t. the source center, just as in the case of the CustomFieldSource, and can then be placed anywhere in the simulation.
To demonstrate, here we make a CustomCurrentSource using our field pattern and the relations \(J = n\times H\) and \(M = -n\times E\) that should recreate the unidirectional excitation (up to the numerical resolution of our source).
[13]:
dataset_JM = td.FieldDataset(
Ey=td.ScalarFieldDataArray(
-scalar_gaussian[None, ..., None] / td.ETA_0,
coords={
"x": [0],
"y": ys,
"z": zs,
"f": [freq0],
},
),
Hz=td.ScalarFieldDataArray(
-scalar_gaussian[None, ..., None],
coords={
"x": [0],
"y": ys,
"z": zs,
"f": [freq0],
},
),
)
custom_src_JM = td.CustomCurrentSource(
source_time=pulse,
center=(-1, 1, 0),
size=(0, 8, 8),
current_dataset=dataset_JM,
)
sim_JM = sim_4.updated_copy(sources=[custom_src_JM])
[14]:
sim_data_JM = web.run(sim_JM, task_name='custom current')
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- e2320d6d-fea6-4986-8e94-133f2bb81e0bv1'.
[15]:
fig, ax = plt.subplots(1)
sim_data.plot_field("field_xy", field_name="Ey", val="abs", ax=ax)
ax.set_xlim([-sim.size[0] / 2, sim.size[0] / 2])
ax.set_ylim([-sim.size[1] / 2, sim.size[1] / 2])
ax.set_title(
f"Flux fwd={float(sim_data['flux_f'].flux):1.2e}, bck={float(sim_data['flux_b'].flux):1.2e}"
)
plt.show()
[ ]: