MaterialGrid and subpixel smoothing

ring2d_harminv.py

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

res = 10                # pixels/μm                                                                               

n = 3.4                 # index of waveguide                                                                      
w = 1                   # width of waveguide                                                                      
pad = 4                 # padding between waveguide and edge of PML                                               
dpml = 2                # thickness of PML                                                                        

for rad in np.arange(1.800,2.001,0.005):
    sxy = 2*(rad+w+pad+dpml)  # cell size                                                                         

    cell_size = mp.Vector3(sxy,sxy)

    design_shape = mp.Vector3(sxy,sxy,0)
    design_region_resolution = 1000
    Nx = int(design_region_resolution*design_shape.x) + 1
    Ny = int(design_region_resolution*design_shape.y) + 1
    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.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=design_params,
                              do_averaging=False,
                              beta=1000,
                              eta=0.5)

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

    # pulse center frequency (from third-order polynomial fit)                                                    
    fcen = -0.018765*rad**3 + 0.137685*rad**2 -0.393918*rad + 0.636202
    # pulse frequency width                                                                                       
    df = 0.02*fcen

    src = [mp.Source(mp.GaussianSource(fcen, fwidth=df),
                     component=mp.Hz,
                     center=mp.Vector3(rad+0.1*w)),
           mp.Source(mp.GaussianSource(fcen, fwidth=df),
                     component=mp.Hz,
                     amplitude=-1,
                     center=mp.Vector3(-(rad+0.1*w)))]

    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)])

    h = mp.Harminv(mp.Hz, mp.Vector3(rad+0.1*w), fcen, df)
    sim.run(mp.after_sources(h),
            until_after_sources=200)

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

res = 10
desres = 1000

f1 = np.genfromtxt('ring_matgrid_simxy_res{}_desres{}_nosubpix.dat'.format(res,desres),delimiter=',')
f2 = np.genfromtxt('ring_matgrid_simxy_res{}_desres{}_subpix.dat'.format(res,desres),delimiter=',')

plt.plot(f1[:,0],f1[:,1],'bo-',label='no smoothing')
plt.plot(f2[:,0],f2[:,1],'ro-',label='subpixel smoothing')

plt.xlabel(r'inner ring radius ($\mu$m)')
plt.ylabel('frequency of resonant mode (300 THz)')
plt.legend()
plt.title('MaterialGrid with resolution {}; Meep resolution {}'.format(desres,res))
plt.savefig('ring_matgrid_res{}_desres{}_freq_vs_res.png'.format(res,desres),bbox_inches='tight',dpi=150)

Smooth the MaterialGrid design region using a Gaussian filter.

    filter_sigma = 1.0
    filtered_design_params = ndimage.gaussian_filter(design_params, filter_sigma, output=np.double)

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