PYTHON MEEP DOCUMENTATION (FDTD SIMULATION)

Primary author of this documentation : EL : Emmanuel.Lambert@intec.ugent.be

This documentation and software is provided as-is with no warranty or guarantee.

Document history :

* EL-19/20/21-08-2009 : document creation
* EL-24-08-2009 : small improvements & clarifications.
* EL-25/26-08-2009 : sections 7 & 8 were added.
* EL-03-04/09/2009 :
       -class "structure_eps_pml" (removed again in v0.8).
       -port to Meep 1.1.1 (class 'volume' was renamed to 'grid_volume' and class 'geometric_volume' to 'volume'
       -minor changes in the bent waveguide sample, to make it more consistent with the Scheme version
* EL-07-08/09/2009 : sections 3.2, 8.2, 8.3 : defining a material function with inline C/C++
* EL-10/09/2009 : additions for MPI mode (multiprocessor)
* EL-21/10/2009 : amplitude factor callback function
* EL-22/10/2009 : keyword arguments for runUntilFieldsDecayed
* EL-01/12/2009 : alignment with version 0.8 - III
* EL-01/12/2009 : release 1.0 / added info about environment variables for inline C/C++
* EL-02/12/2009 : minor change for complex_time function + section 1 | v3
* EL-15/01/2010 : missing parameter library_dirs in pars 3.4.2 and 8.3
* EL-19/01/2010 : example for exploiting mirror symmetry / added information on how to shift the origin
* EL-20/01/2010 : minor change
* EL-22/01/2010 : added paragraph 6.2 explaining the new function 'runWithHarminv' - II
* EL-19/02/2010 : added paragraph 6.3 about custom callback classes
* EL-07/05/2010 : added paragraph 3.3 about using a numpy matrix to interface the material definition with Meep
* EL-10/05/2010 : some restructering/finetuning for the changes with python-meep version 1.3
* EL-28/02/2011 : version 1.4 - polygon callback

0. Introduction

Meep is a free finite-difference time-domain (FDTD) simulation software package developed at MIT to model electromagnetic systems. The default scripting language for Meep is Scheme.

Python-meep is a wrapper around the C++ core of Meep : it allows the scripting of Meep-simulations with Python!

Python-meep can be downloaded at : http://launchpad.net/python-meep

This documentation applies to version 1.4.1 of Python-Meep.

1. The general structure of a python-meep program

In general terms, a python-meep program can be structured as follows :

  • import the python-meep binding :

    from meep import *

    This will load the library _meep.so and Python-files meep.py and python_meep.py from path /usr/local/lib/python2.6/dist-packages/

    If you are running in MPI mode (multiprocessor, see section 9), then you have to import module meep_mpi instead :

    from meep_mpi import *

  • define a computational grid volume

    See section 2 below which explains usage of the grid_volume class.

  • define the waveguide structure (describing the geometry, PML and materials)

    See section 3 below which explains usage of the structure class.

  • create an object which will hold the calculated fields

    See section 4 below which explains usage of the field class.

  • define the sources

    See section 5 below which explains usage of the add_point_source and add_volume_source functions.

  • run the simulation (iterate over the time-steps)

    See section 6 below.

Section 7 gives details about defining and retrieving fluxes.

Section 9 gives some complete examples.

Section 10 outlines some differences between Scheme-Meep and Python-Meep.

2. Defining the computational grid volume

The following set of ‘factory functions’ is provided, aimed at creating a grid_volume object. The first arguments define the size of the computational volume, the last argument is the computational grid resolution (in pixels per distance unit).
  • volcyl(rsize, zsize, resolution)

    Defines a cyclical computational grid volume.

  • volone(zsize, resolution) alternatively called vol1d(zsize, resolution)

    Defines a 1-dimensional computational grid volume along the Z-axis.

  • voltwo(xsize, ysize, resolution) alternatively called vol2d(xsize, ysize, a)

    Defines a 2-dimensional computational grid volumes along the X- and Y-axes

  • vol3d(xsize, ysize, zsize, resolution)

    Defines a 3-dimensional computational grid volume.

e.g.: v = volone(6, 10) defines a 1-dimensional computational volume of length 6, with 10 pixels per distance unit.

3. Defining the waveguide structure

The waveguide structure is defined using the class structure, of which the constructor has the following arguments :

  • required : the computational grid volume (a reference to an object of type grid_volume, see section 2 above)

  • required : a function defining the dielectric properties of the materials in the computational grid volume (thus describing the actual waveguide structure). For all-air, the predefined function ‘one’ can be used (epsilon = constant value 1). See note 3.1 below for more information about defining your own custom material function.

  • optional : a boundary region: this is a reference to an object of type boundary_region. There are a number of predefined functions that can be used to create such an object :
    • no_pml() describing a conditionless boundary region (no PML)
    • pml(thickness) : decribing a perfectly matching layer (PML) of a certain thickness (double value) on the boundaries in all directions.
    • pml(thickness, direction) : decribing a perfectly matching layer (PML) of a certain thickness (double value) in a certain direction (X, Y, Z, R or P).
    • pml(thickness, direction, boundary_side) : describing a PML of a certain thickness (double value), in a certain direction (X, Y, Z, R or P) and on the High or Low side. E.g. if boundary_side is Low and direction is X, then a PML layer is added to the −x boundary. The default puts PML layers on both sides of the direction.
  • optional : a function defining a symmetry to exploit, in order to speed up the FDTD calculation (reference to an object of type symmetry). The following predefined functions can be used to create a symmetry object:
    • identity : no symmetry
    • rotate4(direction, grid_volume) : defines a 90° rotational symmetry with ‘direction’ the axis of rotation.
    • rotate2(direction, grid_volume) : defines a 180° rotational symmetry with ‘direction’ the axis of rotation.
    • mirror(direction, grid_volume) : defines a mirror symmetry plane with ‘direction’ normal to the mirror plane.
    • r_to_minus_r_symmetry : defines a mirror symmetry in polar coordinates

    You can apply a phase factor to the symmetry simply by using the * operator, e.g. if symm is a symmetry object, then you can apply a phase factor of -1 as follows : symm = symm * complex(-1,0).

  • optional: the number of chunks in which to split up the calculated geometry. If you leave this empty, it is auto-configured. Otherwise, you would set this to a factor which is a multiple of the number of processors in your MPI run (for multiprocessor configuration).

e.g. : if v is a variable pointing to the computational grid volume, then :
s = structure(v, one) defines a structure with all air (eps=1),
which is equivalent to:
s = structure(v, one, no_pml(), identity(), 1)

Another example : s = structure(v, EPS, pml(0.1,Y) ) with EPS a custom material function, which is explained in the note below.

An example of exploiting symmetry:

#we define a computational volume of 16x32 with resolution 10
vol = voltwo(16,32,10)
#reposition the origin to the center of the computational volume (symmetry is always defined with relative to the origin)
vol.center_origin()
#create a symmetry through the origin with the Y-axis normal to the mirror plane
symm = mirror(Y, vol)
#apply a phase factor to the symmetry
symm = symm * complex(-1,0)
#define the structure and refer to the symmetry with the 4th parameter
s = structure(vol, EPS, pml(1.0), symm)

In order to describe the geometry of the waveguide, we have to provide a ‘material function’ returning the epsilon value as a function of the position.

In python-meep, we have 3 approaches to define such a material function :
  • by defining a class that inherits from class ``Callback``. Through this inheritance, the core meep library will be able to call back to the Python function which describes the material properties. This method is explained in 3.1, but is rather slow and only suitable for simple examples.
  • by defining a Numpy matrix with the epsilon values for every point (using a class that inherits from ``CallbackMatrix``). Python-meep will interface this matrix to the Meep core and Meep can then access this matrix without doing callbacks to Python, which results in great performance. This technique is explained in paragraph 3.3 and is the preferred approach since Python-Meep version 1.3.
  • by using inline C or C++ : this was the fastest technique in early versions of Python-Meep. It creates quite some overhead for the user in terms of required programming skills. The use of a Numpy array offers comparable performance and therefore the use of inline C/C++ is now discouraged.

3.1. Defining a material function with a pure Python callback class [SLOW PERFORMANCE !]

We define a class that inherits from Callback. The method double_vec(self,vec) must be implemented : given a coordinate in the simulation volume (parameter vec), return a real value for the epsilon.

E.g. :

class epsilon(Callback):  #inherit from Callback for integration with the meep core library
        def __init__(self):
                Callback.__init__(self)
        def double_vec(self,vec): #override of function in the Callback class to set the eps function
                #return the epsilon value for a certain point (indicated by the vector v)
                v = vec
                r = v.x()*v.x() + v.y()*v.y()
                dr = sqrt(r)
                while dr>1:
                    dr-=1
                if dr > 0.7001:
                    return 12.0
                return 1.0

Please note the brackets when referring to the x- and y-components of the vector vec. These are crucial : no error message will be thrown if you refer to it as vec.x or vec.y, but the value will always be zero. So, one should write : vec.x() and vec.y().

The meep-python library has a ‘global’ variable EPS, which is used as a reference for communication between the Meep core library and the Python code. We create an instance of our callback-class en assign it as follows to the global EPS variable :

set_EPS_Callback(epsilon().__disown__())
s = structure(v, EPS, no_pml(), identity())

The call to function __disown__() is for memory management purposes and is absolutely required. At the end of our program, we should call : del_EPS_Callback() in order to clean up the global variable.

As described above, this technique is only useful for simple examples, as it has the slowest performance.

3.2 Eps-averaging

EPS-averaging (anisotrpic averaging) is disabled by default, making this behaviour consistent with the behaviour of the Meep C++ core.

You can enable EPS-averaging using the function use_averaging :

#enable EPS-averaging
use_averaging(True)
...
#disable EPS-averaging
use_averaging(False)
...

Enabling EPS-averaging results in slower performance, but more accurate results.

3.3. Using a NUMPY matrix to interface the material definition with Meep

The approach described in 3.1 lets the Meep core library call back to Python for every query of the epsilon-function : this creates quite some overhead, which results in slow performance. An approach which has a lot better performance is to create a Numpy matrix with the epsilon values. Python-Meep then interfaces this matrix in one shot to the Meep C++ core and no more callback to Python will be required afterwards. However, there are some important remarks to make regarding this approach, as Meep will start from a geometry which is inherently already a staircase approach of the geometry that you intended. Therefore, it will be somehow less accurate. An alternative is to use polygons to describe the geometry, see paragraph 3.4. The polygons are an analytically correct representation of the geometry and Meep will thus start from a representation which does not make any simplication.

Example :

class epsilon(CallbackMatrix):
    def __init__(self, pIsWgBent):
        CallbackMatrix.__init__(self)
        master_printf("Creating the material matrix....\n")
        self.meep_resolution = int(resolution)
        eps_matrix = numpy.zeros([10.0 * resolution, 10.0 * resolution], dtype = float)
        len_x = eps_matrix.shape[0]
        len_y = eps_matrix.shape[1]
        for x in range(0, len_x):
            for y in range(0, len_y):
                if ((y_co >= 4.0) and (y_co <= 5.0)):
                     eps_matrix[x,y] = 12.0;
                else:
                     eps_matrix[x,y] = 1.0;
        master_printf("Setting the material matrix...\n")
        self.setMatrix(eps_matrix)
        self.stored_eps_matrix = eps_matrix #to prevent the garbage collector from cleaning up the matrix...
        master_printf("MeepMaterial object initialized.\n")

The following steps are to be followed :

  • Create a class which intherits from CallbackMatrix
  • Set the attribute self.resolution to the resolution that is used for the simulation.
  • Create a numpy array with all zeros and the correct dimensions (as many elements as points in your simulation grid) : eps_matrix = numpy.zeros([size_X * resolution, size_Y * resolution], dtype = float).
  • Fill this array with the correct values for epsilon.
  • Call the function self.setMatrix with a pointer to the matrix as parameter (this gives the Meep C++ core a reference to the matrix).
  • Now additionally store a pointer to the matrix in a local attribute of your class (this is crucial and prevents the garbage collector from cleaning up the matrix when the __init__ function has finished).

Now, similarly as the approach described in 3.1, create an instance of your class and pass it to the function set_EPS_Callback (EPS is a ‘global’ variable EPS, which is used as a reference for communication between the Meep core library and the Python code).

set_EPS_Callback(epsilon().__disown__())
s = structure(v, EPS, no_pml(), identity())

The call to function __disown__() is for memory management purposes and is absolutely required. At the end of our program, we should call : del_EPS_Callback() in order to clean up the global variable.

3.4. Using POLYGONS to interface the material definition with Meep [PREFERRED TECHNIQUE]

In the previous paragraph, we showed how the geometry could be defined using a Numpy matrix. This has the drawback that Meep will start from a representation that is inherently already a staircase approach of the geometry that you intended. Therefore, it will be somehow less accurate.

In this paragraph, we show how you can define the geometry using polygons. The polygons are an analytically correct representation of the geometry and Meep will thus start from a representation which does not make any simplication.

  • the callback class should inherit from PolygonCallback2D (for 2D-simulation) or PolygonCallback3D (for 3D simulations).

  • using function add_polygon, we specify a list of polygons, each associated with a certain epsilon value (for 2D) or a material stack (for 3D). A polygon defined by N point is represented as a numpy array of dimension N x 2. The function add_polygon returns a reference to the underlying C++ polygon object. Every polygon can have a number of inner polygons (the points inside which are not part of the parent polygon): this allows to model for example a donut shape. Inner polygons can be added by calling the function add_inner_polygon.

  • when Meep does a callback for a certain (x,y) coordinate in the grid, a winding number algorithm is used to check inside which polygon the (x,y) coordinate falls :
    • for 2D, the epsilon value of that polygon is taken
    • for 3D, it is checked in which material layer of the stack the (x,y,z) coordinate falls : the epsilon of that material is then taken
  • if a certain point does not match a polygon, then it is assumed that epsilon = 1.0 (air) [possible improvement : make this parametrizable]

Example for 2D : specification of a Silicon waveguide in a surrounding of air

class epsilon(PolygonCallback2D):
    def __init__(self, pIsWgBent):
        PolygonCallback2D.__init__(self)
        master_printf("Callback function for epsilon contructed. Now defining the polygons...\n")
        #Points outside this polygon are assumed to have epsilon 1.0
        pol_points = numpy.zeros((5,2))
        pol_points[0] = [0.0, 10.0]
        pol_points[1] = [30.0, 10.0]
        pol_points[2] = [30.0, 12.0]
        pol_points[3] = [0.0, 12.0]
        pol_points[4] = [0.0, 10.0]
        self.add_polygon(pol_points, 12.0)
        master_printf("Polygons OK.\n")

Example for 3D : specification of a Si waveguide on top of a SiOx cladding

from meep import *  # make it 'meep_mpi' for MPI-meep and 'meep' for non-MPI meep

grid_size_x = 16.0
grid_size_y = 32.0
grid_size_z = 0.5+0.38

class epsilon(PolygonCallback3D):
    def __init__(self):
        PolygonCallback3D.__init__(self)
        master_printf("Callback function for epsilon contructed.\n")
        master_printf("Now defining the matrial stacks...\n")
        ms = numpy.zeros([4,4])
        #first material stack : 0.5 micron of SiOx (eps = 2.31), with 0.38 micron of air (eps = 1.0) on top
        ms[0] = [1, 0.5, 2.31, 2] #material stack ID, height of material, eps of material, number of lines in this stack
        ms[1] = [1, 0.38, 1.0, 2]
        #second material stack : 0.5 micron of SiOx (eps = 2.31), with 0.38 micron of Si (eps = 12.0) on top
        ms[2] = [2, 0.5, 2.31, 2]
        ms[3] = [2, 0.38, 12.0, 2]
        self.add_material_stacks_from_numpy_matrix(ms, 2)
        master_printf("Now defining the polygons...\n")
        wg_length = 16.0
        pad_size = 4.0
        wg_width = 1.5
        grid_size_y = 32.0
        #polygon 1 : area with material stack 1 (SiOx and air)
        pol1_points = numpy.zeros((5,2))
        pol1_points[0] = [0.0,0.0]
        pol1_points[1] = [wg_length, 0.0]
        pol1_points[2] = [wg_length, pad_size]
        pol1_points[3] = [0.0, pad_size]
        pol1_points[4] = [0.0, 0.0]
        self.add_polygon(pol1_points, 1)
        #polygon 2 : area with material stack 2 (SiOx and Si)
        pol2_points = numpy.zeros((5,2))
        pol2_points[0] = [0.0, pad_size]
        pol2_points[1] = [wg_length, pad_size]
        pol2_points[2] = [wg_length, pad_size+wg_width]
        pol2_points[3] = [0.0, pad_size+wg_width]
        pol2_points[4] = [0.0, pad_size]
        self.add_polygon(pol2_points, 2)
        #polygon 3 : area with material stack 1 (SiOx and air)
        pol3_points = numpy.zeros((5,2))
        pol3_points[0] = [0.0,pad_size+wg_width]
        pol3_points[1] = [wg_length, pad_size+wg_width]
        pol3_points[2] = [wg_length, grid_size_y]
        pol3_points[3] = [0.0, grid_size_y]
        pol3_points[4] = [0.0,pad_size+wg_width]
        self.add_polygon(pol3_points, 1)
        master_printf("Polygons OK....\n")

comp_vol = vol3d(grid_size_x,grid_size_y,grid_size_z,10.0)
material = epsilon()
set_EPS_Callback(material.__disown__())
use_averaging(False)
s = structure(comp_vol, EPS, pml(0.25))
f = fields(s)
eps_file =  prepareHDF5File("3d_waveguide.h5")
f.output_hdf5(Dielectric, comp_vol.surroundings(), eps_file)

3.5. Defining a material function with inline C/C++ [DEPRECATED, NO LONGER TO BE USED]

THIS FUNCTIONALITY IS DEPRECATED SINCE VERSION 1.3 : ONE CAN NOW CREATE A NUMPY MATRIX WITH THE EPSILON-VALUES OR USE POLYGON DESCRIPTIONS, SEE PARAGRAPHS 3.3 AND 3.4 FOR DETAILS.

The approach described in 3.1 lets the Meep core library call back to Python for every query of the epsilon-function. This creates a lot of overhead. An approach which has better performance is to define this epsilon-function with an inline C-function or C++ class.

  • If our epsilon-function needs no other parameters than the position vector (X, Y, Z), then we can suffice with an inline C-function (the geometry dependencies are then typically hardcoded).
  • If our epsilon-function needs to base it’s calculation on a more complex set of parameters (e.g. parameters depending on the geometry), then we have to write a C++ class.

For example, in the bent-waveguide example (section 8.3), we can define a generic C++ class which can return the epsilon-value for both the “bend” and “no bend” case, with variable size parameters. We can also take a simpler approach (section 8.2) and write a function in which the geometry size parameters are hardcoded : we then need 2 inline C-functions : one for the “bend” case and one for the “no bend” case.

Make sure the following environment variables are defined :
  • MEEP_INCLUDE : should point to the path containing meep.hpp (and a subdirectory ‘meep’ with vec.hpp and mympi.hpp), e.g.: /usr/include
  • MEEP_LIB : should point to the path containing libmeep.so, e.g. : /usr/lib
  • PYTHONMEEP_INCLUDE : should point to the path containing custom.hpp, e.g. : /usr/local/lib/python2.6/dist-packages

3.5.1 Inline C-function

[ DEPRECATED SINCE PYTHON-MEEP VERSION 1.3 - USE TECHNIQUE FROM PAR. 3.3 OR 3.4 INSTEAD]

First we create a header file, e.g. “eps_function.hpp” which contains our EPS-function. Not that the geometry dependencies are hardcoded (upperLimitHorizontalWg = 4 and lowerLimitHorizontalWg = 5).

namespace meep
{
        static double myEps(const vec &v) {
                double xCo = v.x();
                double yCo = v.y();
                double upperLimitHorizontalWg = 4;
                double lowerLimitHorizontalWg = 5;
                if ((yCo < upperLimitHorizontalWg) || (yCo > lowerLimitHorizontalWg)){
                        return 1.0;
                }
                else return 12.0;
        }
}

Then, in the Python program, we prepare and set the callback function as shown below. prepareCallbackCfunction returns a pointer to the C-function, which we deliver to the Meep core using set_EPS_CallbackInlineFunction.

def initEPS(isWgBent):
        funPtr = prepareCallbackCfunction("myEps","eps_function.hpp")  #name of your function / name of header file
        set_EPS_CallbackInlineFunction(funPtr)
        print "EPS function successfully set."
        return funPtr

We refer to section 8.2 below for a full example.

3.5.2 Inline C++-class

[ DEPRECATED SINCE PYTHON-MEEP VERSION 1.3 - USE TECHNIQUE FROM PAR. 3.3 OR 3.4 INSTEAD]

A more complex approach is to have a C++ object that can accept more parameters when it is constructed. For example this is the case if want to change the parameters of the geometry from inside Python without touching the C++ code.

We create a header file “eps_class.hpp” which contains the definition of the class (the class must inherit from Callback). In the example below, the parameters upperLimitHorizontalWg and widthWg will be communicated from Python upon construction of the object. If these parameters then change (depending on the geometry), the C++ object will follow automatically.

using namespace meep;

namespace meep
{

class myEpsCallBack : virtual public Callback {

    public:
           myEpsCallBack() : Callback() { };
           ~myEpsCallBack() { cout << "Callback object destructed." << endl; };

           myEpsCallBack(double upperLimitHorizontalWg, double widthWg)  : Callback() {
                _upperLimitHorizontalWg = upperLimitHorizontalWg;
                _widthWg = widthWg;
           };

           double double_vec(const vec &v) { //return the EPS-value, depending on the position vector
                double eps = myEps(v, _upperLimitHorizontalWg, _widthWg);
                return eps;
           };

           complex<double> complex_vec(const vec &x) { return 0; }; //no need to implement
           complex<double> complex_time(const double &t) { return 0; };  //ATTENTION : in python-meep version 0.8, this is complex_time(double t)


     private:
           double _upperLimitHorizontalWg;
           double _widthWg;

           double myEps(const vec &v, double upperLimitHorizontalWg, double widthWg) {
                double xCo = v.x();
                double yCo = v.y();
                if ((yCo < upperLimitHorizontalWg) || (yCo > upperLimitHorizontalWg+widthWg)){
                        return 1.0;
                    }
                }
                return 12.0;
           }

};

}

The syntax in Python is a little bit more complex in this case.

We will need to import the module scipy.weave :

from scipy.weave import *

(this is not required for the previous case of a simple inline function)

First we create a list with the names of the parameters that we want to pass to the C++ class. These variables must be declared in the scope where the inline function call happens (see below).

c_params = ['upperLimitHorizontalWg','widthWg']

Then, we prepare the code snippet, using the function prepareCallbackCObjectCode and passing the class name and parameter names list.

c_code = prepareCallbackCObjectCode("myEpsCallBack", c_params)

Finally, we call the inline function, passing :
  • the code snippet
  • the list of parameter names
  • the inline libraries, include directories and headers (helper functions are provided for this, see below). The call to getInlineHeaders should receive the name of the header file (with the definition of the C++ class) as an argument.

funPtr = inline(c_code,c_params, libraries=getInlineLibraries(), include_dirs = getInlineInclude(), library_dirs=getLibraryDirs(), headers = getInlineHeaders("eps_class.hpp") )

def initEPS():
        #the set of parameters that we want to pass to the Callback object upon construction
        #all these variables must be declared in the scope where the "inline" function call happens
        c_params = ['upperLimitHorizontalWg','widthWg']
        #the C-code snippet for constructing the Callback object
        c_code = prepareCallbackCObjectCode("myEpsCallBack", c_params)
        #do the actual inline C-call and fetch the pointer to the Callback object
        funPtr = inline(c_code,c_params, libraries=getInlineLibraries(), include_dirs = getInlineInclude(), library_dirs=getLibraryDirs(), headers = getInlineHeaders("eps_class.hpp") )
        #set the pointer to the callback object in the Python-meep core
        set_EPS_CallbackInlineObject(funPtr)
        print "EPS function successfully set."
        return

We refer to section 8.3 below for a full example.

4. Defining the initial field

This is optional.

We create an object of type fields, which will contain the calculated field.

We must first create a Python class that inherits from class Callback and that will define the function for initialization of the field. Inheritance from class Callback is required, because the core meep library (written in C++) will have to call back to the Python function. For example, let’s call our initialization class fi.

class fi(Callback):  #inherit from Callback for integration with the meep core library
    def __init__(self):
        Callback.__init__(self)
    def complex_vec(self,v): #override of function in the Callback class to set the field initialization function
         #return the field value for a certain point (indicated by the vector v)
        return complex(1.0,0)

The meep-python library has a ‘global’ variable INIF, that is used to bind the meep core library to our Python field initialization class. To set INIF, we have to use the following statement :

set_INIF_Callback(fi().__disown__())  #link the INIF variable to the fi class

We refer to section 3-note1 for more information about the function __disown__().

E.g.: If s is a variable pointing to the structure and comp denotes the component which we are initializing, then the complete field initialization code looks as follows :

f = fields(s)
comp = Hy
f.initialize_field(comp, INIF)

*Important remark* : at the end of our program, we should then call : del_INIF_Callback() in order to clean up the global variable.

The call to initialize_field is not mandatory. If the initial conditions are zeros for all components, then we can rely on the automatic initialization at creation of the object.

We can additionally define Bloch-periodic boundary conditions over the field. This is done with the use_bloch function of the field class, e.g. :

f.use_bloch(vec(0.0))

to be further elaborated - what is the exact meaning of the vector argument? (not well understood at this time)

5. Defining the sources

The definition of the current sources can be done through 2 functions of the fields class :
  • add_point_source(component, src_time, vec, complex)
  • add_volume_source(component, src_time, volume, complex)

Each require as arguments an electromagnetic component (e.g. Ex, Ey, ... and Hx, Hy, ...) and an object of type src_time, which specifies the time dependence of the source (see below).

For a point source, we must specify the center point of the current source using a vector (object of type vec).

For a volume source, we must specify an object of type volume (to be elablorated).

The last argument is an overall complex amplitude number, multiplying the current source (default 1.0).

The following variants are available :
  • add_point_source(component, double, double, double, double, vec centerpoint, complex amplitude, int is_continuous)
    • This is a shortcut function so that no src_time object must be created. This function is preferably used for point sources.
    • The four real arguments define the central frequency, spectral width, start time and cutoff.
    • Set is_continuous to 1 for a continuous point source, to 0 for a Gaussian point source.
  • add_volume_source(component, src_time, volume)

  • add_volume_source(component, src_time, volume, AMPL)
    • AMPL is a built-in reference to a callback function. Such a callback function returns a factor to multiply the source amplitude with (complex value). It receives 1 parameter, i.e. a vector indicating a position RELATIVE to the CENTER of the source. See the example below.
Three classes, inheriting from src_time, are predefined and can be used off the shelf :
  • gaussian_src_time for a Gaussian-pulse source. The constructor demands 2 arguments of type double :
    • the center frequency ω, in units of 2πc
    • the frequency width w used in the Gaussian
  • continuous_src_time for a continuous-wave source proportional to exp(−iωt). The constructor demands 4 arguments :
    • the frequency ω, in units 2πc/distance (complex number)
    • the temporal width of smoothing (default 0)
    • the start time (default 0)
    • the end time (default infinity = never turn off)
  • custom_src_time for a user-specified source function f(t) with start/end times. The constructor demands 4 arguments :
    • The function f(t) specifying the time-dependence of the source
    • ...(2nd argument unclear, to be further elaborated)...
    • the start time of the source (default -infinity)
    • the end time of the source (default +infinity)

For example, in order to define a continuous line source of length 1, from point (6,3) to point (6,4) in 2-dimensional geometry :

#define a continuous source
srcFreq = 0.125
srcWidth = 20
src = continuous_src_time(srcFreq, srcWidth, 0, infinity)
srcComp = Ez
#make it a line source of size 1 starting on position(6,3)
srcGeo = volume(vec(6,3),vec(6,4))
f.add_volume_source(srcComp, src, srcGeo)

Here is an example of the implementation of a callback function for the amplitude factor :

class amplitudeFactor(Callback):
    def __init__(self):
        Callback.__init__(self)
        master_printf("Callback function for amplitude factor activated.\n")

    def complex_vec(self,vec):
        #BEWARE, these are coordinates RELATIVE to the source center !!!!
        x = vec.x()
        y = vec.y()
        master_printf("Fetching amplitude factor for x=%f - y=%f\n" %(x,y) )
        result = complex(1.0,0)
        return result

...
        #define a continuous source
        srcFreq = 0.125
        srcWidth = 20
        src = continuous_src_time(srcFreq, srcWidth, 0, infinity)
        srcComp = Ez
        #make it a line source of size 1 starting on position(6,3)
        srcGeo = volume(vec(6,3),vec(6,4))
        #create callback object for amplitude factor
        af = amplitudeFactor()
        set_AMPL_Callback(af.__disown__())
        f.add_volume_source(pSrcComp, srcGaussian, srcGeo, AMPL)

6. Running the simulation, retrieving field values and exporting HDF5 files

6.1 General case and use of the function ‘runUntilFieldsDecayed

We can now time-step and retrieve various field values along the way. The actual time step value can be retrieved or set through the variable f.dt.

The default time step in Meep is : Courant factor / resolution (in FDTD, the Courant factor relates the time step size to the spatial discretization: cΔt = SΔx. Default for S is 0.5). If no further parametrization is done, then this default value is used.

To trigger a step in time, you call the function f.step().

To step until the source has fully decayed :

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

The function runUntilFieldsDecayed mimicks the behaviour of ‘stop-when-fields-decayed’ in Meep-Scheme. This will run time steps until the source has decayed to 0.001 of the peak amplitude. After that, by default an additional 50 time steps will be run. The function has 7 arguments, of which 4 are mandatory and 3 are optional keywords arguments :

  • 4 regular arguments : reference to the field, reference to the computational grid volume, the source component, the monitor point.
  • keyword argument pHDF5OutputFile : reference to a HDF5 file (constructed with the function prepareHDF5File); default : None (no ouput to files)
  • keyword argument pH5OutputIntervalSteps : step interval for output to HDF5 (default : 50)
  • keyword argument pDecayedStopAfterSteps : the number of steps to continue after the source has decayed to 0.001 of the peak amplitude at the probing point (default: 50)

We further refer to section 8 below where this function is applied in an example.

A rich functionality is available for retrieving field information. Some examples :

  • f.energy_in_box(v.surroundings())
  • f.electric_energy_in_box(v.surroundings())
  • f.magnetic_energy_in_box(v.surroundings())
  • f.thermo_energy_in_box(v.surroundings())
  • f.total_energy()
  • f.field_energy_in_box(v.surroundings())
  • f.field_energy_in_box(component, v.surroundings()) where the first argument is the electromagnetic component (Ex, Ey, Ez, Er, Ep, Hx, Hy, Hz, Hr, Hp, Dx, Dy, Dz, Dp, Dr, Bx, By, Bz, Bp, Br, Dielectric or Permeability)
  • f.flux_in_box(X, v.surroundings()) where the first argument is the direction (X, Y, Z, R or P)

We can probe the field at certain points by defining a monitor point as follows :

m = monitor_point()
p = vec(2.10) #vector identifying the point that we want to probe
f.get_point(m, p)
m.get_component(Hx)

We can export the dielectric function and e.g. the Ex component of the field to HDF5 files as follows :

#make sure you start your Python session with 'sudo' or write rights to the current path
feps = prepareHDF5File("eps.h5")
f.output_hdf5(Dielectric, v.surroundings(), feps)   #export the Dielectric structure so that we can visually verify it
fex = prepareHDF5File("ex.h5")
while (f.time() < f.last_source_time()):
     f.step()
     f.output_hdf5(Ex, v.surroundings(), fex, 1) #export the Ex component, appending to the file "ex.h5"

If you want to output only part of the field to HDF5, then check the instructions in paragraph 6.3 : it explains how you can take a custom action in every step of the simulation.

6.2 Running a simulation for Harminv analysis (use of the function ‘runWithHarminv‘)

The function ‘runWithHarminv’ is provided for doing Harminv analysis and outputting the fields to HDF5 : .

runWithHarminv(pField, pVol, pComp, pProbingPoint, pFreq, pDf, pMaxbands, pTimeAfterSources, pTimePhase3, pHDF5OutputFile, pH5OutputIntervalSteps)

The execution strategy is as follows :
  • PHASE 1 :
    • run the field pField until all sources have extinct.
    • then run for an extra pTimeAfterSources (default : 300) while collecting data for Harminv analysis at the probing point pProbingPoint (for field component pComp).
    • a custom callback class may be specified : during every step in phase 1, this class will then be called (see paragraph 6.3)
  • PHASE 2 : do Harminv analysis in the frequency band from pFreq - pDf /2 to pFreq - pDf /2 and search for manimum pMaxbands.

  • PHASE 3 : run for an additional pTimePhase3 (default : 1 / pFreq (1 period)) and send HDF5 output every pH5OutputIntervalSteps (default 5) to the HDF5 file specified by pHDF5OutputFile. By default HDF5 files are only generated during phase 3 (every 5 steps, this is parameter pH5OutputIntervalSteps). If you want to change this behaviour, then you can specify a custom callback class : during every step in phase 3, this class will then be called and you can decide wether to output HDF5 or take another specific action. See paragraph 6.3 for detailed instructions on this.

The return value of runWithHarminv is an array containing vectors with the results of the Harminv : frequency, amplitude, quality factor.

Below is an example of a ring in which we excite a Gaussian source. After sources are extinct, we run for an additional 300 time unit and do Harminv analysis. After this phase, we run for another period and generate HDF5-output (which will allow to visualize the resonances). The source code of the example can be found in the /samples directory of the Python-Meep distribution.

from meep import *  # make it 'meep_mpi' for MPI-meep and 'meep' for non-MPI meep

from math import *
import numpy
import sys
import glob

res = 10.0

width = 1.0
innerRadius = 1.0
outerRadius = innerRadius+width
pad = 4.0
dpml = 2.0
gridSizeX = 2 * (innerRadius + width + pad + dpml)
gridSizeY = gridSizeX
centerX = gridSizeX / 2.0
centerY = gridSizeY / 2.0

srcFreqCenter = 0.15 #gaussian source center frequency
srcPulseWidth = 0.10 #gaussian source pulse width
srcComp = Ez #gaussian source component

#this function plots the waveguide material as a function of a vector(X,Y)
class epsilon(Callback):
    def __init__(self):
        Callback.__init__(self)
        master_printf("Callback function for epsilon activated.\n")

    def double_vec(self,vec):
        x = vec.x() - centerX
        y = vec.y() - centerY
        r = sqrt(x*x+y*y)
        if (r>innerRadius) and (r<=outerRadius):
           return 11.56
        else:
           return 1.0

def runSimul():
        master_printf("Starting...\n")
        #create the computational grid
        master_printf("Size of computational volume: %f x %f with resolution %f\n" %(gridSizeX,gridSizeY,res))
        ringWgVol = voltwo(gridSizeX,gridSizeY,res)
        master_printf("Center at : %f - %f\n" %(centerX,centerY))
        #create the field : we create a structure with PML (thickness 'dpml') on all boundaries, in all directions, using the material function 'EPS'
        material = epsilon()
        set_EPS_Callback(material.__disown__())
        s = structure(ringWgVol, EPS, pml(dpml))
        ringField = fields(s)

        #define a gaussian point source
        srcGaussian = gaussian_src_time(srcFreqCenter, srcPulseWidth )
        srcXco = centerX + innerRadius + 0.1
        srcYco = centerY
        ringField.add_point_source(srcComp, srcGaussian, vec(srcXco,srcYco))
        master_printf("Field created...\n")
        filenameEps = "./harminv_Eps.h5"
        filenameComp = "./harminv_Comp.h5"

        #export the dielectric structure (so that we can visually verify the waveguide structure)
        ringDielectricFile =  prepareHDF5File(filenameEps)
        ringField.output_hdf5(Dielectric, ringWgVol.surroundings(), ringDielectricFile)

        #create the file for the field components
        ringWgFileOutputComp = prepareHDF5File(filenameComp)

        master_printf("Calculating...\n")
        ringWgProbingPoint =  vec(srcXco,srcYco)
        #Now run the fields, do Harminv analysis, followed by HDF5 output during 1 period
        results = runWithHarminv(ringField, ringWgVol, srcComp, ringWgProbingPoint, srcFreqCenter, srcPulseWidth, 100, pOutputHDF5OnlyAfterAnalysis = False, pAdditionalTimeAfterHarminv = (1.0 / srcFreqCenter), pTimeAfterSources = 300, pH5OutputIntervalSteps=6, pHDF5OutputFile = ringWgFileOutputComp)
        master_printf("Done..!\n")

        del_EPS_Callback()

        print results


runSimul()

If you want to output only part of the field to HDF5, then check the instructions in paragraph 6.3. That section also explains how you can take a custom action in every step of the simulation.

6.3 Using your own callback class in every step of the simulation (how to output only part of the field)

Both runUntilFieldsDecayed and runWithHarminv allow to specify a custom class that will be called in every step of the simulation. This permits to output only specific parts of the field or take another custom action, such as calculating and storing the total energy, etc.

  • In runUntilFieldsDecayed, replace the parameters pHDF5OutputFile and pH5OutputIntervalSteps with the parameter pStepCallback.
  • In runWithHarminv , replace the parameters pHDF5OutputFile and pH5OutputIntervalSteps with the parameter pStepCallbackPhase3 (for custom callback in phase 3) and add parameter pStepCallbackPhase1 for custom callback in every step of phase 1.

The new parameters should point to a callable Python class (i.e. a class implementing the __call__ method).

Below is an example of a callable class that outputs the HDF5 field every x steps and prints the total energy in the volume :

class myOutputHDF5(object):

        """ HDF5 files are written to HDF5 file every 'pH5OutputIntervalSteps' steps, while also printing the total energy """
        def __init__(self, pHDF5OutputFile, pH5OutputIntervalSteps, pField, pComp, pExportVolume):
                self.stepsCount = 0
                #HDF5 file to which to output the field
                self.pHDF5OutputFile = pHDF5OutputFile
                #interval at which to output the field
                self.pH5OutputIntervalSteps = pH5OutputIntervalSteps
                #reference to the field
                self.pField = pField
                #reference to the component that we want to output
                self.pComp = pComp
                #reference to the volume that we want to ouput
                self.pExportVolume = pExportVolume

        def __call__(self):
                #this function will be called in every step of the simulation
                if (self.stepsCount % self.pH5OutputIntervalSteps == 0):
                        self.doOutput()
                self.stepsCount = self.stepsCount + 1

        def doOutput(self):
                master_printf("Total energy in volume at step %i : %f \n" %(self.stepsCount, self.pField.total_energy()))
                self.pField.output_hdf5(self.pComp, pExportVolume, self.pHDF5OutputFile, 1)

Suppose now that you do not want to output the field of the complete computational volume, but only of a subvolume. Let’s apply this to the example above : the parameter pExportVolume can be assgined a subvolume of the complete computational volume. You can create such a subvolume by specifying 2 vectors :

  • in 2D : if the 2 points defined by these vector are not aligned on one line, they define the rectangle in which the field will be exported.
  • in 3D : if the 2 points defined by these vector are not falling in the same plane, they will define the cube in which the field will be exported.

In every step of the simulation, the callback class will be triggered : every pH5OutputIntervalSteps steps, the field in the subvolume pExportVolume will be exported to HDF5 file.

point1 = vec(10,20)
point2 = vec(25,40)
subVol = volume(point1, point2)
h5file = prepareHDF5File("./example.h5")
customOutput = myOutputHDF5(pHDF5OutputFile = h5file, pH5OutputIntervalSteps = 10, pField = field, pComp = Ez, pExportVolume = subVol)
runUntilFieldsDecayed(... pStepCallback = customOutput)
The /samples directory of the Python-Meep distribution contains 2 examples with custom callback classes :
  • in combination with runUntilFieldsDecayed, the sample script python_meep_straight_wg_customHDF5output.py
  • in combination with runWithHarminv, the sample script python_meep_harminv_ring_customHDF5output.py

7. Defining and retrieving fluxes

First we define a flux plane. This is done through the creation of an object of type volume (specifying 2 vectors as arguments).

Then we apply this flux plane to the field, specifying 4 parameters :
  • the reference to the volume object
  • the minimum frequency (in the example below, this is srcFreqCenter-(srcPulseWidth/2.0))
  • the maximum frequency (in the example below this is srcFreqCenter+(srcPulseWidth/2.0) )
  • the number of discrete frequencies that we want to monitor in the flux (in the example below, this is 100).

After running the simulation, we can retrieve the flux values through the function getFluxData() : this returns a 2-dimensional array with the frequencies and actual flux values.

E.g., if f is the field, then we proceed as follows :

#define the flux plane and flux parameters
fluxplane = volume(vec(1,2),vec(1,3))
flux = f.add_dft_flux_plane(fluxplane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100)

#run the calculation
while (f.time() < f.last_source_time()):
    f.step()

#retrieve the flux data
fluxdata = getFluxData(flux)
frequencies = fluxdata[0]
fluxvalues = fluxdata[1]

8. The “90 degree bent waveguide example in Python

We have ported the “90 degree bent waveguide” example from the Meep-Scheme tutorial to Python.

The original example can be found on the following URL : http://ab-initio.mit.edu/wiki/index.php/Meep_Tutorial (section ‘A 90° bend’).

You can find the source code below and also in the ``/samples/bent_waveguide`` directory of the Python-Meep distribution.

Sections 8.2 contains the same example, but using a Numpy matrix for interfacing the material with Meep (technique described in paragraph 6.4). This is actually the preferred technique, as it has the best performance and requires only programming in Python.

The flavors with inline C/C++ (descibed in paragraphs 3.2 and 3.4) are no longer shown in this tutorial, as the technique is more or less deprecated (the source code is still available in the /samples/bent_waveguide` directory of the Python-Meep distribution).

The following animated gifs can be produced from the HDF5-files (see the script included in directory ‘samples’) :

_images/bentwgNB.gif _images/bentwgB.gif

And here is the graph of the transmission, reflection and loss fluxes, showing the same results as the example in Scheme:

_images/fluxes.png

8.1 With a pure Python class as EPS-function : SLOW PERFORMANCE !

A bottleneck in this version is the epsilon-function, which is written in Python. This means that the Meep core library must do a callback to the Python function, which creates a lot of overhead. An approach which has a much better performance is to use a Numpy matrix : the Meep core library can then directly access this matrix (in the C-space), without overhead (see next paragraph 8.2).

from meep import *
from math import *
from python_meep import *
import numpy
import matplotlib.pyplot
import sys

res = 10.0
gridSizeX = 16.0
gridSizeY = 32.0
padSize = 4.0
wgLengthX = gridSizeX - padSize
wgLengthY = gridSizeY - padSize
wgWidth = 1.0 #width of the waveguide
wgHorYCen = padSize +  wgWidth/2.0 #horizontal waveguide center Y-pos
wgVerXCen = wgLengthX - wgWidth/2.0 #vertical waveguide center X-pos (in case there is a bend)
srcFreqCenter = 0.15 #gaussian source center frequency
srcPulseWidth = 0.1 #gaussian source pulse width
srcComp = Ez #gaussian source component

#this function plots the waveguide material as a function of a vector(X,Y)
class epsilon(Callback):
    def __init__(self, pIsWgBent):
        Callback.__init__(self)
        self.isWgBent = pIsWgBent
    def double_vec(self,vec):
        if (self.isWgBent): #there is a bend
            if ((vec.x()<wgLengthX) and (vec.y() >= padSize) and (vec.y() <= padSize+wgWidth)):
                return 12.0
            elif ((vec.x()>=wgLengthX-wgWidth) and (vec.x()<=wgLengthX) and vec.y()>= padSize ):
                return 12.0
            else:
                return 1.0
        else: #there is no bend
            if ((vec.y() >= padSize) and (vec.y() <= padSize+wgWidth)):
                return 12.0
            else:
                return 1.0

def createField(pCompVol, pWgLengthX, pWgWidth, pIsWgBent, pSrcFreqCenter, pSrcPulseWidth, pSrcComp):
        #we create a structure with PML of thickness = 1.0 on all boundaries,
        #in all directions,
        #using the material function EPS
        material = epsilon(pIsWgBent)
        set_EPS_Callback(material.__disown__())
        s = structure(pCompVol, EPS, pml(1.0) )
        f = fields(s)
        #define a gaussian line source of length 'wgWidth' at X=wgLength/2, Y=padSize
        srcGaussian = gaussian_src_time(pSrcFreqCenter, pSrcPulseWidth )
        srcGeo = volume(vec(1.0,padSize),vec(1.0,padSize+pWgWidth))
        f.add_volume_source(pSrcComp, srcGaussian, srcGeo, 1)
        print "Field created..."
        return f


#FIRST WE WORK OUT THE CASE WITH NO BEND
#----------------------------------------------------------------
print "*1* Starting the case with no bend..."
#create the computational grid
noBendVol = voltwo(gridSizeX,gridSizeY,res)

#create the field
wgBent = 0 #no bend
noBendField = createField(noBendVol, wgLengthX, wgWidth, wgBent, srcFreqCenter, srcPulseWidth, srcComp)

bendFnEps = "./bentwgNB_Eps.h5"
bendFnEz = "./bentwgNB_Ez.h5"
#export the dielectric structure (so that we can visually verify the waveguide structure)
noBendDielectricFile =  prepareHDF5File(bendFnEps)
noBendField.output_hdf5(Dielectric, noBendVol.surroundings(), noBendDielectricFile)

#create the file for the field components
noBendFileOutputEz = prepareHDF5File(bendFnEz)

#define the flux plane for the reflected flux
noBendReflectedfluxPlaneXPos = 1.5 #the X-coordinate of our reflection flux plane
noBendReflectedFluxPlane = volume(vec(noBendReflectedfluxPlaneXPos,wgHorYCen-wgWidth),vec(noBendReflectedfluxPlaneXPos,wgHorYCen+wgWidth))
noBendReflectedFlux = noBendField.add_dft_flux_plane(noBendReflectedFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100)

#define the flux plane for the transmitted flux
noBendTransmfluxPlaneXPos = gridSizeX - 1.5;  #the X-coordinate of our transmission flux plane
noBendTransmFluxPlane = volume(vec(noBendTransmfluxPlaneXPos,wgHorYCen-wgWidth),vec(noBendTransmfluxPlaneXPos,wgHorYCen+wgWidth))
noBendTransmFlux = noBendField.add_dft_flux_plane(noBendTransmFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100 )

print "Calculating..."
noBendProbingpoint = vec(noBendTransmfluxPlaneXPos,wgHorYCen) #the point at the end of the waveguide that we want to probe to check if source has decayed
runUntilFieldsDecayed(noBendField, noBendVol, srcComp, noBendProbingpoint, noBendFileOutputEz)
print "Done..!"

#construct 2-dimensional array with the flux plane data
#see python_meep.py
noBendReflFlux = getFluxData(noBendReflectedFlux)
noBendTransmFlux = getFluxData(noBendTransmFlux)

#save the reflection flux from the "no bend" case as minus flux in the temporary file 'minusflux.h5'
noBendReflectedFlux.scale_dfts(-1);
f = open("minusflux.h5", 'w') #truncate file if already exists
f.close()
noBendReflectedFlux.save_hdf5(noBendField, "minusflux")

del_EPS_Callback()


#AND SECONDLY FOR THE CASE WITH BEND
#----------------------------------------------------------------
print "*2* Starting the case with bend..."
#create the computational grid
bendVol = voltwo(gridSizeX,gridSizeY,res)

#create the field
wgBent = 1 #there is a bend
bendField = createField(bendVol, wgLengthX, wgWidth, wgBent, srcFreqCenter, srcPulseWidth, srcComp)

#export the dielectric structure (so that we can visually verify the waveguide structure)
bendFnEps = "./bentwgB_Eps.h5"
bendFnEz = "./bentwgB_Ez.h5"
bendDielectricFile = prepareHDF5File(bendFnEps)
bendField.output_hdf5(Dielectric, bendVol.surroundings(), bendDielectricFile)

#create the file for the field components
bendFileOutputEz = prepareHDF5File(bendFnEz)

#define the flux plane for the reflected flux
bendReflectedfluxPlaneXPos = 1.5 #the X-coordinate of our reflection flux plane
bendReflectedFluxPlane = volume(vec(bendReflectedfluxPlaneXPos,wgHorYCen-wgWidth),vec(bendReflectedfluxPlaneXPos,wgHorYCen+wgWidth))
bendReflectedFlux = bendField.add_dft_flux_plane(bendReflectedFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100)

#load the minused reflection flux from the "no bend" case as initalization
bendReflectedFlux.load_hdf5(bendField, "minusflux")


#define the flux plane for the transmitted flux
bendTransmfluxPlaneYPos = padSize + wgLengthY - 1.5; #the Y-coordinate of our transmission flux plane
bendTransmFluxPlane = volume(vec(wgVerXCen - wgWidth,bendTransmfluxPlaneYPos),vec(wgVerXCen + wgWidth,bendTransmfluxPlaneYPos))
bendTransmFlux = bendField.add_dft_flux_plane(bendTransmFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100 )

print "Calculating..."
bendProbingpoint = vec(wgVerXCen,bendTransmfluxPlaneYPos) #the point at the end of the waveguide that we want to probe to check if source has decayed
runUntilFieldsDecayed(bendField, bendVol, srcComp, bendProbingpoint, bendFileOutputEz)
print "Done..!"

#construct 2-dimensional array with the flux plane data
#see python_meep.py
bendReflFlux = getFluxData(bendReflectedFlux)
bendTransmFlux = getFluxData(bendTransmFlux)

del_EPS_Callback()

#SHOW THE RESULTS IN A PLOT
frequencies = bendReflFlux[0] #should be equal to bendTransmFlux.keys() or noBendTransmFlux.keys() or ...
Ptrans = [x / y for x,y in zip(bendTransmFlux[1], noBendTransmFlux[1])]
Prefl = [ abs(x / y) for x,y in zip(bendReflFlux[1], noBendTransmFlux[1]) ]
Ploss = [ 1-x-y for x,y in zip(Ptrans, Prefl)]

matplotlib.pyplot.plot(frequencies, Ptrans, 'bo')
matplotlib.pyplot.plot(frequencies, Prefl, 'ro')
matplotlib.pyplot.plot(frequencies, Ploss, 'go' )

matplotlib.pyplot.show()

8.2 With a Numpy matrix for the material definition

See paragraph 3.3 for a description of this technique.

Please beware of the accuracy restrictions, as explained in 3.3.

'''
Variant of the 90 degree BENT WAVEGUIDE SAMPLE which uses a NUMPY matrix

@author: Emmanuel.Lambert@intec.ugent.be
'''

from meep import *  # make it 'meep_mpi' for MPI-meep and 'meep' for non-MPI meep

from math import *
import numpy
import matplotlib.pyplot
import sys

res = 10.0
gridSizeX = 16.0
gridSizeY = 32.0
padSize = 4.0
wgLengthX = gridSizeX - padSize
wgLengthY = gridSizeY - padSize
wgWidth = 1.0 #width of the waveguide
wgHorYCen = padSize +  wgWidth/2.0 #horizontal waveguide center Y-pos
wgVerXCen = wgLengthX - wgWidth/2.0 #vertical waveguide center X-pos (in case there is a bend)
srcFreqCenter = 0.15 #gaussian source center frequency
srcPulseWidth = 0.1 #gaussian source pulse width
srcComp = Ez #gaussian source component


class epsilon(CallbackMatrix):
    def __init__(self, pIsWgBent):
        CallbackMatrix.__init__(self)
        master_printf("Creating the material matrix....\n")
        self.meep_resolution = int(res)
        eps_matrix = numpy.zeros([gridSizeX * res, gridSizeY * res], dtype = float)
        len_x = eps_matrix.shape[0]
        len_y = eps_matrix.shape[1]
        for x in range(0, len_x):
            for y in range(0, len_y):
                y_co = y / res
                x_co = x / res
                if (pIsWgBent): #there is a bend
                    if ((x_co<wgLengthX) and (y_co >= padSize) and (y_co <= padSize+wgWidth)):
                        eps_matrix[x,y] = 12.0;
                    elif ((x_co>=wgLengthX-wgWidth) and (x_co<=wgLengthX) and y_co>=padSize ):
                        eps_matrix[x,y] = 12.0;
                    else:
                        eps_matrix[x,y] = 1.0;
                else: #there is no bend
                        if ((y_co >= padSize) and (y_co <= padSize+wgWidth)):
                             eps_matrix[x,y] = 12.0;
                        else:
                             eps_matrix[x,y] = 1.0;
        master_printf("Setting the material matrix...\n")
        self.setMatrix(eps_matrix)
        self.stored_eps_matrix = eps_matrix #CRUCIAL !! To prevent the garbage collector from cleaning it up...
        master_printf("MeepMaterial object initialized.\n")



class amplitudeFactor(Callback):
    def __init__(self):
        Callback.__init__(self)
        master_printf("Callback function for amplitude factor activated.\n")

    def complex_vec(self,vec):
        #BEWARE, these are coordinates RELATIVE to the source center !!!!
        x = vec.x()
        y = vec.y()
        master_printf("Fetching amplitude factor for x=%f - y=%f\n" %(x,y) )
        result = complex(1.0,0)
        return result

def createField(pCompVol, pWgLengthX, pWgWidth, pIsWgBent, pSrcFreqCenter, pSrcPulseWidth, pSrcComp):
        #we create a structure with PML of thickness = 1.0 on all boundaries,
        #in all directions,
        #using the material function EPS
        material = epsilon(pIsWgBent)
        set_EPS_Callback(material.__disown__())
        #use_averaging(False)  #--> uncomment this line to disable eps-averaging!
        s = structure(pCompVol, EPS, pml(1.0))
        f = fields(s)
        #define a gaussian line source of length 'wgWidth' at X=wgLength/2, Y=padSize
        srcGaussian = gaussian_src_time(pSrcFreqCenter, pSrcPulseWidth )
        srcGeo = volume(vec(1.0,padSize),vec(1.0,padSize+pWgWidth))
        af = amplitudeFactor()
        set_AMPL_Callback(af.__disown__())
        f.add_volume_source(pSrcComp, srcGaussian, srcGeo, AMPL)
        master_printf("Field created...\n")
        return f


master_printf("** Bent waveguide sample, version 17-11-2009 **\n")

master_printf("Running on %d processor(s)...\n",count_processors())

#FIRST WE WORK OUT THE CASE WITH NO BEND
#----------------------------------------------------------------
master_printf("*1* Starting the case with no bend...\n")
#create the computational grid
noBendVol = voltwo(gridSizeX,gridSizeY,res)

#create the field
wgBent = 0 #no bend
noBendField = createField(noBendVol, wgLengthX, wgWidth, wgBent, srcFreqCenter, srcPulseWidth, srcComp)

bendFnEps = "./bentwgNB_Eps.h5"
bendFnEz = "./bentwgNB_Ez.h5"
#export the dielectric structure (so that we can visually verify the waveguide structure)
noBendDielectricFile =  prepareHDF5File(bendFnEps)
noBendField.output_hdf5(Dielectric, noBendVol.surroundings(), noBendDielectricFile)

#create the file for the field components
noBendFileOutputEz = prepareHDF5File(bendFnEz)

#define the flux plane for the reflected flux
noBendReflectedfluxPlaneXPos = 1.5 #the X-coordinate of our reflection flux plane
noBendReflectedFluxPlane = volume(vec(noBendReflectedfluxPlaneXPos,wgHorYCen-wgWidth),vec(noBendReflectedfluxPlaneXPos,wgHorYCen+wgWidth))
noBendReflectedFlux = noBendField.add_dft_flux_plane(noBendReflectedFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100)

#define the flux plane for the transmitted flux
noBendTransmfluxPlaneXPos = gridSizeX - 1.5;  #the X-coordinate of our transmission flux plane
noBendTransmFluxPlane = volume(vec(noBendTransmfluxPlaneXPos,wgHorYCen-wgWidth),vec(noBendTransmfluxPlaneXPos,wgHorYCen+wgWidth))
noBendTransmFlux = noBendField.add_dft_flux_plane(noBendTransmFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100 )

master_printf("Calculating...")
noBendProbingpoint = vec(noBendTransmfluxPlaneXPos,wgHorYCen) #the point at the end of the waveguide that we want to probe to check if source has decayed
runUntilFieldsDecayed(noBendField, noBendVol, srcComp, noBendProbingpoint, pHDF5OutputFile = noBendFileOutputEz, pH5OutputIntervalSteps=50)
master_printf("Done..!")

#construct 2-dimensional array with the flux plane data
#see python_meep.py
noBendReflFlux = getFluxData(noBendReflectedFlux)
noBendTransmFlux = getFluxData(noBendTransmFlux)

#save the reflection flux from the "no bend" case as minus flux in the temporary file 'minusflux.h5'
noBendReflectedFlux.scale_dfts(-1);
f = open("minusflux.h5", 'w') #truncate file if already exists
f.close()
noBendReflectedFlux.save_hdf5(noBendField, "minusflux")

del_EPS_Callback()


#AND SECONDLY FOR THE CASE WITH BEND
#----------------------------------------------------------------
master_printf("*2* Starting the case with bend...\n")
#create the computational grid
bendVol = voltwo(gridSizeX,gridSizeY,res)

#create the field
wgBent = 1 #there is a bend
bendField = createField(bendVol, wgLengthX, wgWidth, wgBent, srcFreqCenter, srcPulseWidth, srcComp)

#export the dielectric structure (so that we can visually verify the waveguide structure)
bendFnEps = "./bentwgB_Eps.h5"
bendFnEz = "./bentwgB_Ez.h5"
bendDielectricFile = prepareHDF5File(bendFnEps)
bendField.output_hdf5(Dielectric, bendVol.surroundings(), bendDielectricFile)

#create the file for the field components
bendFileOutputEz = prepareHDF5File(bendFnEz)

#define the flux plane for the reflected flux
bendReflectedfluxPlaneXPos = 1.5 #the X-coordinate of our reflection flux plane
bendReflectedFluxPlane = volume(vec(bendReflectedfluxPlaneXPos,wgHorYCen-wgWidth),vec(bendReflectedfluxPlaneXPos,wgHorYCen+wgWidth))
bendReflectedFlux = bendField.add_dft_flux_plane(bendReflectedFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100)

#load the minused reflection flux from the "no bend" case as initalization
bendReflectedFlux.load_hdf5(bendField, "minusflux")


#define the flux plane for the transmitted flux
bendTransmfluxPlaneYPos = padSize + wgLengthY - 1.5; #the Y-coordinate of our transmission flux plane
bendTransmFluxPlane = volume(vec(wgVerXCen - wgWidth,bendTransmfluxPlaneYPos),vec(wgVerXCen + wgWidth,bendTransmfluxPlaneYPos))
bendTransmFlux = bendField.add_dft_flux_plane(bendTransmFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100 )

master_printf("Calculating...\n")
bendProbingpoint = vec(wgVerXCen,bendTransmfluxPlaneYPos) #the point at the end of the waveguide that we want to probe to check if source has decayed
runUntilFieldsDecayed(bendField, bendVol, srcComp, bendProbingpoint, pHDF5OutputFile = noBendFileOutputEz)
master_printf("Done..!")

#construct 2-dimensional array with the flux plane data
#see python_meep.py
bendReflFlux = getFluxData(bendReflectedFlux)
bendTransmFlux = getFluxData(bendTransmFlux)

del_EPS_Callback()

#SHOW THE RESULTS IN A PLOT
frequencies = bendReflFlux[0] #should be equal to bendTransmFlux.keys() or noBendTransmFlux.keys() or ...
Ptrans = [x / y for x,y in zip(bendTransmFlux[1], noBendTransmFlux[1])]
Prefl = [ abs(x / y) for x,y in zip(bendReflFlux[1], noBendTransmFlux[1]) ]
Ploss = [ 1-x-y for x,y in zip(Ptrans, Prefl)]

matplotlib.pyplot.plot(frequencies, Ptrans, 'bo')
matplotlib.pyplot.plot(frequencies, Prefl, 'ro')
matplotlib.pyplot.plot(frequencies, Ploss, 'go' )

matplotlib.pyplot.show()

8.3 With POLYGONS for the material definition [PREFERRED APPROACH]

This approach is very accurate and has at the same time a very good performance. See paragraph 3.4 for a description of this technique.

'''
Conversion into Python of the 90 degree BENT WAVEGUIDE SAMPLE from Meep-Scheme tutorial
(see http://ab-initio.mit.edu/wiki/index.php/Meep_Tutorial )

Using the technique of POLYGONS to define the simulation geometry.

@author: Emmanuel.Lambert@intec.ugent.be
'''

from meep import *  # make it 'meep_mpi' for MPI-meep and 'meep' for non-MPI meep

from math import *
import numpy
import matplotlib.pyplot
import sys

res = 10.0
gridSizeX = 16.0
gridSizeY = 32.0
padSize = 4.0
wgLengthX = gridSizeX - padSize
wgLengthY = gridSizeY - padSize
wgWidth = 1.0 #width of the waveguide
wgHorYCen = padSize +  wgWidth/2.0 #horizontal waveguide center Y-pos
wgVerXCen = wgLengthX - wgWidth/2.0 #vertical waveguide center X-pos (in case there is a bend)
srcFreqCenter = 0.15 #gaussian source center frequency
srcPulseWidth = 0.1 #gaussian source pulse width
srcComp = Ez #gaussian source component


#this function plots the waveguide material as a function of a vector(X,Y)
class epsilon(PolygonCallback2D):
    def __init__(self, pIsWgBent):
        PolygonCallback2D.__init__(self)
        master_printf("Callback function for epsilon contructed. Now defining the polygons...\n")
        grid_step = 1.0 / float(res)
        if pIsWgBent:
           master_printf("Starting polygons for case without bend.\n")
           #Points outside this polygon are assumed to have epsilon 1.0
           pol_points = numpy.zeros((7,2))
           pol_points[0] = [0.0,padSize]
           pol_points[1] = [gridSizeX - padSize, padSize]
           pol_points[2] = [gridSizeX - padSize, gridSizeY]
           pol_points[3] = [gridSizeX - padSize - wgWidth, gridSizeY]
           pol_points[4] = [gridSizeX - padSize - wgWidth, padSize + wgWidth]
           pol_points[5] = [0.0, padSize + wgWidth]
           pol_points[6] = [0.0,padSize]
           self.add_polygon(pol_points, 12.0)
        else:
           master_printf("Starting polygons for case without bend.\n")
           #Points outside this polygon are assumed to have epsilon 1.0
           pol_points = numpy.zeros((5,2))
           pol_points[0] = [0.0, padSize]
           pol_points[1] = [gridSizeX, padSize]
           pol_points[2] = [gridSizeX, padSize+wgWidth]
           pol_points[3] = [0.0, padSize+wgWidth]
           pol_points[4] = [0.0, padSize]
           self.add_polygon(pol_points, 12.0)

           master_printf("Polygons OK for case without bend.\n")


class amplitudeFactor(Callback):
    def __init__(self):
        Callback.__init__(self)
        master_printf("Callback function for amplitude factor activated.\n")

    def complex_vec(self,vec):
        #BEWARE, these are coordinates RELATIVE to the source center !!!!
        x = vec.x()
        y = vec.y()
        master_printf("Fetching amplitude factor for x=%f - y=%f\n" %(x,y) )
        result = complex(1.0,0)
        return result

def createField(pCompVol, pWgLengthX, pWgWidth, pIsWgBent, pSrcFreqCenter, pSrcPulseWidth, pSrcComp):
        #we create a structure with PML of thickness = 1.0 on all boundaries,
        #in all directions,
        #using the material function EPS
        material = epsilon(pIsWgBent)
        set_EPS_Callback(material.__disown__())
        #use_averaging(False)  #--> uncomment this line to disable eps-averaging!
        s = structure(pCompVol, EPS, pml(1.0))
        f = fields(s)
        #define a gaussian line source of length 'wgWidth' at X=wgLength/2, Y=padSize
        srcGaussian = gaussian_src_time(pSrcFreqCenter, pSrcPulseWidth )
        srcGeo = volume(vec(1.0,padSize),vec(1.0,padSize+pWgWidth))
        af = amplitudeFactor()
        set_AMPL_Callback(af.__disown__())
        f.add_volume_source(pSrcComp, srcGaussian, srcGeo, AMPL)
        master_printf("Field created...\n")
        return f


master_printf("** Bent waveguide sample using POLYGONS to geometry, version 28-02-2011 **\n")

master_printf("Running on %d processor(s)...\n",count_processors())

#FIRST WE WORK OUT THE CASE WITH NO BEND
#----------------------------------------------------------------
master_printf("*1* Starting the case with no bend...\n")
#create the computational grid
noBendVol = voltwo(gridSizeX,gridSizeY,res)

#create the field
wgBent = 0 #no bend
noBendField = createField(noBendVol, wgLengthX, wgWidth, wgBent, srcFreqCenter, srcPulseWidth, srcComp)

bendFnEps = "./bentwgNB_Eps.h5"
bendFnEz = "./bentwgNB_Ez.h5"
#export the dielectric structure (so that we can visually verify the waveguide structure)
noBendDielectricFile =  prepareHDF5File(bendFnEps)
noBendField.output_hdf5(Dielectric, noBendVol.surroundings(), noBendDielectricFile)

#create the file for the field components
noBendFileOutputEz = prepareHDF5File(bendFnEz)

#define the flux plane for the reflected flux
noBendReflectedfluxPlaneXPos = 1.5 #the X-coordinate of our reflection flux plane
noBendReflectedFluxPlane = volume(vec(noBendReflectedfluxPlaneXPos,wgHorYCen-wgWidth),vec(noBendReflectedfluxPlaneXPos,wgHorYCen+wgWidth))
noBendReflectedFlux = noBendField.add_dft_flux_plane(noBendReflectedFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100)

#define the flux plane for the transmitted flux
noBendTransmfluxPlaneXPos = gridSizeX - 1.5;  #the X-coordinate of our transmission flux plane
noBendTransmFluxPlane = volume(vec(noBendTransmfluxPlaneXPos,wgHorYCen-wgWidth),vec(noBendTransmfluxPlaneXPos,wgHorYCen+wgWidth))
noBendTransmFlux = noBendField.add_dft_flux_plane(noBendTransmFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100 )

master_printf("Calculating...")
noBendProbingpoint = vec(noBendTransmfluxPlaneXPos,wgHorYCen) #the point at the end of the waveguide that we want to probe to check if source has decayed
runUntilFieldsDecayed(noBendField, noBendVol, srcComp, noBendProbingpoint, pHDF5OutputFile = noBendFileOutputEz, pH5OutputIntervalSteps=50)
master_printf("Done..!")

#construct 2-dimensional array with the flux plane data
#see python_meep.py
noBendReflFlux = getFluxData(noBendReflectedFlux)
noBendTransmFlux = getFluxData(noBendTransmFlux)

#save the reflection flux from the "no bend" case as minus flux in the temporary file 'minusflux.h5'
noBendReflectedFlux.scale_dfts(-1);
f = open("minusflux.h5", 'w') #truncate file if already exists
f.close()
noBendReflectedFlux.save_hdf5(noBendField, "minusflux")

del_EPS_Callback()


#AND SECONDLY FOR THE CASE WITH BEND
#----------------------------------------------------------------
master_printf("*2* Starting the case with bend...\n")
#create the computational grid
bendVol = voltwo(gridSizeX,gridSizeY,res)

#create the field
wgBent = 1 #there is a bend
bendField = createField(bendVol, wgLengthX, wgWidth, wgBent, srcFreqCenter, srcPulseWidth, srcComp)

#export the dielectric structure (so that we can visually verify the waveguide structure)
bendFnEps = "./bentwgB_Eps.h5"
bendFnEz = "./bentwgB_Ez.h5"
bendDielectricFile = prepareHDF5File(bendFnEps)
bendField.output_hdf5(Dielectric, bendVol.surroundings(), bendDielectricFile)

#create the file for the field components
bendFileOutputEz = prepareHDF5File(bendFnEz)

#define the flux plane for the reflected flux
bendReflectedfluxPlaneXPos = 1.5 #the X-coordinate of our reflection flux plane
bendReflectedFluxPlane = volume(vec(bendReflectedfluxPlaneXPos,wgHorYCen-wgWidth),vec(bendReflectedfluxPlaneXPos,wgHorYCen+wgWidth))
bendReflectedFlux = bendField.add_dft_flux_plane(bendReflectedFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100)

#load the minused reflection flux from the "no bend" case as initalization
bendReflectedFlux.load_hdf5(bendField, "minusflux")


#define the flux plane for the transmitted flux
bendTransmfluxPlaneYPos = padSize + wgLengthY - 1.5; #the Y-coordinate of our transmission flux plane
bendTransmFluxPlane = volume(vec(wgVerXCen - wgWidth,bendTransmfluxPlaneYPos),vec(wgVerXCen + wgWidth,bendTransmfluxPlaneYPos))
bendTransmFlux = bendField.add_dft_flux_plane(bendTransmFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100 )

master_printf("Calculating...\n")
bendProbingpoint = vec(wgVerXCen,bendTransmfluxPlaneYPos) #the point at the end of the waveguide that we want to probe to check if source has decayed
runUntilFieldsDecayed(bendField, bendVol, srcComp, bendProbingpoint, pHDF5OutputFile = bendFileOutputEz)
master_printf("Done..!")

#construct 2-dimensional array with the flux plane data
#see python_meep.py
bendReflFlux = getFluxData(bendReflectedFlux)
bendTransmFlux = getFluxData(bendTransmFlux)

del_EPS_Callback()

#SHOW THE RESULTS IN A PLOT
frequencies = bendReflFlux[0] #should be equal to bendTransmFlux.keys() or noBendTransmFlux.keys() or ...
Ptrans = [x / y for x,y in zip(bendTransmFlux[1], noBendTransmFlux[1])]
Prefl = [ abs(x / y) for x,y in zip(bendReflFlux[1], noBendTransmFlux[1]) ]
Ploss = [ 1-x-y for x,y in zip(Ptrans, Prefl)]

matplotlib.pyplot.plot(frequencies, Ptrans, 'bo')
matplotlib.pyplot.plot(frequencies, Prefl, 'ro')
matplotlib.pyplot.plot(frequencies, Ploss, 'go' )

matplotlib.pyplot.show()

8.3 With inline C/C++ for the material definition

As described in chapter 3, the use of this technique is discouraged since Python-Meep version 1.3. Therefore, these examples are no longer part of the documentation. The source code can however still be found in the source code distribution, under the samples\bent_waveguide path.

9. Running in MPI mode (multiprocessor configuration)

  • We assume that an MPI implementation is installed on your machine (e.g. OpenMPI).

If you want to run python-meep in MPI mode, then you must must import the meep_mpi module instead of the meep module :

from meep_mpi import *

Then start up the Python script as follows :

mpirun -n 2 ./myscript.py

The -n parameter indicates the number of processors requested.

  • For printing output to the console, use the master_printf statement. This will generate output on the master node only (regular Python print statements will run on all nodes).
  • If you output HDF5 files, make sure your HDF5 library is MPI-enable, otherwise your Python script will stall upon writing HDF5 due to a deadlock.

10. Differences between Python-Meep and Scheme-Meep (libctl)

note: this section does NOT apply to the UGent Intec Photonics Research Group (apart from the coordinate system, the default behaviour for us is made consistent with Scheme-Meep)

The general rule is that Python-Meep has consistent behaviour with the C++ core of Meep.

The default version of Python-Meep has 3 differences compared with the Scheme version of Meep :

  • in Python-meep, eps-averaging is disabled by default (see section 3.2 for details on how to enable eps-averaging)
  • in Python-meep, calculation is done with complex fields by default (in Scheme-Meep v1.1.1, real fields are used by default). You can call function use_real_fields() on your fields-object to enable calculation with real fields only.
  • in Python-meep, the center of the coordinate system is in the lower left corner (in Scheme-Meep v1.1.1, the center of the coordinate system is in the middle of your computational volume). However, by calling the function center_origin() on the computational volume, you can shift the origin to the center. However, the use of this function is discouraged, since it was not sufficiently tested yet in combination with Python-Meep.
#define a computational volume of size 16x32, with resolution 10
v = voltwo(16,32,10)
#shift the origin of the coordinate system to the center of the volume
v.center_origin()

On starting your script, Python-Meep will warn you about these differences. You can suppress these warning by setting the global variables DISABLE_WARNING_COORDINATESYSTEM, DISABLE_WARNING_EPS_AVERAGING and DISABLE_WARNING_REAL_FIELDS to True. You can add site-specific customisations to the file meep-site-init.py : in this script, you can for example suppress the warning, or enable EPS-averaging by default.

Add the following code to meep-site-init.py if you want to calculate with real fields only by default :

#by default enable calculation with real fields only (consistent with Scheme-meep)
def fields(structure, m = 0, store_pol_energy=0):
   f = fields_orig(structure, m, store_pol_energy)
   f.use_real_fields()
   return f

#add a new construct 'fields_complex' that you can use to force calculation with complex fields
def fields_complex(structure, m = 0, store_pol_energy=0):
   master_printf("Calculation with complex fields enalbed.\n")
   return fields_orig(structure, m, store_pol_energy)

global DISABLE_WARNING_REAL_FIELDS
DISABLE_WARNING_REAL_FIELDS = True

Add the following code to meep-site-init.py if you want to enable EPS-averaging by default :

use_averaging(True)

global DISABLE_WARNING_EPS_AVERAGING
DISABLE_WARNING_EPS_AVERAGING = True

11. Feedback and comments

You may post your questions, comments or remarks about this documentation on the general Meep mailinglist : meep-discuss@ab-initio.mit.edu Or send feedback directly to Emmanuel.Lambert@intec.ugent.be and nizamov.shawkat@gmail.com (we would love to hear from you).

_images/python-meep.jpg

Table Of Contents

This Page

web counter