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 define...
Sign up now to simulate .ipynb

How to define fully anisotropic materials in Tidy3D FDTD¶

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

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

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'.                                                          
           View task using web UI at                                                                  webapi.py:188
           'https://tidy3d.simulation.cloud/workbench?taskId=fdve-78d231c2-61bf-4a03-b57b-7f61b6e14d1              
           8v1'.                                                                                                   
Output()


[16:59:02] status = queued                                                                            webapi.py:321
Output()
[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
Output()
[16:59:59] early shutoff detected, exiting.                                                           webapi.py:366


           status = postprocess                                                                       webapi.py:383
Output()
[17:00:25] status = success                                                                           webapi.py:390

Output()


[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'.                                                          
           View task using web UI at                                                                  webapi.py:188
           'https://tidy3d.simulation.cloud/workbench?taskId=fdve-23947891-340a-4ea7-8540-86162eeb55d              
           9v1'.                                                                                                   
Output()


[17:00:45] status = queued                                                                            webapi.py:321
Output()
[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
Output()
[17:01:43] early shutoff detected, exiting.                                                           webapi.py:366


           status = postprocess                                                                       webapi.py:383
Output()
[17:02:09] status = success                                                                           webapi.py:390

Output()


[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'.                                                          
           View task using web UI at                                                                  webapi.py:188
           'https://tidy3d.simulation.cloud/workbench?taskId=fdve-169281b2-d8df-4764-a7eb-1e815dd6b28              
           cv1'.                                                                                                   
Output()


[17:02:29] status = queued                                                                            webapi.py:321
Output()
[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
Output()
[17:03:29] early shutoff detected, exiting.                                                           webapi.py:366


[17:03:30] status = postprocess                                                                       webapi.py:383
Output()
[17:03:49] status = success                                                                           webapi.py:390

Output()


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

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

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

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.


Aerodynamics & Electromagnetic Simulation | Flexcompute
Thanks for subscribing
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
Resources
  • FDTD 101
  • Aircraft CFD 101
  • CFD Essentials
  • Newsroom
  • Blog
  • Python FDTD
  • FDTD
  • FDTD Acceleration
About
  • Vision
  • Story
  • Team
  • Careers
  • Contact Us
  • Terms & Conditions
  • Privacy Policy
  • Trust Guide
  • Copyright © 2023 Flexcompute, inc.

This website uses cookies to ensure you get the best experience. Learn more about our cookie policy.

Accept

Subscribe

Thanks for subscribing

Download Tutorial Presentation Slides

Enter your email address below to receive the presentation slides. In the future, we’ll share you very few emails when we have new tutorial release, development updates, valuable toolkits and technical guidance . You can unsubscribe at any time by clicking the link at the bottom of every email. We’ll never share your information.

By clicking below, you are agreeing to Flexcompute's Privacy Policy.