Add a test case for fix_boundary_sources

[stevengj]

@mawc2019, would be good to add a test case for #1959 (that fails on CI without this PR).

[mawc2019]

Here is an example. Without this PR, the finite-difference and adjoint gradients are inconsistent when the number of processes is large, such as 8. With this PR, the finite-difference and adjoint gradients become consistent.

The code is as follows.

import meep as mp
import meep.adjoint as mpa
import numpy as np
from autograd import numpy as npa
Si, SiO2 = mp.Medium(index=3.4), mp.Medium(index=1.44)

resolution,design_resolution = 20,20
Lx,Ly,pml_size = 2,4,1
cell_size = mp.Vector3(Lx,Ly+2*pml_size)
pml_layers = [mp.PML(pml_size, direction=mp.Y)]
design_width,design_height = Lx-0.4,0.4

fcen,fwidth = 1/1.55,0.01/1.55
source_center,source_size = mp.Vector3(0,-Ly/4),mp.Vector3(Lx,0)

src = mp.GaussianSource(frequency=fcen,fwidth=fwidth)
source = [mp.Source(src,component=mp.Ez,center=source_center,size=source_size)]
Nx,Ny = int(round(design_resolution*design_width)),int(round(design_resolution*design_height))

design_variables = mp.MaterialGrid(mp.Vector3(Nx,Ny),SiO2,Si,grid_type='U_MEAN')
design_region = mpa.DesignRegion(design_variables,volume=mp.Volume(center=mp.Vector3(), size=mp.Vector3(design_width, design_height)))

geometry = [mp.Block(center=design_region.center, size=design_region.size, material=design_variables)] # design region
sim = mp.Simulation(
    cell_size=cell_size,boundary_layers=pml_layers,geometry=geometry,
    sources=source,default_material=SiO2,resolution=resolution)

nf = mpa.FourierFields(sim,mp.Volume(center=mp.Vector3(0,design_height+0.2),size=mp.Vector3(Lx,0)),mp.Ez)
ob_list = [nf]

def J(near_data):
    return npa.sum(npa.abs(near_data)**2)

opt = mpa.OptimizationProblem(
    simulation = sim,objective_functions = J,objective_arguments = ob_list,design_regions = [design_region],
    fcen=fcen,df = 0,nf = 1,decay_by = 1e-5)

x0 = 0.5*np.ones(Nx*Ny)
opt.update_design([x0])
f0, dJ_deps = opt()

db, choose = 1e-4, 5
np.random.seed(20211204)

g_discrete, idx = opt.calculate_fd_gradient(num_gradients=choose,db=db)
g_discrete = np.array(g_discrete).flatten()
g_adjoint = dJ_deps.flatten()

print("Randomly selected indices: ",idx)
print("Finite-difference gradient: ",g_discrete)
print("Adjoint gradient: ",g_adjoint[idx])

[stevengj]

I would suggest a simpler test.

1. For a small 2d problem, compute the adjoint gradient for FourierFields with 1 and 2 processors
2. Compare them: they should be the same up to roundoff errors (but different before fix_boundary_sources).

[mawc2019]

For a small 2d problem, compute the adjoint gradient for FourierFields with 1 and 2 processors

When the number of processes is too small, fix_boundary_sources does not make a difference. For the above setup, the smallest number of processes that allows fix_boundary_sources to make a difference is 7.

Compare them: they should be the same up to roundoff errors (but different before fix_boundary_sources).

The comparison is as follows. In the second figure, the relative difference is defined as 7-process adjoint gradient / 1-process adjoint gradient − 1, which is around 1e-13 with fix_boundary_sources.

The code is as follows, which has the same setup as above.

import meep as mp
import meep.adjoint as mpa
import numpy as np
from autograd import numpy as npa
Si, SiO2 = mp.Medium(index=3.4), mp.Medium(index=1.44)

resolution,design_resolution = 20,20
Lx,Ly,pml_size = 2,4,1
cell_size = mp.Vector3(Lx,Ly+2*pml_size)
pml_layers = [mp.PML(pml_size, direction=mp.Y)]
design_width,design_height = Lx-0.4,0.4

fcen,fwidth = 1/1.55,0.01/1.55
source_center,source_size = mp.Vector3(0,-Ly/4),mp.Vector3(Lx,0)

src = mp.GaussianSource(frequency=fcen,fwidth=fwidth)
source = [mp.Source(src,component=mp.Ez,center=source_center,size=source_size)]
Nx,Ny = int(round(design_resolution*design_width)),int(round(design_resolution*design_height))

design_variables = mp.MaterialGrid(mp.Vector3(Nx,Ny),SiO2,Si,grid_type='U_MEAN')
design_region = mpa.DesignRegion(design_variables,volume=mp.Volume(center=mp.Vector3(), size=mp.Vector3(design_width, design_height)))

geometry = [mp.Block(center=design_region.center, size=design_region.size, material=design_variables)] # design region
sim = mp.Simulation(
    cell_size=cell_size,boundary_layers=pml_layers,geometry=geometry,
    sources=source,default_material=SiO2,resolution=resolution)

nf = mpa.FourierFields(sim,mp.Volume(center=mp.Vector3(0,design_height+0.2),size=mp.Vector3(Lx,0)),mp.Ez)
ob_list = [nf]

def J(near_data):
    return npa.sum(npa.abs(near_data)**2)

opt = mpa.OptimizationProblem(
    simulation = sim,objective_functions = J,objective_arguments = ob_list,design_regions = [design_region],
    fcen=fcen,df = 0,nf = 1,decay_by = 1e-5)

x0 = 0.5*np.ones(Nx*Ny)
opt.update_design([x0])
f0, g_adjoint = opt()

number_of_process = 1
np.savetxt("adjoint_proc"+str(number_of_process)+".dat",g_adjoint.flatten())

[stevengj]

To run this in serial and parallel in the same script, do:

# set up simulation

mp.divide_parallel_processes(mp.count_processes())
# run simulation in serial
mp.end_divide_parallel()

# reset simulation and re-run it in parallel

# compare