Fixing transport equilibrium workflow

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
+end                                                                                                                                                                                                   

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
 

src/ddinit/init.jl

@@ -1,18 +1,26 @@
 """
     init(dd::IMAS.dd, par::Parameters; do_plot = false)
 
-Initialize all IDSs
+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))
     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 +30,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,15 +41,15 @@ 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) || :core_profiles ∈ ods_items
         init_core_profiles(dd, par)
         if do_plot
             display(plot(dd.core_profiles))
         end
     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))

src/parameters_list.jl

@@ -179,4 +179,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
+end

src/workflows/DB5_validation_workflow.jl

@@ -11,13 +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)
+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 
+    IMAS.ohmic_power_steady_state!(dd)
+
     # 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 +64,23 @@ 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)
 
     # 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)
@@ -85,21 +90,27 @@ 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
+
+    failed_df = run_df[failed_runs_ids, :]
 
     # save all input data as well as predicted tau to CSV file
     if !isempty(save_directory)
         CSV.write(joinpath(save_directory, "dataframe.csv"), run_df)
+        CSV.write(joinpath(save_directory, "failed_runs_dataframe.csv"), failed_df)
     end
 
-    failed_df = run_df[failed_runs_ids, :]
-    goof_runs_ids = [k for (k,v) in enumerate(tau_FUSE) if !isnan(v)]
-    run_df = run_df[goof_runs_ids, :]
+
 
     if plot_database
         plot_x_y_regression(run_df, "TAUTH")