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

Incorrect adjoint fields in cylindrical coordinates #2193

Closed
mochen4 opened this issue Aug 11, 2022 · 7 comments · Fixed by #2194
Closed

Incorrect adjoint fields in cylindrical coordinates #2193

mochen4 opened this issue Aug 11, 2022 · 7 comments · Fixed by #2194

Comments

@mochen4
Copy link
Collaborator

mochen4 commented Aug 11, 2022

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:

meep/src/meepgeom.cpp

Lines 2925 to 2927 in 2aa9164

material_grids_addgradient_point(v_local + ng * f_i, vec_to_vector3(p),
scalegrad * cyl_scale, geps, adjoint_c, forward_c,
fwd, adj, frequencies[f_i], gv, du);

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

smartalecH commented Aug 11, 2022

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

@smartalecH smartalecH added the bug label Aug 11, 2022
@mochen4
Copy link
Collaborator Author

mochen4 commented Aug 11, 2022

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

smartalecH commented Aug 12, 2022

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:

self.sim.change_k_point(-1 * self.sim.k_point)

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

if utils._check_if_cylindrical(self.sim):
self.sim.m = -self.sim.m

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

if self.sim.k_point:
self.sim.k_point *= -1


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
Copy link
Collaborator Author

mochen4 commented Aug 12, 2022

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

if self.sim.k_point:
self.sim.k_point *= -1

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
Copy link
Collaborator Author

mochen4 commented Aug 12, 2022

There is a known issue introduced in an earlier PR that the factor should be p.r() instead of 2 * p.r() here:

cyl_scale = (gv.dim == meep::Dcyl) ? 2 * p.r()

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

@smartalecH
Copy link
Collaborator

smartalecH commented Aug 12, 2022

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 smartalecH changed the title Incorrect adjoint fields Incorrect adjoint fields in cylindrical coordinates Aug 12, 2022
@smartalecH
Copy link
Collaborator

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

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

Successfully merging a pull request may close this issue.

2 participants