Unphysical Sawtooth crash effects

[AlexSaperstein]

Sawtooth crashes are having unphysical effects on parameters like Ip, beta_N, tau_E, and li3 (and potentially more), like that shown below. li3 should change as after a sawtooth collapse occurs, but it should transiently decrease, not increase.

[abushanisro]

In sawtooth_solver.py:109–111, it looks like psidot * dt from the pre-crash (peaked) state is added to the flattened psi:

evolved_psi_redistributed_value = (
redistributed_core_profiles.psi.value
+ core_profiles_t.psidot.value * dt
)

Could this mixing of old psidot and new psi explain the unphysical Bp² -> Wpol -> li3 changes?
Also, does Ip_total stay constant through crashes? The plots seem to show small oscillations

Maybe plotting Wpol and Ip before/after the crash could help pinpoint the source

[AlexSaperstein]

Could this mixing of old psidot and new psi explain the unphysical Bp² -> Wpol -> li3 changes?

Potentially, but I don't really know for sure.

Also, does Ip_total stay constant through crashes? The plots seem to show small oscillations

Are you referring to in reality (the black curves) or in the sim (blue curves)? In reality, Ip_total shouldn't really be changing over the sawtooth crash, so it's unrealistic that in the simulation it does. In the simulation, Ip_total is NOT being kept fixed explicitly. Instead, a fixed choice for vloop_lcfs is being used to constrain the total current, allowing Ip_total to evolve with the system.

I have noticed since raising this issue that this bug does seem to be related to using vloop_lcfs to constrain Ip instead of Ip_total. I reproduced the above plot below, but this time plotting the same exact simulation again, with the sole exception that Ip is being constrained by ip_total instead of vloop_lcfs (orange curves), i.e. Ip_total is kept constant throughout the sim. Here you can see that the bug is no longer present. Hopefully this will help!

[jcitrin]

I have noticed since raising this issue that this bug does seem to be related to using vloop_lcfs to constrain Ip instead of Ip_total

Thanks that's a great hint.

We will fix this soon. There are unfortunately some other priorities at the moment.

[AlexSaperstein]

So I've run into a potentially more significant version of this bug while trying to validate collapses on more shots. I was initially under the interpretation that this was a transient, single-time-step bug that resolves itself on the following time-step without any effect on the subsequent simulation. However, after simulating the collapse of shot 1120210022 on C-Mod, I have found that this may not be the case. Below are the time-traces for some relevant collapse parameters for this shot, with black traces corresponding to the real experiment, and blue corresponding to the TORAX simulation of it.

You can see that the ST bug that transiently spikes up following a ST crash persists for many time-steps in this case. During this time, this has a significant effect on the ohmic heating power within the plasma (especially the edge), which transiently spikes the edge Te (see below) and drops Prad.

The change in Ip also appears to be largely localized in the edge (see below). And not only that, the persistent change in Ip is a positive one, as opposed to a negative one like we have seen for all other previous (and benign) manifestations of this bug. Although, this positive persistent change is notably preceded by a negative single-step change, similar to the others.

It is unclear at the moment how prevalent this more significant, persistent manifestation of the bug is, and under what conditions this occurs, instead of the benign versions I've been seeing before. So I don't know if this should really change the priority of this issue at the moment. However, it may be useful for narrowing down the root cause of this bug; once again pointing towards a potential incompatibility between the ST model and the fixed vloop_lcfs constraint.

[theo-brown]

I've been reading some of the JETTO documentation about the Vloop boundary condition, looking for any info on how they navigate the interaction with sawteeth. A couple of quotes (may be obvious to those with more tokamak physics background than me!):

Getting the loop voltage to match experiment is extremely sensitive to the electron temperature near the edge

[with the loop voltage boundary condition,] edge surface currents will be applied each time step in order to maintain the desired total plasma current. These edge surface currents can fluctuate strongly in time if other quantities in the simulation are changing (particulary near the boundary).

These sound pretty consistent with what you're observing in your last two plots.

There's also a little section on how they handle splitting of the flux consumption with sawteeth models, with a reference to Houlberg 1987, but I'm not sure how directly relevant this is.