subpixel smoothing for MaterialGrid causes modified results and field instabilities when used with symmetry

[oskooi]

A follow up to the discussion in #1780 related to subpixel smoothing for the MaterialGrid involving a demonstration of a continuous frequency response from continuous changes in the input geometry. Following @smartalecH's suggestion of increasing the filter radius (to four pixels rather than one), we are now finally able to demonstrate nearly exact agreement for the example involving a resonant H_z-polarized mode of a ring resonator represented using two Cylinder objects and the MaterialGrid. However, the only caveat for producing these results is that symmetries (i.e., the two mirror symmetries in x and y) had to be turned off for the MaterialGrid. This is because either the symmetries produced slightly different results than the no-symmetry case or, for certain values of the ring radius, were causing the fields to blow up. Not being able to use symmetries to reduce the size of the computation (in this case by a factor of 4) when using a MaterialGrid to represent the ring resonator as a level-set function is obviously a major performance limitation.

The next validation test is to demonstrate quadratic convergence of the gradient (using the directional derivative) versus resolution.

MaterialGrid

    design_shape = mp.Vector3(sxy,sxy,0)
    design_region_resolution = 100
    Nx = int(design_region_resolution*design_shape.x) + 1
    Ny = int(design_region_resolution*design_shape.y) + 1
    x = np.linspace(-0.5*design_shape.x,0.5*design_shape.x,Nx)
    y = np.linspace(-0.5*design_shape.y,0.5*design_shape.y,Ny)
    xv, yv = np.meshgrid(x,y)

    weights = np.logical_and(np.sqrt(np.square(xv) + np.square(yv)) > rad,
                             np.sqrt(np.square(xv) + np.square(yv)) < rad+w,
                             dtype=np.double)

    filter_radius = 4/res
    filtered_weights = mpa.conic_filter(weights,
                                        filter_radius,
                                        design_shape.x,
                                        design_shape.y,
                                        design_region_resolution)

    matgrid = mp.MaterialGrid(mp.Vector3(Nx,Ny),
                              mp.air,
                              mp.Medium(index=n),
                              weights=filtered_weights,
                              do_averaging=True,
                              beta=1000,
                              eta=0.5)

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

Cylinder

    geometry = [mp.Cylinder(material=mp.Medium(index=n),
                            radius=rad+w,
                            height=mp.inf,
                            center=mp.Vector3()),
                mp.Cylinder(material=mp.vacuum,
                            radius=rad,
                            height=mp.inf,
                            center=mp.Vector3())]

[smartalecH]

Very nice!

Here are some initial thoughts:

• When implementing simulation symmetries in conjunction with MaterialGrids, I've found you almost always need to also force symmetry with the underlying MaterialGrid Block object by overlaying multiple mirrored versions of itself and using the U_MEAN flag (just like we do with the tutorials). In an ideal world, we shouldn't have to do that, as the underlying design parameters themselves are symmetric. But with all of the interpolation steps and the way the permittivity is sampled from a discrete array, I'm wondering if there are a few more "necessary conditions" that we aren't thinking about (e.g. odd vs even number of samples, where the center lines up, etc.)
• I think it would be good to overlay 4 curves on your figure: (1) the current curve you have which shows the method works; (2) the case where we only filter with a radius of 1 sample; (3) the case where we don't do any filtering; (4) the case where we don't use any subpixel smoothing, and the underlying β=0 (but we manually project at a high β before the interpolation, thereby implementing a classic density method). We should start to see a gradual breakdown in the continuitity with each curve. The final curve should have strong, discrete steps.

[oskooi]

I've found you almost always need to also force symmetry with the underlying MaterialGrid Block object by overlaying multiple mirrored versions of itself and using the U_MEAN flag (just like we do with the tutorials)

Even when explicitly adding symmetry to the MaterialGrid object, the results are noticeably different compared to MaterialGrid with no symmetry and the Cylinder object. Here are results for the resonant frequency and quality factor (Q) for the ring resonator with an inner ring radius of 1.8 μm for the MaterialGrid and Cylinder with and without symmetry:

1. MaterialGrid with symmetry: 0.263043235778326, 8934.286516094802
2. MaterialGrid with no symmetry: 0.26356255323163896, 2910.5877377395595
3. Cylinder with symmetry: 0.26354500641617595, 10892.830346146491
4. Cylinder with no symmetry: 0.2635450062651281, 10892.711997833718

Note that the frequency for (1) differs from (2), (3), and (4) in the fourth decimal digit. The relative error is less than 1%.

MaterialGrid with symmetry

    matgrid = mp.MaterialGrid(mp.Vector3(Nx,Ny),
                              mp.air,
                              mp.Medium(index=n),
                              weights=filtered_weights,
                              grid_type='U_MEAN',
                              do_averaging=True,
                              beta=1000,
                              eta=0.5)

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

    symmetries = [mp.Mirror(mp.X,phase=+1),
                  mp.Mirror(mp.Y,phase=-1)]

    sim = mp.Simulation(cell_size=cell_size,
                        geometry=geometry,
                        sources=src,
                        resolution=res,
                        symmetries=symmetries,
                        boundary_layers=[mp.PML(dpml)])

[oskooi]

Additional results for MaterialGrid with (1) no eps. averaging and (2) no filter which are nearly equivalent and show staircasing artifacts, as expected.

1. MaterialGrid with no eps. averaging

    weights = np.logical_and(np.sqrt(np.square(xv) + np.square(yv)) > rad,
                             np.sqrt(np.square(xv) + np.square(yv)) < rad+w,
                             dtype=np.double)

    filter_radius = 4/res
    filtered_weights = mpa.conic_filter(weights,
                                        filter_radius,
                                        design_shape.x,
                                        design_shape.y,
                                        design_region_resolution)

    beta = 1000
    eta = 0.5
    projected_weights = mpa.tanh_projection(filtered_weights, beta, eta)

    matgrid = mp.MaterialGrid(mp.Vector3(Nx,Ny),
                              mp.air,
                              mp.Medium(index=n),
                              weights=projected_weights,
                              do_averaging=False,
                              beta=0)

2. MaterialGrid with no filter

    weights = np.logical_and(np.sqrt(np.square(xv) + np.square(yv)) > rad,
                             np.sqrt(np.square(xv) + np.square(yv)) < rad+w,
                             dtype=np.double)

    matgrid = mp.MaterialGrid(mp.Vector3(Nx,Ny),
                              mp.air,
                              mp.Medium(index=n),
                              weights=weights,
                              do_averaging=False,
                              beta=0)