emopt.fdtd¶

Implements a simple 3D CW-FDTD (continuous wave finite difference frequenct domain) solver for simulating Maxwell’s equations.

The FDTD method provides a simple, efficient, highly parallelizable, and scalable solution for solving Maxwell’s equations. Typical implementations of FDTD, however, can be tricky to use with highly dispersive materials and are tough to use in conjunction with the adjoint method for inverse electromagnetic design.

In order to overcome these issues, this implementation uses a ramped-CW source which allows us to solve for the frequency domain fields at the exact frequency we care about without doing any form of interpolation. While this eliminates one of FDTD’s advantages (obtaining broadband info with a single simulation), it is easier to implement and plays nice with the adjoint method which is easily derivable for frequency-domain solvers. Furthermore, contrary to explicit frequency domain solvers, CW-FDTD is highly parallelizable which enables it to run quite a bit faster. Furthermore, compared to frequency domain solvers like FDFD, CW-FDTD consumes considerably less RAM, making it useful for optimizing much larger devices at much higher resolutions.

class emopt.fdtd.FDTD(X, Y, Z, dx, dy, dz, wavelength, rtol=1e-06, nconv=None, min_rindex=1.0, complex_eps=False)¶

Bases: emopt.simulation.MaxwellSolver

A 3D continuous-wave finite difference time domain (CW-FDTD) solver.

This class implements a continuous-wave finite difference time domain solver which solves for the freqeuncy-domain fields using the finite difference time domain method. Unlike the typical pulsed FDTD method, this implementation uses a ramped CW source which ensures that the simulation error exhibits convergent behavior (which is useful for optimization). Furthermore, this CW eliminates the need to perform discrete fourier transforms and allows us to compute the fields at the exact frequency/wavelength we desire (without relying on any interpolation which complicates optimizations).

Compared to fdfd.FDFD_3D, FDTD will generally scale to larger/higher resolution problems much better. For such problems, the memory consumption can be an order of magnitude or more lower and it parallelizes significantly better. For very small problems, however, fdfd.FDFD_3D may be faster.

Notes

1. Complex materials are not yet supported, however they will be in the future!

2. Because of how power is computed, you should not do any calculations within 1 grid cells of the PMLs

3. Currently, the adjoint simulation can have some difficulty converging when used for a power-normalized figure of merit. To get around this, you can either increase the relative tolerance, or better yet set a maximum number of time steps using Nmax.

4. Power is currently computed using a transmission box which is placed 1 grid cell within the PML regions and does not account for material absorption.

Parameters:

X (float) – The width of the simulation in the x direction
Y (float) – The width of the simulation in the y direction
Z (float) – The width of the simulation in the z direction
dx (float) – The grid spacing in the x direction
dy (float) – The grid spacing in the y direction
dz (float) – The grid spacing in the z direction
wavelength (float) – The wavelength of the source (and fields)
rtol (float (optional)) – The relative tolerance used to terminate the simulation. (default = 1e-6)
nconv (int (optional)) – The number of field points to check for convergence. (default = 100)
min_rindex (float (optional)) – The minimum refractive index of the simulation. Setting this parameter correctly can speed up your simulations by quite a bit, however you must be careful that its value does not exceed the minimum refractive index in your simulation. (default = 1.0)
complex_eps (boolean (optional)) – Tells the solver if the permittivity is complex-valued. Setting this to False can speed up the solver when run on fewer cores (default = False)

X¶

The width of the simulation in the x direction

Type:	float

Y¶

The width of the simulation in the y direction

Type:	float

Z¶

The width of the simulation in the z direction

Type:	float

dx¶

The grid spacing in the x direction

Type:	float

dy¶

The grid spacing in the y direction

Type:	float

dz¶

The grid spacing in the z direction

Type:	float

Nx¶

The number of grid cells along the x direction

Type:	int

Ny¶

The number of grid cells along the y direction

Type:	int

Nz¶

The number of grid cells along the z direction

Type:	int

courant_num¶

The courant number (fraction of maximum time step needed for stability). This must be <= 1. By default it is 0.95. Don’t change this unless you simulations appear to diverge for some reason.

Type:	float

wavelength¶

The wavelength of the source (and fields)

Type:	float

eps¶

The permittivity distribution.

Type:	emopt.grid.Material3D

mu¶

The permeability distribution.

Type:	emopt.grid.Material3D

rtol¶

The relative tolerance used to terminate the simulation. (default = 1e-6)

Type:	float (optional)

src_min_value¶

The starting value of the ramped source. The source is ramped using a smooth envelope function which never truely goes to zero. This attribute lets us set the “zero” value which will effect the convergence time and minimum error. (default = 1e-5)

Type:	float

src_ramp_time¶

The number of time steps over which the source is ramped. This will affect the convergence time and final error. Note: faster ramps lead to higher error (in the current implementation) and do not gurantee that the simulation will converge faster. (default = 15*Nlambda)

Type:	int

Nmax¶

Max number of time steps

Type:	int

Nlambda¶

Number of spatial steps per wavelength (in minimum index material)

Type:	float

Ncycle¶

Number of time steps per period of oscillation of the fields.

Type:	float

bc¶

The the boundary conditions (PEC, field symmetries, etc)

Type:	str

w_pml¶

The list of pml widths (in real coordinates) in the format [xmin, xmax, ymin, ymax, zmin, zmax]. Use this to change the PML widths.

Type:	list of floats

w_pml_xmin¶

The number of grid cells making up with PML at the minimum x boundary.

Type:	int

w_pml_xmax¶

The number of grid cells making up with PML at the maximum x boundary.

Type:	int

w_pml_ymin¶

The number of grid cells making up with PML at the minimum y boundary.

Type:	int

w_pml_ymax¶

The number of grid cells making up with PML at the maximum y boundary.

Type:	int

w_pml_zmin¶

The number of grid cells making up with PML at the minimum z boundary.

Type:	int

w_pml_zmax¶

The number of grid cells making up with PML at the maximum z boundary.

Type:	int

X_real¶

The width of the simulation excluding PMLs.

Type:	float

Y_real¶

The height of the simulation excluding PMLs.

Type:	float

Ncycle

Nlambda

Nx

Ny

Nz

X

X_real

Y

Y_real

Z

Z_real¶

add_adjoint_sources(src, domain)¶

Add an adjoint source.

Add an adjoint source using a set of 6 current source arrays.

Parameters:	src (list or tuple of numpy.ndarray) – The list/tuple of source arrays in the following order: (Jx, Jy, Jz, Mx, My, Mz) domain (misc.DomainCoordinates) – The domain which specifies the location of the source arrays

bc

build()¶

Assemble the strcutre.

This involves computing the permittiviy and permeability distribution for the simulation.

calc_ydAx(Adiag0)¶

Calculate the product y^T*dA*x.

Parameters:	Adiag0 (tuple of 6 numpy.ndarray) – The ‘initial’ diag[A] obtained from self.get_A_diag()
Returns:	The product y^TdAx
Return type:	complex

clear_adjoint_sources()¶: Clear the adjoint simulation sources.

clear_sources()¶: Clear simulation sources.

courant_num

dx

dy

dz

eps

get_A_diag()¶

Get a representation of the diagonal of the discretized Maxwell’s equations assuming they are assembled in a matrix in the frequency domain.

For the purposes of this implementation, this is just a copy of the permittivity and permeability distribution. In reality, there should be a factor of 1j and -1j for the permittivities and permeabilities, respectively, however we handle these prefactors elsewhere.

Returns:	A copy of the set of permittivity and permeability distributions used internally.
Return type:	tuple of numpy.ndarrays

get_adjoint_field(component, domain=None)¶

Get the adjoint field.

Parameters:	component (str) – The adjoint field component to retrieve. domain (misc.DomainCoordinates (optional)) – The domain in which the field is retrieved. If None, retrieve the field in the whole 3D domain (which you probably should avoid doing for larger problems). (default = None)
Returns:	An array containing the field.
Return type:	numpy.ndarray

get_field(component, domain=None)¶

Get the (raw, uninterpolated) field.

In most cases, you should use :def:`get_field_interp` instead.

Parameters:	component (str) – The field component to retrieve domain (misc.DomainCoordinates (optional)) – The domain in which the field is retrieved. If None, retrieve the field in the whole 3D domain (which you probably should avoid doing for larger problems). (default = None)
Returns:	An array containing the field.
Return type:	numpy.ndarray

get_field_interp(component, domain=None, squeeze=False)¶

Get the desired field component.

Internally, fields are solved on a staggered grid. In most cases, it is desirable to know all of the field components at the same sets of positions. This requires that we interpolate the fields onto a single grid. In emopt, we interpolate all field components onto the Ez grid.

Parameters:	component (str) – The desired field component. domain (misc.DomainCoordinates (optional)) – The domain from which the field is retrieved. (default = None)
Returns:	The interpolated field
Return type:	numpy.ndarray

get_source_power()¶

Get source power.

The source power is the total electromagnetic power radiated by the electric and magnetic current sources.

Returns:	The source power.
Return type:	float

mu

rtol

set_adjoint_sources(src)¶

Set the adjoint sources.

This function works a bit differently than set_sources in order to maintain compatibility with the adjoint_method class. It takes a single argument which is a list of lists containing multiple sets of source arrays and corresponding domains.

Todo

Clean up EMopt interfaces so that this doesn’t feel so hacky.

src : list: List containing two lists: 1 with sets of source arrays and on with DomainCoordinates which correspond to those source array sets.

set_materials(eps, mu)¶

set_sources(src, domain, mindex=0)¶

Set a simulation source.

Simulation sources can be set either using a set of 6 arrays (Jx, Jy, Jz, Mx, My, Mz) or a modes.ModeFullVector object. In either case, a domain must be provided which tells the simulation where to put those sources.

This function operates in an additive manner meaning that multiple calls will add multiple sources to the simulation. To replace the existing sources, simply call :def:`clear_sources` first.

Parameters:	src (tuple or modes.ModeFullVector) – The source arrays or mode object containing source data domain (misc.DomainCoordinates) – The domain which specifies where the source is located mindex (int (optional)) – The mode source index. This is only relevant if using a ModeFullVector object to set the sources. (default = 0)

solve_adjoint()¶

Run an adjoint simulation.

In FDTD, this is a solution to Maxwell’s equations but using a different set of sources than the forward simulation.

solve_forward()¶

Run a forward simulation.

A forward simulation is just a solution to Maxwell’s equations.

src_min_value

src_ramp_time

test_src_func(**kwargs)¶

update(bbox=None)¶

Update the permittivity and permeability distribution.

A small portion of the whole simulation region can be updated by supplying a bounding box which encompasses the desired region to be updated.

Notes

The bounding box accepted by this class is currently different from the FDFD solvers in that it is specified using real-space coordinates not indices. In the future, the other solvers will implement this same functionality.

Todo

Implement 3 point frequency domain calculation. This should be more accurate and produce more consistent results.

Parameters:	bbox (list of floats (optional)) – The region of the simulation which should be updated. If None, then the whole region is updated. Format: [xmin, xmax, ymin, ymax, zmin, zmax]

update_saved_fields()¶

Update the fields contained in the regions specified by self.field_domains.

This function is called internally by solve_forward and should not need to be called otherwise.

w_pml

w_pml_xmax

w_pml_xmin

w_pml_ymax

w_pml_ymin

w_pml_zmax

w_pml_zmin

wavelength

class emopt.fdtd.GhostComm(k0, j0, i0, K, J, I, Nx, Ny, Nz)¶

Bases: future.types.newobject.newobject

update_local_vector(vl)¶

Update the ghost values of a local vector.

Parameters:	vl (np.ndarray) – The local vector

class emopt.fdtd.SourceArray(Jx, Jy, Jz, Mx, My, Mz, i0, j0, k0, I, J, K)¶

Bases: future.types.newobject.newobject

A container for source arrays and its associated parameters.

Parameters:

Jx (numpy.ndarray) – The x-component of the electric current density
Jy (numpy.ndarray) – The x-component of the electric current density
Jz (numpy.ndarray) – The x-component of the electric current density
Mx (numpy.ndarray) – The x-component of the magnetic current density
My (numpy.ndarray) – The x-component of the magnetic current density
Mz (numpy.ndarray) – The x-component of the magnetic current density
i0 (int) – The LOCAL starting z index of the source arrays
j0 (int) – The LOCAL starting y index of the source arrays
k0 (int) – The LOCAL starting k index of the source arrays
I (int) – The z width of the source array
J (int) – The y width of the source array
K (int) – The z width of the source array