Importing custom source data#

Run this notebook in your browser using Binder.

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

[16:42:16] INFO     Using client version: 1.8.0                               __init__.py:112

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

<AxesSubplot: title={'center': 'cross section at z=0.00'}, xlabel='x', ylabel='y'>
<Figure size 432x288 with 1 Axes>
../_images/notebooks_CustomFieldSource_6_2.png
[5]:
sim_data = web.run(sim, task_name="free space gaussian")

[16:42:17] INFO     Using Tidy3D credentials from stored file.                     auth.py:70
[16:42:18] INFO     Authentication successful.                                     auth.py:30
           INFO     Created task 'free space gaussian' with task_id             webapi.py:120
                    '2fdce2d2-30ac-4cfc-976d-193efc544086'.                                  
[16:42:20] INFO     Maximum flex unit cost: 0.03                                webapi.py:248
           INFO     status = queued                                             webapi.py:257
[16:42:23] INFO     status = preprocess                                         webapi.py:269
[16:42:26] INFO     starting up solver                                          webapi.py:273
[16:42:32] INFO     running solver                                              webapi.py:279
[16:42:35] INFO     early shutoff detected, exiting.                            webapi.py:290
           INFO     status = postprocess                                        webapi.py:296
[16:42:37] INFO     status = success                                            webapi.py:302
[16:42:38] INFO     downloading file "output/monitor_data.hdf5" to              webapi.py:585
                    "simulation_data.hdf5"                                                   
[16:42:39] INFO     loading SimulationData from simulation_data.hdf5            webapi.py:407
[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()

<Figure size 432x288 with 2 Axes>
../_images/notebooks_CustomFieldSource_8_1.png

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

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

[16:42:40] INFO     Created task 'zero_x' with task_id                          webapi.py:120
                    'adab473c-eea3-4a24-946e-301a08e7d3a7'.                                  
[16:42:42] INFO     Created task 'nonzero_x' with task_id                       webapi.py:120
                    'a075c71c-7915-4d73-ae07-973074606ad8'.                                  
[16:42:44] Started working on Batch.                                         container.py:361
[16:43:02] Batch complete.                                                   container.py:382

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)

[16:43:03] INFO     downloading file "output/monitor_data.hdf5" to              webapi.py:585
                    "data/adab473c-eea3-4a24-946e-301a08e7d3a7.hdf5"                         
[16:43:04] INFO     loading SimulationData from                                 webapi.py:407
                    data/adab473c-eea3-4a24-946e-301a08e7d3a7.hdf5                           
           INFO     downloading file "output/monitor_data.hdf5" to              webapi.py:585
                    "data/a075c71c-7915-4d73-ae07-973074606ad8.hdf5"                         
[16:43:05] INFO     loading SimulationData from                                 webapi.py:407
                    data/a075c71c-7915-4d73-ae07-973074606ad8.hdf5                           
Text(0.5, 1.0, 'Flux: fwd=1.00e+00, bck=-1.89e-06\nnonzero-size FieldData along x')
<Figure size 720x288 with 4 Axes>
../_images/notebooks_CustomFieldSource_13_12.png

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.

[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=(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()),
)

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

[16:43:06] INFO     Created task 'custom_E' with task_id                        webapi.py:120
                    '66368a70-313c-44b9-bd3b-5249e235d891'.                                  
[16:43:07] INFO     Created task 'custom_EH' with task_id                       webapi.py:120
                    '14ec7668-3ec0-424d-b529-c8bdeffbdf2e'.                                  
[16:43:09] Started working on Batch.                                         container.py:361
[16:43:26] Batch complete.                                                   container.py:382

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)

[16:43:26] INFO     downloading file "output/monitor_data.hdf5" to              webapi.py:585
                    "data/66368a70-313c-44b9-bd3b-5249e235d891.hdf5"                         
[16:43:27] INFO     loading SimulationData from                                 webapi.py:407
                    data/66368a70-313c-44b9-bd3b-5249e235d891.hdf5                           
           INFO     downloading file "output/monitor_data.hdf5" to              webapi.py:585
                    "data/14ec7668-3ec0-424d-b529-c8bdeffbdf2e.hdf5"                         
[16:43:28] INFO     loading SimulationData from                                 webapi.py:407
                    data/14ec7668-3ec0-424d-b529-c8bdeffbdf2e.hdf5                           
Text(0.5, 1.0, 'Flux: fwd=8.53e-03, bck=-2.10e-04\nE and H fields supplied')
<Figure size 720x288 with 4 Axes>
../_images/notebooks_CustomFieldSource_18_12.png

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.

[ ]: