Meep C-plus-plus Tutorial

From AbInitio

(Difference between revisions)
Jump to: navigation, search
Revision as of 18:57, 25 October 2005 (edit)
Ardavan (Talk | contribs)

← Previous diff
Current revision (15:05, 18 April 2017) (edit)
Ardavan (Talk | contribs)

 
Line 1: Line 1:
{{Meep}} {{Meep}}
-Instead of using the libctl/Scheme interface described in the [[Meep tutorial]], one can also access Meep as a C++ library by writing a C++ program that links to it. The C++ interface provides the greatest level of flexibility in designing FDTD simulations.+Instead of using the libctl/Scheme interface described in the [[Meep tutorial]], Meep is also callable as a C++ library by writing a C++ program that links to it. The C++ interface provides the greatest level of flexibility in setting up FDTD simulations. There are 19 examples of such C++ programs in the <code>tests/</code> directory of the Meep codebase which cover a wide range of Meep's functionality. These are a good reference. Finally, we should also note that, while Meep is nominally in C++, it is perhaps better described as "C+". That is, most of the coding style is C-like with a few C++ features.
-(This is, in fact, the original way in which we used Meep, and the way in which all of the Meep test programs, in its <code>tests/</code> directory, are written.)+== Differences from libctl ==
-We should also note that, while Meep is nominally in C++, it is perhaps better described as "C+". That is, most of the coding style is C-like with a few C++ features.+The C++ interface has several differences from the libctl interface besides the obvious difference in syntax.
 + 
 +The most notable difference is that, while the libctl interface puts the origin <math>(0,0,0)</math> at the ''center'' of the computational cell, the C++ interface by default puts the origin at the ''corner'' of the computational cell. That is, an <math>L\times L\times L</math> cell goes from <math>(-L/2,-L/2,-L/2)</math> to <math>(L/2,L/2,L/2)</math> in libctl, but from <math>(0,0,0)</math> to <math>(L,L,L)</math> in C++. (You can change this by calling <code><nowiki>grid_volume::shift_origin</nowiki></code>.)
== Overview == == Overview ==
Line 15: Line 17:
to include the definitions of the Meep routines. Later, when you compile (see below), you must also link to the Meep libraries. to include the definitions of the Meep routines. Later, when you compile (see below), you must also link to the Meep libraries.
-Your main program should then initialize Meep, and will generally then define a computational <code>volume</code> and the associated <code>structure</code> (describing the geometry and materials), initialize the <code>fields</code>, add sources, and then time-step. In short:+Your main program should then initialize Meep, and will generally then define a computational <code>grid_volume</code> and the associated <code>structure</code> (describing the geometry and materials), initialize the <code>fields</code>, add sources, and then time-step. In short:
int main(int argc, char **argv) { int main(int argc, char **argv) {
- initialize mpi(argc, argc); // do this even for non-MPI Meep+ initialize mpi(argc, argv); // do this even for non-MPI Meep
double resolution = 20; // pixels per distance double resolution = 20; // pixels per distance
- volume v = vol2d(5,10, resolution); // 5x10 2d cell+ grid_volume v = vol2d(5,10, resolution); // 5x10 2d cell
structure s(v, eps, pml(1.0)); structure s(v, eps, pml(1.0));
fields f(&s); fields f(&s);
Line 59: Line 61:
The particular Fabry-Perot cavity we will investigate consists of an air region bounded by two distributed Bragg reflectors (quarter-wave stack) of epsilon 12 and 1. We choose the size of the defect to be twice as large as the air thickness in the quarter-wave stack, so that a defect mode is found near midgap. This structure will include N = 5 periods in the Bragg reflector on either side of the defect. We can use a larger N but the quality factor may then be too large to compute. The parameters are set up as follows: The particular Fabry-Perot cavity we will investigate consists of an air region bounded by two distributed Bragg reflectors (quarter-wave stack) of epsilon 12 and 1. We choose the size of the defect to be twice as large as the air thickness in the quarter-wave stack, so that a defect mode is found near midgap. This structure will include N = 5 periods in the Bragg reflector on either side of the defect. We can use a larger N but the quality factor may then be too large to compute. The parameters are set up as follows:
- const double eps1 = 12.0;+ const double eps1 = 12.0; // epsilon of layer 1
- const double eps2 = 1.0;+ const double eps2 = 1.0; // epsilon of layer 2
const double grating_periodicity = 1.0; const double grating_periodicity = 1.0;
- const double d1 = sqrt(eps2) / (sqrt(eps1)+sqrt(eps2));+ const double d1 = sqrt(eps2) / (sqrt(eps1)+sqrt(eps2)); // quarter wave stack dimensions
const double d2 = sqrt(eps1) / (sqrt(eps1)+sqrt(eps2)); const double d2 = sqrt(eps1) / (sqrt(eps1)+sqrt(eps2));
const double half_cavity_width = d2; const double half_cavity_width = d2;
Line 91: Line 93:
initialize mpi(argc,argv); initialize mpi(argc,argv);
const double amicron = 10.0; const double amicron = 10.0;
- const volume vol = vol1d(2*z_center,amicron);+ const grid_volume vol = vol1d(2*z_center,amicron);
structure s(vol,eps,pml(pml_thickness)); structure s(vol,eps,pml(pml_thickness));
fields f(&s); fields f(&s);
-Note the constructor for the volume class in 1D which takes as parameters the size of the cell and the resolution. The sources in Meep are excitations of the polarization vector in <math>\mathbf{D}=\epsilon\mathbf{E}+\mathbf{P}</math>. The polarization can be any one of the six cartesian or cylindrical fields. There are a variety of sources including dipole and current sources, gaussian pulses and a continuous wave sources. For now, we use a gaussian source centered at the midgap frequency with a narrow width, along with the start time and end time of the source.+Note the constructor for the <code>grid_volume</code> class in 1D which takes as parameters the size of the cell and the resolution. The sources in Meep are excitations of the polarization vector in <math>\mathbf{D}=\epsilon\mathbf{E}+\mathbf{P}</math>. The polarization can be any one of the six cartesian or cylindrical fields. There are a variety of sources including dipole and current sources, gaussian pulses and a continuous wave sources. For now, we use a gaussian source centered at the midgap frequency with a narrow width, along with the start time and end time of the source.
const double w_midgap = (sqrt(eps1)+sqrt(eps2))/(4*sqrt(eps1)*sqrt(eps2)); const double w_midgap = (sqrt(eps1)+sqrt(eps2))/(4*sqrt(eps1)*sqrt(eps2));
Line 114: Line 116:
int ttot = int(400/f.dt)+1; int ttot = int(400/f.dt)+1;
complex<double> *p = new complex<double>[ttot]; complex<double> *p = new complex<double>[ttot];
- complex<double> *amps = new complex<double>[maxband];+ complex<double> *amps = new complex<double>[maxbands];
- double *freq_re = new double[maxband];+ double *freq_re = new double[maxbands];
- double *freq_im = new double[maxband];+ double *freq_im = new double[maxbands];
The variable <code>maxbands</code> is the number of mode frequencies we will be looking for while <code>ttot</code> represents the total number of timesteps for which we will monitor the field component at. For this simulation given the large quality factor that we expect, we will time step for a long period of time which in this case is 400 unit seconds. Now we are ready to time step and collect the field information: The variable <code>maxbands</code> is the number of mode frequencies we will be looking for while <code>ttot</code> represents the total number of timesteps for which we will monitor the field component at. For this simulation given the large quality factor that we expect, we will time step for a long period of time which in this case is 400 unit seconds. Now we are ready to time step and collect the field information:
Line 126: Line 128:
} }
-The get_field function does exactly that, obtains the value of the field at a given position within the computational cell using linear interpolation where necessary. Now we have all the information we need and must obtain the frequencies of the modes which we do by invoking a special tool for harmonic inversation, <code>harminv</code>, to extract the real and imaginary frequencies:+The get_field function does exactly that, obtains the value of the field at a given position within the computational cell using linear interpolation where necessary. Now we have all the information we need and must obtain the frequencies of the modes which we do by invoking a special tool for harmonic inversion, <code>harminv</code>, to extract the real and imaginary frequencies:
int num = do_harminv(p, ttot, f.dt, 0.8*w_midgap, 1.2*w_midgap, maxbands, amps, freq_re, freq_im); int num = do_harminv(p, ttot, f.dt, 0.8*w_midgap, 1.2*w_midgap, maxbands, amps, freq_re, freq_im);
-The integer returned by <code>harminv</code> is the number of mode frequencies obtained from the input data <code>p</code>. The particular call to <code>harminv</code> included passing the arrays by value and telling <code>harminv</code> to look for frequencies within 20% of the mid gap frequency up to a maximum of 5 bands. At this point, the necessary information to compute the quality has been stored in the freq_re and freq_im arrays and we compute the quality factor using the formula, Q<math>=-2*\omega_r_e_a_l/omega_i_m_a_g</math>. All that is left to do, is to print these results with the quality factor:+The integer returned by <code>harminv</code> is the number of mode frequencies obtained from the input data <code>p</code>. The particular call to <code>harminv</code> included passing the arrays by value and telling <code>harminv</code> to look for frequencies within 20% of the mid gap frequency up to a maximum of 5 bands. At this point, the necessary information to compute the quality has been stored in the <code>freq_re</code> and <code>freq_im</code> arrays and we compute the quality factor using the formula, Q<math>=-\omega_r/2\omega_i</math>. All that is left to do, is to print these results with the quality factor:
- master_printf("frequency,amplitude,quality factor\n"); + master_printf("frequency,amplitude,quality factor\n");
- for (int i=0; i<num; ++i) + for (int i=0; i<num; ++i)
- master_printf("%0.6g%+0.6gi,%0.5g%+0.5gi,%g\n",freq_re[i],freq_im[i],real(amps[i]),imag(amps[i]),-2*freq_re[i]/freq_im[i]);+ master_printf("%0.6g%+0.6gi,%0.5g%+0.5gi,%g\n",freq_re[i],freq_im[i],real(amps[i]),imag(amps[i]),-0.5*freq_re[i]/freq_im[i]);
- return 0;+ return 0;
} }
-That's it, we are done. The output for this program is <code>1.30435e+07</code>. The large quality factor which is to be expected since our system includes a relatively high number of bragg reflectors (N=5) in each direction and we are exciting the mode at mid gap.+That's it, we are done. The output for this program is <code>3.26087e+06</code>. The large quality factor is to be expected since our system includes a relatively high number of bragg reflectors (N=5) in each direction and we are exciting the mode at mid gap. Suppose we want to see what the resonant mode looks like, say over one period. The following block of code time steps the field for exactly one period while outputting the Ex field for the computational cell at each time step:
 + 
 + double curr_time = f.time();
 + while (f.time() < curr_time + 1/w_midgap) {
 + f.output_hdf5(Ex, vol.surroundings());
 + f.step();
 + }
 + 
 +After we convert the hdf5 files to PNG format and superimpose the images on the dielectric background to produce an animated, we obtain the following:
 + 
 +[[Image:Fabryperot.gif|center|Resonant mode of fabry-perot cavity]]
 + 
 +== Compiling ==
 + 
 +In order to compile your code and link against the Meep library, you must link in several additional libraries that Meep requires. There are three possible ways to do this:
 + 
 +1. After compiling the Meep package following the instructions elsewhere in this wiki, place myfile.cpp (where myfile.cpp is a MEEP C++ file of your choosing) in the <code>tests/</code> subdirectory, cd to the same directory, and invoke:
 + 
 + make myfile.dac
 + 
 +Run the resulting executable via:
 + 
 + ./myfile.dac
 + 
 +2. Use the [[w:pkg-config]] program, which is installed by default on most GNU/Linux systems. Then, to compile a file <code>foo.cpp</code> using Meep, do:
 + 
 + g++ `pkg-config --cflags meep` foo.cpp -o foo `pkg-config --libs meep`
 + 
 +(Naturally, replace <code>g++</code> with the name of your C++ compiler if you are not using the GNU compilers.)
 + 
 +3. Compile with g++, this time invoking each library separately:
 + 
 + g++ -malign-double foo.cpp -o foo -lmeep -lhdf5 -lz -lgsl -lharminv -llapack -lcblas -latlas -lfftw3 -lm
 + 
[[Category:Meep]] [[Category:Meep]]
 +[[Category:Meep examples]]

Current revision

Meep
Download
Release notes
FAQ
Meep manual
Introduction
Installation
Tutorial
Reference
C++ Tutorial
C++ Reference
Acknowledgements
License and Copyright

Instead of using the libctl/Scheme interface described in the Meep tutorial, Meep is also callable as a C++ library by writing a C++ program that links to it. The C++ interface provides the greatest level of flexibility in setting up FDTD simulations. There are 19 examples of such C++ programs in the tests/ directory of the Meep codebase which cover a wide range of Meep's functionality. These are a good reference. Finally, we should also note that, while Meep is nominally in C++, it is perhaps better described as "C+". That is, most of the coding style is C-like with a few C++ features.

Contents

Differences from libctl

The C++ interface has several differences from the libctl interface besides the obvious difference in syntax.

The most notable difference is that, while the libctl interface puts the origin (0,0,0) at the center of the computational cell, the C++ interface by default puts the origin at the corner of the computational cell. That is, an L\times L\times L cell goes from ( − L / 2, − L / 2, − L / 2) to (L / 2,L / 2,L / 2) in libctl, but from (0,0,0) to (L,L,L) in C++. (You can change this by calling grid_volume::shift_origin.)

Overview

We begin with a very brief outline of a Meep C++ program, with minimal explanations, leaving more details for the examples below and the C++ reference page. Your C++ program should begin with:

#include <meep.hpp>
using namespace meep;

to include the definitions of the Meep routines. Later, when you compile (see below), you must also link to the Meep libraries.

Your main program should then initialize Meep, and will generally then define a computational grid_volume and the associated structure (describing the geometry and materials), initialize the fields, add sources, and then time-step. In short:

int main(int argc, char **argv) {
  initialize mpi(argc, argv); // do this even for non-MPI Meep
  double resolution = 20; // pixels per distance
  grid_volume v = vol2d(5,10, resolution); // 5x10 2d cell
  structure s(v, eps, pml(1.0));
  fields f(&s);
  
  f.output_hdf5(Dielectric, v.surroundings());
  
  double freq = 0.3, fwidth = 0.1;
  gaussian_src_time src(freq, fwidth);
  f.add_point_source(Ey, src, vec(1.1, 2.3));
  while (f.time() < f.last_source_time()) {
    f.step();
  }
  
  f.output_hdf5(Hz, v.surroundings());
  
  return 0;
} 

This example, doesn't do much — it just runs a Gaussian source and outputs the \mathbf{H}_z field at the end. The dielectric structure is determined by the user-defined function eps, which has the form:

double eps(const vec &p) {
  if (p.x() < 2 && p.y() < 3)
    return 12.0;
  return 1.0;
}

which returns the dielectric function \varepsilon(\mathbf{x}) (here, just a 2×3 rectangle of ε=12 in the upper-left corner). Unlike in the libctl front-end, by default the origin 0 of the coordinate system is at the corner of the computational cell.

Now that you have the basic flavor, we can proceed to some more specific examples.

Computing the Quality Factor of a Resonator

In this first tutorial, we will write the script to compute the quality factor of a 1D Fabry-Perot cavity. For a 1D system, Meep considers a computational cell along the z coordinate. The control file will be a C++ file (having extension *.cpp). In order to use all the classes and subroutines available in Meep, the first two lines of any control file must be the following:

#include <meep.hpp>
using namespace meep;

The particular Fabry-Perot cavity we will investigate consists of an air region bounded by two distributed Bragg reflectors (quarter-wave stack) of epsilon 12 and 1. We choose the size of the defect to be twice as large as the air thickness in the quarter-wave stack, so that a defect mode is found near midgap. This structure will include N = 5 periods in the Bragg reflector on either side of the defect. We can use a larger N but the quality factor may then be too large to compute. The parameters are set up as follows:

const double eps1 = 12.0; // epsilon of layer 1
const double eps2 = 1.0; // epsilon of layer 2
const double grating_periodicity = 1.0;
const double d1 = sqrt(eps2) / (sqrt(eps1)+sqrt(eps2)); // quarter wave stack dimensions
const double d2 = sqrt(eps1) / (sqrt(eps1)+sqrt(eps2));
const double half_cavity_width = d2;
const int N = 5;

Meep supports perfectly matching layers (PML) as absorbing boundary conditions. The PML begins at the edge of the computational volume and works inwards. Hence, we specify the size of the computational cell as follows (note that z_center is half the cell length):

const double pml_thickness = 1.0;
const double z_center = half_cavity_width + N*grating_periodicity + pml_thickness;

To specify a dielectric structure in Meep, we define a function that takes as input one parameter, a position vector, and returns the value of the dielectric at that position.

double eps(const vec &p) {
 vec r = p - vec(z_center);
 if (abs(r.z()) < half_cavity_width)
   return 1.0;
 else {
   double dz = abs(r.z()) - half_cavity_width;
   while (dz > grating_periodicity)
     dz -= grating_periodicity;
   return (dz < d1) ? eps1 : eps2;
 }
}

We are now ready to set up the computational cell, excite the sources, time step the fields, and compute the resulting quality factor. Here we set up the main part of the control file incorporating some of Meep's classes and sub routines. We will excite the mode at the midgap of the bandgap, which we expect to have the largest quality factor since it will have the largest exponential decay.

int main(int argc, char **argv){ 
 initialize mpi(argc,argv);
 const double amicron = 10.0;
 const grid_volume vol = vol1d(2*z_center,amicron);
 structure s(vol,eps,pml(pml_thickness));
 fields f(&s);

Note the constructor for the grid_volume class in 1D which takes as parameters the size of the cell and the resolution. The sources in Meep are excitations of the polarization vector in \mathbf{D}=\epsilon\mathbf{E}+\mathbf{P}. The polarization can be any one of the six cartesian or cylindrical fields. There are a variety of sources including dipole and current sources, gaussian pulses and a continuous wave sources. For now, we use a gaussian source centered at the midgap frequency with a narrow width, along with the start time and end time of the source.

 const double w_midgap = (sqrt(eps1)+sqrt(eps2))/(4*sqrt(eps1)*sqrt(eps2));
 gaussian_src_time src(w_midgap, 1.0, 0.0, 5.0);
 f.add_point_source(Ex, src, vol.center());

Here we make use of the built in function, vol.center(), that automatically computes the geometric center of the of the computational volume. To see what the unit cell of the dielectric function looks, output as an HDF5 file, we invoke the built in command:

f.output_hdf5(Dielectric, v.surroundings())

This versatile command can be used to visualize all of the field components, field energy when called at a particular time step. The resulting output file will be eps-000000.00.h5 and to view it, we first need to convert it to the portable network graphics (PNG) format with the h5topng tool. We are now ready to begin the time stepping of the fields, and do so in a simple loop:

while (f.time() < f.last_source_time()) f.step();

To calculate the quality factor, we need to monitor a field beyond the time at which our point source has been turned off. We can do this in Meep by establishing an array where all fields at a particular point can be monitored. We set up the area as well as the output arrays we will need later as such:

int maxbands = 5;
int ttot = int(400/f.dt)+1; 
complex<double> *p = new complex<double>[ttot];
complex<double> *amps = new complex<double>[maxbands];
double *freq_re = new double[maxbands];
double *freq_im = new double[maxbands];

The variable maxbands is the number of mode frequencies we will be looking for while ttot represents the total number of timesteps for which we will monitor the field component at. For this simulation given the large quality factor that we expect, we will time step for a long period of time which in this case is 400 unit seconds. Now we are ready to time step and collect the field information:

 int i=0;
 while (f.time() < f.last_source_time() + 400) {
   f.step();
   p[i++] = f.get_field(Ex,vol.center());
 }

The get_field function does exactly that, obtains the value of the field at a given position within the computational cell using linear interpolation where necessary. Now we have all the information we need and must obtain the frequencies of the modes which we do by invoking a special tool for harmonic inversion, harminv, to extract the real and imaginary frequencies:

 int num =  do_harminv(p, ttot, f.dt, 0.8*w_midgap, 1.2*w_midgap, maxbands, amps, freq_re, freq_im);

The integer returned by harminv is the number of mode frequencies obtained from the input data p. The particular call to harminv included passing the arrays by value and telling harminv to look for frequencies within 20% of the mid gap frequency up to a maximum of 5 bands. At this point, the necessary information to compute the quality has been stored in the freq_re and freq_im arrays and we compute the quality factor using the formula, Q = − ωr / 2ωi. All that is left to do, is to print these results with the quality factor:

master_printf("frequency,amplitude,quality factor\n"); 
for (int i=0; i<num; ++i) 
  master_printf("%0.6g%+0.6gi,%0.5g%+0.5gi,%g\n",freq_re[i],freq_im[i],real(amps[i]),imag(amps[i]),-0.5*freq_re[i]/freq_im[i]);
return 0;
}

That's it, we are done. The output for this program is 3.26087e+06. The large quality factor is to be expected since our system includes a relatively high number of bragg reflectors (N=5) in each direction and we are exciting the mode at mid gap. Suppose we want to see what the resonant mode looks like, say over one period. The following block of code time steps the field for exactly one period while outputting the Ex field for the computational cell at each time step:

double curr_time = f.time();
while (f.time() < curr_time + 1/w_midgap) {    
 f.output_hdf5(Ex, vol.surroundings());
 f.step();
}  

After we convert the hdf5 files to PNG format and superimpose the images on the dielectric background to produce an animated, we obtain the following:

Resonant mode of fabry-perot cavity

Compiling

In order to compile your code and link against the Meep library, you must link in several additional libraries that Meep requires. There are three possible ways to do this:

1. After compiling the Meep package following the instructions elsewhere in this wiki, place myfile.cpp (where myfile.cpp is a MEEP C++ file of your choosing) in the tests/ subdirectory, cd to the same directory, and invoke:

make myfile.dac

Run the resulting executable via:

./myfile.dac

2. Use the w:pkg-config program, which is installed by default on most GNU/Linux systems. Then, to compile a file foo.cpp using Meep, do:

g++ `pkg-config --cflags meep` foo.cpp -o foo `pkg-config --libs meep`

(Naturally, replace g++ with the name of your C++ compiler if you are not using the GNU compilers.)

3. Compile with g++, this time invoking each library separately:

g++ -malign-double foo.cpp -o foo -lmeep -lhdf5 -lz -lgsl -lharminv -llapack -lcblas -latlas -lfftw3 -lm
Personal tools