# Importing custom source data#

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.

[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)
src_pos = -4
gaussian_source = td.GaussianBeam(
source_time=pulse,
center=(src_pos, 0, 0),
size=(0, 10, 10),
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 use a monitor with a slightly nonzero size in the propagation direction. This is because of the details associated with the Yee grid used in the FDTD algorithm. 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 related to the numerical grid. 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], 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], name="yz_nonzero_size_x"
)


[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)


[14:27:01] Created task 'free space gaussian' with task_id                      webapi.py:139
'fdve-2ff9d5bf-36ed-49d4-a877-6b8c026fa524v1'.

[14:27:02] status = queued                                                      webapi.py:269

[14:27:06] status = preprocess                                                  webapi.py:263

[14:27:12] Maximum FlexCredit cost: 0.025. Use 'web.real_cost(task_id)' to get    webapi.py:286
the billed FlexCredit cost after a simulation run.

           starting up solver                                                   webapi.py:290

           running solver                                                       webapi.py:300

           status = postprocess                                                 webapi.py:330

[14:27:27] status = success                                                     webapi.py:337

[14:27:28] loading SimulationData from simulation_data.hdf5                     webapi.py:512

[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, 0, 0))
custom_src_2 = sim_data["yz_nonzero_size_x"].to_source(
source_time=pulse,
center=(0, 0, 0),
size=(
0,
8,
8,
),  # by default taken from the data, but it must be reset to zero along x
)
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()),
)


[9]:

batch = web.Batch(simulations={"zero_x": sim_1, "nonzero_x": sim_2}, verbose=True)
batch_results = batch.run(path_dir="data")


[14:28:46] Created task 'zero_x' with task_id                                   webapi.py:139
'fdve-36edde23-9270-4133-9a85-33fe1c7e09f9v1'.

[14:28:47] Created task 'nonzero_x' with task_id                                webapi.py:139
'fdve-5be0ba3d-3d6b-479c-81c2-b783057403a5v1'.

[14:28:48] Started working on Batch.                                         container.py:402

[14:31:09] Batch complete.                                                   container.py:436


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.

[10]:

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()


[14:31:10] loading SimulationData from                                          webapi.py:512
data/fdve-36edde23-9270-4133-9a85-33fe1c7e09f9v1.hdf5

[14:31:10] loading SimulationData from                                          webapi.py:512
data/fdve-5be0ba3d-3d6b-479c-81c2-b783057403a5v1.hdf5


## 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.

[11]:

# 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=(0, 0, 0),
size=(0, 8, 8),
field_dataset=dataset_E,
)
custom_src_4 = td.CustomFieldSource(
source_time=pulse,
center=(0, 0, 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()),
)


[12]:

batch = web.Batch(simulations={"custom_E": sim_3, "custom_EH": sim_4}, verbose=True)
batch_results = batch.run(path_dir="data")


[14:31:11] Created task 'custom_E' with task_id                                 webapi.py:139
'fdve-c7858f32-de82-4896-90a5-15b966f6976cv1'.

[14:31:12] Created task 'custom_EH' with task_id                                webapi.py:139
'fdve-2184dd60-be40-443e-bb6d-53a0f7857e8ev1'.

[14:31:14] Started working on Batch.                                         container.py:402

[14:31:32] Batch complete.                                                   container.py:436


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.

[13]:

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()


[14:31:33] loading SimulationData from                                          webapi.py:512
data/fdve-c7858f32-de82-4896-90a5-15b966f6976cv1.hdf5

[14:31:34] loading SimulationData from                                          webapi.py:512
data/fdve-2184dd60-be40-443e-bb6d-53a0f7857e8ev1.hdf5


## Final notes#

• 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, Hx in 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 E and H fields are provided. An empty normalizing run may be needed to accurately normalize results.

[ ]: