Negative Jacobian

[bendudson]

Can we modify the coordinate system slightly, so that the Jacobian is always positive?

BOUT++ currently rejects inputs with negative Jacobians, but this is what happens when the poloidal field is negative.

In https://bout-dev.readthedocs.io/en/latest/user_docs/coordinates.html can we just remove the factor of sigma (sign of Bpol) from the definition of Z? This would reverse the sign of g^xz and g^yz, change some elements of g.. and mean that J was always positive. It would also mean that:

• x is always outwards from core to SOL
• y is always clockwise in the right hand cut-through of the tokamak
• z is always toroidal angle, anti-clockwise looking from the top

[johnomotani]

I think something like this would be nice to do. I vaguely remember not being sure that all the \sigma_B\theta factors were in the right place, but running out of time to check carefully, when I was looking at the coordinate system derivations to implement hypnotoad2 originally.

Off the top of my head - one way that might be nice to check would be to generate a grid file, then flip the sign of psi and the toroidal field (f_pol?) in the input equilibrium and generate a second grid file. I think then the metric coefficients should be identical (because the magnetic field lines are in exactly the same place, just pointing in the opposite direction), but the curvature vector (or Curl(b/B)) should flip sign. Are there any other signs that ought to change?

[johnomotani]

It might just be a matter of definition, but what sign should Grad_par have? If y is always clockwise but B_pol can change sign, then e_y can be +b or -b. If Grad_par(f)=b.Grad(f) then it would need a sign(B_pol) factor, but would it make more sense to define it as Grad_par(f)=e_y.Grad(f) even if that might be -b.Grad(f)? I guess this one should always cancel with something like the sign(B_pol) factor on a V_par (that would be there for the same reason).

What about the sign of bracket(phi,f)=b.Grad(phi)xGrad(f)/B ?

[bendudson]

Oh no... I think we have to keep Grad_par(f)=b.Grad(f) or we'll go crazy. So then we would have to reverse the direction of the y coordinate in the grid file? Good idea about generating equilibria with reversed B and checking signs. Might even be possible to automate it for a test

This has become more urgent for me because a certain ball-like energy producing design has negative poloidal field...

[johnomotani]

I don't know. I like the idea of keeping the y-coordinate pointing in a fixed direction (like you proposed at the top), but allowing e_y=-b seems likely to be a pain. On the other hand, if the y-coordinate has to go in the opposite direction, we'd have to start the grid at the outer target instead of the inner one - I don't know how much of the code would need a special case to handle that.

One possible work-around that we discussed at one of the STORM meetings some time ago is to rotate the whole equilibrium upside down when running the simulation by messing around with the equilibrium input, and then rotate the outputs back 180 degrees when plotting them. Downside being that this option requires faffing around by the user both when before generating the grid file and on all the outputs.

Looking at this thread, from the point of view of BOUT++ (which is by far the most lines of code) it seems like the simplest (therefore best) thing would be to require e_y=+b always, so the y-coordinate has to go around anticlockwise on the right-hand cut through the tokamak if B_p does. I think we still need to check the signs in the coordinates derivation carefully, because this definition isn't the one that that derivation uses as it's written, I'm just hoping it will all work out OK and not too complicated*. I'd hope the changes to hypnotoad would be confined to the TokamakEquilibrium class, so at worst we'd have two cases to handle through the whole thing and end up doubling the length (with mostly copy-paste) of a ~1500 line file.

* I'm not sure the best definition to start from - one possibility is to say (psi,theta,zeta) is right-handed, but might have theta-clockwise/zeta-anticlockwise or theta-anticlockwise/zeta-clockwise (depending on sign(B_pol)). Or alternatively, we fix theta-clockwise/zeta-anticlockwise and include sign(B_p) in the definitions of both(?) y and z.

[bendudson]

It does seem like e_y = +b might be easiest; this way I think resolves the problem of having to work out whether sqrt(g_22) should be positive or negative. Perhaps it's possible to have this as a small output post-processing step for now, to reverse the y coordinate and flip sign of J?

[johnomotani]

Perhaps it's possible to have this as a small output post-processing step for now, to reverse the y coordinate and flip sign of J?

That could work nicely! One thing to bear in mind - I wouldn't swear that I haven't assumed sign(B_p)>0 in some of the calculations in hypnotoad2. It was something I originally wanted to be consistent with (so there is a signbp included in various places), but couldn't follow through after giving up on checking the sign(B_p) in the coordinates derivation. So watch out for the signs of the metric coefficients and curvature!

[bendudson]

Thanks! Yes there are a lot of traps to fall into...the integrated shear I for example should reverse sign if Bp<0.

[bendudson]

I've made a start, going through the coordinate system derivation in the BOUT++ manual: https://github.com/boutproject/BOUT-dev/pull/2165