Changing the calculation of the normal curvature

[acerfon]

After running some tests, I still think the currently implemented version of the normal curvature is inaccurate in DESC, particularly for vacuum equilibria, where DESC is not particularly good at getting the current density to zero.

If you use other definitions of the normal curvature, including the one with the gradient of the pressure, you get much better consistency at low beta and for vacuum equilibria between the pure vector algebra versions of the quantity and the 'physics'-based versions.

[acerfon]

The issues are deeper than I expected. I am closing the issue to avoid falling into a rabbit hole. Accuracy is very limited for many of these calculations. I guess DESC is not great at calculating derivatives of B, and we just have to live with that weakness, as it is the case for any 3D equilibrium code.

[rahulgaur104]

Hi Antoine, Can you tell us which equilibria you tested or give us a minimal failing example?

Accuracy is very limited for many of these calculations. I guess DESC is not great at calculating derivatives of B, and we just have to live with that weakness, as it is the case for any 3D equilibrium code.

What you are saying may be true for equilibria with a large force error (jxB-grad p) but I don't see why it would be true in general.

[acerfon]

You may just take Daniel Dudt's poloidal.h5 NFP=1 QI and check the different formulas for the normal curvature in Per Helander's review paper. You will not even get 3 significant digits of agreement depending on which formula you are using. When I get a chance, I will cook up other equilibria where you barely have 1.5 digits of accuracy.

However, as I said, I would not worry about that issue if I were you, because like I said, it's just that in DESC, like in any MHD equilibrium code with nested surfaces, the force error jxB-grad p is pointwise not good, 2-3 digits at best. And that hurts for quantities that require derivatives of B.

[acerfon]

Here is the code, if you want to see the bad accuracy for yourself:

filename = "poloidal.h5" eq = desc.io.load(filename)[-1] nfp = eq.NFP s = .5 Mreg = 12eq.M Nreg = 12eq.N rho = jnp.sqrt(s) surf_grid = LinearGrid(rho=rho,M=Mreg,N=Nreg,NFP=nfp,sym=False)

modB = eq.compute("|B|", surf_grid)["|B|"][0] gradB = eq.compute("grad(|B|)", surf_grid)["grad(|B|)"][0] gradpsi = eq.compute("grad(psi)", surf_grid)["grad(psi)"][0] modgradpsi = eq.compute("|grad(psi)|", surf_grid)["|grad(psi)|"][0] gradp = eq.compute("grad(p)", surf_grid)["grad(p)"][0]

Bvec = eq.compute("B", surf_grid)["B"][0] gradBmat = eq.compute("grad(B)", surf_grid)["grad(B)"][0] bunit = Bvec/modB gradbmat = (gradBmat - np.outer(gradB,bunit))/modB kappasecondtest = np.dot(np.dot(bunit,gradbmat),gradpsi/modgradpsi) # Compute normal curvature in the geometric way, as another check print(kappasecondtest)

kappatest = np.dot(gradB/modB+4np.pi1e-7*gradp/modB**2,gradpsi/modgradpsi) # Formula using MHD equilibrium print(kappatest) print(eq.compute("kappa_n", surf_grid)["kappa_n"][0]) # DESC way of determining curvature

Like I said, depending on the different formulas, you have agreement up to at most 2 significant digits.

[rahulgaur104]

Ok, thanks. I'll probably rerun the equilibrium at a higher resolution using solve_continuation_automatic (instead of eq.solve) and than recalculate kappatest again.

[acerfon]

Sounds good. I am curious about your results. As I said, it may not be DESC's fault, but the overall way people have to approach MHD equilibria with nested flux surfaces.

The question, then, would be what do you trust more, curl of B or grad B/B and grad p/B^2 ? My inclination would be for the latter, since we know that in ideal MHD with nested surfaces, curl B is a very funky object, with bad math behavior. And somehow, in DESC, you have decided to trust curl B better for the computation of the normal curvature.

[rahulgaur104]

And somehow, in DESC, you have decided to trust curl B better for the computation of the normal curvature.

Yeah, that's probably not a great choice. However, for stability and gyrokinetic calculations, we actually use the grad(mu_0 p + B^2/2) dot grad psi. I agree with your suggestion and I'll change the calculation of kappa_n in DESC.

[acerfon]

I am not sure I would call my comments a "suggestion", it was more of an open question, which I thought was worth studying. For some equilibria and some flux surfaces it does not make much of a difference. For other equilibria or other flux surfaces, it seems to matter. I really do not know what is best. The expression b\cdot\nabla b \cdot n seems to agree with the present formula DESC is using, but tends to disagree with the formula involving the pressure...