From 9a8a1d5f1476755e00dc516e80a2d4eae3a64d55 Mon Sep 17 00:00:00 2001
From: TimSlendebroek <32385057+TimSlendebroek@users.noreply.github.com>
Date: Wed, 9 Mar 2022 11:34:08 -0800
Subject: [PATCH 01/16] add IMASDD to install and update

---
 Makefile | 20 ++++++++++++++++++--
 1 file changed, 18 insertions(+), 2 deletions(-)

diff --git a/Makefile b/Makefile
index 71d2a568a..4f6f0cac0 100644
--- a/Makefile
+++ b/Makefile
@@ -12,7 +12,7 @@ install: install_FUSE install_IJulia
 	julia -e '\
 using Pkg;\
 Pkg.activate();\
-Pkg.develop(["FUSE", "IMAS", "CoordinateConventions", "FusionMaterials", "AD_GS", "Equilibrium", "AD_TAUENN", "AD_EPEDNN", "AD_TGLFNN", "QED"]);\
+Pkg.develop(["FUSE", "IMAS","IMASDD", "CoordinateConventions", "FusionMaterials", "AD_GS", "Equilibrium", "AD_TAUENN", "AD_EPEDNN", "AD_TGLFNN", "QED"]);\
 Pkg.resolve();\
 try Pkg.upgrade_manifest() catch end;\
 '
@@ -24,7 +24,7 @@ Pkg.add("IJulia");\
 Pkg.build("IJulia");\
 '
 
-install_FUSE: install_IMAS install_FusionMaterials install_AD_GS install_Equilibrium install_TAUENN install_QED
+install_FUSE: install_IMAS install_IMASDD install_FusionMaterials install_AD_GS install_Equilibrium install_TAUENN install_QED
 	if [ ! -d "$(JULIA_PKG_DEVDIR)/FUSE" ]; then ln -s $(CURRENTDIR) $(JULIA_PKG_DEVDIR)/FUSE; fi
 	julia -e '\
 using Pkg;\
@@ -46,6 +46,19 @@ Pkg.resolve();\
 try Pkg.upgrade_manifest() catch end;\
 '
 
+install_IMASDD: install_CoordinateConventions
+	if [ ! -d "$(JULIA_PKG_DEVDIR)/IMASDD" ]; then\
+		julia -e 'using Pkg; Pkg.develop(url="git@github.com:ProjectTorreyPines/IMASDD.jl.git");';\
+	fi
+	julia -e '\
+using Pkg;\
+Pkg.activate("$(JULIA_PKG_DEVDIR)/IMASDD");\
+Pkg.develop(["CoordinateConventions"]);\
+Pkg.resolve();\
+try Pkg.upgrade_manifest() catch end;\
+'
+
+
 install_CoordinateConventions:
 	if [ ! -d "$(JULIA_PKG_DEVDIR)/CoordinateConventions" ]; then\
 		julia -e 'using Pkg; Pkg.develop(url="git@github.com:ProjectTorreyPines/CoordinateConventions.jl.git");';\
@@ -141,6 +154,9 @@ update_FUSE:
 update_IMAS:
 	cd $(JULIA_PKG_DEVDIR)/IMAS; git fetch; git pull
 
+update_IMASDD:
+	cd $(JULIA_PKG_DEVDIR)/IMASDD; git fetch; git pull
+
 update_CoordinateConventions:
 	cd $(JULIA_PKG_DEVDIR)/CoordinateConventions; git fetch; git pull
 

From c9cdfc86fa4720066c88265c65d16732528589ff Mon Sep 17 00:00:00 2001
From: TimSlendebroek <32385057+TimSlendebroek@users.noreply.github.com>
Date: Thu, 31 Mar 2022 16:18:33 -0700
Subject: [PATCH 02/16] bug fix ohmic power for validation database

---
 cases/HDB5.jl                            | 22 ++++++++++++++--------
 src/actors/transport_actor.jl            |  3 +++
 src/ddinit/init.jl                       |  2 +-
 src/parameters_list.jl                   |  3 ++-
 src/workflows/DB5_validation_workflow.jl | 24 ++++++++++++++++++------
 5 files changed, 38 insertions(+), 16 deletions(-)

diff --git a/cases/HDB5.jl b/cases/HDB5.jl
index 80ec048ca..a91dfe882 100644
--- a/cases/HDB5.jl
+++ b/cases/HDB5.jl
@@ -23,8 +23,8 @@ function Parameters(data_row::DataFrames.DataFrameRow)
     par.equilibrium.area = data_row[:AREA]
     par.equilibrium.volume = data_row[:VOL]
     par.equilibrium.ip = abs(data_row[:IP])
-    par.equilibrium.x_point = contains("SN", data_row[:CONFIG]) || contains("DN", data_row[:CONFIG])
-    par.equilibrium.symmetric = !par.equilibrium.x_point || !contains("SN", data_row[:CONFIG])
+    par.equilibrium.x_point = false
+    par.equilibrium.symmetric = true
 
     # Core_profiles parameters
     par.core_profiles.ne_ped = data_row[:NEL] / 1.3
@@ -38,17 +38,23 @@ function Parameters(data_row::DataFrames.DataFrameRow)
     par.core_profiles.impurity = :C
 
     # nbi
-    if data_row[:PNBI] + data_row[:POHM] > 0.0
-        par.nbi.power_launched = data_row[:PNBI] + data_row[:POHM]
-        if data_row[:PNBI] == 0.0
-            par.nbi.beam_energy = 1e9
-        else
+    if data_row[:PNBI] > 0
+        par.nbi.power_launched = data_row[:PNBI]
+        if data_row[:ENBI] > 0
             par.nbi.beam_energy = data_row[:ENBI]
+        else
+            par.nbi.beam_energy = 100e3
         end
         par.nbi.beam_mass = 2
         par.nbi.toroidal_angle = 0.0
     end
 
+    # ohmic 
+
+    if data_row[:POHM] > 0
+        par.oh.ohmic_heating = data_row[:POHM]
+    end
+
     if data_row[:PECRH] > 0
         par.ec.power_launched = data_row[:PECRH]
     end
@@ -76,4 +82,4 @@ function load_hdb5(tokamak::T=:all, extra_signal_names=T[]) where {T<:Union{Stri
         run_df = run_df[run_df.TOK.==String(tokamak), :]
     end
     return run_df
-end
\ No newline at end of file
+end                                                                                                                                                                                                   
\ No newline at end of file
diff --git a/src/actors/transport_actor.jl b/src/actors/transport_actor.jl
index ae61e28ca..598d27b10 100644
--- a/src/actors/transport_actor.jl
+++ b/src/actors/transport_actor.jl
@@ -20,6 +20,9 @@ function TauennActor(dd::IMAS.dd, par::Parameters; do_plot = false, verbose = fa
     if do_plot
         display(plot!(dd.core_profiles))
     end
+    if verbose
+        display(actor.tauenn_parameters)
+    end
     return dd
 end
 
diff --git a/src/ddinit/init.jl b/src/ddinit/init.jl
index 0f9aae5b9..8f3aa1aa0 100644
--- a/src/ddinit/init.jl
+++ b/src/ddinit/init.jl
@@ -33,7 +33,7 @@ function init(dd::IMAS.dd, par::Parameters; do_plot=false)
     end
 
     # initialize core profiles
-    if !ismissing(par.core_profiles, :bulk)
+    if !ismissing(par.core_profiles, :bulk) || par.general.init_from == :ods
         init_core_profiles(dd, par)
         if do_plot
             display(plot(dd.core_profiles))
diff --git a/src/parameters_list.jl b/src/parameters_list.jl
index 2a27dc1e0..88d2ddaa3 100644
--- a/src/parameters_list.jl
+++ b/src/parameters_list.jl
@@ -106,6 +106,7 @@ function Parameters(::Type{Val{:oh}})
     oh = Parameters(nothing)
     oh.technology = Parameters(:coil_technology)
     oh.flattop_duration = Entry(Real, "s", "Duration of the flattop (use Inf for steady-state)")
+    oh.ohmic_heating =  Entry(Real, "W", "Amount of ohmic heating absorbed in the plasma")
     return oh
 end
 
@@ -179,4 +180,4 @@ function Parameters(::Type{Val{:coil_technology}})
     coil_tech.fraction_void = Entry(Real, "", "Fraction of `void` in the coil cross-sectional area. Void is everything (like coolant) that is not structural nor conductor.")
     coil_tech.ratio_SC_to_copper = Entry(Real, "", "Fraction of superconductor to copper cross-sectional areas")
     return coil_tech
-end
\ No newline at end of file
+end
diff --git a/src/workflows/DB5_validation_workflow.jl b/src/workflows/DB5_validation_workflow.jl
index b5fd874da..ee4ce908c 100644
--- a/src/workflows/DB5_validation_workflow.jl
+++ b/src/workflows/DB5_validation_workflow.jl
@@ -11,13 +11,18 @@ import ProgressMeter
 
 Initializes and runs simple equilibrium, core_sources and transport actors and stores the resulting dd in <save_directory>
 """
-function simple_equilibrium_transport_workflow(dd::IMAS.dd, par::Parameters; save_directory::String="", do_plot::Bool=false, warn_nn_train_bounds=true)
+function simple_equilibrium_transport_workflow(dd::IMAS.dd, par::Parameters; save_directory::String="", do_plot::Bool=false, warn_nn_train_bounds=true, transport_model=:tglfnn)
     FUSE.init_equilibrium(dd, par) # already solves the equilibrium once
     FUSE.init_core_profiles(dd, par)
     FUSE.init_core_sources(dd, par)
 
+    # Add ohmic power to core_sources 
+    if !isempty(par.oh.ohmic_heating)
+        IMAS.simple_ohmic_heating_profile!(dd, par.oh.ohmic_heating)
+    end
+
     # run transport actor
-    FUSE.TauennActor(dd, par; transport_model=:tglfnn, warn_nn_train_bounds, verbose=false)
+    FUSE.TauennActor(dd, par; transport_model=transport_model, warn_nn_train_bounds, verbose=false)
 
     # Set beta_normal from equilbrium to the kinetic beta_n
     if !isempty(dd.core_profiles.profiles_1d)
@@ -61,21 +66,24 @@ function transport_validation_workflow(;
     n_samples_per_tokamak::Union{Integer,Symbol}=10,
     save_directory::String="",
     show_dd_plots=false,
-    plot_database=true)
+    plot_database=true,
+    verbose=false)
 
     # load HDB5 database
     run_df = load_hdb5(tokamak)
+    run_df = run_df[run_df.TOK .!= "ASDEX",:]
 
     # pick cases at random
     if n_samples_per_tokamak !== :all
         tok_list = unique(run_df[:, "TOK"])
+        display(@show tok_list)
         run_df = DataFrames.reduce(
             vcat, [run_df[run_df.TOK.==tok, :][Random.shuffle(1:DataFrames.nrow(run_df[run_df.TOK.==tok, :]))[1:minimum([n_samples_per_tokamak, length(run_df[run_df.TOK.==tok, :][:, "TOK"])])], :] for tok in tok_list]
         )
     end
 
     # Run simple_equilibrium_transport_workflow on each of the selected cases
-    run_df[!, "TAUTH_fuse"] = tau_FUSE = Float64[NaN for k in 1:length(run_df[:, "TOK"])]
+    tau_FUSE = zeros(length(run_df[:,"TOK"]))
     tbl = DataFrames.Tables.rowtable(run_df)
     failed_runs_ids = Int[]
     p = ProgressMeter.Progress(length(DataFrames.Tables.rows(tbl)); showspeed=true)
@@ -86,12 +94,17 @@ function transport_validation_workflow(;
             par = Parameters(run_df[idx, :])
             simple_equilibrium_transport_workflow(dd, par; save_directory, do_plot=show_dd_plots, warn_nn_train_bounds=false)
             tau_FUSE[idx] = @ddtime(dd.summary.global_quantities.tau_energy.value)
-        catch
+            if verbose
+                display(println("τ_fuse = $(@ddtime(dd.summary.global_quantities.tau_energy.value)) τ_hdb5 = $(run_df[idx, :TAUTH])"))
+            end
+        catch e
             push!(failed_runs_ids, idx)
+            display(@show e)
         end
         ProgressMeter.next!(p)
     end
     println("Failed runs: $(length(failed_runs_ids)) out of $(length(run_df[:,"TOK"]))")
+    run_df[:,"TAUTH_fuse"] = tau_FUSE
 
     # save all input data as well as predicted tau to CSV file
     if !isempty(save_directory)
@@ -99,7 +112,6 @@ function transport_validation_workflow(;
     end
 
     failed_df = run_df[failed_runs_ids, :]
-    run_df = run_df[DataFrames.completecases(run_df), :]
 
     if plot_database
         plot_x_y_regression(run_df, "TAUTH")

From 0d0857f3cdc87b545848657ca47b71dd90d4f4b2 Mon Sep 17 00:00:00 2001
From: TimSlendebroek <32385057+TimSlendebroek@users.noreply.github.com>
Date: Mon, 4 Apr 2022 16:06:02 -0700
Subject: [PATCH 03/16] take out janky ohmic power

---
 src/parameters_list.jl                   | 1 -
 src/workflows/DB5_validation_workflow.jl | 5 +----
 2 files changed, 1 insertion(+), 5 deletions(-)

diff --git a/src/parameters_list.jl b/src/parameters_list.jl
index 88d2ddaa3..fe1d939b5 100644
--- a/src/parameters_list.jl
+++ b/src/parameters_list.jl
@@ -106,7 +106,6 @@ function Parameters(::Type{Val{:oh}})
     oh = Parameters(nothing)
     oh.technology = Parameters(:coil_technology)
     oh.flattop_duration = Entry(Real, "s", "Duration of the flattop (use Inf for steady-state)")
-    oh.ohmic_heating =  Entry(Real, "W", "Amount of ohmic heating absorbed in the plasma")
     return oh
 end
 
diff --git a/src/workflows/DB5_validation_workflow.jl b/src/workflows/DB5_validation_workflow.jl
index 0bcdec588..4c928a300 100644
--- a/src/workflows/DB5_validation_workflow.jl
+++ b/src/workflows/DB5_validation_workflow.jl
@@ -17,9 +17,7 @@ function simple_equilibrium_transport_workflow(dd::IMAS.dd, par::Parameters; sav
     FUSE.init_core_sources(dd, par)
 
     # Add ohmic power to core_sources 
-    if !isempty(par.oh.ohmic_heating)
-        IMAS.simple_ohmic_heating_profile!(dd, par.oh.ohmic_heating)
-    end
+    IMAS.ohmic_power_steady_state!(dd.equilibrium.time_slice[], dd.core_profiles.profiles_1d[])
 
     # run transport actor
     FUSE.TauennActor(dd, par; transport_model=transport_model, warn_nn_train_bounds, verbose=false)
@@ -71,7 +69,6 @@ function transport_validation_workflow(;
 
     # load HDB5 database
     run_df = load_hdb5(tokamak)
-    run_df = run_df[run_df.TOK .!= "ASDEX",:]
 
     # pick cases at random
     if n_samples_per_tokamak !== :all

From 6d2fe18c25a0c3810ccd4c8bd5d74f755d799e5d Mon Sep 17 00:00:00 2001
From: TimSlendebroek <32385057+TimSlendebroek@users.noreply.github.com>
Date: Tue, 5 Apr 2022 11:24:15 -0700
Subject: [PATCH 04/16] check if dd.X in the ods that you are loading

---
 src/ddinit/init.jl | 15 ++++++++++-----
 1 file changed, 10 insertions(+), 5 deletions(-)

diff --git a/src/ddinit/init.jl b/src/ddinit/init.jl
index 8f3aa1aa0..c751c731b 100644
--- a/src/ddinit/init.jl
+++ b/src/ddinit/init.jl
@@ -1,9 +1,14 @@
 """
     init(dd::IMAS.dd, par::Parameters; do_plot = false)
 
-Initialize all IDSs
+Initialize all IDSs if it is not populated
 """
 function init(dd::IMAS.dd, par::Parameters; do_plot=false)
+    ods_items = []
+    if par.general.init_from == :ods
+        ods_items = keys(IMAS.json2imas(par.ods.filename))
+    end
+
     # initialize equilibrium
     init_equilibrium(dd, par)
     if do_plot
@@ -12,7 +17,7 @@ function init(dd::IMAS.dd, par::Parameters; do_plot=false)
     end
 
     # initialize build
-    if !ismissing(par.build, :vessel) || !ismissing(par.build, :layers)
+    if !ismissing(par.build, :vessel) || !ismissing(par.build, :layers) || :build ∈ ods_items
         init_build(dd, par)
         if do_plot
             plot(dd.equilibrium, color=:gray)
@@ -22,7 +27,7 @@ function init(dd::IMAS.dd, par::Parameters; do_plot=false)
     end
 
     # initialize oh and pf coils
-    if !ismissing(par.pf_active, :n_oh_coils)
+    if !ismissing(par.pf_active, :n_oh_coils) || :pf_active ∈ ods_items
         init_pf_active(dd, par)
         if do_plot
             plot(dd.equilibrium, color=:gray)
@@ -33,7 +38,7 @@ function init(dd::IMAS.dd, par::Parameters; do_plot=false)
     end
 
     # initialize core profiles
-    if !ismissing(par.core_profiles, :bulk) || par.general.init_from == :ods
+    if !ismissing(par.core_profiles, :bulk) || :core_profiles ∈ ods_items
         init_core_profiles(dd, par)
         if do_plot
             display(plot(dd.core_profiles))
@@ -41,7 +46,7 @@ function init(dd::IMAS.dd, par::Parameters; do_plot=false)
     end
 
     # initialize core sources
-    if !ismissing(par.ec, :power_launched) || !ismissing(par.ic, :power_launched) || !ismissing(par.lh, :power_launched) || !ismissing(par.nbi, :power_launched)
+    if !ismissing(par.ec, :power_launched) || !ismissing(par.ic, :power_launched) || !ismissing(par.lh, :power_launched) || !ismissing(par.nbi, :power_launched) || :core_sources ∈ ods_items
         init_core_sources(dd, par)
         if do_plot
             display(plot(dd.core_sources))

From fa604203640e21103eb4850f68594f55b16392bf Mon Sep 17 00:00:00 2001
From: TimSlendebroek <32385057+TimSlendebroek@users.noreply.github.com>
Date: Tue, 5 Apr 2022 15:06:42 -0700
Subject: [PATCH 05/16] equilibrium init fix , ohmic_power fix

---
 src/ddinit/init.jl                       | 7 +++++--
 src/workflows/DB5_validation_workflow.jl | 2 +-
 2 files changed, 6 insertions(+), 3 deletions(-)

diff --git a/src/ddinit/init.jl b/src/ddinit/init.jl
index c751c731b..47f7febf7 100644
--- a/src/ddinit/init.jl
+++ b/src/ddinit/init.jl
@@ -1,16 +1,19 @@
 """
     init(dd::IMAS.dd, par::Parameters; do_plot = false)
 
-Initialize all IDSs if it is not populated
+Initialize all IDSs if there are parameters for it or is initialized from ods
 """
 function init(dd::IMAS.dd, par::Parameters; do_plot=false)
     ods_items = []
+    # Check what is in the ods to load
     if par.general.init_from == :ods
         ods_items = keys(IMAS.json2imas(par.ods.filename))
     end
 
     # initialize equilibrium
-    init_equilibrium(dd, par)
+    if !ismissing(par.equilibrium.B0) || :equilibrium ∈ ods_items
+        init_equilibrium(dd, par)
+    end
     if do_plot
         plot(dd.equilibrium.time_slice[end])
         display(plot!(dd.equilibrium.time_slice[1].boundary.outline.r, dd.equilibrium.time_slice[1].boundary.outline.z))
diff --git a/src/workflows/DB5_validation_workflow.jl b/src/workflows/DB5_validation_workflow.jl
index 4c928a300..12410774e 100644
--- a/src/workflows/DB5_validation_workflow.jl
+++ b/src/workflows/DB5_validation_workflow.jl
@@ -17,7 +17,7 @@ function simple_equilibrium_transport_workflow(dd::IMAS.dd, par::Parameters; sav
     FUSE.init_core_sources(dd, par)
 
     # Add ohmic power to core_sources 
-    IMAS.ohmic_power_steady_state!(dd.equilibrium.time_slice[], dd.core_profiles.profiles_1d[])
+    IMAS.ohmic_power_steady_state!(dd)
 
     # run transport actor
     FUSE.TauennActor(dd, par; transport_model=transport_model, warn_nn_train_bounds, verbose=false)

From d3e8c3bc8b41b3d970496b48a80317a93ed8ef6b Mon Sep 17 00:00:00 2001
From: TimSlendebroek <32385057+TimSlendebroek@users.noreply.github.com>
Date: Wed, 6 Apr 2022 13:30:42 -0700
Subject: [PATCH 06/16] xpoint and symmetric added to validation cases

---
 cases/HDB5.jl          | 36 +++++++++++++++++++++++++-----------
 src/parameters_list.jl | 34 +++++++++++++++++-----------------
 2 files changed, 42 insertions(+), 28 deletions(-)

diff --git a/cases/HDB5.jl b/cases/HDB5.jl
index a91dfe882..944bce591 100644
--- a/cases/HDB5.jl
+++ b/cases/HDB5.jl
@@ -9,7 +9,7 @@ end
 
 function Parameters(data_row::DataFrames.DataFrameRow)
     par = Parameters()
-    par.general.casename = "HDB_$(data_row[:TOK])_$(data_row[:SHOT]))"
+    par.general.casename = "HDB_$(data_row[:TOK])_$(data_row[:SHOT])"
     par.general.init_from = :scalars
 
     # Equilibrium parameters
@@ -23,8 +23,28 @@ function Parameters(data_row::DataFrames.DataFrameRow)
     par.equilibrium.area = data_row[:AREA]
     par.equilibrium.volume = data_row[:VOL]
     par.equilibrium.ip = abs(data_row[:IP])
-    par.equilibrium.x_point = false
-    par.equilibrium.symmetric = true
+
+    # Determine x-points
+    if data_row[:CONFIG] == "SN"
+        # upper single null
+        xpoint = (data_row[:RGEO] * (1 - 1.1 * data_row[:DELTA] * data_row[:AMIN] / data_row[:RGEO]),  data_row[:RGEO] * 1.1 * data_row[:KAPPA] * data_row[:AMIN] / data_row[:RGEO])
+        symmetric = false
+    elseif data_row[:CONFIG] == "SN(L)"
+        # lower single null
+        xpoint = (data_row[:RGEO] * (1 - 1.1 * data_row[:DELTA] * data_row[:AMIN] / data_row[:RGEO]),  -data_row[:RGEO] * 1.1 * data_row[:KAPPA] * data_row[:AMIN] / data_row[:RGEO])
+        symmetric = false
+    elseif data_row[:CONFIG] == "DN"
+        # double null
+        xpoint = (data_row[:RGEO] * (1 - 1.1 * data_row[:DELTA] * data_row[:AMIN] / data_row[:RGEO]),  data_row[:RGEO] * 1.1 * data_row[:KAPPA] * data_row[:AMIN] / data_row[:RGEO])
+        symmetric = true
+    else
+        # no x-points
+        xpoint = false
+        symmetric = true
+    end
+
+    par.equilibrium.x_point = xpoint
+    par.equilibrium.symmetric = symmetric
 
     # Core_profiles parameters
     par.core_profiles.ne_ped = data_row[:NEL] / 1.3
@@ -49,12 +69,6 @@ function Parameters(data_row::DataFrames.DataFrameRow)
         par.nbi.toroidal_angle = 0.0
     end
 
-    # ohmic 
-
-    if data_row[:POHM] > 0
-        par.oh.ohmic_heating = data_row[:POHM]
-    end
-
     if data_row[:PECRH] > 0
         par.ec.power_launched = data_row[:PECRH]
     end
@@ -78,8 +92,8 @@ function load_hdb5(tokamak::T=:all, extra_signal_names=T[]) where {T<:Union{Stri
     run_df = run_df[DataFrames.completecases(run_df), :]
     # some basic filters
     run_df = run_df[(run_df.TOK.!="T10").&(run_df.TOK.!="TDEV").&(run_df.KAPPA.>1.0).&(1.6 .< run_df.MEFF .< 2.2).&(1.1 .< run_df.ZEFF .< 5.9), :]
-    if !(Symbol(tokamak) in [:all,:any])
+    if !(Symbol(tokamak) in [:all, :any])
         run_df = run_df[run_df.TOK.==String(tokamak), :]
     end
     return run_df
-end                                                                                                                                                                                                   
\ No newline at end of file
+end
\ No newline at end of file
diff --git a/src/parameters_list.jl b/src/parameters_list.jl
index fe1d939b5..f91e8a993 100644
--- a/src/parameters_list.jl
+++ b/src/parameters_list.jl
@@ -35,7 +35,7 @@ end
 
 function Parameters(::Type{Val{:material}})
     material = Parameters(nothing)
-    material.wall = Switch(FusionMaterials.available_materials("wall_materials"), "", "Material used for the wall"; default = "Steel, Stainless 316")
+    material.wall = Switch(FusionMaterials.available_materials("wall_materials"), "", "Material used for the wall"; default="Steel, Stainless 316")
     material.blanket = Switch(FusionMaterials.available_materials("blanket_materials"), "", "Material used for blanket coils")
     material.shield = Switch(FusionMaterials.available_materials("shield_materials"), "", "Material used for the shield")
     return material
@@ -45,20 +45,20 @@ function Parameters(::Type{Val{:equilibrium}})
     equilibrium = Parameters(nothing)
     equilibrium.B0 = Entry(Real, IMAS.equilibrium__vacuum_toroidal_field, :b0)
     equilibrium.R0 = Entry(Real, IMAS.equilibrium__vacuum_toroidal_field, :r0)
-    equilibrium.Z0 = Entry(Real, "m", "Z offset of the machine midplane"; default = 0.0)
+    equilibrium.Z0 = Entry(Real, "m", "Z offset of the machine midplane"; default=0.0)
 
-    equilibrium.volume = Entry(Real, "m³", "Scalar volume to match (optional)"; default = missing)
-    equilibrium.area = Entry(Real, "m²", "Scalar area to match (optional)"; default = missing)
+    equilibrium.volume = Entry(Real, "m³", "Scalar volume to match (optional)"; default=missing)
+    equilibrium.area = Entry(Real, "m²", "Scalar area to match (optional)"; default=missing)
 
     equilibrium.ϵ = Entry(Real, "", "Plasma aspect ratio")
     equilibrium.δ = Entry(Real, IMAS.equilibrium__time_slice___boundary, :triangularity)
     equilibrium.κ = Entry(Real, IMAS.equilibrium__time_slice___boundary, :elongation)
     equilibrium.βn = Entry(Real, IMAS.equilibrium__time_slice___global_quantities, :beta_normal)
     equilibrium.ip = Entry(Real, IMAS.equilibrium__time_slice___global_quantities, :ip)
-    equilibrium.x_point = Entry(Bool, IMAS.equilibrium__time_slice___boundary, :x_point)
+    equilibrium.x_point = Entry(Union{NTuple{2},Bool}, IMAS.equilibrium__time_slice___boundary, :x_point)
     equilibrium.symmetric = Entry(Bool, "", "Is plasma up-down symmetric")
-    equilibrium.ngrid = Entry(Int, "", "Resolution of the equilibrium grid"; default = 129)
-    equilibrium.field_null_surface = Entry(Real, "", "ψn value of the field_null_surface. Disable with 0.0"; default = 0.25)#, min=0.0, max=1.0)
+    equilibrium.ngrid = Entry(Int, "", "Resolution of the equilibrium grid"; default=129)
+    equilibrium.field_null_surface = Entry(Real, "", "ψn value of the field_null_surface. Disable with 0.0"; default=0.25)#, min=0.0, max=1.0)
     return equilibrium
 end
 
@@ -70,10 +70,10 @@ function Parameters(::Type{Val{:core_profiles}})
     core_profiles.w_ped = Entry(Real, "", "Pedestal width expressed in fraction of ψₙ")
     core_profiles.zeff = Entry(Real, "", "Effective ion charge")
     core_profiles.rot_core = Entry(Real, IMAS.core_profiles__profiles_1d, :rotation_frequency_tor_sonic)
-    core_profiles.ngrid = Entry(Int, "", "Resolution of the core_profiles grid"; default = 101)
+    core_profiles.ngrid = Entry(Int, "", "Resolution of the core_profiles grid"; default=101)
     core_profiles.bulk = Entry(Symbol, "", "Bulk ion species")
     core_profiles.impurity = Entry(Symbol, "", "Impurity ion species")
-    core_profiles.ejima = Entry(Real, "", "Ejima coefficient"; default = 0.4)
+    core_profiles.ejima = Entry(Real, "", "Ejima coefficient"; default=0.4)
     return core_profiles
 end
 
@@ -85,7 +85,7 @@ function Parameters(::Type{Val{:pf_active}})
         :corners => "like :simple, but PF coils have filaments at the four corners",
         :realistic => "hundreds of filaments per coil (very slow!)",
     ]
-    pf_active.green_model = Switch(options, "", "Model used for the Greens function calculation"; default = :simple)
+    pf_active.green_model = Switch(options, "", "Model used for the Greens function calculation"; default=:simple)
     pf_active.n_oh_coils = Entry(Int, "", "Number of OH coils")
     pf_active.n_pf_coils_inside = Entry(Int, "", "Number of PF coils inside of the TF")
     pf_active.n_pf_coils_outside = Entry(Int, "", "Number of PF coils outside of the TF")
@@ -97,7 +97,7 @@ function Parameters(::Type{Val{:tf}})
     tf = Parameters(nothing)
     tf.n_coils = Entry(Int, "", "Number of TF coils")
     options = [:princeton_D, :rectangle, :triple_arc, :miller, :spline]
-    tf.shape = Switch(options, "", "Shape of the TF coils"; default = :triple_arc)
+    tf.shape = Switch(options, "", "Shape of the TF coils"; default=:triple_arc)
     tf.technology = Parameters(:coil_technology)
     return tf
 end
@@ -111,9 +111,9 @@ end
 
 function Parameters(::Type{Val{:center_stack}})
     center_stack = Parameters(nothing)
-    center_stack.bucked = Entry(Bool, "", "flag for bucked boundary conditions between TF and OH (and center plug, if present"; default = false)
-    center_stack.noslip = Entry(Bool, "", "flag for no slip conditions between TF and OH (and center plug, if present)"; default = false)
-    center_stack.plug = Entry(Bool, "", "flag for center plug"; default = false)
+    center_stack.bucked = Entry(Bool, "", "flag for bucked boundary conditions between TF and OH (and center plug, if present"; default=false)
+    center_stack.noslip = Entry(Bool, "", "flag for no slip conditions between TF and OH (and center plug, if present)"; default=false)
+    center_stack.plug = Entry(Bool, "", "flag for center plug"; default=false)
     return center_stack
 end
 
@@ -122,8 +122,8 @@ function Parameters(::Type{Val{:nbi}})
     nbi = Parameters(nothing)
     nbi.power_launched = Entry(Union{X,Vector{X}} where {X<:Real}, "W", "Beam power")
     nbi.beam_energy = Entry(Union{X,Vector{X}} where {X<:Real}, "eV", "Beam energy")
-    nbi.beam_mass = Entry(Union{X,Vector{X}} where {X<:Real}, "AU", "Beam mass"; default = 2.0)
-    nbi.toroidal_angle = Entry(Union{X,Vector{X}} where {X<:Real}, "rad", "toroidal angle of injection"; default = 0.0)
+    nbi.beam_mass = Entry(Union{X,Vector{X}} where {X<:Real}, "AU", "Beam mass"; default=2.0)
+    nbi.toroidal_angle = Entry(Union{X,Vector{X}} where {X<:Real}, "rad", "toroidal angle of injection"; default=0.0)
     return nbi
 end
 
@@ -159,7 +159,7 @@ function Parameters(::Type{Val{:gasc}})
     gasc = Parameters(nothing)
     gasc.filename = Entry(String, "", "Output GASC .json file from which data will be loaded")
     gasc.case = Entry(Int, "", "Number of the GASC run to load")
-    gasc.no_small_gaps = Entry(Bool, "", "Remove small gaps from the GASC radial build"; default = true)
+    gasc.no_small_gaps = Entry(Bool, "", "Remove small gaps from the GASC radial build"; default=true)
     return gasc
 end
 

From 81b604ccbbfdc72f2f575f4d4d4dbbc168866aa7 Mon Sep 17 00:00:00 2001
From: TimSlendebroek <32385057+TimSlendebroek@users.noreply.github.com>
Date: Wed, 6 Apr 2022 14:03:40 -0700
Subject: [PATCH 07/16] allow passing of xpoint to SolovevEquilibriumActor

---
 src/actors/equilibrium_actor.jl | 51 +++++++++++++++++----------------
 1 file changed, 27 insertions(+), 24 deletions(-)

diff --git a/src/actors/equilibrium_actor.jl b/src/actors/equilibrium_actor.jl
index ec38a0811..eda25f22f 100644
--- a/src/actors/equilibrium_actor.jl
+++ b/src/actors/equilibrium_actor.jl
@@ -10,10 +10,10 @@ mutable struct SolovevEquilibriumActor <: AbstractActor
     S::SolovevEquilibrium
 end
 
-function SolovevEquilibriumActor(dd::IMAS.dd, par::Parameters; verbose = false)
-    actor = SolovevEquilibriumActor(dd.equilibrium, symmetric = par.equilibrium.symmetric)
+function SolovevEquilibriumActor(dd::IMAS.dd, par::Parameters; verbose=false)
+    actor = SolovevEquilibriumActor(dd.equilibrium, xpoint=par.equilibrium.x_point, symmetric=par.equilibrium.symmetric)
     step(actor; verbose)
-    finalize(actor, ngrid = par.equilibrium.ngrid)
+    finalize(actor, ngrid=par.equilibrium.ngrid)
 end
 
 """
@@ -31,9 +31,10 @@ Phys. Plasmas 17, 032502 (2010); https://doi.org/10.1063/1.3328818
 - alpha: Constant affecting the pressure
 """
 function SolovevEquilibriumActor(eq::IMAS.equilibrium;
-    qstar = 1.5,
-    alpha = 0.0,
-    symmetric = true) # symmetric should really be passed/detected through IMAS
+    qstar=1.5,
+    alpha=0.0,
+    xpoint=missing,
+    symmetric=true) # symmetric should really be passed/detected through IMAS
 
     eqt = eq.time_slice[]
     a = eqt.boundary.minor_radius
@@ -43,13 +44,15 @@ function SolovevEquilibriumActor(eq::IMAS.equilibrium;
     ϵ = a / R0
     B0 = @ddtime eq.vacuum_toroidal_field.b0
 
-    if length(eqt.boundary.x_point) > 0
-        xpoint = (eqt.boundary.x_point[1].r, eqt.boundary.x_point[1].z)
-    else
-        xpoint = nothing
+    if ismissing(xpoint)
+        if length(eqt.boundary.x_point) > 0
+            xpoint = (eqt.boundary.x_point[1].r, eqt.boundary.x_point[1].z)
+        else
+            xpoint = nothing
+        end
     end
 
-    S0 = solovev(abs(B0), R0, ϵ, δ, κ, alpha, qstar, B0_dir = sign(B0), Ip_dir = 1, symmetric = symmetric, xpoint = xpoint)
+    S0 = solovev(abs(B0), R0, ϵ, δ, κ, alpha, qstar, B0_dir=sign(B0), Ip_dir=1, symmetric=symmetric, xpoint=xpoint)
 
     SolovevEquilibriumActor(eq, S0)
 end
@@ -60,11 +63,11 @@ end
 Convert IMAS.equilibrium__time_slice to Equilibrium.jl EFIT structure
 """
 function IMAS2Equilibrium(eqt::IMAS.equilibrium__time_slice)
-    dim1 = range(eqt.profiles_2d[1].grid.dim1[1], eqt.profiles_2d[1].grid.dim1[end], length = length(eqt.profiles_2d[1].grid.dim1))
+    dim1 = range(eqt.profiles_2d[1].grid.dim1[1], eqt.profiles_2d[1].grid.dim1[end], length=length(eqt.profiles_2d[1].grid.dim1))
     @assert collect(dim1) ≈ eqt.profiles_2d[1].grid.dim1
-    dim2 = range(eqt.profiles_2d[1].grid.dim2[1], eqt.profiles_2d[1].grid.dim2[end], length = length(eqt.profiles_2d[1].grid.dim2))
+    dim2 = range(eqt.profiles_2d[1].grid.dim2[1], eqt.profiles_2d[1].grid.dim2[end], length=length(eqt.profiles_2d[1].grid.dim2))
     @assert collect(dim2) ≈ eqt.profiles_2d[1].grid.dim2
-    psi = range(eqt.profiles_1d.psi[1], eqt.profiles_1d.psi[end], length = length(eqt.profiles_1d.psi))
+    psi = range(eqt.profiles_1d.psi[1], eqt.profiles_1d.psi[end], length=length(eqt.profiles_1d.psi))
     @assert collect(psi) ≈ eqt.profiles_1d.psi
 
     Equilibrium.efit(Equilibrium.cocos(11), # COCOS
@@ -86,7 +89,7 @@ end
 
 Non-linear optimization to obtain a target `ip` and `beta_normal`
 """
-function step(actor::SolovevEquilibriumActor; verbose = false)
+function step(actor::SolovevEquilibriumActor; verbose=false)
     S0 = actor.S
 
     eqt = actor.eq.time_slice[]
@@ -97,20 +100,20 @@ function step(actor::SolovevEquilibriumActor; verbose = false)
 
     function cost(x)
         # NOTE: Ip/Beta calculation is very much off in Equilibrium.jl for diverted plasmas because boundary calculation is wrong
-        S = solovev(abs(B0), R0, epsilon, delta, kappa, x[1], x[2], B0_dir = sign(B0), Ip_dir = sign(target_ip), symmetric = true, xpoint = nothing)
+        S = solovev(abs(B0), R0, epsilon, delta, kappa, x[1], x[2], B0_dir=sign(B0), Ip_dir=sign(target_ip), symmetric=true, xpoint=nothing)
         beta_cost = (Equilibrium.beta_n(S) - target_beta) / target_beta
         ip_cost = (Equilibrium.plasma_current(S) - target_ip) / target_ip
         c = sqrt(beta_cost^2 + ip_cost^2)
         return c
     end
 
-    res = Optim.optimize(cost, [alpha, qstar], Optim.NelderMead(), Optim.Options(g_tol = 1E-3))
+    res = Optim.optimize(cost, [alpha, qstar], Optim.NelderMead(), Optim.Options(g_tol=1E-3))
 
     if verbose
         println(res)
     end
 
-    actor.S = solovev(abs(B0), R0, epsilon, delta, kappa, res.minimizer[1], res.minimizer[2], B0_dir = sign(B0), Ip_dir = sign(target_ip), symmetric = S0.symmetric, xpoint = S0.xpoint)
+    actor.S = solovev(abs(B0), R0, epsilon, delta, kappa, res.minimizer[1], res.minimizer[2], B0_dir=sign(B0), Ip_dir=sign(target_ip), symmetric=S0.symmetric, xpoint=S0.xpoint)
 
     # @show Equilibrium.beta_t(actor.S)
     # @show Equilibrium.beta_p(actor.S)
@@ -132,9 +135,9 @@ Store SolovevEquilibriumActor data in IMAS.equilibrium format
 """
 function finalize(
     actor::SolovevEquilibriumActor;
-    ngrid::Int = 129,
-    rlims::NTuple{2,<:Real} = (maximum([actor.S.R0 * (1 - actor.S.epsilon * 2), 0.0]), actor.S.R0 * (1 + actor.S.epsilon * 2)),
-    zlims::NTuple{2,<:Real} = (-actor.S.R0 * actor.S.epsilon * actor.S.kappa * 1.7, actor.S.R0 * actor.S.epsilon * actor.S.kappa * 1.7)
+    ngrid::Int=129,
+    rlims::NTuple{2,<:Real}=(maximum([actor.S.R0 * (1 - actor.S.epsilon * 2), 0.0]), actor.S.R0 * (1 + actor.S.epsilon * 2)),
+    zlims::NTuple{2,<:Real}=(-actor.S.R0 * actor.S.epsilon * actor.S.kappa * 1.7, actor.S.R0 * actor.S.epsilon * actor.S.kappa * 1.7)
 )::IMAS.equilibrium__time_slice
 
     tc = transform_cocos(3, 11)
@@ -151,7 +154,7 @@ function finalize(
     empty!(eqt)
     eqt.boundary.geometric_axis.r = actor.S.R0
     eqt.boundary.geometric_axis.z = Z0
-    orig_psi = collect(range(Equilibrium.psi_limits(actor.S)..., length = ngrid))
+    orig_psi = collect(range(Equilibrium.psi_limits(actor.S)..., length=ngrid))
     eqt.profiles_1d.psi = orig_psi * (tc["PSI"] * sign_Ip)
 
     eqt.profiles_1d.pressure = Equilibrium.pressure(actor.S, orig_psi)
@@ -162,8 +165,8 @@ function finalize(
 
     resize!(eqt.profiles_2d, 1)
     eqt.profiles_2d[1].grid_type.index = 1
-    eqt.profiles_2d[1].grid.dim1 = range(rlims[1], rlims[2], length = ngrid)
-    eqt.profiles_2d[1].grid.dim2 = range(zlims[1] + Z0, zlims[2] + Z0, length = Int(ceil(ngrid * actor.S.kappa)))
+    eqt.profiles_2d[1].grid.dim1 = range(rlims[1], rlims[2], length=ngrid)
+    eqt.profiles_2d[1].grid.dim2 = range(zlims[1] + Z0, zlims[2] + Z0, length=Int(ceil(ngrid * actor.S.kappa)))
 
     eqt.profiles_2d[1].psi = [actor.S(rr, zz - Z0) * (tc["PSI"] * sign_Ip) for rr in eqt.profiles_2d[1].grid.dim1, zz in eqt.profiles_2d[1].grid.dim2]
 

From 4735137bf6ef188d1d53d76d8070e61d2a5922a9 Mon Sep 17 00:00:00 2001
From: TimSlendebroek <32385057+TimSlendebroek@users.noreply.github.com>
Date: Thu, 7 Apr 2022 11:25:37 -0700
Subject: [PATCH 08/16] conflict merge

---
 .github/workflows/runtests.yml  |  2 +-
 Makefile                        |  1 +
 src/actors/build_actor.jl       | 15 ++++---
 src/actors/current_actor.jl     | 80 +++++++++++++++++----------------
 src/actors/equilibrium_actor.jl | 22 +++++----
 src/ddinit/init_equilibrium.jl  | 23 +++++-----
 6 files changed, 75 insertions(+), 68 deletions(-)

diff --git a/.github/workflows/runtests.yml b/.github/workflows/runtests.yml
index efc61486c..ae48d7c72 100644
--- a/.github/workflows/runtests.yml
+++ b/.github/workflows/runtests.yml
@@ -35,7 +35,7 @@ jobs:
             using Pkg
             Pkg.activate(".")
             dependencies = PackageSpec[]
-            for package in ["IMAS", "IMASDD", "CoordinateConventions", "FusionMaterials", "VacuumFields", "Equilibrium", "TAUENN", "EPEDNN", "TGLFNN", "QED"]
+            for package in ["IMAS", "IMASDD", "CoordinateConventions", "FusionMaterials", "VacuumFields", "Equilibrium", "TAUENN", "EPEDNN", "TGLFNN", "QED", "FiniteElementHermite"]
               push!(dependencies, PackageSpec(url="https://project-torrey-pines:${{secrets.PTP_READ_TOKEN}}@github.com/ProjectTorreyPines/$package.jl.git"))
             end
             Pkg.add(dependencies)
diff --git a/Makefile b/Makefile
index e255f15bb..1ed8c274f 100644
--- a/Makefile
+++ b/Makefile
@@ -29,6 +29,7 @@ Pkg.build("IJulia");\
 
 update: develop
 	make -j 100 update_FUSE update_IMAS update_IMASDD update_CoordinateConventions update_FusionMaterials update_VacuumFields update_Equilibrium update_TAUENN update_EPEDNN update_TGLFNN update_QED update_FiniteElementHermite
+	julia -e 'using Pkg; Pkg.activate("."); Pkg.precompile()'
 
 update_FUSE:
 	cd $(JULIA_PKG_DEVDIR)/FUSE; git fetch; git pull; julia -e 'using Pkg; Pkg.activate("."); Pkg.resolve()'
diff --git a/src/actors/build_actor.jl b/src/actors/build_actor.jl
index 119c6617c..dacbcf686 100644
--- a/src/actors/build_actor.jl
+++ b/src/actors/build_actor.jl
@@ -22,7 +22,7 @@ function step(flxactor::FluxSwingActor)
     cp1d = cp.profiles_1d[]
 
     bd.flux_swing_requirements.rampup = rampup_flux_requirements(eqt, cp)
-    bd.flux_swing_requirements.flattop = flattop_flux_requirements(cp1d, flxactor.flattop_duration)
+    bd.flux_swing_requirements.flattop = flattop_flux_requirements(eqt, cp1d, flxactor.flattop_duration)
     bd.flux_swing_requirements.pf = pf_flux_requirements(eqt)
 
     oh_peakJ(bd)
@@ -65,12 +65,17 @@ function rampup_flux_requirements(eqt::IMAS.equilibrium__time_slice, cp::IMAS.co
 end
 
 """
-    flattop_flux_requirements(cp1d::IMAS.core_profiles__profiles_1d, flattop_duration)
+    flattop_flux_requirements(eqt::IMAS.equilibrium__time_slice, cp1d::IMAS.core_profiles__profiles_1d, flattop_duration::Real)
 
-Estimate OH flux requirement during flattop 
+Estimate OH flux requirement during flattop (if j_ohmic profile is missing then steady state ohmic profile is assumed)
 """
-function flattop_flux_requirements(cp1d::IMAS.core_profiles__profiles_1d, flattop_duration::Real)
-    return integrate(cp1d.grid.area, cp1d.j_ohmic ./ cp1d.conductivity_parallel .* flattop_duration) # V*s
+function flattop_flux_requirements(eqt::IMAS.equilibrium__time_slice, cp1d::IMAS.core_profiles__profiles_1d, flattop_duration::Real)
+    if ismissing(cp1d, :j_ohmic)
+        j_ohmic = IMAS.j_ohmic_steady_state(eqt, cp1d)
+    else
+        j_ohmic = cp1d.j_ohmic
+    end
+    return integrate(cp1d.grid.area, j_ohmic ./ cp1d.conductivity_parallel .* flattop_duration) # V*s
 end
 
 """
diff --git a/src/actors/current_actor.jl b/src/actors/current_actor.jl
index 4ca2860c7..73adec4d3 100644
--- a/src/actors/current_actor.jl
+++ b/src/actors/current_actor.jl
@@ -16,7 +16,7 @@ function QEDcurrentActor(dd::IMAS.dd)
     QEDcurrentActor(dd, from_imas(dd), η_imas(dd), missing, @ddtime(dd.equilibrium.time), 0.0)
 end
 
-function step(actor::QEDcurrentActor, tmax, Nt, Vedge = nothing, Ip = nothing; resume=false)
+function step(actor::QEDcurrentActor, tmax, Nt, Vedge=nothing, Ip=nothing; resume=false)
     if resume
         if actor.QO !== missing
             actor.QI = actor.QO
@@ -32,53 +32,57 @@ end
 
 function finalize(actor::QEDcurrentActor)
     dd = actor.dd
-    eqt = dd.equilibrium.time_slice[]
 
+    eqt = dd.equilibrium.time_slice[]
     newtime = actor.initial_time + actor.tmax
-    resize!(dd.equilibrium.time_slice, newtime)
-    dd.equilibrium.time_slice[newtime] = new = deepcopy(eqt)
-
-    dΡ_dρ = new.profiles_1d.rho_tor[end]
-    ρ = new.profiles_1d.rho_tor / dΡ_dρ
-
-    new.profiles_1d.q = 1.0 ./ actor.QO.ι.(ρ)
-    new.profiles_1d.j_tor = actor.QO.JtoR.(ρ) ./ new.profiles_1d.gm9
-    new.time = newtime
+    resize!(dd.equilibrium.time_slice, Float64(newtime))
+    dd.equilibrium.time_slice[newtime] = eqt_new = deepcopy(eqt)
+
+    dΡ_dρ = eqt_new.profiles_1d.rho_tor[end]
+    ρ = eqt_new.profiles_1d.rho_tor / dΡ_dρ
+
+    eqt_new.profiles_1d.q = 1.0 ./ actor.QO.ι.(ρ)
+    eqt_new.profiles_1d.j_tor = actor.QO.JtoR.(ρ) ./ eqt_new.profiles_1d.gm9
+    eqt_new.time = newtime
+
+    # if dd has core_profiles set cp1d.j_total from equilibrium (and set j_ohmic as expression)
+    if !isempty(dd.core_profiles.profiles_1d)
+        cp1d_new = cp1d = dd.core_profiles.profiles_1d[]
+        if cp1d.time < newtime
+            resize!(dd.core_profiles.profiles_1d, Float64(newtime))
+            dd.core_profiles.profiles_1d[newtime] = cp1d_new = deepcopy(cp1d)
+        end
+        IMAS.j_total_from_equilibrium!(cp1d_new)
+    end
 
-    new
+    return dd
 end
 
-
 # utils
 function from_imas(dd::IMAS.dd)
     eqt = dd.equilibrium.time_slice[]
-
-    dΡ_dρ = eqt.profiles_1d.rho_tor[end]
-    ρ = eqt.profiles_1d.rho_tor / dΡ_dρ
-
-    B₀ = @ddtime(dd.equilibrium.vacuum_toroidal_field.b0)
-
-    fsa_R⁻² = QED.FE(ρ, eqt.profiles_1d.gm1)
-    F = QED.FE(ρ, eqt.profiles_1d.f)
-
-    # Require dV_dρ=0 on-axis
-    tmp = dΡ_dρ .* eqt.profiles_1d.dvolume_drho_tor
-    tmp[1] = 0.0
-    dV_dρ = QED.FE(ρ, tmp)
-
-    ι = QED.FE(ρ, 1.0 ./ eqt.profiles_1d.q)
-    JtoR = QED.FE(ρ, eqt.profiles_1d.j_tor .* eqt.profiles_1d.gm9)
-
-    JBni = nothing
-    if !ismissing(dd.core_profiles.profiles_1d[], :j_non_inductive)
-        JBni = QED.FE(dd.core_profiles.profiles_1d[].grid.rho_tor_norm, dd.core_profiles.profiles_1d[].j_non_inductive .* B₀)
+    rho_tor = eqt.profiles_1d.rho_tor
+    B0 = @ddtime(dd.equilibrium.vacuum_toroidal_field.b0)
+    gm1 = eqt.profiles_1d.gm1
+    f = eqt.profiles_1d.f
+    dvolume_drho_tor = eqt.profiles_1d.dvolume_drho_tor
+    q = eqt.profiles_1d.q
+    j_tor = eqt.profiles_1d.j_tor
+    gm9 = eqt.profiles_1d.gm9
+
+    if !isempty(dd.core_profiles.profiles_1d) && !ismissing(dd.core_profiles.profiles_1d[], :j_non_inductive)
+        prof1d = dd.core_profiles.profiles_1d[]
+        ρ_j_non_inductive = (prof1d.grid.rho_tor_norm, prof1d.j_non_inductive)
+    else
+        ρ_j_non_inductive = nothing
     end
 
-    return QED.QED_state(ρ, dΡ_dρ, B₀, fsa_R⁻², F, dV_dρ, ι, JtoR, JBni = JBni)
+    return QED.initialize(rho_tor, B0, gm1, f, dvolume_drho_tor, q, j_tor, gm9; ρ_j_non_inductive)
 end
 
-function η_imas(dd::IMAS.dd)
-    rho = dd.core_profiles.profiles_1d[].grid.rho_tor_norm
-    η = 1.0 ./ dd.core_profiles.profiles_1d[].conductivity_parallel
-    return QED.FE(rho, η)
+function η_imas(dd::IMAS.dd; use_log=true)
+    prof1d = dd.core_profiles.profiles_1d[]
+    rho = prof1d.grid.rho_tor_norm
+    η = 1.0 ./ prof1d.conductivity_parallel
+    return QED.η_FE(rho, η; use_log)
 end
\ No newline at end of file
diff --git a/src/actors/equilibrium_actor.jl b/src/actors/equilibrium_actor.jl
index eda25f22f..50e81eb8f 100644
--- a/src/actors/equilibrium_actor.jl
+++ b/src/actors/equilibrium_actor.jl
@@ -11,7 +11,7 @@ mutable struct SolovevEquilibriumActor <: AbstractActor
 end
 
 function SolovevEquilibriumActor(dd::IMAS.dd, par::Parameters; verbose=false)
-    actor = SolovevEquilibriumActor(dd.equilibrium, xpoint=par.equilibrium.x_point, symmetric=par.equilibrium.symmetric)
+    actor = SolovevEquilibriumActor(dd.equilibrium, x_point=par.equilibrium.x_point, symmetric=par.equilibrium.symmetric)
     step(actor; verbose)
     finalize(actor, ngrid=par.equilibrium.ngrid)
 end
@@ -33,7 +33,7 @@ Phys. Plasmas 17, 032502 (2010); https://doi.org/10.1063/1.3328818
 function SolovevEquilibriumActor(eq::IMAS.equilibrium;
     qstar=1.5,
     alpha=0.0,
-    xpoint=missing,
+    x_point=missing,
     symmetric=true) # symmetric should really be passed/detected through IMAS
 
     eqt = eq.time_slice[]
@@ -44,15 +44,13 @@ function SolovevEquilibriumActor(eq::IMAS.equilibrium;
     ϵ = a / R0
     B0 = @ddtime eq.vacuum_toroidal_field.b0
 
-    if ismissing(xpoint)
-        if length(eqt.boundary.x_point) > 0
-            xpoint = (eqt.boundary.x_point[1].r, eqt.boundary.x_point[1].z)
-        else
-            xpoint = nothing
-        end
+    if length(eqt.boundary.x_point) > 0
+        x_point = (eqt.boundary.x_point[1].r, eqt.boundary.x_point[1].z)
+    else
+        x_point = nothing
     end
-
-    S0 = solovev(abs(B0), R0, ϵ, δ, κ, alpha, qstar, B0_dir=sign(B0), Ip_dir=1, symmetric=symmetric, xpoint=xpoint)
+    @show x_point
+    S0 = solovev(abs(B0), R0, ϵ, δ, κ, alpha, qstar, B0_dir=sign(B0), Ip_dir=1, symmetric=symmetric, x_point=x_point)
 
     SolovevEquilibriumActor(eq, S0)
 end
@@ -100,7 +98,7 @@ function step(actor::SolovevEquilibriumActor; verbose=false)
 
     function cost(x)
         # NOTE: Ip/Beta calculation is very much off in Equilibrium.jl for diverted plasmas because boundary calculation is wrong
-        S = solovev(abs(B0), R0, epsilon, delta, kappa, x[1], x[2], B0_dir=sign(B0), Ip_dir=sign(target_ip), symmetric=true, xpoint=nothing)
+        S = solovev(abs(B0), R0, epsilon, delta, kappa, x[1], x[2], B0_dir=sign(B0), Ip_dir=sign(target_ip), symmetric=true, x_point=nothing)
         beta_cost = (Equilibrium.beta_n(S) - target_beta) / target_beta
         ip_cost = (Equilibrium.plasma_current(S) - target_ip) / target_ip
         c = sqrt(beta_cost^2 + ip_cost^2)
@@ -113,7 +111,7 @@ function step(actor::SolovevEquilibriumActor; verbose=false)
         println(res)
     end
 
-    actor.S = solovev(abs(B0), R0, epsilon, delta, kappa, res.minimizer[1], res.minimizer[2], B0_dir=sign(B0), Ip_dir=sign(target_ip), symmetric=S0.symmetric, xpoint=S0.xpoint)
+    actor.S = solovev(abs(B0), R0, epsilon, delta, kappa, res.minimizer[1], res.minimizer[2], B0_dir=sign(B0), Ip_dir=sign(target_ip), symmetric=S0.symmetric, x_point=S0.x_point)
 
     # @show Equilibrium.beta_t(actor.S)
     # @show Equilibrium.beta_p(actor.S)
diff --git a/src/ddinit/init_equilibrium.jl b/src/ddinit/init_equilibrium.jl
index 50645c33d..fdf4a1442 100644
--- a/src/ddinit/init_equilibrium.jl
+++ b/src/ddinit/init_equilibrium.jl
@@ -27,7 +27,7 @@ function init_equilibrium(
     δ::Real,
     βn::Real,
     ip::Real,
-    x_point::Union{AbstractVector,NTuple{2},Bool} = false)
+    x_point::Union{AbstractVector,NTuple{2},Bool}=false)
 
     eqt = resize!(eq.time_slice)
     eqt.boundary.minor_radius = ϵ * R0
@@ -71,21 +71,20 @@ function init_equilibrium(dd::IMAS.dd, par::Parameters)
         # init equilibrium
         init_equilibrium(
             dd.equilibrium;
-            B0 = par.equilibrium.B0,
-            R0 = par.equilibrium.R0,
-            Z0 = par.equilibrium.Z0,
-            ϵ = par.equilibrium.ϵ,
-            κ = par.equilibrium.κ,
-            δ = par.equilibrium.δ,
-            βn = par.equilibrium.βn,
-            ip = par.equilibrium.ip,
-            x_point = par.equilibrium.x_point)
+            B0=par.equilibrium.B0,
+            R0=par.equilibrium.R0,
+            Z0=par.equilibrium.Z0,
+            ϵ=par.equilibrium.ϵ,
+            κ=par.equilibrium.κ,
+            δ=par.equilibrium.δ,
+            βn=par.equilibrium.βn,
+            ip=par.equilibrium.ip,
+            x_point=par.equilibrium.x_point)
 
         # equilibrium
         SolovevEquilibriumActor(dd, par)
     end
 
-
     # field null surface
     if par.equilibrium.field_null_surface > 0.0
         pushfirst!(dd.equilibrium.time_slice, field_null_surface(dd.equilibrium.time_slice[], par.equilibrium.field_null_surface))
@@ -102,7 +101,7 @@ end
 
 Return field null surface by scaling an existing equilibrium time_slice
 """
-function field_null_surface(eqt::IMAS.equilibrium__time_slice, scale::Real = 0.25, abs_psi_boundary::Real = 0.1)
+function field_null_surface(eqt::IMAS.equilibrium__time_slice, scale::Real=0.25, abs_psi_boundary::Real=0.1)
     eqb = IMAS.equilibrium__time_slice()
     eqb.global_quantities.psi_boundary = sign(eqt.profiles_1d.psi[1] - eqt.profiles_1d.psi[end]) * abs_psi_boundary
     eqb.boundary.outline.r, eqb.boundary.outline.z, _ = IMAS.flux_surface(eqt, eqt.profiles_1d.psi[1] * (1 - scale) + eqt.profiles_1d.psi[end] * scale)

From e802f4604fb4abfdfba00618051018ba8ad8be56 Mon Sep 17 00:00:00 2001
From: TimSlendebroek <32385057+TimSlendebroek@users.noreply.github.com>
Date: Thu, 7 Apr 2022 13:45:12 -0700
Subject: [PATCH 09/16] Symmetric and x_point resolution

Allowing for upper/lower and double null
---
 cases/HDB5.jl                            | 17 ++++++++--------
 sample/fpp_gasc_59_step.json             |  4 ----
 src/actors/equilibrium_actor.jl          | 26 ++++++++++++++++--------
 src/ddinit/init_equilibrium.jl           | 14 ++++++++++---
 src/workflows/DB5_validation_workflow.jl | 13 ++++++------
 5 files changed, 44 insertions(+), 30 deletions(-)

diff --git a/cases/HDB5.jl b/cases/HDB5.jl
index 944bce591..88a9bce23 100644
--- a/cases/HDB5.jl
+++ b/cases/HDB5.jl
@@ -13,7 +13,7 @@ function Parameters(data_row::DataFrames.DataFrameRow)
     par.general.init_from = :scalars
 
     # Equilibrium parameters
-    par.equilibrium.B0 = abs(data_row[:BT])
+    par.equilibrium.B0 = data_row[:BT]
     par.equilibrium.R0 = data_row[:RGEO]
     par.equilibrium.Z0 = 0.0
     par.equilibrium.ϵ = data_row[:AMIN] / data_row[:RGEO]
@@ -22,28 +22,30 @@ function Parameters(data_row::DataFrames.DataFrameRow)
     par.equilibrium.βn = 1.0
     par.equilibrium.area = data_row[:AREA]
     par.equilibrium.volume = data_row[:VOL]
-    par.equilibrium.ip = abs(data_row[:IP])
+    par.equilibrium.ip = data_row[:IP]
+
+    x_point = (data_row[:RGEO] * (1 - 1.1 * data_row[:DELTA] * data_row[:AMIN] / data_row[:RGEO]), data_row[:RGEO] * 1.1 * data_row[:KAPPA] * data_row[:AMIN] / data_row[:RGEO])
 
     # Determine x-points
     if data_row[:CONFIG] == "SN"
         # upper single null
-        xpoint = (data_row[:RGEO] * (1 - 1.1 * data_row[:DELTA] * data_row[:AMIN] / data_row[:RGEO]),  data_row[:RGEO] * 1.1 * data_row[:KAPPA] * data_row[:AMIN] / data_row[:RGEO])
+        x_point = (x_point[1], x_point[2])
         symmetric = false
     elseif data_row[:CONFIG] == "SN(L)"
         # lower single null
-        xpoint = (data_row[:RGEO] * (1 - 1.1 * data_row[:DELTA] * data_row[:AMIN] / data_row[:RGEO]),  -data_row[:RGEO] * 1.1 * data_row[:KAPPA] * data_row[:AMIN] / data_row[:RGEO])
+        x_point = (x_point[1], -x_point[2])
         symmetric = false
     elseif data_row[:CONFIG] == "DN"
         # double null
-        xpoint = (data_row[:RGEO] * (1 - 1.1 * data_row[:DELTA] * data_row[:AMIN] / data_row[:RGEO]),  data_row[:RGEO] * 1.1 * data_row[:KAPPA] * data_row[:AMIN] / data_row[:RGEO])
+        x_point = (x_point[1], x_point[2])
         symmetric = true
     else
         # no x-points
-        xpoint = false
+        x_point = false
         symmetric = true
     end
 
-    par.equilibrium.x_point = xpoint
+    par.equilibrium.x_point = x_point
     par.equilibrium.symmetric = symmetric
 
     # Core_profiles parameters
@@ -76,7 +78,6 @@ function Parameters(data_row::DataFrames.DataFrameRow)
     if data_row[:PICRH] > 0
         par.ic.power_launched = data_row[:PICRH]
     end
-
     return par
 end
 
diff --git a/sample/fpp_gasc_59_step.json b/sample/fpp_gasc_59_step.json
index f1f5e22ae..1dd7e7936 100644
--- a/sample/fpp_gasc_59_step.json
+++ b/sample/fpp_gasc_59_step.json
@@ -22534,10 +22534,6 @@
 },
 "x_point": [
 {
-"r": NaN,
-"z": NaN
-},
-{
 "r": 5.029669610853508,
 "z": -0.001802959285504403
 }
diff --git a/src/actors/equilibrium_actor.jl b/src/actors/equilibrium_actor.jl
index 6cd254f9f..398f75fad 100644
--- a/src/actors/equilibrium_actor.jl
+++ b/src/actors/equilibrium_actor.jl
@@ -11,7 +11,7 @@ mutable struct SolovevEquilibriumActor <: AbstractActor
 end
 
 function SolovevEquilibriumActor(dd::IMAS.dd, par::Parameters; verbose=false)
-    actor = SolovevEquilibriumActor(dd.equilibrium, x_point=par.equilibrium.x_point, symmetric=par.equilibrium.symmetric)
+    actor = SolovevEquilibriumActor(dd.equilibrium)
     step(actor; verbose)
     finalize(actor, ngrid=par.equilibrium.ngrid)
 end
@@ -32,9 +32,7 @@ Phys. Plasmas 17, 032502 (2010); https://doi.org/10.1063/1.3328818
 """
 function SolovevEquilibriumActor(eq::IMAS.equilibrium;
     qstar=1.5,
-    alpha=0.0,
-    x_point=missing,
-    symmetric=true) # symmetric should really be passed/detected through IMAS
+    alpha=0.0)
 
     eqt = eq.time_slice[]
     a = eqt.boundary.minor_radius
@@ -44,13 +42,19 @@ function SolovevEquilibriumActor(eq::IMAS.equilibrium;
     ϵ = a / R0
     B0 = @ddtime eq.vacuum_toroidal_field.b0
 
+    # check number of x_points
+    if mod(length(eqt.boundary.x_point), 2) == 0
+        symmetric = true
+    else
+        symmetric = false
+    end
+
     if length(eqt.boundary.x_point) > 0
-        x_point = (eqt.boundary.x_point[1].r, eqt.boundary.x_point[1].z)
+        x_point = (eqt.boundary.x_point[1].r, -abs(eqt.boundary.x_point[1].z))
     else
         x_point = nothing
     end
-    @show x_point
-    S0 = solovev(abs(B0), R0, ϵ, δ, κ, alpha, qstar, B0_dir=sign(B0), Ip_dir=1, symmetric=symmetric, x_point=x_point)
+    S0 = Equilibrium.solovev(abs(B0), R0, ϵ, δ, κ, alpha, qstar, B0_dir=Int64(sign(B0)), Ip_dir=1, x_point=x_point, symmetric=symmetric)
 
     SolovevEquilibriumActor(eq, S0)
 end
@@ -144,7 +148,12 @@ function finalize(
     eqt = eq.time_slice[]
     sign_Ip = sign(eqt.global_quantities.ip)
     sign_Bt = sign(eqt.profiles_1d.f[end])
+
     Z0 = eqt.boundary.geometric_axis.z
+    flip_z = 1.0
+    if mod(length(eqt.boundary.x_point), 2) == 1 && eqt.boundary.x_point[1].z > Z0
+        flip_z = -1.0
+    end
 
     eq.vacuum_toroidal_field.r0 = actor.S.R0
     @ddtime eq.vacuum_toroidal_field.b0 = actor.S.B0 * sign_Bt
@@ -166,8 +175,7 @@ function finalize(
     eqt.profiles_2d[1].grid.dim1 = range(rlims[1], rlims[2], length=ngrid)
     eqt.profiles_2d[1].grid.dim2 = range(zlims[1] + Z0, zlims[2] + Z0, length=Int(ceil(ngrid * actor.S.kappa)))
 
-    eqt.profiles_2d[1].psi = [actor.S(rr, zz - Z0) * (tc["PSI"] * sign_Ip) for rr in eqt.profiles_2d[1].grid.dim1, zz in eqt.profiles_2d[1].grid.dim2]
-
+    eqt.profiles_2d[1].psi = [actor.S(rr, flip_z * (zz - Z0)) * (tc["PSI"] * sign_Ip) for rr in eqt.profiles_2d[1].grid.dim1, zz in eqt.profiles_2d[1].grid.dim2]
     IMAS.flux_surfaces(eqt)
 
     return eqt
diff --git a/src/ddinit/init_equilibrium.jl b/src/ddinit/init_equilibrium.jl
index fdf4a1442..72fd2a126 100644
--- a/src/ddinit/init_equilibrium.jl
+++ b/src/ddinit/init_equilibrium.jl
@@ -13,7 +13,8 @@
         δ::Real,
         βn::Real,
         ip::Real,
-        x_point::Union{Vector,NTuple{2},Bool} = false)
+        x_point::Union{Vector,NTuple{2},Bool} = false,
+        symmetric::Bool=true)
 
 Initialize equilibrium IDS based on some basic Miller geometry parameters
 """
@@ -27,7 +28,8 @@ function init_equilibrium(
     δ::Real,
     βn::Real,
     ip::Real,
-    x_point::Union{AbstractVector,NTuple{2},Bool}=false)
+    x_point::Union{AbstractVector,NTuple{2},Bool}=false,
+    symmetric::Bool=true)
 
     eqt = resize!(eq.time_slice)
     eqt.boundary.minor_radius = ϵ * R0
@@ -44,6 +46,11 @@ function init_equilibrium(
         resize!(eqt.boundary.x_point, 1)
         eqt.boundary.x_point[1].r = x_point[1]
         eqt.boundary.x_point[1].z = x_point[2]
+        if symmetric
+            resize!(eqt.boundary.x_point, 2)
+            eqt.boundary.x_point[2].r = x_point[1]
+            eqt.boundary.x_point[2].z = -x_point[2]
+        end
     end
     eq.vacuum_toroidal_field.r0 = R0
     @ddtime eq.vacuum_toroidal_field.b0 = B0
@@ -79,7 +86,8 @@ function init_equilibrium(dd::IMAS.dd, par::Parameters)
             δ=par.equilibrium.δ,
             βn=par.equilibrium.βn,
             ip=par.equilibrium.ip,
-            x_point=par.equilibrium.x_point)
+            x_point=par.equilibrium.x_point,
+            symmetric=par.equilibrium.symmetric)
 
         # equilibrium
         SolovevEquilibriumActor(dd, par)
diff --git a/src/workflows/DB5_validation_workflow.jl b/src/workflows/DB5_validation_workflow.jl
index 12410774e..22f4ad0bb 100644
--- a/src/workflows/DB5_validation_workflow.jl
+++ b/src/workflows/DB5_validation_workflow.jl
@@ -16,8 +16,10 @@ function simple_equilibrium_transport_workflow(dd::IMAS.dd, par::Parameters; sav
     FUSE.init_core_profiles(dd, par)
     FUSE.init_core_sources(dd, par)
 
-    # Add ohmic power to core_sources 
-    IMAS.ohmic_power_steady_state!(dd)
+    # Set j_ohmic to steady_state ohmic current 
+    IMAS.j_ohmic_steady_state!(dd.equilibrium.time_slice[], dd.core_profiles.profiles_1d[])
+    # Add ohmic heating
+    IMAS.ohmic_source!(dd)
 
     # run transport actor
     FUSE.TauennActor(dd, par; transport_model=transport_model, warn_nn_train_bounds, verbose=false)
@@ -73,7 +75,6 @@ function transport_validation_workflow(;
     # pick cases at random
     if n_samples_per_tokamak !== :all
         tok_list = unique(run_df[:, "TOK"])
-        display(@show tok_list)
         run_df = DataFrames.reduce(
             vcat, [run_df[run_df.TOK.==tok, :][Random.shuffle(1:DataFrames.nrow(run_df[run_df.TOK.==tok, :]))[1:minimum([n_samples_per_tokamak, length(run_df[run_df.TOK.==tok, :][:, "TOK"])])], :] for tok in tok_list]
         )
@@ -95,7 +96,9 @@ function transport_validation_workflow(;
             end
         catch e
             push!(failed_runs_ids, idx)
-            display(@show e)
+            if verbose
+                display(@show e)
+            end
         end
         ProgressMeter.next!(p)
     end
@@ -110,8 +113,6 @@ function transport_validation_workflow(;
         CSV.write(joinpath(save_directory, "failed_runs_dataframe.csv"), failed_df)
     end
 
-
-
     if plot_database
         plot_x_y_regression(run_df, "TAUTH")
     end

From 261aa37af044974e907dc36a5f031e173468ac46 Mon Sep 17 00:00:00 2001
From: TimSlendebroek <32385057+TimSlendebroek@users.noreply.github.com>
Date: Thu, 7 Apr 2022 14:05:13 -0700
Subject: [PATCH 10/16] merge conflict bugs....

---
 src/actors/equilibrium_actor.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/actors/equilibrium_actor.jl b/src/actors/equilibrium_actor.jl
index 398f75fad..6212c1b18 100644
--- a/src/actors/equilibrium_actor.jl
+++ b/src/actors/equilibrium_actor.jl
@@ -102,7 +102,7 @@ function step(actor::SolovevEquilibriumActor; verbose=false)
 
     function cost(x)
         # NOTE: Ip/Beta calculation is very much off in Equilibrium.jl for diverted plasmas because boundary calculation is wrong
-        S = solovev(abs(B0), R0, epsilon, delta, kappa, x[1], x[2], B0_dir = sign(B0), Ip_dir = 1, symmetric = true, xpoint = nothing)
+        S = solovev(abs(B0), R0, epsilon, delta, kappa, x[1], x[2], B0_dir = sign(B0), Ip_dir = 1, symmetric = true, x_point = nothing)
         beta_cost = (Equilibrium.beta_n(S) - target_beta) / target_beta
         ip_cost = (Equilibrium.plasma_current(S) - target_ip) / target_ip
         c = sqrt(beta_cost^2 + ip_cost^2)

From 206d799ad3ab8c0af622926b0ec470fdf11af4dc Mon Sep 17 00:00:00 2001
From: TimSlendebroek <32385057+TimSlendebroek@users.noreply.github.com>
Date: Thu, 7 Apr 2022 15:01:25 -0700
Subject: [PATCH 11/16] j_ohmic_steady_state! function changed

---
 src/workflows/DB5_validation_workflow.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/workflows/DB5_validation_workflow.jl b/src/workflows/DB5_validation_workflow.jl
index 22f4ad0bb..9c80a8d6f 100644
--- a/src/workflows/DB5_validation_workflow.jl
+++ b/src/workflows/DB5_validation_workflow.jl
@@ -17,7 +17,8 @@ function simple_equilibrium_transport_workflow(dd::IMAS.dd, par::Parameters; sav
     FUSE.init_core_sources(dd, par)
 
     # Set j_ohmic to steady_state ohmic current 
-    IMAS.j_ohmic_steady_state!(dd.equilibrium.time_slice[], dd.core_profiles.profiles_1d[])
+    IMAS.j_ohmic_steady_state!(dd.core_profiles.profiles_1d[])
+
     # Add ohmic heating
     IMAS.ohmic_source!(dd)
 

From f0cd2861c424b3036bb844f6646dece2400c17f8 Mon Sep 17 00:00:00 2001
From: TimSlendebroek <32385057+TimSlendebroek@users.noreply.github.com>
Date: Thu, 7 Apr 2022 15:30:46 -0700
Subject: [PATCH 12/16] !ismissing mistake fixed

---
 src/ddinit/init.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/ddinit/init.jl b/src/ddinit/init.jl
index 7dc57e9ec..f051d308b 100644
--- a/src/ddinit/init.jl
+++ b/src/ddinit/init.jl
@@ -11,7 +11,7 @@ function init(dd::IMAS.dd, par::Parameters; do_plot=false)
     end
 
     # initialize equilibrium
-    if !ismissing(par.equilibrium.B0) || :equilibrium ∈ ods_items
+    if !ismissing(par.equilibrium, :B0) || :equilibrium ∈ ods_items
         init_equilibrium(dd, par)
     end
     if do_plot

From 6e1536b735f9ab68d843d086119a02f334eecf7b Mon Sep 17 00:00:00 2001
From: Orso Meneghini <orso82@gmail.com>
Date: Thu, 7 Apr 2022 21:39:13 -0700
Subject: [PATCH 13/16] better finding of x-points

---
 src/actors/current_actor.jl              |  2 +-
 src/ddinit/init_build.jl                 |  2 +-
 src/physics.jl                           |  2 +-
 src/utils.jl                             | 19 -------------------
 src/workflows/DB5_validation_workflow.jl |  9 +++------
 5 files changed, 6 insertions(+), 28 deletions(-)

diff --git a/src/actors/current_actor.jl b/src/actors/current_actor.jl
index 73adec4d3..544cc85db 100644
--- a/src/actors/current_actor.jl
+++ b/src/actors/current_actor.jl
@@ -52,7 +52,7 @@ function finalize(actor::QEDcurrentActor)
             resize!(dd.core_profiles.profiles_1d, Float64(newtime))
             dd.core_profiles.profiles_1d[newtime] = cp1d_new = deepcopy(cp1d)
         end
-        IMAS.j_total_from_equilibrium!(cp1d_new)
+        IMAS.j_total_from_equilibrium!(eqt_new, cp1d_new)
     end
 
     return dd
diff --git a/src/ddinit/init_build.jl b/src/ddinit/init_build.jl
index db16d579d..0045d3205 100644
--- a/src/ddinit/init_build.jl
+++ b/src/ddinit/init_build.jl
@@ -264,7 +264,7 @@ function wall_from_eq(bd::IMAS.build, eqt::IMAS.equilibrium__time_slice; diverto
         if sign(pz[1] - Z0) != sign(pz[end] - Z0)
             # open flux surface does not encicle the plasma
             continue
-        elseif minimum_distance_two_shapes(pr, pz, rlcfs, zlcfs) > (maximum(zlcfs) - minimum(zlcfs)) / 20
+        elseif IMAS.minimum_distance_two_shapes(pr, pz, rlcfs, zlcfs) > (maximum(zlcfs) - minimum(zlcfs)) / 20
             # secondary Xpoint far away
             continue
         elseif (sum(pz) - Z0) < 0
diff --git a/src/physics.jl b/src/physics.jl
index dbb523f94..2ffc29679 100644
--- a/src/physics.jl
+++ b/src/physics.jl
@@ -139,7 +139,7 @@ function optimize_shape(r_obstruction, z_obstruction, target_clearance, func, r_
         cost_inside = sum(inpoly) / cost_area
 
         # target clearance  O(1)
-        minimum_distance = minimum_distance_two_shapes(R, Z, r_obstruction, z_obstruction)
+        minimum_distance = IMAS.minimum_distance_two_shapes(R, Z, r_obstruction, z_obstruction)
         cost_min_clearance = (minimum_distance - target_clearance) / target_clearance
 
         # favor up/down symmetric solutions
diff --git a/src/utils.jl b/src/utils.jl
index 6a261c873..853780460 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -28,25 +28,6 @@ function atan_eq(r, z, r0, z0)
     return r, z, θ
 end
 
-"""
-    minimum_distance_two_shapes(R_obj1, Z_obj1, R_obj2, Z_obj2)
-
-Returns minimum distance between two shapes
-"""
-function minimum_distance_two_shapes(R_obj1, Z_obj1, R_obj2, Z_obj2)
-    R_obj1, Z_obj1, R_obj2, Z_obj2 = promote(R_obj1, Z_obj1, R_obj2, Z_obj2)
-    distance = Inf
-    for k1 in 1:length(R_obj1)
-        for k2 in 1:length(R_obj2)
-            @inbounds d = (R_obj1[k1] - R_obj2[k2])^2 + (Z_obj1[k1] - Z_obj2[k2])^2
-            if distance > d
-                distance = d
-            end
-        end
-    end
-    return sqrt(distance)
-end
-
 struct GASC
     filename::String
     case::Int
diff --git a/src/workflows/DB5_validation_workflow.jl b/src/workflows/DB5_validation_workflow.jl
index 9c80a8d6f..291363753 100644
--- a/src/workflows/DB5_validation_workflow.jl
+++ b/src/workflows/DB5_validation_workflow.jl
@@ -11,19 +11,16 @@ import ProgressMeter
 
 Initializes and runs simple equilibrium, core_sources and transport actors and stores the resulting dd in <save_directory>
 """
-function simple_equilibrium_transport_workflow(dd::IMAS.dd, par::Parameters; save_directory::String="", do_plot::Bool=false, warn_nn_train_bounds=true, transport_model=:tglfnn)
+function simple_equilibrium_transport_workflow(dd::IMAS.dd, par::Parameters; save_directory::String="", do_plot::Bool=false, warn_nn_train_bounds=true, transport_model=:tglfnn, verbose=false)
     FUSE.init_equilibrium(dd, par) # already solves the equilibrium once
     FUSE.init_core_profiles(dd, par)
     FUSE.init_core_sources(dd, par)
 
     # Set j_ohmic to steady_state ohmic current 
-    IMAS.j_ohmic_steady_state!(dd.core_profiles.profiles_1d[])
-
-    # Add ohmic heating
-    IMAS.ohmic_source!(dd)
+    IMAS.j_ohmic_steady_state!(dd.equilibrium.time_slice[], dd.core_profiles.profiles_1d[])
 
     # run transport actor
-    FUSE.TauennActor(dd, par; transport_model=transport_model, warn_nn_train_bounds, verbose=false)
+    FUSE.TauennActor(dd, par; transport_model=transport_model, warn_nn_train_bounds, verbose=verbose)
 
     # Set beta_normal from equilbrium to the kinetic beta_n
     if !isempty(dd.core_profiles.profiles_1d)

From 0106b30d09ac0c2b0db761960de42dd902c34bf6 Mon Sep 17 00:00:00 2001
From: Orso Meneghini <orso82@gmail.com>
Date: Thu, 7 Apr 2022 21:40:17 -0700
Subject: [PATCH 14/16] init_workflow --> init

---
 test/runtests_workflows.jl | 20 ++++++++++----------
 1 file changed, 10 insertions(+), 10 deletions(-)

diff --git a/test/runtests_workflows.jl b/test/runtests_workflows.jl
index bc215f1e5..96c3781d1 100644
--- a/test/runtests_workflows.jl
+++ b/test/runtests_workflows.jl
@@ -2,36 +2,36 @@ using Revise
 using FUSE
 using Test
 
-@testset "init_workflows" begin
+@testset "init" begin
 
     @testset "ITER_ods" begin
         dd = IMAS.dd()
-        par = FUSE.Parameters(:ITER; init_from = :ods)
-        FUSE.init_workflow(dd, par)
+        par = FUSE.Parameters(:ITER; init_from=:ods)
+        FUSE.init(dd, par)
     end
 
     @testset "ITER_scalars" begin
         dd = IMAS.dd()
-        par = FUSE.Parameters(:ITER; init_from = :scalars)
-        FUSE.init_workflow(dd, par)
+        par = FUSE.Parameters(:ITER; init_from=:scalars)
+        FUSE.init(dd, par)
     end
 
     @testset "CAT" begin
         dd = IMAS.dd()
         par = FUSE.Parameters(:CAT)
-        FUSE.init_workflow(dd, par)
+        FUSE.init(dd, par)
     end
 
     @testset "FPP_gasc" begin
         dd = IMAS.dd()
-        par = FUSE.Parameters(:FPP, init_from = :gasc)
-        FUSE.init_workflow(dd, par)
+        par = FUSE.Parameters(:FPP, init_from=:gasc)
+        FUSE.init(dd, par)
     end
 
     @testset "FPP_ods" begin
         dd = IMAS.dd()
-        par = FUSE.Parameters(:FPP, init_from = :ods)
-        FUSE.init_workflow(dd, par)
+        par = FUSE.Parameters(:FPP, init_from=:ods)
+        FUSE.init(dd, par)
     end
 
 end
\ No newline at end of file

From 843b716697ad1ad1eff3657eb7f53801d5f55f80 Mon Sep 17 00:00:00 2001
From: Orso Meneghini <orso82@gmail.com>
Date: Thu, 7 Apr 2022 22:15:25 -0700
Subject: [PATCH 15/16] recalculate flux_surfaces when loading from ODS

---
 src/ddinit/init_equilibrium.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/ddinit/init_equilibrium.jl b/src/ddinit/init_equilibrium.jl
index 72fd2a126..e5b021cb8 100644
--- a/src/ddinit/init_equilibrium.jl
+++ b/src/ddinit/init_equilibrium.jl
@@ -69,6 +69,7 @@ function init_equilibrium(dd::IMAS.dd, par::Parameters)
         if !ismissing(dd1.equilibrium, :time) && length(keys(dd1.equilibrium.time)) > 0
             dd.global_time = max(dd.global_time, maximum(dd1.equilibrium.time))
             dd.equilibrium = dd1.equilibrium
+            IMAS.flux_surfaces(dd.equilibrium.time_slice[])
         else
             init_from = :scalars
         end

From ad7f714a750ddcb54217e12658a995f49595c1e0 Mon Sep 17 00:00:00 2001
From: Orso Meneghini <orso82@gmail.com>
Date: Thu, 7 Apr 2022 22:29:21 -0700
Subject: [PATCH 16/16] minor renaming of sources

---
 src/actors/sources_actor.jl | 36 ++++++++++++++++++++++++++++--------
 1 file changed, 28 insertions(+), 8 deletions(-)

diff --git a/src/actors/sources_actor.jl b/src/actors/sources_actor.jl
index 3add817b3..5f82df4dc 100644
--- a/src/actors/sources_actor.jl
+++ b/src/actors/sources_actor.jl
@@ -44,10 +44,15 @@ function step(actor::SimpleNBIactor)
         ne_vol = integrate(volume_cp, cp1d.electrons.density) / volume_cp[end]
         j_parallel = actor.current_efficiency / eqt.boundary.geometric_axis.r / (ne_vol / 1e19) * power_launched
 
-        isource = resize!(cs.source, "identifier.name" => "beam_$idx")
+        name = "nbi"
+        if length(actor.dd.ic_antennas.antenna) > 1
+            name = "nbi_$idx"
+        end
+
+        isource = resize!(cs.source, "identifier.name" => name)
         gaussian_source_to_dd(
             isource,
-            "beam_$idx",
+            name,
             2,
             rho_cp,
             volume_cp,
@@ -103,10 +108,15 @@ function step(actor::SimpleECactor)
         ne_vol = integrate(volume_cp, cp1d.electrons.density) / volume_cp[end]
         j_parallel = actor.current_efficiency / eqt.boundary.geometric_axis.r / (ne_vol / 1e19) * power_launched
 
-        isource = resize!(cs.source, "identifier.name" => "ec_launcher_$idx")
+        name = "ec"
+        if length(actor.dd.ic_antennas.antenna) > 1
+            name = "ec_$idx"
+        end
+
+        isource = resize!(cs.source, "identifier.name" => name)
         gaussian_source_to_dd(
             isource,
-            "ec_launcher_$idx",
+            name,
             3,
             rho_cp,
             volume_cp,
@@ -160,10 +170,15 @@ function step(actor::SimpleICactor)
         ne_vol = integrate(volume_cp, cp1d.electrons.density) / volume_cp[end]
         j_parallel = actor.current_efficiency / eqt.boundary.geometric_axis.r / (ne_vol / 1e19) * power_launched
 
-        isource = resize!(cs.source, "identifier.name" => "ic_antenna_$idx")
+        name = "ic"
+        if length(actor.dd.ic_antennas.antenna) > 1
+            name = "ic_$idx"
+        end
+
+        isource = resize!(cs.source, "identifier.name" => name)
         gaussian_source_to_dd(
             isource,
-            "ic_antenna_$idx",
+            name,
             5,
             rho_cp,
             volume_cp,
@@ -217,10 +232,15 @@ function step(actor::SimpleLHactor)
         ne_vol = integrate(volume_cp, cp1d.electrons.density) / volume_cp[end]
         j_parallel = actor.current_efficiency / eqt.boundary.geometric_axis.r / (ne_vol / 1e19) * power_launched
 
-        isource = resize!(cs.source, "identifier.name" => "lh_antenna_$idx")
+        name = "lh"
+        if length(actor.dd.ic_antennas.antenna) > 1
+            name = "lh_$idx"
+        end
+
+        isource = resize!(cs.source, "identifier.name" => name)
         gaussian_source_to_dd(
             isource,
-            "lh_antenna_$idx",
+            name,
             4,
             rho_cp,
             volume_cp,