Meep C-plus-plus Tutorial

From AbInitio

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

← Previous diff
Revision as of 19:54, 23 October 2005 (edit)
Ardavan (Talk | contribs)

Next diff →
Line 57: Line 57:
using namespace meep; using namespace meep;
-The particular cavity we will investigate is a 1D waveguide bounded by two distributed bragg reflectors whose parameters we set up as follows:+The particular cavity we will investigate is an air cavity bounded by two distributed bragg reflectors (quarter-wave stack) whose parameters we set up as follows:
const double half_cavity_width = 0.5*0.68; const double half_cavity_width = 0.5*0.68;
Line 88: Line 88:
} }
-We are now ready to set up the computational area, excite the sources, and compute the resulting quality factor. Here we set up the main part of the control file, as with any C/C++, incorporating classes and sub routines.+We are now ready to set up the computational area, excite the sources, and compute the resulting quality factor. Here we set up the main part of the control file, as with any C/C++, incorporating 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 within the 1D mirrors.
int main(int argc, char **argv) { int main(int argc, char **argv) {
Line 102: Line 102:
f.add_point_source(Ez, freq, width, 0, 5.0, v.center()); f.add_point_source(Ez, freq, width, 0, 5.0, v.center());
-The parameters to the add_point_source function of the fields class are, in order, the frequency, width, start time, peak time and the location of the excitation in the computational cell. Here we make use of the built in function, v.center(), that automatically computes the geometric center of the of the computational volume. We are now ready to begin the time stepping of the fields, and do so in a simple loop:+To see what the unit cell of the dielectric function looks, output has an HDF5 file, we invoke the built in command:
 + 
 + f.output_hdf5(Dielectric, v.surroundings())
 + 
 +This 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. The parameters to the add_point_source function of the fields class are, in order, the frequency, width, start time, peak time and the location of the excitation in the computational cell. Here we make use of the built in function, v.center(), that automatically computes the geometric center of the of the computational volume. 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(); while (f.time() < f.last_source_time()) f.step();
Line 114: Line 118:
} }
-Now that we have the values of the fields at many time steps, we need to compute the quality factor using the formula, Q=\oemga_real/2\pi*omega_imag. To determine the frequency of the resonant mode, we invoke a built in function specifically designed for the purpose of harmonic inversion, harminv. +Now that we have the values of the fields at many time steps, we need to compute the quality factor using the formula, Q=\oemga_real/2\pi*omega_imag. To determine the frequency of the resonant mode, we invoke a built in function specifically designed for the purpose of harmonic inversion, harminv. The arguments to harminv are passed by reference, and so we declare as such:
- +
- do_harminv(p, fmin, fmax, maxbands, &amps, &freqs_real, &freq_imag)+
- +
 + complex<double> *amps, *freqs;
 + int num;
 + p->harminv(Ez, &amps, &freqs, &num, 0.8*freq, 1.2*freq, 5);
 +
 +At this point, we have the complex amplitudes and frequencies of the modes within 20% of our mid gap frequency. All that is left to do, is to print these results with the quality factor:
 +
 + master_printf("amplitude, frequency, quality factor\n");
 + for (int i=0; i<num; ++i)
 + master_printf("%g%+gi, %g%+gi, %g\n",real(amps[i]),imag(amps[i]),real(freqs[i]),imag(freqs[i]),-real(freqs[i])/2*imag(freqs[i]));
 + return 0;
 + }
 +
 +That's it, we are done. The output for this program is. You will note the large quality factor which is to be expected since our system includes a relatively high number of bragg reflectors (N=15) in each direction and we are exciting the mode at mid gap.
[[Category:Meep]] [[Category:Meep]]

Revision as of 19:54, 23 October 2005

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, 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.

(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 tests/ directory, are written.)

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.

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 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, argc); // do this even for non-MPI Meep
  double resolution = 20; // pixels per distance
  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 cavity in 2D cartesian coordinates. 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 cavity we will investigate is an air cavity bounded by two distributed bragg reflectors (quarter-wave stack) whose parameters we set up as follows:

const double half_cavity_width = 0.5*0.68;
const double air_slit_width = 0.38;
const double grating_periodicity = 0.48;
const double half_waveguide_width = 1.0;
const double num_air_slits = 15.0;
const double high_dielectric = 12.0;
const double low_dielectric = 11.5;

Meep supports perfectly matching layers (PML) and 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:

const double pml_thickness = 1.0;
const double x_center = 7.7 + pml_thickness;
const double y_center = 10.5 + 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) {
  const vec r = p - vec(x_center, y_center);
  double dx = fabs(r.x() - half_cavity_width);
  if (dx > 0 && dx < num_air_slits*grating_periodicty) {
    while (dx > grating_periodicity) dx -= grating_periodicity;
    if (dx < air_slit_width) return 1.0; 
  }
  if (fabs(r.y()) < half_waveguide_width) 
    return high_dielectric;
  else
    return low_dielectric;
 }

We are now ready to set up the computational area, excite the sources, and compute the resulting quality factor. Here we set up the main part of the control file, as with any C/C++, incorporating 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 within the 1D mirrors.

int main(int argc, char **argv) {
 double amicron = 20;
 volume vol = vol2d(2*x_center,2*y_center, amicron);
 structure s(vol,eps,pml_thickness);
 fields f(&s);
 const double wavelength = 1.72;
 const double freq = 1.0/wavelength;

The sources in Meep are excitations of the polarization vector in D=\epsilonE+P. The polarization can be any one of the six cartesian or cylindrical fields. There are a variety of sources including dipoles, gaussian pulses and a continuous wave source. For now, we use the default dipole source for the TM mode consisting of the Ez, Hy, Hx fields.

f.add_point_source(Ez, freq, width, 0, 5.0, v.center());

To see what the unit cell of the dielectric function looks, output has an HDF5 file, we invoke the built in command:

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

This 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. The parameters to the add_point_source function of the fields class are, in order, the frequency, width, start time, peak time and the location of the excitation in the computational cell. Here we make use of the built in function, v.center(), that automatically computes the geometric center of the of the computational volume. 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 a monitor point where all fields at a particular point can be monitored. The monitor point is in the form of a linked list that is reset at each time step in the loop, as follows:

monitor_point *p = NULL;
while (f.time() <= f.last_source_time()+100.0) {
  f.step;
  p = f.get_new_point(v.center(), p);
}

Now that we have the values of the fields at many time steps, we need to compute the quality factor using the formula, Q=\oemga_real/2\pi*omega_imag. To determine the frequency of the resonant mode, we invoke a built in function specifically designed for the purpose of harmonic inversion, harminv. The arguments to harminv are passed by reference, and so we declare as such:

complex<double> *amps, *freqs;
int num;   
p->harminv(Ez, &amps, &freqs, &num, 0.8*freq, 1.2*freq, 5);

At this point, we have the complex amplitudes and frequencies of the modes within 20% of our mid gap frequency. All that is left to do, is to print these results with the quality factor:

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

That's it, we are done. The output for this program is. You will note the large quality factor which is to be expected since our system includes a relatively high number of bragg reflectors (N=15) in each direction and we are exciting the mode at mid gap.

Personal tools