Defining fully anisotropic materials#

Tidy3D’s capabilities include modeling of non-dispersive fully anisotropic materials. In this tutorial we explain how to set up and run a simulation containing such materials. Specifically, we will consider scattering of a plane wave from fully anisotropic dielectric spheres and compare results to an equivalent setup containing only diagonally anisotropic materials.

Simulation setup#

[1]:
# standard python imports
import numpy as np
import xarray as xr
import matplotlib.pylab as plt

# tidy3D import
import tidy3d as td
from tidy3d.constants import C_0
import tidy3d.web as web

We start with defining the simulation domain as a cube with side length of 1.6 um, and selecting 11 sample frequencies in the vicinity of the wavelength of interest \(\lambda = 1\) um.

[2]:
sim_size = (1.6, 1.6, 1.6)

wvl_um = 1
freq0 = td.C_0 / wvl_um
fwidth = freq0 / 10
num_freqs = 11
freqs = np.linspace(freq0 - fwidth / 2, freq0 + fwidth / 2, num_freqs)
find0 = num_freqs // 2

run_time = 10 / fwidth

Next we define main permittivity and conductivity components of medium we would like to simulate.

[3]:
permittivity_main = (15, 10, 5)
conductivity_main = (0.01, 0.02, 0.03)

We will create 3 different anisotropic media to use and compare in 3 different simulations. First, a diagonally anisotropic medium that will be used in a reference simulation is created using the familiar AnisotropicMedium class.

[4]:
# diagonally anisotropic material
medium_diag = td.AnisotropicMedium(
    xx=td.Medium(permittivity=permittivity_main[0], conductivity=conductivity_main[0]),
    yy=td.Medium(permittivity=permittivity_main[1], conductivity=conductivity_main[1]),
    zz=td.Medium(permittivity=permittivity_main[2], conductivity=conductivity_main[2]),
)

Next, we would like to create two fully anisotropic media obtained as rotations of the above diagonally anisotropic medium around x and z axes, respectively. Fully anisotropic materials are created using class FullyAnisotropicMedium, which requires full permittivity and conductivity tensor during initialization. Tidy3D provides a convenience class RotationAroundAxis that can be used for transformation of material property tensors. Using it we can define a fully anisotropic medium in two ways. First, one can use this class directly to perform rotation of material property tensors, and then use the resulting tensors to initialize a FullyAnisotropicMedium as following.

[5]:
# Fully anisotropic media obtained as a rotation around x axis by angle phi
phi = np.pi / 3

# create diagonal tensor
permittivity_phi = np.diag(permittivity_main)
conductivity_phi = np.diag(conductivity_main)

# create rotation around axis
rotation_around_x = td.RotationAroundAxis(axis=(1, 0, 0), angle=phi)

# transform material property tensors
permittivity_phi = rotation_around_x.rotate_tensor(permittivity_phi)
conductivity_phi = rotation_around_x.rotate_tensor(conductivity_phi)

# define a fully anisotropic medium
medium_phi = td.FullyAnisotropicMedium(
    permittivity=permittivity_phi, conductivity=conductivity_phi
)

The alternative way is to use a shortcut class method FullyAnisotropicMedium.from_diagonal() that performs the necessary tensor transformations internally. Note that this class method follows the signature of AnisotropicMedium class, that is, it expects main components of material properties defined as Medium’s.

[6]:
# Fully anisotropic media obtained as a rotation around z axis by angle theta
theta = np.pi / 7

# create rotation around axis
rotation_around_z = td.RotationAroundAxis(axis=(0, 0, 1), angle=theta)

# define a fully anisotropic medium as a transformation of diagonally anisotropic medium
medium_theta = td.FullyAnisotropicMedium.from_diagonal(
    xx=td.Medium(permittivity=permittivity_main[0], conductivity=conductivity_main[0]),
    yy=td.Medium(permittivity=permittivity_main[1], conductivity=conductivity_main[1]),
    zz=td.Medium(permittivity=permittivity_main[2], conductivity=conductivity_main[2]),
    rotation=rotation_around_z,
)

One can confirm that the three created media have identical main components by visualizing their complex refractive indexes using the built-in function .plot(). Note that in case of fully anisotropic materials the direction of each main component is provided in figure legend as they are no longer aligned with Cartesian directions.

[7]:
# Visualize material properties
fig, ax = plt.subplots(1, 3, figsize=(14, 4))
medium_diag.plot(freqs, ax=ax[0])
medium_phi.plot(freqs, ax=ax[1])
medium_theta.plot(freqs, ax=ax[2])
plt.show()
../_images/notebooks_FullyAnisotropic_14_0.png

Create three spheres of the same size but made of different anisotropic materials.

[8]:
sphere_diag = td.Structure(geometry=td.Sphere(center=(0, 0, 0), radius=0.3), medium=medium_diag)
sphere_phi = sphere_diag.updated_copy(medium=medium_phi)
sphere_theta = sphere_diag.updated_copy(medium=medium_theta)

To set up three equivalent simulations we use TFSF sources injecting plane waves whose propagation fronts are rotated in the same way as the material properties of the three anisotropic media.

[9]:
source_diag = td.TFSF(
    center=(0, 0, 0),
    size=(1.1, 1.1, 1.1),
    source_time=td.GaussianPulse(freq0=freq0, fwidth=fwidth),
    injection_axis=0,
    direction="+",
    name="tfsf",
    angle_phi=0,
    angle_theta=0,
    pol_angle=0,
)

source_phi = source_diag.updated_copy(angle_phi=phi)
source_theta = source_diag.updated_copy(angle_theta=theta)

For a quantitative comparison of results of the three equivalent simulation we monitor total scattered flux from the sphere and field distribution in two orthogonal planes.

[10]:
mnt_flux = td.FluxMonitor(center=(0, 0, 0), size=(1.3, 1.3, 1.3), freqs=freqs, name="flux")
mnt_field_yz = td.FieldMonitor(center=(0, 0, 0), size=(0, td.inf, td.inf), freqs=freqs, name="field_yz")
mnt_field_xy = td.FieldMonitor(center=(0, 0, 0), size=(td.inf, td.inf, 0), freqs=freqs, name="field_xy")

PML boundary conditions are placed on each side of the simulation domain and a simple uniform grid is used for grid discretization.

[11]:
boundary_spec = td.BoundarySpec.all_sides(boundary=td.PML())
grid_spec = td.GridSpec.uniform(dl=0.005)

Combining all the components defined above we create three equivalent simulations which are rotations of each other.

[12]:
sim_diag = td.Simulation(
    center=(0, 0, 0),
    size=sim_size,
    grid_spec=grid_spec,
    structures=[sphere_diag],
    sources=[source_diag],
    monitors=[mnt_flux, mnt_field_yz, mnt_field_xy],
    run_time=run_time,
    boundary_spec=boundary_spec,
)

sim_phi = sim_diag.updated_copy(structures=[sphere_phi], sources=[source_phi])
sim_theta = sim_diag.updated_copy(structures=[sphere_theta], sources=[source_theta])

Visualize simulation setups for confirmation.

[13]:
fig, ax = plt.subplots(1, 3, figsize=(10, 3))
sim_diag.plot(z=0, ax=ax[0])
sim_phi.plot(z=0, ax=ax[1])
sim_theta.plot(z=0, ax=ax[2])
plt.tight_layout()
plt.show()
../_images/notebooks_FullyAnisotropic_26_0.png

Results#

Submit and run simulations on the server.

[14]:
sim_data_diag = web.run(simulation=sim_diag, task_name="fully_anisotropic_diag", path="data/simulation_data_diag.hdf5", verbose=True)
sim_data_phi = web.run(simulation=sim_phi, task_name="fully_anisotropic_phi", path="data/simulation_data_phi.hdf5", verbose=True)
sim_data_theta = web.run(simulation=sim_theta, task_name="fully_anisotropic_theta", path="data/simulation_data_theta.hdf5", verbose=True)
[16:58:59] Created task 'fully_anisotropic_diag' with task_id                                         webapi.py:186
           'fdve-78d231c2-61bf-4a03-b57b-7f61b6e14d18v1'.                                                          
[16:59:02] status = queued                                                                            webapi.py:321
[16:59:09] status = preprocess                                                                        webapi.py:315
[16:59:16] Maximum FlexCredit cost: 0.436. Use 'web.real_cost(task_id)' to get the billed FlexCredit  webapi.py:338
           cost after a simulation run.                                                                            
           starting up solver                                                                         webapi.py:342
           running solver                                                                             webapi.py:352
[16:59:59] early shutoff detected, exiting.                                                           webapi.py:366
           status = postprocess                                                                       webapi.py:383
[17:00:25] status = success                                                                           webapi.py:390
[17:00:43] loading SimulationData from data/simulation_data_diag.hdf5                                 webapi.py:568
[17:00:43] Created task 'fully_anisotropic_phi' with task_id                                          webapi.py:186
           'fdve-23947891-340a-4ea7-8540-86162eeb55d9v1'.                                                          
[17:00:45] status = queued                                                                            webapi.py:321
[17:00:52] status = preprocess                                                                        webapi.py:315
[17:00:58] Maximum FlexCredit cost: 0.706. Use 'web.real_cost(task_id)' to get the billed FlexCredit  webapi.py:338
           cost after a simulation run.                                                                            
           starting up solver                                                                         webapi.py:342
           running solver                                                                             webapi.py:352
[17:01:43] early shutoff detected, exiting.                                                           webapi.py:366
           status = postprocess                                                                       webapi.py:383
[17:02:09] status = success                                                                           webapi.py:390
[17:02:27] loading SimulationData from data/simulation_data_phi.hdf5                                  webapi.py:568
[17:02:28] Created task 'fully_anisotropic_theta' with task_id                                        webapi.py:186
           'fdve-169281b2-d8df-4764-a7eb-1e815dd6b28cv1'.                                                          
[17:02:29] status = queued                                                                            webapi.py:321
[17:02:36] status = preprocess                                                                        webapi.py:315
[17:02:43] Maximum FlexCredit cost: 0.706. Use 'web.real_cost(task_id)' to get the billed FlexCredit  webapi.py:338
           cost after a simulation run.                                                                            
           starting up solver                                                                         webapi.py:342
           running solver                                                                             webapi.py:352
[17:03:29] early shutoff detected, exiting.                                                           webapi.py:366
[17:03:30] status = postprocess                                                                       webapi.py:383
[17:03:49] status = success                                                                           webapi.py:390
[17:04:09] loading SimulationData from data/simulation_data_theta.hdf5                                webapi.py:568

To confirm that the three simulation produces close results we first consider the total scattered flux from the sphere. The slight difference in obtained results is attributed to the fact that currently Tidy3D does not support subpixel averaging for fully anisotropic materials.

[15]:
# retrive flux values from flux monitor
flux_diag = sim_data_diag["flux"].flux
flux_phi = sim_data_phi["flux"].flux
# scale power by injection angle
flux_theta = sim_data_theta["flux"].flux * np.cos(theta)

# visualize
flux_diag.plot()
flux_phi.plot()
flux_theta.plot()
plt.legend(["diag", "phi", "theta"])
plt.show()
../_images/notebooks_FullyAnisotropic_31_0.png

Now we look at the angular field intensity distributions around the sphere at several distances from the sphere’s center.

[16]:
# field sample locations
sample_radius = [0.1, 0.2, 0.4, 0.5]
sample_angle = np.linspace(0, 2 * np.pi, 100)

When comparing the reference (diagonally anisotropic) simulation with the simulation obtained by rotation around \(x\) axis we use the fact that field distributions must coincide in the \(yz\) plane when appropriately rotated.

[17]:
# get field intensity distribution at central frequency
int_diag = sim_data_diag.get_intensity("field_yz").isel(x=0, f=find0)
int_full = sim_data_phi.get_intensity("field_yz").isel(x=0, f=find0)

fig, ax = plt.subplots(len(sample_radius), 1, constrained_layout=True, figsize=(8, 8))

for i, r in enumerate(sample_radius):
    # sample points in reference simulation
    y_diag = xr.DataArray(r * np.cos(sample_angle), dims="u")
    z_diag = xr.DataArray(r * np.sin(sample_angle), dims="u")

    # rotated sample points in fully anisotropic simulation
    y_full = xr.DataArray(r * np.cos(sample_angle + phi), dims="u")
    z_full = xr.DataArray(r * np.sin(sample_angle + phi), dims="u")

    # interpolate at sample points
    int_sampled_diag = int_diag.interp(y=y_diag, z=z_diag)
    int_sampled_full = int_full.interp(y=y_full, z=z_full)

    # plot comparisons
    ax[i].plot(sample_angle, int_sampled_diag.data)
    ax[i].plot(sample_angle, int_sampled_full.data)
    ax[i].set_xlabel("Angle [rad]")
    ax[i].legend(["Reference", "Fully Anisotropic"])
    ax[i].set_title(f"Electric field intensity at r = {r}")

plt.show()
../_images/notebooks_FullyAnisotropic_35_0.png

When comparing the reference (diagonally anisotropic) simulation with the simulation obtained by rotation around \(z\) axis we use the fact that field distributions must coincide in the \(xy\) plane when appropriately rotated.

[18]:
# get field intensity distribution at central frequency
int_diag = sim_data_diag.get_intensity("field_xy").isel(z=0, f=find0)
# scale power by the injection angle
int_full = sim_data_theta.get_intensity("field_xy").isel(z=0, f=find0) * np.cos(theta)

fig, ax = plt.subplots(len(sample_radius), 1, constrained_layout=True, figsize=(8, 8))

for i, r in enumerate(sample_radius):
    # sample points in reference simulation
    x_diag = xr.DataArray(r * np.cos(sample_angle), dims="u")
    y_diag = xr.DataArray(r * np.sin(sample_angle), dims="u")

    # rotated sample points in fully anisotropic simulation
    x_full = xr.DataArray(r * np.cos(sample_angle + theta), dims="u")
    y_full = xr.DataArray(r * np.sin(sample_angle + theta), dims="u")

    # interpolate at sample points
    int_sampled_diag = int_diag.interp(x=x_diag, y=y_diag)
    int_sampled_full = int_full.interp(x=x_full, y=y_full)

    # plot comparisons
    ax[i].plot(sample_angle, int_sampled_diag.data)
    ax[i].plot(sample_angle, int_sampled_full.data)
    ax[i].set_xlabel("Angle [rad]")
    ax[i].legend(["Reference", "Fully Anisotropic"])
    ax[i].set_title(f"Electric field intensity at r = {r}")

plt.show()
../_images/notebooks_FullyAnisotropic_37_0.png

Again, we see that equivalent simulations produce very close results, and the small differences are attributed to the lack of subpixel averaging for fully anisotropic materials.