Feature blog aims to introduce additional important features of SOLPS-ITER which did not fit in the main step-by-step tutorials. Familiarity with the basic workflow is assumed.

Gas puffing and pumping

Motivation

Adding explicit pumps (particle sinks) and fuelling puff valves (particle sources) is required for proper simulation of the plasma and neutral density. While there may be some implicit sinks and sources provided by some of the fluid boundary conditions that help to keep your particle balance, it is often not enough for a realistic simulation.

Note

Puffing and pumping is best realized in a coupled B2.5 + Eirene run, since it mainly involves the neutral particles. Unless stated otherwise, this tutorial focuses on the Eirene-based features, which should be preferred whenever possible. The alternatives for B2.5 standalone runs are discussed more briefly in a note at the end of each section.

Particle sources and sinks

The main particle sources and sinks in the tokamak edge plasma are:

Recycling. Plasma ions impinge on the divertor target (or limiter, in the limiter configuration), grab a free electron and become neutral particles. For an indefinite period, they chill at the target and enjoy chemical reactions, until they randomly desorb due to thermal motion or they are struck free by another impinging particle. The desorbed neutral can be an atom (D⁰) or a molecule (D₂, CH₄). The neutral collides with a plasma electron, the molecule is broken up and the atoms are ionised. Depending on the plasma temperature and density, the ionisation can take place just above the target, in the divertor volume, or even in the confined plasma (especially with an open divertor). This entire process from a plasma ion to a freshly ionised ion is called recycling. Recycling may and may not preserve the total number of ions. One defines the recycling coefficient $R$ as the ratio between the ion flux leaving the solid surface and impinging on the solid surface.

\[R = \frac{\Gamma_{\text{neutrals from target}}}{\Gamma_{\text{ions to target}}} \approx \text{0.8-1}\]

If the recycling coefficient is unity ($R=1$), every neutralised particle is released from the target again. This is what happens in steady state. There is a constant inventory of neutrals sticking to the solid surface, and the stream of incoming ions and outgoing neutrals balance each other out. If the recycling coefficient is below unity, usually $R=0.8-0.99$, some percentage of neutralised particles are "absorbed" into the target. This usually happens on the first wall of tokamaks with a short pulse. The wall was pumped clean of particles prior to the discharge, and the plasma doesn't last long enough to saturate the first wall with neutrals. Hence, the first wall acts similarly to a vacuum pump during the discharge. More on this below.
Gas puff. Besides the initial neutral gas pressure, plasma density in a tokamak is controlled by a real-time neutral gas injection, called gas puff. Gas puff is typically implemented as a small tube ending inside a tokamak port, with a piezovalve at the end. The piezovalve reacts to electric field by opening and letting gas flow out. The electric field is created with the piezovalve voltage (MARTE_NODE.FeedbackDataCollection.GasPuffAnalogueOutput in the experimental COMPASS database). The number of atoms/molecules released per second $\Gamma_{puff}$ [s^-1] at a certain valve voltage can be estimated calibrating the piezovalve.
Vacuum pumps. Vacuum pumps are, put simply, surfaces with a recycling coefficient below unity, $R<1$. (They are usually designed to pump neutrals, not ions, but the logic stands. If a neutral particle impinges on the pump surface, it has the chance $1-R$ that it will not desorb again.) The number of particles (atoms or molecules, depending on the pumped gas) absorbed per second $\Gamma_{pump}$ [s^-1] depends on the recycling coefficient $R$, the pump surface $A_{pump}$, its surface temperature $T$ [K] and the particle mass $m$:

$\Gamma_{pump} = 3.638 A_{pump} (1-R) \sqrt{k_BT/m}$

Sourced from the EIRENE manual, section 2.6.2 Effective pumping speed. Pumping surfaces don't necessarily have to be vacuum pumps; they can also be an entrance to neutral beam injectors or other adjoint vacuum spaces which are not modelled by EIRENE and where neutrals may be lost or wander forever.
Neutral beam injection, pellet injection etc.
Plasma-neutral interaction. Though recycling usually overshadows every other process of ion-neutral conversion, volumetric plasma-neutral interaction is important for the conversion of energy. The dominant process is usually charge exchange, with volumetric recombination playing a smaller role.

Considering all of these particle sources and sinks, one should arrive at a particle balance in steady state:

\[\Gamma_{sources} \, [\text{s}^{-1}] = \Gamma_{sinks} \, [\text{s}^{-1}]\]

SOLPS-ITER is inherently a steady-state code, so achieving particle balance is the main hallmark of a converged simulation:

\[\frac{\Gamma_{sources} - \Gamma_{sinks}}{\Gamma_{sources}} \lesssim 0.01\]

It should be noted, however, that experimental plasma is not necessarily in the steady state regarding particle balance. The primary actor is the tokamak chamber, which can fuel the plasma both with main ions and impurities.

The status of the tokamak chamber is a major uncertainty in tokamak operation. It depends on the time of the day, as the wall status evolves during the experimental campaign, and on chamber conditioning. In the COMPASS tokamak, the helium glow discharge, performed before every discharge, is a major source of helium impurities. Chamber baking is performed... in the morning? Boronisation, coating the plasma-facing components with a thin layer of boron oxide, mitigates especially oxygen impurities, and wears off over a few weeks or months.

Besides affecting the impurity content, the main chamber can be a major particle sink of the main ion species. In a short plasma discharge, a wall swept clean of neutrals with the glow discharge can be more important in terms of particles removed per second than even a cryopump, owing to its large surface and proximity to the plasma. This is the case of the COMPASS tokamak. In a long discharge, the main chamber will become saturated with neutrals and its recycling coefficient will approach unity, so that real vacuum pumps can impact the particle balance. This is the case of the ITER tokamak, and possibly of COMPASS Upgrade.

construction TODO: How to calculate the time scales of particle balance and assess the relative importance of the first wall pumping?

When setting up gas puffing and pumping in a SOLPS-ITER simulation, consider these time scales. If the discharge is short, chances are that you don't need to implement the vacuum pumps at all, and should rather choose $R<1$ for the main chamber.

Pumping

Instructions

To set up a case with Eirene pumping surfaces, you'll need a simulation (ideally converged) without any pumps.

The Goal

We want to configure one or more surfaces, either existing or new ones, with non-zero absorption near the location of the actual pump or its inlet. In case the wall pumping is the dominating effect, select the whole wall instead. Consider if the targets should have different absorption.

Terminology: DivGeo recognizes elements, which are used to define corresponding Eirene surfaces (namely the "additional surfaces").

Inside the SOLPS work environment, create a new case directory and its baserun and run subdirectories.
```
cd $SOLPSTOP/runs
mkdir case_with_pumping
cd case_with_pumping
mkdir baserun
mkdir run
```

Copy the .dg file from your pump-less simulation over to the new baserun.

cd baserun
cp ../../case_without_pumping/baserun/compass.dg compass.dg

Open DivGeo.
```
setenv DEVICE compass
dg compass.dg &
```
In DivGeo, open Variables>General Surface Data. Mark every wall element which you would like to designate as pump. If you prefer, you may also define your pump as a new free-floating surface, see more about those in the "Gas puff" section. More advanced configuration of the General Surface Data may be needed, refer to the corresponding surface properties in the EIRENE manual.
Set Absorption to something small, like 0.01. This corresponds to the recycling coefficient $R=0.99$, which should not impact the simulation much. Don't forget to hit Set.
Commands>Check variables, Commands>Rebuild Carre objects, File>Save, File>Output. Close DivGeo and re-run the standard carre - & triang workflow.
- only the Eirene input file should get modified

If the simulation runs smoothly, congrats, you have added a pump into the simulation.

Changing the recycling coefficient

Assuming a run with coupled B2.5 and EIRENE, you can change the recycling coefficient in the input file. Find the block 6a. General data for reflection model and the corresponding surface model (SURFMOD) with PUMP in its name. There should be at least one for each group of surfaces with different absorption value configured in DivGeo, e.g.:

Formatted (input.dat)JSON (eirene.input.json)

SURFMOD_4_CARBON_SPT_PUMP_580K
    1     2     0     0     1
1.20600E+03-5.00000E-02 0.00000E+00 0.00000E+00 0.00000E+00-1.00000E+00
1.00000E+00 9.00000E-01 1.00000E+00 1.00000E+00 5.00000E-01 1.00000E+00
1.00000E+00 1.00000E+00 0.00000E+00 0.00000E+00 1.00000E+00

{
  // ...
  "REFLECTION_MODELS": {
    "MANUAL": "http://www.eirene.de/eirene.pdf#subsection.2.6.1",
    // ...
    "SURFMODS": [
      {
        "NAME": "SURFMOD_4_CARBON_SPT_PUMP_580K",
        // ...
        "RECYCT": 9.0E-01,
        // ...
      }
    ]
  }
}

Here, the recycling coefficient is 9.00000E-01 (2nd number in the row for Formatted input file) and specifies recycling value $R=0.9$. In the EIRENE manual, the recycling coefficient is named RECYCT (chapter 2.6.1).

B2.5 Standalone Case

(Untested) Authors have not attempted to implement pumping in B2.5 standalone runs. The key might be either in configuring boundary conditions for continuity equation (BCCON) in the b2.boundary.parameters file or in configuring recycling for fluid neutrals in b2.neutrals.parameters. Please refer to other sources for more information.

Main ion species gas puff

If one has implemented a particle pump, the lost particles must be compensated for in the particle balance. Under default boundary conditions, where the plasma density is controlled by the electron density at the core boundary $n_{e,core}$ (BCCON = 1), this refuelling is done mainly through the particle flux from the core. The flux automatically adjusted to such a value that matches the density required in CONPAR. However, a more advanced way to control the plasma density is through the gas puff intensity. That's the topic of this tutorial.

Before you implement a gas puff into your simulation, consider if you need it. In low density discharges of the COMPASS tokamak, gas puff was actually off during the flat top. The piezovalve was open before the discharge and during the current ramp-up, but once steady state was established, the density feedback system found that no gas puffing was needed. The plasma density was sustained purely on the particles present in the chamber. The question how, or whether, to implement gas puff in such a case has not been resolved yet. In the following, we will assume a moderate to high density discharge, where the gas puff was on.

The COMPASS tokamak had three piezovalves for gas injection: one on the LFS just under the outer midplane, which was the most active one in deuterium density control, one on the HFS which was used for deuterium density control only rarely, and one at the divertor target (near LFS strikepoint), which was used for nitrogen, neon and other impurity puffing. If you are making an interpretative simulation of the COMPASS tokamak, chances are, the LFS gas puff was active. Check this with your favourite vacuum technician and the operator who performed the discharge, as the COMPASS experimental database does not hold records of the gas puff location. For calibration of the gas puff (conversion from piezovalve voltage to puffed particles per second), contact Jordan Cavalier. The puffed particle fluxes were of the order $\Gamma_{puff} \sim 10^{20}-10^{21}$ D₂ molecules/s.

Instructions

To set up main ion species gas puff in a SOLPS-ITER simulation, you will require a puff-less simulation, ideally converged and with a semi-realistic pump already implemented.

Goal

For a single gas puff, we want to create or select one surface and create a new gas puffing Eirene stratum (= particle source). The surface is only used as a location reference for the stratum. It's possible to use one of the existing wall surfaces, but not the target surfaces. Otherwise, a new free-floating surface must be created and configured as invisible to the particles.

Terminology: DivGeo recognizes elements, which are used to define corresponding Eirene surfaces (namely the "additional surfaces").

Inside the SOLPS work environment, create a new case directory and its baserun and run subdirectories.
```
cd $SOLPSTOP/runs
mkdir case_with_gas_puff
cd case_with_gas_puff
mkdir baserun
mkdir run
```
Copy the .dg file from your gas-puff-less simulation over to the new baserun.
```
cd baserun
cp ../../case_without_gas_puff/baserun/compass.dg compass.dg
```
Open DivGeo.
```
setenv DEVICE compass
dg compass.dg &
```
In DivGeo, create or choose one surface at or near the pump's location.
- The easy way is to choose one of the existing wall surfaces. But it must not be one of the target surfaces (i.e. in contact with B2 mesh), otherwise an obscure error will occur later when starting the simulation.
- The honest way is to create a new free-floating surface just inside the vessel wall (Edit>Create>Point, then Add element). Its normal vector must point outside, like the wall surfaces. Finally, make it invisible to the particles by specifying Variables>General Surface Data>Surface Type = 0 and ...>Surface Side = 0. Don't forget to hit Set.
In DivGeo, select Variables>Add>Gas puff. Ctrl+U to unselect all elements, then mark the element where you wish the gas puff to be located and click Puffing slot > Set.
Make sure that the Gas species shows the species which you want. Usually D2 for a deuterium fueling gas puff in a coupled simulation with D plasma species.
Set Puffed flux to something small, like 1.0e19, and hit Set. This is a fairly small gas puff which should not impact your simulation significantly.
Commands>Check variables, Commands>Rebuild Carre objects, File>Save, File>Output. Close DivGeo and re-run the standard carre - & triang workflow.
- Eirene input file should get modified
- b2.neutrals.parameters should get modified

If the simulation runs smoothly, congrats, you have added a gas puff into the simulation. This, however, is not the end of the story. To truly control the plasma density with the gas puff intensity, you must set the core boundary condition to zero flux and find the ideal gas puff particle flux. To do this:

Open b2.boundary.parameters and set the core (S) continuity equation boundary condition (BCCON) to 8, meaning "prescribe the total particle flux over this boundary". Set the corresponding CONPAR (where you were likely to keep the desired $n_{e,core}$ until now) to 0.0. This will turn off any deuterium ion flux from the core.
To change the gas puff intensity, open b2.neutrals.parameters. One of the first lines contains the list of strata types:
b2.neutrals.parameters
```
...
crcstra= 'W', 'E', 'S', 'S', 'N', 'C', 'V', 'T',
...
```
This is the list of stratas (see Questions and answers: What is a stratum?, short answer: particle source in Eirene). The stratum C corresponds to your gas puff. A bit lower in the file, you will find the userfluxparm variable:
b2.neutrals.parameters
```
...
userfluxparm(1,1)=  0.00    ,  0.00    ,  0.00    ,  0.00    ,  0.00    , 1.000E+19,  0.00    ,
...
```
The nonzero value 1.000E+19 near the end, corresponding to the stratum C in the list of strata types, is the amount of particles puffed per second. If you have multiple gas puffs, you may have several C strata. Edit the particle flux value to something semi-realistic. In the case of COMPASS, 5.000E+20 may suffice for starters.
Run the simulation according to Running SOLPS-ITER. I'd recommend running it inside an interactive QSUB job,
```
qsub -IX -l select=1:host=soroban-node-03
```
...and checking its progress using 2dt nesepm. If your separatrix density starts falling steeply, soft-land the simulation,
```
touch b2mn.exe.dir/.quit
```
...and increase the gas puff intensity. If the separatrix density starts rising, decrease the gas puff in the same way. Pure-deuterium, moderate-density COMPASS simulations take this treatment quite well and do not attempt to crash. Settle on a gas puff value which yields your desired separatrix density.

B2.5 Standalone Case

Gas puff for a standalone case can be implemented using boundary conditions for the continuity equation BC (BCCON) in the b2.boundary.parameters file. Use BCCON = 8 to prescribe a total particle flux with constant flux density in the CONPAR value, particles/s - that is your gas puff intensity. Apply this BC to your chosen neutral fluid at an outer radial boundary (type N or S in SOL or PFR regions). There are two options for choosing the boundary:

Implement a simplified "global" puff using the existing far SOL boundary (type N), spanning most of the wall outline. This is the easy way, where you only need to change BCCON and CONPAR for an already existing boundary. Depending on the order of fluid species (here index 0) and boundaries (here index 6), the end result could look like this:
b2.boundary.parameters
```
BCCON(0,6)= 8,  9,
! ...
CONPAR(0,6)= 1.000E+19,  0.0,
```
(Untested) Implement a localized puff by adding a new boundary definition. This requires some work and you should know what you are doing. Increment NBC, specify BCCHAR, and use e.g. the BC_LIST_SIZE, BC_LIST_X, BC_LIST_Y to select your gas puff cell. Append values to the BCCON and CONPAR arrays as needed. Leave the BCs for other equations without modifications, which will result in the "no condition" default behaviour.

Impurity gas puff

Once you've managed one gas puff, there is actually nothing stopping you from adding another one. The procedure is the same as for the main ion species gas puff, there is only a few things to consider:

In DivGeo, add another puffing stratum via Variables>Add>Gas puff.
Assuming you have already added your new Plasma species into the simulation, you should be able to configure Gas species for the gas puff as your new impurity species. Note that while D2 molecule is always automatically available, it is not necessarily the case for other species. You may need to use the atomic species instead.
Set Puffed flux to something 1-2 orders of magnitude smaller than the main species.
The strata will appear in b2.neutrals.parameters in the order that corresponds to their numbering in DivGeo. Look for the second C-type stratum.

The hardest part is to provide a sensible value for the impurity gas puff intensity. In contrast to the main species, there is no easy reference in terms of the total plasma density, since the impurity should only represent a small fraction of the plasma and it should be localized in the divertor region (assuming standard seeding experiment setup). More specialized diagnostics data are often needed, including neutral pressure in divertor, spectroscopy, etc. Furthermore, the second gas puff represents another degree of freedom in the model and typically makes it even more difficult to optimize against the experimental data. A good idea is to employing a feedback control scheme for one or both of the gas puffs, as described in the following section.

B2.5 Standalone Case

Follow the same procedure as for the main ion species gas puff in standalone case, but modify the boundary condition BCCON(is, ib) and parameter CONPAR(is, ib) for the impurity species neutral fluid (index is).

Feedback schemes (gas puff)

The problem with setting the gas puff intensity manually is that your separatrix density now depends on two free, unpredictable parameters: the diffusion coefficient $D_n$ (possibly in a profile) and the gas puff intensity (where previously $n_{e,core}$ could be kept constant). This is why it's a good idea to tune a case with gas puff using feedback schemes.

A feedback scheme is a mechanism where an active actuator (e.g. gas puff) is being adjusted in time to achieve a desired requirement (e.g. a value of separatrix density). In SOLPS-ITER, you construct your feedback scheme by independently configuring three things:

the requirement (NA_FEEDBACK_CHOICE, see below)
the actuator (NA_FEEDBACK_ACTUATOR, see below)
the feedback formula (NA_FEEDBACK_OPTION, see below), which takes the difference or ratio between requirement and current value, and calculates the next actuator value

Currently, all of the implemented actuators are focused on density changes in one way or the other (mainly gas puff and pumping). The implemented requirements cover a wider range, including total density of species at different locations, particle flux, relative concentration, radiated power, and more. In this case, we will focus solely on the gas puff actuator and the separatrix density requirement.

Examples of gas puff feedback schemes among the SOLPS examples include:

ITER_2588_D+He+N: feedback on the deuterium gas puff intensity to preserve "the total particle content for that species" (not clear whether deuterium ions or neutrals) summed over a given rectangle of B2.5 cells
ITER_2308_Honly_20MW: feedback on the core boundary hydrogen particle (not clear whether neutrals or ions) flux to preserve the neutral hydrogen particle flux through the core boundary

All available feedback schemes are documented primarily in the description of switches specified in the b2.feedback_control.parameters file (refer to the B2 input doc viewer open_in_new).

The NEW and OLD feedback scheme switches

Historically, there are two ways how to set up feedback schemes. You might run into a number of switches in b2mn.dat, which are documented as "feedback switches", e.g. b2stbc_isfeedback - those are the old-style switches and they are redundant in 3.0.8+. The new-style configuration of feedback is done almost entirely in the b2.feedback_control.parameters file.

Instructions

Set up a feedback scheme acting on the gas puff intensity (gas puff located in an outer midplane port) to reach a desired outer midplane separatrix density $n_{e,sep}$:

Procure a run with manual gas puff and make a copy of it. You needn't touch DivGeo this time.

In the run folder, create the main feedback configuration file b2.feedback_control.parameters with the following contents:

b2.feedback_control.parameters

&FEEDBACK_CONTROL
NA_FEEDBACK_CHOICE=     3,  ! the feedback aims for a specific OMP separatrix electron density
NA_FEEDBACK_TARGET=     2.0E+19,  ! requested OMP separatrix electron density

NA_FEEDBACK_ACTUATOR=   1,  ! the feedback controls the gas puff intensity
NA_FEEDBACK_IB=        -6,  ! actuating gas puff Eirene stratum 6 ('C'-type in b2.neutrals.parameters)
NA_FEEDBACK_ACTUATOR_MIN=   1.000E+19,  ! minimum gas puff intensity in particles/s
NA_FEEDBACK_ACTUATOR_MAX=   5.000E+21,  ! maximum gas puff intensity in particles/s

NA_FEEDBACK_OPTION=     1,  ! choose a particular formula for calculating next time step's gas puff intensity
NA_FEEDBACK_ALPHA=      0.5,  ! choose a mid-range speed for the gas puff intensity changes
/

Add these lines to b2mn.dat to enabled the feedback scheme:

b2mn.dat

'b2stbc_feedback'       '1'  ! turn on the feedback scheme
'eirene_ionising_core'  '1'  ! neutrals entering the core should return as ions through the core boundary

The userfluxparm value in b2.neutrals.parameters now specifies the initial gas puff intensity, until a b2.feedback_save.parameters file is written on the next simulation run.
Set the requested elapsed time in b2mn.dat to 60 s and run the simulation, starting from the previous manual gas puff solution.
```
rm b2mn.prt ; cp b2fstate b2fstati ; correct_b2yt_timestamps
b2run b2mn >& run.log &
```
If the simulation runs alright with no errors, you are halfway done.

B2.5 Standalone Case

Assuming a gas puff for the standalone case is implemented (see the previous sections), the procedure is almost identical to setting up a feedback scheme in Eirene-coupled case. Follow the same steps, then make two extra adjustments before running the simulation:

Change the gas puff BC to BCCON = 11 (feedback-enabled boundary condition). CONPAR now defines the initial gas puff intensity, until a b2.feedback_save.parameters file is written on the next simulation run. For example:
b2.boundary.parameters
```
BCCON(0,6)= 11,  9,
! ...
CONPAR(0,6)= 1.000E+19,  0.0,
```
Change the actuator target index NA_FEEDBACK_IB to a positive value that corresponds to the BC index of the gas puff in b2.boundary.parameters. For example:
b2.feedback_control.parameters
```
NA_FEEDBACK_IB= 6,  ! actuating gas puff boundary condition 6 (BCCON=11 in b2.boundary.parameters)
```

The rub with controlling plasma density through feedback on gas puff intensity is that getting it to work is only half the battle. The real challenge is to get the plasma solution you want. The main danger is that you may get an oscillatory solution, where the plasma density and total particle content oscillate with no end in sight. To learn our way around this problem, let's delve a little deeper into how the feedback scheme we've just set up works.

Feedback formula

Using NA_FEEDBACK_OPTION = 1 in b2.feedback_control.parameters translates, according to the B2.5 input documentationopen_in_new and my intuition, into the following formula:

\[\Gamma_{puff,new} = \Gamma_{puff,old} \cdot \frac{1 + \alpha \cdot \frac{\text{target } n_{e,sep}}{\text{current } n_{e,sep}}}{1 + \alpha}\]

The $\alpha$ parameter is set using NA_FEEDBACK_ALPHA in b2.feedback_control.parameters, and it controls the speed of the feedback control changes. The smaller it is, the more the rescaling factor approaches unity, and the slower $\Gamma_{puff}$ reacts to the ratio between the target and current $n_{e,sep}$. Its default value is $\alpha=0.001$, which is ridiculously tiny. For my COMPASS simulations, I have reached good results with $\alpha = 1.0$. In a complicated simulation where you can afford long run times, however, lower $\alpha$ presumably brings in more stability. In an oscillatory feedback solution, see below, $\alpha$ controls the period and magnitude of the oscillations.

Managing the oscillations

Oscillatory plasma solutions: theoretical introduction. As the article Kukushkin & Krasheninnikov, Bifurcations and oscillations in divertor plasma describes, a tokamak flux tube can act really weird if you inject enough power into it. As we know, SOL physics is pretty non-linear, especially around the target when $T_e$ drops to several eV. Small changes in the input parameters (or upstream parameters) can result in large changes in the solution (or target parameters). Such non-linearities can result in systems, such as the aforementioned flux tube with a high input power, which can sustain two particular solutions but nothing in-between. In the simplified 1D case presented in the article, the crucial quantity is the neutral pressure above the target $p_n$. Say, you are trying to reach a certain $p_n$ by gas puff feedback, which increases or decreases the number of particles $N_D$ in the flux tube, and that the desired $p_n$ is in the "forbidden zone". You start from a low-density, low-$p_n$ solution and increase the gas puff. Before the simulation reaches the desired $p_n$, the system non-linearity will cause a sudden plasma cool-down. $p_n$ will jump up way higher than you wanted it to. So you decrease the gas puff (or number of particles in the flux tube $N_D$), and $p_n$ obediently starts falling. But then, before you get close to the $p_n$ you want, the non-linearity kicks in again in the opposite direction. The target plasma heats up and $p_n$ drops below your desired value. So you increase the gas puff... Rinse, repeat. You now have sustained oscillations in the plasma solution.

The exact mechanism of the non-linearity which causes this behaviour is not explained by Kukushkin and Krasheninnikov. I assume it has to do with power and pressure loss processes in the divertor plasma, which kick in when target $T_e$ drops below 10 eV or so, and then leave just as abruptly when the target plasma heats up. However, I have seen sustained oscillatory solutions in SOLPS feedback cases where the target remained safely attached all the way through. K&K explain that in a 2D plasma, self-sustained oscillations can appear even without feedback, just by the virtue of neighbouring flux tubes sharing neutrals based on their $p_n$. Such solutions would oscillate all on their own, even without gas puff feedback. K&K claim that such oscillations, on the time scale of neutral diffusion, have indeed been observed at JET. So maybe it's that. Or maybe it's Maybelinne. The important thing is that K&K's explanation matches the feedback oscillations I have observed:

Solutions with low $n_{e,sep} \sim 2 \times 10^{18}$ m^-3 are stable and gas puff feedback converges correctly.
Solutions with high $n_{e,sep} \sim 2 \times 10^{19}$ m^-3 are likewise stable and gas puff feedback converges correctly.
Solutions with intermediate $n_{e,sep}$, even when started from a converged solution with the correct plasma parameters, begin oscillating.

Oscillatory plasma solutions: how they look and behave. After running a case with density control via gas puff intensity feedback, this is what you want to see (requested $n_{e,sep}=2 \times 10^{19}$ m^-3):

Damped oscillations

Nice and damped oscillations, converging toward the desired value .

This is what you do not want to see (requested $n_{e,sep}=1.3 \times 10^{19}$ m^-3):

Sustained oscillations

Grim sustained oscillations, centered around the desired value.

The value of the $\alpha$ parameter affects the oscillation period and magnitude.

Changing $\alpha = 0.1 \rightarrow 0.01$	Changing $\alpha = 0.1 \rightarrow 0.5$

Lower $\alpha$ increases the oscillations magnitude and decreases the period. I tried going as high as $\alpha=2.0$, at which point the oscillations of the solution were comparable to the simulation noise, but they were still present. It's up in the air if this sort of solution is as good as a converged one. The separatrix density is correct, and you can smooth out the oscillation with tools within SOLPS-ITER built just for these cases. But it feels risky. Wouldn't it be better to get a really converged solution?

How to force an oscillating solution to converge. Start from a high-density plasma solution. Simple as that. (Shut up, future me, I am allowed some optimism.) Props to Haosheng Wu, who mentioned this solution at the SOLPS Slack (channel q_and_a, April 2024). Set the desired $n_{e,sep}$ to a high value, let the simulation converge, then decrease $n_{e,sep}$ step-wise. I've gone from $n_{e,sep}$ 2e19 to 1.7e9 to 1.3e19 within two days, resolving oscillations which had stopped my work for two months.

Feedback adjustment

Three consequent simulations, gradually decreasing requested separatrix electron density. At 2.2 ms, $\alpha$ was decreased from 0.1 (I think) to 1.0 to prompt faster convergence.

I am not sure how this trick gets around the non-linearity described by Kukushkin and Krasheninnikov, but I'm just happy it works. On the flip side, it appears that relaunching a converged low-density simulation with different parameters may provoke renewed oscillations. In the worst case, always start from a high-density case.

Final word of hope. In COMPASS plasmas with low density (line-averaged density $\overline{n}_e < 6 \times 10^{19}$ m^-3), gas puff was usually turned off in the experiment. Initial working gas pressure was sufficient to sustain the required plasma density. It appears that such SOLPS plasmas work well with simple boundary conditions of fixed $n_{e,core}$, and don't require gas puff. Conversely, COMPASS plasmas with high density required additional gas puff in the experiment. This effect should be included in a realistic interpretative SOLPS-ITER simulation. Luckily, the separatrix density was so large that gas puff feedback schemes should have no problems converging. On other tokamaks, who knows. But I have faith in you.