The Jacobian orientation test is not invariant under the poloidal rotation applied to an asymmetric boundary

[CharlesCNorton]

An asymmetric boundary is rotated in the poloidal angle before the orientation of the Jacobian is tested. The rotation is at readin.f90:218,

delta = ATAN2(rbs(0,1) - zbc(0,1), rbc(0,1) + zbs(0,1))

which puts the boundary into the representation with RBS(m=1) = ZBC(m=1), and the test is at readin.f90:367,

rtest = SUM(rbcc(1:ntor1,mj))
ztest = SUM(zbsc(1:ntor1,mj))
lflip=(rtest*ztest .lt. zero)

rtest*ztest is not invariant under that rotation. Under theta -> theta + delta the m=1 coefficients mix as

rbcc' = rbcc cos + rbsc sin      zbsc' = zbsc cos - zbcc sin
rbsc' = rbsc cos - rbcc sin      zbcc' = zbcc cos + zbsc sin

so the same boundary in two gauges can give products of different magnitude, of different sign, or zero. The combination that is invariant is the m=1 determinant rbcc*zbsc - rbsc*zbcc, which reduces to rtest*ztest when rbsc and zbcc vanish, that is for a stellarator-symmetric boundary.

Reproduction

Relabel input.up_down_asym_current by theta -> pi - theta. With ntor = 0 that is

RBC(m) ->  (-1)^m RBC(m)     RBS(m) -> -(-1)^m RBS(m)
ZBS(m) -> -(-1)^m ZBS(m)     ZBC(m) ->  (-1)^m ZBC(m)

The surface is unchanged as a point set and only the poloidal direction reverses, so the orientation test must fire. Printing delta, rtest, ztest and lflip:

input	delta	rtest	ztest	lflip
as given	0.4636476	0.80498447	0.53665631	F
relabelled	1.5707963	0.60000000	3.6739e-17	F

The rotation drove ztest to zero, so the product is zero, the test cannot fire, and the wrongly oriented boundary is passed to the solver with no m=1 sine component in Z. That run takes 7 Jacobian resets and stops at the 2001 iteration cap with fsqr = 2.07e-02, against 1311 iterations to fsqr = 5.07e-12 for the same surface as given.

That the rotation is the cause

Forcing delta = 0 makes both inputs converge in the same 1519 iterations to the same fsqr = 3.588e-12, and their wout files agree to 1e-8 in every variable.

The same relabelling of a stellarator-symmetric case is absorbed exactly as it should be: input.solovev and its relabelled twin both converge in 204 iterations to the same fsqr = 9.36e-13, and the wout files differ only in the sign of the poloidally oriented quantities, iotas, iotaf, chi, chipf, q_factor, buco, bsubumnc, bsupumnc, bsubsmns, jcurv, jdotb, currvmnc and ctor. So flip_theta is right and the fault is in the test that decides whether to call it.

PARVMEC readin.f:663 and VMEC++ boundaries.cc:244 carry the same test after the same rotation. VMEC++ already records the doubt next to it: a TODO proposing the signed polygon area of the boundary instead, of which the m=1 determinant is the leading term.