Comparing experiment to model
Interpretative modelling of tokamak experiments comes with careful comparison of modelling and experimental results. Here such a comparison is showcased on the COMPASS discharge #17588. More comprehensive information is given in Katka's PhD thesis study, chapter 3 Interpretative modelling results.
Scouting a discharge
Before you commit to modelling a particular discharge, scout it in advance. This involves judging:
-
Is there a part of the discharge where the parameters are reasonably stable and which captures the physics you wish to study (L-mode, H-mode, impurity seeding, isotope mix...)?
-
Is the discharge covered with diagnostics which are pertinent to the physics you wish to study and to the overall simulation quality?
Creating a simulation takes a lot of time and effort, even if it's a simple pure D run. So scout carefully and pick the discharge wisely.
Time traces of basic plasma parameters
First, consult the COMPASS logbook entry. Typically, you will model a steady-state plasma with SOLPS-ITER, so make sure the real plasma is also steady-state. This means the following parameters must be stable:
- plasma current, loop voltage
- electron density: line-averaged as given by the interferometer, peak electron density as given by the Thomson scattering diagnostic, gas puff
- plasma shape (a video from the standard EFIT reconstruction is linked in the logbook entry), strike point position, outer midplane separatrix position, plasma top separatrix position
- \(H_\alpha\) line emission intensity (especially crucial for investigating inter-ELM H-mode)
- heating power: NBI power, Ohmic power, energy content in the plasma
- plasma composition: impurity seeding presence/magnitude
If you like, plot the time traces yourself and investigate them up close.
This will give you an idea whether the shot is low- or high-density, what kind of heating power it has and which time period you should like to model. In shot #17588, I decided to model time 1100 ms, just before NBI 2 comes online (magenta line).
To determine the power crossing the separatrix \(P_{SOL}\), on COMPASS one may use the signals P_sep_norm
(power crossing the separatrix calculated using a more sophisticated algorithm than \(P_{ohm}-P_{rad}\), normalised to the separatrix surface) and A_sep
(separatrix surface). The NBI contribution is more tricky as the absorbed power is smaller than the nominal power; ask Klára Bogár (klara.bogar@ipp.cas.cz) for a more nuanced estimate.
Find available diagnostics
As stated in the Table, sheet Diagnostics, the crucial diagnostics for SOLPS-ITER interpretative simulations are:
- the Thomson scattering diagnostic (electron temperature \(T_e\), electron density \(n_e\)) -
Te
,ne
and their stray-light-corrected variants - the "new" divertor probe array (plasma potential \(\Phi\), ion saturation current \(I_{sat}\) and electron temperature \(T_e\)) -
DIVBPP01
-DIVBPP55
,DIVLPA01
-DIVLPA54
andDIVLPB01
-DIVLPB54
- the IR camera (divertor heat flux \(q_\parallel\)) -
FIRCAM_heat_flux
- the AXUV diode (bolometer) array (line-averaged radiated power \(P_{rad}\)) -
Prad
and many signals startingAXUV
- the horizontal reciprocating probe (plasma potential \(\Phi\), ion saturation current \(I_{sat}\) and electron temperature \(T_e\)) -
rcp_position_horizontal
and an array of signals depending on the probe head, such asBPP1_floating
orLP2_Isat_current
Figure: Essential COMPASS diagnostics in the poloidal cross-section. Example SOLPS-ITER grids included.
To find whether these diagnostics measure, search the CDB (COMPASS DataBase) for the corresponding test signals:
- Thomson scattering diagnostic -
Te
- "new" divertor array -
DIVBPP01
- IR camera -
FIRCAM_RAW
- AXUV diode (bolometer) array -
Prad
- horizontal reciprocating probe -
rcp_position_horizontal
If these signals are present, the diagnostic was most likely collecting data. If they aren't present, the diagnostic was offline and you might want to switch to another discharge. The reciprocating probe is not entirely essential, but is is very helpful for gauging the equilibrium reconstruction accuracy at the outer midplane.
Check experimental data quality
Here I'll list a couple of criteria I know to judge the shape of the diagnostics in question. For more information, ask the responsible diagnosticians (listed in the Table, sheet Diagnostics).
Horizontal reciprocating probe
Located on the outer midplane, the horizontal reciprocating probe (HRCP) provides some of the best "upstream" measurements. Its measurements of the plasma potential (ball-pen probes) are unique. Its uncertainties/drawbacks/issues involve:
- The probe has a pronounced effect on the plasma. Sometimes you can see both on the Thomson scattering profile and on the divertor that as the probe is inserted, the whole edge plasma cools down. It is an open question whether this means the measurement of the probe tips are untrustworthy, as the majority of the cooling is probably done by the probe head and the tips are deeper in the plasma.
- HRCP often doesn't measure at all, so its measurements aren't routinely available.
- When HRCP does measure, it might not penetrate deep enough to hit the separatrix.
- Over the years, multiple probe heads were installed on the horizontal reciprocating manipulator. You must find out which probe head was installed in your shot and what the signals from individual pins were called.
- The probe position might be given inaccurately in the database because of an undetected systematic shift of the probe drive.
The primary quantities measured by the horizontal reciprocating probe are:
- plasma potential \(\Phi\) - signal
BPP1_floating
or something similar - ion saturation current \(I_{sat}\) -
LP2_Isat_current
or something similar - floating potential \(V_{fl}\) -
LP1_floating
or something similar.
I'll come back to this and make pictures of how the signals should look when I have the time.
The "new" divertor array
Located on the divertor, the "new" divertor array provides useful "downstream" measurements. One of its biggest selling points is measuring with high temporal resolution, which is sadly not pertinent to SOLPS-ITER simulations since those tend to describe the steady-state and the measurements are averaged over a time period anyway. What is fortunate is that it measures almost in every shot, so its measurements are routinely available (in D-shaped discharges, obviously). Its uncertainties/drawbacks/issues involve:
- The profile often shows "teeth", where one probe measures a systematically and strangely lower/higher signal than its neighbour. Be careful about them. Some say they're the effect of drift and drift-enabled SOLPS-ITER modelling would reproduce them, but they might also be bad probes.
- The divertor can be shadowed by one of the reciprocating probes. It isn't entirely clear whether it then gives bonkers values (the profile can exhibit strange spikes which are there only when the probe's deep enough) or correct measurements.
- All the HFS (inner target), probe measurements are untrusted because they yield negative \(T_e\). It is not known why, but it is known that reversing the fields (and thus the direction of the \(E \times B\) drift) alleviates this problems. Thus it is thought to be a magnetic shadowing problem.
The primary quantities measured by the "new" divertor array are:
- plasma potential \(\Phi\) - signals
DIVBPP01
toDIVBPP55
- ion saturation current \(I_{sat}\) -
DIVLPB01
toDIVLPB24
andDIVLPA25
toDIVLPA54
- floating potential \(V_{fl}\) -
DIVLPA01
toDIVLPA24
andDIVLPB25
toDIVLPB54
To check the data quality, plot the \(T_e = (V_{BPP} - V_{LP}/1.4)\) time trace for the whole array, compare it to the whole shot time traces and look for badly measuring probes. Decide where you'll want to average the divertor quantity values (red lines in the picture).
Figure: Reciprocating probes movement vs divertor \(T_e\) measurements in COMPASS discharge #19418. Only the horizontal reciprocating probe is on. Note how the divertor profiles change after it passes the velocity shear layer (\(\Phi_{max}\), denoted by dotted lines). Also note how the strike point moves around \(t=1100\) ms.
IR camera
Viewing a special divertor tile, the IR camera provides "downstream" measurements of the target heat flux \(q_\parallel\). It's our most trusted measurement of the target heat flux (and it works across the whole target to boot, unlike the "new" divertor array). Its uncertainties/drawbacks/issues involve:
- Calculating the electron temperature from the raw camera data is a long process which involves a simulation using the THEODOR suite. As a result, there might be raw data available but no processed data available. In that situation, if you're sure you really want to model this discharge and you need the IR camera data, ask Petr Vondráček (vondracek@ipp.cas.cz) nicely to process the data for you.
- Background target heat flux must be subtracted from the data. In the zeroth approximation, you can just estimate its value as constant.
The primary quantity measured by the IR camera is the perpendicular target heat flux \(q_\perp\) - FIRCAM_heat_flux_profile
. Note that you'll have to convert this to the parallel heat flux using the sine of the divertor tile tilt angle (3 degrees) and the magnetic line impact angle (depends on the equilibrium and may be calculated using PLEQUE).
The Thomson scattering diagnostic
Measuring along a vertical chord on the plasma top, the Thomson scattering diagnostic provides invaluable "upstream" measurements. Unlike the horizontal reciprocating probe, it doesn't disturb the plasma and it is more or less a routine measurement. Its uncertainties/drawbacks/issues involve:
- The data needs to be corrected for stray light (light which wasn't scattered from the plasma but is reflecting off of walls). The corrected signal is stored under the name
Te/THOMSON:shot_number:stray_corrected
(and similar for the density and the errorbars). The correction is especially important in the SOL, so always strive to use the corrected data. If they aren't available and you are certain you want to model this discharge and you need the Thomson data, ask Miroslav Šos (sos@ipp.cas.cz) nicely. - The edge plasma can be cooled by one of the reciprocating probes. Check for this and don't model the times when the TS signals drop because of the probe.
- It is recommended to disregard the datapoints which have more than a 50% errorbar. Sadly, the physical meaning of the errorbars is not known (95% confidence interval? standard deviation?), though they do carry qualitative information.
The primary quantities measured by the Thomson scattering diagnostic are:
- the electron temperature \(T_e\) -
Te
orTe/THOMSON:shot_number:stray_corrected
- the electron temperature errorbars -
Te_err
orTe_err/THOMSON:shot_number:stray_corrected
- the electron density \(n_e\) -
ne
orne/THOMSON:shot_number:stray_corrected
- the electron density errorbars -
ne_err
orne_err/THOMSON:shot_number:stray_corrected
To check the data quality, plot the \(T_e\) and \(n_e\) time trace together with the separatrix position at the plasma top and compare it to the whole shots time traces. Decide where you'll want to collect the profiles (red lines in the picture).
The AXUV diode (bolometer) array
Consult Martin Imríšek (imrisek@ipp.cas.cz). The total radiated power is available in the Prad
signal, while the individual diode signals are named AXUV_A_1
to AXUV_F_20
.
Pick time to model
In steady-state plasma interpretative modelling, it is good practice to model one particular time instance of the discharge. This will inform the equilibrium reconstruction you will use as grid basis, time interval of collecting and averaging diagnostic data, the investigated physics (NBI heating, gas puff magnitude), the edge plasma parameters (density in a discharge with density ramp) and other aspects of the modelling. When choosing the modelled time instance, consider the following criteria.
- The plasma is steady-state. If possible, don't model a time instance where NBI has just been turned on, the plasma is going through MHD instabilities, RMP is being used or an ELM is underway. Pick a quiescent discharge phase (flat top), a long inter-ELM period or a spot on the density ramp where the plasma parameters evolve slowly and calmly.
- Diagnostic measurements are available. For instance, an inter-ELM time instance may be chosen where the Thomson scattering diagnostic takes a measurement.
- The plasma is not disturbed by diagnostic measurements. This applies particularly to reciprocating probes, which can cool the edge plasma if inserted close to the separatrix. Keep track of the probe movement, plot the upstream and target temperature evolution in time and model time instances when the probe is still far from the separatrix. (Ironically, if you're going to use the reciprocating probe data for experiment-model comparison, you will take probe data from the disturbed part of the discharge. This is an unavoidable, intrinsic feature of the probe measurements. If you are certain the probe cools the edge plasma substantially, don't use this experimental data. In the worst case, pick a different discharge to model.)
When averaging diagnostic data for experiment-model comparison, keep in mind that every diagnostic is different. Probe data are collected with high time resolution, which you don't need because your model averages over turbulence. Consequently, you can average them over a duration of almost any choice, as long as you process the errorbars in a sensible manner. The Thomson scattering diagnostic might make a measurement directly at the time you model, but even though it has usable errorbars, you will be better informed of its uncertainty if you use several consecutive measurements during a steady-state part of the discharge. A reciprocating probe may be measuring in a slightly different plasma than the time instance you're modelling and the plasma state may change as it travels; take this into account when gauging its measurement uncertainty. Don't average over ELMs when you model inter-ELM plasma. Discard outliers only after you have a sound justification for it. Consider spatial uncertainties as well as uncertainties in the measured values; the reciprocating probe position may be badly calibrated, the bolometer lines of sight might be slightly misplaced, the trustworthiness of anything (including, e.g. tomography reconstruction) is greatly reduced if it uses a flawed equilibrium reconstruction as input.
Separatrix position: the ever-present problem of code-experiment matching
As discussed at length in Katka's PhD thesis study, Katka's Master's thesis and [Jirakova 2018], the magnetic equilibrium reconstructions performed by EFIT++ at COMPASS are inaccurate by 1-2 cm when it comes to the separatrix position. Unfortunately, accurate separatrix position is absolutely crucial for experiment-code data comparison. If EFIT misplaces the upstream separatrix by 2 cm, you won't match the Thomson scattering profiles with realistic boundary conditions and diffusion coefficients for the love of God. The only solution - supported theoretically by the notion that you know the EFIT separatrix, and therefore the SOLPS-ITER grid built on top of it, is misplaced - is to shift the profiles radially by an ad hoc amount until they match nicely.
"Correcting the separatrix position is witchcraft." - Katka Hromasova
On the one hand, make no mistake: this magical solution supported by hand-waving is an embarrassment for anyone attempting interpretative simulations. On the other hand, there is often no other way. My suggestion is the following. Do what you can, make your best guesses and then disclose everything including the uncertainties. Being upfront about the spells you employed is better than mystifying your reviewers and colleagues with your clairvoyance. Possible ways to get a better hold of the separatrix position include:
- Scout several equilibrium variants for the time you're modelling. The standard CDB variant, the
V4_std_O
variant (described in [Kovanda 2019]), a reconstruction using experimental plasma pressure profile or your own reconstructions (described in Creating an equilibrium file), whatever you can get your hands on. These will inform you of the uncertainties in the separatrix position. If you're lucky, you'll even come across a reconstruction which requires little to no separatrix position correction. Then, supposing you can justify why this reconstruction is better than the others (e.g. it involves more constraining input, its input has higher quality, it matches divertor probe data...), you might be able to avoid the separatrix uncertainty issue altogether. Cool! - Check time traces of the strike point positions
R_strike_points
, OMP separatrix positionR_mid_out
and plasma top separatrix positionZ_TOP_TS
. Their temporal variation will inform you of the reconstruction convergence and uncertainty. Do this for all the EFIT variants you can get your hands on; some converge worse than others. - If reciprocating probe data (horizontal or vertical) is available and the probe is inserted deeply enough, check the velocity shear layer (VSL) position. It can most easily be detected as the maximum in the ball-pen probe potential, plus minus a few mm due to the \(T_e\) contribution. If no ball-pen probe is available, try cross-correlating two poloidally spaced \(I_{sat}\) Langmuir probes to find the poloidal plasma velocity \(v_p\) and locate the point where it changes sign. It is not known whether the VSL forms at the separatrix (though Jiří Adámek argues it does), but I'm pretty sure it forms systematically near it, probably a bit outside it. [Jirakova 2018]
- Use common physics sense. Separatrix electron temperature can be 200 eV in ITER, but the same result in COMPASS it means EFIT places the separatrix way too inward. If the strike point position is misplaced in regard to the target \(T_e\) peak by 1 cm even though the shot is firmly attached, the reconstruction isn't probably that good. If the equilibrium reconstruction says the plasma-wall clearance is really small between the OMP and the plasma top, and you do see a big reduction in far SOL temperature on the plasma top and inner target, then that reconstruction might be on to something.
Figure: Comparison of four equilibrium variants in H-mode discharge #16908 (two reconstructed for the whole discharge, two only for the modelled time) to the strike point position recorded by the "new" divertor array and the IR camera (bottom right). Notice the inner strike point disagreement. Is this caused by the probes (which sometimes measure strange values on the HFS, probably due to currents flowing into the grounded divertor tiles), or by the equilibrium reconstruction inaccuracy? Is the "Kovanda FLMC on" variant, which includes input from the divertor Mirnov coils and therefore could arguably yield more accurate results at the divertor, better than the others? It is known that the Mirnov coil data is noisy and causes bad convergence of the EFIT reconstruction. What is the best choice?
Figure: The four equilibrium reconstruction variants at \(t=1130\) ms. Notice the spread of OMP separatrix position, and how every variant other than the standard CDB variant (blue) places the separatrix more outward.
- Learn EFIT++ and work on fixing its systematic errors. This is the hardest way, but if you make it, you will be hailed a God by the edge community. You can start where Ondřej Kovanda left off - ask around about his work and search TOSEP for his presentations.
- Lastly, and this is a work in progress, I believe you can backtrack and use a simple SOLPS-ITER simulation to validate the reconstruction. As argued in Katka's PhD thesis study, if you build a SOLPS-ITER simulation atop a particular equilibrium reconstruction and need no radial shifting or extravagant boundary conditions and diffusion coefficient profiles to match the experimental data, this is solid evidence that the reconstruction is accurate. Furthermore, if you do have to result to radially shifting the profiles but can achieve good experiment-code match without extravagant boundary conditions and diffusion coefficient profiles, the employed shift indicates by how much EFIT misplaced the separatrix. SOLPS-ITER is a physically very sound code. It's hard to use, but it's one of the best at gauging the overall edge plasma state. Making a SOLPS-ITER simulation for every equilibrium variant is lengthy, but if it can be automatised, I believe it can become a convincing argument why certain reconstructions are better than others.
Comparing code output and experimental data
If comparing simulation and experimental data requires data post-processing, it is generally better to perform the transformations on the modelled data and process the experimental data as little as possible. For instance, it is better to integrate plasma radiation along a line of sight than to perform tomographic reconstructions from bolometer signals.
In this figure, we compare the upstream \(T_e\) (Thomson scattering), upstream \(n_e\) (Thomson scattering), upstream and target \(p_e=e T_e n_e\) (Thomson scattering, divertor probes), upstream and target \(T_e\) (Thomson scattering, divertor probes), upstream and target \(q_\parallel\) (divertor probes, IR camera) and the target \(I_{sat}\) (divertor probes).
Upstream \(T_e\) and \(n_e\)
Measured by the Thomson scattering diagnostic. The raw input of the diagnostic is scattered light intensity in several wavelength channels, but the data processing is automatised and well-documented to the point where comparing directly \(T_e\) and \(n_e\) is the most viable option.
Outer target ion saturation current \(I_{sat}\)
Measured by the "new" divertor probe array.
\(I_{sat} = feA_{probe}c_{s}n_{e}\)
where \(f = 0.5\) is a factor due to density fall inside the sheath, \(e = 1.6 \times 10^{- 19}C\) is the elementary charge, \(A_{\text{probe}}\) is the effective probe collecting area, \(c_{s}\) is the sound speed at the sheath entrance and \(n_{e}\) is the density at the sheath entrance.
The factor \(f = 0.5\) effectively says that because ions are accelerated to the sound speed within the sheath, the density of the plasma at the probe surface is only half of the density at the sheath entrance.
The effective collecting area used in the roof-shaped divertor probes is their parallel geometric projection along the field lines, \(A_{probe} = 2.8\) mm2. This means that the Larmor radius is considered infinitely small (thus particles cannot arrive “from the side” but only parallel to the magnetic field) and sheath expansion is considered non-existent.
The sound speed in its most basic definition is
\(c_s = \sqrt{\frac{\gamma_e e T_e+\gamma_i Z e T_i}{m_e+m_i}}\)
where \(\gamma\) is the polytropic coefficient for each species, \(m_{i} = 3.34 \times 10^{- 27}\) kg is the deuteron mass and the temperatures are given in electronvolts. Since \(m_{e} \ll m_{i}\), the electron mass in the denominator can be safely neglected. The numerator is more complicated, as \(\gamma\) is a mystery and \(T_{i}\) is typically not even measured. According to David Tskhakaya, \(\gamma_{e} = \gamma_{i} = 1\) (isothermal process) is an appropriate assumption. As to the temperature, experiments and modelling suggest that \(T_{e} < T_{i}\) in the SOL as electrons cool down faster, but this can change with a number of factors: collisionality along the flux tube, relation between ion and electron temperature at upstream, heating type… If no better data is available, one often assumes \(T_{e} = T_{i}\) (thermal equilibrium). Employing all of the above assumptions, we write \(c_{s} = \sqrt{2eT_e/m_i}\).
The said relations are derived from the collisionless sheath theory.
Parallel heat flux \(q_\parallel\)
In B2.5 the electron and ion heat flux are composed of the following parts:
- conductive heat flux, proportional to the heat conductivity and temperature gradient of the respective species
- convective heat flux, proportional to the species temperature and particle flux
- kinetic heat flux, proportional to the species mass and squared velocity (neglected for electrons)
However, at the target the temperature gradient typically goes to zero (I'm not sure if that happens naturally, is assumed within the isothermal sheath theory or is enforced by the boundary condition, but there you go), and so only the latter two parts remain. According to Eq. (14) in [Kotov & Reiter, 2009], the heat flux is
\(q_\parallel = {\gamma_{e}T_{e}\Gamma_{e} + \frac{5}{2} \sum_\alpha n_{\alpha}u_{\alpha}\ + \frac{1}{2} \sum_\alpha m_{\alpha}n_{\alpha}u_{\alpha}u_{\alpha}^{\parallel 2}}_{}T_{e}\)
where \(\gamma_{e}\) is the electron sheath heat transmission coefficient. Comparing this to actual B2.5 boundary conditions yields… maybe a somewhat qualitative agreement? Here’s me delving into the electron heat flux boundary condition BCENE=3
:
You can see that there are two terms, out of which the first could be the convective heat flux (it even has the BCENE
boundary condition value \(p_{h_{e},1}\) in it). The second one is maybe zero, maybe non-zero, and \(p_{h_{e},1}\) definitely is not equal to \(\gamma_{e}\) from the previous equation. The situation is even more complicated with ions, where the prescription for BCENI=3
is:
Why is the kinetic energy flux (the only term with masses and squared velocities) zero unless drifts are on? Why are there two recommended values for \(p_{h_{i},1}\)? (Actually I might know that one. In B2.5 equations the ion heat flux is defined with a factor \(\frac{3}{2}\) before the convective part, while other sources, such as the abovementioned Kotov & Reiter, demand that it is \(\frac{5}{2}\). So I guess this parameter isn’t as free as you’d think.) Who is Boris? (I found this out years later by accident. BoRiS is a 3D transport code for modelling stellarators.) Why is there electron temperature in that term and what does that term even mean?
In the end, as Kotov and Reiter suggest, it’s easier to just sum the incident heat fluxes and divide them by \(\Gamma_{e}T_{e}\). That way you’ll get the total sheath heat transmission coefficient, which they define as
\(\gamma = \gamma_{e} + \frac{5}{2}\frac{t}{z} + \frac{1}{2}(1 + \frac{t}{z})M^{2}\)
where \(t = \frac{T_{i}}{T_{e}}\), \(z\) is the effective ion charge, \(M\) is the Mach number and it is assumed all ion species have the same velocity at the target. On the other hand, a similar definition of \(\gamma\) based on B2.5 equations would be a terrible mess.
On the experiment side, the situation is delightfully simple. The total target heat flux is
\(q_\parallel = \gamma fec_{s}n_{e}T_{e}\)
where this time \(\gamma\) is the sheath heat transmission coefficient. Otherwise the same definitions and considerations as in the ion saturation current apply.
The sheath heat transmission coefficient is usually assumed to be \(\gamma = 7\); however, to match the outer divertor array and the IR camera, \(\gamma = 11\) must be used. This value can be retrieved when one uses non-ambipolar heat flux to the target, which is the case when the target tiles are grounded. We use the latter value.
Formulas for the divertor probes measurements
The directly measured quantities are the ball-pen and Langmuir probe floating potentials, \(V_{\text{BPP}}\) and \(V_{\text{LP}}\) respectively, and the ion saturation current \(I_{\text{sat}}\). From them, all the rest can be calculated as:
Electron temperature: \(T_{e} = \frac{V_{\text{BPP}} - V_{\text{LP}}}{1.4}\)
The coefficient is 1.4 for the new divertor array; it is different for the reciprocating probes. By subtracting the raw potential data, \(T_{e}
\) can be calculated including its fluctuations.
Electron density: \(n_{e} = I_{\text{sat}}/feA_{\text{probe}} \sqrt{\frac{2eT_e}{m_i}}\)
Since the expression diverges for \(T_{e} = 0\), the electron density cannot be reliably calculated including its fluctuations. The second best is to limit the raw \(T_{e}\) data by some threshold, for instance 0.5 eV.
Parallel heat flux: \(q_\parallel = \gamma\frac{I_{\text{sat}}T_{e}}{A_{\text{probe}}}\)
The heat flux can be calculated including its fluctuations.
Bolometry
b2plot
contains routines for integrating plasma radiation along a given line of sight.
Relevant manual sections:
- I.5.3-I.5.8: general description of chord files and examples of chord files shipped with SOLPS-ITER (great for reverse engineering)
- I.7.3: example commands to visualise your chords (good for reverse engineering) and plotting the radiation intensity collected among many chords (the x axis in figure I.9(b) is chord number)
Bolometers measure the sum of radiation intensity over a broad spectral range (see the COMPASS wiki page on bolometry) integrated along the bolometer line of sight (so-called chord). COMPASS has six AXUV bolometer arrays: A, B, C, D, E and F, each consisting of 20 bolometers. Their position and chord numbering are on the wiki. Lines of sight of each bolometer are recorded in one of the Excel sheets on the wiki; they are given as the detector \([R, Z]\) position and its poloidal viewing angle in degrees (viewing outward is 0, viewing upward is 90 etc.).
The b2plot input files, however, require a different format of chords: the Cartesian coordinates of the start and end of each chord. To convert \([R_{detector}, Z_{detector}\), viewing angle] to \([R_{start}, Z_{start}, R_{end}, Z_{end}]\), one may use the following Python code:
angle = np.radians(angle)
chord_length = 1 #m
R_end = R_start + chord_length*np.cos(angle)
Z_end = Z_start + chord_length*np.sin(angle)
b2plot
in the format specified in the manual:
‘Name of the array’ number_of_points_sampled_along_the_chord
X1_start Y1_start Z1_start X2_end Y2_end Z2_end chord_number
...
All COMPASS bolometers look straight toward the centre of the tokamak (toward the \(Z\) axis). So if you’re using a single line of sight for each bolometer, using \(Y=0\) for all chords suffices. However, when you wish to take into account that the detector actually collects radiation from a cone by averaging over a number of chords spanning its entire view cone, you will have to dip into the \(Y\) direction as well.
Note that most of the bolometer signal comes from the core plasma. This means that if a chord crosses the core, you can’t use this chord for code-experiment comparison. Instead one must use the chords which travel exclusively through the area simulated by SOLPS-ITER. To check which chords fulfil this requirement, use the example given in the I.7.3 section of the manual. Based on my analysis, the bolometer array most relevant to SOLPS modelling is the C array since it overlooks the divertor area, where most of the edge radiation can be expected to occur, with a grand total of 11 chords.
Diffusion coefficient profiles
Ahh, anomalous diffusion coefficients. Physically nonsensical and lacking in meaning, but so beloved and so often studied. Setting the particle diffusion coefficient \(D_{n,a}\) for every ion species \(a\) including impurities and neutrals, electron heat diffusivity \(\chi_e\) and ion heat diffusivity \(\chi_i\) will be necessary to achieve your experiment-code match. Information on how to set them and lamentation on how to interpret them can be found in Katka's PhD thesis study, in Questions and answers: How do I give the diffusion coefficient a radial profile?, and in Processing SOLPS-ITER output: Checking the perpendicular transport coefficients. Tips for setting them:
- Don't overfit. Yes, it's very nice when your SOLPS-ITER results follow the experimental profile point by point, but what is the actual physical value of such a result if you set a different diffusion coefficient for every flux tube and their values are all over the place? 10 \% deviations are considered a great match in edge plasma modelling, a wonder of modern science, and I don't think it's worth the time to pursue a closer match.
- According to David Coster, simulations of L-mode can often do with a flat profile of \(\chi\) and \(D_n\). H-mode will probably need a profile of one or both to match the pedestal shape.
- Theoretically, much of the perpendicular transport is carried by turbulence, which carries particles and energy inside a single structure. Consequently, the values of \(D_n\) and \(\chi\) should be close to one another. My SOLPS-ITER modelling has not supported this so far, with heat diffusivity being up to 10x higher than particle diffusivity, but it's something to keep in mind. A mismatch in \(\chi\) and \(D_n\) by two orders of magnitude should probably make you wonder if something is wrong.
- \(D_n\) can be realistically set for D\(^{1+}\), as this makes the largest contribution to the electron density profile. \(\chi_e\) can be sort of realistically set as well thanks to the \(T_e\) profile measurements, as long as you're confident about what fraction of \(P_{SOL}\) goes to the electrons and to the ions. However, \(D_n\) of other ion species (neutrals, impurities) and \(\chi_i\) are completely up in the air (to my knowledge). Lacking any foothold, do what every edge physicist does: assume equality. Set equal \(D_n\) for all ion species and \(\chi_e=\chi_i\). Yes, this may not reflect reality at all, it might overestimate carbon transport to the core, it might distort \(T_e/T_i\) in the SOL. But it's probably the best guess you have, and pretending otherwise will get you accused of overfitting (by me).