Running a simulation is as simple as using the grid.run method.

Just like for the lengths in the grid, the total_time of the simulation can be specified as an integer (number of time_steps) or as a float (in seconds).

Let's visualize the grid. This can be done with the grid.visualize method:

This method will by default visualize all objects in the grid, as well as the field intensity at the current time_step at a certain x, y OR z-plane. By setting show=False, one can disable the immediate visualization of the matplotlib image.

An as quick as possible explanation of the FDTD discretization of the Maxwell equations.

An electromagnetic FDTD solver solves the time-dependent Maxwell Equations

These two equations are called Ampere's Law and Faraday's Law respectively.

In these equations, e and u are the relative permittivity and permeability tensors respectively. e0 and u0 are the vacuum permittivity and permeability and their square root can be absorbed into E and H respectively, such that E := e0*E and H := u0*H.

Doing this, the Maxwell equations can be written as update equations:

The electric and magnetic field can then be discretized on a grid with interlaced Yee-coordinates, which in 3D looks like this:

According to the Yee discretization algorithm, there are inherently two types of fields on the grid: E-type fields on integer grid locations and H-type fields on half-integer grid locations.

The beauty of these interlaced coordinates is that they enable a very natural way of writing the curl of the electric and magnetic fields: the curl of an H-type field will be an E-type field and vice versa.

This way, the curl of E can be written as

this can be written efficiently with array slices (note that the factor (1/du) was left out):

The curl for H can be obtained in a similar way (note again that the factor (1/du) was left out):

The update equations can now be rewritten as

The number (c*dt/du) is a dimensionless parameter called the Courant number sc. For stability reasons, the Courant number should always be smaller than 1/D, with D the dimension of the simulation. This can be intuitively be understood as the condition that information should always travel slower than the speed of light through the grid. In the FDTD method described here, information can only travel to the neighboring grid cells (through application of the curl). It would therefore take D time steps to travel over the diagonal of a D-dimensional cube (square in 2D, cube in 3D), the Courant condition follows then automatically from the fact that the length of this diagonal is 1/D.

This yields the final update equations for the FDTD algorithm:

This is also how it is implemented:

Ampere's Law can be updated to incorporate a current density:

Making again the usual substitutions sc := c*dt/du, E := e0*E and H := u0*H, the update equations can be modified to include the current density:

Making one final substitution Es := -dt*inv(e)*J/e0 allows us to write this in a very clean way:

Where we defined Es as the electric field source term.

It is often useful to also define a magnetic field source term Hs, which would be derived from the magnetic current density if it were to exist. In the same way, Faraday's update equation can be rewritten as

When a material has a electric conductivity s, a conduction-current will ensure that the medium is lossy. Ampere's law with a conduction current becomes

Making the usual substitutions, this becomes:

This update equation depends on the electric field on a half-integer time step (a magnetic field time step). We need to substitute E(t+dt/2)=(E(t)+E(t+dt))/2 to interpolate the electric field to the correct time step.

Which, yield the new update equations:

Note that the more complicated the permittivity tensor e is, the more time consuming this algorithm will be. It is therefore sometimes a nice hack to transfer the absorption to the magnetic domain by introducing a (nonphysical) magnetic conductivity, because the permeability tensor u is usually just equal to one:

The electromagnetic energy density can be given by

making the above substitutions, this becomes in simulation units:

The Poynting vector is given by

Which in simulation units becomes

The energy introduced by a source Es can be derived from tracking the change in energy density

This could also be derived from Poyntings energy conservation law:

where the first term just describes the redistribution of energy in a volume and the second term describes the energy introduced by a current density.

Note: although it is unphysical, one could also have introduced a magnetic source. This source would have introduced the following energy:

Since the u-tensor is usually just equal to one, using a magnetic source term is often more efficient.

Similarly, one can also keep track of the absorbed energy due to an electric conductivity in the following way:

or if we want to keep track of the absorbed energy by magnetic a magnetic conductivity:

The electric term and magnetic term in the energy density are usually of the same size. Therefore, the same amount of energy will be absorbed by introducing a magnetic conductivity sm as by introducing a electric conductivity s if:

Assuming we want periodic boundary conditions along the X-direction, then we have to make sure that the fields at Xlow and Xhigh are the same. This has to be enforced after performing the update equations:

Note that the electric field E is dependent on curl_H, which means that the first indices of E will not be updated through the update equations. It's those indices that need to be set through the periodic boundary condition. Concretely: E[0] needs to be set to equal E[-1]. For the magnetic field, the inverse is true: H is dependent on curl_E, which means that its last indices will not be set. This has to be done by the boundary condition: H[-1] needs to be set equal to H[0]:

a Perfectly Matched Layer (PML) is the state of the art for introducing absorbing boundary conditions in an FDTD grid. A PML is an impedance-matched absorbing area in the grid. It turns out that for a impedance-matching condition to hold, the PML can only be absorbing in a single direction. This is what makes a PML in fact a nonphysical material.

Consider Ampere's law for the Ez component, where we use the following substitutions: E := e0*E, H := u0*H and s := inv(e)*s/e0 are already introduced:

This becomes in the frequency domain:

We can split this equation in a x-propagating wave and a y-propagating wave:

We can define the S-operators as follows

In general, we prefer to add a stability factor au and a scaling factor ku to Su:

Summing the two equations for Ez back together after dividing by the respective S-operator gives

Converting this back to the time domain gives

where sx denotes the inverse Fourier transform of (1/Sx) and [*] denotes a convolution. The expression for su can be proven [after some derivation] to look as follows:

where d(t) denotes the Dirac delta function and C(t) an exponentially decaying function given by:

Plugging this in gives:

This can be written as an update equation:

Where we defined Feu as

and Pseuv as the convolution updating the component Eu by taking the derivative of Hw in the v direction:

This can be rewritten [after some derivation] as an update equation in itself:

Where the constants bu and cu are derived to be:

The final PML algorithm for the electric field now becomes:

1. Update Fe=[Fex, Fey, Fez] by using the update equation for the Ps-components.
2. Update the electric fields the normal way
3. Add Fe to the electric fields.

or as python code:

The same has to be applied for the magnetic field.

These update equations for the PML were based on Schneider, Chap. 11.

As a bare FDTD library, this is dimensionally agnostic for any unit system you may choose. No conversion factors are applied within the library API; this is left to the user. (The code used to calculate the Courant limit may be a sticking point depending on the time scale involved).

However, as noted above (H := u0*H), it is generally good numerical practice to scale all values to get the maximum precision from floating-point types.

In particular, a scaling scheme detailed in "Novel architectures for brain-inspired photonic computers", Chapters 4.1.2 and 4.1.6, is highly recommended.

A set of conversion functions to and from reduced units are available for users in conversions.py.

You can run a linter in the root using pylint fdtd.

(c) Floris laporte - MIT License