edge artifacts of `MaterialGrid`

[oskooi]

There seems to be unwanted artifacts along the outer edges of a MaterialGrid which could introduce errors into a simulation. This can be demonstrated by using plot2D to visualize a MaterialGrid in a 2d cell consisting of a linear interpolation along a single direction of medium1 (index n1) and medium2 (index n2). The edge artifacts are independent of the resolution. The edge artifacts seem to disappear when the MaterialGrid is defined as the default_medium (rather than a Block object) which suggests there is probably a bug in the interpolation for grid voxels near the edges.

(Possibly related to #1399.)

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

resolution = 50

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

Nx = 101
Ny = 101
lin_interp = np.linspace(0,1,Nx)
weights = np.tile(lin_interp,(Ny,1))
mat_interp = mp.MaterialGrid(mp.Vector3(Nx,Ny,0),
                             mp.Medium(index=1.5),
                             mp.Medium(index=3.5),
                             weights=weights)

geometry = [mp.Block(material=mat_interp,
                     center=mp.Vector3(),
                     size=mp.Vector3(1,1,0))]

sim = mp.Simulation(cell_size=cell_size,
                    resolution=resolution,
                    geometry=geometry)

sim.plot2D()
plt.title('resolution = {}, (Nx,Ny) = ({},{})'.format(resolution,Nx,Ny))
plt.savefig('matgrid_graded_res{}_{}_{}.png'.format(resolution,Nx,Ny),
            bbox_inches='tight',
            dpi=150)

[smartalecH]

Couldn't this be from the get_epsilon routine called by the plotting function? (Remember there's quite a bit of "averaging" that occurs due to the yee grid).

Note that our new approach (#1951) does a better job at handling the edges... so it could indeed be an underlying issue with the current averaging code...

[oskooi]

plot2D calls get_epsilon_grid which loops over a set of grid points: any grid point within a MaterialGrid will use bilinear interpolation to obtain an averaged value from its four neighbors whereas any point outside will not use interpolation. However, the images above show artifacts on the outer edges which do not seem to be interpolated because they are discontinuous.

[oskooi]

I suspect there could be a bug during the linear interpolation somewhere in these lines:

meep/src/fields.cpp

Lines 734 to 773 in 35b10e5

	/* implement mirror boundary conditions for i outside 0..n-1: */
	int mirrorindex(int i, int n) { return i >= n ? 2 * n - 1 - i : (i < 0 ? -1 - i : i); }
	/* map the cell coordinates into the range [0,1].
	anything outside [0,1] is *mirror* reflected into [0,1] */
	void map_coordinates(double rx, double ry, double rz,
	int nx, int ny, int nz,
	int &x1, int &y1, int &z1,
	int &x2, int &y2, int &z2,
	double &dx, double &dy, double &dz,
	bool do_fabs) {
	/* mirror boundary conditions for r just beyond the boundary */
	rx = rx < 0.0 ? -rx : (rx > 1.0 ? 1.0 - rx : rx);
	ry = ry < 0.0 ? -ry : (ry > 1.0 ? 1.0 - ry : ry);
	rz = rz < 0.0 ? -rz : (rz > 1.0 ? 1.0 - rz : rz);
	/* get the point corresponding to r in the epsilon array grid: */
	x1 = mirrorindex(int(rx * nx), nx);
	y1 = mirrorindex(int(ry * ny), ny);
	z1 = mirrorindex(int(rz * nz), nz);
	/* get the difference between (x,y,z) and the actual point */
	dx = rx * nx - x1 - 0.5;
	dy = ry * ny - y1 - 0.5;
	dz = rz * nz - z1 - 0.5;
	/* get the other closest point in the grid, with mirror boundaries: */
	x2 = mirrorindex(dx >= 0.0 ? x1 + 1 : x1 - 1, nx);
	y2 = mirrorindex(dy >= 0.0 ? y1 + 1 : y1 - 1, ny);
	z2 = mirrorindex(dz >= 0.0 ? z1 + 1 : z1 - 1, nz);
	/* take abs(d{xyz}) to get weights for {xyz} and {xyz}2: */
	if (do_fabs) {
	dx = fabs(dx);
	dy = fabs(dy);
	dz = fabs(dz);
	}
	}

[smartalecH]

When I run your example using #1951 (which uses the same interpolation scheme as master) I get the following image (I don't seem to see any artifacts):

As mentioned earlier, I think the issue is with the boundary smoothing (which is fixed in #1951).