Reimplement helicity constraint

[zhisong]

I did a naive attempt to enable helicity constraint. Without modifying the already contained subroutines, I just reversed the changes that disabled the helicity constraint. This is in branch "helicity-fix". One needs to set the option LBeltrami= 2.

I have also uploaded some random example to calculate the Taylor state in a large aspect ratio torus, in Inputfiles/Verification/helicity . There are two input files with the same helicity but different initial guess of mu and enforced symmetry, thus different configurations.

Frankly speaking, now I can have very little control on choosing the global minimum energy state or converging to a local one. Even setting the initial guess of mu to exactly the result after the computation does not guarantee the convergence to the same mu. I think we need to work further for a rigorous benchmark.

[jloizu]

Thanks Zhisong for opening up this issue.

I don't understand how Lconstraint=0 (which presumably fixes mu and fluxes) and LBeltrami=2 (which presumably treats mu as an independent degree of freedom) are compatible and if they really allow for calculating the field at fixed helicity and not fixed mu.

[zhisong]

It seems that Lconstraint only works for the linear solver (LBeltrami=4). For LBeltrami=2, the option Lconstraint is not used.

[zhisong]

I have benchmarked the helicity constraint with the Taylor analytic solution in the cylindrical geometry. It seems to match quite well.

The Mathematica file that generates this plot is attached. One will need to change the extension back to .m from .txt.
I have updated the helicity-fix branch. The input files that generate this result can be found in InputFiles/Verification/helicity/Taylor

One trick I did was to randomize the initial guess. This helps to converge better to the helical state. If we give zero to the initial guess of the field, the solution will always be the axisymmetric one. Please see the modification in preset.h.

taylor.nb.txt
taylor.pdf

[jloizu]

I tested the helicity-fix branch (which compiles and executes) in slab geometry to see what is the Hessian obtained from SPEC in a case where I have the analytical solution for its elements. It seems that when LBeltrami=2, the helicity is indeed constrained when calculating a Beltrami field, but the Hessian looks exactly like the one computed at fixed mu (and not at fixed helicity).

Could it be that during the calculation of the force derivative (i.e. the hessian) with respect to geometry SPEC still uses the linear solver for the Beltrami fields and thus considers mu as given?

[jloizu]

The reason I am saying this is because in dforce, when LComputederivatives=1, it looks like the matrices dMA, dMB, etc. are used to invert the Beltrami solution, but if I understand correctly that should only be done if the Beltrami field is calculated using the linear method, right?

[zhisong]

My understand is, dMA, dMB are used for any methods. The different methods are different algorithms to minimize the energy functional in the quadratic form (or cubic for \mu * helicity).
To properly implement helicity constraint in Hessian, one will need to look at how the parameter mu changes if the interfaces are moved, and add the contribution to the Hessian. This is non-trivial.

[zhisong]

And my worry is that, the matrix dMA - \mu dMB in bifurcation cases (e.g. helical Taylor states) will be singular and non-invertible, since \mu is the eigenvalue. Therefore, we will not able to use the same way as the linear method to construct the Hessian.

[zhisong]

I have attached a derivation of how the Hessian should be found. Please see if you agree or not. @jloizu @SRHudson
Hessian.pdf

[jloizu]

Here is a funny discovery:

Running SPEC with the same input file twice gives different results! This is for the cases with LBeltrami=2 where the helicity is constrained. For example, the test case uplodaded by @zhisong,

SPEC/InputFiles/Verification/helicity/Taylor/

I understand that the Newton method is sensitive to the initial guess and small differences in the guess can produce different results. But here it is by running the exact same input file that the results are different! One has to run several times to actually see a different result appearing. So it is relatively robust, but not quite! This is very disturbing...please try it and let me know! (you may set Wma02aa=T in order to better see the obtained mu)

[zhucaoxiang]

@jloizu I have seen this in other codes before. I guess this might be caused by the ill-conditioning of the Hessian matrix. To examine it, we could do a SVD or eigenvalue decomposition on the Hessian matrix and compare the largest and smallest singular value / eigenvalue.

[jloizu]

@zhucaoxiang what do you mean here by the Hessian? Because this case is for Nvol=1, namely no force-balance required. Do you mean the Jacobian used by the Newton method for the calculation of the Beltrami field?

[zhucaoxiang]

@jloizu Yeah, I guess the Newton method will invert the Jacobian to find descent direction. Correct me if I was wrong. The machine round-off error will be amplified when inverting an ill-conditioning matrix.

[zhisong]

What @zhucaoxiang suggested could be one of the problem.
The source of randomness may also come from the fact that I randomized all the field variables(initial guess) to some small number. If I leave them zero, the Newton iteration will seldom converge to nonaxisymmetric solutions.
I agree the randomness is not elegant. We should find another more reliable way to initialize.

[jloizu]

@zhisong in which commit did you randomize the field initial guess? Does this apply when Linitgues=1? I cannot find this. Of course that would already explain why a result may not be reproducible. It does not sound like a good strategy, right?

@zhucaoxiang I am surprised that in the very simple slab axisymmetric case I sent the matrix to invert is so ill-posed. I would naively think that there is a problem in the approach itself, because this should be a very simple system to solve and there are no bifurcations.

[zhisong]

@jloizu That was in the latest commit of branch helicity-fix.
I have add the following lines in preset.h

call random_seed()

call random_number(Ate(vvol,ideriv,ii)%s)
call random_number(Aze(vvol,ideriv,ii)%s)
if (.not. YESstellsym) then
    call random_number(Ato(vvol,ideriv,ii)%s)
    call random_number(Azo(vvol,ideriv,ii)%s)
endif

[zhisong]

Following what @SRHudson suggested, I have made the random initialization an option rather than compulsory. Now if we set Linitgues=3, the initial field will be randomized. The maximum of the random number is controlled by "maxrndgues" in the namelist "locallist".

I have pushed this commit to the branch helicity-fix.

[jloizu]

I tried the new version of the branch helicity-fix with Linitgues=1 and now the results are reproducible (as expected from the lack of random guess) and also the code successfully finds Beltrami solutions with constrained helicity, at least for simple slab cases I was having trouble with before. Thanks @zhisong !

[arunav2123]

Helicity in hessian matrix has been implemented for LBeltrami = 2 i.e Newon's Method. A separate branch named "SPEC-helicity-hessian" has opened for review. Results looks promising.....

[arunav2123]

First plot , using "SPEC-helicity-fix" with LBeltrami =2, Lconstraint=0

Second plot, using SPEC-helicity-hessian with LBeltrami=2, Lconstraint=0

For both case: Input file "G1V05L0Fi.001.sp.h5' has been used....

[arunav2123]

Three Input files has been uploaded i.e .
For Igeometry = 2 : Cylinder_vol2.sp, cylinder_vol2_compare.sp
For Igeometry =3 : G3V02L1Fi.001.sp

[jloizu]

I have pulled the branch helicity-hessian and verified to machine precision that the Hessian at fixed helicity is properly calculated. I have used a 3-volume slab case (one with a stable current sheet, one with an unstable current sheet) that I have studied analytically in detail and for which I have exact predictions for the Hessian at fixed helicity, see paper below

https://aip.scitation.org/doi/10.1063/1.5091765

Well done @zhisong and @arunav2123 ! I will upload the verification cases with some description.

[jloizu]

You can find now the verification of the Hessian to machine precision here:

SPEC/InputFiles/Verification/hessian/fixed_helicity/

Notice this is as of now only in the helicity-hessian branch.

[zhisong]

We could merge this branch. However, the hessian for toroidal geometry needs further development.

[jloizu]

@zhisong Yes, I agree. I would suggest that you add a warning if LBeltrami=2 and Igeometry=3 and Nvol>1, saying something like "Hessian needs to be revised/verified". Then, we could merge.

[zhisong]

Done. I have also screened warning messages for the Newton's method. I suggest we continue this thread for future work of Igeometry=3.