search

NanoComp / meep

free finite-difference time-domain (FDTD) software for electromagnetic simulations
GNU General Public License v2.0
1.25k stars 626 forks source link

Incorrect adjoint fields in cylindrical coordinates #2193

Closed mochen4 closed 2 years ago

mochen4 commented 2 years ago

There is an issue with the gradients after PR #1855. Specifically, if I plot the forward fields, adjoint fields, and coordinates of the points from this line: https://github.com/NanoComp/meep/blob/2aa9164580bdad950ba834df95039eb1e6dbcce7/src/meepgeom.cpp#L2925-L2927

The forward fields are the same as before the PR, yet the adjoint fields are different compared to before. cc @smartalecH

Here is a simple test case (in cylindrical coordinates) I used:

import meep as mp
import meep.adjoint as mpa
import numpy as np
from autograd import numpy as npa
from matplotlib import pyplot as plt

np.random.rand(240)

mp.verbosity(0)
resolution = 20
Si = mp.Medium(epsilon=13)

design_region_width = 0.1
design_region_height = 0.1
pml_size = 1.0
padl, padr, padz = 0.5 , 0.5, 0.5
Sr = design_region_width + pml_size + padl + padr
Sz = 2*pml_size + design_region_height + 2*padz + 1
cell_size = mp.Vector3(Sr,0,Sz)
pml_layers = [mp.PML(pml_size)]

nf = 2
fcen = 1/1.55
width = 0.1
fwidth = width * fcen
source_center  = mp.Vector3(padl+design_region_width/2,0,-padz-design_region_height/2)
source_size    = mp.Vector3(design_region_width,0,0)
src = mp.GaussianSource(frequency=fcen,fwidth=fwidth)
source = [mp.Source(src,component=mp.Ep,center=source_center, size=source_size)]

design_region_resolution = 20
Nx = int(design_region_resolution*design_region_width)
Ny = int(design_region_resolution*design_region_height)
design_variables = mp.MaterialGrid(mp.Vector3(Nx,1, Ny),mp.air, Si, do_averaging=False)
design_region = mpa.DesignRegion(design_variables,volume=mp.Volume(center=mp.Vector3(padl+design_region_width/2, 0, 0), size=mp.Vector3(design_region_width, 0, design_region_height)))

geometry = [mp.Block(center=design_region.center, size=design_region.size, material=design_variables)] # design region
dimensions = mp.CYLINDRICAL
sim = mp.Simulation(cell_size=cell_size, boundary_layers=pml_layers, geometry=geometry, sources=source, default_material=mp.air,
                    resolution=resolution, dimensions=dimensions, m=-1, force_complex_fields=True)

far_x = [mp.Vector3(0.1,0,20)]
NearRegions = [mp.Near2FarRegion(center=mp.Vector3(padl+(design_region_width)/2,0,padz+design_region_height/2), 
                size=mp.Vector3(0.5,0,0), weight=+1)]
FarFields = mpa.Near2FarFields(sim, NearRegions ,far_x)
ob_list = [FarFields]

def J(ff):
    return npa.abs(ff[0,0,1])**2
obj_fs = [J]
opt = mpa.OptimizationProblem(
    simulation = sim,
    objective_functions = obj_fs,
    objective_arguments = ob_list,
    design_regions = [design_region],
    fcen = fcen,
    nf=1,
    df=0,
    minimum_run_time = 2000,
    maximum_run_time = 2001
)

x = 0.5*np.ones(Nx*Ny)
f0, adj = opt([x])
#fd, idx = opt.calculate_fd_gradient(num_gradients=4, db = 1e-4)

print(f0)
print("adj",adj)

The fields before PR:

fwd c er, fwd 0.645056+i0.153555, adj c er, adj 0.000002+i0.000011, (0.475000, 0.000000,0.000000, -0.050000) 
fwd c er, fwd 0.643162+i0.383013, adj c er, adj 0.000005+i0.000010, (0.475000, 0.000000,0.000000, 0.000000) 
fwd c er, fwd 0.602949+i0.602070, adj c er, adj 0.000008+i0.000007, (0.475000, 0.000000,0.000000, 0.050000) 
fwd c er, fwd 0.405817+i0.497857, adj c er, adj 0.000007+i0.000005, (0.475000, 0.000000,0.000000, 0.100000) 
fwd c er, fwd 0.651281+i-0.039765, adj c er, adj -0.000000+i0.000015, (0.525000, 0.000000,0.000000, -0.050000) 
fwd c er, fwd 0.638741+i0.398589, adj c er, adj 0.000005+i0.000010, (0.525000, 0.000000,0.000000, 0.000000) 
fwd c er, fwd 0.734200+i0.989445, adj c er, adj 0.000012+i0.000007, (0.525000, 0.000000,0.000000, 0.050000) 
fwd c er, fwd 0.230273+i0.344623, adj c er, adj 0.000005+i0.000003, (0.525000, 0.000000,0.000000, 0.100000) 
fwd c er, fwd 1.118304+i0.707820, adj c er, adj 0.000006+i0.000014, (0.575000, 0.000000,0.000000, -0.050000) 
fwd c er, fwd 0.652416+i0.453383, adj c er, adj 0.000006+i0.000012, (0.575000, 0.000000,0.000000, 0.000000) 
fwd c er, fwd 0.561694+i0.367975, adj c er, adj 0.000008+i0.000016, (0.575000, 0.000000,0.000000, 0.050000) 
fwd c er, fwd 0.188903+i0.293060, adj c er, adj 0.000005+i0.000004, (0.575000, 0.000000,0.000000, 0.100000) 
fwd c er, fwd 0.630064+i0.389666, adj c er, adj 0.000004+i0.000010, (0.625000, 0.000000,0.000000, -0.050000) 
fwd c er, fwd 0.528196+i0.387201, adj c er, adj 0.000005+i0.000010, (0.625000, 0.000000,0.000000, 0.000000) 
fwd c er, fwd 0.420290+i0.369146, adj c er, adj 0.000006+i0.000010, (0.625000, 0.000000,0.000000, 0.050000) 
fwd c er, fwd 0.256575+i0.360613, adj c er, adj 0.000006+i0.000006, (0.625000, 0.000000,0.000000, 0.100000) 
fwd c ep, fwd 0.092786+i-0.482891, adj c ep, adj -0.000010+i0.000002, (0.500000, 0.000000,0.000000, -0.050000) 
fwd c ep, fwd 0.221215+i-0.499613, adj c ep, adj -0.000010+i0.000004, (0.500000, 0.000000,0.000000, 0.000000) 
fwd c ep, fwd 0.332522+i-0.473199, adj c ep, adj -0.000008+i0.000005, (0.500000, 0.000000,0.000000, 0.050000) 
fwd c ep, fwd 0.378766+i-0.394397, adj c ep, adj -0.000006+i0.000006, (0.500000, 0.000000,0.000000, 0.100000) 
fwd c ep, fwd 0.375797+i-2.848804, adj c ep, adj -0.000069+i0.000011, (0.550000, 0.000000,0.000000, -0.050000) 
fwd c ep, fwd 1.271099+i-3.059732, adj c ep, adj -0.000065+i0.000024, (0.550000, 0.000000,0.000000, 0.000000) 
fwd c ep, fwd 2.042060+i-2.922385, adj c ep, adj -0.000056+i0.000034, (0.550000, 0.000000,0.000000, 0.050000) 
fwd c ep, fwd 0.335374+i-0.340894, adj c ep, adj -0.000006+i0.000006, (0.550000, 0.000000,0.000000, 0.100000) 
fwd c ep, fwd 0.011707+i-0.289568, adj c ep, adj -0.000008+i0.000001, (0.600000, 0.000000,0.000000, -0.050000) 
fwd c ep, fwd 0.114356+i-0.312847, adj c ep, adj -0.000007+i0.000002, (0.600000, 0.000000,0.000000, 0.000000) 
fwd c ep, fwd 0.207578+i-0.305733, adj c ep, adj -0.000006+i0.000004, (0.600000, 0.000000,0.000000, 0.050000) 
fwd c ep, fwd 0.279757+i-0.274860, adj c ep, adj -0.000005+i0.000005, (0.600000, 0.000000,0.000000, 0.100000) 
fwd c ep, fwd -0.002345+i-0.243960, adj c ep, adj -0.000007+i0.000001, (0.650000, 0.000000,0.000000, -0.050000) 
fwd c ep, fwd 0.085468+i-0.258471, adj c ep, adj -0.000007+i0.000002, (0.650000, 0.000000,0.000000, 0.000000) 
fwd c ep, fwd 0.166449+i-0.251492, adj c ep, adj -0.000006+i0.000003, (0.650000, 0.000000,0.000000, 0.050000) 
fwd c ep, fwd 0.234422+i-0.224610, adj c ep, adj -0.000005+i0.000004, (0.650000, 0.000000,0.000000, 0.100000) 
fwd c ez, fwd -0.112498+i-0.320694, adj c ez, adj -0.000003+i0.000006, (0.500000, 0.000000,0.000000, -0.075000) 
fwd c ez, fwd -0.135251+i-0.123785, adj c ez, adj -0.000000+i0.000001, (0.500000, 0.000000,0.000000, -0.025000) 
fwd c ez, fwd -0.144965+i-0.156319, adj c ez, adj -0.000001+i0.000001, (0.500000, 0.000000,0.000000, 0.025000) 
fwd c ez, fwd -0.295754+i-0.590018, adj c ez, adj -0.000005+i0.000002, (0.500000, 0.000000,0.000000, 0.075000) 
fwd c ez, fwd -0.372778+i-0.470543, adj c ez, adj -0.000004+i0.000002, (0.550000, 0.000000,0.000000, -0.075000) 
fwd c ez, fwd -0.661254+i-1.214330, adj c ez, adj -0.000010+i0.000009, (0.550000, 0.000000,0.000000, -0.025000) 
fwd c ez, fwd -0.455460+i-1.192297, adj c ez, adj -0.000009+i0.000011, (0.550000, 0.000000,0.000000, 0.025000) 
fwd c ez, fwd -0.076187+i-0.446885, adj c ez, adj -0.000003+i0.000006, (0.550000, 0.000000,0.000000, 0.075000) 
fwd c ez, fwd -0.539074+i-0.550588, adj c ez, adj -0.000005+i-0.000002, (0.600000, 0.000000,0.000000, -0.075000) 
fwd c ez, fwd -0.099551+i-0.277188, adj c ez, adj -0.000003+i0.000001, (0.600000, 0.000000,0.000000, -0.025000) 
fwd c ez, fwd 0.001548+i-0.236501, adj c ez, adj -0.000002+i0.000003, (0.600000, 0.000000,0.000000, 0.025000) 
fwd c ez, fwd 0.127514+i-0.251088, adj c ez, adj -0.000000+i0.000008, (0.600000, 0.000000,0.000000, 0.075000) 
fwd c ez, fwd -0.313957+i-0.418766, adj c ez, adj -0.000004+i0.000000, (0.650000, 0.000000,0.000000, -0.075000) 
fwd c ez, fwd -0.143470+i-0.340765, adj c ez, adj -0.000003+i0.000002, (0.650000, 0.000000,0.000000, -0.025000) 
fwd c ez, fwd -0.032632+i-0.303925, adj c ez, adj -0.000002+i0.000003, (0.650000, 0.000000,0.000000, 0.025000) 
fwd c ez, fwd 0.049530+i-0.298294, adj c ez, adj -0.000001+i0.000005, (0.650000, 0.000000,0.000000, 0.075000)

The fields after PR:

fwd c dr, fwd 0.645056+i0.153555, adj c dr, adj 0.000004+i-0.000002, (0.475000, 0.000000,0.000000, -0.050000) 
fwd c dr, fwd 0.643162+i0.383013, adj c dr, adj 0.000004+i-0.000001, (0.475000, 0.000000,0.000000, 0.000000) 
fwd c dr, fwd 0.602949+i0.602070, adj c dr, adj 0.000005+i-0.000001, (0.475000, 0.000000,0.000000, 0.050000) 
fwd c dr, fwd 0.405817+i0.497857, adj c dr, adj 0.000004+i-0.000001, (0.475000, 0.000000,0.000000, 0.100000) 
fwd c dr, fwd 0.651281+i-0.039765, adj c dr, adj 0.000004+i-0.000004, (0.525000, 0.000000,0.000000, -0.050000) 
fwd c dr, fwd 0.638741+i0.398589, adj c dr, adj 0.000005+i-0.000001, (0.525000, 0.000000,0.000000, 0.000000) 
fwd c dr, fwd 0.734200+i0.989445, adj c dr, adj 0.000009+i0.000001, (0.525000, 0.000000,0.000000, 0.050000) 
fwd c dr, fwd 0.230273+i0.344623, adj c dr, adj 0.000003+i-0.000001, (0.525000, 0.000000,0.000000, 0.100000) 
fwd c dr, fwd 1.118304+i0.707820, adj c dr, adj 0.000006+i0.000004, (0.575000, 0.000000,0.000000, -0.050000) 
fwd c dr, fwd 0.652416+i0.453383, adj c dr, adj 0.000005+i-0.000000, (0.575000, 0.000000,0.000000, 0.000000) 
fwd c dr, fwd 0.561694+i0.367975, adj c dr, adj 0.000004+i-0.000005, (0.575000, 0.000000,0.000000, 0.050000) 
fwd c dr, fwd 0.188903+i0.293060, adj c dr, adj 0.000003+i-0.000001, (0.575000, 0.000000,0.000000, 0.100000) 
fwd c dr, fwd 0.630064+i0.389666, adj c dr, adj 0.000004+i0.000002, (0.625000, 0.000000,0.000000, -0.050000) 
fwd c dr, fwd 0.528196+i0.387201, adj c dr, adj 0.000004+i-0.000000, (0.625000, 0.000000,0.000000, 0.000000) 
fwd c dr, fwd 0.420290+i0.369146, adj c dr, adj 0.000004+i-0.000002, (0.625000, 0.000000,0.000000, 0.050000) 
fwd c dr, fwd 0.256575+i0.360613, adj c dr, adj 0.000004+i-0.000001, (0.625000, 0.000000,0.000000, 0.100000) 
fwd c dp, fwd 0.092786+i-0.482891, adj c dp, adj 0.000000+i0.000002, (0.500000, 0.000000,0.000000, -0.050000) 
fwd c dp, fwd 0.221215+i-0.499613, adj c dp, adj 0.000001+i0.000002, (0.500000, 0.000000,0.000000, 0.000000) 
fwd c dp, fwd 0.332522+i-0.473199, adj c dp, adj 0.000001+i0.000001, (0.500000, 0.000000,0.000000, 0.050000) 
fwd c dp, fwd 0.378766+i-0.394397, adj c dp, adj 0.000002+i0.000002, (0.500000, 0.000000,0.000000, 0.100000) 
fwd c dp, fwd 0.375797+i-2.848804, adj c dp, adj 0.000005+i0.000017, (0.550000, 0.000000,0.000000, -0.050000) 
fwd c dp, fwd 1.271099+i-3.059732, adj c dp, adj 0.000010+i0.000017, (0.550000, 0.000000,0.000000, 0.000000) 
fwd c dp, fwd 2.042060+i-2.922385, adj c dp, adj 0.000014+i0.000016, (0.550000, 0.000000,0.000000, 0.050000) 
fwd c dp, fwd 0.335374+i-0.340894, adj c dp, adj 0.000002+i0.000002, (0.550000, 0.000000,0.000000, 0.100000) 
fwd c dp, fwd 0.011707+i-0.289568, adj c dp, adj 0.000001+i0.000003, (0.600000, 0.000000,0.000000, -0.050000) 
fwd c dp, fwd 0.114356+i-0.312847, adj c dp, adj 0.000001+i0.000003, (0.600000, 0.000000,0.000000, 0.000000) 
fwd c dp, fwd 0.207578+i-0.305733, adj c dp, adj 0.000002+i0.000003, (0.600000, 0.000000,0.000000, 0.050000) 
fwd c dp, fwd 0.279757+i-0.274860, adj c dp, adj 0.000002+i0.000003, (0.600000, 0.000000,0.000000, 0.100000) 
fwd c dp, fwd -0.002345+i-0.243960, adj c dp, adj 0.000001+i0.000003, (0.650000, 0.000000,0.000000, -0.050000) 
fwd c dp, fwd 0.085468+i-0.258471, adj c dp, adj 0.000002+i0.000003, (0.650000, 0.000000,0.000000, 0.000000) 
fwd c dp, fwd 0.166449+i-0.251492, adj c dp, adj 0.000002+i0.000003, (0.650000, 0.000000,0.000000, 0.050000) 
fwd c dp, fwd 0.234422+i-0.224610, adj c dp, adj 0.000003+i0.000003, (0.650000, 0.000000,0.000000, 0.100000) 
fwd c dz, fwd -0.112498+i-0.320694, adj c dz, adj 0.000000+i-0.000005, (0.500000, 0.000000,0.000000, -0.075000) 
fwd c dz, fwd -0.135251+i-0.123785, adj c dz, adj -0.000001+i-0.000002, (0.500000, 0.000000,0.000000, -0.025000) 
fwd c dz, fwd -0.144965+i-0.156319, adj c dz, adj -0.000002+i-0.000002, (0.500000, 0.000000,0.000000, 0.025000) 
fwd c dz, fwd -0.295754+i-0.590018, adj c dz, adj -0.000006+i-0.000004, (0.500000, 0.000000,0.000000, 0.075000) 
fwd c dz, fwd -0.372778+i-0.470543, adj c dz, adj -0.000001+i-0.000005, (0.550000, 0.000000,0.000000, -0.075000) 
fwd c dz, fwd -0.661254+i-1.214330, adj c dz, adj -0.000005+i-0.000013, (0.550000, 0.000000,0.000000, -0.025000) 
fwd c dz, fwd -0.455460+i-1.192297, adj c dz, adj -0.000007+i-0.000013, (0.550000, 0.000000,0.000000, 0.025000) 
fwd c dz, fwd -0.076187+i-0.446885, adj c dz, adj -0.000004+i-0.000005, (0.550000, 0.000000,0.000000, 0.075000) 
fwd c dz, fwd -0.539074+i-0.550588, adj c dz, adj -0.000002+i-0.000005, (0.600000, 0.000000,0.000000, -0.075000) 
fwd c dz, fwd -0.099551+i-0.277188, adj c dz, adj -0.000001+i-0.000002, (0.600000, 0.000000,0.000000, -0.025000) 
fwd c dz, fwd 0.001548+i-0.236501, adj c dz, adj -0.000001+i-0.000002, (0.600000, 0.000000,0.000000, 0.025000) 
fwd c dz, fwd 0.127514+i-0.251088, adj c dz, adj -0.000001+i-0.000005, (0.600000, 0.000000,0.000000, 0.075000) 
fwd c dz, fwd -0.313957+i-0.418766, adj c dz, adj -0.000000+i-0.000004, (0.650000, 0.000000,0.000000, -0.075000) 
fwd c dz, fwd -0.143470+i-0.340765, adj c dz, adj -0.000000+i-0.000003, (0.650000, 0.000000,0.000000, -0.025000) 
fwd c dz, fwd -0.032632+i-0.303925, adj c dz, adj -0.000001+i-0.000003, (0.650000, 0.000000,0.000000, 0.025000) 
fwd c dz, fwd 0.049530+i-0.298294, adj c dz, adj -0.000001+i-0.000004, (0.650000, 0.000000,0.000000, 0.075000)
smartalecH commented 2 years ago

In your master_printf() statement, can you use %e for the fields, rather than %f? That way it prints them in scientific notation (thus allowing for maximum digits of precision).

mochen4 commented 2 years ago 
printf("fwd c %s, fwd %e+i%e, adj c %s, adj %e+i%e, (%f, %f) \n", meep::component_name(forward_c), real(fwd), imag(fwd), meep::component_name(adjoint_c), real(adj), imag(adj), p.r(), p.z())

Here is the results with the above printf statement: Before:

fwd c er, fwd 6.450557e-01+i1.535550e-01, adj c er, adj 2.066160e-06+i1.116971e-05, (0.475000, -0.050000) 
fwd c er, fwd 6.431620e-01+i3.830131e-01, adj c er, adj 5.180422e-06+i9.543417e-06, (0.475000, 0.000000) 
fwd c er, fwd 6.029485e-01+i6.020699e-01, adj c er, adj 8.090047e-06+i7.381786e-06, (0.475000, 0.050000) 
fwd c er, fwd 4.058174e-01+i4.978575e-01, adj c er, adj 7.433247e-06+i5.209937e-06, (0.475000, 0.100000) 
fwd c er, fwd 6.512812e-01+i-3.976512e-02, adj c er, adj -1.294129e-07+i1.540280e-05, (0.525000, -0.050000) 
fwd c er, fwd 6.387413e-01+i3.985892e-01, adj c er, adj 5.281752e-06+i9.943139e-06, (0.525000, 0.000000) 
fwd c er, fwd 7.342002e-01+i9.894447e-01, adj c er, adj 1.214755e-05+i6.780654e-06, (0.525000, 0.050000) 
fwd c er, fwd 2.302731e-01+i3.446227e-01, adj c er, adj 5.329128e-06+i3.413410e-06, (0.525000, 0.100000) 
fwd c er, fwd 1.118304e+00+i7.078195e-01, adj c er, adj 6.246399e-06+i1.380167e-05, (0.575000, -0.050000) 
fwd c er, fwd 6.524161e-01+i4.533830e-01, adj c er, adj 5.761601e-06+i1.228389e-05, (0.575000, 0.000000) 
fwd c er, fwd 5.616938e-01+i3.679747e-01, adj c er, adj 8.247404e-06+i1.638200e-05, (0.575000, 0.050000) 
fwd c er, fwd 1.889034e-01+i2.930596e-01, adj c er, adj 4.787147e-06+i3.973302e-06, (0.575000, 0.100000) 
fwd c er, fwd 6.300638e-01+i3.896657e-01, adj c er, adj 3.585567e-06+i1.018880e-05, (0.625000, -0.050000) 
fwd c er, fwd 5.281961e-01+i3.872010e-01, adj c er, adj 4.764728e-06+i1.039714e-05, (0.625000, 0.000000) 
fwd c er, fwd 4.202898e-01+i3.691465e-01, adj c er, adj 5.890000e-06+i1.025350e-05, (0.625000, 0.050000) 
fwd c er, fwd 2.565745e-01+i3.606132e-01, adj c er, adj 5.717705e-06+i6.239493e-06, (0.625000, 0.100000) 
fwd c ep, fwd 9.278563e-02+i-4.828911e-01, adj c ep, adj -1.020829e-05+i2.073916e-06, (0.500000, -0.050000) 
fwd c ep, fwd 2.212145e-01+i-4.996127e-01, adj c ep, adj -9.601048e-06+i3.892076e-06, (0.500000, 0.000000) 
fwd c ep, fwd 3.325218e-01+i-4.731990e-01, adj c ep, adj -8.277467e-06+i5.442711e-06, (0.500000, 0.050000) 
fwd c ep, fwd 3.787657e-01+i-3.943971e-01, adj c ep, adj -6.472010e-06+i6.195760e-06, (0.500000, 0.100000) 
fwd c ep, fwd 3.757966e-01+i-2.848804e+00, adj c ep, adj -6.873723e-05+i1.146188e-05, (0.550000, -0.050000) 
fwd c ep, fwd 1.271099e+00+i-3.059732e+00, adj c ep, adj -6.526817e-05+i2.388106e-05, (0.550000, 0.000000) 
fwd c ep, fwd 2.042060e+00+i-2.922385e+00, adj c ep, adj -5.582382e-05+i3.418362e-05, (0.550000, 0.050000) 
fwd c ep, fwd 3.353738e-01+i-3.408944e-01, adj c ep, adj -6.132979e-06+i5.579830e-06, (0.550000, 0.100000) 
fwd c ep, fwd 1.170709e-02+i-2.895680e-01, adj c ep, adj -7.915985e-06+i9.528919e-07, (0.600000, -0.050000) 
fwd c ep, fwd 1.143558e-01+i-3.128470e-01, adj c ep, adj -7.444631e-06+i2.365050e-06, (0.600000, 0.000000) 
fwd c ep, fwd 2.075779e-01+i-3.057332e-01, adj c ep, adj -6.487342e-06+i3.624752e-06, (0.600000, 0.050000) 
fwd c ep, fwd 2.797569e-01+i-2.748600e-01, adj c ep, adj -5.485223e-06+i4.749444e-06, (0.600000, 0.100000) 
fwd c ep, fwd -2.345337e-03+i-2.439603e-01, adj c ep, adj -7.118196e-06+i5.744760e-07, (0.650000, -0.050000) 
fwd c ep, fwd 8.546788e-02+i-2.584715e-01, adj c ep, adj -6.717534e-06+i1.829815e-06, (0.650000, 0.000000) 
fwd c ep, fwd 1.664485e-01+i-2.514918e-01, adj c ep, adj -5.988234e-06+i2.998626e-06, (0.650000, 0.050000) 
fwd c ep, fwd 2.344216e-01+i-2.246099e-01, adj c ep, adj -5.034933e-06+i4.038350e-06, (0.650000, 0.100000) 
fwd c ez, fwd -1.124977e-01+i-3.206942e-01, adj c ez, adj -2.775077e-06+i5.874407e-06, (0.500000, -0.075000) 
fwd c ez, fwd -1.352510e-01+i-1.237850e-01, adj c ez, adj -4.689497e-07+i1.333520e-06, (0.500000, -0.025000) 
fwd c ez, fwd -1.449648e-01+i-1.563195e-01, adj c ez, adj -7.041813e-07+i9.195791e-07, (0.500000, 0.025000) 
fwd c ez, fwd -2.957536e-01+i-5.900177e-01, adj c ez, adj -5.229288e-06+i1.640336e-06, (0.500000, 0.075000) 
fwd c ez, fwd -3.727776e-01+i-4.705428e-01, adj c ez, adj -4.346507e-06+i1.991257e-06, (0.550000, -0.075000) 
fwd c ez, fwd -6.612543e-01+i-1.214330e+00, adj c ez, adj -9.958379e-06+i8.513755e-06, (0.550000, -0.025000) 
fwd c ez, fwd -4.554603e-01+i-1.192297e+00, adj c ez, adj -8.769191e-06+i1.109616e-05, (0.550000, 0.025000) 
fwd c ez, fwd -7.618686e-02+i-4.468849e-01, adj c ez, adj -2.688492e-06+i5.516855e-06, (0.550000, 0.075000) 
fwd c ez, fwd -5.390735e-01+i-5.505883e-01, adj c ez, adj -5.173560e-06+i-1.949995e-06, (0.600000, -0.075000) 
fwd c ez, fwd -9.955087e-02+i-2.771879e-01, adj c ez, adj -2.842984e-06+i1.322933e-06, (0.600000, -0.025000) 
fwd c ez, fwd 1.547694e-03+i-2.365006e-01, adj c ez, adj -2.087622e-06+i2.885026e-06, (0.600000, 0.025000) 
fwd c ez, fwd 1.275135e-01+i-2.510875e-01, adj c ez, adj -1.721241e-08+i8.444322e-06, (0.600000, 0.075000) 
fwd c ez, fwd -3.139573e-01+i-4.187660e-01, adj c ez, adj -3.878146e-06+i3.632693e-07, (0.650000, -0.075000) 
fwd c ez, fwd -1.434703e-01+i-3.407647e-01, adj c ez, adj -2.926732e-06+i1.823297e-06, (0.650000, -0.025000) 
fwd c ez, fwd -3.263224e-02+i-3.039245e-01, adj c ez, adj -2.089580e-06+i3.380786e-06, (0.650000, 0.025000) 
fwd c ez, fwd 4.952949e-02+i-2.982943e-01, adj c ez, adj -1.141573e-06+i5.424936e-06, (0.650000, 0.075000)

After:

fwd c dr, fwd 6.450557e-01+i1.535550e-01, adj c dr, adj 3.742905e-06+i-1.983668e-06, (4.750000e-01, -5.000000e-02) 
fwd c dr, fwd 6.431620e-01+i3.830131e-01, adj c dr, adj 4.439404e-06+i-1.339612e-06, (4.750000e-01, 0.000000e+00) 
fwd c dr, fwd 6.029485e-01+i6.020699e-01, adj c dr, adj 5.375141e-06+i-7.010624e-07, (4.750000e-01, 5.000000e-02) 
fwd c dr, fwd 4.058174e-01+i4.978575e-01, adj c dr, adj 3.884652e-06+i-1.367128e-06, (4.750000e-01, 1.000000e-01) 
fwd c dr, fwd 6.512812e-01+i-3.976512e-02, adj c dr, adj 4.410885e-06+i-4.106187e-06, (5.250000e-01, -5.000000e-02) 
fwd c dr, fwd 6.387413e-01+i3.985892e-01, adj c dr, adj 4.739418e-06+i-1.192076e-06, (5.250000e-01, 0.000000e+00) 
fwd c dr, fwd 7.342002e-01+i9.894447e-01, adj c dr, adj 8.781568e-06+i1.211187e-06, (5.250000e-01, 5.000000e-02) 
fwd c dr, fwd 2.302731e-01+i3.446227e-01, adj c dr, adj 2.828184e-06+i-1.012647e-06, (5.250000e-01, 1.000000e-01) 
fwd c dr, fwd 1.118304e+00+i7.078195e-01, adj c dr, adj 6.344741e-06+i4.432403e-06, (5.750000e-01, -5.000000e-02) 
fwd c dr, fwd 6.524161e-01+i4.533830e-01, adj c dr, adj 4.835983e-06+i-2.795558e-07, (5.750000e-01, 0.000000e+00) 
fwd c dr, fwd 5.616938e-01+i3.679747e-01, adj c dr, adj 3.881450e-06+i-5.397551e-06, (5.750000e-01, 5.000000e-02) 
fwd c dr, fwd 1.889034e-01+i2.930596e-01, adj c dr, adj 2.682128e-06+i-1.183209e-06, (5.750000e-01, 1.000000e-01) 
fwd c dr, fwd 6.300638e-01+i3.896657e-01, adj c dr, adj 4.374641e-06+i1.568593e-06, (6.250000e-01, -5.000000e-02) 
fwd c dr, fwd 5.281961e-01+i3.872010e-01, adj c dr, adj 4.298628e-06+i-9.556317e-08, (6.250000e-01, 0.000000e+00) 
fwd c dr, fwd 4.202898e-01+i3.691465e-01, adj c dr, adj 3.976348e-06+i-1.841314e-06, (6.250000e-01, 5.000000e-02) 
fwd c dr, fwd 2.565745e-01+i3.606132e-01, adj c dr, adj 3.680750e-06+i-1.355337e-06, (6.250000e-01, 1.000000e-01) 
fwd c dp, fwd 9.278563e-02+i-4.828911e-01, adj c dp, adj 4.373042e-07+i1.529266e-06, (5.000000e-01, -5.000000e-02) 
fwd c dp, fwd 2.212145e-01+i-4.996127e-01, adj c dp, adj 8.957709e-07+i1.540579e-06, (5.000000e-01, 0.000000e+00) 
fwd c dp, fwd 3.325218e-01+i-4.731990e-01, adj c dp, adj 1.353846e-06+i1.478889e-06, (5.000000e-01, 5.000000e-02) 
fwd c dp, fwd 3.787657e-01+i-3.943971e-01, adj c dp, adj 1.598014e-06+i1.691261e-06, (5.000000e-01, 1.000000e-01) 
fwd c dp, fwd 3.757966e-01+i-2.848804e+00, adj c dp, adj 5.213380e-06+i1.660870e-05, (5.500000e-01, -5.000000e-02) 
fwd c dp, fwd 1.271099e+00+i-3.059732e+00, adj c dp, adj 9.614529e-06+i1.706159e-05, (5.500000e-01, 0.000000e+00) 
fwd c dp, fwd 2.042060e+00+i-2.922385e+00, adj c dp, adj 1.380111e-05+i1.615548e-05, (5.500000e-01, 5.000000e-02) 
fwd c dp, fwd 3.353738e-01+i-3.408944e-01, adj c dp, adj 2.120063e-06+i2.203443e-06, (5.500000e-01, 1.000000e-01) 
fwd c dp, fwd 1.170709e-02+i-2.895680e-01, adj c dp, adj 8.922617e-07+i2.739605e-06, (6.000000e-01, -5.000000e-02) 
fwd c dp, fwd 1.143558e-01+i-3.128470e-01, adj c dp, adj 1.484191e-06+i2.820291e-06, (6.000000e-01, 0.000000e+00) 
fwd c dp, fwd 2.075779e-01+i-3.057332e-01, adj c dp, adj 2.079399e-06+i2.754291e-06, (6.000000e-01, 5.000000e-02) 
fwd c dp, fwd 2.797569e-01+i-2.748600e-01, adj c dp, adj 2.402376e-06+i2.600257e-06, (6.000000e-01, 1.000000e-01) 
fwd c dp, fwd -2.345337e-03+i-2.439603e-01, adj c dp, adj 1.026560e-06+i2.859971e-06, (6.500000e-01, -5.000000e-02) 
fwd c dp, fwd 8.546788e-02+i-2.584715e-01, adj c dp, adj 1.526597e-06+i2.954450e-06, (6.500000e-01, 0.000000e+00) 
fwd c dp, fwd 1.664485e-01+i-2.514918e-01, adj c dp, adj 2.048087e-06+i2.951919e-06, (6.500000e-01, 5.000000e-02) 
fwd c dp, fwd 2.344216e-01+i-2.246099e-01, adj c dp, adj 2.525493e-06+i2.882366e-06, (6.500000e-01, 1.000000e-01) 
fwd c dz, fwd -1.124977e-01+i-3.206942e-01, adj c dz, adj 4.122350e-07+i-4.952361e-06, (5.000000e-01, -7.500000e-02) 
fwd c dz, fwd -1.352510e-01+i-1.237850e-01, adj c dz, adj -8.163604e-07+i-2.481619e-06, (5.000000e-01, -2.500000e-02) 
fwd c dz, fwd -1.449648e-01+i-1.563195e-01, adj c dz, adj -1.729374e-06+i-2.412992e-06, (5.000000e-01, 2.500000e-02) 
fwd c dz, fwd -2.957536e-01+i-5.900177e-01, adj c dz, adj -5.991521e-06+i-4.215360e-06, (5.000000e-01, 7.500000e-02) 
fwd c dz, fwd -3.727776e-01+i-4.705428e-01, adj c dz, adj -8.314765e-07+i-5.089111e-06, (5.500000e-01, -7.500000e-02) 
fwd c dz, fwd -6.612543e-01+i-1.214330e+00, adj c dz, adj -4.764110e-06+i-1.316858e-05, (5.500000e-01, -2.500000e-02) 
fwd c dz, fwd -4.554603e-01+i-1.192297e+00, adj c dz, adj -6.846978e-06+i-1.314016e-05, (5.500000e-01, 2.500000e-02) 
fwd c dz, fwd -7.618686e-02+i-4.468849e-01, adj c dz, adj -3.991132e-06+i-5.086486e-06, (5.500000e-01, 7.500000e-02) 
fwd c dz, fwd -5.390735e-01+i-5.505883e-01, adj c dz, adj -1.922884e-06+i-4.626941e-06, (6.000000e-01, -7.500000e-02) 
fwd c dz, fwd -9.955087e-02+i-2.771879e-01, adj c dz, adj -6.277213e-07+i-1.938816e-06, (6.000000e-01, -2.500000e-02) 
fwd c dz, fwd 1.547694e-03+i-2.365006e-01, adj c dz, adj -7.059988e-07+i-1.983497e-06, (6.000000e-01, 2.500000e-02) 
fwd c dz, fwd 1.275135e-01+i-2.510875e-01, adj c dz, adj -1.357829e-06+i-5.064835e-06, (6.000000e-01, 7.500000e-02) 
fwd c dz, fwd -3.139573e-01+i-4.187660e-01, adj c dz, adj -4.985871e-07+i-3.646495e-06, (6.500000e-01, -7.500000e-02) 
fwd c dz, fwd -1.434703e-01+i-3.407647e-01, adj c dz, adj -4.297963e-07+i-3.020549e-06, (6.500000e-01, -2.500000e-02) 
fwd c dz, fwd -3.263224e-02+i-3.039245e-01, adj c dz, adj -6.346266e-07+i-3.137108e-06, (6.500000e-01, 2.500000e-02) 
fwd c dz, fwd 4.952949e-02+i-2.982943e-01, adj c dz, adj -1.156082e-06+i-4.039807e-06, (6.500000e-01, 7.500000e-02)
smartalecH commented 2 years ago

Hmm is it possible the m number is not being properly flipped during the adjoint run?

I remember for #1855 I had to modify the way I flipped the k_point since we no longer need to rebuild the simulation:

https://github.com/NanoComp/meep/blob/ea43f83c2f8c3e74aab2f8dae3e2ed96b7347777/python/adjoint/optimization_problem.py#L233

I'm not sure the current code to change the m number actually changes it (since we don't reinitialize)

https://github.com/NanoComp/meep/blob/ea43f83c2f8c3e74aab2f8dae3e2ed96b7347777/python/adjoint/optimization_problem.py#L263-L264

(Note that we still have the following line... which I don't think does anything either... if it does we need to remove it!)

https://github.com/NanoComp/meep/blob/ea43f83c2f8c3e74aab2f8dae3e2ed96b7347777/python/adjoint/optimization_problem.py#L267-L268

EDIT: I'm guessing this is the issue, as my gradients for a simulation with m=0 match the finite differences really well. I'll see if I can properly update the m number.

mochen4 commented 2 years ago

EDIT: I'm guessing this is the issue, as my gradients for a simulation with m=0 match the finite differences really well. I'll see if I can properly update the m number.

I don't think this is the main issue, at least not the only issue. I ran the test with m=0 and the gradients didn't match. You also mentioned there were gradients mismatch in Cartesian cases? Perhaps we should debug m=0 or 2D Cartesian first?

(Note that we still have the following line... which I don't think does anything either... if it does we need to remove it!) https://github.com/NanoComp/meep/blob/ea43f83c2f8c3e74aab2f8dae3e2ed96b7347777/python/adjoint/optimization_problem.py#L267-L268

The purposes were to reset the parameter to the value in the forward simulation, so that we could subsequently run calculate_fd_gradient. Perhaps change it to something like self.sim.change_k_point(-1 * self.sim.k_point). Similarly for m values .

mochen4 commented 2 years ago

There is a known issue introduced in an earlier PR that the factor should be p.r() instead of 2 * p.r() here: https://github.com/NanoComp/meep/blob/ea43f83c2f8c3e74aab2f8dae3e2ed96b7347777/src/meepgeom.cpp#L2923

Before #1855, the adjoint gradients would be off by a factor of 2.

smartalecH commented 2 years ago

I slightly modified your example:

import meep as mp
import meep.adjoint as mpa
import numpy as np
from autograd import numpy as npa
from matplotlib import pyplot as plt

np.random.rand(240)

mp.verbosity(0)
resolution = 20
Si = mp.Medium(epsilon=13)

design_region_width = 1.0
design_region_height = 1.0
pml_size = 1.0
padl, padr, padz = 0.5 , 0.5, 0.5
Sr = design_region_width + pml_size + padl + padr
Sz = 2*pml_size + design_region_height + 2*padz + 1
cell_size = mp.Vector3(Sr,0,Sz)
pml_layers = [mp.PML(pml_size)]

nf = 2
fcen = 1/1.55
width = 0.1
fwidth = width * fcen
source_center  = mp.Vector3(padl+design_region_width/2,0,-padz-design_region_height/2)
source_size    = mp.Vector3(design_region_width,0,0)
src = mp.GaussianSource(frequency=fcen,fwidth=fwidth)
source = [mp.Source(src,component=mp.Ep,center=source_center, size=source_size)]

design_region_resolution = 20
Nx = int(design_region_resolution*design_region_width)
Ny = int(design_region_resolution*design_region_height)
design_variables = mp.MaterialGrid(mp.Vector3(Nx,1, Ny),mp.air, Si, do_averaging=False)
design_region = mpa.DesignRegion(design_variables,volume=mp.Volume(center=mp.Vector3(padl+design_region_width/2, 0, 0), size=mp.Vector3(design_region_width, 0, design_region_height)))

geometry = [mp.Block(center=design_region.center, size=design_region.size, material=design_variables)] # design region
dimensions = mp.CYLINDRICAL
sim = mp.Simulation(cell_size=cell_size, boundary_layers=pml_layers, geometry=geometry, sources=source, default_material=mp.air,
                    resolution=resolution, dimensions=dimensions, m=0, force_complex_fields=True)

far_x = [mp.Vector3(0.1,0,20)]
NearRegions = [mp.Near2FarRegion(center=mp.Vector3(padl+(design_region_width)/2,0,padz+design_region_height/2), 
                size=mp.Vector3(0.5,0,0), weight=+1)]
FarFields = mpa.Near2FarFields(sim, NearRegions ,far_x)
ob_list = [FarFields]

def J(ff):
    return npa.abs(ff[0,0,1])**2
obj_fs = [J]
opt = mpa.OptimizationProblem(
    simulation = sim,
    objective_functions = obj_fs,
    objective_arguments = ob_list,
    design_regions = [design_region],
    frequencies = [fcen]
)

x = 0.5*np.ones(Nx*Ny)

## ensure reproducible results
np.random.seed(314159)
deps = 1e-5
dp = deps*np.random.rand(Nx*Ny)

# adjoint simulation
f0, adj = opt([x])
adjoint_dd = (dp@adj).flatten() # adjoint directional derivative

# finite difference approximation
opt.update_design([x+dp/2])
f_pp, _ = opt(need_gradient=False)
opt.update_design([x-dp/2])
f_pm, _ = opt(need_gradient=False)
fd_grad = f_pp-f_pm

# errors
abs_error = np.abs(adjoint_dd-fd_grad)
rel_error = abs_error / np.abs(fd_grad)

print("Directional derivative -- adjoint solver: {}, FD: {}|| rel_error = {} || abs_error = {}".format(adjoint_dd,fd_grad,rel_error,abs_error))

Using #2194 and m=0, I get the following results (I have not yet changed the factor of 2):

Directional derivative -- adjoint solver: [1.89824615e-13], FD: 1.898245466881937e-13|| rel_error = [3.61207742e-07] || abs_error = [6.85660959e-20]

Which is pretty accurate... I modified some of the above parameters (e.g. the size of the design region) just to stress test it a bit. It seems good so far, but there are definitely more tests we need to do.

However, this doesn't work if $m\neq0$. For example, for m=-1 (your original example) I get the following:

Directional derivative -- adjoint solver: [2.34424655e-10], FD: 2.7949400497988275e-11|| rel_error = [7.3874663] || abs_error = [2.06475254e-10]

As I mention in #2194, I get the same result even if I don't change the m number before the adjoint run... I also have to force a max runtime (maximum_run_time=2000), as we've seen before. I bet the lack of convergence is related to #2182...

You also mentioned there were gradients mismatch in Cartesian cases?

I spent more time looking into this. The Cartesian cases actually seem rather robust. I still need to look at the cases introduced by @mawc2019... but I have a feeling those are more "edge" cases.

Note that #2194 modified a few lines of the recombination step which seems to help with accuracy.

smartalecH commented 2 years ago

OK @mochen4 my latest commit in #2194 seems to have fixed it. I can get gradients within about 2% (sometimes much better) using the example above with various m numbers (try changing the size of the design region too).