Measure the discretization by the method of manufactured solutions

[CharlesCNorton]

Reference tests compare against Fortran VMEC, which cannot say whether either code converges to the continuum problem. This measures it: an analytic mapping is made the exact solution of a modified problem by installing the negative of its continuum ideal-MHD force as a source, so the distance from the solver to the mapping is discretization error.

Change

IdealMhdModel::SetForceSource adds a fixed spectral force to the decomposed force, before the m = 1 gauge rotation. It is refused unless the input sets the new enable_force_source, and it is reinstalled after every InitializeRadial, which an axis reguess and each multi-grid step trigger. VmecModel.set_state_from_fourier goes through FourierGeometry::InitFromState, the routine a hot restart already uses, so the basis conventions keep one implementation.

Measured

examples/manufactured_solution.py: helical axis, rotating ellipse, non-trivial lambda; nfp = 3, aspect ratio 18, beta(0) = 2 per cent, mpol = 5, ntor = 3, 24 x 24. Errors are relative L2 over s >= 0.1.

ns	energy	R, Z force	lambda force	solved R, Z	solved lambda
25	2.858e-08	2.590e-04	1.715e-05	7.980e-05	1.755e-03
51	6.586e-09	5.866e-05	3.999e-06	1.116e-05	1.810e-04
101	1.646e-09	1.426e-05	9.586e-07	2.463e-06	7.014e-05
201	4.116e-10	3.525e-06	2.347e-07	4.593e-07	1.614e-05
order	2.00	2.02	2.03	2.42	2.12

Angular error falls by 9.4 for every two grid points, from 3.2e-06 at 16 x 16 to 4.2e-11 at 26 x 26. The source drops the force at the mapping by 5.4e+04 at ns = 25 and 7.9e+06 at ns = 201, and the solve reaches fsq = 1.5e-16 at every resolution.

The m = 0 lambda is the exception. Its force is second order, 8.5e-09 at ns = 201, while its state error sits at 1e-02 and does not fall, so this is conditioning rather than consistency. FourierGeometry::extrapolateTowardsAxis copies that component from the first interior surface onto the axis, and the force transform does not transpose the copy.

Asymmetric and constrained current

The mapping gains a sine series in R and a cosine series in Z and lambda at m <= 1, obeying RSC_{1n} = ZCC_{1n}: 2.01 for the R, Z force and 2.04 for lambda, with the source cancelling by 4.8e+07 at ns = 201. The symmetric mapping run through the asymmetric path reproduces the symmetric energy error exactly at every resolution.

With ncurr = 1 the solver solves chi' from currH(s) = <B_u>(s) instead of reading an iota profile. Prescribing the mapping's own <B_u> as a cubic_spline_i profile returns its iota and leaves the order unchanged at 2.01.

Free boundary

only_coils takes the vacuum field from the mgrid with no Laplace solve, and a field bilinear in (R, Z) is interpolated exactly, so the vacuum pressure is analytic. At the last radial node only the half-grid cell inside the boundary exists, so the force there is (1/ds) dw/dx_s + (1/2) EL(w) + edge, whose first and third terms are the pressure jump (1/ds) R (p_total - p_vac) (z_u, -r_u).

ns	volume	last node	fsq with source
25	4.034e-05	1.837e-05	3.463e-11
51	9.319e-06	4.232e-06	4.794e-12
101	2.277e-06	1.058e-06	1.542e-12
201	5.758e-07	2.645e-07	1.134e-12
order	2.00	2.00

The edge condition is second order, like the volume. There is no free-boundary counterpart to the solved-state column above: the free-boundary descent limit-cycles around fsq = 1e-09 for 40000 iterations, where the fixed-boundary one reaches 1.5e-16, so the state it stops at is not a solution of the modified problem.

Tests

Thirteen Python tests in 8 seconds pin the orders, the rejections and the free-boundary edge; //vmecpp/vmec/vmec:force_source_test holds the source's effect on the assembled force to within 1e-12 of the source itself across 1, 2, 4 and 8 radial threads, which a source landing on the wrong surfaces would miss by its own size. sympy is a test extra and the module skips without it; scipy is needed only by --study fit.

[CharlesCNorton]

fixing ci and other things