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

WIP: loop_in_chunks over owned points #1895

Merged
merged 161 commits into from
Mar 10, 2022
Merged

Conversation

mochen4
Copy link
Collaborator

@mochen4 mochen4 commented Jan 12, 2022

Fixes #1555.

@review-notebook-app
Copy link

Check out this pull request on  ReviewNB

See visual diffs & provide feedback on Jupyter Notebooks.


Powered by ReviewNB

@smartalecH
Copy link
Collaborator

You might want to rebase (or even pull and force push).

@smartalecH
Copy link
Collaborator

FYI I tried merging this branch with #1855 and ran a simple adjoint test case both with and without mpi (no symmetries). The fields blew up in both cases.

@mochen4
Copy link
Collaborator Author

mochen4 commented Jan 24, 2022

I am still having some difficulties. For concreteness, let's focus on the first test case in symmetry.cpp: test_1d_periodic_mirror. The grid volume is from z=0 to z=32, with Z mirror symmetry (mirror plane at z=16) and periodic BC. The where volume is 0 to 32 for Ex and Ey, -1 to 33 for Ez.

To make sure we are on the same page, here is the situation:
Without using symmertry:

  1. Ex and Ey chunks are 2 to 32. The overlap is 2 to 32. Then the periodic shift changes z=32 to 0, which finds another correct overlap.
  2. Ez chunks are 1 to 31. The overlaps is 1 to 31. Then one shift correctly overlaps 31 with -1, the other correctly overlaps 1 with 33.

Using symmetry:

  1. Ex and Ey chunks are 2 to 16, then the symmetry transform is 16 to 30. The transform has z=16 double counted. In addition, boundary points aren't captured by periodic shifts, which don't produce any overlaps.
  2. Ez chunk is 1 to 15, its transform is 17 to 31. Then one shift correctly overlaps 31 with -1, the other correctly overlaps 1 with 33. The result is correct. Just like the case without symmetry.

Issue:
The Ex and Ey chunks starts at z=2 instead of 0 isn't because periodic boundary, but because the little_owned_corner doesn't include the boundary by definition; on the other hand, big_owned_corner includes the boundary.

In grid_volume::halve, if I keep the other half 16 to 32 of the grid volume, the owned ivecs under consideration in loop_in_chunks becomes 18 to 32. The symmetric transform is 0 to 14. z=16 won't be captured, while the boundary points are "doubly double counted": both z=0 and z=32 will be counted twice, but z=0 and 32 are equivalent. In the correct scenario without symmetry, only z=32 is counted twice, z=0 is not included.

If I change the start location and size of the halved volume, either the boundaries are still "doubly double counted", or not counted at all. Moreover, the halved chunk isn't really half the original chunk, which may be potentially problematic.

@mochen4
Copy link
Collaborator Author

mochen4 commented Jan 24, 2022

By the way, using metallic BC in the same test case, without using symmetry, the overlap will be 2 to 32. Although both chunk and where is from 0 to 32, z=0 is not owned by the chunk, and there is no periodic shifts that overlaps 32 to 0. I am not sure this is what we want?

@stevengj
Copy link
Collaborator

stevengj commented Jan 26, 2022

As discussed: when you nontrivially (sn != 0) transform chunks[i]->gv, but before you shift it, then in any flipped direction you should remove the mirror plane and you should remove 0 if that direction is periodic

@stevengj
Copy link
Collaborator

As discussed today, you only want to remove the symmetry-plane overlap from a transformed chunk in a direction where that transformed gv is the "lower half" of the original user_volume. For any halved direction, this is determined by whether the is (lower corner) corresponds to where the upper half starts in grid_volume::halve.

@mochen4
Copy link
Collaborator Author

mochen4 commented Feb 16, 2022

It is almost working now. Here are the current issues when running CI:

  1. With double-precision (py3.7 MPI False, py3.7 MPI True, py3.10 MPI False): all c++ tests passed. Two python tests failed: test_adjoint_solver.py and test_dump_load.py.
    1.1 For test_adjoint_solver.py, sim.get_eigenmode_coefficients returns NaN (in the forward run). It is weird that many other tests also call sim.get_eigenmode_coefficients and they worked well.
    1.2 For test_dump_load.py, the error message from the log is:
HDF5-DIAG: Error detected in HDF5 (1.10.4) thread 140033161258816:
  #000: ../../../src/H5Dio.c line 336 in H5Dwrite(): can't write data
    major: Dataset
    minor: Write failed
  #001: ../../../src/H5Dio.c line 736 in H5D__write(): no output buffer
    major: Invalid arguments to routine
    minor: Bad value
munmap_chunk(): invalid pointer
  1. With single-precision (py3.10 MPI True), two C++ tests failed: aniso_disp and h5test. CI didn't generate logs so I tried them locally.
    2.1 For aniso_disp, the real part is good (0.4156 vs 0.41562), but the imaginary part is -6.10321e-07 rather than 4.8297e-7.
    2.2 h5test passed locally, so I don't know why it failed on CI.

@mochen4
Copy link
Collaborator Author

mochen4 commented Feb 24, 2022

@mawc2019 if you still have the original code of #1555, you may try whether this patch fixes the issue when you have a chance.

@mawc2019
Copy link
Contributor

mawc2019 commented Mar 2, 2022

I test your branch using an example almost the same as that in #1555. The results show that the issue still exists and is even more obvious. The code is as follows, which is based on your example code.

import meep as mp
import meep.adjoint as mpa
import numpy as np
from autograd import numpy as npa
from autograd import tensor_jacobian_product, grad

from scipy import special, signal
Si = mp.Medium(index=3.4)
SiO2 = mp.Medium(index=1.44)

design_region_width = 6
design_region_height = 2

pml_size = 1.0
resolution = 30

Sx = design_region_width
Sy = 2*pml_size + design_region_height + 5
cell_size = mp.Vector3(Sx,Sy)

nf = 2
frequencies = np.array([1/1.52, 1/1.58])
design_region_resolution = int(1*resolution)
pml_layers = [mp.PML(pml_size, direction=mp.Y)]

fcen = 1/1.55
width = 0.2
fwidth = width * fcen
source_center  = [0,-(design_region_height/2 + 1.5),0]
source_size    = mp.Vector3(Sx,0,0)
src = mp.GaussianSource(frequency=fcen,fwidth=fwidth)
source = [mp.EigenModeSource(src,eig_band = 1,size = source_size,center=source_center)]

Nx = int(design_region_resolution*design_region_width)
Ny = int(design_region_resolution*design_region_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_region_width, design_region_height, 0)))

geometry = [
    mp.Block(center=design_region.center, size=design_region.size, material=design_variables),
    mp.Block(center=design_region.center, size=design_region.size, material=design_variables, e1=mp.Vector3(x=-1))] # design region
kpoint = mp.Vector3()
sim = mp.Simulation(cell_size=cell_size,boundary_layers=pml_layers,k_point=kpoint,geometry=geometry,
                    sources=source,default_material=SiO2,symmetries=[mp.Mirror(direction=mp.X)],resolution=resolution)

far_x = [mp.Vector3(0,40,0)]
NearRegions = [mp.Near2FarRegion(center=mp.Vector3(0,design_region_height/2+1.5), size=mp.Vector3(design_region_width+1.5,0), weight=+1)]
FarFields = mpa.Near2FarFields(sim, NearRegions ,far_x)
ob_list = [FarFields]

def J(alpha):
    return npa.sum(npa.abs(alpha[0,:,2])**2)

opt = mpa.OptimizationProblem(
    simulation = sim,
    objective_functions = J,
    objective_arguments = ob_list,
    design_regions = [design_region],
    frequencies=frequencies,
    decay_by = 1e-6,
    maximum_run_time = 2000
)

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

g_adjoint = np.sum(dJ_deps, axis=1).flatten()
np.savetxt("adjoint_far_0302.out",g_adjoint)

I am sorry that I did not test it earlier. My only computer broke down and was left for repair in the last few days.

src/loop_in_chunks.cpp Outdated Show resolved Hide resolved
src/meep/vec.hpp Outdated Show resolved Hide resolved
src/meep/vec.hpp Outdated Show resolved Hide resolved
src/meep/vec.hpp Outdated Show resolved Hide resolved
@mochen4
Copy link
Collaborator Author

mochen4 commented Mar 8, 2022

@mawc2019 There seems to be an issue with overlapping design regions #1984. The adjoint gradients in the code are simply incorrect. The gradients are correct after removing this linemp.Block(center=design_region.center, size=design_region.size, material=design_variables, e1=mp.Vector3(x=-1)). And my testing showed that in that case, the issue in #1555 seems to be fixed.

@stevengj stevengj merged commit 26737ab into NanoComp:master Mar 10, 2022
@mawc2019
Copy link
Contributor

This commit sometimes causes incorrect adjoint gradients in the FourierFields adjoint solver. Using the example code below, the master branch before this commit yields

Finite-difference gradient:  [-119.00431169 -347.72326993  -78.66792763    7.16398801 -358.82263881]
Adjoint gradient:  [-119.00437279 -347.72554719  -78.66516035    7.16401012 -358.82190825]

but after this commit yields

Finite-difference gradient:  [-119.00431171 -347.72326992  -78.66792763    7.16398801 -358.82263881]
Adjoint gradient:  [ -93.28070124 -390.82490261  -76.1459804    17.70988613 -372.10113594]
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

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,yee_grid=True)
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])

@mochen4
Copy link
Collaborator Author

mochen4 commented Mar 14, 2022

This commit sometimes causes incorrect adjoint gradients in the FourierFields adjoint solver. Using the example code below, the master branch before this commit yields

Finite-difference gradient:  [-119.00431169 -347.72326993  -78.66792763    7.16398801 -358.82263881]
Adjoint gradient:  [-119.00437279 -347.72554719  -78.66516035    7.16401012 -358.82190825]

but after this commit yields

Finite-difference gradient:  [-119.00431171 -347.72326992  -78.66792763    7.16398801 -358.82263881]
Adjoint gradient:  [ -93.28070124 -390.82490261  -76.1459804    17.70988613 -372.10113594]

I didn't get same results with your exact code, but the gradients are indeed different. And in my case, the finite-difference gradients also changed somehow.
before:

Finite-difference gradient:  [-107.421875 -361.328125  -87.890625   19.53125  -351.5625  ]
Adjoint gradient:  [-119.00410995 -347.72466942  -78.66489225    7.16400295 -358.82102172]

after:

Finite-difference gradient:  [ -97.65625  -361.328125  -87.890625    0.       -351.5625  ]
Adjoint gradient:  [ -93.28038962 -390.8236593   -76.14582934   17.7098952  -372.10023175]

However, I personally don't think there is an issue with this PR. For near2far, the accuracy is actually improved overall (the only outlier is when the value itself is close to 0):
before:

Finite-difference gradient:  [-6.08377480e-04  2.94305584e-03  8.40114568e-03 -1.19700100e-05
 -7.30430429e-04]
Adjoint gradient:  [-5.70741537e-04  2.49842873e-03  7.46726086e-03 -6.45785975e-06
 -6.90577202e-04]

after:

Finite-difference gradient:  [-6.08377480e-04  2.94305584e-03  8.40114568e-03 -1.19700099e-05
 -7.30430429e-04]
Adjoint gradient:  [-5.98187898e-04  2.94184305e-03  8.41811757e-03 -9.69859556e-07
 -7.30216065e-04]

@mawc2019
Copy link
Contributor

mawc2019 commented Mar 15, 2022

I am using Python 3.7. Which version of Python are you using? Did you run the code on your own computer or Supercloud?

@mochen4
Copy link
Collaborator Author

mochen4 commented Mar 15, 2022

I am using Python 3.7. Which version of Python are you using? Did you run the code on your own computer or Supercloud?

Python 3.9, run locally on my laptop.

@mochen4
Copy link
Collaborator Author

mochen4 commented Mar 28, 2022

I ran a few more tests and plotted the the finite-difference gradients vs. the adjoint gradients (on some randomly selected pixels). When the design region does not touch the cell boundaries or when PML is added to all directions, results are basically consistent before and after this PR.

Setup Before this PR After this PR
Wenchao's setup image image
design not touching boundaries image image
add PML

@mawc2019
Copy link
Contributor

mawc2019 commented Mar 29, 2022

I reran the test with Python 3.8 in SuperCloud, and the results are exactly the same as my previous results using Python 3.7. You may also rerun the test in Supercloud when you have a chance.

In your previous results, one of the finite-difference gradients become exactly zero after this PR. Maybe that is not something usual:

Finite-difference gradient: [ -97.65625 -361.328125 -87.890625 0. -351.5625 ]

@mochen4
Copy link
Collaborator Author

mochen4 commented Mar 30, 2022

I reran the test with Python 3.8 in SuperCloud, and the results are exactly the same as my previous results using Python 3.7. You may also rerun the test in Supercloud when you have a chance.

I just did, and the results are actually similar to what I have locally. I tried other cases as well, and the conclusion is consistent with my previous comment: except in your original test case, the gradients don't change much before and after this commit. At least from my perspective, there isn't any issue with this PR. You should try other cases to see whether the adjoint gradients are always worse for you.

Here are the plots for your reference:
image
image
image

@mawc2019
Copy link
Contributor

mawc2019 commented Mar 30, 2022

I just reran the test with the same setup, but this time, I output the gradients with respect to all design pixels and calculated the relative error as ||adjoint_gradient-finite-difference_gradient||/||finite-difference_gradient||, where || || denotes the 2-norm.
As shown in the figures below, when the design region does not touch the cell boundary, the results before and after this PR are almost the same. However, when the design region touches the cell boundary, the results are nearly perfect before this PR but become significantly worse after this PR. Moreover, in the worst case shown in the bottom left panel, the pattern displays a mirror symmetry with respect to the line of x=y.
image

I tried another setup, and the results are similar.
image

Besides, the adjoint gradients in the Near2Far and FourierFields solvers still vary with the number of processes even with this PR. For example, we can slightly change the previous setup from FourierFields to Near2Far, with

far_pt = [mp.Vector3(0,20,0)]
NearRegions = [mp.Near2FarRegion(center=mp.Vector3(0,design_height+0.2), size=mp.Vector3(Lx,0), weight=+1)]
FarFields = mpa.Near2FarFields(sim, NearRegions ,far_pt)
ob_list = [FarFields]

def J(alpha):
    return npa.sum(npa.abs(alpha[0,:,2])**2)

The relative difference between 1 process and 8 processes is calculated as 8-process adjoint gradient / 1-process adjoint gradient − 1, just like in #1555. The result is as follows, which has a similar pattern as that in #1555:
image

In addition, sometimes the adjoint gradient blows up. Here is an example where the design region does not touch the cell boundary.

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 = 3,3,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 = 0.8,0.008
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)),1

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.35),size=mp.Vector3(Lx,0)),mp.Ez,yee_grid=False)
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])

When running with 16 processes, the code yields the adjoint gradients as
Adjoint gradient: [1.63547167e+153 3.12609241e+153 5.74916256e+153 2.16949324e+153 1.72648562e+153]
with the warning message
RuntimeWarning: invalid value encountered in double_scalars return (change/closure['maxchange']) <= tol and _sim.round_time() >= minimum_run_time
The adjoint gradient does not blow up and the warning message does not appear before this PR with the above example.

@mochen4
Copy link
Collaborator Author

mochen4 commented Mar 31, 2022

Besides, the adjoint gradients in the Near2Far and FourierFields solvers still vary with the number of processes even with this PR.

I tested this as well previously with your code. The adjoint gradients actually agreed with different numbers of processes after this PR but not before. In your plot, notice that the relative differences are now at 1e-5, while in #1555 they went up to almost 1e+2.

In addition, sometimes the adjoint gradient blows up. Here is an example where the design region does not touch the cell boundary

I tried the code on Supercloud and actually didn't see any blowup in gradients or any warning.

@mawc2019
Copy link
Contributor

notice that the relative differences are now at 1e-5

This difference is indeed small but much larger than the machine precision, so something unusual still seems to be in the code. I also have seen a FourierFields case where the variation is large but I have not posted it here.

I tried the code on Supercloud and actually didn't see any blowup in gradients or any warning.

Did you use 16 processes? This problem probably does not appear with a small number of processes.

Maybe we need to confirm that our results are the same at first.

@mochen4
Copy link
Collaborator Author

mochen4 commented Mar 31, 2022

This difference is indeed small but much larger than the machine precision, so something unusual still seems to be in the code. I also have seen a FourierFields case where the variation is large but I have not posted it here.

I would argue that the differences are generally much smaller than before. There might be other places that contribute to the differences, but not necessarily an issue with this PR.

Did you use 16 processes? This problem probably does not appear with a small number of processes.

I did, and it worked fine.

@smartalecH
Copy link
Collaborator

This difference is indeed small but much larger than the machine precision

We don't expect the difference to reach machine precision. A relative error of 1e-5 is actually pretty good.

@mawc2019
Copy link
Contributor

mawc2019 commented Mar 31, 2022

There might be other places that contribute to the differences, but not necessarily an issue with this PR.

That is possible. Probably this PR itself is perfect, but still does not suffice to completely resolve the problem.

We don't expect the difference to reach machine precision. A relative error of 1e-5 is actually pretty good.

Some other examples may give much more obvious differences. For example, in the previous case, the adjoint gradient is normal at 1 process but diverges at 16 processes:

When running with 16 processes, the code yields the adjoint gradients as
Adjoint gradient: [1.63547167e+153 3.12609241e+153 5.74916256e+153 2.16949324e+153 1.72648562e+153]

However, Mo did not see the same thing.

@mawc2019
Copy link
Contributor

mawc2019 commented Apr 1, 2022

It seems that Meep still has the same issue as #1555, although it may become moderate. Here is another Near2Far example, which is slightly changed from a previous FourierFields example. In this Near2Far example, for some design pixels, the relative difference in the adjoint gradient can be more than 15% between 1 and 2 processes. The relative difference here is defined as 2-process adjoint gradient / 1-process adjoint gradient − 1.
image

Among 2, 4, 8, and 16 processes, the adjoint gradients show negligible differences. But at 32 processes, the differences become obvious again. Here is the relative difference between 16 and 32 processes, i.e., 32-process adjoint gradient / 16-process adjoint gradient − 1.
image

For finite-difference gradients, the relative differences among different numbers of processes are negligible, which are around 1e-10.

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 = 3,3,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 = 0.8,0.008
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)),1

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)

far_pt = [mp.Vector3(0,20,0)]
NearRegions = [mp.Near2FarRegion(center=mp.Vector3(0,design_height+0.35), size=mp.Vector3(Lx,0), weight=+1)]
FarFields = mpa.Near2FarFields(sim, NearRegions ,far_pt)
ob_list = [FarFields]

def J(alpha):
    return npa.sum(npa.abs(alpha[0,:,2])**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, Nx*Ny
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()

numb_proc = 1 # number of processes
np.savetxt("finite_far_"+str(numb_proc)+"proc_0401.out",g_discrete)
np.savetxt("adjoint_far_"+str(numb_proc)+"proc_0401.out",g_adjoint)
print("Randomly selected indices: ",idx)
print("Finite-difference gradient: ",g_discrete)
print("Adjoint gradient: ",g_adjoint[idx])

@mawc2019
Copy link
Contributor

mawc2019 commented Nov 23, 2022

It seems that this PR also works for eigenmode coefficients, but the results between different numbers of processes can still be noticeably different. Some 3d tests for EigenmodeCoefficient were performed on the latest version of Meep, namely, v1.25.0. The comparisons of results are shown here. The relative differences in the results between different numbers of processes can be a few percent to 80 percent. As such an issue was already fixed by this PR, some bug might somehow be reintroduced after this PR.

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

Successfully merging this pull request may close these issues.

Adjoint gradients in Near2Far change with the number of processes
8 participants