segmentation fault for 3D design region of adjoint solver

[oskooi]

It seems #1855 introduced a bug in the adjoint solver for the case of a 3D design region. (I discovered this while working on a unit test for #2054.) This can be demonstrated using the following test which involves a 3D unit cell of a cylindrical grating on a substrate with 2D periodicity. The MaterialGrid is actually 2D but the DesignRegion is 3D with extrusion in the z direction. The test involves an objective function based on the mode coefficient of a transmitted mode.

With the current master branch (commit: 6c975ac), this test crashes with a segmentation fault that is traced using gdb to the following lines introduced in #1855:

Calculating gradient...
*** Process received signal ***
Signal: Segmentation fault (11)
Signal code:  (128)
Failing at address: (nil)
*** Process received signal ***
Signal: Segmentation fault (11)
Signal code: Address not mapped (1)
Failing at address: 0x101

meep/src/meepgeom.cpp

Lines 3023 to 3028 in 6c975ac

	// clear all the dft data structures
	for (int i=0;i<3;i++){
	for (int ii=0;i<adjoint_dft_chunks[i].size();i++){
	delete adjoint_dft_chunks[i][ii];
	delete forward_dft_chunks[i][ii];
	}

With the commit just prior to #1855 (75ced52), the test runs fine and the gradient is evaluated without any errors:

Calculating gradient...
dir-deriv:, [-1.32462946e-15]

import meep as mp
import meep.adjoint as mpa
import math
import numpy as np
from autograd import numpy as npa


resolution = 25  # pixels/μm                                                                                      

nSi = 3.45
Si = mp.Medium(index=nSi)
nSiO2 = 1.45
SiO2 = mp.Medium(index=nSiO2)

wvl = 0.5  # wavelength                                                                                           
fcen = 1/wvl

dpml = 1.0  # PML thickness                                                                                       
dsub = 3.0  # substrate thickness                                                                                 
dair = 3.0  # air padding                                                                                         
hcyl = 0.5  # cylinder height                                                                                     
rcyl = 0.2  # cylinder radius                                                                                     

sx = 1.4
sy = 0.8
sz = dpml+dsub+hcyl+dair+dpml

cell_size = mp.Vector3(sx,sy,sz)

design_region_size = mp.Vector3(sx,sy)
design_region_resolution = int(2*resolution)
Nx = int(design_region_size.x*design_region_resolution)
Ny = int(design_region_size.y*design_region_resolution)

# ensure reproducible results                                                                                     
rng = np.random.RandomState(59387385)

xcoord = np.linspace(-0.5*design_region_size.x,0.5*design_region_size.x,Nx)
ycoord = np.linspace(-0.5*design_region_size.y,0.5*design_region_size.y,Ny)
xv, yv = np.meshgrid(xcoord,ycoord)

# cylindrical grating                                                                                             
p = np.sqrt(np.square(xv) + np.square(yv)) < rcyl
p = p.flatten(order='F')

# random epsilon perturbation for design region                                                                   
deps = 1e-5
dp = deps*rng.rand(Nx*Ny)

boundary_layers = [mp.PML(thickness=dpml,direction=mp.Z)]

# periodic boundary conditions                                                                                    
k_point = mp.Vector3()

# plane of incidence: XZ                                                                                          
# polarization: (1) S/TE: Ey or (2) P/TM: Ex                                                                      
src_cmpt = mp.Ex
sources = [mp.Source(src=mp.GaussianSource(fcen,fwidth=0.2*fcen),
                     size=mp.Vector3(sx,sy,0),
                     center=mp.Vector3(0,0,-0.5*sz+dpml),
                     component=src_cmpt)]

tran_pt = mp.Vector3(0,0,0.5*sz-dpml)
stop_cond = mp.stop_when_fields_decayed(10,src_cmpt,tran_pt,1e-7)

substrate = [mp.Block(size=mp.Vector3(mp.inf,mp.inf,dpml+dsub),
                      center=mp.Vector3(0,0,-0.5*sz+0.5*(dpml+dsub)),
                      material=SiO2)]

def adjoint_solver(design_params):
    matgrid = mp.MaterialGrid(mp.Vector3(Nx,Ny,0),
                              mp.air,
                              Si,
                              weights=np.ones((Nx,Ny)))

    matgrid_region = mpa.DesignRegion(matgrid,
                                      volume=mp.Volume(center=mp.Vector3(0,0,-0.5*sz+dpml+dsub+0.5*hcyl),
                                                       size=mp.Vector3(design_region_size.x,
                                                                       design_region_size.y,
                                                                       hcyl)))

    matgrid_geometry = [mp.Block(center=matgrid_region.center,
                                 size=matgrid_region.size,
                                 material=matgrid)]

    geometry = substrate + matgrid_geometry

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

    obj_list = [mpa.EigenmodeCoefficient(sim,
                                         mp.Volume(center=tran_pt,
                                                   size=mp.Vector3(sx,sy,0)),
                                         1,
                                         eig_parity=mp.ODD_Y)]

    def J(mode_mon):
        return npa.power(npa.abs(mode_mon),2)

    opt = mpa.OptimizationProblem(simulation=sim,
                                  objective_functions=J,
                                  objective_arguments=obj_list,
                                  design_regions=[matgrid_region],
                                  frequencies=[fcen],
                                  maximum_run_time=100)

    f, dJ_du = opt([design_params])

    sim.reset_meep()

    return f, dJ_du


adjsol_obj, adjsol_grad = adjoint_solver(p)

if adjsol_grad.ndim < 2:
    adjsol_grad = np.expand_dims(adjsol_grad,axis=1)
adj_scale = (dp[None,:]@adjsol_grad).flatten()
print("dir-deriv:, {}".format(adj_scale))

[smartalecH]

Odd since this is just the cleanup step after all the gradient calculations are performed on each chunk... I'll take a closer look. Nice catch.

[mawc2019]

A 2d example showing the segmentation fault was posted in #1855.
In addtion, in a test for the FourierFields adjoint solver that was also posted there, both finite-difference and adjoint gradients change obviously once this PR is used. These finite-difference and adjoint gradients become inconsistent after this change.