L to H and H to L transition in dynamic plasma actor with ActorPedestal

src/actors/pedestal/EPED_actor.jl

@@ -13,6 +13,7 @@ Base.@kwdef mutable struct FUSEparameters__ActorEPED{T<:Real} <: ParametersActor
     T_ratio_pedestal::Entry{T} =
         Entry{T}("-", "Ratio of ion to electron temperatures (or rho at which to sample for that ratio, if negative; or rho_nml-(rho_ped-rho_nml) if 0.0)"; default=1.0)
     ped_factor::Entry{T} = Entry{T}("-", "Pedestal height multiplier (width scaled by sqrt of this factor)"; default=1.0)
+    density_factor::Entry{T} = Entry{T}("-", "Scale input density by given factor"; default=1.0)
     only_powerlaw::Entry{Bool} = Entry{Bool}("-", "EPED-NN uses power-law pedestal fit (without NN correction)"; default=false)
     #== data flow parameters ==#
     ip_from::Switch{Symbol} = switch_get_from(:ip)
@@ -62,7 +63,7 @@ function _step(actor::ActorEPED{D,P}) where {D<:Real,P<:Real}
     par = actor.par
 
     cp1d = dd.core_profiles.profiles_1d[]
-    sol = run_EPED(dd, actor.inputs, actor.epedmod; par.ne_ped_from, par.zeff_ped_from, par.βn_from, par.ip_from, par.only_powerlaw, par.warn_nn_train_bounds)
+    sol = run_EPED(dd, actor.inputs, actor.epedmod; par.ne_ped_from, par.zeff_ped_from, par.βn_from, par.ip_from, par.only_powerlaw, par.warn_nn_train_bounds, par.density_factor)
 
     if sol.pressure.GH.H < 1.1 * cp1d.pressure_thermal[end] / 1e6
         actor.pped = 1.1 * cp1d.pressure_thermal[end] / 1E6
@@ -125,11 +126,12 @@ function run_EPED(
     βn_from::Symbol,
     ip_from::Symbol,
     only_powerlaw::Bool,
-    warn_nn_train_bounds::Bool)
+    warn_nn_train_bounds::Bool,
+    density_factor::Float64)
 
     inputs = EPEDNN.InputEPED()
     epedmod = EPEDNN.loadmodelonce("EPED1NNmodel.bson")
-    return run_EPED(dd, inputs, epedmod; ne_ped_from, zeff_ped_from, βn_from, ip_from, only_powerlaw, warn_nn_train_bounds)
+    return run_EPED(dd, inputs, epedmod; ne_ped_from, zeff_ped_from, βn_from, ip_from, only_powerlaw, warn_nn_train_bounds, density_factor)
 end
 
 """
@@ -155,7 +157,8 @@ function run_EPED(
     βn_from::Symbol,
     ip_from::Symbol,
     only_powerlaw::Bool,
-    warn_nn_train_bounds::Bool)
+    warn_nn_train_bounds::Bool,
+    density_factor::Float64)
 
     cp1d = dd.core_profiles.profiles_1d[]
     eqt = dd.equilibrium.time_slice[]
@@ -165,7 +168,7 @@ function run_EPED(
         @warn "EPED-NN is only trained on m_effective = 2.0 & 2.5 , m_effective = $m"
     end
 
-    neped = IMAS.get_from(dd, Val{:ne_ped}, ne_ped_from, nothing)
+    neped = IMAS.get_from(dd, Val{:ne_ped}, ne_ped_from, nothing) * density_factor
     zeffped = IMAS.get_from(dd, Val{:zeff_ped}, zeff_ped_from, nothing)
     βn = IMAS.get_from(dd, Val{:βn}, βn_from)
     ip = IMAS.get_from(dd, Val{:ip}, ip_from)

src/actors/pedestal/WPED_actor.jl

@@ -8,6 +8,7 @@
     #== actor parameters ==#
     ped_to_core_fraction::Entry{T} = Entry{T}("-", "Ratio of edge (@rho=0.9) to core stored energy [0.05 for L-mode, 0.3 for neg-T plasmas, missing keeps original ratio]")
     rho_ped::Entry{T} = Entry{T}("-", "Defines rho at which the edge region starts")
+    density_factor::Entry{T} = Entry{T}("-", "Scale input density by given factor"; default=1.0)
     #== data flow parameters ==#
     ne_ped_from::Switch{Symbol} = switch_get_from(:ne_ped)
     zeff_ped_from::Switch{Symbol} = switch_get_from(:zeff_ped)
@@ -55,10 +56,10 @@
     par = actor.par
 
     cp1d = dd.core_profiles.profiles_1d[]
-
     summary_ped = dd.summary.local.pedestal
+
     @ddtime summary_ped.position.rho_tor_norm = par.rho_ped
-    @ddtime summary_ped.n_e.value = IMAS.get_from(dd, Val{:ne_ped}, par.ne_ped_from, par.rho_ped)
+    @ddtime summary_ped.n_e.value = IMAS.get_from(dd, Val{:ne_ped}, par.ne_ped_from, par.rho_ped) * par.density_factor
     @ddtime summary_ped.zeff.value = IMAS.get_from(dd, Val{:zeff_ped}, par.zeff_ped_from, par.rho_ped)
     IMAS.blend_core_edge(:L_mode, cp1d, summary_ped, NaN, par.rho_ped; what=:densities)
 

src/actors/pedestal/pedestal_actor.jl

@@ -11,12 +11,17 @@
     rho_nml::Entry{T} = Entry{T}("-", "Defines rho at which the no man's land region starts")
     rho_ped::Entry{T} = Entry{T}("-", "Defines rho at which the pedestal region starts") # rho_nml < rho_ped
     density_match::Switch{Symbol} = Switch{Symbol}([:ne_line, :ne_ped], "-", "Matching density based on ne_ped or line averaged density"; default=:ne_ped)
-    model::Switch{Symbol} = Switch{Symbol}([:EPED, :WPED, :auto, :none], "-", "Pedestal model to use"; default=:EPED)
+    model::Switch{Symbol} = Switch{Symbol}([:EPED, :WPED, :dynamic, :none], "-", "Pedestal model to use"; default=:EPED)
     #== data flow parameters ==#
     ip_from::Switch{Symbol} = switch_get_from(:ip)
     βn_from::Switch{Symbol} = switch_get_from(:βn)
     ne_from::Switch{Symbol} = switch_get_from(:ne_ped)
     zeff_ped_from::Switch{Symbol} = switch_get_from(:zeff_ped)
+
+    #== L to H and H to L transition model ==#
+    tau_t::Entry{T} = Entry{T}("[s]", "pedestal temperature LH transition tanh evolution time (95% of full transition)")
+    tau_n::Entry{T} = Entry{T}("[s]", "pedestal density LH transition tanh evolution time (95% of full transition)")
+    density_ratio_L_over_H::Entry{T} = Entry{T}("[-]", "n_Lmode / n_Hmode")
     #== display and debugging parameters ==#
     do_plot::Entry{Bool} = act_common_parameters(; do_plot=false)
 end
@@ -29,6 +34,10 @@
     none_actor::ActorNoOperation{D,P}
     eped_actor::ActorEPED{D,P}
     wped_actor::ActorWPED{D,P}
+    state::Vector{Symbol}
+    t_lh::Float64
+    t_hl::Float64
+    cp1d_transition::IMAS.core_profiles__profiles_1d{D}
 end
 
 """
@@ -49,7 +58,7 @@
     eped_actor = ActorEPED(dd, act.ActorEPED; ne_ped_from=par.ne_from, par.zeff_ped_from, par.βn_from, par.ip_from, par.rho_nml, par.rho_ped)
     wped_actor = ActorWPED(dd, act.ActorWPED; ne_ped_from=par.ne_from, par.zeff_ped_from, par.rho_ped, par.do_plot)
     none_actor = ActorNoOperation(dd, act.ActorNoOperation)
-    return ActorPedestal(dd, par, act, none_actor, none_actor, eped_actor, wped_actor)
+    return ActorPedestal(dd, par, act, none_actor, none_actor, eped_actor, wped_actor, Symbol[], -Inf, -Inf, IMAS.core_profiles__profiles_1d())
 end
 
 """
@@ -69,25 +78,106 @@
         actor.ped_actor = actor.none_actor
     elseif par.model == :EPED
         actor.ped_actor = actor.eped_actor
-        mode = :H_mode
+        run_selected_pedstal_model(actor)
     elseif par.model == :WPED
         actor.ped_actor = actor.wped_actor
-        mode = :L_mode
-    elseif par.model == :auto
-        if IMAS.satisfies_h_mode_conditions(dd)
-            actor.ped_actor = actor.eped_actor
-            mode = :H_mode
+        run_selected_pedstal_model(actor)
+    elseif par.model == :dynamic
+        @assert par.ne_from == :pulse_schedule "Dyanmic requires ne_from pulse schedule"
+
+        if IMAS.satisfies_h_mode_conditions(dd; threshold_multiple=2.0)
+            push!(actor.state,:H_mode)
+        elseif !IMAS.satisfies_h_mode_conditions(dd; threshold_multiple=1.0)
+            push!(actor.state,:L_mode)
         else
+            push!(actor.state,actor.state[end])
+        end
+
+        if length(actor.state) >= 2 && actor.state[end-1:end] == [:L_mode, :H_mode] 
+            # L to H
+            actor.t_lh = dd.global_time
+        elseif length(actor.state) >= 2 && actor.state[end-1:end] == [:H_mode, :L_mode] 
+            # H to L
+            actor.t_hl = dd.global_time
+
+            # The L to H and H to L model are triggered when hmode conditions are triggered so therefore pulse_schedule needs to be updated based on this model
+            if par.density_match == :ne_ped
+                @ddtime(dd.pulse_schedule.density_control.n_e_pedestal.reference =  @ddtime(dd.summary.local.pedestal.n_e.value))
+            elseif par.density_match == :ne_line
+                @ddtime(dd.pulse_schedule.density_control.n_e_line.reference =  IMAS.geometric_midplane_line_averaged_density(eqt, cp1d))
+            end
+        end
+
+        if actor.state[end] == :L_mode
+            α_t = LH_tanh_response(par.tau_t, dd.global_time,actor.t_hl;time_steps=100)
+            α_n = LH_tanh_response(par.tau_n, dd.global_time,actor.t_hl;time_steps=100)
+
             actor.ped_actor = actor.wped_actor
-            mode = :L_mode
-            if eqt.boundary.triangularity < 0.0
-                actor.wped_actor.par.ped_to_core_fraction = 0.3
-            else
-                actor.wped_actor.par.ped_to_core_fraction = 0.05
+            actor.ped_actor.par.density_factor = 1. - α_n * (1 - par.density_ratio_L_over_H)
+
+            run_selected_pedstal_model(actor)
+
+            Te_ongoing = (1 .- α_t) .* actor.cp1d_transition.electrons.temperature .+  α_t .* dd.core_profiles.profiles_1d[].electrons.temperature
+            Ti_ongoing = (1 .- α_t) .* actor.cp1d_transition.ion[1].temperature .+  α_t .* dd.core_profiles.profiles_1d[].ion[1].temperature
+
+            cp1d.electrons.temperature = Te_ongoing
+            for ion in cp1d.ion
+                ion.temperature = Ti_ongoing
+            end
+        else
+            # H mode
+            α_t = LH_tanh_response(par.tau_t,dd.global_time,actor.t_lh;time_steps=100)
+            α_n = LH_tanh_response(par.tau_n,dd.global_time, actor.t_lh;time_steps=100)
+            actor.ped_actor = actor.eped_actor
+
+            actor.ped_actor.par.density_factor =  1 + (1 / par.density_ratio_L_over_H - 1) * α_n
+            run_selected_pedstal_model(actor)
+
+            Te_ongoing = (1 .- α_t) .* actor.cp1d_transition.electrons.temperature .+  α_t .* dd.core_profiles.profiles_1d[].electrons.temperature
+            Ti_ongoing = (1 .- α_t) .* actor.cp1d_transition.ion[1].temperature .+  α_t .* dd.core_profiles.profiles_1d[].ion[1].temperature
+
+            cp1d.electrons.temperature = Te_ongoing
+            for ion in cp1d.ion
+                ion.temperature = Ti_ongoing
             end
         end
     end
 
+    return actor
+end
+
+"""
+    LH_tanh_response(τ::Float64,t_now::Float64, t_LH::Float64; time_steps::Int=100)
+
+Returns a parameter that follows a tanh like response where τ represent the value of 0.95 @ τ time starting from t_LH
+"""
+function LH_tanh_response(τ::Float64,t_now::Float64, t_LH::Float64; time_steps::Int=100)
+    if isinf(t_LH)
+        return 0.0
+    end
+    t, α = LH_tanh_response(τ, t_LH; time_steps)
+    return maximum([minimum([IMAS.interp1d(t,α).(t_now),1.0]),0.0])
+end
+
+function LH_tanh_response(τ::Float64,t_LH::Float64; time_steps::Int=100)
+    t = (0:τ/time_steps:τ) .+ t_LH
+    α = tanh.((2pi.*(t .- t_LH .- τ / 4.0)) ./ τ)
+    α = ((α .- α[1]) ./ (α[end]-α[1]))
+    return t, α
+end
+
+"""
+    run_selected_pedstal_model(actor::ActorPedestal)
+
+Runs selected pedestal model this prevents code duplication for using different par.model settings
+"""
+function run_selected_pedstal_model(actor::ActorPedestal)
+    dd = actor.dd
+    par = actor.par
+
+    eq = dd.equilibrium
+    eqt = eq.time_slice[]
+    cp1d = dd.core_profiles.profiles_1d[]
     if par.density_match == :ne_ped
         finalize(step(actor.ped_actor))
 
@@ -125,10 +215,5 @@
         # turn pressures back into expressions
         IMAS.empty!(cp1d, :pressure_thermal)
         IMAS.empty!(cp1d, :pressure)
-
-    else
-        error("act.ActorPedestal.density_match can be either one of [:ne_ped, :ne_line]")
     end
-
-    return actor
-end
+end

src/optimization.jl

@@ -32,7 +32,11 @@
 
     # create working directory
     original_dir = pwd()
-    savedir = abspath(joinpath(save_folder, "$(generation)__$(Dates.now())__$(getpid())"))
+    if Sys.iswindows()
+        savedir = abspath(joinpath(save_folder, "$(generation)__$(join(split(string(Dates.now()),":"),"-"))__$(getpid())"))        
+    else
+        savedir = abspath(joinpath(save_folder, "$(generation)__$(Dates.now())__$(getpid())"))
+    end
     mkdir(savedir)
     cd(savedir)
 

src/utils_begin.jl

@@ -324,19 +324,20 @@
 function localhost_memory()
     if Sys.isapple()
         cmd = `sysctl hw.memsize` # for OSX
-        mem_size = parse(Int, match(r"\d+", readchomp(cmd)).match) / 1024^3
+        mem_size = parse(Int, match(r"\d+", readchomp(cmd)).match) / 1024^3 # GiB
     elseif Sys.isunix()
         # General Unix command (including macOS and Linux)
         cmd = `free -b` # get memory in bytes
-        mem_size = parse(Int, match(r"\d+", readchomp(cmd)).match) / 1024^3
+        mem_size = parse(Int, match(r"\d+", readchomp(cmd)).match) / 1024^3 # GiB
     elseif Sys.iswindows()
         # Windows command
-        cmd = `wmic ComputerSystem get TotalPhysicalMemory`
-        mem_size = parse(Int, match(r"\d+", readchomp(cmd)).match) / 1024^3
+        cmd = `powershell -Command "Get-CimInstance Win32_ComputerSystem | Select-Object -ExpandProperty TotalPhysicalMemory"`
+        mem_bytes = parse(Int, readchomp(cmd))
+        mem_size = mem_bytes / 1024^3 # GiB
     elseif Sys.islinux()
         # Linux-specific command
         cmd = `grep MemTotal /proc/meminfo`
-        mem_size = parse(Int, match(r"\d+", readchomp(cmd)).match) / 1024^2 # Linux reports in KB
+        mem_size = parse(Int, match(r"\d+", readchomp(cmd)).match) / 1024^2 # GiB
     else
         error("couldn't determine the mem_size")
     end

src/utils_end.jl

@@ -808,19 +808,28 @@
         [], []
     end
 end
-
 """
     get_julia_process_memory_usage()
 
 Returns memory used by current julia process
 """
 function get_julia_process_memory_usage()
-    pid = getpid()
-    mem_info = read(`ps -p $pid -o rss=`, String)
-    mem_usage_kb = parse(Int, strip(mem_info))
+    if Sys.iswindows()
+        pid = getpid()
+        # Use PowerShell to get the current process's WorkingSet (memory in bytes)
+        cmd = `powershell -Command "(Get-Process -Id $pid).WorkingSet64"`
+        mem_bytes_str = readchomp(cmd)
+        mem_bytes = parse(Int, mem_bytes_str)
+        return mem_bytes        
+    else
+        pid = getpid()
+        mem_info = read(`ps -p $pid -o rss=`, String)
+        mem_usage_kb = parse(Int, strip(mem_info))
+    end
     return mem_usage_kb * 1024
 end
 
+
 """
     save(memtrace::MemTrace, filename::String="memtrace.txt")