Loading geometries

[theo-brown]

I'm aware that the geometry portion of the code is under active development, with work being done to move towards allowing time-dependent geometries.

However, in the mean time I am trying to get to grips with how to load geometries, supported file types, requirements, etc.

1. What is the .mat2cols file format and how is it constructed from a CHEASE equilibrium?

My original guess was that it was an ASCII-formatted file output by running save(...., '-ascii'); from within MATLAB; this seems to be backed up by the PINT commit where it was added.
However, I can't load it back in to MATLAB using load(..., '-ascii') due to the presence of the header line, so I would have to write a script to load it in to a structure in the correct format. Presumably this already has been done by TORAX/PINT developers somewhere!

🎯 If there are custom scripts used in creating the mat2cols from CHEASE/MATLAB, it would be great if they could be provided given that the most mature support for equilibria in TORAX is CHEASE.

2. What is the best way (currently) of loading a non-CHEASE equilibrium?

Answering my own question here for completeness: I believe it is to create a StandardGeometryIntermediates and then pass it to build_standard_geometry.

🎯 It might be neat to have more constructor classmethods for StandardGeometryIntermediates that load different common equilibrium formats, such as GEQDSK. This might involve writing code for manual flux-surface-averaging, or other methods for removing the dependency on CHEASE.

3. What is the planned future best way of loading a non-CHEASE equilibrium, if different?

Open discussion on this point, really. Is the above suggestion the best method for loading geometry?

[theo-brown]

Ok, from digging around the CHEASE repo, it looks like the tool for creating the file is either the Perl script o.chease_to_cols or the MATLAB script cheasedata_mat_to_cols.m. I presume the example equilibrium file provided by TORAX is created by the latter. It's a custom ASCII/txt format.

If you have a CHEASE .mat file, you can use the following commands to turn it into a .mat2cols from Python:

import h5py
import numpy as np


def chease_mat_to_mat2cols(input_file, output_file):
    with h5py.File(input_file, "r") as h5:
        labels = []
        label_references = h5["chease_data"]["profiledata"]["labels"][0]
        for label_reference in label_references:
            labels.append("".join(chr(c.item()) for c in h5[label_reference][:]))

        np.savetxt(
            output_file,
            h5["chease_data"]["profiledata"]["data"][:].T,
            header=" ".join(labels),
            comments="% ",
            fmt="%13.5e",
        )

[theo-brown]

Great point about IMAS. I'm not super familiar with the equilibrium IDS, would it include flux-surface-averaged quantities? If so, I'll probably jump over GEQDSK and head straight to IMAS.
If not, that means that flux-surface-averaging will still need to be performed somewhere between whatever equilibrium code and TORAX, and we can try and hunt for an existing tool for that.

[theo-brown]

I'm aware of Ben Dudson's FreeQDSK, but that's just for reading/writing.

PSFC has eqtools which I think might support some sort of flux-volume-average.

Not entirely sure - @allen-adastra have you used this / does it still exist?

[jcitrin]

megpy may also be useful

[allen-adastra]

I'm aware of Ben Dudson's FreeQDSK, but that's just for reading/writing.

PSFC has eqtools which I think might support some sort of flux-volume-average.

Not entirely sure - @allen-adastra have you used this / does it still exist?

I did not know about this, thanks for sharing!

[theo-brown]

Sorry for the @, thought you might be aware as it's PSFC!

[fcasson]

For solving the transport equations, there is a shortcut to provide the 1D geometric g0,1,2,3 quantities directly from JETTO and avoiding eqdsk or CHEASE (all transport solvers need the same things here)

The proof for the relation <grad(rho)> = SURF/DVEQ can be obtained from equations 3 and 5 in the ASTRA manual. If f is set to grad(V) in equation 3:

the relation in equation 5 is obtained:

which is valid for any magnetic flux surface coordinate a. If a is set to rho, one obtains S = SURF = dV/drho * <|grad(rho)|> = DVEQ * <|grad(rho)|>.

When we get to models that need 2D geometry information (transport models like TGLF and source models like GRAY), g0,1,2,3 will not be enough. I'm not yet clear if TORAX already carries some 2D information around in order to provide it to other models (this goes beyond the paper), but it seems to load more than g01,2,3 from mat2cols? Some more definitions will be needed to write 2D interfaces, and here it might make sense to restrict interfaces to IMAS definitions.

[jcitrin]

Indeed, if you already have a set of solved equilibria from ESCO in JETTO and can utilize jsp output to get what you need, that is a viable path. Would need to write a jsp loader in geometry_loader and from_jsp classmethod in StandardGeometryIntermediates in geometry.py, etc. Ultimately coupling to the IMAS equilibrium IDS is the solution, so writing a jsp loader will be a bit of a duplication. On the other hand it's something that can be done right now, and will still be useful for cross-verification with the future IMAS equilibrium coupling.

Note that part of the code is a moving target at the moment since we're in the middle of generalizing to time-dependent equilibria.

Regarding 2D geometry information. It isn't currently carried around but can extend when needed.

[theo-brown]

I think much of the hard part of the JSP loading stuff exists already in jetto-pythontools, for example, which means it might not be too tricky in the short term.

[theo-brown]

It's a bit tricky to decipher which of the geometry parameters are needed. Some of them are CHEASE outputs that only are used to derive a different geometry parameter, that another code might compute directly.

Other than the g_0 to g_3 parameters listed in the paper, what geometry terms are required and what are their definitions?
Their in-code names are visible here but it's a bit hard to decipher some of them.

[jcitrin]

StandardGeometryIntermediates has everything that's needed from CHEASE which should have analogues in other codes.

What is needed in TORAX in the end is

• rho (toroidal flux coordinate)
• g0,g1,g2,g3
• F (standard definition of Bphi*R flux function)
• volume (array of flux surfaces volumes)
• area (array of flux surfaces areas)
• delta (standard definition of triangularity, can also be from averaged delta_upper and delta_lower)
• Rin (array of min R on flux surfaces)
• Rout (array of max R on flux surfaces)
• Rmaj (geometric major radius)
• B0 (vacuum on-axis B-field)

Everything else can be calculated from them. I may have forgotten something, but it's all in build_standard_geometry

[theo-brown]

Thanks for extracting this! Got a bit overwhelmed by the amount of terms flying around in build_standard_geometry, so it's very helpful to have it clearly spelt out.

[fcasson]

Just to be absolutely sure about different axis locations...

Rmaj = geometric axis for given surface = (Rin + Rout) / 2
B0 = B(Rmaj(rho=0)) ??

[jcitrin]

Correct. Rmaj is the major radius at the vacuum field magnetic axis (geometric axis). B0 is the vacuum magnetic field on that axis.

[theo-brown]

TORAX/IMAS/JETTO geometry terms:

Variable	IMAS IDS	JETTO
rho	equilibrium.time_slice[:].profiles_1d.rho_tor	JSP/XRHO
g0 =	equilibrium.time_slice[:].profiles_1d.dvolume_drho_tor * equilibrium.time_slice[:].profiles_1d.gm7 †	JSP/SURF †
g1 = <grad V^2>	equilibrium.time_slice[:].profiles_1d.dvolume_drho_tor^2 * equilibrium.time_slice[:].profiles_1d.gm3 †	JSP/DVEQ^2 * JSP/GROQ †
g2 = <grad V^2/R^2>	?	JSP/K
g3 = <1/R^2>	equilibrium.time_slice[:].profiles_1d.gm1	JSP/A
F	equilibrium.time_slice[:].profiles_1d.f	JSP/F
flux surface volumes	equilibrium.time_slice[:].profiles_1d.volume	JSP/VOL
flux surface areas	equilibrium.time_slice[:].profiles_1d.surface	?
delta (triangularity)	(equilibrium.time_slice[:].profiles_1d.triangularity_upper + equilibrium.time_slice[:].profiles_1d.triangularity_lower)/2	JSP/TRI
Rin	equilibrium.time_slice[:].profiles_1d.r_inboard	JSP/RI
Rout	equilibrium.time_slice[:].profiles_1d.r_outboard	JSP/R
Rmaj	?	?
B0	?	?