Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
355 changes: 9 additions & 346 deletions README.md

Large diffs are not rendered by default.

403 changes: 403 additions & 0 deletions docs/Contribute.md

Large diffs are not rendered by default.

35 changes: 16 additions & 19 deletions docs/feature_blog/Diffusion_coefficients.md

Large diffs are not rendered by default.

10 changes: 5 additions & 5 deletions docs/feature_blog/Drifts.md
Original file line number Diff line number Diff line change
Expand Up @@ -43,7 +43,7 @@ V_\perp^{(dia)} = -\frac{1}{enB} \frac{1}{h_y} \frac{\partial (nT_i)}{\partial y
V_y^{(dia)} = \frac{B_z}{enB^2} \frac{1}{h_x} \frac{\partial (nT_i)}{\partial x}.
$$

This term is stupid when you put it inside a divergence, because first you differentiate the ion pressure $nT_i$ and then you have to differentiate it *again*. This causes numerical instabilities. Luckily, as [[Bufferand, 2016]](https://www.sciencedirect.com/science/article/pii/S2352179116301946)<span class="material-symbols-outlined">open_in_new</span> explains in formulas (11) and (12), the diamagnetic flux can be decomposed into two components: the so-called magnetization flux and the grad-B drift of the ion guiding centres. The magnetization flux divergence is, by definition, zero, and the second term in grad-B drift is usually negligible. This means that *inside the divergence*, the diamagnetic drift can be replaced by a part of the **grad-B drift**:
This term is problematic when you put it inside a divergence, because first you differentiate the ion pressure $nT_i$ and then you have to differentiate it *again*. This causes numerical instabilities. Luckily, as [[Bufferand, 2016]](https://www.sciencedirect.com/science/article/pii/S2352179116301946)<span class="material-symbols-outlined">open_in_new</span> explains in formulas (11) and (12), the diamagnetic flux can be decomposed into two components: the so-called magnetization flux and the grad-B drift of the ion guiding centres. The magnetization flux divergence is, by definition, zero, and the second term in grad-B drift is usually negligible. This means that *inside the divergence*, the diamagnetic drift can be replaced by a part of the **grad-B drift**:

$$
\tilde{V}_\perp^{(dia)} = \frac{T_iB_z}{eb_z} \frac{1}{h_y} \frac{\partial}{\partial y} \left( \frac{1}{B^2} \right)
Expand All @@ -55,13 +55,13 @@ $$

This is great news! The magnetic field $B$ derivatives don't change between time steps because SOLPS-ITER is an electrostatic code, so you can just carefully precalculate them and get rid of the double differentiation. Or so they thought... until they tried modelling H-mode.

In formulas (7-9), [[Rozhansky, 2009]](https://iopscience.iop.org/article/10.1088/0029-5515/49/2/025007)<span class="material-symbols-outlined">open_in_new</span> describes the final transformation. It uses physics of the magnetic equilibrium. In a tokamak, the grad-B drift is vertical and of opposite direction for electrons and ions, so it causes vertical charge separation. The resulting electric field drives electric currents along the magnetic field lines. These are the **Pfirsch-Schlütter currents**, and their divergence is the same (taken negatively) as the divergence of the grad-B drift.
In formulas (7-9), [[Rozhansky, 2009]](https://iopscience.iop.org/article/10.1088/0029-5515/49/2/025007)<span class="material-symbols-outlined">open_in_new</span> describes the final transformation which was needed to facilitate drift simulations. It uses physics of the magnetic equilibrium. In a tokamak, the grad-B drift is vertical and of opposite direction for electrons and ions, so it causes vertical charge separation. The resulting electric field drives electric currents along the magnetic field lines. These are the **Pfirsch-Schlütter currents**, and their divergence is the same (taken negatively) as the divergence of the grad-B drift.

$\nabla( \textbf{j}_{grad-B} + \textbf{j}_{PS}) = 0$

Replacing the grad-B drift with Pfirsch-Schlütter currents inside divergences improves convergence further. This is why you'll encounter terms with "PSch" labels in the B2.5 equations, and why there are tildes ($\tilde{\Gamma}$) all over the place.
Replacing the grad-B drift with Pfirsch-Schlütter currents inside divergences improves convergence further. This is why you'll encounter terms with "PSch" labels in the B2.5 equations, and why there are tildes ($\tilde{\Gamma}$) all over the place. The more tildes, the more numerical tricks have been applied to a quantity.

Agonisingly, the physical nature of these three phenomena is all different. Diamagnetic flows are cross-field and fluid, the grad-B drift is cross-field and for the guiding centre, and P-Sch currents are parallel. If you calculate them on their own, their magnitude and direction is all different. The only thing they have in common is that their *divergence* is the same. So, keep in mind: inside divergences, diamagnetic, grad-B and Pfirsch-Schlütter are interchangeable, but if you're calculating real particle fluxes (like for energy convection), you have to use the entire diamagnetic drift formula.
Agonisingly, the physical nature of these three phenomena is all different. Diamagnetic flows are cross-field and fluid, the grad-B drift is cross-field and for the guiding centre, and Pfirsch-Schlütter currents are parallel. If you calculate them on their own, their magnitude and direction is all different. The only thing they have in common is that their *divergence* is the same (save for the sign). So, keep in mind: inside divergences, diamagnetic, grad-B and Pfirsch-Schlütter are interchangeable, but if you're calculating real particle fluxes (like for energy convection), you have to use the entire diamagnetic drift formula.



Expand Down Expand Up @@ -117,7 +117,7 @@ It's hard to say what are sensible values (1.5, 1.1, 1.0001) for any of the ramp

## Modify boundary conditions in `b2.boundary.parameters`

This is the most complicated step, but at the same time, it might not be important to make the simulation stable (Honza's conjecture). It is, however, important for the simulation to be physically correct. Most of the information here is based on the documentation of the boundary conditions and on the official example `ITER_2588_Donly_standalone_drifts` from the `solps-iter/examples` directory. Refer to that example for more details. See the [B2.5 switches](/solps-doc/extras/b2input) for a documentation on the boundary conditions. The main idea is:
This is the most complicated step, but at the same time, it might not be important to make the simulation stable (Honza's conjecture). It is, at any rate, important for the simulation to be physically correct. Most of the information here is based on the documentation of the boundary conditions and on the official example `ITER_2588_Donly_standalone_drifts` from the `solps-iter/examples` directory. Refer to that example for more details. See the [B2.5 switches](/solps-doc/extras/b2input) for a documentation on the boundary conditions. The main idea is:

- There are special versions of the sheath boundary conditions that are modified to properly account for drifts.
- It is advisable to use leakage conditions instead of decay lengths for the radial boundaries. But I'm not sure if that is important for drifts or just a good idea in general.
Expand Down
6 changes: 4 additions & 2 deletions docs/feature_blog/Gas_puffing_and_pumping.md
Original file line number Diff line number Diff line change
Expand Up @@ -343,7 +343,8 @@ Assuming a gas puff for the standalone case is implemented (see the previous sec
At low plasma densities, enabling density control using gas puff feedback can lead to oscillations in the plasma solution:

![Sustained oscillations](../img/feedback_sustained_oscillations.png)
<p align="center"><i>Sustained oscillations of the separatrix electron density, centered around the desired value.</i></p>

<center><i>Sustained oscillations of the separatrix electron density, centered around the desired value.</i></center>

This section discusses how such oscillations come to be and how to control them.

Expand All @@ -366,7 +367,7 @@ The exact mechanism of the non-linearity which causes this behaviour is not expl
**How oscillatory plasma solutions 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<sup>-3</sup>):

![Damped oscillations](../img/feedback_damped_oscillations.png)
<p align="center"><i>Nice and damped oscillations, converging toward the desired value .</i></p>
<center><i>Nice and damped oscillations, converging toward the desired value.</i></center>

The sustained oscillations shown in the beginning of this section are what you do not want to see (requested $n_{e,sep}=1.3 \times 10^{19}$ m<sup>-3</sup>).

Expand All @@ -381,6 +382,7 @@ Lower $\alpha$ increases the oscillations magnitude and decreases the period. At
**How to force an oscillating solution to converge?** Start from a high-density plasma solution. Simple as that. 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.

![Feedback adjustment](../img/feedback_adjustment.png)

*Three consequent simulations, gradually decreasing requested separatrix electron density. At 2.2 ms, $\alpha$ was increased from 0.1 to 1.0 to prompt faster convergence.*

I am not sure how (or if) this trick gets around the non-linearity described by Kukushkin and Krasheninnikov. On occasion, relaunching a converged low-density simulation with different parameters renews the oscillations. In the worst case, always start from a high-density case.
6 changes: 3 additions & 3 deletions docs/feature_blog/Impurities.md
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@ This feature blog entry covers:
**Do you really need impurities?** Each ionisation state of the added impurity becomes an ion species to keep track of, increasing output file size, simulation complexity and runtime. Think carefully if you need the impurity.

/// hint | Example impurity: tungsten
Tungsten is a key impurity in core physics. Even a small concentration can kill a fusion reactor by decreasing the energy confinement time. The source of tungsten in tokamaks is sputtering from the walls, so one would think modelling tungsten with SOLPS-ITER is key. Not so. Tungsten concentrations in the SOL are typically so small that they don't affect the plasma solution. Tungsten has 75 charge states including neutrals, so it's quite impossible to model in SOLPS without charge bundling. And finally, tungsten concentrations/densities are so small in tokamaks that it's in the kinetic regime rather than highly collisional, so Braginskii equations are terrible at describing its transport. You are much better off describing tungsten sputtering and transport with a dedicated kinetic code such as ERO2.0. SOLPS-ITER can provide ERO with the background plasma, which does not contain tungsten, and ERO can provide SOLPS with the power radiated by tungsten in the core, which translates into decreased power crossing the separatrix $P_{sep}$.
Tungsten is a key impurity in core physics. Even a small concentration can kill a fusion reactor by decreasing the energy confinement time. The source of tungsten in tokamaks is sputtering from the walls, so one would think modelling tungsten with SOLPS-ITER is key. Not so. Tungsten concentrations in the SOL are typically so small that they don't affect the plasma solution. Tungsten has 75 charge states including neutrals, so it's quite impossible to model in SOLPS without charge bundling. And finally, tungsten concentrations/densities are so small in tokamaks that it's in the kinetic regime rather than highly collisional, so Braginskii equations are inadequate at describing its transport. You are much better off describing tungsten sputtering and transport with a dedicated kinetic code such as ERO2.0. SOLPS-ITER can provide ERO with the background plasma, which does not contain tungsten, and ERO can provide SOLPS with the power radiated by tungsten in the core, which translates into decreased power crossing the separatrix $P_{sep}$.
///

**Baby steps.** SOLPS-ITER works better (and faster) if the simulation complexity is increased gradually. Even when you're aiming to make a simulation with impurities from the start, begin by building a pure-D simulation. Let it converge, implement self-consistent particle balance using [gas puffing and pumping](Gas_puffing_and_pumping.md), let it converge, and only then add impurities. You'll save runtime and headache with divergence.
Expand All @@ -31,7 +31,7 @@ You can use the initial plasma solution `b2fstati` without an impurity (but othe


/// hint | More resources
- [SOLPS Conversion Tutorial](https://iterorganization.sharepoint.com/:w:/r/sites/SOLPS-ITER/Shared%20Documents/General/SOLPS-ITER/Manuals%20and%20Documentation/SOLPS_Conversion_tutorial.docx?d=wef2a123c427245728cf76c7b4a0b3542&csf=1&web=1&e=a7aEwQ)<span class="material-symbols-outlined">open_in_new</span>, section *Adding new species to a SOLPS-ITER run*
- [SOLPS Conversion Tutorial](https://iterorganization.sharepoint.com/:w:/r/sites/SOLPS-ITER/Shared%20Documents/General/Tutorials/2017_Bonnin_Conversion_tutorial.docx?d=w665bc959a164494eacab09fe6b4fb948&csf=1&web=1&e=CU7cNc)<span class="material-symbols-outlined">open_in_new</span>, section *Adding new species to a SOLPS-ITER run*
///


Expand Down Expand Up @@ -89,7 +89,7 @@ Import the grids into your DivGeo file and check that they look ok. They don't n
/// warning | Use the `stencil` boundary conditions files
Copy over the `stencil` files generated in the `baserun` rather than boundary condition files from your old, carbon-less simulation. The number of ion species is different, and though checking the text files line by line is a drag, typing in all the new ion species is even a bigger drag.

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...](../supplementary/Questions_and_answers.md#transport-parameters-12s)), change the number before `*` accordingly. We're adding carbon here, which is 7 extra species (C<sup>0</sup>, C<sup>1+</sup>, ..., C<sup>6+</sup>), so `2` goes to `9`.
If you're using the `2*0.3` notation in `b2.transport.parameters` (see Q&A: [I copied `b2.transport.parameters` from the manual...](../supplementary/Questions_and_answers.md#transport-parameters-12s)), change the number before `*` accordingly. We're adding carbon here, which is 7 extra species (C<sup>0</sup>, C<sup>1+</sup>, ..., C<sup>6+</sup>), so `2` goes to `9`.
///

Make a dry run of SOLPS-ITER to check that everything is working. If no errors (especially regarding boundary condition files) appear, you're good to go and regulate the sputtering source.
Expand Down
Loading