Releasing v.0.2.0

Project.toml

@@ -1,18 +1,20 @@
 name = "SD4SOLPS"
 uuid = "d8db6f1b-e564-4c04-bba3-ef399f78c070"
 authors = ["David Eldon <[email protected]>"]
-version = "0.1.0"
+version = "0.2.0"
 
 [deps]
 ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
 Contour = "d38c429a-6771-53c6-b99e-75d170b6e991"
 EFIT = "cda752c5-6b03-55a3-9e33-132a441b0c17"
 GGDUtils = "b7b5e640-9b39-4803-84eb-376048795def"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
+JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 OMAS = "91cfaa06-6526-4804-8666-b540b3feef2f"
 PhysicalConstants = "5ad8b20f-a522-5ce9-bfc9-ddf1d5bda6ab"
 PlotUtils = "995b91a9-d308-5afd-9ec6-746e21dbc043"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+PolygonOps = "647866c9-e3ac-4575-94e7-e3d426903924"
 SOLPS2IMAS = "09becab6-0636-4c23-a92a-2b3723265c31"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"

apply_format.sh

@@ -0,0 +1,2 @@
+#!/bin/bash
+julia -e 'using JuliaFormatter; format(".", verbose=true)'

data/.gitignore

@@ -0,0 +1,2 @@
+*.mesh_ext.json
+

src/SD4SOLPS.jl

@@ -7,7 +7,7 @@ using Interpolations: Interpolations
 #import GGDUtils
 
 export find_files_in_allowed_folders
-export geqdsk_to_imas
+export geqdsk_to_imas!
 
 include("$(@__DIR__)/supersize_profile.jl")
 include("$(@__DIR__)/repair_eq.jl")
@@ -63,12 +63,12 @@ function find_files_in_allowed_folders(
 end
 
 """
-    geqdsk_to_imas()
+    geqdsk_to_imas!()
 
 Transfers the equilibrium reconstruction in an EFIT-style gEQDSK file into
 the IMAS DD structure.
 """
-function geqdsk_to_imas(eqdsk_file, dd; time_index=1)
+function geqdsk_to_imas!(eqdsk_file, dd; time_index=1)
     # https://github.com/JuliaFusion/EFIT.jl/blob/master/src/io.jl
     g = EFIT.readg(eqdsk_file)
     # Copying ideas from OMFIT: omfit/omfit_classes/omfit_eqdsk.py / to_omas()
@@ -122,6 +122,39 @@ function geqdsk_to_imas(eqdsk_file, dd; time_index=1)
         gq.q_min.rho_tor_norm = g.rhovn[qmin_idx]
     end
 
+    # X-points
+    xrs, xzs, xpsins, xseps = EFIT.x_points(g; within_limiter_only=false)
+    if length(xrs) > 0
+        bx = eqt.boundary.x_point
+        resize!(bx, length(xrs))
+        for i ∈ length(xrs)
+            bx[i].r = xrs[i]
+            bx[i].z = xzs[i]
+        end
+        nprim = sum(xseps .== 1)
+        if nprim > 0
+            bsx = eqt.boundary_separatrix.x_point
+            resize!(bsx, nprim)
+            xrprim = xrs[xseps.==1]
+            xzprim = xzs[xseps.==1]
+            for i ∈ nprim
+                bsx[i].r = xrprim[i]
+                bsx[i].z = xzprim[i]
+            end
+        end
+        nsec = sum(xseps .== 2)
+        if nsec > 0
+            bssx = eqt.boundary_secondary_separatrix.x_point
+            resize!(bssx, nsec)
+            xrsec = xrs[xseps.==2]
+            xzsec = xzs[xseps.==2]
+            for i ∈ nsec
+                bssx[i].r = xrsec[i]
+                bssx[i].z = xzsec[i]
+            end
+        end
+    end
+
     # Boundary / LCFS
     eqt.boundary.outline.r = g.rbbbs
     eqt.boundary.outline.z = g.zbbbs
@@ -134,7 +167,8 @@ function geqdsk_to_imas(eqdsk_file, dd; time_index=1)
     limiter.type.description = "first wall"
     resize!(limiter.unit, 1)
     limiter.unit[1].outline.r = g.rlim
-    return limiter.unit[1].outline.z = g.zlim
+    limiter.unit[1].outline.z = g.zlim
+    return
 end
 
 """
@@ -238,14 +272,15 @@ function preparation(
     println("    eqdsk = ", eqdsk)
 
     dd = SOLPS2IMAS.solps2imas(b2fgmtry, b2time, gridspec, b2mn)
-    geqdsk_to_imas(eqdsk, dd)
+    geqdsk_to_imas!(eqdsk, dd)
     println("Loaded input data into IMAS DD")
 
     fill_in_extrapolated_core_profile!(dd, "electrons.density"; method=core_method)
     fill_in_extrapolated_core_profile!(dd, "electrons.temperature"; method=core_method)
     # ... more profiles here as they become available in b2time
     println("Extrapolated core profiles")
 
+    cached_mesh_extension!(dd, eqdsk_file, b2fgmtry)
     fill_in_extrapolated_edge_profile!(dd, "electrons.density"; method=core_method)
     # ... more profiles here
     println("Extrapolated edge profiles (but not really (placeholder only))")

src/repair_eq.jl

@@ -48,7 +48,7 @@ function add_rho_to_equilibrium!(dd::OMAS.dd)
         psi = eqt.profiles_1d.psi
         n = length(psi)
         if (length(eqt.profiles_1d.rho_tor_norm) > 0)
-            if max(eqt.profiles_1d.rho_tor_norm) > 0
+            if maximum(eqt.profiles_1d.rho_tor_norm) > 0
                 println("Slice #", it, " already has a rho_tor_norm profile; skipping")
                 continue
             end