enabling support for subpixel smoothing of MaterialGrid

[oskooi]

The MaterialGrid feature added by @smartalecH in #1242 (which is not documented in the User Interface and is intended mainly for use with the adjoint solver) provides support for importing the material geometry as a NumPy array. One limitation of the current implementation is that it is based on a bilinear interpolation of the user-specified pixel grid onto Meep's Yee grid which is therefore only first-order accurate (rather than second order) with respect to the resolution:

meep/src/meepgeom.cpp

Lines 693 to 701 in 7067536

	// material grid: interpolate onto user specified material grid to get properties at r
	case material_data::MATERIAL_GRID:
	meep::realnum u;
	int oi;
	geom_box_tree tp;
	tp = geom_tree_search(p, restricted_tree, &oi);
	u = matgrid_val(p, tp, oi, md); // interpolate onto material grid
	epsilon_material_grid(md, u); // interpolate material from material grid point

The accuracy of the MaterialGrid feature can be improved by replacing the bilinear interpolation with anisotropic subpixel smoothing using a quadrature scheme similar to what already exists for the material function but without the slow callbacks to Python since it would be implemented entirely in C/C++.

As a demonstration of the convergence properties of the current implementation of the MaterialGrid, we can compute the relative error in the resonant mode of a 2d photonic crystal using a unit cell comprised of a cylinder with n=3.5 in air and non-zero Bloch-periodic boundaries (so that the mode is propagating at an oblique angle to the Cartesian grid). The source polarization is H_z (i.e., in-plane electric field) in order to generate non-parallel fields along the discontinuous dielectric interface. (This type of convergence test was originally demonstrated in Optics Letters, vol. 31, pp. 2972-74 (2006)). For comparison, we can also compute the same result using a Cylinder geometric object which already supports subpixel smoothing.

The following schematic shows the H_z profile for the Cylinder and MaterialGrid geometry representations at three different resolutions. The discrete artifacts of the MaterialGrid are slightly more apparent than the Cylinder representation. (These images were generated using plot2D and the red dots are the source locations.)

The simulation involves using Harminv to compute the frequency of the resonant mode over a range of resolutions (20-200). The "exact' result is computed at a resolution of 300.

import numpy as np
import argparse
import meep as mp
import matplotlib

parser = argparse.ArgumentParser()
parser.add_argument('-res',
                    type=int,
                    default=20,
                    help='resolution (default: 20 pixels/um)')
parser.add_argument('-geom_type',
                    type=int,
                    choices=[1,2],
                    default=1,
                    help='type of geometry: 1: Cylinder (default), 2: material grid')
args = parser.parse_args()

cell_size = mp.Vector3(1,1,0)

rad = 0.301943

if args.geom_type == 1:
    geometry = [mp.Cylinder(radius=rad,
                            center=mp.Vector3(),
                            height=mp.inf,
                            material=mp.Medium(index=3.5))]
else:
    design_shape = mp.Vector3(1,1,0)
    design_region_resolution = int(2*args.res)
    Nx = int(design_region_resolution*design_shape.x)
    Ny = int(design_region_resolution*design_shape.y)
    x = np.linspace(-0.5*cell_size.x,0.5*cell_size.x,Nx)
    y = np.linspace(-0.5*cell_size.y,0.5*cell_size.y,Ny)
    xv, yv = np.meshgrid(x,y)
    design_params = np.sqrt(np.square(xv) + np.square(yv)) < rad

    matgrid = mp.MaterialGrid(mp.Vector3(Nx,Ny),
                              mp.air,
                              mp.Medium(index=3.5),
                              design_parameters=design_params,
                              grid_type='U_SUM')

    geometry = [mp.Block(center=mp.Vector3(),
                         size=mp.Vector3(design_shape.x,design_shape.y,0),
                         material=matgrid)]

fcen = 0.3
df = 0.2*fcen
sources = [mp.Source(mp.GaussianSource(fcen,fwidth=df),
                     component=mp.Hz,
                     center=mp.Vector3(-0.1057,0.2094,0))]

k_point = mp.Vector3(0.3892,0.1597,0)

sim = mp.Simulation(resolution=args.res,
                    cell_size=cell_size,
                    geometry=geometry,
                    sources=sources,
                    k_point=k_point)

h = mp.Harminv(mp.Hz, mp.Vector3(0.3718,-0.2076), fcen, df)
sim.run(mp.after_sources(h),
        until_after_sources=200)

for m in h.modes:
    print("harminv:, {}, {}, {}".format(args.res,m.freq,m.Q))

A plot of the relative error vs. resolution for the two methods demonstrates that the MaterialGrid is indeed first-order accurate whereas the Cylinder is second order.

To verify that anisotropic subpixel smoothing has been correctly set up for the MaterialGrid, we can use this same convergence test to verify that it is second-order accurate.