Merge pull request #210 from ProjectTorreyPines/eqtransport_convergence

Replaced user iterations of EqulibriumTransport actor with convergenc…

cases/ARC.jl

@@ -9,6 +9,7 @@ function case_parameters(::Type{Val{:ARC}})::Tuple{ParametersAllInits,Parameters
     ini.general.casename = "ARC"
     ini.general.init_from = :scalars
 
+    ini.equilibrium.boundary_from = :scalars
     ini.equilibrium.R0 = 3.45
     ini.equilibrium.ϵ = 0.27 #a/R0
     ini.equilibrium.κ = 1.80 #kappa_a = 1.60, kappa_sep = 1.80

cases/CAT.jl

@@ -9,6 +9,7 @@ function case_parameters(::Type{Val{:CAT}})::Tuple{ParametersAllInits,Parameters
 
     ini.general.casename = "CAT"
     ini.general.init_from = :ods
+    ini.equilibrium.boundary_from = :ods
 
     ini.ods.filename = joinpath(@__DIR__, "..", "sample", "CAT_eq_ods.json")
 

cases/D3D.jl

@@ -9,6 +9,7 @@ function case_parameters(::Type{Val{:D3D}})::Tuple{ParametersAllInits,Parameters
 
     ini.general.casename = "D3D"
     ini.general.init_from = :ods
+    ini.equilibrium.boundary_from = :ods
 
     ini.ods.filename = joinpath(@__DIR__, "..", "sample", "D3D_eq_ods.json")
 

cases/HDB5.jl

@@ -24,6 +24,7 @@ function case_parameters(data_row::DataFrames.DataFrameRow)
     ini.general.init_from = :scalars
 
     # Equilibrium parameters
+    ini.equilibrium.boundary_from = :scalars
     ini.equilibrium.B0 = data_row[:BT]
     ini.equilibrium.R0 = data_row[:RGEO]
     ini.equilibrium.Z0 = 0.0

cases/ITER.jl

@@ -21,12 +21,16 @@ function case_parameters(::Type{Val{:ITER}}; init_from::Symbol)::Tuple{Parameter
         ini.ods.filename = joinpath(@__DIR__, "..", "sample", "ITER_eq_ods.json")
         act.ActorCXbuild.rebuild_wall = false
         act.ActorHFSsizing.fixed_aspect_ratio = true
+
+        ini.equilibrium.boundary_from = :MXH_params
+        ini.equilibrium.xpoints_number = 1
     else
         ini.equilibrium.B0 = -5.3
         ini.equilibrium.ip = 15e6
         ini.equilibrium.pressure_core = 0.643e6
 
         ini.equilibrium.boundary_from = :MXH_params
+        ini.equilibrium.xpoints_number = 1
 
         R0 = 6.2
         Z0 = 0.0
@@ -48,8 +52,6 @@ function case_parameters(::Type{Val{:ITER}}; init_from::Symbol)::Tuple{Parameter
                 asin(δ), -ζ, -0.05597, -0.01655, 0.00204, 0.00306]
         end
 
-        ini.equilibrium.xpoints_number = 1
-
         act.ActorCXbuild.rebuild_wall = true
         act.ActorHFSsizing.fixed_aspect_ratio = true
         # act.ActorEquilibrium.model = :CHEASE

cases/SPARC.jl

@@ -9,6 +9,7 @@ function case_parameters(::Type{Val{:SPARC}})::Tuple{ParametersAllInits,Paramete
     ini.general.casename = "SPARC"
     ini.general.init_from = :scalars
 
+    ini.equilibrium.boundary_from = :scalars
     ini.equilibrium.R0 = 1.85
     ini.equilibrium.ϵ = 0.308 #a/R0
     ini.equilibrium.κ = 1.97 #kappa_a = 1.75 (kappa_lower-single-null = 1.65), kappa_sep = 1.97

docs/src/develop.md

@@ -147,6 +147,10 @@ end
 
 ## Tips and more
 
+### Revise.jl
+Use [Revise.jl](https://github.com/timholy/Revise.jl) to modify code and use the changes without restarting Julia.
+We recommend adding `import Revise` to your `~/.julia/config/startup.jl`.
+
 ### Development in VSCode
 
 [VSCode](https://code.visualstudio.com) is an excellent development environment for Julia.

src/actors/chease_actor.jl

@@ -12,7 +12,7 @@ end
 Base.@kwdef struct FUSEparameters__ActorCHEASE{T} <: ParametersActor where {T<:Real}
     free_boundary = Entry(Bool, "", "Convert fixed boundary equilibrium to free boundary one"; default=true)
     clear_workdir = Entry(Bool, "", "Clean the temporary workdir for CHEASE"; default=true)
-    rescale_eq_to_ip = Entry(Bool, "", "Scale equilibrium to match Ip"; default=false)
+    rescale_eq_to_ip = Entry(Bool, "", "Scale equilibrium to match Ip"; default=true)
 end
 
 """

src/actors/compound_actors.jl

@@ -12,7 +12,8 @@ end
 
 Base.@kwdef struct FUSEparameters__ActorEquilibriumTransport{T} <: ParametersActor where {T<:Real}
     do_plot = Entry(Bool, "", "plot"; default=false)
-    iterations = Entry(Int, "", "transport-equilibrium iterations"; default=1)
+    max_iter = Entry(Int, "", "max number of transport-equilibrium iterations"; default=5)
+    convergence_error = Entry(Float64, "", "Convergence error threshold"; default=1E-2)
 end
 
 """
@@ -57,18 +58,53 @@ function _step(actor::ActorEquilibriumTransport)
     # Set j_ohmic to steady state
     finalize(step(actor.actor_jt))
 
-    for iteration in 1:par.iterations
+    # set CHEASE switches specific to this workflow
+    act_chease = deepcopy(act.ActorCHEASE)
+    actor.actor_eq.par = act_chease
+    act_chease.free_boundary = false
+    act_chease.rescale_eq_to_ip = true
+
+    iter = 1
+    total_error = 0.0
+    while (iter == 1) || (total_error > par.convergence_error)
+        # get current and pressure profiles before updating them
+        j_tor_before = dd.core_profiles.profiles_1d[].j_tor
+        pressure_before = dd.core_profiles.profiles_1d[].pressure
+
         # run transport actor
         finalize(step(actor.actor_tr))
 
+        # Set j_ohmic to steady state
+        finalize(step(actor.actor_jt))
+
         # prepare equilibrium input based on transport core_profiles output
         prepare(dd, :ActorEquilibrium, act)
 
         # run equilibrium actor with the updated beta
         finalize(step(actor.actor_eq))
 
-        # Set j_ohmic to steady state
-        finalize(step(actor.actor_jt))
+        # evaluate change in current and pressure profiles after the update
+        j_tor_after = dd.core_profiles.profiles_1d[].j_tor
+        pressure_after = dd.core_profiles.profiles_1d[].pressure
+        error_jtor = sum((j_tor_after .- j_tor_before) .^ 2) / sum(j_tor_before .^ 2)
+        error_pressure = sum((pressure_after .- pressure_before) .^ 2) / sum(pressure_before .^ 2)
+        total_error = sqrt(error_jtor + error_pressure) / 2.0
+
+        if act.ActorEquilibrium.model == :Solovev
+            total_error = 0 # temporary fix to force Solovev to run exactly once
+        end
+
+        if iter == par.max_iter
+            @warn "Max number of iterations ($max_iter) has been reached with convergence error of $(round(total_error,digits = 3)) compared to threshold of $par.convergence_error"
+            break
+        end
+        iter += 1
+
+    end
+
+    if act.ActorCHEASE.free_boundary
+        act_chease.free_boundary = true
+        finalize(step(actor.actor_eq))
     end
 
     if par.do_plot

src/ddinit/gasc.jl

@@ -101,6 +101,7 @@ end
 Convert equilibrium information in GASC solution to FUSE `ini` and `act` parameters
 """
 function gasc_2_equilibrium(gasc::GASC, ini::ParametersAllInits, act::ParametersAllActors)
+    ini.equilibrium.boundary_from = :scalars
     ini.equilibrium.B0 = gasc.inputs["conductors"]["magneticFieldOnAxis"]
     ini.equilibrium.R0 = gasc.inputs["radial build"]["majorRadius"]
     ini.equilibrium.Z0 = 0.0

src/ddinit/init_equilibrium.jl

@@ -46,6 +46,7 @@ function init_equilibrium(dd::IMAS.dd, ini::ParametersAllInits, act::ParametersA
             else
                 mxh = IMAS.MXH(pr, pz, 4)
             end
+            ini.equilibrium.xpoints_number = 0 # to use the x-points from ODS
 
         else
             eqt = resize!(dd.equilibrium.time_slice)

src/parameters_inits.jl

@@ -28,9 +28,9 @@ Base.@kwdef mutable struct FUSEparameters__equilibrium{T} <: ParametersInit wher
     xpoints_number::Entry{Integer} = Entry(Integer, "", "Number of x-points")
     ngrid::Entry{Int} = Entry(Int, "", "Resolution of the equilibrium grid"; default=129)
     field_null_surface::Entry{T} = Entry(T, "", "ψn value of the field_null_surface. Disable with 0.0"; default=0.5)
-    boundary_from::Switch{Symbol} = Switch(Symbol, [:scalars, :MXH_params, :rz_points, :ods], "", "The starting r, z boundary taken from"; default=:scalars)
-    MXH_params::Entry{Union{Nothing,Vector{<:T}}} = Entry(Union{Nothing,Vector{<:T}}, "", "Vector of MXH flats", default=missing)
-    rz_points::Entry{Union{Nothing,Vector{Vector{<:T}}}} = Entry(Union{Nothing,Vector{Vector{<:T}}}, "m", "R_Z boundary as Vector{Vector{<:Real}}} : r = rz_points[1], z = rz_points[2]", default=missing)
+    boundary_from::Switch{Symbol} = Switch(Symbol, [:scalars, :MXH_params, :rz_points, :ods], "", "The starting r, z boundary taken from")
+    MXH_params::Entry{Vector{<:T}} = Entry(Vector{<:T}, "", "Vector of MXH flats")
+    rz_points::Entry{Vector{Vector{<:T}}} = Entry(Vector{Vector{<:T}}, "m", "R_Z boundary as Vector{Vector{<:Real}}} : r = rz_points[1], z = rz_points[2]")
 end
 
 Base.@kwdef mutable struct FUSEparameters__core_profiles{T} <: ParametersInit where {T<:Real}