Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

segmentation fault for 3D design region of adjoint solver #2055

Closed
oskooi opened this issue May 3, 2022 · 2 comments · Fixed by #2066
Closed

segmentation fault for 3D design region of adjoint solver #2055

oskooi opened this issue May 3, 2022 · 2 comments · Fixed by #2066
Labels

Comments

@oskooi
Copy link
Collaborator

oskooi commented May 3, 2022

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))
@oskooi oskooi added the bug label May 3, 2022
@smartalecH
Copy link
Collaborator

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
Copy link
Contributor

mawc2019 commented May 7, 2022

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

Successfully merging a pull request may close this issue.

3 participants