user-specified chunk layouts

[oskooi]

Closes #1510.

This PR adds a new BinaryPartition class to the Python API following the outline proposed by @stevengj in #1510 (comment) which can be used to specify a cell partition comprised of arbitrary sized chunks. A BinaryPartition class object is then passed as the chunk_layouts parameter of the Simulation class object.

A new binary_partition class is added in C++ which is designed to be used with a new split_by_binarytree routine to divide the cell into chunks similar to the existing split_by_cost routine.

This PR is missing the SWIG typemaps (or a glue routine) to convert the BinaryPartition object from Python into a binary_partition object in C++. It is therefore not yet ready to be merged. This PR, when ready, will also include a unit test.

[stevengj]

For a SWIG wrapper you need something like the following in meep.i:

%typemap(in) binary_partition * {
    $1 = py_bp_to_bp($input);
    if(!$1) {
        SWIG_fail;
    }
}

%typemap(arginit) binary_partition * {
    $1 = NULL;
}

%typemap(freearg) binary_partition * {
    delete $1;
}

where you have defined a function:

binary_partition *py_bp_to_bp(PyObject *bp) {
    ....
}

that converts the Python object to the C++ object.

[stevengj]

For py_bp_to_bp you have to use the Python C API, something like:

binary_partition *py_bp_to_bp(PyObject *bp) {
    binary_partition *bp = NULL;
    if (bp == Py_None) return bp;

    PyObject *id = PyObject_GetAttrString(bp, "id");
    PyObject *split_dir = PyObject_GetAttrString(bp, "split_dir");
    PyObject *split_pos = PyObject_GetAttrString(bp, "split_pos");
    PyObject *left = PyObject_GetAttrString(bp, "left");
    PyObject *right = PyObject_GetAttrString(bp, "right");

    if (!id || !split_dir || !split_pos || !left || !right) {
        // error....
    }

    if (PyLong_Check(id)) {
         bp = new binary_partition(PyLong_AsLong(id));
    else {
         bp = new binary_partition(direction(PyLong_AsLong(split_dir)), PyFloat_AsDouble(split_pos));
         bp->left = py_bp_to_bp(left);
         bp->right = py_bp_to_bp(right);
    }

    Py_XDECREF(id);
    Py_XDECREF(split_dir);
    Py_XDECREF(split_pos);
    Py_XDECREF(left);
    Py_XDECREF(right);
    return bp; 
}

[stevengj]

We could also add a

%typecheck (SWIG_TYPECHECK_POINTER) binary_partition * {
    $1 = PyObject_IsInstance($input, py_binary_partition_object());
}

where py_binary_partition_type is defined in typemap_utils.cpp as:

static PyObject *get_meep_mod() {
  // Return value: Borrowed reference
  static PyObject *meep_mod = NULL;
  if (meep_mod == NULL) { meep_mod = PyImport_ImportModule("meep"); }
  return meep_mod;
}

static PyObject *py_binary_partition_object() {
  // Return value: Borrowed reference
  static PyObject *bp_type = NULL;
  if (bp_type == NULL) {
    bp_type = PyObject_GetAttrString(get_meep_mod(), "BinaryPartition");
  }
  return bp_type;
}

[oskooi]

The following error appears during make with the above changes applied to meep.i and typemap_utils.cpp:

meep-python.cxx: In function ‘PyObject* _wrap_py_bp_to_bp(PyObject*, PyObject*)’:
meep-python.cxx:130336:42: error: ‘py_bp_to_bp’ was not declared in this scope
       result = (meep::binary_partition *)py_bp_to_bp(arg1);

[oskooi]

Moving py_bp_to_bp from meep.i into typemap_utils.cpp and running make produces this error:

typemap_utils.cpp: In function ‘meep::binary_partition* py_bp_to_bp(PyObject*)’:
typemap_utils.cpp:1034:29: error: declaration of ‘meep::binary_partition* bp’ shadows a parameter
     meep::binary_partition *bp = NULL;

[stevengj]

Alternatively, you could just write a function

def binarypartition(list):
    ....

that just constructs the C++ meep::binary_partition object directly (by calling its SWIG-wrapped constructors), which you can then pass directly to C++ (with no additional typemaps or conversion functions).

[oskooi]

This feature seems to be working now. The current setup involves creating a BinaryPartition class object in Python and passing it as the chunk_layout parameter of the Simulation constructor.

A unit test has been added which verifies that the chunk ids and volumes specified in a given BinaryPartition class object (based on the example from #1510) are equivalent to the actual chunks created during the structure class initialization.

As another verification, the chunk layout for the 2d test case can be visualized using the visualize_chunks routine:

import meep as mp
import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt

chunk_layout = mp.BinaryPartition(data=[ (mp.X,-2.0), 0, [ (mp.Y,1.5), [ (mp.X,3.0), 1, [ (mp.Y,-0.5), 4,3 ] ], 2 ] ])

cell_size = mp.Vector3(10.0,5.0,0)

sim = mp.Simulation(cell_size=cell_size,
		    resolution=10,
                    chunk_layout=chunk_layout)

sim.visualize_chunks()
plt.savefig('bp_chunk_layout.png',dpi=150,bbox_inches='tight')

Executing this script via mpirun -np 5 python3 chunks.py produces this image:

The only differences in this figure compared to that from #1510 (comment) are the cell origin which is (0,0) rather than (5,2.5) and the five chunk ids which are in the range [0,4] rather than [1,5].

[oskooi]

An example demonstrating the use of this feature (based on the results shown above) has been added as a new section to the Parallel Meep page of the user manual.

[stevengj]

Key point: this is not a "chunk ID", it is a "process ID" which must be between 0 and #processes–1 (inclusive). Note also that the same process ID can be assigned to as many chunks as you want, which just means that process timesteps multiple chunks.

[stevengj]

One option would be to change split_by_const etcetera so that they actually return a binary_partition and then we call split_by_binarytree. That would make it easy to export the tree if we want, and would ensure that the binary_partition code gets further testing (because now all splitting would pass through it).

[stevengj]

No need for num_chunks, but it won't actually use more than 5 processes.

[stevengj]

I think the binary_partition needs to be passed to the structure constructor, rather than choosing the chunk division and then overwriting it with the new chunk division in load_chunk_layout as you are doing now. (e.g. you shouldn't have to pass num_chunks as an additional parameter — it should be determined from the partition.)

[oskooi]

Modifying the structure class constructor to accept a binary_partition parameter and revamping split_by_cost to return a binary_partition object would require a lot of changes which should probably be a separate PR. For now, I just added the following comment to the doc strings for BinaryPartition and the tutorial example:

The `num_chunks` parameter of the `Simulation` constructor must be set to the
number of chunks specified by the binary tree.

As long as the num_chunks parameter is set to the number of chunks in the binary tree, the simulation will run with any number of MPI processes (i.e., run using not just calling visualize_chunks but actually time stepping the fields).