Meep C-plus-plus Tutorial

From AbInitio

Revision as of 04:56, 23 October 2005; Stevenj (Talk | contribs)
(diff) ←Older revision | Current revision | Newer revision→ (diff)
Jump to: navigation, search
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 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 a 1D waveguide bounded by two distributed bragg reflectors 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_sluts = 15.0;
const double high_dielectric = 12.0;
const double low_dielectric = 11.5;

Meep supports periodic 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.

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());

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

Personal tools