Minor bug fix in writeg and parse_gfile_header

README.md

@@ -5,9 +5,9 @@
 [EFIT (Equilibrium Fitting)](https://fusion.gat.com/theory/Efit) is a computer code developed to translate measurements from plasma diagnostics into useful information like plasma geometry, stored energy, and current profiles.
 
 
-EFIT.jl provides basic functionality for reading [EFIT GEQDSK](https://fusion.gat.com/theory/Efitgeqdsk) files.
+`EFIT.jl` provides basic functionality for reading [EFIT GEQDSK](https://fusion.gat.com/theory/Efitgeqdsk) files.
 
-```
+```julia-repl
 julia> using EFIT
 
 julia> g = readg(EFIT.test_gfile);
@@ -41,3 +41,24 @@ julia> minor_radius(g)
 julia> aspect_ratio(g)
 2.6031891587258285
 ```
+
+
+`EFIT.jl` can convert `dd` into GEQDSK files as in the following example:
+```julia
+import EFIT
+import EFIT.IMASdd
+
+# Load example dd
+filename = joinpath(dirname(dirname(pathof(EFIT))), "test", "test_dd_eq.json")
+dd = IMASdd.json2imas(filename);
+
+# First, convert equilibrium time_slices in dd to a list of `GEQDSK` struct
+ggs = EFIT.imas2geqdsk(dd)
+
+# Then, write into files
+for gg in ggs
+    # default file name has "g00000.xxxxx" format, where xxxxx is time in `ms`
+    file_name = gg.file
+    EFIT.writeg(gg, file_name)
+end
+```

src/geqdsk_imas.jl

@@ -475,22 +475,32 @@ function imas2geqdsk(
         zbbbs = eqt.boundary.outline.z ./ tc["Z"]
         nbbbs = length(rbbbs)
 
-        # 1D profiles
-        p1 = eqt.profiles_1d
-        psi = p1.psi ./ tc["PSI"]
-        qpsi = p1.q ./ tc["Q"]
-        pres = p1.pressure ./ tc["P"]
-        pprime = p1.dpressure_dpsi ./ tc["PPRIME"]
-        fpol = p1.f ./ tc["F"]
-        ffprim = p1.f_df_dpsi ./ tc["F_FPRIME"]
-        rhovn = p1.rho_tor_norm
-
         # 2D flux map
-        p2 = eqt.profiles_2d[1]
+        p2 = findfirst(:rectangular, eqt.profiles_2d)
         r = p2.grid.dim1 ./ tc["R"]
         z = p2.grid.dim2 ./ tc["Z"]
         psirz = p2.psi ./ tc["PSI"]
 
+        # 1D profiles
+        function interpolate_1d_profile(ori_1D::AbstractVector{<:Real}, N::Int)
+            if length(ori_1D) == N
+                return ori_1D
+            else
+                itp=IMASdd.interp1d(p1.psi_norm, ori_1D, :cubic)
+                return itp.(range(0,1,N))
+            end
+        end
+
+        p1 = eqt.profiles_1d
+        nw = length(r)
+        psi = interpolate_1d_profile(p1.psi ./ tc["PSI"], nw)
+        qpsi = interpolate_1d_profile(p1.q ./ tc["Q"], nw)
+        pres = interpolate_1d_profile(p1.pressure ./ tc["P"], nw)
+        pprime = interpolate_1d_profile(p1.dpressure_dpsi ./ tc["PPRIME"], nw)
+        fpol = interpolate_1d_profile(p1.f ./ tc["F"], nw)
+        ffprim = interpolate_1d_profile(p1.f_df_dpsi ./ tc["F_FPRIME"], nw)
+        rhovn = interpolate_1d_profile(p1.rho_tor_norm, nw)
+
         bcentr = IMASdd.@ddtime(dd.equilibrium.vacuum_toroidal_field.b0) ./ tc["B"]
         time = eqt.time
         @debug "eqt.time = $(eqt.time), time out = $time"

src/io.jl

@@ -114,6 +114,9 @@ function parse_gfile_header(headerline::String; set_time=nothing)
         try
             time /= 1000.0  # Assume it was written in ms and convert to s
         catch e
+            if time isa Expr
+                time = string(time)
+            end
             if time isa String
                 rg2 = r"([[:digit:]]+|(?:[[:punct:]]|[[:blank:]])+)"
                 try
@@ -273,24 +276,25 @@ function writeg(g::GEQDSKFile, filename::String;
         @warn("Warning: A file with the name '$filename' already exists. It will be overwritten.")
     end
 
+    # Note: GEQDSKFile.time is converted to ms (to follow readg's notation)
+    # Remove all newline characters in description
+    clean_desc = replace(desc, "\n" => "")
+    description = join([clean_desc, shot, time], "   ")
+    if length(description) > 48
+        error("""
+            Description too long (length = $(length(description)) > max 48).
+            Current description: '$description'
+            Please shorten one or more of the following kwargs:
+                desc = "$desc"
+                shot = "$shot"
+                time = "$time"
+            """)
+    end
+
     try
         # Open the target file for writing
         open(filename, "w") do f
             # Write header
-            # Note: GEQDSKFile.time is converted to ms (to follow readg's notation)
-            # Remove all newline characters in description
-            clean_desc = replace(desc, "\n" => "")
-            description = join([clean_desc, shot, time], "   ")
-
-            if length(description) > 48
-                error("Description too long (length = $(length(description)) > max 48).\n" *
-                      "Current description: '$description'\n" *
-                      "Please shorten one or more of the following kwargs:\n" *
-                      "  desc = \"$desc\"\n" *
-                      "  shot = \"$shot\"\n" *
-                      "  time = \"$time\"")
-            end
-
             @printf(f,"%-48s%4d%4d%4d\n",description, 0, g.nw, g.nh)
 
             @printf(f,"%16.9E%16.9E%16.9E%16.9E%16.9E\n", g.rdim, g.zdim, g.rcentr, g.rleft, g.zmid)
@@ -318,7 +322,8 @@ function writeg(g::GEQDSKFile, filename::String;
         @info "Successfully wrote GEQDSKFile to '$filename'."
         return true
     catch e
-        @error "An error occurred while writing to file '$filename': $e"
+        @error """An error occurred while writing to file '$filename':
+             $(sprint(showerror, e))"""
         return false
     end
 end

test/runtests.jl

@@ -1,8 +1,8 @@
 using Test
 using EFIT
 import EFIT
-import IMASdd
-import CoordinateConventions
+import EFIT.IMASdd
+import EFIT.CoordinateConventions
 
 g = readg(EFIT.test_gfile)
 g2 = readg(EFIT.test_gfile2)
@@ -117,6 +117,10 @@ end
         @test new_g == ori_g
         @test isequal(new_g, ori_g)
 
+        # keywords test
+        @test writeg(new_g, tmp_filename; desc="my description", shot="my shot 12345", time="1000ms")
+        @test_throws Exception writeg(ori_g, tmp_filename; desc="This is a long description that will cause an error.")
+
         # Delete temporary file
         rm(tmp_filename, force=true)
     end
@@ -174,8 +178,14 @@ end
     # Add one geqdsk
     println("test file 1 info: filename = $(gs[1].file), time = $(gs[1].time)")
     dd = IMASdd.dd()
-    eqt = resize!(dd.equilibrium.time_slice, length(gs))[1]
     dd.equilibrium.time = [g.time for g in gs]
+
+    resize!(dd.equilibrium.time_slice, length(gs))
+    for (k, eqt) in pairs(dd.equilibrium.time_slice)
+        eqt.time = gs[k].time
+    end
+
+    eqt = dd.equilibrium.time_slice[1]
     EFIT.geqdsk2imas!(gs[1], eqt; wall=dd.wall, cocos_clockwise_phi=cocos_clockwise_phi, add_derived=true)
     @test length(eqt.profiles_2d[1].grid.dim1) > 1
     @test eqt.time == gs[1].time
@@ -256,3 +266,55 @@ end
     gg3.file = gs[2].file
     gcompare(gg3, gs[2])
 end
+
+@testset "imas2geqdsk & geqdsk2imas!" begin
+    ori_dd = IMASdd.json2imas("test_dd_eq.json")
+
+    gg=EFIT.imas2geqdsk(ori_dd)
+
+    new_dd = IMASdd.dd()
+    EFIT.geqdsk2imas!(gg, new_dd)
+
+    # compare wall
+    @test new_dd.wall == ori_dd.wall
+    @test new_dd.equilibrium.vacuum_toroidal_field == ori_dd.equilibrium.vacuum_toroidal_field
+
+    ori_eqt = ori_dd.equilibrium.time_slice[]
+    new_eqt = new_dd.equilibrium.time_slice[]
+
+    # compare boundary
+    @test new_eqt.boundary.geometric_axis == ori_eqt.boundary.geometric_axis
+    @test new_eqt.boundary.outline == ori_eqt.boundary.outline
+
+    # compare global_quantities
+    ori_gq = new_eqt.global_quantities
+    new_gq = new_eqt.global_quantities
+
+    @test ori_gq.ip == new_gq.ip
+    @test ori_gq.magnetic_axis == new_gq.magnetic_axis
+    @test ori_gq.psi_axis == new_gq.psi_axis
+    @test ori_gq.psi_boundary == new_gq.psi_boundary
+
+    # compare profiles_1d
+    ori_p1d = ori_eqt.profiles_1d
+    new_p1d = new_eqt.profiles_1d
+
+    target_fields = filter(x -> x ∉ IMASdd.private_fields, fieldnames(typeof(new_p1d)))
+    for field in target_fields
+        if !ismissing(new_p1d, field) && !isempty(new_p1d, field)
+            ori1D_vec = getproperty(ori_p1d, field)
+            new1D_vec = getproperty(new_p1d, field)
+
+            itp=IMASdd.interp1d(ori_p1d.psi_norm, ori1D_vec, :cubic)
+            @test isapprox(new1D_vec, itp.(range(0,1, length(new_p1d.psi))))
+        end
+    end
+
+    # compare profiles_2d
+    ori_p2d = findfirst(:rectangular, ori_eqt.profiles_2d)
+    new_p2d = findfirst(:rectangular, new_eqt.profiles_2d)
+
+    @test new_p2d.grid == ori_p2d.grid
+    @test isapprox(new_p2d.psi, ori_p2d.psi)
+
+end