Aerodynamics & Electromagnetic Simulation | Flexcompute
  • Sign up for Free
  • Sign up for Free
  • Sign in
  • Talk to an Expert
Aerodynamics & Electromagnetic Simulation | Flexcompute
  • Flow360 (CFD)
    • Flow360 Solver
    • Flow360 Documentation
    • Publications
    • Technical Support
    • Contact Us
  • Tidy3D (EM)
    • Tidy3D Solver
    • Tidy3D Documentation
    • Tidy3D Knowledge Base
    • Learning FDTD through Examples
    • White Paper
    • Publications
    • Technical Support
    • Contact Us
  • Solution
    • Flow360 Solution
  • Resources
    • FDTD 101
    • Aircraft CFD 101
    • CFD Essentials
    • Newsroom
    • Blog
    • Python FDTD
    • FDTD
    • FDTD Acceleration
  • About
    • Vision
    • Story
    • Team
    • Careers
    • Contact Us
  • Login
    • Flow360
    • Tidy3D
Learning FDTD through Examples  / How to perform...
Sign up now to simulate .ipynb

How to perform near field to far field projections in Tidy3D FDTD¶

This tutorial will show you how to use field projections to obtain electromagnetic field data far away from a structure with knowledge of only near-field data.

When projecting fields, geometric approximations can be invoked to allow computing fields far away from the structure quickly and with good accuracy, but in Tidy3D we can also turn these approximations off when projecting fields at intermediate distances away, which gives a lot of flexibility.

These field projections are particularly useful for eliminating the need to simulate large regions of empty space around a structure.

In this notebook, we will

  • show how to compute projected fields on your local machine after a simulation is run, or on our servers during the simulation run.

  • show how to extract various quantities related to projected fields such as fields in different coordinate systems, power, and radar cross section.

  • demonstrate how, when far field approximations are used, the fields can dynamically be re-projected to new distances without having to run a new simulation.

  • study when geometric far field approximations should and should not be invoked, depending on the projection distance and the geometry of the structure.

  • show how to set up projections for finite-sized objects (e.g., scattering at a sphere) vs. thin but large-area structures (e.g., metasurfaces).

Table of contents¶

  1. Simulation setup
  2. Far-field projector setup
  3. Server-side far field projection
  4. Coordinate system conversion, power computation
  5. Re-projection to a new far field distance
  6. Exact field projections without making the far-field approximation
  7. Projection to a grid defined in reciprocal space
  8. Some final notes
[1]:
# standard python imports
import numpy as np
import matplotlib.pyplot as plt

# tidy3d imports
import tidy3d as td
import tidy3d.web as web

Far Field for a Uniformly Illuminated Aperture ¶

First, we will consider the simple case of an aperture in a perfect electric conductor sheet illuminated by a plane wave. The far fields in this case are known analytically, which allows for a straightforward comparison to Tidy3D's field projection functionality. We will show how to compute the far fields both on your local machine, and on the server. The geometry is shown below.

Geometry setup¶

[2]:
# size of the aperture (um)
width = 1.5
height = 2.5

# free space central wavelength (um)
wavelength = 0.75
# center frequency
f0 = td.C_0 / wavelength

# Define materials
air = td.Medium(permittivity=1)
pec = td.PECMedium()

# PEC plate thickness
thick = 0.2

# FDTD grid resolution
min_cells_per_wvl = 20

# create the PEC plate
plate = td.Structure(
    geometry=td.Box(size=[td.inf, thick, td.inf], center=[0, 0, 0]), medium=pec
)

# create the aperture in the plate
aperture = td.Structure(
    geometry=td.Box(size=[width, 1.5 * thick, height], center=[0, 0, 0]), medium=air
)

# make sure to append the aperture to the plate so that it overrides that region of the plate
geometry = [plate, aperture]

# define the boundaries as PML on all sides
boundary_spec = td.BoundarySpec.all_sides(boundary=td.PML())

# set the total domain size in x, y, and z
sim_size = [width * 2, 2, height * 2]

Source setup¶

For our incident field, we create a plane wave incident from the left, with the electric field polarized in the -z direction.

[3]:
# bandwidth in Hz
fwidth = f0 / 10.0

# time dependence of source
gaussian = td.GaussianPulse(freq0=f0, fwidth=fwidth)

# place the source to the left, propagating in the +y direction
offset_src = -0.3
source = td.PlaneWave(
    center=(0, offset_src, -0),
    size=(td.inf, 0, td.inf),
    source_time=gaussian,
    direction="+",
    pol_angle=np.pi / 2,
)

# Simulation run time
run_time = 50 / fwidth

Create monitor¶

First, we'll see how to do field projections using your machine after you've downloaded near fields from a Tidy3D simulation.

We create a surface FieldMonitor just to the right of the aperture to capture the near field data in the frequency domain.

[4]:
offset_mon = 0.3
monitor_near = td.FieldMonitor(
    center=[0, offset_mon, 0], size=[td.inf, 0, td.inf], freqs=[f0], name="near_field"
)

Create Simulation¶

Now we can put everything together and define the simulation with a simple uniform mesh, and then we'll visualize the geometry to make sure everything looks right.

[5]:
sim = td.Simulation(
    size=sim_size,
    center=[0, 0, 0],
    grid_spec=td.GridSpec.uniform(dl=wavelength / min_cells_per_wvl),
    structures=geometry,
    sources=[source],
    monitors=[monitor_near],
    run_time=run_time,
    boundary_spec=boundary_spec,
)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 3))
sim.plot(x=0, ax=ax1)
sim.plot(y=0, ax=ax2)
plt.show()

Run simulation¶

[6]:
sim_data = web.run(
    sim, task_name="aperture_1", path="data/aperture_1.hdf5", verbose=True
)
[16:32:44] Created task 'aperture_1' with task_id 'fdve-e7550611-d9c7-480a-be68-0f844034d6a4v1'.      webapi.py:139
           View task using web UI at                                                                  webapi.py:141
           'https://tidy3d.simulation.cloud/workbench?taskId=fdve-e7550611-d9c7-480a-be68-0f844034d6a              
           4v1'.                                                                                                   
Output()


[16:32:47] status = queued                                                                            webapi.py:271
Output()
[16:32:50] status = preprocess                                                                        webapi.py:265

[16:32:54] Maximum FlexCredit cost: 0.025. Use 'web.real_cost(task_id)' to get the billed FlexCredit  webapi.py:288
           cost after a simulation run.                                                                            
           starting up solver                                                                         webapi.py:292
           running solver                                                                             webapi.py:302
Output()
[16:33:00] early shutoff detected, exiting.                                                           webapi.py:316


[16:33:01] status = postprocess                                                                       webapi.py:333
Output()
[16:33:04] status = success                                                                           webapi.py:340

Output()


[16:33:06] loading SimulationData from data/aperture_1.hdf5                                           webapi.py:512

Far field points ¶

Now, we'll define the set of observation angles far away from the source at which we'd like to measure the far fields.

[7]:
# radial distance away from the origin at which to project fields
r_proj = 50 * wavelength

# theta and phi angles at which to observe fields - part of the half-space to the right
theta_proj = np.linspace(np.pi / 10, np.pi - np.pi / 10, 100)
phi_proj = np.linspace(np.pi / 10, np.pi - np.pi / 10, 100)

Now, we define a far-field monitor, FieldProjectionAngleMonitor, which stores the information regarding the far field projection grid, and then we define the object that does the actual projections, FieldProjector.

[8]:
# far field projection monitor
monitor_far = td.FieldProjectionAngleMonitor(
    center=[
        0,
        offset_mon,
        0,
    ],  # the monitor's center defined the local origin - the projection distance
    # and angles will all be measured with respect to this local origin
    size=[td.inf, 0, td.inf],
    # the size and center of any far field monitor should indicate where the *near* fields are recorded
    freqs=[f0],
    name="far_field",
    phi=list(phi_proj),
    theta=list(theta_proj),
    proj_distance=r_proj,
    far_field_approx=True,  # we leave this to its default value of 'True' because we are interested in fields sufficiently
    # far away that geometric far field approximations can be invoked to speed up the calculation
)

# helper functin to call the projector
def get_proj_fields(sim_data, monitor_near, monitor_far, pts_per_wavelength=10):
    # object that does projections is constructed using the near-field monitor, because those are the fields to be projected
    projector = td.FieldProjector.from_near_field_monitors(
        sim_data=sim_data,
        near_monitors=[monitor_near],
        normal_dirs=["+"],  # we are projecting along the + direction
        pts_per_wavelength=pts_per_wavelength,  # to speed up calculations, the fields on the near-field monitor can be downsampled to these
        # many points per wavelength (default is already 10)
    )
    return projector.project_fields(monitor_far)


# execute the projector, with the far field monitor as input, to do the projection
# let's also time this, for later use
import time

t0 = time.perf_counter()
projected_field_data = get_proj_fields(sim_data, monitor_near, monitor_far)
t1 = time.perf_counter()
proj_time = t1 - t0
Output()


Analytical solution¶

Before we plot and analyze the results, we need reference data with which to perform comparisons. In our simple aperture example, an analytical expression for the far fields is already available, so we'll simply implement the analytic formula here at the observation points of interest.

[9]:
def analytic_fields_aperture(
    proj_monitor, sim_size, aperture_height, aperture_width, r_proj
):
    """Compute the far fields analytically."""
    # in Tidy3D, the plane wave source is normalized so that a total flux of 1 is injected into the simulation domain,
    # which corresponds to an electric field strength that is inversely proportional to the square root of the in-plane domain area
    thetas_ext = np.array(proj_monitor.theta)[None, :, None, None]
    phis_ext = np.array(proj_monitor.phi)[None, None, :, None]
    f = np.array(proj_monitor.freqs)[None, None, None, :]
    E0 = np.sqrt(2.0 * td.ETA_0 / sim_size[0] / sim_size[2])
    k = 2.0 * np.pi * f / td.C_0
    ux = k * np.sin(thetas_ext) * np.cos(phis_ext) * aperture_width / 2.0
    uz = k * np.cos(thetas_ext) * aperture_height / 2.0
    Etheta = (
        -k
        / 2.0
        / np.pi
        / r_proj
        * E0
        * np.sin(thetas_ext)
        * np.exp(1j * k * r_proj)
        * aperture_height
        * aperture_width
        * np.sinc(ux / np.pi)
        * np.sinc(uz / np.pi)
    )
    Hphi = Etheta / td.ETA_0

    # for convenience, let's encapsulate the data into one of Tidy3D's native data structures designed for
    # storing far fields - this is the same format in which data will be returned when using Tidy3D's
    # 'FieldProjector', so comparisons will be easier to make
    coords = dict(
        r=np.array([r_proj]),
        theta=np.array(proj_monitor.theta),
        phi=np.array(proj_monitor.phi),
        f=np.array(proj_monitor.freqs),
    )
    Etheta_data = td.FieldProjectionAngleDataArray(Etheta, coords=coords)
    Hphi_data = td.FieldProjectionAngleDataArray(Hphi, coords=coords)
    Er_data = td.FieldProjectionAngleDataArray(np.zeros_like(Etheta), coords=coords)
    Ephi_data = td.FieldProjectionAngleDataArray(np.zeros_like(Etheta), coords=coords)
    Hr_data = td.FieldProjectionAngleDataArray(np.zeros_like(Etheta), coords=coords)
    Htheta_data = td.FieldProjectionAngleDataArray(np.zeros_like(Etheta), coords=coords)
    return td.FieldProjectionAngleData(
        monitor=proj_monitor,
        Er=Er_data,
        Etheta=Etheta_data,
        Ephi=Ephi_data,
        Hr=Hr_data,
        Htheta=Htheta_data,
        Hphi=Hphi_data,
        projection_surfaces=proj_monitor.projection_surfaces,
    )


analytic_field_data = analytic_fields_aperture(
    monitor_far, sim_size, height, width, r_proj
)

Plot and compare¶

Now we can compare the analytic fields to those computed via Tidy3D's FieldProjector, and also compute the root mean squared error between the two.

[10]:
def make_field_plot(phi, theta, vals1, vals2):
    n_plots = 2
    fig, ax = plt.subplots(1, n_plots, tight_layout=True, figsize=(8, 3.8))
    im1 = ax[0].pcolormesh(
        phi * 180 / np.pi,
        theta * 180 / np.pi,
        np.real(vals1),
        cmap="RdBu",
        shading="auto",
    )
    im2 = ax[1].pcolormesh(
        phi * 180 / np.pi,
        theta * 180 / np.pi,
        np.real(vals2),
        cmap="RdBu",
        shading="auto",
    )
    fig.colorbar(im1, ax=ax[0])
    fig.colorbar(im2, ax=ax[1])
    ax[0].set_title("Analytic")
    ax[1].set_title("Field projection")
    for _ax in ax:
        _ax.set_xlabel("$\phi$ (deg)")
        _ax.set_ylabel("$\\theta$ (deg)")


# RMSE
def rmse(array_ref, array_test):
    error = array_test - array_ref
    rmse = np.sqrt(np.mean(np.abs(error.flatten()) ** 2))
    nrmse = rmse / np.abs(np.max(array_ref.flatten()) - np.min(array_ref.flatten()))
    return nrmse


# plot Etheta
Etheta_analytic = analytic_field_data.Etheta.isel(f=0, r=0)
Etheta_proj = projected_field_data.Etheta.isel(f=0, r=0)
make_field_plot(phi_proj, theta_proj, Etheta_analytic, Etheta_proj)

# print the normalized RMSE
print(
    f"Normalized root mean squared error: {rmse(Etheta_analytic.values, Etheta_proj.values) * 100:.2f} %"
)

plt.show()
Normalized root mean squared error: 4.79 %

We obtain good agreement to analytical results. Now let's see if we can repeat this simulation but compute the far fields on the server, during the simulation run.

Server-side field projection ¶

All we have to do is provide the FieldProjectionAngleMonitor monitor as an input to the Tidy3D Simulation object as one of its monitors. Now, we no longer need to provide a separate near-field FieldMonitor - the near fields will automatically be recorded based on the size and location of the FieldProjectionAngleMonitor.

[11]:
sim2 = td.Simulation(
    size=sim_size,
    center=[0, 0, 0],
    grid_spec=td.GridSpec.uniform(dl=wavelength / min_cells_per_wvl),
    structures=geometry,
    sources=[source],
    monitors=[
        monitor_far
    ],  # just provide the far field FieldProjectionAngleMonitor as the input monitor
    run_time=run_time,
    boundary_spec=boundary_spec,
)

Run the new simulation.

[12]:
sim_data2 = web.run(
    sim2, task_name="aperture_2", path="data/aperture_2.hdf5", verbose=True
)
[16:33:09] Created task 'aperture_2' with task_id 'fdve-92acd7c3-ea08-42b3-a358-26d27de5ca57v1'.      webapi.py:139
           View task using web UI at                                                                  webapi.py:141
           'https://tidy3d.simulation.cloud/workbench?taskId=fdve-92acd7c3-ea08-42b3-a358-26d27de5ca5              
           7v1'.                                                                                                   
Output()


[16:33:11] status = queued                                                                            webapi.py:271
Output()
[16:33:12] status = preprocess                                                                        webapi.py:265

[16:33:17] Maximum FlexCredit cost: 0.025. Use 'web.real_cost(task_id)' to get the billed FlexCredit  webapi.py:288
           cost after a simulation run.                                                                            
           starting up solver                                                                         webapi.py:292
           running solver                                                                             webapi.py:302
Output()
[16:33:23] early shutoff detected, exiting.                                                           webapi.py:316


           status = postprocess                                                                       webapi.py:333
Output()
[16:33:26] status = success                                                                           webapi.py:340

Output()


[16:33:32] loading SimulationData from data/aperture_2.hdf5                                           webapi.py:512

Now the projected fields are already contained in the returned sim_data2 object - all we have to do is access it as follows, and then plot and compare to analytical results as before.

[13]:
# extract the computed projected fields
projected_field_data_server = sim_data2[monitor_far.name]

# plot Etheta
Etheta_proj_server = projected_field_data_server.Etheta.isel(f=0, r=0)
make_field_plot(phi_proj, theta_proj, Etheta_analytic, Etheta_proj_server)

# print the normalized RMSE
print(
    f"Normalized root mean squared error: {rmse(Etheta_analytic.values, Etheta_proj_server.values) * 100:.2f} %"
)

plt.show()
Normalized root mean squared error: 4.78 %

We obtain nearly identical results, except that they are computed much faster on our servers. Note also that in some cases, the server-side computations may be slightly more accurate than client-side ones, because on the server, the near fields are not downsampled at all.

To see the performance gains of using server-side computations, let's compare the time taken in each case.

[14]:
# use the simulation log to find the time taken for server-side computations
server_time = float(
    sim_data2.log.split("Field projection time (s):    ", 1)[1].split("\n", 1)[0]
)
print(f"Client-side field projection took {proj_time:.2f} s")
print(f"Server-side field projection took {server_time:.2f} s")
Client-side field projection took 1.65 s
Server-side field projection took 0.60 s

As expected, the server computes far fields faster than the local CPU-based computation, though it's a relatively small gain in this case. The gains in computation time are expected to be greater for larger and more complex setups.

Other far field quantities and coordinate systems ¶

So far, we've been looking at the electric field in spherical coordinates. However, we can also look at the fields in other coordinate systems, e.g., E_x, E_y, E_z, and the radiated power, as follows:

[15]:
def make_cart_plot(phi, theta, vals1, vals2, vals3):
    n_plots = 3
    fig, ax = plt.subplots(1, n_plots, tight_layout=True, figsize=(12.6, 3.8))
    im1 = ax[0].pcolormesh(
        phi * 180 / np.pi,
        theta * 180 / np.pi,
        np.real(vals1),
        cmap="RdBu",
        shading="auto",
    )
    im2 = ax[1].pcolormesh(
        phi * 180 / np.pi,
        theta * 180 / np.pi,
        np.real(vals2),
        cmap="RdBu",
        shading="auto",
    )
    im3 = ax[2].pcolormesh(
        phi * 180 / np.pi,
        theta * 180 / np.pi,
        np.real(vals3),
        cmap="RdBu",
        shading="auto",
    )
    fig.colorbar(im1, ax=ax[0])
    fig.colorbar(im2, ax=ax[1])
    fig.colorbar(im3, ax=ax[2])
    ax[0].set_title("Ex")
    ax[1].set_title("Ey")
    ax[2].set_title("Ez")
    for _ax in ax:
        _ax.set_xlabel("$\phi$ (deg)")
        _ax.set_ylabel("$\\theta$ (deg)")


# get the fields in Cartesian coordinates from the projected data we already computed above
fields_cartesian = projected_field_data.fields_cartesian.isel(f=0, r=0)

# plot Ex, Ey, Ez
make_cart_plot(
    phi_proj, theta_proj, fields_cartesian.Ex, fields_cartesian.Ey, fields_cartesian.Ez
)

# get the power
power = projected_field_data.power.isel(f=0, r=0)

# plot the power
fig, ax = plt.subplots(1, 1, tight_layout=True, figsize=(4.3, 3.8))
im = ax.pcolormesh(
    phi_proj * 180 / np.pi,
    theta_proj * 180 / np.pi,
    power,
    cmap="inferno",
    shading="auto",
)
fig.colorbar(im, ax=ax)
_ = ax.set_title("Power")
_ = ax.set_xlabel("$\phi$ (deg)")
_ = ax.set_ylabel("$\\theta$ (deg)")

plt.show()

Re-projection to a different distance ¶

We can re-project the already-computed far fields to a different distance away from the structure - we neither need to run another simulation nor re-run the FieldProjector. Instead, the fields can simply be renormalized as shown below.

Note that by default, if no proj_distance was provided in the FieldProjectionAngleMonitor, the fields are projected to a distance of 1m.

[16]:
# new projection distance
r_proj_new = 20 * wavelength

# re-project our far field data above to this new distance
reprojected_field_data = projected_field_data.renormalize_fields(r_proj_new)

# now all the fields stored in 'projected_field_data' correspond to this new distance
# compare to the analytical fields at this new distance
analytic_field_data_new = analytic_fields_aperture(
    monitor_far, sim_size, height, width, r_proj_new
)

# plot Etheta
Etheta_analytic = analytic_field_data_new.Etheta.isel(f=0, r=0)
Etheta_proj = reprojected_field_data.Etheta.isel(f=0, r=0)
make_field_plot(phi_proj, theta_proj, Etheta_analytic, Etheta_proj)

# print the normalized RMSE
print(
    f"Normalized root mean squared error: {rmse(Etheta_analytic.values, Etheta_proj.values) * 100:.2f} %"
)

plt.show()
Normalized root mean squared error: 4.79 %

More accurate field projections ¶

In the field projections used above, the far field approximation is used: it is assumed that the fields are measured at a distance much greater than the size of our simulation in the transverse direction. Accordingly, geometric approximations are invoked, and any quantity whose magnitude drops off as 1/r^2 or faster is ignored. The advantages of these approximations are:

  • the projections are computed relatively fast
  • the projections are cast in a simple mathematical form which allows re-projecting the fields to different distance without the need to re-run a simulation or to re-run the FieldProjector.

However, in some cases we may want to project to intermediate distances where the far field approximation is no longer valid. Tidy3D's field projection functionality allows doing this very easily: simply flip a switch when defining the FieldProjectionAngleMonitor! The resulting computations will be a bit slower, but the results will be significantly more accurate.

Note: when the far field approximations are turned off, we can no longer simply use renormalize_fields to re-project the fields at a new distance. Instead, we would need to re-run the FieldProjector.

Below, we will demonstrate this feature by looking at fields only a few wavelengths away from the aperture. Note that our analytical results also made far field approximations, so here we'll make our simulation domain a bit larger and measure the actual fields on a monitor, so that we can compare these actual fields to those computed by the FieldProjector.

Also, this time we'll use the FieldProjectionCartesianMonitor, which is the counterpart to the FieldProjectionAngleMonitor where the observation grid is defined in Cartesian coordinates, not angles.

[17]:
# project fields only a few wavelengths away from the aperture
r_proj_intermediate = 4 * wavelength

# create a field monitor to measure these fields at the intermediate projection distance,
# so that we have something to which we can compare the 'FieldProjector' results
monitor_intermediate = td.FieldMonitor(
    center=[0, offset_mon + r_proj_intermediate, 0],
    size=[td.inf, 0, td.inf],
    freqs=[f0],
    name="inter_field",
)

# make a larger simulation along y to accommodate the plane at which the intermediate fields need to be measured
shift = 1.2 * r_proj_intermediate
sim_size3 = [sim_size[0], sim_size[1] + shift, sim_size[2]]
# move the sim center
sim_center = [0, (sim_size[1] + shift) / 2 - sim_size[1] / 2, 0]
sim3 = td.Simulation(
    size=sim_size3,
    center=sim_center,
    grid_spec=td.GridSpec.uniform(dl=wavelength / min_cells_per_wvl),
    structures=geometry,
    sources=[source],
    monitors=[
        monitor_near,
        monitor_intermediate,
    ],  # provide both near field and intermediate field monitors
    run_time=run_time,
    boundary_spec=boundary_spec,
)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 3))
sim3.plot(x=0, ax=ax1)
sim3.plot(y=0, ax=ax2)
plt.show()

Run the new simulation.

[18]:
sim_data3 = web.run(
    sim3, task_name="aperture_3", path="data/aperture_3.hdf5", verbose=True
)
[16:33:33] Created task 'aperture_3' with task_id 'fdve-3cdf602c-25c6-44d7-b01d-95b5386f2c2bv1'.      webapi.py:139
           View task using web UI at                                                                  webapi.py:141
           'https://tidy3d.simulation.cloud/workbench?taskId=fdve-3cdf602c-25c6-44d7-b01d-95b5386f2c2              
           bv1'.                                                                                                   
Output()


[16:33:35] status = queued                                                                            webapi.py:271
Output()
[16:33:37] status = preprocess                                                                        webapi.py:265

[16:33:42] Maximum FlexCredit cost: 0.025. Use 'web.real_cost(task_id)' to get the billed FlexCredit  webapi.py:288
           cost after a simulation run.                                                                            
           starting up solver                                                                         webapi.py:292
           running solver                                                                             webapi.py:302
Output()
[16:33:48] early shutoff detected, exiting.                                                           webapi.py:316


           status = postprocess                                                                       webapi.py:333
Output()
[16:33:52] status = success                                                                           webapi.py:340

Output()


[16:33:54] loading SimulationData from data/aperture_3.hdf5                                           webapi.py:512

Now let's create the FieldProjectionCartesianMonitor, this time turning off the far field approximations, and then run the FieldProjector again.

The FieldProjectionCartesianMonitor's xy observation grid is defined in a local coordinate system whose z axis points in the direction along which we want to project fields, in this case the +y axis. The mapping between local and global coordinates is as follows:

  • proj_axis=0: local x = global y, local y = global z
  • proj_axis=1: local x = global x, local y = global z
  • proj_axis=2: local x = global x, local y = global y
[19]:
# make the projection monitor which projects fields without approximations
xs = np.linspace(-sim_size[0] / 2, sim_size[0] / 2, 100)
ys = np.linspace(-sim_size[1] / 2, sim_size[1] / 2, 100)
monitor_intermediate_proj = td.FieldProjectionCartesianMonitor(
    center=[
        0,
        offset_mon,
        0,
    ],  # the monitor's center defined the local origin - the projection distance
    # and angles will all be measured with respect to this local origin
    size=[td.inf, 0, td.inf],
    freqs=[f0],
    name="inter_field_proj",
    proj_axis=1,  # because we're projecting along the +y axis
    x=list(xs),  # local x coordinates - corresponds to global x in this case
    y=list(ys),  # local y coordinates - corresponds to global z in this case
    proj_distance=r_proj_intermediate,
    far_field_approx=False,  # turn off the far-field approximation (is 'True' by default)
)

# execute the projector, with the projection field monitor as input, to do the projection
# let's also time this, for later use
import time

t0 = time.perf_counter()
projected_field_data_noapprox = get_proj_fields(
    sim_data3, monitor_near, monitor_intermediate_proj
)
t1 = time.perf_counter()
proj_time_new = t1 - t0
Output()


Now let's compare the following three results:

  • Directly-measured fields at the projection distance
  • Projected fields with approximations turned off
  • Projected fields with approximations turned on (just to compare the accuracy)
[20]:
# Compute projected fields *with* far field approximations, to facilitate an accuracy comparison
monitor_intermediate_proj_approx = td.FieldProjectionCartesianMonitor(
    center=[
        0,
        offset_mon,
        0,
    ],  # the monitor's center defined the local origin - the projection distance
    # and angles will all be measured with respect to this local origin
    size=[td.inf, 0, td.inf],
    freqs=[f0],
    name="inter_field_proj_approx",
    proj_axis=1,  # because we're projecting along the +y axis
    x=list(xs),  # local x coordinates - corresponds to global x in this case
    y=list(ys),  # local y coordinates - corresponds to global z in this case
    proj_distance=r_proj_intermediate,
    far_field_approx=True,  # turn on the far-field approximation
)

# execute the projector, with the projection field monitor as input, to do the projection
# let's also time this, for later use
import time

t0 = time.perf_counter()
projected_field_data_approx = get_proj_fields(
    sim_data3, monitor_near, monitor_intermediate_proj_approx
)
t1 = time.perf_counter()
proj_time_new_approx = t1 - t0
Output()


[21]:
# let's see how long this took compared to the previous case when the approximations were turned on
print(
    f"Client-side field projection *with approximations on* took {proj_time_new_approx:.2f} s"
)
print(
    f"Client-side field projection *with approximations off* took {proj_time_new:.2f} s"
)
Client-side field projection *with approximations on* took 3.68 s
Client-side field projection *with approximations off* took 12.03 s

As expected, when the approximations are turned off, the projections take longer. Now let's see if it was worth it!

[22]:
# Helper function to plot fields
def make_cart_plot(phi, theta, vals1, vals2, vals3):
    n_plots = 3
    fig, ax = plt.subplots(1, n_plots, tight_layout=True, figsize=(9, 3))
    im1 = ax[0].pcolormesh(ys, xs, np.real(vals1), cmap="RdBu", shading="auto")
    im2 = ax[1].pcolormesh(ys, xs, np.real(vals2), cmap="RdBu", shading="auto")
    im3 = ax[2].pcolormesh(ys, xs, np.real(vals3), cmap="RdBu", shading="auto")
    fig.colorbar(im1, ax=ax[0])
    fig.colorbar(im2, ax=ax[1])
    fig.colorbar(im3, ax=ax[2])
    ax[0].set_title("Ex")
    ax[1].set_title("Ey")
    ax[2].set_title("Ez")
    for _ax in ax:
        _ax.set_xlabel("$y$ (micron)")
        _ax.set_ylabel("$x$ (micron)")
[23]:
# plot the actual measured fields
fields_meas = sim_data3[monitor_intermediate.name].colocate(x=xs, z=ys)
make_cart_plot(
    ys,
    xs,
    fields_meas.Ex.isel(f=0, y=0),
    fields_meas.Ey.isel(f=0, y=0),
    fields_meas.Ez.isel(f=0, y=0),
)
plt.suptitle("Measured fields")

# projected field without approximations - get them in Cartesian coords
fields_proj_noapprox = projected_field_data_noapprox.fields_cartesian
make_cart_plot(
    ys,
    xs,
    fields_proj_noapprox.Ex.isel(f=0, y=0),
    fields_proj_noapprox.Ey.isel(f=0, y=0),
    fields_proj_noapprox.Ez.isel(f=0, y=0),
)
plt.suptitle("Projected, no approximations")

# projected field with approximations - get them in Cartesian coords
fields_proj_approx = projected_field_data_approx.fields_cartesian
make_cart_plot(
    ys,
    xs,
    fields_proj_approx.Ex.isel(f=0, y=0),
    fields_proj_approx.Ey.isel(f=0, y=0),
    fields_proj_approx.Ez.isel(f=0, y=0),
)
_ = plt.suptitle("Projected, with far field approximations")

# RMSE
Emag_meas = np.sqrt(
    np.abs(fields_meas.Ex) ** 2
    + np.abs(fields_meas.Ey) ** 2
    + np.abs(fields_meas.Ez) ** 2
)
Emag_proj_noapprox = np.sqrt(
    np.abs(fields_proj_noapprox.Ex) ** 2
    + np.abs(fields_proj_noapprox.Ey) ** 2
    + np.abs(fields_proj_noapprox.Ez) ** 2
)
Emag_proj_approx = np.sqrt(
    np.abs(fields_proj_approx.Ex) ** 2
    + np.abs(fields_proj_approx.Ey) ** 2
    + np.abs(fields_proj_approx.Ez) ** 2
)
print(
    f"Normalized RMSE for |E|, no far field approximation: {rmse(Emag_meas.values, Emag_proj_noapprox.values) * 100:.2f} %"
)
print(
    f"Normalized RMSE for |E|, with far field approximation: {rmse(Emag_meas.values, Emag_proj_approx.values) * 100:.2f} %"
)

plt.show()
Normalized RMSE for |E|, no far field approximation: 0.64 %
Normalized RMSE for |E|, with far field approximation: 24.03 %

Without approximations, the projected fields match the measured ones extremely well! Instead, when approximations are used, the match is very poor. Thus, the accurate field projections can be extremely useful when the projection distance is not large compared to the structure size, but one still wants to avoid simulating all the empty space around the structure.

We should also note that this more accurate version of field projections can be run on the server in exactly the same way as before: just supply the projection monitor with its far_field_approx field set to False into the simulation's list of monitors as before. Everything else remains exactly the same, as shown below.

[24]:
sim4 = td.Simulation(
    size=sim_size,
    center=[0, 0, 0],
    grid_spec=td.GridSpec.uniform(dl=wavelength / min_cells_per_wvl),
    structures=geometry,
    sources=[source],
    monitors=[monitor_intermediate_proj],  # only need to supply the projection monitor
    run_time=run_time,
    boundary_spec=boundary_spec,
)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 3))
sim4.plot(x=0, ax=ax1)
sim4.plot(y=0, ax=ax2)
plt.show()
[25]:
# run the simulation
sim_data4 = web.run(
    sim4, task_name="aperture_4", path="data/aperture_4.hdf5", verbose=True
)
[16:34:11] Created task 'aperture_4' with task_id 'fdve-ca4f3cbb-15ee-4295-87b7-0d11590a73a7v1'.      webapi.py:139
           View task using web UI at                                                                  webapi.py:141
           'https://tidy3d.simulation.cloud/workbench?taskId=fdve-ca4f3cbb-15ee-4295-87b7-0d11590a73a              
           7v1'.                                                                                                   
Output()


[16:34:13] status = queued                                                                            webapi.py:271
Output()
[16:34:16] status = preprocess                                                                        webapi.py:265

[16:34:23] Maximum FlexCredit cost: 0.025. Use 'web.real_cost(task_id)' to get the billed FlexCredit  webapi.py:288
           cost after a simulation run.                                                                            
           starting up solver                                                                         webapi.py:292
           running solver                                                                             webapi.py:302
Output()
[16:34:30] early shutoff detected, exiting.                                                           webapi.py:316


           status = postprocess                                                                       webapi.py:333
Output()
[16:34:35] status = success                                                                           webapi.py:340

Output()


[16:34:40] loading SimulationData from data/aperture_4.hdf5                                           webapi.py:512
[26]:
# extract the projected fields as before and plot them
projected_field_data_server = sim_data4[monitor_intermediate_proj.name]

# plot the actual measured fields from the previous simulation
fields_meas = sim_data3[monitor_intermediate.name].colocate(x=xs, z=ys)
make_cart_plot(
    ys,
    xs,
    fields_meas.Ex.isel(f=0, y=0),
    fields_meas.Ey.isel(f=0, y=0),
    fields_meas.Ez.isel(f=0, y=0),
)
plt.suptitle("Measured fields")

# projected field without approximations computed on the server
fields_proj_noapprox = projected_field_data_server.fields_cartesian
make_cart_plot(
    ys,
    xs,
    fields_proj_noapprox.Ex.isel(f=0, y=0),
    fields_proj_noapprox.Ey.isel(f=0, y=0),
    fields_proj_noapprox.Ez.isel(f=0, y=0),
)
plt.suptitle("Projected, no approximations, computed on the server")

# RMSE
Emag_proj_server = np.sqrt(
    np.abs(fields_proj_noapprox.Ex) ** 2
    + np.abs(fields_proj_noapprox.Ey) ** 2
    + np.abs(fields_proj_noapprox.Ez) ** 2
)
print(
    f"Normalized RMSE for |E|, no far field approximation, computed on the server: {rmse(Emag_meas.values, Emag_proj_server.values) * 100:.2f} %\n"
)

# use the simulation log to find the time taken for server-side computations
server_time = float(
    sim_data4.log.split("Field projection time (s):    ", 1)[1].split("\n", 1)[0]
)
print(
    f"Client-side field projection *without approximations* took {proj_time_new:.2f} s"
)
print(f"Server-side field projection *without approximations* took {server_time:.2f} s")

plt.show()
Normalized RMSE for |E|, no far field approximation, computed on the server: 0.91 %

Client-side field projection *without approximations* took 12.03 s
Server-side field projection *without approximations* took 1.04 s