Ensure the βMHD / βkinetic error is below a pre-determined threshold

[orso82]

Asked by @daveweisberg

• What is the best practice to ensure the βMHD / βkinetic error is below a pre-determined threshold (e.g. 1%)?​
• Is this a different process than converging on pressure / current error?​

[orso82]

@daveweisberg please note that the ActorEquilibriumTransport controls the convergence error with two parameters https://github.com/ProjectTorreyPines/FUSE.jl/blob/d51f84bd73f380fcb1d0c45a07a7a2cefe7ce94b/src/actors/compound/equilibrium_transport_actor.jl#L8-L9

I suspect that the cases you do not see converging to tolerance are terminated sooner because they reach the maximum number of iterations. It would be good to isolate one of the problematic cases to understand why this happens. One option may choose to raise an error if the desired convergence is not reached by the maximum number of iterations.

We may also want to look at more sophisticated iteration schemes, or even just add relaxation to current scheme to allow damping of oscillations in the solutions.

[daveweisberg]

I've published a new development branch with my optimization workflow that shows this problem. Please see https://github.com/ProjectTorreyPines/FUSE.jl/issues/373 for more details. In a nutshell, I have not been able to reduce the size of this betaN error either by changing the EquilibriumTransport convergence_error or max_iter, or by changing the equilibrium actor used in the workflow.

[TimSlendebroek]

@daveweisberg

Thanks for uploading the notebook with outputs to the PR as this makes it easy to look at :)

The Equilibrium-transport actor converges on Jtor(rho) and the kinetic pressure(rho) there is no option to convergence on /betaN

The convergence error was set to 0.01 which was reached in 4 iterations for the case you showcased. Which looking at the digest is actually pretty good (pressure and Jtor profiles as well as the fluxsurfaces):

I am not entirely sure what isn't working as both betas converged to <1% ? maybe it wasn't that exact case you had an issue with but during optimization?

[daveweisberg]

@TimSlendebroek @bclyons12 Yes, the pressure & current profiles have no problem converging quickly. However, what I'm concerned about is the betaN_MHD / betaN_tot error. In the case you posted betaN_MHD = 4.71 and betaN_tot = 4.56, which is an error of about 3%. Other cases can have errors of more than 10%. I don't understand why the equilibrium solver cannot exactly match the target betaN, or how to force it to converge to some specified betaN tolerance.

[TimSlendebroek]

The equilibrium solvers (Cheese and tequila) take in P, J and the boundary shape. There is no betaN normal target or tolerance in the equilibrium solvers nor in the transport part, could you point me to where you found this?

BetaN is more a derived quantity hence it's set as an expression in IMAS but I am also a bit confused as BetaN should at least be a proxy for convergence as it depends on quantities that should convergence however taking ratios might obfuscate that. BetaMHD and betaN from kinetic pressure might be slightly different depending on the last actor but for the cases where you see 10% difference between the last and second to last iteration with P and J converging is very concerning. Can you share that case's notebook?

[orso82]

The issue is that the Solovev actor does not match the required beta_n but only the core pressure https://github.com/ProjectTorreyPines/FUSE.jl/blob/master/src/actors/equilibrium/solovev_actor.jl#L96

One possibility could be to make beta_n with the Solovev linear pressure be equal to beta_n using the pressure profile defined in eqt.profiles_1d.pressure vector. For this we need a way to take the volume average of an vector in the MXHequilibrium package. @lstagner suggested writing a function g(r,theta,zeta) = f_spline(rho_spline(r)) and then use average(g, S::PlasmaShape,:volume).

[daveweisberg]

@TimSlendebroek If you replace cell 21 with the following inputs, you will get a lower betaN solution with about the same magnitude betaN error: Input:

ini.equilibrium.B0 = 5.5
ini.equilibrium.ip = 8e6
ini.equilibrium.R0 = 5.0
ini.core_profiles.greenwald_fraction = 1.2
ini.core_profiles.greenwald_fraction_ped = 0.85
ini.equilibrium.pressure_core = 0.5e6

act.ActorHFSsizing.do_plot = true
act.ActorHFSsizing.verbose = true
act.ActorEquilibrium.model = :TEQUILA

Output:

Jtor0_before = 2.09 MA
P0_before = 467.27 kPa
βn_MHD = 1.72
βn_tot = 1.51
Te_ped = 1.75e+03 eV
rho_ped = 0.9500

Jtor0_after = 1.88 MA
P0_after = 467.27 kPa
βn_MHD = 1.40
βn_tot = 1.54
Te_ped = 2.48e+03 eV
rho_ped = 0.9288
[ Info: Iteration = 1 , convergence error = 0.04364, threshold = 0.01

Jtor0_after = 1.95 MA
P0_after = 467.27 kPa
βn_MHD = 1.40
βn_tot = 1.54
Te_ped = 2.49e+03 eV
rho_ped = 0.9294
[ Info: Iteration = 2 , convergence error = 0.00613, threshold = 0.01

Here, the betaN error is ~10%.

[daveweisberg]

@orso82 I understand the limitations of Solovev, but this issue also occurs if I use TEQUILA or CHEASE. For example, the case I just posted above uses TEQUILA and has a betaN error of ~10%.

[daveweisberg]

@TimSlendebroek @bclyons12 @orso82 Update: I've checked the profiles from the MHD (dd.equilibrium.time_slice[].profiles_1d) and kinetic (dd.core_profiles.profiles_1d[]) parts of the data structure and have noticed several discrepancies:

The MHD and kinetic pressure profiles don't align when using the rho_tor_norm coordinate, while the MHD and kinetic volume profiles don't align when using the psi_norm coordinate. Is this why the betaN values are different?

[orso82]

This is definitely pointing to a bug. I am thinking perhaps something in the IMAS dynamic expressions...

[bclyons12]

I haven't dug into this, but some general thoughts of where pitfalls might be:

1. The equilibrium IDS does the silly thing where it uses psi, not psi_norm as the radial coordinate. psi gets updated after the equilibrium step, so we need to be careful that things are properly mapped using psi_norm.
2. After the equilibrium update, rho_tor_norm will change as a function of psi and psi_norm. That means in the core_profiles IDS, we need to make sure that psi is being updated appropriately.

My guess is that some combination of those two leads to the inconsistency in the coordinates.

[orso82]

Indeed this was an issue with how we handle one-time expressions: https://github.com/ProjectTorreyPines/IMAS.jl/blob/master/src/expressions/onetime.jl

These expressions are frozen the first time they are accessed. This is necessary to ensure that core_profiles, core_sources, and core_transport grids do not change after changing the equilibrium. The idea is that we want to freeze in the core_profiles, core_sources, and core_transport grids the rho, psi, volume, area, ... info that were used when those IDSs were filled. While this is generally ok, this is not desirable when iterating the equilibrium solver with other actors. In this case, at each iteration we want core_profiles, core_sources, and core_transport to take the grids from the latest equilibrium. https://github.com/ProjectTorreyPines/FUSE.jl/blob/c6b74326cd89a93ac80485cce36134cc1b1ff288/src/actors/compound/equilibrium_transport_actor.jl#L89

It seems to work ok now: