Spec dofs

[abaillod]

I finally found the time to document and improve my code. In this branch I implemented some additional features for optimizing SPEC equilibria:

• Can now use any input profile as dofs, including the toroidal flux profile (i.e. the position of SPEC interfaces). See for example the test in tests/mhd/spec_test/SpecTests.test_optimize_net_toroidal_current()
• Can now use the Vns (i.e. B.n on the computational boundary) as dofs, to perform quasi-free-boundary optimization. These harmonics are stored in a NormalField object, implemented in src/simsopt/field/normal_field.py
• Spec now stores the latest converged equilibrium, and uses it as initial guess for the next iteration. Interfaces are stored as SurfaceRZFourier objects.

These new features have been used to obtain the results presented in this publication. I am still working on the coupling of SPEC with booz_xform and the Redl formulas for the computation of the bootstrap current. These changes will come in a future PR.

Please let me know if you have any question or if anything needs to be improved!

[abaillod]

@landreman, I think I implemented all your suggestions. For some reasons, there are however some tests that are not successful. I will have a look at this tomorrow

[abaillod]

@landreman This is ready to be reviewed. I don't know why the docker container / docker check is failing. It seems a password is missing?

[mbkumar]

Don't worry about the docker container test failing. Bharat Medasani Engineer Princeton Plasma Physics Lab (PPPL)

…

On Thu, May 4, 2023 at 10:05 AM abaillod ***@***.***> wrote: @landreman <https://github.com/landreman> This is ready to be reviewed. I don't know why the docker container / docker check is failing. It seems a password is missing? — Reply to this email directly, view it on GitHub <#304 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AA62VEHKAA4G6BJKU3XVJETXEOZT5ANCNFSM6AAAAAAXFE2AIA> . You are receiving this because your review was requested.Message ID: ***@***.***>

[landreman]

Looks good, just left a couple questions. The remaining issue is why the spec integrated tests are failing in the CI.

[landreman]

Are these new lines necessary?

[abaillod]

In the old version of spec.py, the coordinate axis given as input to SPEC, i.e. rac and zas, were ignored. Instead, these arrays were filled with zeros and SPEC constructed its own guess for the coordinate axis.

I modified spec.py to read the initial guess in SPEC input file. See these lines and these lines.

When using the initial guess for the axis in QH-residues.sp, the residue1 differs by about 0.5% from the value in the test. Another fix would be to modify tests/test_files/QH-residue.sp and remove the axis from it. A similar fix had to be made to examples/2_Intermediate/inputs/QH-residues.sp so that the example examples/2_Intermediate/eliminate_magnetic_islands.py still works. Let me know what you think is best!

[mbkumar]

Why is the residue1 value changing by 0.5% based on the initial guess? Is it settling at a different local minima? Is that acceptable?

[landreman]

When using the initial guess for the axis in QH-residues.sp, the residue1 differs by about 0.5% from the value in the test. Another fix would be to modify tests/test_files/QH-residue.sp and remove the axis from it. A similar fix had to be made to examples/2_Intermediate/inputs/QH-residues.sp so that the example examples/2_Intermediate/eliminate_magnetic_islands.py still works. Let me know what you think is best!

Got it. I'm indifferent then between 3 options: (1) including those lines to set the axis to 0, (2) editing QH-residue.sp so the axis is 0 there, or (3) loosening the tolerance of the test so it passes without setting the axis to 0. For option (1), it might be useful to state in a comment why those lines are there.

[abaillod]

@mbkumar since a finite resolution is used when constructing the SPEC equilibrium, changing the coordinate axis position (i.e. and also the coordinate grid) changes the equilibrium, and thus also the value of the residue.

I will change the input file QH-residue.sp to have no initial guess for the coordinate axis, for consistency with the file in ./examples/

[landreman]

To what extent does the default update_guess=True confuse finite differencing and make optimization depend on the number of MPI processes? If the initial guesses for the axis and internal surfaces are not the same for every function evaluation, then the objective function depends slightly on the history of previous objective function evaluations on that process. That history in turn depends on how many MPI processes are used. For finite difference derivatives, where tiny differences in the objective function are particularly important, this could be an issue.

[abaillod]

Good point - I haven't thought about this.

My first approach would be to set update_guess=False whenever a worker calls spec.run() to evaluate a finite-difference derivative, and only call with update_guess=True when a step is taken by the optimizer. I guess this means we would have to modify least_square_mpi_solve._f_proc0 and least_square_mpi_solve._constrained_mpi_workers_task?

I don't think this is ideal - I would prefer to only modify spec.py. Any idea?

[abaillod]

I thought about this and concluded that it probably has no impact on the evaluation of the finite differences.

When evaluating finite differences, only small deviations dx of each degree of freedom x are considered around x_0. The equilibrium obtained by SPEC at x_0+dx is thus very close to the one at x_0, and can also be used as an initial guess. The way I understand it is that the objective function depends on SPEC equilibrium, but does not depend on SPEC initial guess. What matters is that the initial guess is good enough for SPEC to find the equilibrium. If this is the case when the initial guess is the solution at x_0, it is also certainly the case if the initial guess is the solution at x_0+dx.

The only possible issue is if there are bifurcations, and SPEC solution (and therefore also the objective function) depends on whether the initial guess is the solution at x_0 or x_0+dx. In this case, however, the problem is probably ill-posed (in the sense that SPEC force-minimization landscape is not smooth).

Does that make sense?

[landreman]

When evaluating finite differences, only small deviations dx of each degree of freedom x are considered around x_0. The equilibrium obtained by SPEC at x_0+dx is thus very close to the one at x_0, and can also be used as an initial guess. The way I understand it is that the objective function depends on SPEC equilibrium, but does not depend on SPEC initial guess.

Ideally, if the equations are solved in spec to very tight tolerance, the objective should not depend on the initial guess. But for finite tolerances, the initial guess may affect the result a bit, and it only takes a small change to the objective to cause an O(1) change to a derivative computed by finite differencing. As an indication that we might want to worry about this, consider the issue discussed in this PR of the residue changing by 0.5% when the coordinate axis shape changed. Also note that when there are multiple worker groups, when an MPI process in group > 0 starts evaluating the objective at x_j + dx for a Jacobian, the initial condition will come from a run from the previous Jacobian evaluation, at x_{j-1}, which may not be super close to x_j.

It's fine with me to have this default update_guess = True, since it's working for you. Just wanted to flag this in case it affects your calculations.

[abaillod]

You are right, we do not necessarily solve SPEC equations to machine precision, which might cause issues when evaluating finite differences. As I see it, we have two choices:

(1) Leave it as it is, with the risk that the different initial guesses might affect the Jacobian evaluation.

(2) Enforce the initial guess to be the same for each finite-difference evaluation. The only (simple) way to do it is to modify the calls to the objective function from the solver and give an additional argument that contains the information on whether or not this call is for a finite-difference evaluation.

For example, in finite_difference.py,

options = dict('finite_diff'=True)
out = np.asarray(self.fn( options ))

and then, in spec.py,

def run(self, options):
    if options['finite_diff']:
        update_guess=False
    else:
        update_guess=True

The downside is that all objective functions must accept an options argument... which might require some changes in different files.

While option (2) is probably better, I think leaving it for another PR is best. What do you think?

[landreman]

I agree that this can be addressed later in a separate pull request.

[mbkumar]

Hello Antoine,

Thanks for the PR with the code improvements and additional features. My comments are minor, but please address them.

[mbkumar]

Why is the residue1 value changing by 0.5% based on the initial guess? Is it settling at a different local minima? Is that acceptable?

[landreman]

I think this is ready to be merged as soon as the failing integrated tests are fixed. The error message is different in https://github.com/hiddenSymmetries/simsopt/actions/runs/4916745286/jobs/8780906806?pr=304 vs https://github.com/hiddenSymmetries/simsopt/actions/runs/4916745286/jobs/8780906678?pr=304 but it appears the failing example is stellarator_benchmarks/1DOF_circularCrossSection_varyAxis_targetIota_spec.py in both cases. @abaillod are you able to fix this?

[abaillod]

I think this is ready to be merged as soon as the failing integrated tests are fixed. The error message is different in https://github.com/hiddenSymmetries/simsopt/actions/runs/4916745286/jobs/8780906806?pr=304 vs https://github.com/hiddenSymmetries/simsopt/actions/runs/4916745286/jobs/8780906678?pr=304 but it appears the failing example is stellarator_benchmarks/1DOF_circularCrossSection_varyAxis_targetIota_spec.py in both cases. @abaillod are you able to fix this?

Ok, this is now fixed. The problem was again related to SPEC coordinate axis; I fixed it by adding these lines and by removing the axis from there.