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.

Multi-species plasma and impurities

construction This section is under construction.

Motivation

Real plasma in a tokamak includes small concentration of impurities from the divertor targets, limiters, and the wall due to erosion (C, W, Be, ...). Impurities may be also intentionally injected in large quantity to seed divertor detachment (N, Ne, Ar, ...). Reactor plasma will include complete fuel mixture (D and T) and reaction products (He). There are plenty of scenarios which require additional species in the simulation and plenty of complication to be introduced in the process.

General advice

It is generally better to set up and let converge a deuterium simulation before adding impurities to it. This is because SOLPS-ITER works better (and faster) in "baby steps", where the complexity is gradually increased. This applies in a variety of situations. When aiming for a simulation on a very fine grid, one may first set up a simulation on a coarse grid and then use b2yt to import the converged results to the new grid, where they will run more stably and faster than if one started from the "flat profiles" solution. When aiming for a simulation with impurities and drifts, one may first set up a drift-less deuterium simulation (with a high time step, such as 10^-4 s), then add impurities (reducing the time step to 10^-5 s) and finally switch on drifts (reducing the time step to 10^-7 s). For this reason, it is assumed in the present section Add impurities that you already have a converged D simulation and want to a run a simulation on the same grid but with more ion species.

Some information on adding impurities may be found in the SOLPS Conversion Tutorial hosted at the ITER SharePoint, section Adding new species to a SOLPS-ITER run. Here we will take a more beginner approach, which basically amounts to stealing the DivGeo file from your converged D simulation and building a whole new simulation.

Carbon sputtering

The divertor and limiter tiles of COMPASS (and many other tokamaks) are made of graphite, a low-Z material which conducts heat well and does not melt. Carbon is, however, sputtered, both physically and chemically.

During physical sputtering, ions or high energy neutrals hit the solid surface and dash out one or more atoms or ions from the surface. The yield of physical sputtering (number of atoms/ions sputtered per incident particle) depends on the energy of the incident particle, and is described by the Roth-Bogdansky formula in SOLPS-ITER [Roth & Carcia-Rosales 1996] [EIRENE manual].

During chemical sputtering, complex chemical bonds are formed between the target atoms and the hydrogen/deuterium/tritium atoms recycling at the target. For example, carbon and hydrogen may form methane, which then evaporates from the target and is quickly broken down and ionised in the plasma. It is difficult to describe this process analytically, so SOLPS-ITER uses a constant chemical sputtering yield $\gamma_{chem} \sim 0.01$ (number of atoms sputtered per incident ion/neutral).

Here's how to add physically and chemically sputtered carbon to a SOLPS-ITER simulation:

First, build and run the initial D simulation using the entirety of solps-doc. This will include:
1. setting up the grid,
2. choosing somewhat realistic boundary conditions (mainly $P_{SOL}$ and $n_{core}$),
3. running the simulation from flat profiles (or another pure-D simulation solution) until it converges (see convergence criteria in Running SOLPS-ITER: Checking if the run has converged), and
4. tweak the boundary conditions (mainly $D_\perp$ and $\chi_\perp$) until experiment-model agreement is achieved.
After carbon is added, the SOL will probably cool down due to increased radiation energy losses, so aim only for a rough experiment-model match. The desired output of this step is a b2fstate file containing a converged solution which is pretty close to the experiment. In particular, don't bother fine-tuning the diffusion coefficient profiles in H-mode. The profiles will fall apart after carbon is included.

Create a new simulation folder and its baserun and run subfolders.

cd $SOLPSTOP/runs
mkdir compass-D+C
cd compass-D+C
mkdir baserun
mkdir run

Copy the DivGeo file from your deuterium simulation into the carbon simulation.
```
cd baserun
cp ../../compass-D/baserun/compass.dg compass.dg
```
Make sure you've selected a device name (it's not needed when you don't run DivGeo) and run DivGeo in baserun.
```
setenv DEVICE compass
dg &
```
File>Open the copied DivGeo file.

Select everything (Ctrl+A). Go to Variables>General surface data and browse all of the variables. In particular, you'll be interested in the physical and chemical sputtering settings. To turn on both physical and chemical sputtering of carbon, set:

Variable	Value	Meaning
Wall Material	C	The first wall is made of carbon.
Wall Temperature	-0.026	Particles desorbed from the wall will have a Maxwellian velocity distribution with the mean energy 0.026 eV = 28 °C. (Discuss the exact value with your experimentalists.)
Phys. sput. model	2	Use the standard physical sputtering model.
Phys. sput. factor	1	Don't artificially enhance or reduce physical sputtering yield.
Sputtered atom	-1	Sputter atomic carbon from the wall and trace the sputtered atoms so they may enter the plasma.
Chem. sputter.	1	Use the standard chemical sputtering model.
Chem. sput. factor	0.01	(Probably) Set the chemical sputtering yield to 0.01.

Don't forget to click Set after changing every parameter! If you change the parameter values first and then hit Set on one of them, all the other parameter values will reset. It's a good check to close the dialogue window, reopen it and check all the parameter values are as you like them.

Do the same for the divertor targets. Unselect all (Ctrl+U) and go to Variables>Target specification>#1 and then #2 and set the variables to the same values as above. The target temperature is tricky; currently I use -0.05 = 307 °C. There are switches in b2mn.dat which can calculate the target temperature self-consistently from its heat load, thickness, material and temperature on the cooled side, but the default setting is that the target/wall temperature is supplied by the user.
Go to Variables>Add>Plasma species and set Species to C. Leave the boundary conditions as they are, you'll change them later.
Commands>Check variables, Commands>Rebuild Carre objects, File>Save, File>Output. You can close DivGeo now.

Set up links to this DivGeo file in baserun.

lns compass           # without the .dg part!

Run Carre.
```
carre -
```
For instructions, tips and troubleshooting, see the Creating a new run guide.
Run triang.
```
triang
```
For instructions, tips and troubleshooting, see the section Creating a new run guide.
Import the grids into your DivGeo file and check that they look ok. They don't need to exactly match the grids of the deuterium simulation. In theory, if you used the same input (for instance setting pasmin in Carre), the grids should be identical. In practice, how am I supposed to remember how I built a grid two years ago? (There should be some way to reuse the old grids but mehhhh.)

Set up the boundary conditions files in baserun as described in the Creating a new run guide. In particular,

cp b2ah.dat.stencil b2ah.dat
cp b2ai.dat.stencil b2ai.dat     #may be left out if b2ai.dat is already present
cp b2ar.dat.stencil b2ar.dat     #may be left out if b2ar.dat is already present

b2run b2ah
b2run b2ai
b2run b2ar

Set up boundary conditions in run. You can copy b2.transport.parameters and b2mn.dat from the deuterium simulation.

cd ../run
cp ../baserun/b2.boundary.parameters.stencil b2.boundary.parameters
cp ../baserun/b2.neutrals.parameters.stencil b2.neutrals.parameters
cp ../baserun/input.dat input.dat
cp /some_source/b2.transport.parameters b2.transport.parameters
cp /some_source/b2mn.dat b2mn.dat

Set realistic boundary conditions. If uncertain, refer to example runs which come shipped with SOLPS and check their boundary condition files line by line.
1. In b2.boundary.parameters, set the input power $P_{SOL}$ (BCENE(0) = BCENI(0) = 8, ENEPAR(1,1,0) = ..., ENIPAR(1,1,0) = ...) the same as in the deuterium simulation.
2. In b2.boundary.parameters, set the main ion $n_{core}$ (BCCON(0,1,1) = 1, CONPAR(0,1,1,1) = ...) the same as in the deuterium simulation.
3. In b2.transport.parameters, set the perpendicular transport parameters the same as in the deuterium simulation (for starters). Take note that there is now a different number of species in your simulation, so if you're using the 2*0.3 notation in b2.transport.parameters (see Questions and answers: I copied b2.transport.parameters from the manual...), change the number before * accordingly. We're adding carbon here, which is 7 extra species (C0, C1+, ..., C6+), so 2 goes to 9.
4. Set any additional deuterium/overall boundary conditions you might have used in the deuterium conditions (like the SOL fall-off length, heat flux limiters or target Mach numbers).
5. In b2.boundary.parameters, set the impurity BCCON. For starters, you may use BCCON(0,1,...) = 8 and CONPAR (0,1,1,...) = 0.0 - the core boundary condition is "set the particle flux value" and its value is zero. In experiment, there is probably an influx of C⁶⁺ from the core and outflux of C⁴⁺ and C⁵⁺ into the core, but this will do for now. At the divertor targets, leave the boundary conditions at 3, sheath conditions.
6. (Optional) Refer to an example run, such as the AUG D+He+C run, and copy the boundary conditions which seem reasonable.
7. Inspect b2mn.dat and set the switches as you like. In particular:
  1. Make sure to use the "physics" boundary conditions (switches b2stbc_boundary_namelist, b2tqna_transport_namelist and b2stbr_neutrals_namelist). If you're using a diffusion coefficient profile, the switches include b2tqna_inputfile.
  2. Set the time step (switch b2mndr_dtim).
  3. Set the number of iterations or simulation runtime (switches b2mndr_ntim and b2mndr_elapsed). I propose to start with a very short simulation (60 s) to get some results quickly and check that everything is alright.
  4. Set whether new tallies (used in the 2dt command and others) should be appended to old ones (switch b2mndr_stim).
  5. Choose between a coupled and a standalone run (switches b2mndr_eirene and b2mndr_rescale_neutrals_sources).
  6. Produce tallies and power+particle balance files (switches tallies_netcdf and balance_netcdf).
8. Check that all deuterium/overall boundary conditions match in your D and D+C simulation.
Set up carbon sputtering. In short, assuming a coupled (B2.5+EIRENE) run:

Open input.dat, find the section *** 6a. General data for reflection model and in lines which look like
Formatted (input.dat)JSON (eirene.input.json)
SURFMOD_1_CARBON_SPT_580K 1 12 -1 0 1
{ // ... "REFLECTION_MODELS": { "MANUAL": "http://www.eirene.de/eirene.pdf#subsection.2.6.1", // ... "SURFMODS": [ { "NAME": "SURFMOD_1_CARBON_SPT_580K", // ... "ISRC": 0, // ... } ] } }
...change the value of ISRC from 0 to -1. Do this for each of the surface models (SURFMOD) that have CARBON_SPT in their name.
Formatted (input.dat)JSON (eirene.input.json)
SURFMOD_1_CARBON_SPT_580K 1 12 -1 -1 1
{ // ... "REFLECTION_MODELS": { "MANUAL": "http://www.eirene.de/eirene.pdf#subsection.2.6.1", // ... "SURFMODS": [ { "NAME": "SURFMOD_1_CARBON_SPT_580K", // ... "ISRC": -1, // ... } ] } }
Mind the spaces! The numbers must stay vertically aligned.

Now comes the long explanation. My current (August 2021) state of knowledge is that both B2.5 and EIRENE can handle sputtering, and that these processes are controlled by two independent sets of switches.

Assuming your simulation is standalone, the appropriate switch set in b2mn.dat is (probably):
b2mn.dat
```
! Common sputtering switches
'b2stbr_sput_dst' '2' ! sputtered ion species index is 2 = neutral carbon

! Physical sputtering switches
'b2stbr_sput_phys' '1.00' ! turn on physical sputtering and set the yield to 1.00

! Chemical sputtering switches
'b2stbr_sput_chem_model' '0' ! use the empirical formula for yield calculation
'b2stbr_sput_frc' '0.01' ! probably chemical sputtering yield; overridden by chemical_sputter_yield in b2.neutrals.parameters
```
Additionally, in b2.neutrals.parameters you may specify chemical sputtering yield individually for each structure element with the chemical_sputter_yield array, which will probably come predefined in b2.neutrals.parameters.stencil. The array index corresponds to the element numbers in DivGeo. Look for further information in $SOLPSTOP/modules/B2.5/src/documentation/b2input.xml. I have unfortunately found no other sputtering documentation.

If your simulation is coupled, the appropriate switches are in input.dat. Contrary to what you might have heard, this file is human-readable, just very particular about each individual character. (I've heard stories of looking for badly placed spaces from too many people not to be afraid.) Find the section 6a. General data for reflection model, which may look something like this:
```
*** 6a. General data for reflection model
T
cpath modules/Eirene/Database/Surfacedata/TRIM/
D_on_C
C_on_C
1.00000E+00 1.00000E+00
1.00000E+00
1.00000E+00
1.00000E+00 1.66667E-01 1.66667E-01 1.66667E-01 1.66667E-01 1.66667E-01
1.66667E-01 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00
1.00000E+00 5.00000E+01 1.00000E-01 0.00000E+00 0.00000E+00 0.00000E+00
SURFMOD_1_CARBON_SPT_580K
 1    12    -1     0     1
1.20600E+03-5.00000E-02 0.00000E+00 0.00000E+00 0.00000E+00-1.00000E+00
1.00000E+00 1.00000E+00 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
SURFMOD_2_CARBON_SPT_1160K
 1     2     0     0     1
1.20600E+03-1.00000E-01 0.00000E+00 0.00000E+00 0.00000E+00-1.00000E+00
1.00000E+00 1.00000E+00 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
SURFMOD_3_CARBON_SPT_302K
 1    12    -1     0     1
1.20600E+03-2.60000E-02 0.00000E+00 0.00000E+00 0.00000E+00-1.00000E+00
1.00000E+00 1.00000E+00 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
```
Your interest lies with the three paragraphs beginning with SURFMOD. These define three groups of structure elements with common sputtering properties. In this case, the first one are the divertor targets (temperature approx 300 °C), the second one is a mystery and the third one is the first wall (temperature approx 30 °C). If you've done everything correct, there should only be two paragraphs, 580K and 302K. If one of them is missing, check the DivGeo surface data and build the grids again. What the numbers inside these paragraphs mean is (mostly) described in the EIRENE manual, section 2.6 Input Data for Surface Interaction Models. Every number in a SURFMOD paragraph denotes one variable, whose names are, in turn:
```
SURFMOD_modname
ILREF ILSPT ISRS ISRC
ZNML EWALL EWBIN TRANSP(1 ,N) TRANSP(2 ,N) FSHEAT
RECYCF RECYCT RECPRM EXPPL EXPEL EXPIL
RECYCS RECYCC SPTPRM ESPUTS ESPUTC
```
At first glance (in the August 2021 version of the EIRENE manual and SOLPS-ITER), you will notice that there are actually 5 numbers in the first paragraph line, not 4. In such cases, the source code is the ultimate judge. As can be found in the source code of EIRENE, file $SOLPSTOP/modules/Eirene/src/startup-routines/input.f, the fifth variable is called LCHSPNWL and it concerns switching between SOLPS4.3 and SOLPS-ITER version 3.0.1 input format. Otherwise, the variable values should match some realistic physical picture of the tokamak first wall or divertor target.

Sputtering is mainly controlled by the variables ILSPT, ISRS, ISRC, RECYCS and RECYCC.
- ILSPT (here 12) is composed of two numbers, M = 1 and N = 2. It controls the sputtering model for chemical sputtering (M = 1 stands for "constant chemical sputtering rate") and physical sputtering (N = 2 stands for "modified Roth–Bogdansky formula for sputter yield").
- ISRS (here -1) controls physical sputtering. Its value here means, put simply, that physical sputtering is on and produces atoms of the same type as the wall/target material - carbon.
- ISRC (here 0) controls chemical sputtering. Its value here means that chemical sputtering is off. This happens despite our previous effort to turn chemical sputtering on in DivGeo, and was the source of much anguish to me once. Set this number to -1 in the first and third paragraph to turn chemical sputtering on and set it to produce atoms of the same type as the wall/target material - carbon. Important: substitute 0 for -1, including the space at the beginning. The reason is that input.dat is very peculiar about its blank spaces. "Shifting" the rest of the line by one character will result in the run.log error Fortran runtime error: Bad value during integer read, as Fortran reads - as the next input value instead of -1.
- RECYCS (here 1.00000E+00) can mean different things depending on which physical sputtering model was chosen in ILSPT. In our case, RECYCS is (probably) the physical sputtering yield, an additional multiplier to the atom flux calculated by the Roth–Bogdansky formula.
- RECYCC (here 1.00000E+00) can mean different things depending on which chemical sputtering model was chosen in ILSPT. In our case, RECYCC is (probably) the chemical sputtering yield.
My current knowledge of sputtering and the appropriate choice of sputtering yields is limited. Until further update, use RECYCS = 1.0 and RECYCC = 0.01 and consult the EIRENE manual and references therein (such as [Roth 1997] and [Roth 1999]) for further details.
Go to the simulation folder and correct the time stamps.
```
cd ..
correct_baserun_timestamps
```
Go to the run folder and link the baserun files here.
```
cd run
setup_baserun_eirene_links
```
Check if the simulation will run.
```
b2run -n b2mn
```
Inspect the output and look for error messages.

Copy the results of the deuterium simulation as your starting state.

cp ../../compass-D/run/b2fstate b2fstati_deuterium
cp b2fstati_deuterium b2fstati

Correct the time stamps.
```
correct_b2yt_timestamps
```
Run the simulation.
```
b2run b2mn >& run.log &
```
After 60 seconds, check that the simulation has produced output. If everything went ok, congrats, you now have a simulation with sputtered carbon!