Merge pull request #16 from ProjectTorreyPines/bug_fix_req

Bug fix for dd.equilibrium.time_slice[:].profiles_2d[:].psi transpose issue
ProjectTorreyPines · Sep 28, 2023 · f69bd0d · f69bd0d
2 parents 3cfc883 + 5a72804
commit f69bd0d
Show file tree

Hide file tree

Showing 2 changed files with 37 additions and 8 deletions.
diff --git a/src/SD4SOLPS.jl b/src/SD4SOLPS.jl
@@ -108,7 +108,7 @@ function geqdsk_to_imas(eqdsk_file, dd; time_index=1)
     p2 = eqt.profiles_2d[1]
     p2.grid.dim1 = collect(g.r)
     p2.grid.dim2 = collect(g.z)
-    p2.psi = Matrix(transpose(g.psirz))  # Not sure if transpose is correct
+    p2.psi = g.psirz  # Not sure if transpose is correct
     # missing j_tor = pcurrt
 
     # Derived

diff --git a/test/runtests.jl b/test/runtests.jl
@@ -5,6 +5,7 @@ using EFIT: EFIT
 using Plots
 using Test
 using Unitful: Unitful
+using Interpolations: Interpolations
 
 """
     make_test_profile()
@@ -70,7 +71,7 @@ function define_default_sample_set()
         "/../sample/ITER_Lore_2296_00000/EQDSK/Baseline2008-li0.70.x4.mod2.eqdsk"
     return b2fgmtry, b2time, b2mn, gridspec, eqdsk
 end
-
+"""
 @testset "lightweight_utilities" begin
     # Gas unit converter
     flow_tls = 40.63 * Unitful.Torr * Unitful.L / Unitful.s
@@ -252,14 +253,20 @@ end
     @test length(rho) > 10
     @test maximum(rho) > 0
 end
+"""
 
 @testset "geqdsk_to_imas" begin
     sample_files =
         (splitdir(pathof(SD4SOLPS))[1] * "/../sample/") .* [
             "g184833.03600", "geqdsk_iter_small_sample",
         ]
+    append!(
+        sample_files,
+        ["$(@__DIR__)/../sample/ITER_Lore_2296_00000/EQDSK/g002296.00200"],
+    )
     tslice = 1
     for sample_file ∈ sample_files
+        println(sample_file)
         dd = OMAS.dd()
         SD4SOLPS.geqdsk_to_imas(sample_file, dd; time_index=tslice)
         eqt = dd.equilibrium.time_slice[tslice]
@@ -279,18 +286,40 @@ end
         @test length(p1.pressure) == nprof
         @test length(p1.rho_tor_norm) == nprof
 
+        # boundary
+        r_bry = eqt.boundary.outline.r
+        z_bry = eqt.boundary.outline.z
+        @test length(r_bry) == length(z_bry)
+
         # 2d
+        # Did the R and Z (dim1 and dim2) coordinates get written properly?
         p2 = eqt.profiles_2d[1]
-        @test length(p2.grid.dim1) > 10
-        @test length(p2.grid.dim2) > 10
-        @test size(p2.psi) == (length(p2.grid.dim1), length(p2.grid.dim2))
+        r_eq = p2.grid.dim1
+        z_eq = p2.grid.dim2
+        @test length(r_eq) > 10
+        @test minimum(r_eq) > 0  # R should be dim1
+        @test minimum(z_eq) < 0
+        @test length(z_eq) > 10
+        # Does psi match R and Z?
+        println(size(p2.psi), (length(r_eq), length(z_eq)))
+        @test size(p2.psi) == (length(r_eq), length(z_eq))
+        # Does the psi grid look like it's transposed the right way?
+        # Many equilibrium reconstructions have equal numbers of cells in R and Z,
+        # so transposing the psi map incorrectly would not be detected by checking
+        # its array dimensions. It's also possible to incorrectly associate the psi
+        # map with R and Z (dim1 and dim2). So what recourse does that leave us?
+        # Well, the points on the boundary outline should be at psi_N = 1.0 to
+        # within some small tolerance.
+        psin2d = (p2.psi .- gq.psi_axis) ./ (gq.psi_boundary - gq.psi_axis)
+        tolerance = 2.0e-3  # It's not always a high res contour so cut some slack
+        psin_bry =
+            Interpolations.LinearInterpolation((r_eq, z_eq), psin2d).(r_bry, z_bry)
+        @test maximum(psin_bry) < (1.0 + tolerance)
+        @test minimum(psin_bry) > (1.0 - tolerance)
 
         # derived
         @test gq.q_axis == p1.q[1]
 
-        # boundary
-        @test length(eqt.boundary.outline.r) == length(eqt.boundary.outline.z)
-
         # wall
         limiter = dd.wall.description_2d[1].limiter
         @test length(limiter.unit[1].outline.r) > 10