Bug in SPEC (introduced by Stuart) about spec HOT 14 CLOSED

I have been looking at the jo00aa.h routine very closely. I do understand every part of it now – i.e., what the routine is doing and why.

Comparing the version of jo00aa.h from master and that from SPEC1.0, I could see no bug introduced in the master version, and so I conclude that the problem may arise from either of these differences: (i) some inut paramaters and variables (halfmm vs regumm, lrad vs Lrad, cheby vs chebyshev); (ii) the way the vector potential is stored in the master version; (iii) the routine for returning coordinates, metrics, etc. is coords.h vs co01aa.h and is called differently.

Also, it seems to me that the problem (beltrami error insensitive to resolution) may only exist in cartesian geometry...

from spec.

I have made some progress in solving this issue. I realized that the code is (in slab geometry) not solving the Beltrami equation, curl(B) = muB, but its negative! Namely, it is solving cur(B) = - muB.

In fact, if I modify the file jo00aa.h, in which the error E=curl(b)-mu*B is calculated, replacing

jerror(ii) = jerror(ii) + weight(jquad) * sum( abs( gJu(1:Ntz,ii) - mu(lvol) * gBu(1:Ntz,ii,0) ) )

with

jerror(ii) = jerror(ii) + weight(jquad) * sum( abs( gJu(1:Ntz,ii) + mu(lvol) * gBu(1:Ntz,ii,0) ) )

then I get very good convergence of the error with radial resolution Lrad!

This also implies that this issue is directly related to the other issue I recently opened, called "sign error in the poloidal flux".

from spec.

If I modify the file mp00ac.h, namely the Beltrami linear solver routine, by replacing

matrix(1:NN,1:NN) = dMA(1:NN,1:NN) - lmu * dMD(1:NN,1:NN)

with

matrix(1:NN,1:NN) = dMA(1:NN,1:NN) + lmu * dMD(1:NN,1:NN)

then I get very good convergence of the error with radial resolution Lrad (and the poloidal flux also has the expected sign).

from spec.

If I modify the file mp00ac.h, namely the Beltrami linear solver routine, by replacing

matrix(1:NN,1:NN) = dMA(1:NN,1:NN) - lmu * dMD(1:NN,1:NN)

with

matrix(1:NN,1:NN) = - dMA(1:NN,1:NN) - lmu * dMD(1:NN,1:NN)

the results also appear correct and converge!

from spec.

The problem is present in any geometry. It was not visible in, e.g., the l2stell case (G3V02L1Fi.001.sp) because the Beltrami coefficient was basically zero. For equilibria in which mu is non-zero, the error from jo00aa.h is not converging, unless one makes either of the two modifications mentioned above.

from spec.

This could explain why I was having trouble with my ITER equilibria (mu non-zero). This would be consistent as a very old version of the code produced meaningful equilibria. Is there a branch I can test?

from spec.

You could create a branch called, e.g. "beltrami_fix", from the master branch, then make the modification that I suggested, e.g. modify the file mp00ac.h, namely the Beltrami linear solver routine, by replacing

matrix(1:NN,1:NN) = dMA(1:NN,1:NN) - lmu * dMD(1:NN,1:NN)

with

matrix(1:NN,1:NN) = dMA(1:NN,1:NN) + lmu * dMD(1:NN,1:NN)

and see if the run you had problems with now runs well.

I am trying to see where the origin of the problem lies. Not sure yet. But clearly the beltrami equation is not solved. Its negative is.

from spec.

@jloizu Quick update. The axisymmetric version computes an equilibria, with the sign change. I'm trying the n=4 case now will take a few hours at least.

from spec.

@SRHudson : I have been trying to compare the dMA and dMD matrices that are constructed to solve the Beltrami equation, and I can see that they are different from the matrices computed in the old (SPEC1.0) version of the code - for the same equilibrium calculation.

More striking is the fact that the sizes are quite different. These matrices have size NN x NN, where NN=NAdof is the number of degrees of freedom for the field. In the GIT version of SPEC, in slab geometry and stellarator symmetry,

NAdof(vvol) = 2 * ( mn ) * ( Lrad(vvol)+1 ) + mn + mn + mn-1 + 1 + 1

while in the SPEC1.0 version, we had

Nmagneticdof(vvol) = ( mn ) * 2 * ( Lrad(vvol)-1 ) + ( mn-1 ) = 2

For example, for Nvol=1, Mpol=Ntor=0, Lrad=2, then mn=1 and we have

NAdof=10
Nmagneticdof=2

so the matrices are very different...

from spec.

Hi, At some time I changed the construction of the Beltrami matrices. Originally I had enforced the boundary conditions on the interfaces that the normal magnetic field be zero by constraining the representation of the vector potential. This made the construction of the matrices very complicated, as some harmonics of the vector potential depended on others. Later, I choose to treat all the harmonics of the vector potential as independent degrees of freedom and to enforce all the constraints (even the contraints on the enclosed toroidal and poloidal fluxes) using Lagrange multipliers. This does mean that there are more independent degrees of freedom, but the source is so much easier to read and debug that I think it is worth the very small increase in computational cost.
Does this make sense?
P.S. I am finishing up my topic with FOCUS. I want to finish this first because I will use FOCUS to get a good set of coils for free-boundary SPEC, and I am updating OCULUS, which has a routine for constructing the input for SPEC given the coils. 2018 will be a HUGE year for free-boundary SPEC.

from spec.

OK, I understand that the matrices have different sizes then. The strict comparison is thus meaningless.
However, I believe that the current (GIT) version of SPEC has a problem. It is NOT solving cur(B)=muB.
I don't understand where the problem originates, but clearly the code is solving cur(B)= - muB ! @SRHudson : Do you have any clue?

PS: yes, I foresee a good year 2018! If you say that OCULUS has a routine for constructing the input for SPEC given the coils, then I should not bother learning how to use a Biot-Savart solver given coils, right? I had started this with Carolin...but maybe I should just wait for OCULUS?

from spec.

The bug must be coming from either of these 3 possibilities:

1. The matrix dMA has the wrong sign (in mp00ac),
2. The matrix dMD has the wrong sign (in mp00ac),
3. The matrix to be inverted should be M = dMA + mu * dMD and not M = dMA - mu * dMD
   Looking at how dMA and dMD are calculated, I cannot see any sign error so far...

@SRHudson : could it be then (3)?

This is a very important problem to solve, and free-boundary SPEC cannot be claimed to be verified until this is resolved.

from spec.

After the last commits of @SRHudson , this problem has been fixed! It was option (2) then, namely there was an error in the sign of the dMD matrix. Now, the code really solves curl(B)=mu*B!

I tested compilation of release and debug versions, in both ipp clusters, and tested execution of the code, which successfully runs in both test cases (G3V02L1Fi.001 and G1V02L0Fi.001). I also checked that the resulting interface geometry is the same as before (essentially to machine precision). And now, the right sign of the Beltrami solution is obtained. I also checked that now (unlike before) the error (computed in jo00aa) converges well to zero with radial resolution.

Problem solved!

from spec.

Thanks Stuart.
Closing issue.

from spec.