material_grid API

[stevengj]

Might be nice to have something like the material_grid feature in MPB where you can have a material function based on a numpy array of ε or u ∈ [0,1] values and it would linearly interpolate onto the Yee grid (implemented in C++ for performance) and then interpolate between two given materials as a function of u.

We might even implement the feature in MPB where you can have overlapping material grids, which is useful for projecting onto symmetry groups during optimization (http://arxiv.org/abs/arXiv:1405.4350) and potentially other applications (e.g. I initially implemented it to model overlapping lithography steps).

[stevengj]

We'll need to design this API so that it can be used by the adjoint solver. That is, it has to expose the interpolation weights or do the product of interpolation weights with fields/currents in C++.

[stevengj]

Given such a grid feature, we might want to implement a simpler subpixel averaging feature that just uses an MxMxM subgrid.

[smartalecH]

I'm trying to map out the API.

It seems like we can adopt an approach that is very similar to the material_file material type. For example, we can create the grid using something like:

material_type make_material_grid(int Nx, int Ny, material_data *m1, material_data *m2);

which is similar to

meep/src/meepgeom.cpp

Lines 1580 to 1584 in 65b69ef

	// this routine subsumes the content of the old
	// 'read_epsilon_file' routine
	material_type make_file_material(const char *eps_input_file) {
	material_data *md = new material_data();
	md->which_subclass = material_data::MATERIAL_FILE;

The update call would also be similar, only with the proper interpolations:

meep/src/meepgeom.cpp

Lines 281 to 299 in 65b69ef

	void epsilon_file_material(material_data *md, vector3 p) {
	default_material = (void *)md;
	if (md->which_subclass != material_data::MATERIAL_FILE)
	meep::abort("epsilon-input-file only works with a type=file default-material");
	if (!(md->epsilon_data)) return;
	medium_struct *mm = &(md->medium);
	double rx =
	geometry_lattice.size.x == 0 ? 0 : 0.5 + (p.x - geometry_center.x) / geometry_lattice.size.x;
	double ry =
	geometry_lattice.size.y == 0 ? 0 : 0.5 + (p.y - geometry_center.y) / geometry_lattice.size.y;
	double rz =
	geometry_lattice.size.z == 0 ? 0 : 0.5 + (p.z - geometry_center.z) / geometry_lattice.size.z;
	mm->epsilon_diag.x = mm->epsilon_diag.y = mm->epsilon_diag.z =
	meep::linear_interpolate(rx, ry, rz, md->epsilon_data, md->epsilon_dims[0],
	md->epsilon_dims[1], md->epsilon_dims[2], 1);
	mm->epsilon_offdiag.x.re = mm->epsilon_offdiag.y.re = mm->epsilon_offdiag.z.re = 0;
	}

Which will need a proper call in get_material_pt:

meep/src/meepgeom.cpp

Lines 516 to 533 in 65b69ef

	// the goal of this routine is to fill in the 'medium' field
	// within the material structure as appropriate for the
	// material properties at r.
	void geom_epsilon::get_material_pt(material_type &material, const meep::vec &r) {
	vector3 p = vec_to_vector3(r);
	boolean inobject;
	material =
	(material_type)material_of_unshifted_point_in_tree_inobject(p, restricted_tree, &inobject);
	material_data *md = material;
	switch (md->which_subclass) {
	// material read from file: interpolate to get properties at r
	case material_data::MATERIAL_FILE:
	if (md->epsilon_data)
	epsilon_file_material(md, p);
	else
	material = (material_type)default_material;
	return;

For the adjoint solver, we can do an implicit vector-jacobian-product (similar to the MPB routine). Obviously we would need proper SWIG hooks too, which shouldn't be too bad thanks to Chris.

Does this sound right? Also, can you explain each of the MPB material_grid enumerations:

https://github.com/NanoComp/mpb/blob/0935b95c362873c2dd871646c5f12899e48af1c6/mpb/mpb.h#L98

[stevengj]

Definitely it should be 3d (Nx x Ny x Nz).

[stevengj]

That enumeration is to support combining overlapping material grids from different objects. U_PROD (and U_MIN) are designed to represent overlapping lithography steps like in this paper. U_SUM (which takes the mean of overlapping grids) is useful for imposing symmetries as in this paper. See here and here in the code.