NanoComp / meep

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

Minimal Optimization at Higher Resolutions #1825

Closed WillE9a closed 2 years ago

WillE9a commented 2 years ago

Hi all,

I've run into an issue where the dielectric distribution does not change much at all when the resolution is increased. I've included a small 2D example code to maximize coupling to a single wave guide that demonstrates the problem. I am not sure if it's just an issue with the optimization finding a worse local maxima when using the higher resolution/larger parameter space or something else, and I was hoping someone could explain what exactly is causing this behavior. Here is the code, results summarized below: X_dipole_single_wg.zip initial_structure

Starting with a uniform initial permittivity distribution of 1/2, the optimization is run at 20 resolution (design resolution is the same as overall resolution) and the resulting structure has significant changes to that distribution. Then, when I run the same optimization at 30 resolution, there is very minimal change to the structure and the final permittivity is very close to the initial distribution.

Here is the 20 resolution run: eps_20res FOM_20res

Here is the 30 resolution run: eps_30res FOM_30res

There are a few additional things I've checked/would like to check:

I appreciate any feedback on if this behavior is expected or abnormal, and if there is any advice or things to watch out for when rescaling the 2D permittivity grid.

Thanks in advance!

smartalecH commented 2 years ago

From your posted code, it looks like you aren't using any filters.

Filters are an important step within density methods, as they help regularize (and even condition) the optimization problem. There are quite a few papers that discuss this (maybe look at Sigmund's review paper).

The tutorials also discuss this in depth (although they haven't been updated in awhile and may have bit-rotted...)

oskooi commented 2 years ago

There is also another aspect of your design that may impact the performance of the optimizer. The Cylinder object at the center of the design region is effectively "blocking" the design voxels underneath. As a test, try removing this Cylinder and see whether it affects the results as you vary the resolution.

geometry = [
    mp.Block(center=mp.Vector3(y=Sy/4), material=GaAs, size=mp.Vector3(waveguide_width, Sy/2, 0)),  # vertical waveguide
    mp.Block(center=design_region.center, size=design_region.size, material=design_variables), # design region
    mp.Block(center=design_region.center, size=design_region.size, material=design_variables,
             e1=mp.Vector3(x=-1), e2=mp.Vector3(y=1)), # enforce mirror symmetry about x=0 plane
    mp.Cylinder(center=design_region.center, material=GaAs, radius=buffer_radius) # buffer region to protect QD
]
WillE9a commented 2 years ago

Thanks for the feedback @smartalecH @oskooi.

I have been trying to implement the filtering process that the examples use but keep running into an issue where the materialGrid dimensions are off by one (i.e., for the 'filtering and thresholding design parameters' example it tries to reshape the fixed 1D array of size 5625 into a 2D array with shape (76, 76)).

Inside the conic_filter function, I found that it is adding 1 to my Nx & Ny dimensions, but in a previous version, there was no change to the materialGrid dimension. Should the filtering functions have this change to material grid dimension? I am using pymeep 1.21 currently and unable to run that example.

smartalecH commented 2 years ago

I have been trying to implement the filtering process that the examples use but keep running into an issue where the materialGrid dimensions are off by one

Yes, the tutorials have bit-rotted, unfortunately. They need to be updated (the adjoint package changes rather frequently, which is why there isn't really any documentation yet).

Try simply adding 1 to your Nx and Ny variables (such that the corresponding design grid expands by 1 in each dimension).

Something like

design_region_resolution = resolution*2
design_region_width = 2
design_region_height = 2
Nx = int(design_region_resolution*design_region_width) + 1
Ny = int(design_region_resolution*design_region_height) + 1

You can always look at the adjoint test for reference (since that should never bit-rot...)