NanoComp / meep

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

Adjoint optimisation involving LDOS objective #2799

Closed zekishaw closed 7 months ago

zekishaw commented 7 months ago

Hi, Im having some trouble with the abovee.

Particularly, in any scenario where I have LDOS included as an objective I get a value error in nlopt:

Full traceback is:

File ~/.miniforge3/envs/meep1p27_spy/lib/python3.11/site-packages/spyder_kernels/py3compat.py:356 in compat_exec exec(code, globals, locals)

File ~/Code/MEEP/adjoint/adjoint_int.py:1050 opt = run_opt(x0, maxeval, len_constraint=True)

File ~/Code/MEEP/adjoint/adjoint_int.py:334 in run_opt x[:] = des_opt.optimize(x)

File ~/.miniforge3/envs/meep1p27_spy/lib/python3.11/site-packages/nlopt.py:328 in optimize return _nlopt.opt_optimize(self, *args)

ValueError: nlopt invalid argument ..

I have done some debugging and checking, "in lcf, x is: ((47089,), None), eta is 0.5, beta is 10", and when trying to output the filtered field as desired, I get very low numbers. (x*10^-18).

I can provide a minimum worked example, both before and after trying a length constraint, but I wandered if there may be something immediately obvious...

I also just want to make sure I understand what is going a little, so to check:

the field filtering is taking the design weights after an optimisation run, and is then mapping the design weights based on any length constraint you have put in? (so its not neccessarily the dft field?)

i.e. it makes sure that x weights are appropriately distributed such that the line constraint is obeyed.

The projected field takes this information, and according to how high you are thresholding (beta) applies these weights as weights within the optimisation geometry.

des_opt.set_max_objective(lambda a,g: lcf(a,g, beta)) -- this lambda function here takes in a, g, which is the output of an nlopt run, where a is the x-array, g is the gradients output?

f0, dJ_du = opt([mapping(v, eta_i, beta)]) -- and then this line here simply takes the value of the objective function, and the gradient, from the mapping of v [which is the x value taken into the cost function from nlopt], eta_i (interpolation point where you could be either material) and the beta value?

Let me know if a minimum worked example will be of use, I have been using the metalens far fields as a simple example to play with, so its essentially that example, optimising for LDOS instead. [As there is no example for LDOS at the minute and help here will be of such massive use to me I would also happily write a jupyter notebook as a tutorial, in return]

zekishaw commented 7 months ago

@mochen4 -- I've just tagged yourself as you seem to have masterminded the LDOS adjoint objective.

mochen4 commented 7 months ago

ValueError: nlopt invalid argument ..

This error message usually has nothing to do with meep adjoint. Instead, you probably didn't set up the NLopt optimization correctly; for example, the size and dimension of arrays don't match, or not in the right format etc. You can just add print statement to relevant places, and see which has problems. All meep function calls should work properly, so for debugging purpose you can just have a small maximum_run_time since accuracy isn't an issue for debugging. Debugging may be a little tedious, but there shouldn't be anything technical.

I have done some debugging and checking, "in lcf, x is: ((47089,), None), eta is 0.5, beta is 10", and when trying to output the filtered field as desired, I get very low numbers. (x*10^-18).

I don't know what you were printing here, but if "None" was for gradient, it looked like the it was not assigned. I also don't know what you mean by filtering here.

the field filtering is taking the design weights after an optimisation run, and is then mapping the design weights based on any length constraint you have put in? (so its not neccessarily the dft field?)

i.e. it makes sure that x weights are appropriately distributed such that the line constraint is obeyed.

The projected field takes this information, and according to how high you are thresholding (beta) applies these weights as weights within the optimisation geometry.

The filtering step takes the raw design weights and "smeared them out" to some new weights, which are then binarized by projection. The filtering doesn't really make sure the length constraint is obeyed; it only softly spread out the design weights according to the filter radius; beta specifies the strength of projection.

Filtering has nothing to do with dft fields, which are, roughly speaking, Fourier transformed electromagnetic fields stored in Meep.

des_opt.set_max_objective(lambda a,g: lcf(a,g, beta)) -- this lambda function here takes in a, g, which is the output of an nlopt run, where a is the x-array, g is the gradients output?

a,g are input for functions used by NLopt. NLopt generally takes functions of format that has variable and gradient array as input, and the function will modify gradient in place, and return some function value.

f0, dJ_du = opt([mapping(v, eta_i, beta)]) -- and then this line here simply takes the value of the objective function, and the gradient, from the mapping of v [which is the x value taken into the cost function from nlopt], eta_i (interpolation point where you could be either material) and the beta value?

mapping takes the raw design weights and process them through filtering and projection; the returned processed weights are used in the actual simulations; and opt(...) computes the function value and gradients for that processed weights.

Let me know if a minimum worked example will be of use, I have been using the metalens far fields as a simple example to play with, so its essentially that example, optimising for LDOS instead. [As there is no example for LDOS at the minute and help here will be of such massive use to me I would also happily write a jupyter notebook as a tutorial, in return]

As I said, I believe the issue is more related to your usage of NLopt. Checking the NLopt documentation and make sure you understand how it works might be more useful.

zekishaw commented 7 months ago

It appears it was an nlopt issue all this time!

I changed the output from

"np.real(f0)" to "np real(f0)[0]".

Thank you for you answers and for the advice