Open eldond opened 1 month ago
I think
https://github.com/ProjectTorreyPines/IMAS.jl/blob/b367a8ea1606f120c7c190aa278cd816cba4123d/src/physics/thermal_loads.jl#L501
should be if psi_sign == 1
...
ping @ggdose
@eldond could you please provide me with a dd to test the case?
@ggdose please see examples/wall_heat_flux_demo.jl
in https://github.com/ProjectTorreyPines/SynthDiag.jl/pull/35 .
It works with the proposed change of psi_sign == -1
to psi_sign == 1
@ggdose please see
examples/wall_heat_flux_demo.jl
in ProjectTorreyPines/SynthDiag.jl#35 . It works with the proposed change ofpsi_sign == -1
topsi_sign == 1
BTW it's necessary to use dvc to pull our samples. See instructions in SOLPSTestSamples
I failed in getting your dd. https://github.com/ProjectTorreyPines/SOLPSTestSamples leads me to a non-existent folder.
From the code alone, the bug seems to be that psi_sign (which should tell you whether psi in the sol is decreasing or increasing) is computed always after sorting the psi_level vector. This means that psi_sign is always +1. this is obviously wrong.
For decreasing psi (like the case you described: [0, ..., -blah, -0.9, -1.]) psi_sign must be -1. The function find_levels_from_P always orders psi_levels sep to wall ([1] is sep, [end] is wall).
Therefore the fix of the bug would be to move the line 499 https://github.com/ProjectTorreyPines/IMAS.jl/blob/db1e833e830024a67d12b90b10d90868adab2066/src/physics/thermal_loads.jl#L499 to line 497 (or before line 498, but after line 495) https://github.com/ProjectTorreyPines/IMAS.jl/blob/db1e833e830024a67d12b90b10d90868adab2066/src/physics/thermal_loads.jl#L490-L511
Please let me know if this fixes your issue. I apologize that I could not be able to load your case
Our samples are too big to load in git; we load instructions for pulling them with dvc into git. Git checks out the dvc instructions, then dvc pulls the actual files.
made a PR with the fix
Levels are being chosen in thermal_loads.jl in
mesher_HF()
https://github.com/ProjectTorreyPines/IMAS.jl/blob/db1e833e830024a67d12b90b10d90868adab2066/src/physics/thermal_loads.jl#L490-L511 The key I think is that the levels are sorted and possibly reversed. In my case, where psi_axis>0 and psi_boundary=0, the chosen levels start <0 (out in the SOL) and count up to 0 (the LCFS). However, line 509 assigns the LCFS flux value to element 1, whereas the separatrix is the end of the array. When this set of levels is used later in sol.jlsol()
, it trips a "levels_is_not_monotonic_in_Ip_direction
" exception https://github.com/ProjectTorreyPines/IMAS.jl/blob/db1e833e830024a67d12b90b10d90868adab2066/src/physics/sol.jl#L190-L195 because the array goes something like[0, -1., -0.9, -blah, ..., 0]
(not those exact values but you get the idea).Either the levels array is going in the wrong direction or the wrong element is being replaced in thermal_loads.jl, or my geqdsk is messed up in terms of Ip and sign of my psi flux map. My geqdsk has positive current and psi_axis > psi_boundary. It's not a real empirical EFIT result but a test/design equilibrium for ITER, so maybe it's configured weirdly and should've been sanitized while loading.