Creating an equilibrium file

The computational grids of both B2.5 and EIRENE follow the shape of the flux surfaces, so a SOLPS-ITER simulation will require a tokamak magnetic equilibrium reconstruction as basis. This document explains how to extract the required equilibrium file (.equ) from EFIT output in an arbitrary COMPASS discharge. COMPASS-Upgrade equilibrium reconstructions are analogous, only they are not loaded from CDB (COMPASS DataBase) but from the CUDB (COMPASS Upgrade DataBase).

In short, you procure an EQDSK (read equidisk) file either by conversion via PLEQUE, or by running your own EFIT reconstruction via a script. Then you copy this file into your SOLPS-ITER simulation folder. Finally, you convert it into the desired .equ format using a simple command. Thankfully, it’s quite straightforward.

Many thanks to Lukáš Kripner (kripner@ipp.cas.cz), who created the PLEQUE Python package for EFIT equilibrium reconstruction manipulation and visualisation. Without it, these tutorials wouldn't be possible.

B2.5 and EIRENE grids

B2.5 and EIRENE grids, aligned to the flux surfaces.

Before you start

The files we will be working with do not contain any information but the equilibrium itself. Notably, they tell you nothing of which shot they’re from or which time they’re from. Believe me, it isn't fun to search your command history desperately trying to find which discharge phase you’re modelling. Be sure to keep this metadata elsewhere! File names work well, for example.

Load an equilibrium reconstruction from the C(U)DB

Log into Soroban.

ssh username@soroban # add tok.ipp.cas.cz if you are outside the IPP network

Load the modules compass/2021 (Python, CDB access and other) and pleque/master.
```
module load compass/2021 pleque/master
```
In case of module clashes, try module purge first.
Open an iPython terminal.
```
ipython
```

Import the cdb function from PLEQUE and load the desired equilibrium.

from pleque.io.compass import cdb
eq = cdb(17588, 1100)  # shot, time in ms

This code loads the standard CDB equilibrium reconstruction. This is the reconstruction that is loaded by default and which is, presumably, usually used for experimental data analysis and processing. There can, however, be multiple reconstruction variants present for a single discharge, described at the COMPASS wiki. You can view them in the WebCDB. For instance, there are 5 reconstruction variants available in the discharge #17588:

default: standard CDB reconstruction
EFIT_HRTS: reconstruction with the plasma pressure taken as $p=2p_e$ from Thomson scattering measurements, see below)
EFIT_LCFS_2PM: reconstruction constrained by the LCFS position calculated using the two-point model
V4_std_O: reconstruction with corrected IPR coils positions
v7_std_hires: standard CDB reconstruction with high resolution

To load a particular reconstruction variant using PLEQUE:

eq = cdb(17588, 1100, variant='EFIT_HRTS')

To access COMPASS Upgrade equilibrium reconstructions:

from pleque.io.compass import cudb
eq = compass.cudb(24300, 1.5)  # scenario number, time in s

Exploring the equilibrium with PLEQUE

The eq object is the central entity of PLEQUE. Try this with it:

# Plot an overview of the equilibrium
eq.plot_overview()

# Plot the separatrix, LCFS (last closed flux surface) and the first wall
eq.separatrix.plot()
eq.lcfs.plot()
eq.first_wall.plot()

# Obtain the strike point and magnetic axis position
R_strike_points = eq.strike_points.R
Z_magnetic_axis = eq.magnetic_axis.Z

More examples can be found in the PLEQUE documentation.

To save a PLEQUE equilibrium in the EQDSK format:

eq.to_geqdsk('standard_CDB_reconstruction_16908_1100ms.eqdsk')

The .eqdsk extension is not needed since it's just a simple text file, but it is instructive. If you need finer or coarser grid, use the parameters nx and ny. Their default values are nx=64 and ny=128.

Reconstruction variant: corrected IPR coils positions, flux loops and Mirnov loops

For this reconstruction variant, we thank Ondřej Kovanda, who has sadly left IPP in search of greener pastures.

As elaborated in [Kovanda 2019] and [Jirakova 2019], the standard CDB reconstruction is provably inaccurate because one of its main inputs, magnetic field measured by the Inner Partial Rogowski (IPR) coils, is not measured where the standard CDB reconstruction thinks it's measured. The IPR coils are very slightly misaligned, and correcting this misalignment results in demonstrably better equilibrium reconstructions. Unfortunately, it also results in somewhat worse convergence of the EFIT equilibrium reconstruction algorithm, which prevented its system-wide deployment as the new and better standard CDB reconstruction. Nevertheless, the reconstruction variant with corrected IPR coils positions was calculated en masse and is available as the v4_std_O variant.

Another reconstruction variant that Ondra was working on concerned adding more constraining input to EFIT, in particular the measurements of flux loops and divertor Mirnov coils. Flux loops don't constrain the solution very strongly and the divertor Mirnov coils data is extremely noisy, but they are present in most discharges and can be used for equilibrium reconstruction. The algorithm fails even more often than only with corrected IPR coils positions, but it can also result in a much more accurate equilibrium reconstruction, especially in the divertor region thanks to the divertor Mirnov coils.

With the following tutorial, you can make your own equilibrium reconstruction with Ondra's corrections. It is based upon the EFIT interface page from the COMPASS wiki, written by Ondra himself.

Log in to Abacus:
```
ssh -oHostKeyAlgorithms=+ssh-rsa -X username@abacus01
```
You might be prompted for your password.

Load the compass/dev module and set the working directory for output of equilibrium:

module load compass/dev
cd /compass/Shared/Common/IT/projects/EFIT_interface
export EFIT_RESULTS_PATH="/compass/Shared/Users/jirakova/"

In case of modules clash, try module purge first.

Run the interface script efit.py, described on the COMPASS wiki.
```
./efit.py 13929 1.150 -t -e -x --floff --mirnoff
```
This creates a reconstruction of the discharge #13929 at $t=1.150$ ms with corrected IPR coils positions. It should be the same one as the v4_std_O variant stored in the CDB. The arguments meaning is the following:
- {shot number} - selects appropriate shot
- {time in seconds} -t - selects the time in the shot
- -e - create EQDSK file as an output
- -x - do not generate an EFIT video (we do not need it)
- --floff - do not use the flux loops as reconstruction input
- --mirnoff - do not use the Mirnov coils as reconstruction input

Important note: By default, the script uses data from flux loops and Mirnov coils as constraining input for the reconstruction. However, in most shots these two diagnostics don’t produce reliable signals, so their inclusion doesn’t improve the reconstruction quality.

This will create a directory 13929 inside $EFIT_RESULTS_PATH, which contains the subdirectory run_0 (or run_1 etc. if you are running the same shot multiple times), which in turn contains the subdirectory eqdsk, which finally contains the equilibrium file (it’s a plain text file). In other words, the path to the equilibrium file will look approximately like this:

/compass/Shared/Users/jirakova/13929/run\_0/eqdsk/g\_p13929_t1.15000

The name of the file indicates the shot number and time in seconds at which the equilibrium was reconstructed. The EQDSK file may then be loaded with PLEQUE, see above.

If you forget which run is which, each contains a runcmd file with the command you used to launch efit.py.

To create a reconstruction which uses corrected IPR coils positions and flux loops and divertor Mirnov coils signals, simple leave out the --floff and --mirnoff options:

./efit.py 13929 1.150 -t -e -x

Reconstruction variant: using Thomson scattering profiles of $p_e$

The total pressure profile $p(\Psi)$, where $\Psi$ is the flux function, is one of the free parameters of the Grad-Shafranov equation, and therefore of all Grad-Shafranov solvers such as EFIT. Usually it is assumed that the total pressure profile is quadratic and falls to zero at the separatrix, which is inaccurate and introduces errors into the magnetic equilibrium reconstruction. Thanks to Michael Komm and Tomáš Markovič, the electron pressure profile measured by the Thomson scattering diagnostic can be used as a better approximation, theoretically yielding better results. The details are given in Michael's presentation.

To create your own reconstruction using the TS electron pressure profile:

Connect to Abacus and start an interactive job.
```
ssh -oHostKeyAlgorithms=+ssh-rsa -X username@abacus01   # tok.ipp.cas.cz
qsub -I
```
The -oHostKeyAlgorithms=+ssh-rsa may be left out depending on your SSH key setup. You may be prompted for your password.

Load the appropriate modules.

module purge
module load compass
module switch pycdb/0.5-dev
module unload idl
module load idl/8.2
module load hdf5
module load idam

Go to the folder where the script is stored.

cd /compass/Shared/Common/IT/projects/idam-sdas-utils

Run the script. It will take about 30 minutes to complete.
```
./create_efit_video.pl -p 2 -y TS_2pe_EFIT 16908
```
The argument -p 2 tells the script to use the total pressure profile $p = p_e . 2$ (for lack of knowledge about $p_i$), -y TS_2pe_EFIT names the calculated signal variant, and 16908 is the shot number.

The equilibrium reconstruction variant can then be loaded from CDB using PLEQUE, see above.

Reconstruction variant: using separatrix location calculated by the two-point model

The equilibrium reconstruction algorithm can be constrained not only by "magnetic field at [R,Z,$\Phi$] is 0.35 T", but also by the presumable separatrix position. It is all explained in Michael Komm and Miroslav Šos's presentation. Here's how it works in a nutshell.

The outer strike point electron temperature $T_{e,t}$ is gauged as the peak $T_e$ measured by the combined ("new") divertor probe array. These measurements are routinely available in diverted discharges.
The upstream separatrix $T_{e,u}$ is calculated from $T_{e,t}$ using the two-point model. The model is simple and, for the typical conditions of COMPASS, might as well be substituted by the assumption that $T_{e,u} = T_{e,t}$.
The upstream $T_e$ profile measured by the Thomson scattering diagnostic is reasonably smoothed and the location where $T_e = T_{e,u}$ is found. This is the separatrix location $Z_{sep}$.
The separatrix location [$R$=0.557 m, $Z = Z_{sep}$] is fed to EFIT as a constraint. This theoretically improves the reconstruction accuracy.

To create your own reconstruction using the two-point model separatrix location, kindly ask Michael Komm (komm@ipp.cas.cz) for help.

Reconstruction variant: using custom separatrix location

If you have reliable data where the separatrix is (for now, only at one spot inside the vacuum vessel), you can feed it as constraint to a modified EFIT reconstruction script written by Michael Komm (komm@ipp.cas.cz).

Connect to Abacus and start an interactive job.
```
ssh -oHostKeyAlgorithms=+ssh-rsa -X username@abacus01   # tok.ipp.cas.cz
qsub -I
```
The -oHostKeyAlgorithms=+ssh-rsa may be left out depending on your SSH key setup. You may be prompted for your password.

Load the appropriate modules.

module purge
module load compass
module switch pycdb/0.5-dev
module unload idl
module load idl/8.2
module load hdf5
module load idam

Go to the folder where the script is stored.

cd /compass/Shared/Common/IT/projects/idam-sdas-utils

Run the script. It will take about 30 minutes to complete.
```
./create_efit_video.pl -f 0.557-0.285-1100 -y "t1100_ZsepTS_285" -q 17588
```
Arguments: - -f 0.557-0.285-1100: Constrain the separatrix position at [$R=0.557$ m, $Z=0.285$ m] at $t=1100$ ms. - -y "t1100_ZsepTS_285": Name of the variant. The name shouldn't contain blank spaces, shouldn't begin with a number and should be shorter than 20 letters. - -q: Do not carry out post-processing, such as calculating $P_{rad}$. This is optional. - 17588: Shot number.

The equilibrium reconstruction variant can then be loaded from CDB using PLEQUE, see above.

Comparing reconstructions

If you're wondering which of these reconstruction variants to use, load PLEQUE and see the difference for yourself.

from pleque.io.readers import read_geqdsk
import pyplot as plt

eq1 = read_geqdsk('/path/to/eqdsk1')
eq2 = read_geqdsk('/path/to/eqdsk2')
eq1.first_wall.plot(c='k')
eq1.separatrix.plot(label='eq1')
eq2.separatrix.plot(label='eq2')
plt.gca().set_aspect('equal')
plt.legend()

This is an example of such a comparison for discharge #16908, time 1130 ms (inter-ELM H-mode). Featured are the standard CDB reconstruction ("CDB"), Ondra's reconstruction with flux loops and Mirnov coils off ("Kovanda") and Ondra's reconstruction with flux loops and Mirnov coils on ("Kovanda FLMC"). On the outer midplane ($Z=0$ m), the separatrix position is 71.6 cm, 72.1 cm and 72.9 cm respectively. That's a big difference when it comes to edge modelling. Which of these reconstructions should one use? Frankly, I don't know. It is discussed in Katka's PhD thesis study, chapters 1.4 and 3.3.

eq1

Converting the EQDSK file for DivGeo

DivGeo is the UI (User Interface) which allows users to prepare grids for SOLPS-ITER. As basic input, it requires an equilibrium file in the format .equ, which is different from the EFIT output format, EQDSK. Luckily, SOLPS-ITER knows how to make the conversion.

Copy the equilibrium file into `baserun/`

At this point, we switch from the Abacus command line to wherever your SOLPS-ITER is installed - say, Soroban.

ssh -X username@soroban-front-01

(The -X option is mandatory here, as it allows for opening windows outside the terminal, such as DivGeo.) Alternatively, you may set up the NX client X2Go (described in Remote access: X2Go) to open a virtual desktop.

A new equilibrium usually means a whole new SOLPS-ITER simulation. One can either start solving the plasma state from scratch (so-called "flat profiles", although they are far from flat) or import a previous solution into the new magnetic geometry using the b2yt command (see chapter 3.14 of the manual). By the conventions of SOLPS-ITER, a new case requires you to create a new folder for the simulation. We assume you've already created it in $SOLPSTOP/runs/ as described in Creating a new SOLPS-ITER run.

If your SOLPS-ITER installation is located on the same server as your EQDSK file:

cp /path/to/eqdsk/g_pNNNNN_tTTTTT /path/to/baserun/g_pNNNNN_tTTTTT

If your SOLPS-ITER installed, for instance, on Marconi Gateway, you can use SCP to transfer the files:

scp /path/to/eqdsk/g_pNNNNN_tTTTTT username@login.eufus.eu:/path/to/baserun/g_pNNNNN_tTTTTT

Converting and refining the equilibrium file

Finally, switch to the SOLPS-ITER terminal again. Assuming that the simulation directory is ready and the SOLPS-ITER environment is set up, you may now convert the file to the DivGeo .equ format:

e2d g\_pNNNNN\_tTTTTT g\_pNNNNN\_tTTTTT.equ

Check the equilibrium by opening DivGeo

dg &

and File > Import > Equilibrium. Ctrl+P to zoom in, File > Import > Topology > SN to display the separatrix. You might find that the equilibrium is too coarse, especially around the X-point.

Fortunately, there is an option to interpolate the grid to twice the grid density:

d2d g\_pNNNNN\_tTTTTT.equ g\_pNNNNN\_tTTTTT.x2.equ

Rinse and repeat this command about a few times (3-4x), interpolating the equilibrium to a higher density. This will yield something like the following lovely resolution:

Important note: It is recommended to increase the resolution no more than 32x, otherwise you can get the error "The grid you have specified seems to be too fine".

Experiment until you have a suitable .equ file to build grids upon. If you encounter problems, refer to the DivGeo tutorial or Creating a new SOLPS-ITER run.