fix: cosmui, sinmui have wrong normalization factors for lasym=T

[jonathanschilling]

This is an issue to track work related to the case of an (I think) incorrectly initialized dnorm factor in fixaray.f.
The code changes happen on the branch fix_cosmui.

I found this problem when trying to understand the Fourier transform of the forces in VMEC.
The problem appears only when running VMEC with lasym=T.

Any (periodic) function can be decomposed into an even-parity contribution and an odd-parity contribution.
Some force components (armn, brmn, ..., clmn, frcon, fzcon) have even parity, some have odd parity when symmetry is assumed.
In the asymmetric case, the forces are decomposed into definite-parity contributions
and then Fourier-transformed over the poloidal half-interval [0, pi] (1, ..., ntheta2).
The decomposition into definite-parity contributions is done in symforce.
Those-parity contributions that would be the only ones in symmetric mode are Fourier-transformed by tomnsps.
The other-parity contributions are Fourier-transformed using tomnspa.
Concluding, the Fourier integrals computed to get the Fourier coefficients of the forces are only ever computed
over the poloidal half-interval [0, pi] and never over the full-poloidal interval [0, 2pi[ .

The Fourier integrals in tomnsp* make use of the Fourier basis functions
stored in cosmui, sinmui, cosmumi and sinmumi.
These are filled in Sources/Initialization_Cleanup/fixaray.f.
Here, the normalization of the Fourier integrals is being done using the dnorm factor.
Note that this enters also in the flux-surface average weighting factor wint through cosmui3.

dnorm gets changed to 1/(nzeta*ntheta3) if lasym=T.
In that case, ntheta3 equals ntheta1, which refers to the full poloidal interval [0, 2 pi[.
The Fourier integrals done in tomnsp* use this dnorm factor to normalize the discrete integrals there,
which however are still only being computed over half the poloidal range.
I think this is incorrect.

Instead, dnorm as used for cosmui, ..., sinmumi should always be 1/(nzeta*ntheta2-1) as in the symmetric case.
Furthermore, cosmui3 and in turn the weighting factors wint still need to adjust to the changed poloidal interval.
Thus, I propose to introduce dnorm3 in fixaray, which shall get adjusted to reflect the actual number of poloidal grid points.

The actual errors resulting from this mixup seems to be that in the asymmetric case,
the Fourier coefficients of the forces are too low by a factor of ntheta3/(ntheta2-1) which is approx. 2.
VMEC seems to be able to cope with that by just using more iterations,
but I think this should be fixed anyway.

The expected results of fixing this issue are:

• for originally symmetric cases, when switching to lasym=T, the number of iterations should approximately stay constant
• for asymmetric cases, the number of iterations should go down

The SOLOVEV case from the DESC repository is used as a simple test case here:

&INDATA
  ! based on input.SOLOVEV from github.com/PlasmaControl/DESC/examples 
 
  LASYM = F
  
  NFP =    1
  MPOL =   16
  NTOR =   0
 
  NS_ARRAY    = 16    256
  NITER_ARRAY = 1000  20000 
  FTOL_ARRAY  = 1E-16 1E-16 
 
  DELT =   9.00E-01
  NSTEP =  250
 
  GAMMA =   0.000000E+00
  PHIEDGE =   1.00000000000000E+00
  CURTOR =   0.00000000000000E+00
  PMASS_TYPE = "power_series"
  AM =  0.125 -0.125
  NCURR =    0
  AI = 1.0
  
  RAXIS =   3.999
  ZAXIS =   0.00000000000000E+00
  
  RBC( 0,0) =  3.999     ZBS( 0,0) =  0.000
  RBC( 0,1) =  1.026     ZBS( 0,1) =  1.580
  RBC( 0,2) = -0.068     ZBS( 0,2) =  0.010
/
&END

Running this case with the current master branch, I get the following output:

  NS =   16 NO. FOURIER MODES =   16 FTOLV =  1.000E-16 NITER =   1000
  PROCESSOR COUNT - RADIAL:    1

  ITER    FSQR      FSQZ      FSQL    RAX(v=0)    DELT       WMHD

    1  9.56E-02  2.88E-03  3.25E-02  3.999E+00  9.00E-01  2.6381E+00
  250  2.15E-15  6.16E-16  4.15E-16  3.990E+00  9.00E-01  2.5484E+00
  272  8.89E-17  3.66E-17  2.15E-17  3.990E+00  9.00E-01  2.5484E+00

  NS =  256 NO. FOURIER MODES =   16 FTOLV =  1.000E-16 NITER =  20000
  PROCESSOR COUNT - RADIAL:    1

  ITER    FSQR      FSQZ      FSQL    RAX(v=0)    DELT       WMHD

    1  9.94E-04  1.52E-04  1.43E-07  3.990E+00  9.00E-01  2.5484E+00
  250  2.37E-10  5.61E-11  6.49E-13  3.989E+00  5.93E-01  2.5484E+00
  500  4.14E-12  2.00E-13  1.92E-15  3.989E+00  5.93E-01  2.5484E+00
  750  9.16E-13  2.13E-14  3.32E-18  3.989E+00  5.93E-01  2.5484E+00
 1000  1.81E-13  3.04E-15  3.45E-20  3.989E+00  5.93E-01  2.5484E+00
 1250  1.66E-14  2.26E-16  3.81E-22  3.989E+00  5.93E-01  2.5484E+00
 1500  8.72E-16  1.17E-17  1.30E-23  3.989E+00  5.93E-01  2.5484E+00
 1698  1.00E-16  1.41E-18  4.74E-24  3.989E+00  5.93E-01  2.5484E+00

Then lasym=T is put in the input file and the output now reads:

  NS =   16 NO. FOURIER MODES =   16 FTOLV =  1.000E-16 NITER =   1000
  PROCESSOR COUNT - RADIAL:    1

  ITER    FSQR      FSQZ      FSQL    RAX(v=0)   ZAX(v=0)    DELT       WMHD

    1  2.39E-02  7.20E-04  8.13E-03  3.999E+00  0.000E+00  9.00E-01  2.6381E+00
  250  8.19E-12  3.05E-13  4.01E-13  3.990E+00  2.528E-16  9.00E-01  2.5484E+00
  500  4.96E-16  6.97E-18  1.57E-18  3.990E+00 -5.771E-17  9.00E-01  2.5484E+00
  533  9.77E-17  1.49E-18  1.68E-19  3.990E+00 -8.564E-16  9.00E-01  2.5484E+00

  NS =  256 NO. FOURIER MODES =   16 FTOLV =  1.000E-16 NITER =  20000
  PROCESSOR COUNT - RADIAL:    1

  ITER    FSQR      FSQZ      FSQL    RAX(v=0)   ZAX(v=0)    DELT       WMHD

    1  2.59E-04  4.06E-05  4.06E-08  3.990E+00 -8.564E-16  9.00E-01  2.5484E+00
  250  5.93E-12  9.04E-13  1.48E-14  3.989E+00 -2.752E-16  9.00E-01  2.5484E+00
  500  8.62E-13  9.57E-14  3.46E-17  3.989E+00 -2.291E-15  9.00E-01  2.5484E+00
  750  2.95E-13  2.13E-14  3.01E-19  3.989E+00 -1.379E-15  9.00E-01  2.5484E+00
 1000  7.49E-14  2.57E-15  3.11E-20  3.989E+00 -1.565E-15  9.00E-01  2.5484E+00
 1250  1.84E-14  3.96E-16  3.08E-21  3.989E+00 -1.876E-15  9.00E-01  2.5484E+00
 1500  2.75E-15  4.85E-17  5.63E-22  3.989E+00 -4.524E-15  9.00E-01  2.5484E+00
 1750  3.40E-16  6.61E-18  8.85E-23  3.989E+00 -3.454E-15  9.00E-01  2.5484E+00
 1906  9.91E-17  2.29E-18  4.60E-23  3.989E+00 -1.048E-15  9.00E-01  2.5484E+00

The number of iterations has increased:

• for ns=16 from 272 to 533
• for ns=256 from 1698 to 1906

With the changes on this branch, the symmetric case runs in the same number of iterations
and the asymmetric case runs in 312 (for ns=16) and 1685 (for ns=256) iterations.

TODO:

• find openly-available asymmetric test case and demonstrate that number of iterations goes down (has been tested locally, but I don't know if I am allowed to post the input file here)
• see where else in the code cosmui, ... are used and see if this has implications there

[jonathanschilling]

Further code changes

The following lines needed to be commented out in bcovar.f:

At line 450 in bcovar_par:

IF (lasym) THEN
   tcon = p5*tcon
END IF

At line 881 in bcovar:

IF (lasym) tcon = p5*tcon

Testing with the SOLOVEV case

For the SOLOVEV case listed above, the numbers of iterations until convergence are as follows with PARVMEC=T:

• symmetric
  • ns=16: 272
  • ns=256: 1698
• asymmetric
  • ns=16: 272
  • ns=256: 1711

and with PARVMEC=F:

• symmetric
  • ns=16: 272
  • ns=256: 1697
• asymmetric
  • ns=16: 272
  • ns=256: 1691

Testing with the CTH-like Stellarator test case

I tested this now also with the test.vmec case of the CTH-like Stellarator used for unit-testing PARVMEC.

Here are the results using the current master branch:
PARVMEC=T:

• symmetric
  • ns=15: 1550
• asymmetric
  • ns=15: currently still stalls at ~1e-15

and with PARVMEC=F:

• symmetric
  • ns=15: 1550
• asymmetric
  • ns=15: currently still stalls at ~1e-15

It is noted that in the asymmetric case, the initial R and Z force residuals printed to the screen
are vastly different to the symmetric case:

symmetric case:

  ITER    FSQR      FSQZ      FSQL    RAX(v=0)    DELT        WMHD      DEL-BSQ
    1  2.33E-01  1.72E-02  5.22E-02  7.566E-01  7.00E-01  5.3447E-02  1.000E+00

asymmetric case:

  ITER    FSQR      FSQZ      FSQL    RAX(v=0)   ZAX(v=0)    DELT        WMHD      DEL-BSQ
    1  5.79E-02  4.15E-03  1.30E-02  7.566E-01  0.000E+00  7.00E-01  5.3447E-02  1.000E+00

I then re-ran this case with the changes on fix_cosmui.
First the results with PARVMEC=T:

• symmetric
  • ns=15: 1550
• asymmetric
  • ns=15: currently still stalls at ~1e-15

and with PARVMEC=F:

• symmetric
  • ns=15: 1550
• asymmetric
  • ns=15: currently still stalls at ~1e-15

However, with these, changes, the force residual printout is the same:

symmetric case:

  NS =   15 NO. FOURIER MODES =   41 FTOLV =  1.000E-20 NITER =  25000
  PROCESSOR COUNT - RADIAL:    1  VACUUM:    1

  ITER    FSQR      FSQZ      FSQL    RAX(v=0)    DELT        WMHD      DEL-BSQ

    1  2.33E-01  1.72E-02  5.22E-02  7.566E-01  7.00E-01  5.3447E-02  1.000E+00

  In VACUUM, np =  5  mf =  6  nf =  4 nu = 16  nv =   36
  2*pi * a * -BPOL(vac) =   4.32E-02 TOROIDAL CURRENT =   4.32E-02
  R * BTOR(vac) =  -4.51E-01 R * BTOR(plasma) =  -4.52E-01

  VACUUM PRESSURE TURNED ON AT   53 ITERATIONS

  200  3.19E-08  5.45E-09  2.93E-09  7.440E-01  7.00E-01  5.0582E-02  1.267E-03
  400  3.88E-10  1.09E-10  1.73E-10  7.439E-01  7.00E-01  5.0581E-02  1.236E-03
  600  6.93E-12  1.61E-12  2.65E-12  7.440E-01  7.00E-01  5.0581E-02  1.343E-03
  800  5.27E-14  1.46E-14  3.20E-14  7.440E-01  7.00E-01  5.0581E-02  1.356E-03
 1000  1.26E-15  3.31E-16  8.18E-16  7.440E-01  7.00E-01  5.0581E-02  1.357E-03
 1200  3.11E-17  8.37E-18  1.99E-17  7.440E-01  7.00E-01  5.0581E-02  1.357E-03
 1400  4.70E-19  1.44E-19  2.86E-19  7.440E-01  7.00E-01  5.0581E-02  1.357E-03
 1550  9.62E-21  3.92E-21  5.00E-21  7.440E-01  7.00E-01  5.0581E-02  1.357E-03

 EXECUTION TERMINATED NORMALLY

asymmetric case:

  NS =   15 NO. FOURIER MODES =   41 FTOLV =  1.000E-20 NITER =  25000
  PROCESSOR COUNT - RADIAL:    1  VACUUM:    1

  ITER    FSQR      FSQZ      FSQL    RAX(v=0)   ZAX(v=0)    DELT        WMHD      DEL-BSQ

    1  2.33E-01  1.72E-02  5.22E-02  7.566E-01  0.000E+00  7.00E-01  5.3447E-02  1.000E+00

  In VACUUM, np =  5  mf =  6  nf =  4 nu = 16  nv =   36
  2*pi * a * -BPOL(vac) =   4.32E-02 TOROIDAL CURRENT =   4.32E-02
  R * BTOR(vac) =  -4.51E-01 R * BTOR(plasma) =  -4.52E-01

  VACUUM PRESSURE TURNED ON AT   53 ITERATIONS

  200  3.19E-08  5.45E-09  2.93E-09  7.440E-01 -2.655E-09  7.00E-01  5.0582E-02  1.267E-03
  400  3.88E-10  1.09E-10  1.73E-10  7.439E-01 -2.859E-08  7.00E-01  5.0581E-02  1.236E-03
  600  6.93E-12  1.61E-12  2.65E-12  7.440E-01 -1.629E-07  7.00E-01  5.0581E-02  1.343E-03
  800  5.13E-14  1.46E-14  3.08E-14  7.440E-01 -7.909E-07  7.00E-01  5.0581E-02  1.356E-03
 1000  6.71E-16  3.73E-16  3.24E-16  7.440E-01 -4.109E-06  7.00E-01  5.0581E-02  1.357E-03
 1200  7.93E-15  5.57E-15  3.87E-15  7.440E-01 -1.951E-05  7.00E-01  5.0581E-02  1.357E-03
 1400  1.77E-13  1.25E-13  8.66E-14  7.440E-01 -9.253E-05  7.00E-01  5.0581E-02  1.357E-03
 1600  1.98E-15  4.71E-16  1.24E-15  7.440E-01 -3.165E-06  6.80E-01  5.0581E-02  1.357E-03
 1800  4.45E-15  2.47E-15  1.96E-15  7.440E-01 -1.473E-05  6.80E-01  5.0581E-02  1.357E-03
 2000  9.42E-14  5.27E-14  4.15E-14  7.440E-01 -6.818E-05  6.80E-01  5.0581E-02  1.357E-03
^C

We note that somewhere around iteration 800, the asymmetric case starts to break down.
Also, the time step gets reduced in the asymmetric case, which does not happen in the symmetric case.

[jonathanschilling]

Going a little bit deeper, the issue seems to come from the zero-current algorithm in add_fluxes.
bsupu already has deviations wrt. the symmetric case on the order of 2e-14 in the first evaluation and things get successively worse with every iteration.

[jonathanschilling]

As a test, I switched to constrained-iota with an iota profile fitted to the equilibrium profile obtained from the free-boundary symmetric run:

NCURR = 0
piota_type = "power_series"
ai = 1.27485, -0.366667, 2.85963, -5.95724, 3.0703

However, the convergence issues in the asymmetric case still persist.
Thus, the error cannot be only in add_fluxes.

[jonathanschilling]

cosmui, sinmui are also used in alias to perform the Fourier-space bandpass filtering.
There, a forward Fourier transform of gcon (the constraint force) is performed using only m=1, ..., (mpol-2),
followed by inverse Fourier transform using the reduced poloidal mode number set.

The scaling-down of tcon in bcovar was presumably done as a workaround to the errornous normalization being applied in alias. Now that cosmui, ... are fixed, it thus seems ok to also remove this workaround.

[jonathanschilling]

Other places where cosmui, sinmui, cosmui3 are used:

• alias: ok, integration only performed over ntheta2, side effects of changes fixed by removing tcon = 0.5_dp * tcon

• angle_constraints: used to initialized cos2u, sin2u arrays, but these are not used

• tomnsp*: ok, intended effect

• fixaray: ok, source of changes

• profil3d: ok, used to initialize wint and pwint arrays --> handled by dnorm3 introduced in fixaray

• jxbforce: integrals only over ntheta2; "fix" handled by removing the following line:

  dnorm1 = 2*dnorm1              !SPH012314 in FixAray
  

• wrout: integrals only over ntheta2; "fix" handled by removing the following line:

  !Changed integration norm in fixaray, SPH012314
           tmult = 2*tmult
  

This should conclude the changes in other parts of VMEC.
A unit test should be set up now to verify that when VMEC is run for a symmetric input,
the outputs of the symmetric code path (lasym=F) agree with the outputs from the asymmetric code path (lasym=T).

[jonathanschilling]

The following two sketches illustrate how the inverse Fourier transforms of the geometry to real space in totzsp* are setup
in the asymmetric case (with symrzl to extend quantities to the full poloidal interval in the asymmetric case):

The following two sketches illustrate how the forward Fourier transforms/integrals of the forces to Fourier space in tomnsp are setup in the asymmetric case (with symforce to decompose the forces on the full poloidal interval into definite-parity contributions on the reduced poloidal interval):