From faae82599fe719150b53d37b6cf42cf9adf2ba94 Mon Sep 17 00:00:00 2001 From: jmcclena Date: Thu, 12 Dec 2024 10:46:35 -0800 Subject: [PATCH 1/4] initial setup for raising q>1 for stationary transport actor --- src/actors/current/current_actor.jl | 16 +++++++------- src/actors/current/qed_actor.jl | 34 +++++++++++++++++++++++++++-- 2 files changed, 40 insertions(+), 10 deletions(-) diff --git a/src/actors/current/current_actor.jl b/src/actors/current/current_actor.jl index 45ac1839d..21943252c 100644 --- a/src/actors/current/current_actor.jl +++ b/src/actors/current/current_actor.jl @@ -5,7 +5,7 @@ Base.@kwdef mutable struct FUSEparameters__ActorCurrent{T<:Real} <: ParametersAc _parent::WeakRef = WeakRef(nothing) _name::Symbol = :not_set _time::Float64 = NaN - model::Switch{Symbol} = Switch{Symbol}([:SteadyStateCurrent, :QED, :none], "-", "Current actor to run"; default=:SteadyStateCurrent) + model::Switch{Symbol} = Switch{Symbol}([:SteadyStateCurrent, :QED, :none], "-", "Current actor to run"; default=:QED) allow_floating_plasma_current::Entry{Bool} = Entry{Bool}("-", "Zero loop voltage if non-inductive fraction exceeds 100% of the target Ip"; default=true) #== data flow parameters ==# ip_from::Switch{Symbol} = switch_get_from(:ip) @@ -57,9 +57,9 @@ function _step(actor::ActorCurrent) # freeze jboot and j_non_inductive before updating johmic cp1d = dd.core_profiles.profiles_1d[] - for field in [:j_bootstrap, :j_non_inductive] - IMAS.refreeze!(cp1d, field) - end + #for field in [:j_bootstrap, :j_non_inductive] + # IMAS.refreeze!(cp1d, field) + #end step(actor.jt_actor) @@ -79,13 +79,13 @@ function _finalize(actor::ActorCurrent) # freeze cp1d j_total and j_tor after johmic update # important to freeze first j_total and then j_tor cp1d = dd.core_profiles.profiles_1d[] - for field in (:j_total, :j_tor) - IMAS.refreeze!(cp1d, field) - end + #for field in (:j_total, :j_tor) + # IMAS.refreeze!(cp1d, field) + #end # similarly, freeze parallel plasma current eqt = dd.equilibrium.time_slice[] - IMAS.refreeze!(eqt.profiles_1d, :j_parallel) + #IMAS.refreeze!(eqt.profiles_1d, :j_parallel) # update core_sources related to current IMAS.bootstrap_source!(dd) diff --git a/src/actors/current/qed_actor.jl b/src/actors/current/qed_actor.jl index 3b04d9887..9fb7e5965 100644 --- a/src/actors/current/qed_actor.jl +++ b/src/actors/current/qed_actor.jl @@ -120,8 +120,33 @@ function _step(actor::ActorQED) Ip = nothing Vedge = IMAS.get_from(dd, Val{:vloop}, par.vloop_from) end + B0 = eqt.global_quantities.vacuum_toroidal_field.b0 + actor.QO = QED.steady_state(actor.QO, η_imas(dd.core_profiles.profiles_1d[]); Vedge, Ip) + + j_total = QED.JB(actor.QO; ρ=cp1d.grid.rho_tor_norm) ./ B0 + p2 = plot(j_total) + display(p2) + if true + rho = cp1d.grid.rho_tor_norm + qval = 1.0 ./ abs.(actor.QO.ι.(rho)) + #println(qval) + qdes = 1.0 + i_qeq1 = findlast(qval .< qdes) + print(i_qeq1) + #η = η_imas(dd.core_profiles.profiles_1d[],i_qeq1) + actor.QO = QED.steady_state(actor.QO, η_imas(dd.core_profiles.profiles_1d[],i_qeq1=i_qeq1); Vedge, Ip) + #p = plot(qval) + #display(p) + #p2 = plot!(1.0./abs.(actor.QO.ι.(rho))) + #display(p2) + + j_total = QED.JB(actor.QO; ρ=cp1d.grid.rho_tor_norm) ./ B0 + p2 = plot!(j_total) + display(p2) + #print(qval-1.0./abs.(actor.QO.ι.(rho)),qval,1.0./abs.(actor.QO.ι.(rho))) + end end return actor @@ -135,7 +160,8 @@ function _finalize(actor::ActorQED) cp1d = dd.core_profiles.profiles_1d[] j_total = QED.JB(actor.QO; ρ=cp1d.grid.rho_tor_norm) ./ B0 - + p = plot!(j_total) + display(p) if ismissing(cp1d, :j_non_inductive) cp1d.j_ohmic = j_total else @@ -191,8 +217,12 @@ function qed_init_from_imas(eqt::IMAS.equilibrium__time_slice, cp1d::IMAS.core_p return QED.initialize(rho_tor, B0, gm1, f, dvolume_drho_tor, q, j_tor, gm9; ρ_j_non_inductive, ρ_grid) end -function η_imas(cp1d::IMAS.core_profiles__profiles_1d; use_log::Bool=true) +function η_imas(cp1d::IMAS.core_profiles__profiles_1d; use_log::Bool=true,i_qeq1=nothing) rho = cp1d.grid.rho_tor_norm η = 1.0 ./ cp1d.conductivity_parallel + if i_qeq1 != nothing + η[1:i_qeq1] .= η[i_qeq1] + #println(i_qeq1,η) + end return QED.η_FE(rho, η; use_log) end From ffc8ee9303fa7fffdd9f77036548e8588c9409f4 Mon Sep 17 00:00:00 2001 From: jmcclena Date: Thu, 12 Dec 2024 14:08:26 -0800 Subject: [PATCH 2/4] add back in refreezes --- src/actors/current/current_actor.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/actors/current/current_actor.jl b/src/actors/current/current_actor.jl index 21943252c..e3eb7b9c1 100644 --- a/src/actors/current/current_actor.jl +++ b/src/actors/current/current_actor.jl @@ -57,9 +57,9 @@ function _step(actor::ActorCurrent) # freeze jboot and j_non_inductive before updating johmic cp1d = dd.core_profiles.profiles_1d[] - #for field in [:j_bootstrap, :j_non_inductive] - # IMAS.refreeze!(cp1d, field) - #end + for field in [:j_bootstrap, :j_non_inductive] + IMAS.refreeze!(cp1d, field) + end step(actor.jt_actor) @@ -79,13 +79,13 @@ function _finalize(actor::ActorCurrent) # freeze cp1d j_total and j_tor after johmic update # important to freeze first j_total and then j_tor cp1d = dd.core_profiles.profiles_1d[] - #for field in (:j_total, :j_tor) - # IMAS.refreeze!(cp1d, field) - #end + for field in (:j_total, :j_tor) + IMAS.refreeze!(cp1d, field) + end # similarly, freeze parallel plasma current eqt = dd.equilibrium.time_slice[] - #IMAS.refreeze!(eqt.profiles_1d, :j_parallel) + IMAS.refreeze!(eqt.profiles_1d, :j_parallel) # update core_sources related to current IMAS.bootstrap_source!(dd) From c7086054d991d9c2c72bdc6502bd11abaa1ca965 Mon Sep 17 00:00:00 2001 From: jmcclena Date: Thu, 12 Dec 2024 15:38:24 -0800 Subject: [PATCH 3/4] clean cruft --- src/actors/current/qed_actor.jl | 31 +++++++++---------------------- 1 file changed, 9 insertions(+), 22 deletions(-) diff --git a/src/actors/current/qed_actor.jl b/src/actors/current/qed_actor.jl index 9fb7e5965..22918566c 100644 --- a/src/actors/current/qed_actor.jl +++ b/src/actors/current/qed_actor.jl @@ -14,6 +14,7 @@ Base.@kwdef mutable struct FUSEparameters__ActorQED{T<:Real} <: ParametersActor{ #== data flow parameters ==# ip_from::Switch{Symbol} = switch_get_from(:ip) vloop_from::Switch{Symbol} = switch_get_from(:vloop) + q0_desired::Entry{Float64} = Entry{Float64}("-", "Desired minimum q-profile"; default=1.0) end mutable struct ActorQED{D,P} <: SingleAbstractActor{D,P} @@ -126,26 +127,14 @@ function _step(actor::ActorQED) actor.QO = QED.steady_state(actor.QO, η_imas(dd.core_profiles.profiles_1d[]); Vedge, Ip) j_total = QED.JB(actor.QO; ρ=cp1d.grid.rho_tor_norm) ./ B0 - p2 = plot(j_total) - display(p2) - if true + + if par.q0_desired > 0 rho = cp1d.grid.rho_tor_norm qval = 1.0 ./ abs.(actor.QO.ι.(rho)) - #println(qval) - qdes = 1.0 - i_qeq1 = findlast(qval .< qdes) - print(i_qeq1) - #η = η_imas(dd.core_profiles.profiles_1d[],i_qeq1) - actor.QO = QED.steady_state(actor.QO, η_imas(dd.core_profiles.profiles_1d[],i_qeq1=i_qeq1); Vedge, Ip) - #p = plot(qval) - #display(p) - #p2 = plot!(1.0./abs.(actor.QO.ι.(rho))) - #display(p2) + i_qdes = findlast(qval .< par.q0_desired) + actor.QO = QED.steady_state(actor.QO, η_imas(dd.core_profiles.profiles_1d[],i_qdes=i_qdes); Vedge, Ip) j_total = QED.JB(actor.QO; ρ=cp1d.grid.rho_tor_norm) ./ B0 - p2 = plot!(j_total) - display(p2) - #print(qval-1.0./abs.(actor.QO.ι.(rho)),qval,1.0./abs.(actor.QO.ι.(rho))) end end @@ -160,8 +149,7 @@ function _finalize(actor::ActorQED) cp1d = dd.core_profiles.profiles_1d[] j_total = QED.JB(actor.QO; ρ=cp1d.grid.rho_tor_norm) ./ B0 - p = plot!(j_total) - display(p) + if ismissing(cp1d, :j_non_inductive) cp1d.j_ohmic = j_total else @@ -217,12 +205,11 @@ function qed_init_from_imas(eqt::IMAS.equilibrium__time_slice, cp1d::IMAS.core_p return QED.initialize(rho_tor, B0, gm1, f, dvolume_drho_tor, q, j_tor, gm9; ρ_j_non_inductive, ρ_grid) end -function η_imas(cp1d::IMAS.core_profiles__profiles_1d; use_log::Bool=true,i_qeq1=nothing) +function η_imas(cp1d::IMAS.core_profiles__profiles_1d; use_log::Bool=true,i_qdes=nothing) rho = cp1d.grid.rho_tor_norm η = 1.0 ./ cp1d.conductivity_parallel - if i_qeq1 != nothing - η[1:i_qeq1] .= η[i_qeq1] - #println(i_qeq1,η) + if i_qdes != nothing + η[1:i_qdes] .= η[i_qdes] end return QED.η_FE(rho, η; use_log) end From 5c5bd96cf6ac9f9eb9786117f9f613d2d9dfae65 Mon Sep 17 00:00:00 2001 From: Orso Meneghini Date: Mon, 16 Dec 2024 12:07:30 -0800 Subject: [PATCH 4/4] small updates to qed jardin sawteeth --- src/actors/current/qed_actor.jl | 22 ++++++++++------------ 1 file changed, 10 insertions(+), 12 deletions(-) diff --git a/src/actors/current/qed_actor.jl b/src/actors/current/qed_actor.jl index 22918566c..90fe50d1c 100644 --- a/src/actors/current/qed_actor.jl +++ b/src/actors/current/qed_actor.jl @@ -14,7 +14,7 @@ Base.@kwdef mutable struct FUSEparameters__ActorQED{T<:Real} <: ParametersActor{ #== data flow parameters ==# ip_from::Switch{Symbol} = switch_get_from(:ip) vloop_from::Switch{Symbol} = switch_get_from(:vloop) - q0_desired::Entry{Float64} = Entry{Float64}("-", "Desired minimum q-profile"; default=1.0) + qmin_desired::Entry{Float64} = Entry{Float64}("-", "Desired minimum magnitude of q-profile"; default=1.0) end mutable struct ActorQED{D,P} <: SingleAbstractActor{D,P} @@ -123,17 +123,16 @@ function _step(actor::ActorQED) end B0 = eqt.global_quantities.vacuum_toroidal_field.b0 - actor.QO = QED.steady_state(actor.QO, η_imas(dd.core_profiles.profiles_1d[]); Vedge, Ip) - j_total = QED.JB(actor.QO; ρ=cp1d.grid.rho_tor_norm) ./ B0 - if par.q0_desired > 0 - rho = cp1d.grid.rho_tor_norm - qval = 1.0 ./ abs.(actor.QO.ι.(rho)) - i_qdes = findlast(qval .< par.q0_desired) - actor.QO = QED.steady_state(actor.QO, η_imas(dd.core_profiles.profiles_1d[],i_qdes=i_qdes); Vedge, Ip) + # using Jardin's model for stationary sawteeth to raise q>1 + # requires running QED twice (the first time to evaluate q, the second to relax using updated q profile) + if par.qmin_desired > 0 + qval = 1.0 ./ abs.(actor.QO.ι.(cp1d.grid.rho_tor_norm)) + i_qdes = findlast(qval .< par.qmin_desired) + actor.QO = QED.steady_state(actor.QO, η_imas(dd.core_profiles.profiles_1d[]; i_qdes); Vedge, Ip) j_total = QED.JB(actor.QO; ρ=cp1d.grid.rho_tor_norm) ./ B0 end end @@ -205,11 +204,10 @@ function qed_init_from_imas(eqt::IMAS.equilibrium__time_slice, cp1d::IMAS.core_p return QED.initialize(rho_tor, B0, gm1, f, dvolume_drho_tor, q, j_tor, gm9; ρ_j_non_inductive, ρ_grid) end -function η_imas(cp1d::IMAS.core_profiles__profiles_1d; use_log::Bool=true,i_qdes=nothing) - rho = cp1d.grid.rho_tor_norm +function η_imas(cp1d::IMAS.core_profiles__profiles_1d; use_log::Bool=true, i_qdes::Int=0) η = 1.0 ./ cp1d.conductivity_parallel - if i_qdes != nothing + if i_qdes != 0 η[1:i_qdes] .= η[i_qdes] end - return QED.η_FE(rho, η; use_log) + return QED.η_FE(cp1d.grid.rho_tor_norm, η; use_log) end