Big differences from CHEASE after FUSE.init, and hangs on EquilibriumTransport

[daveweisberg]

On a typical FPP optimization run, the following case is an example of how TEQUILA and CHEASE produce very different results after the FUSE.init step. While CHEASE produces a reasonable initialized core solution and converges after 4 EquilibriumTransport iterations, TEQUILA produces a very high-temp initialized core solution and then hangs indefinitely during the EquilibriumTransport actor.

Executed code:

# load case and reduce plasma current
load_dir = "/fusion/ga/projects/ird/ptp/weisbergd/julia/dev/FUSE/playground/weisberg_FPP_opt/very_exploratory/opt_FPP_Solovev_Kr_HTS0.1_AspectRatio_fk0.88_fLH_ECaccess_qpol2.5/opt_runs/2023-06-05T16:27:22.052__2437505"
ini = FUSE.json2ini(joinpath(load_dir, "ini.json"))
act = FUSE.json2act(joinpath(load_dir, "act.json"))
dd0 = IMAS.json2imas(joinpath(load_dir, "dd.json"));
ini.equilibrium.ip = 8e6

# run with CHEASE
act.ActorEquilibriumTransport.do_plot = true
act.ActorEquilibrium.model = :CHEASE
act.ActorEquilibriumTransport.max_iter = 5
act.ActorCXbuild.rebuild_wall = true
act.ActorCHEASE.free_boundary = true
dd = IMAS.dd();
FUSE.init(dd, ini, act);
FUSE.digest(dd);

FUSE.ActorWholeFacility(dd, act);
FUSE.digest(dd);

[ Info: Iteration = 1 , convergence error = 0.486, threshold = 0.05
[ Info: Iteration = 2 , convergence error = 0.195, threshold = 0.05
[ Info: Iteration = 3 , convergence error = 0.055, threshold = 0.05
[ Info: Iteration = 4 , convergence error = 0.014, threshold = 0.05

# run with TEQUILA 
act.ActorEquilibriumTransport.do_plot = true
act.ActorEquilibrium.model = :TEQUILA
act.ActorEquilibriumTransport.max_iter = 5
act.ActorCXbuild.rebuild_wall = true
dd.ActorTEQUILA.free_boundary = true
dd = IMAS.dd();
FUSE.init(dd, ini, act);
FUSE.digest(dd);

FUSE.ActorWholeFacility(dd, act);
FUSE.digest(dd);

This does not produce any output, seemingly hanging before the 1st EquilibriumTransport iteration is completed.

[bclyons12]

This is likely related to #42 and was "fixed" by 95063fb. I suspect this case initializes a very high pressure profile and very high current-density profile. This needs to be scaled down to match the target Ip.

CHEASE appears to do a Bateman scaling to keep the flux surfaces unchanged, so it scales the pressure down too. This results in less Shafranov shift and easier convergence, but the pressure is not what is input to the equilibrium actor.

TEQUILA only scales the current profile, maintaining whatever pressure profile is requested. This is huge compared to the target Ip, so the Shafranov shift is very large. TEQUILA struggles to do a full Picard iteration with such a high Shafranov shift, so I introduced an optional relaxation parameter with 95063fb. I don't think this has been coupled to FUSE yet.

I'll double check that that's what's happening in this case, but if so, it's not an issue for TEQUILA, but for FUSE in my opinion. I'll leave the issue open until I confirm though.

[daveweisberg]

So the solution is to reduce ini.equilibrium.pressure_core? I've tried this, and it doesn't seem to help. In fact, the resulting P0 produced by FUSE.init using TEQUILA is 2.6x the ini value.

Note that the CHEASE solution above converges to P0 = 1.29MPa, so I'll reduce ini.equilibrium.pressure_core to that value before running FUSE.init :

load_dir = "/fusion/ga/projects/ird/ptp/weisbergd/julia/dev/FUSE/playground/weisberg_FPP_opt/very_exploratory/opt_FPP_Solovev_Kr_HTS0.1_AspectRatio_fk0.88_fLH_ECaccess_qpol2.5/opt_runs/2023-06-05T16:27:22.052__2437505"
ini = FUSE.json2ini(joinpath(load_dir, "ini.json"))
act = FUSE.json2act(joinpath(load_dir, "act.json")) 
ini.equilibrium.ip = 8e6
ini.equilibrium.pressure_core = 1.3e6

# run with TEQUILA
act.ActorEquilibriumTransport.do_plot = true
act.ActorEquilibrium.model = :TEQUILA
act.ActorEquilibriumTransport.max_iter = 5
act.ActorCXbuild.rebuild_wall = true
act.ActorCHEASE.free_boundary = true
dd = IMAS.dd();
FUSE.init(dd, ini, act);
FUSE.digest(dd);

[daveweisberg]

One more update: If I further reduce ini.equilibrium.pressure_core, then FUSE.init gives a reasonable core solution. However, EquilibriumTransport still hangs with TEQUILA when it runs with this initial input. Why?

load_dir = "/fusion/ga/projects/ird/ptp/weisbergd/julia/dev/FUSE/playground/weisberg_FPP_opt/very_exploratory/opt_FPP_Solovev_Kr_HTS0.1_AspectRatio_fk0.88_fLH_ECaccess_qpol2.5/opt_runs/2023-06-05T16:27:22.052__2437505"
ini = FUSE.json2ini(joinpath(load_dir, "ini.json"))
act = FUSE.json2act(joinpath(load_dir, "act.json")) 
ini.equilibrium.ip = 8e6
ini.equilibrium.pressure_core = 0.6e6

# run with TEQUILA
act.ActorEquilibriumTransport.do_plot = true
act.ActorEquilibrium.model = :TEQUILA
act.ActorEquilibriumTransport.max_iter = 5
act.ActorCXbuild.rebuild_wall = true
act.ActorCHEASE.free_boundary = true
dd = IMAS.dd();
FUSE.init(dd, ini, act);
FUSE.digest(dd);

FUSE.ActorEquilibriumTransport(dd, act);
FUSE.digest(dd);

which still hangs indefinitely...

[bclyons12]

Trying to investigate, but I'm getting some runtime errors. @daveweisberg are you running on any special branches?

act = FUSE.json2act(joinpath(load_dir, "act.json"))

runs but reports

Error: reading act.ActorDivertors.impurities : TypeError(:setfield!, "", Union{Missing, Vector{Symbol}}, Float64[])
└ @ SimulationParameters ~/.julia/dev/SimulationParameters/src/utils.jl:92

Then

FUSE.init(dd, ini, act);

errors with

[ Info: Equilibrium
[ Info:  CHEASE
[ Info: HCD
[ Info:  ECsimple
[ Info:  ICsimple
[ Info:  LHsimple
[ Info:  NBsimple
[ Info: SteadyStateCurrent
MethodError: no method matching j_ohmic_steady_state!(::IMASDD.equilibrium__time_slice{Float64}, ::IMASDD.core_profiles{Float64}, ::Bool)

Closest candidates are:
  j_ohmic_steady_state!(::IMASDD.equilibrium__time_slice, ::IMASDD.core_profiles__profiles_1d)
   @ IMAS ~/.julia/dev/IMAS/src/physics/currents.jl:35


Stacktrace:
  [1] _step(actor::FUSE.ActorSteadyStateCurrent{Float64, Float64})
    @ FUSE ~/.julia/dev/FUSE/src/actors/current/steadycurrent_actor.jl:42
  [2] macro expansion
    @ ~/.julia/dev/FUSE/src/actors/abstract_actors.jl:37 [inlined]
  [3] macro expansion
    @ ~/.julia/packages/TimerOutputs/RsWnF/src/TimerOutput.jl:237 [inlined]
  [4] step(::FUSE.ActorSteadyStateCurrent{Float64, Float64}; kw::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ FUSE ~/.julia/dev/FUSE/src/actors/abstract_actors.jl:34
  [5] step(::FUSE.ActorSteadyStateCurrent{Float64, Float64})
    @ FUSE ~/.julia/dev/FUSE/src/actors/abstract_actors.jl:29
  [6] FUSE.ActorSteadyStateCurrent(dd::IMASDD.dd{Float64}, act::FUSE.ParametersActors{Float64}; kw::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ FUSE ~/.julia/dev/FUSE/src/actors/current/steadycurrent_actor.jl:32
  [7] FUSE.ActorSteadyStateCurrent(dd::IMASDD.dd{Float64}, act::FUSE.ParametersActors{Float64})
    @ FUSE ~/.julia/dev/FUSE/src/actors/current/steadycurrent_actor.jl:30
  [8] macro expansion
    @ ~/.julia/dev/FUSE/src/ddinit/init_currents.jl:12 [inlined]
  [9] macro expansion
    @ ~/.julia/packages/TimerOutputs/RsWnF/src/TimerOutput.jl:237 [inlined]
 [10] init_currents(dd::IMASDD.dd{Float64}, ini::FUSE.ParametersInits{Float64}, act::FUSE.ParametersActors{Float64})
    @ FUSE ~/.julia/dev/FUSE/src/ddinit/init_currents.jl:8
 [11] macro expansion
    @ ~/.julia/dev/FUSE/src/ddinit/init.jl:47 [inlined]
 [12] macro expansion
    @ ~/.julia/packages/TimerOutputs/RsWnF/src/TimerOutput.jl:237 [inlined]
 [13] init(dd::IMASDD.dd{Float64}, ini::FUSE.ParametersInits{Float64}, act::FUSE.ParametersActors{Float64}; do_plot::Bool)
    @ FUSE ~/.julia/dev/FUSE/src/ddinit/init.jl:12
 [14] init(dd::IMASDD.dd{Float64}, ini::FUSE.ParametersInits{Float64}, act::FUSE.ParametersActors{Float64})
    @ FUSE ~/.julia/dev/FUSE/src/ddinit/init.jl:10
 [15] top-level scope
    @ In[20]:1

[bclyons12]

I switched FUSE back to master and the second error disappears, so there's something inconsistent with this run case and the ProjectTorreyPines/FUSE.jl#348. Perhaps @TimSlendebroek knows why

[daveweisberg]

@bclyons12 I'm using the master branch for everything, including TEQUILA.

[bclyons12]

So the TEQUILA equilibrium from init looks fine to me. I tried

# run with TEQUILA 
act.ActorEquilibriumTransport.do_plot = true
act.ActorEquilibrium.model = :TEQUILA
act.ActorEquilibriumTransport.max_iter = 5
act.ActorCXbuild.rebuild_wall = true
act.ActorTEQUILA.free_boundary = true
ini.equilibrium.ip = 8e6
ini.equilibrium.pressure_core = 1.0e6
dd = IMAS.dd();
FUSE.init(dd, ini, act);
FUSE.digest(dd);

and get

so the pressure on-axis matches ini.equilibrium.pressure_core. The digests reports a P0 much higher though... I'm not sure why:

[bclyons12]

Without the core pressure setting, I run this:

load_dir = "/fusion/ga/projects/ird/ptp/weisbergd/julia/dev/FUSE/playground/weisberg_FPP_opt/very_exploratory/opt_FPP_Solovev_Kr_HTS0.1_AspectRatio_fk0.88_fLH_ECaccess_qpol2.5/opt_runs/2023-06-05T16:27:22.052__2437505"
ini = FUSE.json2ini(joinpath(load_dir, "ini.json"))
act = FUSE.json2act(joinpath(load_dir, "act.json"))
dd0 = IMAS.json2imas(joinpath(load_dir, "dd.json"));
ini.equilibrium.ip = 8e6

# run with TEQUILA 
act.ActorEquilibriumTransport.do_plot = true
act.ActorEquilibrium.model = :TEQUILA
act.ActorEquilibriumTransport.max_iter = 5
act.ActorCXbuild.rebuild_wall = true
act.ActorTEQUILA.free_boundary = true
dd = IMAS.dd();
FUSE.init(dd, ini, act);
FUSE.digest(dd);

which gives

Some things to note:

1. The core pressure on the equilibrium is ~1.5 MPa, but the digest report 4.28.
2. Note how over driven the current is. The noninductive current is ~14 MA, but the request is for 8 MA. Thus, you get this huge negative ohmic current. TEQUILA won't work if the current density changes sign since psi would be non-monotonic and rho_poloidal is a non-sensical parameter. Formerly TEQUILA would hang in such situations. Now I force it to check and throw an error:
3. The CHEASE run avoids this issue because it scales the pressure down, so the bootstrap current isn't so high and the ohmic current remains positive.

[orso82]

@bclyons12 the pressure in the digest is evaluated here:
https://github.com/ProjectTorreyPines/IMAS.jl/blob/525f9e3f5340aa4a1f4480b992a8a6c16e5f03fc/src/extract.jl#L60
which is simply:

P0 = dd.core_profiles.profiles_1d[].pressure[1] / 1E6

[bclyons12]

So there's a big inconsistency after init between the core_profiles and equilibrium pressures. I guess that makes sense before you do the EquilibriumTransport iterations, but it means that the reported P0 shouldn't be compared to the requested one at that step

[orso82]

The logic for setting the equilibrium pressure when initializing dd is here
https://github.com/ProjectTorreyPines/FUSE.jl/blob/master/src/ddinit/init_equilibrium.jl#L117-L136

You can see that P0 is provided externally from ini.equilibrium.pressure_core.
Instead of relying on an external input we could try to guess P0, maybe using some scaling law?

[bclyons12]

Right, I just think @daveweisberg was expecting the digest P0 to match ini.equilibrium.pressure_core, but it won't until the iterations converge.