FDTD Plus

This is a Wiki for a finite-difference time-domain (FDTD) simulation software package written in C++ and developed by Sunil Sandhu of Fan Group at Stanford University.

The purpose of this Wiki is to

• document the software package.
• allow group members to post feedback on the software.

Overview

The name of the FDTD software package is FDTD Plus. A parallel version of FDTD Plus executable is available to Fan group members at several servers,

• XSEDE Stampede /home1/00056/tg457106/share/bin/fdtd_plus_mpi
• Hera /usr/local/bin/fdtd_plus_mpi
• XSEDE COMET /users/ux455370/share/fdtd_plus_mpi
• corn.stanford.edu /afs/ir.stanford.edu/users/c/e/centaur/share/fdtd_plus_mpi
• The executables above do not have the CellParam::PrecisionType = float option.
• To run FDTD Plus:

   fdtd_plus_mpi <directory location> <input file name> <x-proc> <y-proc> <z-proc>

• <directory location> is the directory where the input file is located and is the location where all FDTD Plus outputs are performed.
• <x-proc>, <y-proc> and <z-proc> are the number of processors in the 3 directions.
  • In the next release, the allocation of processors in the 3 directions will be done internally.
• During the run, errors, warnings, status of the run etcetera are output to a file called output.dat.
• Email centaur(at)stanford.edu if you can't access any of the above sites or for any other questions.

Convenience script (not implemented in Hera)

There is a convenience script available for submitting FDTD Plus jobs into any queue in the above platforms using a minimum set of arguments. The usage of the script at the command line is

             qsub_fdtd_plus <options>

• To obtain information on the available options of qsub_fdtd_plus type

             qsub_fdtd_plus --help

Example usage

• An example use is

             qsub_fdtd_plus -o=$SCRATCH/run01 -f=$HOME/input/input01.txt -t=00:30:00 -x=2

which will submit a FDTD Plus run with the following settings

• user input file path $HOME/input/input01.txt (this file will be copied to the working directory indicated in the -o option),
• working directory (i.e. where all outputs are done) at $SCRATCH/run01 (note: run01 directory is created if it does not exist),
• for a wall time of 30min,
• using 2 processors in the x direction (default value for -x,-y,-z directions is 1 processor).

Access to script

For Shanhuiserve, you can use it right away at the command line. As for other platforms, the following steps have to be performed before qsub_fdtd_plus can be used at the command line.

• For SDSC, add the following at the end of the .soft file (.soft is located in your $HOME directory),

         PATH += /users/ux455370/share

User input file

The FDTD Plus simulation is set up based on instructions written in a user input file.

Input file syntax

The input file consists of instructions enclosed in parentheses (round brackets).

            syntax of an instruction:
                          (instruction)

            error! no closing parenthesis:
                          (instruction_1)
                          (instructio_2     <-instruction_2 missing closing parenthesis
                          (instruction_3)

Certain instructions introduced in later sections have nested instructions.

           nested instruction example:
                           (instruction_1 (instruction_2) )

The main rule in the user input file is never have two parentheses together.

     this will result in an error: 
                           (instruction_1 (instruction_2))  <-closed parentheses together
            
     this fixes the above error: 
                           (instruction_1 (instruction_2) )

An instruction can be disabled by adding a slash "\" to the open parenthesis.

           disable instruction: 
                            (\instruction)

           only instruction 2 is disabled: 
                            (instruction_1 (\instruction_2) )

Having "(\" together with another parenthesis (e.g. "(\(" ) will result in an error too.

A block of instructions can be disabled by enclosing them in "(\" and ")".

          instructions 3-5 are disabled:
                            (instruction_1)   
                            (instruction_2) 
                       (\            
                            (instruction_3)   
                            (instruction_4)   
                            (instruction_5) 
                        )
                            (instruction_6)   
                            (instruction_7)

Comments can be added anywhere outside instructions. If parentheses (i.e. round brackets) need to be used for comments, the open parenthesis must have a slash attached to it, "(\".

• Having a single "(\" without an enclosing ")" will result in an error too.

         this will result in an error:
                       This comment is fine.
                       #This is another example of a comment that is fine too
                       (\This is a comment that is fine.. I guess you get the idea!)
                       (\This comment has an error..  <-no closing parenthesis 
                                   (instruction)

Available instructions

There are five groups of instructions, namely Options, Materials, Sources, Objects, and Outputs .

• All instructions except for Objects can be placed in any order. Source Objects and Material Objects are processed separately according to the order they are written.
  • The order of Source Objects have no relation with the order of Material Objects (see comment in Sources).

The following sections describe the different possible user input instructions available.

Example file

Click here for an example input file.

Options

This user input sets up the dimensions, boundaries and cell type for the simulation.

        syntax:
                (Option OptionType <parameters> )

• In the definition of the different options below, any option that is not marked "(optional)" is a required option.

CellParam

Defines the type of cell used in the FDTD computation including the field type.

         (Option CellParam (CellType      value) 
                           (FieldType     value)
                           (PrecisionType value) 
         )

CellType

The type of cells used to represent the materials in the computation region,
$value\in \{$linear$\}$,

• linear is a field independent material.
  • A linear cell only communicates E- and H-fields to neighboring cells.

FieldType

The mathematical field type of the cell physical fields that are communicated to other cells,
$value\in \{$complex, real$\}$.

• Most simulations can be performed using FieldType = real. If a material type or output type requires FieldType = complex, this will be highlighted in its description below.

PrecisionType (optional)

The floating point precision of fields, coefficients etc in the cell,
$value\in \{$double, float$\}$.

• Default value = double.
• To include the PrecisionType = float during compilation of FDTD Plus, use the FDTD_PLUS_PRECISION_TYPE_FLOAT flag.
  • The FDTD Plus executable available in the servers above do not have the PrecisionType = float option.

DimensionParam

Defines the spatial and time dimensions of the compute region.

     (Option DimensionParam (LengthScale   value)         
                            (Center        x_value y_value z_value) 
                            (Size          x_value y_value z_value) 
                            (Resolution    x_value y_value z_value)
                            (TimeParam     (Start value) (End value) ) 
                            (CourantNumber value)  
     )

LengthScale (optional)

The length scale of the compute region,
$value\in \mathbb {R} _{>0}$.

• Unit = meters.
• LengthScale does not need to be defined for lossless, non-dispersive, linear materials.
  • Since the Maxwell equations are scale invariant in the case of lossless, non-dispersive, linear materials, a length scale (i.e. LengthScale) can be used to represent a unit of distance in these cases.
  • In the description of material types below, materials requiring the definition of LengthScale will be highlighted.

Center

The center coordinates of the compute region,
$x\_value,y\_value,z\_value\in \mathbb {R}$

• Unit = LengthScale.

Size

The size of the compute region,
$x\_value,y\_value,z\_value\in \mathbb {R} _{>0}$

• Unit = LengthScale.

Resolution

The spatial increment in the compute region,
$x\_value,y\_value,z\_value\in \mathbb {R} _{>0}$

• Unit = LengthScale.

TimeParam

The time parameters for the FDTD run.

Start (optional)

The start time of computation
$value\in \mathbb {R}$

• Unit = LengthScale / c, where c = speed of light in vacuum.
• Default value = 0.

End

The end time of computation
$value\in \mathbb {R} _{\geq {Start}}$

• Unit = LengthScale / c.

CourantNumber (optional)

The courant number,
$value\in \mathbb {R} _{>0}$.

• Default value = $0.5\cdot {\sqrt {3}}\approx {0.866}$.
• Using CourantNumber, the FDTD time increment is calculated as

                  $\Delta {t}={\frac {S}{\sqrt {{\frac {1}{\Delta {x}^{2}}}+{\frac {1}{\Delta {y}^{2}}}+{\frac {1}{\Delta {z}^{2}}}}}}$

• where,
  • $\Delta {t}$ has unit = LengthScale / c,
  • $\Delta {x}$,$\Delta {y}$,$\Delta {z}$ = Resolution,
  • S = CourantNumber.
• In order to ensure numerical stability, CourantNumber should be slightly less than $n_{min}$,
  • $n_{min}$ = minimum refractive index in the compute region.

BoundaryParam

Defines the boundaries of the compute region.

     (Option BoundaryParam (Type        x_value y_value z_value) 
                           (PmlParam   (Size (X value1 value2) (Y value1 value2) (Z value1 value2) )
                                       (SigmaFactor (X value1 value2) (Y value1 value2) (Z value1 value2) ) 
                                       (KappaFactor (X value1 value2) (Y value1 value2) (Z value1 value2) )    
                                       (Type value)
                                       (DecayConstant value)
                           )    
                           (BlochK     (X value) (Y value) (Z value) )
                           (Symmetry   (X value) (Y value) (Z value) )
     )

Type

The type of boundaries in the compute region,
$x\_value,y\_value,z\_value\in \{$pml, bloch, none$\}$,

• pml is a Perfectly Matched Layer (PML) boundary,
• bloch is a Bloch boundary,
• none results in the boundary that includes the origin of the compute region to be of type PEC (E-fields parallel to boundary set to zero), and the opposite boundary of the compute region to be of type PMC (H-fields parallel to boundary set to zero).

PmlParam (optional)

The parameters for PML boundaries.

• PmlParam is only required to be defined if any BoundaryParam::Type == pml.

Size

The size of PML boundaries.

• Size is included in DimensionParam::Size.
• Size is only required to be defined for directions with BoundaryParam::Type == pml.

value1

The PMLboundary size for boundaries that include the origin of the compute region,
$value1\in \mathbb {R} _{\geq {0}}$

• Unit = DimensionParam::LengthScale.

value2 (optional)

The PMLboundary size for boundaries that do not include the origin of the compute region,
$value2\in \mathbb {R} _{\geq {0}}$

• Unit = DimensionParam::LengthScale.
• Default value = $value1$.

SigmaFactor (optional)

This controls the imaginary part of the $s$parameter of UPML (see Taflove). Controls the normal incidence reflection R(0) at the PMLboundaries. Generally, $R(0)=Exp\left[-\gamma \right]$ where $\gamma =SigmaFactor\times {1.6}\times {\frac {PmlParams::Size}{DimensionParams::Resolution}}$

• This parameter adds flexibility to designing your PML boundaries.
• By default SigmaFactor = 0.8. This results in $\gamma =12.8$ for 10 PML cells. But for 20 PML cells, $\gamma =25.6$ (which can cause numerical issues). Be careful or you may get undesirable reflections or numerical issues. Read Taflove for more information.
• You should read Taflove's Computational Electrodynamics (within chapter on PML, section on Theoretical Performance of the PML) to get a better idea of this parameter.

value1

The SigmaFactor for PMLboundaries that include the origin of the compute region,
$value1\in \mathbb {R} _{>0}$

• Unit = Dimensionless.
• Default value = $0.8$.

value2 (optional)

The SigmaFactor for PMLboundaries that do not include the origin of the compute region,
$value2\in \mathbb {R} _{>0}$

• Unit = Dimensionless.
• Default value = $value1$.

KappaFactor (optional)

This controls the real part of the $s$ parameter of UPML (see Taflove).

• KappaFactor is $\kappa _{max}$ in Taflove.

• This parameter adds flexibility to designing your PML boundaries.
• By default KappaFactor = 1.
• Only use this parameter if you REALLY know what you are doing.

value1

The KappaFactor for PMLboundaries that include the origin of the compute region,
$value1\in \mathbb {R} _{>0}$

• Unit = Dimensionless.
• Default value = $1$.

value2 (optional)

The KappaFactor for PMLboundaries that do not include the origin of the compute region,
$value2\in \mathbb {R} _{>0}$

• Unit = Dimensionless.
• Default value = $value1$.

Type

The type of PMLboundary,
$value\in \{$upml$\}$,

• upml is a Uniaxial PML. For Uniaxial PML,
  • PmlParam::DecayConstant corresponds to the reflection error R(0).
  • The order of pml variation (m) in a particular boundary is chosen to be m = 4 when number of pml cells < 50, and m = 2 for other cases.

DecayConstant (optional)

The PMLdecay constant,
$value\in \mathbb {R} _{>0}$

• Default value = $e^{-16}$.

BlochK (optional)

The Bloch k-vector for Bloch boundaries,
$value\in \mathbb {R}$

• Unit = $2\pi$ / DimensionParam::Size.
• BlochK is only required to be defined for directions with BoundaryParam::Type == bloch.
• When BoundaryParam::Type == bloch, only Bloch k-vector value = 0 allows CellParam::FieldType = real.

Symmetry (optional)

For a compute region that has mirror symmetry in the x-, y- or z-direction, this option enforces the symmetry at the center plane normal to that direction (i.e. each symmetry results in a reduction by half of the corresponding DimensionParam::Size used for FDTD computation),
$value\in \{$pec, pmc$\}$,

• pec results in a Perfect Electric Conductor (PEC) center plane where the E-fields parallel to the plane are set to zero.
• pmc results in a Perfect Magnetic Conductor (PMC) center plane where the H-fields parallel to the plane are set to zero.

OutputStructure (optional)

Prints the structure within the compute region into a hdf5 file,

     (Option OutputStructure (FileName value) )

• This option is disabled by default.

FileName (optional)

The name of the output file,
$value\in {s}tring$.

• Default value = "structure".
• The suffix ".h5" is appended to FileName internally.

DefaultMaterial (optional)

Sets the default material in the compute region,

     (Option DefaultMaterial value)

$value\in {s}tring$.

• value must be a defined MaterialName (see Material below).
• Default value = "vacuum". The "vacuum" material type is generated internally.

CheckInput (optional)

Option to check the user input file. This option runs the simulation up to the point before cell construction.

     (Option CheckInput)

• This option is disabled by default.
• If OutputStructure is enabled, CheckInput will run OutputStructure.
• Read comment in Sources for a bug in CheckInput.

Materials

This user input defines a material type.

   syntax:
          (Material MaterialType (MaterialName value) 
                                 (MaterialToken value)
                                  <other parameters> 
           )

Common material parameters

MaterialName

The name of the material,
$value\in {s}tring$.

• Duplicate MaterialName values are not allowed.
• "vacuum" is a reserved MaterialName.

MaterialToken

The index of the material (used by Option::OutputStructure to print material),
$value\in \mathbb {I} _{>0}$.

• Duplicate MaterialToken values are not allowed.
• 0 is a reserved MaterialToken for "vacuum".

Linear materials

These are the material types defined for CellParam::CellType = linear.

Lossless

Defines a lossless MaterialType.

         (Material Lossless (MaterialName    value) 
                            (MaterialToken   value)
                            (RelPermittivity value)
                            (RelPermeability value) 
         )

RelPermittivity

The relative permittivity of the material,
$value\in \mathbb {R}$.

RelPermeability

The relative permeability of the material,
$value\in \mathbb {R}$.

Lossy

Defines a lossy MaterialType.

         (Material Lossy    (MaterialName    value) 
                            (MaterialToken   value)
                            (RelPermittivity value)
                            (RelPermeability value) 
                            (Conductivity    value) 
         )

• The parameters for Lossy are the same as Lossless except for the additional conductivity parameter.
• Option::DimensionParam::LengthScale must be defined for Lossy.

Conductivity

The conductivity of the material,
$value\in \mathbb {R}$.

• Unit = siemens/meter.

PRP

Defines an arbitrary dispersive MaterialType based on general complex-conjugate pole-residue pairs (PRP).

         (Material Prp      (MaterialName     value) 
                            (MaterialToken    value)
                            (RelPermeability  value)
                            (EpsilonInfinity  value) 
                            (PoleResiduePair  (Pair1 residue_value pole_value )  
                                              (Pair2 residue_value pole_value ) 
                                                              |
                                                              |
                                                              |
                                              (PairN residue_value pole_value )
                             )  
         )

• Option::DimensionParam::LengthScale must be defined for PRP.
• PRP can only be used with Option::CellParam::FieldType = real.

EpsilonInfinity

The relative permittivity at infinite frequency,
$value\in \mathbb {R}$.

PoleResiduePair

The N pole-residue pairs,
$residue\_value,\,pole\_value\,\in \,\mathbb {C}$.

• Unit = radians/second.
• If your Xth (residue, pole) pair is (a + j b, c + j d), the format for entry into the PoleResiduePair parameter is (PairX (a, b) (c, d) ), where residue_value = (a, b) and pole_value = (c, d) in the above syntax for PRP definition.

Example PRP parameters

• Click here for gold (Au) parameters.
• Click here for silver (Ag) parameters.

Sources

This user input defines a source type.

   syntax:
          (Source SourceType (SourceName value) 
                             (Amplitude  value)
                             (Duration   value) 
                             (Center     value) 
                             (Width      value)
                             (Direction  x_value y_value z_value) 
                              <other parameters> 
          )

• In FDTD Plus, sources and materials are populated using separate arrays during initialization.
  • Sources and materials are treated separately, up to the point of cell construction in the compute region, where a source from the source array is attached to a material in the material array.
    • The order in which source objects are input in the user input file have no relation with the order of material objects.
• Currently, only the Material::PRP material cannot be attached to a source.
  • If this is attempted, Option::CheckInput will not catch the error. The program will abort during cell construction. This bug will be fixed in the next release.

Common source parameters

SourceName

The name of the source,
$value\in {s}tring$.

• Duplicate SourceName values are not allowed.

Amplitude

The amplitude of the excitation,
$value\in \mathbb {R}$.

• Unit = amperes (if Option::DimensionParam::LengthScale is defined).
  • if Option::DimensionParam::LengthScale is not defined, value must be normalized to Option::DimensionParam::LengthScale by user $($i.e. if source amplitude is $I_0$ (unit = amperes), it should be entered as $value=I_{0}/Option::DimensionParam::LengthScale)$.
  • In SPP, the unit of Amplitude was not amperes. If your Amplitude value of the source in SPP was 1, enter it as 1 for FDTD Plus when Option::DimensionParam::LengthScale is not defined; if Option::DimensionParam::LengthScale is defined, enter it as the value of Option::DimensionParam::LengthScale.

Duration

The duration of the excitation,
$value\in \mathbb {R} _{\geq {0}}$.

• Unit = Option::DimensionParam::LengthScale / c.

Center

The center of the excitation,
$value\in \mathbb {R}$.

• Unit = Option::DimensionParam::LengthScale / c.

Width

The width of the excitation,
$value\in \mathbb {R}$.

• Unit = Option::DimensionParam::LengthScale / c.
• A rule of thumb is to have Width $\le$ Center / 4. This is to allow a smooth turn on of the source.

Direction

The direction (polarization) of the source,
$x\_value,y\_value,z\_value\in \mathbb {R}$.

• Unit = Option::DimensionParam::LengthScale.
• This direction vector is normalized internally during initialization.

GaussTime

Defines a Gaussian source.

         (Source GaussTime (SourceName value) 
                           (Duration   value)
                           (Center     value)
                           (Width      value) 
                           (Frequency  value)
                           (Phase      value)
                           (Amplitude  value)
                           (Direction  x_value y_value z_value)   
        )

• The Gaussian excitation is defined as,

             ${\frac {A_{0}}{\Delta _{0}\Delta _{1}}}\cdot {E}xp\left[-{\frac {(t-t_{0})^{2}}{2\cdot \Delta _{t}^{2}}}\right]\cdot \sin \left[2\pi {f}(t-t_{0})-\phi \right]\cdot \left(u\left[t-(t_{0}-{\frac {T}{2}})\right]-u\left[t-(t_{0}+{\frac {T}{2}})\right]\right)$

• t = time (Unit = Option::DimensionParam::LengthScale / c).
• $A_0$ = Amplitude.
• $\Delta _{0},\Delta _{1}$ = Option::DimensionParam::Resolution.
  • $\Delta _{0}\cdot \Delta _{1}$ is the area of the unit cell plane perpendicular to the direction of the current source.
• $\Delta _{t}$ = Width.
• $t_0$ = Center.
• $f$ = Frequency.
• $\phi$ = Phase.
• $T$ = Duration.
• $u\left[t\right]$ = unit step function.

Frequency

The center frequency of the excitation,
$value\in \mathbb {R}$.

• Unit = $2\pi \,c$ / Option::DimensionParam::LengthScale.

Phase

The phase of the excitation,
$value\in \mathbb {R}$.

• Unit = radians.

Rectangle

Defines a sinusoidal source with a Gaussian turn on/off.

         (Source Rectangle (SourceName value) 
                           (Duration   value)
                           (Center     value)
                           (Width      value) 
                           (Frequency  value)
                           (Phase      value)
                           (Amplitude  value)
                           (Direction  x_value y_value z_value)  
                           (Start        value) 
        )

• The Rectangular excitation pulse is defined as,

             ${\frac {A_{0}}{\Delta _{0}\Delta _{1}}}\left[\left(1-u(t-t_{start})\right)\cdot {E}xp\left[-{\frac {(t-t_{start})^{2}}{2\cdot \Delta _{t}^{2}}}\right]+\left(u(t-t_{start})-u(t-t_{end})\right)+u(t-t_{end})\cdot {E}xp\left[-{\frac {(t-t_{end})^{2}}{2\cdot \Delta _{t}^{2}}}\right]\right]\cdot \sin \left[2\pi {f}(t-t_{0})-\phi \right]\cdot \left(u\left[t-(t_{0}-{\frac {T}{2}})\right]-u\left[t-(t_{0}+{\frac {T}{2}})\right]\right)$

• t = time (Unit = Option::DimensionParam::LengthScale / c).
• $A_0$ = Amplitude.
• $\Delta _{0},\Delta _{1}$ = Option::DimensionParam::Resolution.
  • $\Delta _{0}\cdot \Delta _{1}$ is the area of the unit cell plane perpendicular to the direction of the current source.
• $\Delta _{t}$ = Width of the Gaussian turn on/off.
• $t_0$ = Center.
• $f$ = Frequency.
• $\phi$ = Phase.
• $T$ = Duration.
• $u\left[t\right]$ = unit step function.
• $t_{start}$ = Start of the constant envelope or peak/end of the Gaussian turn on.
• $t_{end}$ = Center + Center - $t_{start}$ = End of the constant envelope or peak/start of the Gaussian turn off.

Width

The width of the Gaussian turn on/off pulse,
$value\in \mathbb {R}$.

• Unit = Option::DimensionParam::LengthScale / c.
• A rule of thumb is to have Width $\le$ (Start - (Center - Duration / 2)) / 4. This is to allow a smooth turn on/off of the source.

Frequency

The center frequency of the excitation,
$value\in \mathbb {R}$.

• Unit = $2\pi \,c$ / Option::DimensionParam::LengthScale.

Phase

The phase of the excitation,
$value\in \mathbb {R}$.

• Unit = radians.

Start

Determines start and end of the constant envelope (or the peaks of the Gaussian turn on/off),
$value\in \mathbb {R}$.
The end of the constant envelope is calculated as Center + Center - Start.

• Unit = Option::DimensionParam::LengthScale / c.
• Default value = Center - Duration / 2 + 4 * Width.

Objects

This user input defines the objects that can be built in the compute region.

   syntax:
          (Object ObjectType (Center x_value y_value z_value) 
                             (MaterialName/SourceName value)  
                                   <other parameters> 
                             (Array (Axis0    x_value y_value z_value)
                                    (Axis1    x_value y_value z_value)
                                    (Axis2    x_value y_value z_value)
                                    (Period   value0  value1  value2) 
                                    (NoPeriod value0  value1  value2)
                             )
           )

• Either a Source or Material can be attached to an object.

Common object parameters

Center

The center coordinate of the object,
$x\_value,y\_value,z\_value\in \mathbb {R}$.

• Unit = Option::DimensionParam::LengthScale.

MaterialName/SourceName

The MaterialName (SourceName) of a defined Material (Source),
$value\in {s}tring$.

• Defining both MaterialName and SourceName for the same object will throw an initialization error. If you want to define a source and a material at the same point, you have to use separate objects.
  • Remember, in FDTD Plus a Source is attached to a Material.

Array (optional)

Creates an array of the object.

• In the description of different object types, this optional parameter will be indicated as <optional array param>.

AxisN (N = 0, 1, 2)

The Nth axis along which the object is cloned,
$x\_value,y\_value,z\_value\in \mathbb {R}$.

• Unit = Option::DimensionParam::LengthScale.
• Currently, if axes other than the scalar multiple of the x-, y-, z-unit axes are chosen, the result of object array population is undefined. This will be fixed in the next release.

Period

The spatial increment of the object center along AxisN,
$value0,value1,value2\in \mathbb {R}$.

• valueN is the spatial increment along AxisN.
• Unit = Option::DimensionParam::LengthScale.

NoPeriod

The number of periods along AxisN,
$value0,value1,value2\in \mathbb {I} _{>0}$.

• valueN is the number of periods along AxisN.

Point

Builds a point (i.e. populates one cell).

          (Object Point (Center x_value y_value z_value)
                        (MaterialName/SourceName value) 
                            <optional array param>   
          )

Plane

Builds a plane,

          (Object Plane (Center   x_value y_value z_value)
                        (MaterialName/SourceName value)
                        (Axis0    x_value y_value z_value)
                        (Axis1    x_value y_value z_value)
                        (Size     value0  value1)
                            <optional array param>
          )

AxisN (N = 0, 1)

The two axes that are parallel to any two non-parallel edges of the plane,
$x\_value,y\_value,z\_value\in \mathbb {R}$.

• Unit = Option::DimensionParam::LengthScale.
• Currently, if axes other than the scalar multiple of the x-, y-, z-unit axes are chosen, the result of the plane population is undefined. This will be fixed in the next release.

Size

The size along AxisN,
$value0,value1,value2\in \mathbb {R}$.

• valueN is the size along AxisN.
• Unit = Option::DimensionParam::LengthScale.

Block

Builds a parallelepiped,

          (Object Block (Center   x_value y_value z_value)
                        (MaterialName/SourceName value)
                        (Axis0    x_value y_value z_value)
                        (Axis1    x_value y_value z_value)
                        (Axis2    x_value y_value z_value)
                        (Size     value0  value1  value2)
                            <optional array param>
          )

AxisN (N = 0, 1, 2)

The three axes that are parallel to any three non-parallel edges of the parallelepiped,
$x\_value,y\_value,z\_value\in \mathbb {R}$.

• Unit = Option::DimensionParam::LengthScale.
• Currently, if axes other than the scalar multiple of the x-, y-, z-unit axes are chosen, the result of the block population is undefined. This will be fixed in the next release.

Size

The size along AxisN,
$value0,value1,value2\in \mathbb {R}$.

• valueN is the size along AxisN.
• Unit = Option::DimensionParam::LengthScale.

Cylinder

Builds a cylinder,

          (Object Cylinder (Center x_value y_value z_value)
                           (MaterialName/SourceName value)
                           (Axis   x_value y_value z_value)
                           (Radius value)
                           (Height value)
                            <optional array param>
          )

Axis

The axis of the cylinder,
$x\_value,y\_value,z\_value\in \mathbb {R}$.

• Unit = Option::DimensionParam::LengthScale.

Height

The size along Axis,
$value\in \mathbb {R}$.

• Unit = Option::DimensionParam::LengthScale.

Radius

The radius of the cylinder,
$value\in \mathbb {R}$.

• Unit = Option::DimensionParam::LengthScale.

Sphere

Builds a sphere,

          (Object Sphere (Center x_value y_value z_value)
                           (MaterialName/SourceName value)
                           (Radius value)
                            <optional array param>
          )

Radius

The radius of the sphere,
$value\in \mathbb {R}$.

• Unit = Option::DimensionParam::LengthScale.

Outputs

This user input defines the outputs of the FDTD simulation.

   syntax:
          (Output OutputType (StartTime value) 
                             (EndTime   value)
                             (TimeStep  value)
                             (FileName  value)   
                             <other parameters> 
           )

• The leap frog algorithmn of FDTD results in field components being stored at different time points and different spatial points.
• Before any output computation below, each field component is spatially averaged (with the corresponding field components in its neighboring cells) to the same spatial point wihin the cell; in addition, time averaging is performed to obtain all field components at the same time point.

Common output parameters

StartTime

The start time of the output,
$value\in \mathbb {R}$

• Unit = Option::DimensionParam::LengthScale / c.

EndTime

The end time of the output,
$value\in \mathbb {R} _{>StartTime}$

• Unit = Option::DimensionParam::LengthScale / c.

TimeStep

The time step of the output,
$value\in \mathbb {R}$

• Unit = Option::DimensionParam::LengthScale / c.

FileName

The name of the file,
$value\in {s}tring$.

FieldTime

Output selected fields at a particular point,

          (Output FieldTime (Position    x_value y_value z_value)
                            (OutputField value)  
                            (StartTime   value) 
                            (EndTime     value)
                            (TimeStep    value)
                            (FileName    value)   
          )

• Click here for help on post-processing of output data.

Position

The position coordinate at which output is performed,
$x\_value,y\_value,z\_value\in \mathbb {R}$.

OutputField

These are the fields requested. Possible values are dependent on value of CellParam::CellType.

• See below for material dependent output definitions.

FieldSpace

Output selected fields in the whole compute region,

          (Output FieldSpace (OutputField value)  
                             (StartTime   value) 
                             (EndTime     value)
                             (TimeStep    value)
                             (FileName    value)   
          )

• The output is performed into a hdf5 file.
  • The suffix ".h5" is appended to FileName internally.
• If Option::CellParam::FieldType=complex, both real and imaginary parts of field are output.
• todo (in future release):
  • To allow the output of the field in a plane instead of the whole region.

FluxTimePlane

At each requested TimeStep, output the sum of the time average Poynting vector over the plane,

          (Output FluxTimePlane (Center          x_value y_value z_value)
                                (NormalDirection value)
                                (Size            value0 value1)  
                                (StartTime       value) 
                                (EndTime         value)
                                (TimeStep        value)
                                (FileName        value)   
          )

• Click here for help on post-processing of output data.
• Output as Power / $\left(Option::DimensionParam::LengthScale\right)^{2}$ where Power unit = Watts.

• The sum of the time average Poynting vector over the plane is defined as,

          ${\frac {1}{2}}\Delta _{1}\Delta _{2}\sum _{plane}{\hat {n}}\cdot {R}e[E(t,r)^{*}\times {H(t,r)}]$

• $\Delta _{1},\Delta _{2}$ are the Option::DimensionParam::Resolution along the two perpendicular axes parallel to the plane (see Size below).
• $\hat n$ is the unit normal axis of the plane (i.e. NormalDirection).
• FluxTimePlane requires CellParam::FieldType = complex.
  • since the calculation of the time average Poynting vector requires E(t,r) and H(t,r) to be complex fields.

Center

The center coordinate of the plane,
$x\_value,y\_value,z\_value\in \mathbb {R}$.

• Unit = Option::DimensionParam::LengthScale.

NormalDirection

The normal direction of the plane,
$value\in \{\,$positive_x, positive_y,positive_z,negative_x,negative_y,negative_z$\,\}$.

Size

The size along the two perpendicular axis parallel to the plane,
$value0,value1\in \mathbb {R}$.

• Unit = Option::DimensionParam::LengthScale.
• For NormalDirection = positive_x, negative_x, value0 corresponds to size along y-axis and value1 corresponds to size along z-axis.
• For NormalDirection = positive_y, negative_y, value0 corresponds to size along x-axis and value1 corresponds to size along z-axis.
• For NormalDirection = positive_z, negative_z, value0 corresponds to size along x-axis and value1 corresponds to size along y-axis.

• To do
  • Implement time average Poynting vector through a box.

FluxSpectraPlane

Output the transmitted power spectra through a plane,

          (Output FluxSpectraPlane (Center          x_value y_value z_value)
                                   (NormalDirection value)
                                   (Size            value0 value1)  
                                   (StartTime       value) 
                                   (EndTime         value)
                                   (TimeStep        value)   
                                   (StartFrequency  value) 
                                   (EndFrequency    value)
                                   (FrequencyStep   (Value value) (Force ) )
                                   (FileName        value)   
          )

• Click here for help on post-processing of output data.
• Output as Power / $\left(Option::DimensionParam::LengthScale\right)^{2}$ where Power unit = Watts.
• The first step in calculating the transmitted power spectra involves taking the Discrete Fourier Transform (DFT) of each field parallel to the plane from StartTime to EndTime at time increment TimeStep. The DFT for field F(t,r) at spatial point r and frequency $\omega$ is computed as

          $F(\omega ,r)=\sum _{t'=StartTime}^{EndTime}Re\left[F(t',r)\right]\cdot {e}^{j\omega {t'}}$

• At EndTime, the transmitted power at each frequency $\omega$ through the plane is obtained by summing the complex Poynting vector across the plane,

          $T(\omega )={\frac {1}{2}}\Delta _{0}\Delta _{1}\sum _{plane}{\hat {n}}\cdot \left[E(\omega ,r)^{*}\times {H}(\omega ,r)\right]$

• $\Delta _{0},\Delta _{1}$ are the Option::DimensionParam::Resolution along the two perpendicular axes parallel to the plane (see Size in FluxTimePlane).
• $\hat n$ is the unit normal axis of the plane (i.e. NormalDirection).
• FluxSpectraPlane can be performed with CellParam::FieldType = real.
  • since E(t,r) and H(t,r) in the DFT summation above can be real fields.

• The plane parameters Center, NormalDirection and Size are the same as FluxTimePlane.

EndTime, StartTime

• Additional restrictions (besides those specified above)
  • If EndTime > Option::DimensionParam::TimeParam::End or EndTime <= StartTime, error is thrown during initialization.
  • If StartTime < Option::DimensionParam::TimeParam::Start, error is thrown during initialization.

TimeStep (optional)

• See above for definition.
• TimeStep is optional for FluxSpectraPlane.
• Default value = FDTD time increment.

StartFrequency

The start frequency of the flux spectra calculation,
$value\in \mathbb {R} _{>0}$

• Unit = $2\pi c$ / Option::DimensionParam::LengthScale.

EndFrequency

The end frequency of the flux spectra calculation,
$value\in \mathbb {R} _{\leq {f}_{max}}$

• $f_{max}\,=\,1/\left(2\cdot {T}imeStep\right)$.
• Unit = $2\pi \,c$ / Option::DimensionParam::LengthScale.
• value is automatically adjusted to $f_{max}$ if $value\,>\,f_{max}$ (in order to avoid aliasing).

FrequencyStep

Value

The frequency step of the flux spectra calculation,
$value\in \mathbb {R} _{>0}$

• Unit = $2\pi \,c$ / Option::DimensionParam::LengthScale.

Force (optional)

• If Force is not specified, Value is automatically adjusted to 1 / (EndTime - StartTime).
• Force is disabled by default.
• Currently, you need to leave a space after Force (i.e. (Force ) is fine but (Force) will throw an error). This bug will be fixed in the next release.

Material specific output definitions

Linear material outputs

This section defines the material dependent parameters and computation related to outputs for CellParam::CellType = linear.

OutputField

The fields requested,
$value\in \{$all, ex,ey,ez,hx,hy,hz$\}$.

• OutputField can be any combination of above values. For example, (OutputField ex ey hz) will only output Ex-, Ey- and Hz-fields.
• H-fields are output as $H/{\sqrt {c\cdot \epsilon _{0}}}$, where H unit = amperes / meter.
• E-fields are output as $E/{\sqrt {c\cdot \mu _{0}}}$, where E unit = volts / meter.

Example user input file 1

Field image

Below is an image of the Ey field generated from running the following input file.

User input

  # This is a 2D simulation of a waveguide in the $\Gamma {-}K$ direction of
  # a hexagonal lattice photonic crystal consisting of air holes in Silicon.
  # The outer holes of the waveguide have a reduced radius.
  # 2D here means only one cell is used in the z-direction with Bloch kz = 0.

   (Option CellParam      (CellType linear)
                          (FieldType complex)
   )
   (Option BoundaryParam  (Type bloch pml bloch)
                          (PmlParam (Type upml) (Size (Y 1.2) )
                           // SigmaFactor here gives R[0]=Exp[-10] 
                          (SigmaFactor (Y 0.26) ) )
                          (BlochK (X 0.3660) (Z 0) )
   )
   (Option DimensionParam (Center 0 0 0)
                          (Resolution  0.05 0.05 0.125)
                          (Size 1 20.4 0.125)
                          (TimeParam (End 3500) )
   )
   // To enable CheckInput, remove "\".
   (\Option CheckInput)
   // File will be output to structure.h5 by default.
   (Option OutputStructure)
   // RelPermittivity here is from effective index of  Si.
   (Material Lossless (MaterialName Si)
                      (MaterialToken 3)
                      (RelPermittivity 8.6696892)
                      (RelPermeability 1)
   )
   (Source GaussTime (SourceName Source1)
                     (Duration 1200)
                     (Center 600)
                     (Width 125)
                     (Frequency 0.251)
                     (Phase 0)
                     (Amplitude 1)
                     (Direction 0 1 0)
   )
   (Output FieldTime (Position 0 0 0)
                     (StartTime 0)
                     (EndTime 3500)
                     (TimeStep  0.19)
                     (FileName field_time_center.dat)
                     (OutputField ey hz)
   )
   (Output FieldSpace (StartTime 1500)
                      (EndTime 1500.3)
                      (TimeStep 0.2)
                      (FileName field_space)
                      (OutputField ey)
    )
    // Build a source.
    (Object Point (SourceName Source1) (Center 0 0 0) )

    // Build  tri-lattice air-holes in Si.
    (Object  Block (MaterialName Si)
                   (Center 0 0 0)
                   (Axis0  1 0 0)
                   (Axis1  0 1 0)
                   (Axis2  0 0 1)
                   (Size 1 20.0 0.125)
    )
    //Note : vacuum is the default material in FDTD+ (it is generated internally).
    (Object Cylinder (MaterialName vacuum)
                     (Center 0 -8.66025403784 0)
                     (Axis 0 0 1)
                     (Radius 0.3)
                     (Height 0.125)
                     (Array  (NoPeriod 1 11 1)
                             (Period 1 1.73205080757 1)
                             (PeriodAxis0 1 0 0)
                             (PeriodAxis1 0 1 0)
                             (PeriodAxis2 0 0 1)
                     )
    )
    (Object Cylinder (MaterialName vacuum)
                     (Center -0.5 -9.52627944163 0)
                     (Axis 0 0 1)
                     (Radius 0.3)
                     (Height 0.125)
                     (Array (NoPeriod 2 12 1)
                            (Period 1 1.73205080757 1)
                            (PeriodAxis0 1 0 0)
                            (PeriodAxis1 0 1 0)
                            (PeriodAxis2 0 0 1)
                      )
    )
     // Create a waveguide.
    (Object Cylinder (MaterialName Si)
                     (Center 0 0 0)
                     (Axis 0 0 1)
                     (Radius 0.35)
                     (Height 0.125)
    )
    // These two instructions reduce the hole size of the outer row of holes
    // surrounding the waveguide.
    (Object Cylinder (MaterialName Si)
                     (Center -0.5 -0.866025403784 0)
                     (Axis 0 0 1)
                     (Radius 0.35)
                     (Height 0.125)
                     (Array  (NoPeriod 2 2 1)
                             (Period 1 1.73205080757 0)
                             (PeriodAxis0 1 0 0)
                             (PeriodAxis1 0 1 0)
                             (PeriodAxis2 0 0 1)
                      )
    )
    (Object Cylinder (MaterialName vacuum)
                     (Center -0.5 -0.866025403784 0)
                     (Axis 0 0 1)
                     (Radius 0.265)
                     (Height 0.125)
                     (Array  (NoPeriod 2 2 1)
                             (Period 1 1.73205080757 0)
                             (PeriodAxis0 1 0 0)
                             (PeriodAxis1 0 1 0)
                             (PeriodAxis2 0 0 1)
                      )
    )

Example Material::Prp parameters

Gold (Au) Prp parameters

      (Material Prp (MaterialName Au)   // MaterialName value is user defined.
                    (MaterialToken 8)   // MaterialToken value is user defined.
                    (EpsilonInfinity 1)
                    (RelPermeability 1)
                    (PoleResiduePair
                        (Pair1 (2.2217e16, -6.3742e18) (-1.4982e14, 3.0413e11) )
                        (Pair2 (-6.4292e15, 1.3091e18) (-2.8617e13, -3.4173e13) )
                        (Pair3 (2.1046e16, 1.4907e15)  (-3.2597e15, -4.3949e15) )
                        (Pair4 (1.0807e15, 5.889e14)   (-4.4078e14, -3.9961e15) )
                        (Pair5 (-7.173e15, 2.5124e18)  (-3.3045e12, -1.0357e13) )
                        (Pair6 (-5.3709e15, 2.517e15)  (-1.9179e15, 7.2146e15) )
                        (Pair7 (-7.8496e14, 3.3901e18) (2.39e13, -5.5125e12) )
                        (Pair8 (-7.1488e15, 5.1946e18) (4.5398e13, 1.2772e12) )
                    )
      )

Silver (Ag) Prp parameters

      (Material Prp (MaterialName Ag)   // MaterialName value is user defined.
                    (MaterialToken 9)   // MaterialToken value is user defined.
                    (EpsilonInfinity 1)
                    (RelPermeability 1)
                    (PoleResiduePair
                          (Pair1  (2.004847623185898e16, 6.373518292877005e18)
                                  (-1.450986957494932e13, -1.471286657124635e13)
                          )
                          (Pair2 (4.419387050893835e15, -2.455715772083395e16)
                                 (-1.617050856766036e15, -7.484019360934898e14)
                          )
                          (Pair3  (-1.253762773346362e16, 1.111152157244156e18)
                                  (-1.621233619602248e16, -2.197813817662167e15)
                          )
                          (Pair4 (2.397213142810919e15, -2.295684023806528e16)
                                 (-3.358496690509831e15, -4.204829676212065e15)
                          )
                          (Pair5 (-1.287525821691552e16, -2.726588225501935e15)
                                 (-1.826237040321377e14, -1.915550886065391e14)
                          )
                          (Pair6 (3.258726396082003e15, -1.140962611954672e17)
                                 (-9.157642928770140e15, -1.161886436310815e16)
                          )
                          (Pair7 (1.485343112935528e15, -4.055619609190444e14)
                                 (-4.591136080183809e14, -6.063844504839434e15)
                          )
                          (Pair8 (3.256535035230075e13, -8.975107704148206e13)
                                 (-2.459619405093500e14, -6.564514351528399e15)
                          )
                    )
      )

Silicon (Si) Prp parameters

      // Fit for wavelengths up to ~800nm.
      (Material Prp (MaterialName Si)   // MaterialName value is user defined.
                    (MaterialToken 2)   // MaterialToken value is user defined.
                    (EpsilonInfinity 1)
                    (RelPermeability 1)
                    (PoleResiduePair
                       (Pair1 (-1.2393e17, -1.5487e17)
                               (-6.3301e14, -6.4586e15)
                       )
                       (Pair2 (-3.6772e16, 7.6757e16)
                               (-7.3714e14, -6.1149e15)
                       )
                       (Pair3 (1.7097e17, 1.0258e17)
                               (-5.666e14, -6.3804e15)
                       )
                       (Pair4 (-1.2185e16, 7.215e15)
                               (-3.3952e15, -8.0514e15)
                       )
                       (Pair5 (3.2464e15, -3.7598e15)
                               (-1.952e14, 5.126e15)
                       )
                       (Pair6 (-1.0e11, 7e13)
                               (-1e15, 3e15)
                       )
                       (Pair7 (-5000000, 2.93e12)
                               (-3e14, 1.9435e15)
                       )
                       (Pair8 (-5000, -2.5e12)
                               (-1.98e14, 1.1e15)
                       )
                       (Pair9 (-5000, -4.6e12)
                               (-2.5e14, 2.5e15)
                       )
                       (Pair10 (-5000, -1e12)
                               (-2.98e14, 4.75e14)
                       )
                    )
      )

Reading output data with Mathematica

• This patch of Mathematica commands can help in reading real/complex output data from Output::FieldTime, Output::FluxTimePlane and Output::FluxSpectraPlane.
• In FDTD Plus, complex data is output by C++ std::complex::operator<< in the cartesian format (i.e. "a + j b" is output as "(a, b)").

   (* Set the directory location. *)
   In[5]:= SimNo = "sim68";
   In[6]:= DirName = "C:\Research\FdtdPlusSimulations\\" <> SimNo ; SetDirectory[DirName]
   Out[6]= "C:\\Research\\FdtdPlusSimulations\\sim68"
   
   (* File name.  *)
   In[7]:= FileName = "field_time.dat";
   (* Read the first few lines (optional step; helps you see the headers). *)
   In[8]:= ReadList[FileName, Word, 4, RecordLists -> True, WordSeparators -> {"\n"}]
   Out[8]= {{"time_value	hz	ey	"}, {"0	(0,0)	(0,0)	"}, {"0.189038	\
                (-0.000608598,0.000132897)	(0.000814581,-0.000388972)	"}, {"0.378076	\
                (-0.000812379,1.68512e-05)	(0.00168293,-0.000526529)	"}}

   (* Read the file. *)
   (* Data is stored in list MyData. *)
   (* Note header above (i.e. "time value hz ey") is removed. *)
   In[9]:= MyString = Import[FileName, "Text"]; MyString =
   StringReplace[MyString, {"(" -> "", ")" -> "", "," -> "\t"}]; MyStream =
   StringToStream[MyString]; MyData = {}; Read[MyStream, String]; MyData =
   ReadList[MyStream, Real, RecordLists -> True,
   RecordSeparators -> {"\n"}]; Close[MyStream]; Clear[MyString];

   (* Output a line from MyData. *)   
   In[10]:= MyData$[[2]]$
   Out[10]= {0.189038, -0.000608598, 0.000132897, 0.000814581, -0.000388972}

Reading output data with Matlab

• This Matlab command can help in reading complex output data from Output::FieldTime, Output::FluxTimePlane and Output::FluxSpectraPlane.
• In FDTD Plus, complex data is output by C++ std::complex::operator<< in the cartesian format (i.e. "a + j b" is output as "(a, b)").

   [t eyr eyi hzr hzi] = textread('field_time_center.dat',
                                  '%f%*c%f%*c%f%*c%*c%f%*c%f%*c','headerlines',1);