June 30th, 2011TutorialsJuan Sebastian Totero Gongora 0 Comments

The first simulation

We begin by considering the simplest situation of a 2D dipole source propagating in the free space. From the Nanocpp root directory, move in the default directory:

cd newsim/

There are several files:

newsim $ ls
newsim $ angles.in· CVS· geometry.in· input.in· makefile· myfile.cc· spectr.m

for this simulation, we need to edit the file input.in, which contains all the high level options for the Nanocpp code. Open your favorite editor (we are linux friendly and assume the use of emacs for this tutorial) and edit input.in.

As you see, the file is divided into two main sections:

######### COMPILE-TIME PARAMETERS ##########


########## RUN-TIME PARAMETERS ###########

the former contains all the parameters that need to be known at compile time (e.g., optimizations, dimensions, source types, …), while the latter has all the data for running a Nanocpp simulation (geometry, number of cells, …). One of the strength of Nanocpp, which allows it to scale over hundreds of thousands of processors, lies in its templated+macro structure. Nanocpp, in fact, self-generates the C++ code from the parameters of the compile-time section and produces an executable that is always optimized and free from any if statement.

We begin by selecting the dimensionality of the problem:

N=2···························· # dimensionality of the problem (2,3)

and the source type:

spatial=point+Ex······································· # + ·

Nanocpp embodies a high level C++ parser class that allows to combine all the various primitives with a high-level syntax. The data will then be interpreted by the parser and traduced into low-level optimized C++ instructions. In this case, to specify a point source we use “point”, which stands for point source, followed by the fields we want to activate, in this case only Ex. The order of this instructions is irrelevant, and the instructions are case insensitive. This give the maximum degree of flexibility.

We complete the compile-time section by selecting the measurements that we want to enable in the code. In this case, we choose just a point probe:

measurements=probes+Ex················· # probes (Ei and Hi) + farfield (only in 3D) + energy + poynting

We now type


and move to the run-time section of input.in. As a first step, we define the computational domain of our system:

number_of_cells=100 100 60············· # x y z number of unit cells
box_size=2 2 2························· # physical size of simulation box along x y z

The number of cells gives the number of YEE cells of our simulation, while the Box size if the physical size of the box. In Nanocpp, SI units can be used, as well as adimensional units (like in this example). These are defined as follows:


being a an arbitrary length scale. We then move to the Time & I/O parameters:

##### TIME & I/O
base_dir=./res················· # directory where to save results
CFL=2.························· # CFL Condition
total_steps=10000······················ # steps used for performing final averages (total steps = equilibration + measure)
io_screen=1000················· # I/O on screen
chkpoint_every=10000··········· # checkpoint step
chk_load=0····················· # 1 true, 0 false
chk_numb=0····················· # chkpoint number: 0 or 1

Parameters total_steps, io_screen and chk_every are the total number of steps, the frequency of I/O on screen and the checkpoint, respectively, while CFL is a scaling factor to determine the timestep dt:


with dt_{max} the maximum time step allowed for a stable leapfrog time-marching algorithm. The standard value CFL=2 is equivalent to a Courant stability factor of 0.5, as conventionally used in FDTD simulations (see, e.g., here).

The parameter base_dir is the root directory where to save all results (the dir is automatically created if does not exist). The checkpoint is enabled by the Boolean variable chk_load, and the checkpoint number to be loaded is selected by chk_numb, which has been introduced as a mechanism of redundancy to avoid· crashing problem during the writing of the checkpoint. In particular, Nanocpp uses a toroidal system of two checkpoints, which are labeled as “0” and “1”. In the unlikely event of a system crash during a checkpoint (e.g., checkpoint 0), the simulation can be safely restarted by using the previous file (checkpoint 1), which has been written before.

We can now select the time waveform of the source:


waveform=cw+0.25+3+1e6 ·················· # + ·
# gaussian + wavelength + waist + peak time position
# gaussian2 + wavelength + waist in wl
# spl + wavelength + starting time
# cw + wavelength + m-duration + n-duration

that is specified by selecting a type (gaussian, single-pulse or spl, and continuous or cw) followed by the· source parameters as explained in the comments. For this simulation· we use a cw source, with wavelength 0.25, 3 rising peaks and 1e6 peaks in the steady-state condition (before switching it off).

The last parameter we need to specify is the· location of the probes:

probe_list=1.2+1.3 # probe point list (x y z) separated by +; probe are saved as p1_<field>.bin,p2_<field>.bin,...

# probes owned by one processor are grouped in p_<0>_<1>_..._<N>_<field>.bin and saved in [0-t0,1-t0,...,N-t0;0-t1,1-t1,...,N-t1; ... ]

This is a list of (x,z) positions separated by a ‘+’. More than one probe can be added to this list; probes are grouped by processors and saved in different files whose name and structure is explained in the comments. Check that the probe point list contains coordinates inside the box_size. All external coordinates are simply ignored by Nanocpp and will not produce any output.

We are now ready for running our first simulation. Type:

mpiexec -n 8 ./nano input.in

and Nanocpp will start its computations. Once done, open the probe file inside the res/probes/ directory with your favorite numerical program (we use matlab) and plot the result:


Energy and Poynting

To enable the calculation of energy and poynting vectors, we first add them in the measurement section of the compile-time parameters part of input.in:


we then re-compile the code:

make clean;make

and finally select the frequency of the energy & poynting I/O operations:

erg_io=10000··································· # frequency of io save

in this case every 10000 simulation steps. For nondispersive materials, the energy calculated by Nanocpp is the electromagnetic energy density u:


For energy calculations in dispersive materials, please see the section on dispersive simulations on this page.

We can then launch Nanocpp (now employing 16 processors):

mpiexec -n 16 ./nanocpp input.in

and look at the files “erg_<num>.bin” for the energy and “px_<num>.bin”, “pz_<num>.bin” for the poynting vector \mathbf{S}=\mathbf{E} \times \mathbf{H}. The next figure shows the results visualized with Matlab (pcolor for energy and quiver plot for poynting):



Other than just point sources, NANOCPP supports TFSF sources for plane waves and Gaussian beams. To enable a plane TFSF source with Ex polarization, go to the sources section in the compile-time parameters and type:


then, in the tfsf source specifiers, put the parameters of the box containing the plane wave source:

tfsf_box=0.6 1.4 0.5 1.5 0.5 0.5······· # xsta xend ysta yend zsta zend separated by just one space
tfsf_power=1e-3	     	     	     	  # average power of the source [W] at central frequency

The box is defined by the edges along x,y,z directions. The average power, conversely, represents the input power of the TFSF source in SI. After re-compiling, we run the simulation. Here is a plot of the energy obtained:


NANOCPP uses a slight modification to the standard TFSF formulation, which allows the plane wave to be numerically exact, up to double precision. In order to use a different time-waveform for the source, and pick up a Gaussian pulse:


which selects a Gaussian envelope with·\lambda=0.25 and bandwidth \Delta\lambda=0.05. We then increase the TFSF box size along z:

tfsf_box=0.6 1.4 0.2 1.8

to better see the tail of the energy field of the source. This is the energy distribution after 600 timesteps:


Localized TFSF sources are available, as well. To enable them, add the single specifier in the sources data list:


and then, in the TFSF box specifiers, select the size along x and the starting position along z:

tfsf_box=0.9 1.1 0.3 0.1 0.5 0.5

in this case, we selected a source at z=0.3, whose extent along x goes from 0.9 to 1.1. Then re-compile and run. This is the energy after 600 timesteps:


When using a TFSF source, NANOCPP assumes that the source propagates in the background material of the simulation, which is general is air. To use the TFSF in the case where there is more than one geometry, and the source is embedded into a material of dielectric constant \epsilon_r, you can use the function:

void set_tfsf_background(double epsilon_r)

which sets the value of dielectric constant for the TFSF source.

Periodic Boundaries

By deafult, NANOCPP assumes non-periodic boundaries terminated with UPML. Periodic boundaries are enabled by adding the directives periodic_<direction> (direction=x,y, or z) in the problem_type compile-time parameter in input.in:


different periodic directives periodic_x, periodic_y and periodic_z can be added together to specify more than one periodic boundary at the same time. In order to enable the periodicity, set also the corresponding UPML thickness to 0, e.g., along x:

upml=0 20

To enable a TFSF plane wave source, use the:


directive, and specify a TFSF window from 0 to L-dx, being L the size of the computational domain along the chosen periodic direction and dx the grid spatial resolution.
The following input.in and geometry.in enable a 2D TFSF plane wave source periodic along x in air:


Geometric Shapes

Material parameters and shapes are contained in the file specified by the media_filename keyword:

number_of_cells=200 200 100············ # x y z number of unit cells
box_size=2 2 2························· # physical size of simulation box along x y z media_filename=geometry.in

The file geometry.in is divided into two main sections:· material parameters· and geometry. The standard file comes with two materials only, air and Silicon, which is enough for this example. Concerning the geometry section, the following shape primitives are supported:

### sphere····· material id· (center······ x y z)· (radius) ### cylinder··· material id·· (center······ x y z) (length)· (radius) (orientation) ### cuboid····· material id·· (lower corner x y z)· (upper corner x y z) ### elipsoid··· material id · (center······ x y z)· (axis········ x y z)

The first parameter of every shape is the material id. This is an integer number that connects the shape to a material defined in the material parameters section. The other geometrical parameters are self-explanatory. In the geometry.in file there is only one shape, a sphere centered in (1 1) which fills the whole space with air (material id 0).

#background sphere 0 1 1 100

In this example we will add an ellipse:

ellipsoid 1 1 1 0.3 0.2

made of Silicon (material id 1), centered at (1 1) and with axes 0.3 along x and 0.2 along z. We then run NANOCPP again, and look at the output of NANOCPP:

filename: geometry.in
number of shapes: 1
sphere:·· type (0);······ center (1, 1);······ radius (100)
materials properties:
type: 0, eps_inf: 1
type: 1, eps_inf: 12

saved (z,x) 201 200 int on type_Ex.bin
saved (z,x) 200 201 int on type_Ez.bin

after every initialization, NANOCPP writes the spatial distribution of the materials id in the files (rootdir)/type_E(axis).bin. Data is saved as binary integers. Since the YEE grid uses different spatial positions for each electromagnetic field, a different file for each electromagnetic component is created and saved. In what follows we display the Ex media distribution, obtained from the 201 x 200 integer data of res/type_Ex.bin:


Geometrical operations on shapes are also supported. In NANOCPP, in particular, the rotation primitive is available:

### rotation··· (center· x y z) (rotation angles xz yz xy)

and it applies to all shapes defined below it. A rotation overwrites all the previous rotation directives. This allows to apply this primitive to single and multiple shapes, thus achieving the maximum degree of flexibility when defining a geometry. In the following example, we rotate the ellipse by 45 degrees along xz:

rotation 1. 1. 0.7854 
ellipsoid 1 1 1 0.3 0.2

Re-run the simulation again. The file type_Ex.bin now shows a rotated geometry:


Thanks to Pritam Roy for the FDTD simulations

Supercontinuum sources

A supercontinuum source is represented in NANOCPP as a frequency comb:


characterized by N equispaced frequencies centered at \omega_0. To use this sources, select the fcb option in the sources list:

# fcb + wl0 + wl1 + N + m-duration + n-duration 
# (fill wl0, wl1 with N mnm cw)

this is the output provided by these
input.in and geometry.in files, where a frequency comb of 5 wavelengths in the range \lambda\in[1,1.5] has been selected:



Parallel Far-Field calculation is available in Nanocpp by specifying the compile-time parameter farfield in the MEASUREMENT section of input.in:

##### MEASUREMENT measurements=probes+Ex+farfield+energy 
# probes (Ei and Hi) + farfield + energy + poynting

Additional optional parameters bscatt and fscatt can be added to calculate backward and forward total electromagnetic far-fields. Please remember to compile the code with:

make clean;make

if not done before. Nanocpp implements the Time-Domain Near-To-Far-Field (NTFF) algorithm, explained in Chap. 8, Section 8.6 of the Taflove’s Book. This allows to calculate the time evolution of the scattered electromagnetic fields at specific angles selected by the user. Several original optimization strategies have been developed in order to optimize both memory requirements and execution time.

In order to set up the Far-Field computation routine, the user needs to specify a cuboid surface, to be used as a source for the surface equivalence theorem (see Chapter 8 of Taflove’s for more details):

farfield_box=0.2 1.8 0.2 1.8 0.5 0.5
# far field box xsta xend ysta yend zsta zend separated by just one space

Angles [\theta,\phi], where the scattered fields have to be calculated, should be saved in a binary file that needs to be specified in the file input.in:

n_angles=100        # number of angles (theta,phi) to process
ff_file=angles.in   # binary filename containing angles [theta,phi, theta,phi, ...]

In this case the programs looks for the angles in the file angles.in. At the end of the simulation, Nanocpp saves the scattered electric field E_\theta and E_\phi in the files <resdir>/eth.bin and <resdir>/ephi.bin, respectively. These files input.in, geometry.in and angles.in calculate the scattered far field from an electric dipole radiator oscillating parallel to \mathbf{x}, with a current J\propto\frac{\sin(\omega t)}{\lambda}. For this type of source, (see Chaper 9 of Jackson’s Book) the far field intensity is I\propto\frac{\sin\theta^2}{\lambda^4}. Figures below show the results of Nanocpp:



Dispersive Materials

NANOCPP allows to add materials with an arbitrary dispersion, thanks to a new formulation based on an ultrafast unconditionally stable algorithm with the same veocity of explicit schemes. In NANOCPP, the dispersion \epsilon(\omega) of a generic material is expressed in fractional form:

  \epsilon(\omega)=\frac{\sum_n^N b_n(i\omega)^n}{\sum_n^N a_n(i\omega)^n}

with arbitrary N and polynomial coefficients a and b. Once a material is defined, NANOCPP checks if the behavior beyond the critical frequency of the simulation gives instability. In NANOCPP, we use the convention that negative values of the imaginary part of  \epsilon(\omega) mean absorption, while positive values mean gain.

NANOCPP is also equipped by a function:

set_index(int id,double w0,double er,double ei,double der,int iter=1000);

which assigns to the material id, at the frequency w0, the values of real er and imaginary ei permittivity specified by the user. Optionally, the user can also specify der=\frac{\partial er(\omega)}{\partial\omega}|_{\omega_0}, and a number of iterations iter. If this data is given, the function uses a Newton-Raphson algorithm to match both the dispersion values and the dispersion derivative, providing a good fit for a larger range of wavelengths. The function set_index uses a specific form of Lorentz oscillator that can represent any material, including metals with negative values of the dielectric constant, without giving any sort of instability. When using the function set_index, NANOCPP prompts at the output the a and b coefficients of the updated material.

The dispersive \epsilon_r(\omega) and \epsilon_i(\omega) are saved in the res directory for frequencies varying \pm 1 decade with respect to the \omega_0 of the source.

Dispersive parameters are specified in NANOCPP by writing the number of poles, which is equal to N-1, and coefficients b and a in the geometry file. Please refer to the standard geometry.in for more information about the possible syntax variations. Here is a 4 poles representation of Silicon from 300nm to 1000nm:

@ 1 12 4 1.2791e+64 7.2843e+47 4.9912e+32 1.1920e+16 1.0000e+00 1.0913e+63 6.1040e+46 6.8576e+31 2.0394e+15 1.0000e+00

Dispersive materials are enabled by adding the word “dispersive” in the compile time parameters of input.in:


the following input.in and geometry.in simulate a classical problem of reflection and absorption in a Si substrate.

In dispersive materials, the dielectric displacement D(\omega) in the frequency domain can be written as follows:


The Maxwell equation for the time dynamics of the displacement is then:

  \nabla\times H(\omega)=i\omega\epsilon_0\epsilon_\infty E(\omega)+J(\omega),

with J(\omega)=i\omega\epsilon_0\Delta\epsilon(\omega)E(\omega). The Pointing theorem in this case read as follows:

  \partial_t u+\nabla\cdot S=-j(t)\cdot e(t),

with u=\frac{1}{2}\left[\epsilon_0\epsilon_\infty e(t)^2+\mu_0 h(t)^2 \right] being the energy density of the field in the dispersive case, and j(t) the inverse Fourier transform of the current J. In NANOCPP, the value \epsilon_\infty is automatically computed by the program and updated after each call to a set_index function, prior to any calculation of the energy density u(t).

Harmonic Analysis of Power Flow

In many cases is desirable to measure the reflection/transmission at different frequencies and obtain the corresponding spectra. One approach is to launch different simulations by using cw sources at different frequencies. This solution takes typically time especially when high spectral resolution is needed. When the dispersion curve of the material is known, a much more effective approach can be developed by using a broadband SPL input source, which behaves as a dirac delta in time an probe the density of states of the system in a continnum set of frequencies. The main issue of this approach is how to measure the energy, which is a nonlinear operator involving products of \mathbf{E} and \mathbf{H} fields. The presence of products require to store in memory the time evolution each component of the field in each point of integration. The corresponding data can be quite large, especially in 3D. To avoid this problem an enable an efficient analysis in NANOCPP, we took inspiration from JPEG compression of data and use an orthogonal, complete basis of cosine functions:

  c_n=A_n\cos(n\pi x/L)

with normalization constants  A_n defined on the integration length  L, to store all temporal information of the electromagnetic fields needed for power flow computations. NANOCPP automatically calculates the time expansion coefficients and save them in an appropriate matrix file, which can be then easily analyzed. This main advantage of this approach is that it typically converge with just a few cosine frequencies (less than 10 for accuracy < 1%). The following files:





solves a simple problem of a fabry-perot resonator and calculates the transmission spectrum. For periodic systems, 1 mode is enough for convergence.