Changes in geqdsk.read() from 0.7.0 and 0.8.0

[jhillairet]

We've found an issue in reading GEQDSK files using geqdsk.read() between the versions 0.7.0 (or 0.6.1) and 0.8. With the same file and machine description, version 0.7.0 converges properly, while 0.8 not.

Before giving a minimal example, is there something evident we should be aware of between these two versions?

[arthuradriaens-code]

In the releases "What's changed" it says: "Switch to external FreeQDSK for *qdsk reading/writing by @ZedThree in #85" I guess the switch to FreeQDSK might be giving problems

[ZedThree]

@arthuradriaens-code That's my immediate worry too!

@jhillairet Could you also try version 0.7 and see whether that works or not?

[jhillairet]

@jhillairet Could you also try version 0.7 and see whether that works or not?

@ZedThree it also works with 0.7.0. I'll update the issue title.

[ZedThree]

Are you able to upload your eqdsk file? If not, could you see if there are any differences in the output of geqdsk.read() between the two versions?

[jhillairet]

Sure. Attached is the GEQDSK file and below the output of geqdsk_freegs.read(). I use rtol=0.1 to make it faster, otherwise it tries to converge during many minutes.

Python 3.10 + freegs 0.7:

  nx = 65, ny = 65
370
CURRENT:  -487377.89145023254
Plasma current: -487377.89145023254 Amps, input: -499000.0 Amps

==========================
A : Coil(R=0.738, Z=0, current=-3753.3, turns=195, control=True)
Bh : Coil(R=1.118, Z=1.75, current=-10836.9, turns=176, control=True)
Dh : Coil(R=2.885, Z=1.92, current=1158.7, turns=95, control=True)
Eh : Coil(R=3.77, Z=1.523, current=373.7, turns=96, control=True)
Fh : Coil(R=4.375, Z=0.637, current=1077.3, turns=96, control=True)
Fb : Coil(R=4.375, Z=-0.637, current=27.1, turns=96, control=True)
Eb : Coil(R=3.77, Z=-1.523, current=2167.0, turns=96, control=True)
Db : Coil(R=2.885, Z=-1.92, current=-179.0, turns=95, control=True)
Bb : Coil(R=1.118, Z=-1.75, current=-6876.8, turns=176, control=True)
divh : Circuit([("Xh1", Coil(R=2.0115, Z=0.7432, current=-5536.6, turns=4, control=True), 1.0), ("Xh2", Coil(R=2.0755, Z=0.7772, current=-5536.6, turns=4, control=True), 1.0), ("Xh3", Coil(R=2.1665, Z=0.8058, current=-5536.6, turns=4, control=True), 1.0), ("Xh4", Coil(R=2.2305, Z=0.8398, current=-5536.6, turns=4, control=True), 1.0)], current=-5536.629971103787, control=True)
divb : Circuit([("Xb1", Coil(R=2.0115, Z=-0.7403, current=-11115.4, turns=4, control=True), 1.0), ("Xb2", Coil(R=2.0755, Z=-0.7743, current=-11115.4, turns=4, control=True), 1.0), ("Xb3", Coil(R=2.1665, Z=-0.8029, current=-11115.4, turns=4, control=True), 1.0), ("Xb4", Coil(R=2.2305, Z=-0.8369, current=-11115.4, turns=4, control=True), 1.0)], current=-11115.363358316634, control=True)
==========================
Plasma current: -487470.1155736676 Amps, input: -499000.0 Amps
Plasma pressure on axis: 37545.6775 Pascals
==========================
A : Coil(R=0.738, Z=0, current=-2643.2, turns=195, control=True)
Bh : Coil(R=1.118, Z=1.75, current=-9650.8, turns=176, control=True)
Dh : Coil(R=2.885, Z=1.92, current=1051.1, turns=95, control=True)
Eh : Coil(R=3.77, Z=1.523, current=519.6, turns=96, control=True)
Fh : Coil(R=4.375, Z=0.637, current=1052.2, turns=96, control=True)
Fb : Coil(R=4.375, Z=-0.637, current=46.9, turns=96, control=True)
Eb : Coil(R=3.77, Z=-1.523, current=2271.6, turns=96, control=True)
Db : Coil(R=2.885, Z=-1.92, current=-266.5, turns=95, control=True)
Bb : Coil(R=1.118, Z=-1.75, current=-5692.8, turns=176, control=True)
divh : Circuit([("Xh1", Coil(R=2.0115, Z=0.7432, current=-5628.8, turns=4, control=True), 1.0), ("Xh2", Coil(R=2.0755, Z=0.7772, current=-5628.8, turns=4, control=True), 1.0), ("Xh3", Coil(R=2.1665, Z=0.8058, current=-5628.8, turns=4, control=True), 1.0), ("Xh4", Coil(R=2.2305, Z=0.8398, current=-5628.8, turns=4, control=True), 1.0)], current=-5628.816824346162, control=True)
divb : Circuit([("Xb1", Coil(R=2.0115, Z=-0.7403, current=-11205.0, turns=4, control=True), 1.0), ("Xb2", Coil(R=2.0755, Z=-0.7743, current=-11205.0, turns=4, control=True), 1.0), ("Xb3", Coil(R=2.1665, Z=-0.8029, current=-11205.0, turns=4, control=True), 1.0), ("Xb4", Coil(R=2.2305, Z=-0.8369, current=-11205.0, turns=4, control=True), 1.0)], current=-11205.043114638209, control=True)
==========================

Python 3.12 + freegs 0.8:

CURRENT:  -471226.0535814483
Plasma current: -471226.0535814483 Amps, input: -499000.0 Amps

==========================
A : Coil(R=0.738, Z=0, current=-235309.2, turns=195, control=True)
Bh : Coil(R=1.118, Z=1.75, current=-243582.1, turns=176, control=True)
Dh : Coil(R=2.885, Z=1.92, current=18906.4, turns=95, control=True)
Eh : Coil(R=3.77, Z=1.523, current=-28088.1, turns=96, control=True)
Fh : Coil(R=4.375, Z=0.637, current=2685.2, turns=96, control=True)
Fb : Coil(R=4.375, Z=-0.637, current=5089.3, turns=96, control=True)
Eb : Coil(R=3.77, Z=-1.523, current=-29567.6, turns=96, control=True)
Db : Coil(R=2.885, Z=-1.92, current=20142.1, turns=95, control=True)
Bb : Coil(R=1.118, Z=-1.75, current=-247293.7, turns=176, control=True)
divh : Circuit([("Xh1", Coil(R=2.0115, Z=0.7432, current=15920.9, turns=4, control=True), 1.0), ("Xh2", Coil(R=2.0755, Z=0.7772, current=15920.9, turns=4, control=True), 1.0), ("Xh3", Coil(R=2.1665, Z=0.8058, current=15920.9, turns=4, control=True), 1.0), ("Xh4", Coil(R=2.2305, Z=0.8398, current=15920.9, turns=4, control=True), 1.0)], current=15920.898761407567, control=True)
divb : Circuit([("Xb1", Coil(R=2.0115, Z=-0.7403, current=18318.3, turns=4, control=True), 1.0), ("Xb2", Coil(R=2.0755, Z=-0.7743, current=18318.3, turns=4, control=True), 1.0), ("Xb3", Coil(R=2.1665, Z=-0.8029, current=18318.3, turns=4, control=True), 1.0), ("Xb4", Coil(R=2.2305, Z=-0.8369, current=18318.3, turns=4, control=True), 1.0)], current=18318.298365373674, control=True)
==========================
psi_relchange: 0.19061389112449315
bndry_relchange: 1.0
bndry_change: -16.859387105599563


psi_relchange: 0.32051298791337046
bndry_relchange: 0.04267078169742521
bndry_change: -0.751468996225757


psi_relchange: 0.2983775159686134
bndry_relchange: 0.03817278379428058
bndry_change: -0.6989357247127224


psi_relchange: 0.133005691676382
bndry_relchange: 0.016061816118604674
bndry_change: -0.29888921306791616


psi_relchange: 0.15832647503865838
bndry_relchange: 0.004518759536850448
bndry_change: -0.08446985387369566


psi_relchange: 0.15960668669158853
bndry_relchange: 0.009803900237463513
bndry_change: -0.1850802952339201


psi_relchange: 0.18743843793976184
bndry_relchange: 0.030239223818907826
bndry_change: -0.5886637944552149


psi_relchange: 0.20253618482208216
bndry_relchange: 0.008577584619557899
bndry_change: -0.16842360673689427


psi_relchange: 0.1484529785446875
bndry_relchange: 0.008411016420651408
bndry_change: -0.16655387445740288


psi_relchange: 0.062140961302465623
bndry_relchange: 0.00209731544129227
bndry_change: -0.04161805908396232


Plasma current: -893916.1695133409 Amps, input: -499000.0 Amps
Plasma pressure on axis: 37545.6775 Pascals
==========================
A : Coil(R=0.738, Z=0, current=463126.0, turns=195, control=True)
Bh : Coil(R=1.118, Z=1.75, current=397286.5, turns=176, control=True)
Dh : Coil(R=2.885, Z=1.92, current=3199.2, turns=95, control=True)
Eh : Coil(R=3.77, Z=1.523, current=9775.9, turns=96, control=True)
Fh : Coil(R=4.375, Z=0.637, current=20278.2, turns=96, control=True)
Fb : Coil(R=4.375, Z=-0.637, current=13876.6, turns=96, control=True)
Eb : Coil(R=3.77, Z=-1.523, current=12072.8, turns=96, control=True)
Db : Coil(R=2.885, Z=-1.92, current=4507.5, turns=95, control=True)
Bb : Coil(R=1.118, Z=-1.75, current=393851.6, turns=176, control=True)
divh : Circuit([("Xh1", Coil(R=2.0115, Z=0.7432, current=-14900.7, turns=4, control=True), 1.0), ("Xh2", Coil(R=2.0755, Z=0.7772, current=-14900.7, turns=4, control=True), 1.0), ("Xh3", Coil(R=2.1665, Z=0.8058, current=-14900.7, turns=4, control=True), 1.0), ("Xh4", Coil(R=2.2305, Z=0.8398, current=-14900.7, turns=4, control=True), 1.0)], current=-14900.705543363762, control=True)
divb : Circuit([("Xb1", Coil(R=2.0115, Z=-0.7403, current=-13026.7, turns=4, control=True), 1.0), ("Xb2", Coil(R=2.0755, Z=-0.7743, current=-13026.7, turns=4, control=True), 1.0), ("Xb3", Coil(R=2.1665, Z=-0.8029, current=-13026.7, turns=4, control=True), 1.0), ("Xb4", Coil(R=2.2305, Z=-0.8369, current=-13026.7, turns=4, control=True), 1.0)], current=-13026.680592987348, control=True)
==========================

plasma current

WEST56898_time6_Equilibrium.zip

[ZedThree]

Thanks! Sorry, I had forgotten there was an extra layer under geqdsk.read() for just reading the geqsk file. I can confirm that that is identical between 0.7.0 and 0.8.0 (using freeqdsk) at least, so the difference is due to something else.

It looks like there were only three changes to geqdsk.read() since 0.7.0 (excluding the file reader and some formatting changes), so it's likely one of these is the culprit:

• 9e1c6ef Constrained geqdsk reading working
• 734265d Coil current constraints in geqdsk reader
• 0eba3d8 FIxes to geqdsk reader for limited plasmas