Swap plane Poisson solver with slab Poisson solver in GK IWL simulations & 2x LBO energy conservation

[manauref]

This goes with PR 93 of gkylcas

Problem

Currently the IWL simulations are using the Poisson solver on planes (fem_poisson_deflated). This has been shown to break energy conservation. One reason we have continued to do that is that:
a) We didn't have the ability to bias the limiter corner in the slab Poisson solver (fem_poisson_perp).
b) We didn't have a working workflow that used the slab Poisson solver and twist-shift BCs (TS BCs).

Here we address both of these points, and deprecate the use of fem_poisson_deflated in favor of fem_poisson_perp.

Solution

Biasing

We implemented the ability to bias the limiter corner in fem_poisson_perp. Some info is in the DR #898 .

Unlike biasing in fem_poisson, where can bias whole planes in the solver, in fem_poisson_perp we bias lines. That's at the moment we just need to bias the lines at (x,z)=(x_LCFS,z_min) and (x,z)=(x_LCFS,z_max).

Tests with biasing were added to the fem_poisson_perp unit test.

Workflow with TS BCs

We implemented a new set of steps to make sure the potential is twist-shift periodic with this new slab solver. These are:

1. Smooth the charge density, without BCs (fem_parproj).
2. Solve the perpendicular Poisson problem (fem_poisson_perp).
3. In the core, apply TS BC to phi at the lower z-boundary, and a ghost-from-skin-surf BC at the upper z-boundary.
4. Smooth phi, with Dirichlet BCs that take the Dirichlet value from the ghost in the core, and witih Dirichlet BCs that take the Dirichlet value from the skin in the SOL.

Commentary on alternatives

• Regarding step 1: It's unclear whether BCs are needed here, or what they should be. In SOL simulations we don't apply BCs, and that's the option that makes the operator self-adjoint and preserves energy. Perhaps indirectly using TS BCs as done in steps 3-4 would be worthwhile, but doesn't seem necessary so far.
• Note that steps 3-4 only apply TS BC at the lower boundary. This is as it's done in main now. We tried, both in main and in this branch, but applying TS BC only at the upper boundary or in combination with the lower, does not work.

Additional perks

We noticed that the LBO collision operator was not conserving energy in 2x. We fixed that in this branch.

This fix is not the final one, though. Through more testing I found scenarios (large gradients + normNu) in which the LBO still doesn't conserve energy. I think we need to compute nu-weighted moments to truly conserve energy, which we may pursue in another branch.

Tests

This branch is valgrind and compute-sanitizer clean.

Regression tests

Regression tests run, but they are expected to yield slightly different solutions since we are using a different algorithm now.

rt_gk_d3d_iwl_2x2v_p1

rt_gk_tcv_iwl_adapt_source_2x2v_p1

rt_gk_d3d_iwl_3x2v_p1

rt_gk_tcv_iwl_adapt_source_3x2v_p1

Production simulation

We ran a TCV case provided by @Antoinehoff in both main and this branch. Here are some snapshots.

[Antoinehoff]

Here are one-to-one comparison of production like simulations between the current state of main (5506bbe) and this branch (8fee1f5). The simulation setup is GK TCV iwl 3x2v, see rt_gk_tcv_iwl_adapt_source_3x2v_p1.c, for a coarse resolution (24x16x12x12x6) as presented in Hoffmann et al. 2026. The only parameter alterations are an increase of the input power, 0.5MW instead of ~0.25MW, to accelerate turbulence development, and a reduction of the collisionality scaling factor, $\nu_{frac}=0.5$ instead of $1.0$, to increase the time step.
The results from the main and this branch are labelled main and poissonperp, respectively.

Edit: I updated the plots and analysis with longer runs as a noticeable difference rise when approaching the quasi-steady state.

Note: The poisson perp run crashed at t~1320mus, the main run did not reach that time yet so it is unclear if the crash is solely due to the solver. I've observed that reducing nu_frac can destabilize the simulation at longer time.

General numerics comments

The poissonperp run presents a very stable time step ~4.6ns against a fluctuating one for the main run, ~3 to 4.6ns. The poissonperp simulation is consequently faster, achieving ~400mus against 330mus for main, in 6h and 4 GPUs.

Energy

It seems that the poissonperp leads to a lower confined energy state.
We see a strong reduction of the ion integrated Hamiltonian because of the global reduction of the potential value. This does not affect the total Hamiltonian. However, the electron Hamiltonian is lower in poissonperp than in main despite the stronger negativity of the potential, which indicates that the thermal energy of the electrons is strongly reduced in the poissonperp.
main

poissonperp

Potential

Poisson perp saturates into a lower potential values everywhere, increasing the ExB shear at the inner radial boundary and decreasing tit at the outer one.
main

poissonperp

Density

Fairly similar, no comments.
main

poissonperp

Temperature

The electron quasi-steady state temperature is remarkably reduced with the Poisson perp solver as it is expected from the energy analysis above. The SOL electron temperature matches better with experiment but it is unclear for the core (recall that the initial conditions of this simulation is ~ the experimental profile measurement). For the ions, a stronger temperature "well" is observed in the Poisson perp simulation.
main

poissonperp

Density fluctuations

We see that the fluctuations at the inner radial boundary are reduced in the poisson perp case, which may be due to the increase ExB shear in that region (potential well is larger in poisson perp).
main

poissonperp

Temperature fluctuations

It seems that the y-flows are stronger in the main simulation.
main

poissonperp

[Antoinehoff]

Thanks for this great work!!

[Antoinehoff]

Longer runs show significant discrepancies in the quasi steady state, who's right? The above analysis is updated.

[manauref]

Longer run are a bit worrying: everything looks like it is loosing energy. I think we need to see the steady state of this simulation Also I need help to plot the energy balance.

can you edit the comment adding the equivalent plots with the simulation run in main? i can message you privately about the energy balance problem

[Antoinehoff]

Longer run are a bit worrying: everything looks like it is loosing energy. I think we need to see the steady state of this simulation Also I need help to plot the energy balance.

can you edit the comment adding the equivalent plots with the simulation run in main? i can message you privately about the energy balance problem

I just restarted it, it is not as advanceed yet but here is how the integrated moments evolved for now. I'll update previous comment tomorrow.

[manauref]

@Antoinehoff one concern you had is that your energy conservation time trace looked worse in this branch, right? But the adaptive sources aim to conserve M2 energy, not H energy, right?