Finish VMEC I/O

[f0uriest]

Currently we can read VMEC inputs (i think it automatically detects?), but it would be good if we could also read/write vmec outputs and interpolate between the two bases (for comparisons, and for using DESC outputs in places that expect VMEC format)

To Do:

• [x] fix bug in input conversion
• [x] load vmec solution as initial guess (ie convert vmec to desc)
• [x] save desc solution as vmec .nc format

[ddudt]

Specifically, we want to output a wout_FILENAME.nc file that follows the same format as VMEC so we can pass it to other codes that rely on VMEC solutions

[ddudt]

We also want to be able to load an initial guess from a VMEC output

[f0uriest]

The utility to convert a vmec input file to desc gives the wrong boundary, as seen here:

However, if we load a previously solved vmec case, that seems to work fine:

[f0uriest]

However, when loading a previous vmec output (ie to use as an initial guess) it sets the pressure/iota profiles based on the spectral coefficients saved in the cetcdf, which may be zero if the vmec equilibria was solved with a current profile instead of iota. Would it be possible to read the discrete iota values calculated in that case and fit them to a power series? otherwise using the vmec solution as an initial guess will be weird because as far as desc knows, the equilibrium has no rotational transform

[ddudt]

The input conversion bug should be fixed by 6a0a474e341dc252253fe4e7bf132f1a922f115c.

Keeping this issue open because we still need to check the output conversion from a DESC solution to the NetCDF file in VMEC format. We have solutions saved in this format, but haven't tested if those files run with external programs that expect a VMEC solution as input.

[ddudt]

Saving a DESC solution in a VMEC-compatible NetCDF format has now been added and tested.

[patil3]

According to this discussion, it seems that desc should be able to use a vmec output file (.nc) as an initial guess. However, when I tried this with desc -o dshape --vmec wout_DSHAPE.nc -vv DSHAPE, I got an error (see disp.txt). disp.txt Am I using this functionality incorrectly? I have attached my input file 'DSHAPE' and the guess file 'wout_DSHAPE.nc' is the same as the one in the repository. DSHAPE.txt

[f0uriest]

You do seem to be using it correctly, I get the same error on master. It seems like there's some inconsistency in how the constraints are passed around. I cleaned up a lot of that in #107 and using vmec as an initial guess seems to work on that branch which should be merged in soon. One other thing to be aware of when using vmec as an initial guess is that it still uses the continuation parameters from the DESC input file, so if you are starting from a vacuum solution but the vmec one is finite beta it might not be a very good guess.

[patil3]

WIth the geometry branch, the code works as expected with vmec output as the initial guess. There's just one small issue that I wanted to bring up: when such a guess is specified, for some reason, the code uses L_grid = L, M_grid - M, N_grid = N. When a guess is not specified, L_grid etc. are used correctly from the input file. I have attached guess.txt and noguess.txt which show the respective outputs for the input file DSHAPE.txt. The guess I used is wout_DSHAPE.nc from repository. DSHAPE.txt guess.txt noguess.txt

[patil3]

I was able to fix this. In the original version, in main .py, eq_fam[0] is created from VMECIO.load which created an equilibrium object with L_grid = L etc. When eq_fam[0] = ir.inputs[0] is done, it changes the eq_fam[0].inputs but not eq_fam[0].L_grid etc. We need the equilibrium object eq_fam[0] to have attributes as per the input file. I think a much better way to read equilibrium would be to read the VMEC output into a temporary equilibrium eq_from_vmec and then assign only its solution data to eq_fam[0]:

    if ir.args.vmec is not None:
        # copy equilibrium from vmec output
        eq_from_vmec = VMECIO.load(
            ir.args.vmec,
            L=ir.inputs[0]["L"],
            M=ir.inputs[0]["M"],
            N=ir.inputs[0]["N"],
            spectral_indexing=ir.inputs[0]["spectral_indexing"],
        )
        # assign only the solution data to equil_fam[0]
        equil_fam[0].R_lmn = eq_from_vmec.R_lmn
        equil_fam[0].Z_lmn = eq_from_vmec.Z_lmn
        equil_fam[0].L_lmn = eq_from_vmec.L_lmn
        equil_fam[0].x0 = eq_from_vmec.x0
        equil_fam[0].x = eq_from_vmec.x

I used this and it seems to have resolved the issue.

[patil3]

When a solution is saved using VMECIO.save and then it's read using VMECIO.load(), it seems to have large force imbalance errors. This may be crucial for using VMEC solutions (.nc) as initial guesses. Since starting every run from scratch is time-consuming for non-axisymmetric cases, I tried saving an already obtained DESC solution as an .nc file and then using it as a guess using the --vmec flag. But it seems that the bug in VMECIO prevents DESC from using the guess effectively. I am currently working on using the desc .h5 file itself as a guess which should make resuming/restarting/using a guess easier. To reproduce an example (files attached): DSHAPE_vmecio_bug.zip

desc DSHAPE: solve an equilibrium In python terminal: eq1=EquilibriaFamily.load(load_from='DSHAPE_output.h5')[-1]: load the solution plot_section(eq1,"|F|", norm_F=True, log=True): check the actual error

VMECIO.save(eq1,'DSHAPE_output.nc'): save as .nc eq2=VMECIO.load('DSHAPE_output.nc'): load from the saved .nc VMECIO.plot_vmec_comparison(eq1, 'DSHAPE_output.nc'),VMECIO.plot_vmec_comparison(eq2, 'DSHAPE_output.nc'): check if the surfaces match

plot_section(eq2,"|F|", norm_F=True, log=True): check the error for the solution loaded from .nc

[f0uriest]

@ddudt can probably give better insights, but I think this is somewhat expected for two reasons:

1) the default for VMECIO.save is only 128 flux surfaces which is pretty low. It's enough that the flux surfaces look the same by eye but force balance can still be pretty high. Increasing the number of surfaces when saving might help.

2) In general, the conversion from a desc solution to a vmec one is "exact" since we can evaluate all the vmec quantities at any point in the domain. However, going from vmec to desc requires fitting discrete data to global basis functions, and the surface spacing that vmec uses isn't very good for fitting zernike polynomials, so there will also be some error there. Adding more surfaces can help a bit, but because of the spacing vmec uses the least squares solution becomes poorly conditioned if you go too high and you'll end up with numerical issues in the fit

[ddudt]

Reason 2 that Rory gave is the main issue. Converting a DESC solution into the VMEC format is exact* because of our spectral basis, but fitting a VMEC solution into that basis introduces numerical errors. There isn't much we can do to improve the fit because the VMEC solutions are always given on surfaces linearly spaced in toroidal flux, and those quadrature points result in a poorly conditioned system. The VMECIO class was primarily intended to allow DESC equilibria to be used with external codes that expect a VMEC input. Loading VMEC solutions should be used with caution.

*The conversion is exact only for R, Z, and lambda. The coefficients of other "derived" quantities such as the magnetic field components must be Fourier transformed from the "real space" data on each surface. This transform is done with a higher resolution than the primary quantities and seems to be very accurate.