From 48f4c579074f1d6b4b07bc1e250653714fe7315e Mon Sep 17 00:00:00 2001 From: Orso Meneghini Date: Thu, 12 Sep 2024 22:10:37 -0700 Subject: [PATCH 01/19] FRESCO instead of Solovev --- Makefile | 21 +-- Project.toml | 1 + src/FUSE.jl | 2 +- src/actors/equilibrium/equilibrium_actor.jl | 10 +- src/actors/equilibrium/fresco_actor.jl | 118 ++++++++++++ src/actors/equilibrium/solovev_actor.jl | 199 -------------------- 6 files changed, 130 insertions(+), 221 deletions(-) create mode 100644 src/actors/equilibrium/fresco_actor.jl delete mode 100644 src/actors/equilibrium/solovev_actor.jl diff --git a/Makefile b/Makefile index 09cf12982..64e6f8cb0 100644 --- a/Makefile +++ b/Makefile @@ -15,7 +15,7 @@ else FUSE_LOCAL_BRANCH=$(shell echo $(GITHUB_REF) | sed 's/refs\/heads\///') endif -FUSE_PACKAGES_MAKEFILE := ADAS BoundaryPlasmaModels CHEASE CoordinateConventions EPEDNN FiniteElementHermite Fortran90Namelists FuseUtils FusionMaterials FuseExchangeProtocol IMAS IMASdd MXHEquilibrium MeshTools MillerExtendedHarmonic NEO NNeutronics QED RABBIT SimulationParameters TEQUILA TGLFNN TJLF VacuumFields XSteam ThermalSystemModels +FUSE_PACKAGES_MAKEFILE := ADAS BoundaryPlasmaModels CHEASE CoordinateConventions EPEDNN FiniteElementHermite Fortran90Namelists FRESCO FuseUtils FusionMaterials FuseExchangeProtocol IMAS IMASdd MXHEquilibrium MeshTools MillerExtendedHarmonic NEO NNeutronics QED RABBIT SimulationParameters TEQUILA TGLFNN TJLF VacuumFields XSteam ThermalSystemModels FUSE_PACKAGES_MAKEFILE := $(sort $(FUSE_PACKAGES_MAKEFILE)) FUSE_PACKAGES := $(shell echo '$(FUSE_PACKAGES_MAKEFILE)' | awk '{printf("[\"%s\"", $$1); for (i=2; i<=NF; i++) printf(", \"%s\"", $$i); print "]"}') DEV_PACKAGES := $(shell find ../*/.git/config -exec grep ProjectTorreyPines \{\} /dev/null \; | cut -d'/' -f 2) @@ -113,6 +113,9 @@ FiniteElementHermite: Fortran90Namelists: $(call clone_pull_repo,$@) +FRESCO: + $(call clone_pull_repo,$@) + CHEASE: $(call clone_pull_repo,$@) @@ -158,13 +161,6 @@ help: header help_info # =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%= # -# @# @")), "# @"); \ - open(file*"", "w") do f \ - write(f, strip(header)); \ - write(f, "\n\n# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%= #\n\n"); \ - write(f, targets); \ - end;\ - '# @"*join(sort!(split(strip(targets), "# @devs blank_examples:clean_examples # Clean and convert examples to md without executing cd docs; julia notebooks_to_md.jl @@ -454,13 +450,6 @@ header: @printf " \033[1;30m╚═╝ ╚═════╝ ╚══════╝╚══════╝\033[0m\n" @printf " Project Torrey Pines (PTP)\n" -# Update Makefile targets -sort_targets: - @julia -e 'using DelimitedFiles; \ - file = "Makefile"; \ - lines = readlines(file); \ - header, targets = split(join(lines, "\n"), "# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%= #"; limit=2); \ - targets = "# @devs init_expressions: # Generates init_expressions.json file, which lists entries that are # always expected to be expressions when coming out of init() @@ -546,7 +535,7 @@ list_open_compats: ) # @devs -make_apache: error_missing_repo_var +apache: error_missing_repo_var # Update LICENSE, NOTICE.md, github workflows, docs, juliaformatter and gitignore in preparation of public release # The starting information is taken from IMASdd.jl and moved to the target repo # >> make apache repo=CHEASE diff --git a/Project.toml b/Project.toml index e11f3f924..2658571fc 100644 --- a/Project.toml +++ b/Project.toml @@ -16,6 +16,7 @@ DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b" EPEDNN = "e64856f0-3bb8-4376-b4b7-c03396503991" FFMPEG = "c87230d0-a227-11e9-1b43-d7ebe4e7570a" +FRESCO = "f2db6eb9-2c54-4c18-912b-ce467796777a" FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" FuseUtils = "c8e4ee58-0bec-4a6a-b552-9dd24d6374ac" diff --git a/src/FUSE.jl b/src/FUSE.jl index 24082c24a..b515d654b 100644 --- a/src/FUSE.jl +++ b/src/FUSE.jl @@ -71,7 +71,7 @@ include("actors.jl") include(joinpath("actors", "noop_actor.jl")) -include(joinpath("actors", "equilibrium", "solovev_actor.jl")) +include(joinpath("actors", "equilibrium", "fresco_actor.jl")) include(joinpath("actors", "equilibrium", "chease_actor.jl")) include(joinpath("actors", "equilibrium", "tequila_actor.jl")) include(joinpath("actors", "equilibrium", "equilibrium_actor.jl")) diff --git a/src/actors/equilibrium/equilibrium_actor.jl b/src/actors/equilibrium/equilibrium_actor.jl index db72df608..280c2a203 100644 --- a/src/actors/equilibrium/equilibrium_actor.jl +++ b/src/actors/equilibrium/equilibrium_actor.jl @@ -6,7 +6,7 @@ Base.@kwdef mutable struct FUSEparameters__ActorEquilibrium{T<:Real} <: Paramete _name::Symbol = :not_set _time::Float64 = NaN #== actor parameters ==# - model::Switch{Symbol} = Switch{Symbol}([:Solovev, :CHEASE, :TEQUILA], "-", "Equilibrium actor to run"; default=:TEQUILA) + model::Switch{Symbol} = Switch{Symbol}([:FRESCO, :CHEASE, :TEQUILA], "-", "Equilibrium actor to run"; default=:TEQUILA) symmetrize::Entry{Bool} = Entry{Bool}("-", "Force equilibrium up-down symmetry with respect to magnetic axis"; default=false) #== data flow parameters ==# j_p_from::Switch{Symbol} = Switch{Symbol}([:equilibrium, :core_profiles], "-", "Take j_tor and pressure profiles from this IDS"; default=:core_profiles) @@ -20,7 +20,7 @@ mutable struct ActorEquilibrium{D,P} <: CompoundAbstractActor{D,P} dd::IMAS.dd{D} par::FUSEparameters__ActorEquilibrium{P} act::ParametersAllActors - eq_actor::Union{Nothing,ActorSolovev{D,P},ActorCHEASE{D,P},ActorTEQUILA{D,P}} + eq_actor::Union{Nothing,ActorFRESCO{D,P},ActorCHEASE{D,P},ActorTEQUILA{D,P}} end """ @@ -38,14 +38,14 @@ end function ActorEquilibrium(dd::IMAS.dd, par::FUSEparameters__ActorEquilibrium, act::ParametersAllActors; kw...) logging_actor_init(ActorEquilibrium) par = par(kw...) - if par.model == :Solovev - eq_actor = ActorSolovev(dd, act.ActorSolovev; par.ip_from) + if par.model == :FRESCO + eq_actor = ActorFRESCO(dd, act.ActorFRESCO; par.ip_from) elseif par.model == :CHEASE eq_actor = ActorCHEASE(dd, act.ActorCHEASE; par.ip_from) elseif par.model == :TEQUILA eq_actor = ActorTEQUILA(dd, act.ActorTEQUILA; par.ip_from) else - error("ActorEquilibrium: model = `$(par.model)` can only be `:Solovev` or `:CHEASE`") + error("ActorEquilibrium: model = `$(par.model)` can only be `:TEQUILA`, `:FRESCO`, `:CHEASE`") end return ActorEquilibrium(dd, par, act, eq_actor) end diff --git a/src/actors/equilibrium/fresco_actor.jl b/src/actors/equilibrium/fresco_actor.jl new file mode 100644 index 000000000..3d902bb48 --- /dev/null +++ b/src/actors/equilibrium/fresco_actor.jl @@ -0,0 +1,118 @@ +import FRESCO + +#= =========== =# +# ActorFRESCO # +#= =========== =# +Base.@kwdef mutable struct FUSEparameters__ActorFRESCO{T<:Real} <: ParametersActor{T} + _parent::WeakRef = WeakRef(nothing) + _name::Symbol = :not_set + _time::Float64 = NaN + #== actor parameters ==# + control::Switch{Symbol} = Switch{Symbol}([:vertical, :shape], "-", ""; default=:shape) + relax::Entry{Float64} = Entry{Float64}("-", "Relaxation on the Picard iterations"; default=0.5) + #tolerance::Entry{Float64} = Entry{Float64}("-", "Tolerance for terminating iterations"; default=1e-4) + #== data flow parameters ==# + ip_from::Switch{Symbol} = switch_get_from(:ip) + # fixed_grid::Switch{Symbol} = Switch{Symbol}([:poloidal, :toroidal], "-", "Fix P and Jt on this rho grid"; default=:toroidal) + #== display and debugging parameters ==# + do_plot::Entry{Bool} = act_common_parameters(; do_plot=false) + debug::Entry{Bool} = Entry{Bool}("-", "Print debug information withing FRESCO solve"; default=false) + #== IMAS psi grid settings ==# + nR::Entry{Int} = Entry{Int}("-", "Grid resolution along R") + nZ::Entry{Int} = Entry{Int}("-", "Grid resolution along Z") +end + +mutable struct ActorFRESCO{D,P} <: SingleAbstractActor{D,P} + dd::IMAS.dd{D} + par::FUSEparameters__ActorFRESCO{P} + canvas::Union{Nothing,FRESCO.Canvas} +end + +""" + ActorFRESCO(dd::IMAS.dd, act::ParametersAllActors; kw...) + +Runs the Fixed boundary equilibrium solver FRESCO +""" +function ActorFRESCO(dd::IMAS.dd, act::ParametersAllActors; kw...) + actor = ActorFRESCO(dd, act.ActorFRESCO; kw...) + step(actor) + finalize(actor) + return actor +end + +function ActorFRESCO(dd::IMAS.dd{D}, par::FUSEparameters__ActorFRESCO{P}; kw...) where {D<:Real,P<:Real} + logging_actor_init(ActorFRESCO) + par = par(kw...) + return ActorFRESCO(dd, par, nothing) +end + +""" + _step(actor::ActorFRESCO) + +Runs FRESCO on the r_z boundary, equilibrium pressure and equilibrium j_tor +""" +function _step(actor::ActorFRESCO) + dd = actor.dd + par = actor.par + eqt = dd.equilibrium.time_slice[] + eq1d = eqt.profiles_1d + + # Ip_target = eqt.global_quantities.ip + # + # if par.fixed_grid === :poloidal + # rho = sqrt.(eq1d.psi_norm) + # rho[1] = 0.0 + # P = (FRESCO.FE(rho, eq1d.pressure), :poloidal) + # # don't allow current to change sign + # Jt = (FRESCO.FE(rho, [sign(j) == sign(Ip_target) ? j : 0.0 for j in eq1d.j_tor]), :poloidal) + # Pbnd = eq1d.pressure[end] + # elseif par.fixed_grid === :toroidal + # rho = eq1d.rho_tor_norm + # P = (FRESCO.FE(rho, eq1d.pressure), :toroidal) + # # don't allow current to change sign + # Jt = (FRESCO.FE(rho, [sign(j) == sign(Ip_target) ? j : 0.0 for j in eq1d.j_tor]), :toroidal) + # Pbnd = eq1d.pressure[end] + # end + + psi_norm = eq1d.psi_norm + gpp = IMAS.interp1d(psi_norm, eq1d.dpressure_dpsi, :cubic) + pprime = x -> gpp(x) + gffp = IMAS.interp1d(psi_norm, eq1d.f_df_dpsi, :cubic) + ffprime = x -> gffp(x) + profile = FRESCO.PprimeFFprime(pprime, ffprime) + + actor.canvas = FRESCO.Canvas(dd, par.nR, par.nZ); + + FRESCO.solve!(actor.canvas, profile, 100, 3; relax=par.relax, par.debug, par.control) # forward mode + + # using Plots + # p1 = Plots.heatmap(Rs, Zs, psi', aspect_ratio=:equal, xrange=(0,12.5), yrange=(-9,9), size=(400,500)) + # Plots.plot!(p1, coils, label=nothing) + # Plots.contour!(p1, Rs, Zs, psi', levels=[C.Ψbnd], color=:white) + # display(p1) + + return actor +end + +# finalize by converting FRESCO canvas to dd.equilibrium +function _finalize(actor::ActorFRESCO) + canvas = actor.canvas + dd = actor.dd + eq = dd.equilibrium + eqt = eq.time_slice[] + eq1d = eqt.profiles_1d + eq2d = resize!(eqt.profiles_2d, 1)[1] + + Raxis, Zaxis, _ = FRESCO.find_axis(canvas) + eqt.global_quantities.magnetic_axis.r = Raxis + eqt.global_quantities.magnetic_axis.z = Zaxis + eq1d.psi = range(canvas.Ψaxis, canvas.Ψbnd, length(eq1d.psi)) + # p, p', f, ff' don't change + + eq2d.grid_type.index = 1 + eq2d.grid.dim1 = canvas.Rs + eq2d.grid.dim2 = canvas.Zs + eq2d.psi = canvas.Ψ + + return actor +end diff --git a/src/actors/equilibrium/solovev_actor.jl b/src/actors/equilibrium/solovev_actor.jl deleted file mode 100644 index 8ce4d0d58..000000000 --- a/src/actors/equilibrium/solovev_actor.jl +++ /dev/null @@ -1,199 +0,0 @@ -import MXHEquilibrium -import Optim - -#= ============ =# -# ActorSolovev # -#= ============ =# -Base.@kwdef mutable struct FUSEparameters__ActorSolovev{T<:Real} <: ParametersActor{T} - _parent::WeakRef = WeakRef(nothing) - _name::Symbol = :not_set - _time::Float64 = NaN - #== actor parameters ==# - ngrid::Entry{Int} = Entry{Int}("-", "Grid size (for R, Z follows proportionally to plasma elongation)"; default=129) - qstar::Entry{T} = Entry{T}("-", "Initial guess of kink safety factor"; default=1.5) - alpha::Entry{T} = Entry{T}("-", "Initial guess of constant relating to pressure"; default=0.0) - #== data flow parameters ==# - ip_from::Switch{Symbol} = switch_get_from(:ip) - #== display and debugging parameters ==# - verbose::Entry{Bool} = act_common_parameters(; verbose=false) -end - -mutable struct ActorSolovev{D,P} <: SingleAbstractActor{D,P} - dd::IMAS.dd{D} - par::FUSEparameters__ActorSolovev{P} - S::Union{Nothing,MXHEquilibrium.SolovevEquilibrium} -end - -""" - ActorSolovev(dd::IMAS.dd, act::ParametersAllActors; kw...) - -Solovev equilibrium actor, based on: -“One size fits all” analytic solutions to the Grad–Shafranov equation -Phys. Plasmas 17, 032502 (2010); https://doi.org/10.1063/1.3328818 -""" -function ActorSolovev(dd::IMAS.dd, act::ParametersAllActors; kw...) - actor = ActorSolovev(dd, act.ActorSolovev; kw...) - step(actor) - finalize(actor) - return actor -end - -function ActorSolovev(dd::IMAS.dd, par::FUSEparameters__ActorSolovev; kw...) - logging_actor_init(ActorSolovev) - par = par(kw...) - return ActorSolovev(dd, par, nothing) -end - -""" - _step(actor::ActorSolovev) - -Non-linear optimization to obtain a target `ip` and `pressure_core` -""" -function _step(actor::ActorSolovev) - dd = actor.dd - par = actor.par - - # initialize eqt from pulse_schedule and core_profiles - eq = dd.equilibrium - eqt = eq.time_slice[] - - # magnetic field - B0 = @ddtime eq.vacuum_toroidal_field.b0 - target_ip = eqt.global_quantities.ip - target_pressure_core = eqt.profiles_1d.pressure[1] - - # plasma shape as MXH - mxh = IMAS.MXH(eqt.boundary.outline.r, eqt.boundary.outline.z, 4) - Z0off = mxh.Z0 # Solovev has a bug for Z!=0.0 - mxh.Z0 -= Z0off - plasma_shape = MXHEquilibrium.MillerExtendedHarmonicShape(mxh.R0, mxh.Z0, mxh.ϵ, mxh.κ, mxh.c0, mxh.c, mxh.s) - # check number of x_points to infer symmetry - if mod(length(eqt.boundary.x_point), 2) == 0 - symmetric = true - else - symmetric = false - end - - # add x_point info - if length(eqt.boundary.x_point) > 0 - x_point = (eqt.boundary.x_point[1].r, -abs(eqt.boundary.x_point[1].z) - Z0off) - else - x_point = nothing - end - - # first run of Solovev - # display(plot(mxh)) - # @show abs(B0) - # @show plasma_shape - # @show par.alpha - # @show par.qstar - # @show x_point - # @show symmetric - S0 = MXHEquilibrium.solovev(abs(B0), plasma_shape, par.alpha, par.qstar; B0_dir=Int64(sign(B0)), Ip_dir=1, x_point, symmetric) - - # optimize `alpha` and `qstar` to get target_ip and target_pressure_core - function cost(x) - S = MXHEquilibrium.solovev(abs(B0), plasma_shape, x[1], x[2]; B0_dir=Int64(sign(B0)), Ip_dir=1, symmetric, x_point) - psimag, psibry = MXHEquilibrium.psi_limits(S) - pressure_cost = (MXHEquilibrium.pressure(S, psimag) - target_pressure_core) / target_pressure_core - ip_cost = (MXHEquilibrium.plasma_current(S) - target_ip) / target_ip - return sqrt(pressure_cost^2 + 100.0 * ip_cost^2) - end - res = Optim.optimize(cost, [S0.alpha, S0.qstar], Optim.NelderMead(), Optim.Options(; g_tol=1E-3)) - if par.verbose - println(res) - end - - # final Solovev solution - actor.S = MXHEquilibrium.solovev(abs(B0), plasma_shape, res.minimizer[1], res.minimizer[2]; B0_dir=Int64(sign(B0)), Ip_dir=1, symmetric, x_point) - - return actor -end - -""" - _finalize(actor::ActorSolovev) - -Store ActorSolovev data in IMAS.equilibrium format -""" -function _finalize(actor::ActorSolovev) - dd = actor.dd - par = actor.par - mxh_eq = actor.S - - eqt = dd.equilibrium.time_slice[] - - target_ip = eqt.global_quantities.ip - target_psi_norm = getproperty(eqt.profiles_1d, :psi_norm, missing) - target_pressure = getproperty(eqt.profiles_1d, :pressure, missing) - target_j_tor = getproperty(eqt.profiles_1d, :j_tor, missing) - - MXHEquilibrium_to_dd!(dd.equilibrium, mxh_eq, par.ngrid; cocos_in=3) - - # force total plasma current to target_ip to avoid drifting after multiple calls of SolovevActor - eqt2d = findfirst(:rectangular, eqt.profiles_2d) - eqt2d.psi = (eqt2d.psi .- eqt.profiles_1d.psi[end]) .* (target_ip / eqt.global_quantities.ip) .+ eqt.profiles_1d.psi[end] - eqt.profiles_1d.psi = (eqt.profiles_1d.psi .- eqt.profiles_1d.psi[end]) .* (target_ip / eqt.global_quantities.ip) .+ eqt.profiles_1d.psi[end] - # match entry target_pressure and target_j_tor as if Solovev could do this - if !ismissing(target_pressure) - eqt.profiles_1d.pressure = IMAS.interp1d(target_psi_norm, target_pressure).(eqt.profiles_1d.psi_norm) - end - if !ismissing(target_j_tor) - eqt.profiles_1d.j_tor = IMAS.interp1d(target_psi_norm, target_j_tor, :cubic).(eqt.profiles_1d.psi_norm) - end - IMAS.p_jtor_2_pprime_ffprim_f!(eqt.profiles_1d, mxh_eq.S.R0, mxh_eq.B0) - - # record optimized values of qstar and alpha in `act` for subsequent calls to the same actor - actor.par.qstar = actor.S.qstar - actor.par.alpha = actor.S.alpha - - return actor -end - -function MXHEquilibrium_to_dd!(eq::IMAS.equilibrium, mxh_eq::MXHEquilibrium.AbstractEquilibrium, ngrid::Int; cocos_in::Int=11) - tc = MXHEquilibrium.transform_cocos(cocos_in, 11) - - rlims = MXHEquilibrium.limits(mxh_eq.S; pad=0.3)[1] - zlims = MXHEquilibrium.limits(mxh_eq.S; pad=0.3)[2] - - eqt = eq.time_slice[] - Ip = eqt.global_quantities.ip - sign_Ip = sign(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 = mxh_eq.S.R0 - @ddtime eq.vacuum_toroidal_field.b0 = mxh_eq.B0 * sign_Bt - - t0 = eqt.time - empty!(eqt) - eqt.time = t0 - - eqt.global_quantities.ip = Ip - eqt.boundary.geometric_axis.r = mxh_eq.S.R0 - eqt.boundary.geometric_axis.z = Z0 - orig_psi = collect(range(MXHEquilibrium.psi_limits(mxh_eq)..., ngrid)) - eqt.profiles_1d.psi = orig_psi * (tc["PSI"] * sign_Ip) - - eqt.profiles_1d.pressure = MXHEquilibrium.pressure.(mxh_eq, orig_psi) - eqt.profiles_1d.dpressure_dpsi = MXHEquilibrium.pressure_gradient.(mxh_eq, orig_psi) ./ (tc["PSI"] * sign_Ip) - - eqt.profiles_1d.f = MXHEquilibrium.poloidal_current.(mxh_eq, orig_psi) .* (tc["F"] * sign_Bt) - eqt.profiles_1d.f_df_dpsi = - MXHEquilibrium.poloidal_current.(mxh_eq, orig_psi) .* MXHEquilibrium.poloidal_current_gradient.(mxh_eq, orig_psi) .* (tc["F_FPRIME"] * sign_Bt * sign_Ip) - - eqt2d = resize!(eqt.profiles_2d, 1)[1] - eqt2d.grid_type.index = 1 - eqt2d.grid.dim1 = range(rlims[1], rlims[2], ngrid) - eqt2d.grid.dim2 = range(zlims[1] + Z0, zlims[2] + Z0, Int(ceil(ngrid * mxh_eq.S.κ))) - eqt2d.psi = [mxh_eq(rr, flip_z * (zz - Z0)) * (tc["PSI"] * sign_Ip) for rr in eqt2d.grid.dim1, zz in eqt2d.grid.dim2] - - fw = IMAS.first_wall(IMAS.top_dd(eqt).wall) - IMAS.flux_surfaces(eqt, fw.r, fw.z) - - return -end \ No newline at end of file From dbd3e13ff41dfb7139537daa90b0eb644a1b3d90 Mon Sep 17 00:00:00 2001 From: Orso Meneghini Date: Fri, 13 Sep 2024 15:06:44 -0700 Subject: [PATCH 02/19] switch_get_from with default --- src/actors.jl | 18 +++++++++--------- src/actors/equilibrium/equilibrium_actor.jl | 2 +- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/src/actors.jl b/src/actors.jl index 668da39cc..97e3d5a71 100644 --- a/src/actors.jl +++ b/src/actors.jl @@ -50,24 +50,24 @@ end # switch_get_from # #= =============== =# """ - switch_get_from(quantity::Symbol) + switch_get_from(quantity::Symbol; default::Union{Symbol,Missing}=missing)::Switch{Symbol} -Switch to pick form which IDS `quantity` comes from +Switch to pick the IDS that `quantity` comes from """ -function switch_get_from(quantity::Symbol)::Switch{Symbol} +function switch_get_from(quantity::Symbol; default::Union{Symbol,Missing}=missing)::Switch{Symbol} txt = "Take $quantity from this IDS" if quantity == :ip - swch = Switch{Symbol}([:core_profiles, :equilibrium, :pulse_schedule], "-", txt) + swch = Switch{Symbol}([:core_profiles, :equilibrium, :pulse_schedule], "-", txt; default) elseif quantity == :vacuum_r0_b0 - swch = Switch{Symbol}([:equilibrium, :pulse_schedule], "-", txt; default=:pulse_schedule) + swch = Switch{Symbol}([:equilibrium, :pulse_schedule], "-", txt; default) elseif quantity == :vloop - swch = Switch{Symbol}([:core_profiles, :equilibrium, :pulse_schedule, :controllers__ip], "-", txt) + swch = Switch{Symbol}([:core_profiles, :equilibrium, :pulse_schedule, :controllers__ip], "-", txt; default) elseif quantity == :βn - swch = Switch{Symbol}([:core_profiles, :equilibrium], "-", txt) + swch = Switch{Symbol}([:core_profiles, :equilibrium], "-", txt; default) elseif quantity == :ne_ped - swch = Switch{Symbol}([:core_profiles, :summary, :pulse_schedule], "-", txt) + swch = Switch{Symbol}([:core_profiles, :summary, :pulse_schedule], "-", txt; default) elseif quantity == :zeff_ped - swch = Switch{Symbol}([:core_profiles, :summary, :pulse_schedule], "-", txt) + swch = Switch{Symbol}([:core_profiles, :summary, :pulse_schedule], "-", txt; default) else error("`$quantity` not supported in switch_get_from()") end diff --git a/src/actors/equilibrium/equilibrium_actor.jl b/src/actors/equilibrium/equilibrium_actor.jl index 280c2a203..ccd7c2c64 100644 --- a/src/actors/equilibrium/equilibrium_actor.jl +++ b/src/actors/equilibrium/equilibrium_actor.jl @@ -11,7 +11,7 @@ Base.@kwdef mutable struct FUSEparameters__ActorEquilibrium{T<:Real} <: Paramete #== data flow parameters ==# j_p_from::Switch{Symbol} = Switch{Symbol}([:equilibrium, :core_profiles], "-", "Take j_tor and pressure profiles from this IDS"; default=:core_profiles) ip_from::Switch{Symbol} = switch_get_from(:ip) - vacuum_r0_b0_from::Switch{Symbol} = switch_get_from(:vacuum_r0_b0) + vacuum_r0_b0_from::Switch{Symbol} = switch_get_from(:vacuum_r0_b0; default=:pulse_schedule) #== display and debugging parameters ==# do_plot::Entry{Bool} = act_common_parameters(; do_plot=false) end From ee6c65e140519c04a8956fa9fb310da800c41f5c Mon Sep 17 00:00:00 2001 From: Orso Meneghini Date: Fri, 13 Sep 2024 15:08:16 -0700 Subject: [PATCH 03/19] fresco p',ff' from p,j and previous eqt --- src/actors/equilibrium/chease_actor.jl | 34 +++++++++--------- src/actors/equilibrium/equilibrium_actor.jl | 37 +++++++++++++++---- src/actors/equilibrium/fresco_actor.jl | 39 ++++++++++---------- src/actors/equilibrium/tequila_actor.jl | 40 ++++++++++----------- 4 files changed, 86 insertions(+), 64 deletions(-) diff --git a/src/actors/equilibrium/chease_actor.jl b/src/actors/equilibrium/chease_actor.jl index acb896130..99f3ca395 100644 --- a/src/actors/equilibrium/chease_actor.jl +++ b/src/actors/equilibrium/chease_actor.jl @@ -50,7 +50,7 @@ function _step(actor::ActorCHEASE) # initialize eqt from pulse_schedule and core_profiles eqt = dd.equilibrium.time_slice[] - eq1d = eqt.profiles_1d + eqt1d = eqt.profiles_1d # boundary pr = eqt.boundary.outline.r @@ -69,9 +69,9 @@ function _step(actor::ActorCHEASE) ϵ = eqt.boundary.minor_radius / r_geo # pressure and j_tor - psin = eq1d.psi_norm - j_tor = [sign(j) == sign(Ip) ? j : 0.0 for j in eq1d.j_tor] - pressure = eq1d.pressure + psin = eqt1d.psi_norm + j_tor = [sign(j) == sign(Ip) ? j : 0.0 for j in eqt1d.j_tor] + pressure = eqt1d.pressure rho_pol = sqrt.(psin) pressure_sep = pressure[end] @@ -164,8 +164,8 @@ function gEQDSK2IMAS(g::CHEASE.EFIT.GEQDSKFile, eq::IMAS.equilibrium) tc = MXHEquilibrium.transform_cocos(1, 11) # chease output is cocos 1 , dd is cocos 11 eqt = eq.time_slice[] - eq1d = eqt.profiles_1d - eq2d = resize!(eqt.profiles_2d, 1)[1] + eqt1d = eqt.profiles_1d + eqt2d = resize!(eqt.profiles_2d, 1)[1] @ddtime(eq.vacuum_toroidal_field.b0 = g.bcentr) eq.vacuum_toroidal_field.r0 = g.rcentr @@ -176,17 +176,17 @@ function gEQDSK2IMAS(g::CHEASE.EFIT.GEQDSKFile, eq::IMAS.equilibrium) eqt.global_quantities.magnetic_axis.z = g.zmaxis eqt.global_quantities.ip = g.current - eq1d.psi = g.psi .* tc["PSI"] - eq1d.q = g.qpsi - eq1d.pressure = g.pres - eq1d.dpressure_dpsi = g.pprime .* tc["PPRIME"] - eq1d.f = g.fpol .* tc["F"] - eq1d.f_df_dpsi = g.ffprim .* tc["F_FPRIME"] - - eq2d.grid_type.index = 1 - eq2d.grid.dim1 = g.r - eq2d.grid.dim2 = g.z - eq2d.psi = g.psirz .* tc["PSI"] + eqt1d.psi = g.psi .* tc["PSI"] + eqt1d.q = g.qpsi + eqt1d.pressure = g.pres + eqt1d.dpressure_dpsi = g.pprime .* tc["PPRIME"] + eqt1d.f = g.fpol .* tc["F"] + eqt1d.f_df_dpsi = g.ffprim .* tc["F_FPRIME"] + + eqt2d.grid_type.index = 1 + eqt2d.grid.dim1 = g.r + eqt2d.grid.dim2 = g.z + eqt2d.psi = g.psirz .* tc["PSI"] return nothing end diff --git a/src/actors/equilibrium/equilibrium_actor.jl b/src/actors/equilibrium/equilibrium_actor.jl index ccd7c2c64..7cc0e3c83 100644 --- a/src/actors/equilibrium/equilibrium_actor.jl +++ b/src/actors/equilibrium/equilibrium_actor.jl @@ -70,6 +70,31 @@ function _step(actor::ActorEquilibrium) # initialize eqt for equilibrium actors prepare(actor) + # FRESCO needs p' and f'' to run + # We can estimate of p' and ff' based on equilibrium at previous time-slice + # if that's not available, we run TEQUILA first + if par.model == :FRESCO + eqt = dd.equilibrium.time_slice[] + idxeqt = IMAS.index(eqt) + if idxeqt != 1 && !ismissing(dd.equilibrium.time_slice[idxeqt-1].global_quantities, :ip) + eqt1d__1 = dd.equilibrium.time_slice[idxeqt-1].profiles_1d + psi__1 = IMAS.interp1d(eqt1d__1.psi_norm, eqt1d__1.psi).(eqt1d.psi_norm) + R__1 = IMAS.interp1d(eqt1d__1.psi_norm, eqt1d__1.gm8).(eqt1d.psi_norm) + one_R__1 = IMAS.interp1d(eqt1d__1.psi_norm, eqt1d__1.gm9).(eqt1d.psi_norm) + one_R2__1 = IMAS.interp1d(eqt1d__1.psi_norm, eqt1d__1.gm1).(eqt1d.psi_norm) + b0 = eqt.global_quantities.vacuum_toroidal_field.b0 + r0 = eqt.global_quantities.vacuum_toroidal_field.r0 + dpressure_dpsi, f_df_dpsi, f = IMAS.calc_pprime_ffprim_f(psi__1, R__1, one_R__1, one_R2__1, r0, b0; eqt1d.pressure, eqt1d.j_tor) + eqt1d.dpressure_dpsi = dpressure_dpsi + eqt1d.f_df_dpsi = f_df_dpsi + eqt1d.f = f + else + ActorTEQUILA(dd, actor.act; par.ip_from) + fw = IMAS.first_wall(dd.wall) + IMAS.flux_surfaces(eqt, fw.r, fw.z) + end + end + # step selected equilibrium actor step(actor.eq_actor) @@ -173,7 +198,7 @@ function prepare(actor::ActorEquilibrium) # add/clear time-slice eqt = resize!(dd.equilibrium.time_slice) resize!(eqt.profiles_2d, 1) - eq1d = dd.equilibrium.time_slice[].profiles_1d + eqt1d = dd.equilibrium.time_slice[].profiles_1d # scalar quantities eqt.global_quantities.ip = ip @@ -225,11 +250,11 @@ function prepare(actor::ActorEquilibrium) end # set j_tor and pressure, forcing zero derivative on axis - eq1d = dd.equilibrium.time_slice[].profiles_1d - eq1d.psi = psi0 - eq1d.rho_tor_norm = rho_tor_norm0 - eq1d.j_tor = j_itp.(sqrt.(eq1d.psi_norm)) - eq1d.pressure = p_itp.(sqrt.(eq1d.psi_norm)) + eqt1d = dd.equilibrium.time_slice[].profiles_1d + eqt1d.psi = psi0 + eqt1d.rho_tor_norm = rho_tor_norm0 + eqt1d.j_tor = j_itp.(sqrt.(eqt1d.psi_norm)) + eqt1d.pressure = p_itp.(sqrt.(eqt1d.psi_norm)) return dd end diff --git a/src/actors/equilibrium/fresco_actor.jl b/src/actors/equilibrium/fresco_actor.jl index 3d902bb48..f12ecf211 100644 --- a/src/actors/equilibrium/fresco_actor.jl +++ b/src/actors/equilibrium/fresco_actor.jl @@ -18,8 +18,8 @@ Base.@kwdef mutable struct FUSEparameters__ActorFRESCO{T<:Real} <: ParametersAct do_plot::Entry{Bool} = act_common_parameters(; do_plot=false) debug::Entry{Bool} = Entry{Bool}("-", "Print debug information withing FRESCO solve"; default=false) #== IMAS psi grid settings ==# - nR::Entry{Int} = Entry{Int}("-", "Grid resolution along R") - nZ::Entry{Int} = Entry{Int}("-", "Grid resolution along Z") + nR::Entry{Int} = Entry{Int}("-", "Grid resolution along R"; default=129) + nZ::Entry{Int} = Entry{Int}("-", "Grid resolution along Z"; default=129) end mutable struct ActorFRESCO{D,P} <: SingleAbstractActor{D,P} @@ -55,35 +55,32 @@ function _step(actor::ActorFRESCO) dd = actor.dd par = actor.par eqt = dd.equilibrium.time_slice[] - eq1d = eqt.profiles_1d + eqt1d = eqt.profiles_1d # Ip_target = eqt.global_quantities.ip # # if par.fixed_grid === :poloidal - # rho = sqrt.(eq1d.psi_norm) + # rho = sqrt.(eqt1d.psi_norm) # rho[1] = 0.0 - # P = (FRESCO.FE(rho, eq1d.pressure), :poloidal) + # P = (FRESCO.FE(rho, eqt1d.pressure), :poloidal) # # don't allow current to change sign - # Jt = (FRESCO.FE(rho, [sign(j) == sign(Ip_target) ? j : 0.0 for j in eq1d.j_tor]), :poloidal) - # Pbnd = eq1d.pressure[end] + # Jt = (FRESCO.FE(rho, [sign(j) == sign(Ip_target) ? j : 0.0 for j in eqt1d.j_tor]), :poloidal) + # Pbnd = eqt1d.pressure[end] # elseif par.fixed_grid === :toroidal - # rho = eq1d.rho_tor_norm - # P = (FRESCO.FE(rho, eq1d.pressure), :toroidal) + # rho = eqt1d.rho_tor_norm + # P = (FRESCO.FE(rho, eqt1d.pressure), :toroidal) # # don't allow current to change sign - # Jt = (FRESCO.FE(rho, [sign(j) == sign(Ip_target) ? j : 0.0 for j in eq1d.j_tor]), :toroidal) - # Pbnd = eq1d.pressure[end] + # Jt = (FRESCO.FE(rho, [sign(j) == sign(Ip_target) ? j : 0.0 for j in eqt1d.j_tor]), :toroidal) + # Pbnd = eqt1d.pressure[end] # end - psi_norm = eq1d.psi_norm - gpp = IMAS.interp1d(psi_norm, eq1d.dpressure_dpsi, :cubic) - pprime = x -> gpp(x) - gffp = IMAS.interp1d(psi_norm, eq1d.f_df_dpsi, :cubic) - ffprime = x -> gffp(x) - profile = FRESCO.PprimeFFprime(pprime, ffprime) + gpp = IMAS.interp1d(eqt1d.psi_norm, eqt1d.dpressure_dpsi, :cubic) + gffp = IMAS.interp1d(eqt1d.psi_norm, eqt1d.f_df_dpsi, :cubic) + profile = FRESCO.PprimeFFprime(x -> gpp(x), x -> gffp(x)) - actor.canvas = FRESCO.Canvas(dd, par.nR, par.nZ); + actor.canvas = FRESCO.Canvas(dd, par.nR, par.nZ) - FRESCO.solve!(actor.canvas, profile, 100, 3; relax=par.relax, par.debug, par.control) # forward mode + FRESCO.solve!(actor.canvas, profile, 100, 3; relax=par.relax, par.debug, par.control) # using Plots # p1 = Plots.heatmap(Rs, Zs, psi', aspect_ratio=:equal, xrange=(0,12.5), yrange=(-9,9), size=(400,500)) @@ -100,13 +97,13 @@ function _finalize(actor::ActorFRESCO) dd = actor.dd eq = dd.equilibrium eqt = eq.time_slice[] - eq1d = eqt.profiles_1d + eqt1d = eqt.profiles_1d eq2d = resize!(eqt.profiles_2d, 1)[1] Raxis, Zaxis, _ = FRESCO.find_axis(canvas) eqt.global_quantities.magnetic_axis.r = Raxis eqt.global_quantities.magnetic_axis.z = Zaxis - eq1d.psi = range(canvas.Ψaxis, canvas.Ψbnd, length(eq1d.psi)) + eqt1d.psi = range(canvas.Ψaxis, canvas.Ψbnd, length(eqt1d.psi)) # p, p', f, ff' don't change eq2d.grid_type.index = 1 diff --git a/src/actors/equilibrium/tequila_actor.jl b/src/actors/equilibrium/tequila_actor.jl index 6713de039..90c1175c6 100644 --- a/src/actors/equilibrium/tequila_actor.jl +++ b/src/actors/equilibrium/tequila_actor.jl @@ -62,27 +62,27 @@ function _step(actor::ActorTEQUILA) dd = actor.dd par = actor.par eqt = dd.equilibrium.time_slice[] - eq1d = eqt.profiles_1d + eqt1d = eqt.profiles_1d # BCL 5/30/23: ψbound should be set time dependently, related to the flux swing of the OH coils - # For now setting to zero as initial eq1d.psi profile from prepare() can be nonsense + # For now setting to zero as initial eqt1d.psi profile from prepare() can be nonsense actor.ψbound = 0.0 Ip_target = eqt.global_quantities.ip if par.fixed_grid === :poloidal - rho = sqrt.(eq1d.psi_norm) + rho = sqrt.(eqt1d.psi_norm) rho[1] = 0.0 - P = (TEQUILA.FE(rho, eq1d.pressure), :poloidal) + P = (TEQUILA.FE(rho, eqt1d.pressure), :poloidal) # don't allow current to change sign - Jt = (TEQUILA.FE(rho, [sign(j) == sign(Ip_target) ? j : 0.0 for j in eq1d.j_tor]), :poloidal) - Pbnd = eq1d.pressure[end] + Jt = (TEQUILA.FE(rho, [sign(j) == sign(Ip_target) ? j : 0.0 for j in eqt1d.j_tor]), :poloidal) + Pbnd = eqt1d.pressure[end] elseif par.fixed_grid === :toroidal - rho = eq1d.rho_tor_norm - P = (TEQUILA.FE(rho, eq1d.pressure), :toroidal) + rho = eqt1d.rho_tor_norm + P = (TEQUILA.FE(rho, eqt1d.pressure), :toroidal) # don't allow current to change sign - Jt = (TEQUILA.FE(rho, [sign(j) == sign(Ip_target) ? j : 0.0 for j in eq1d.j_tor]), :toroidal) - Pbnd = eq1d.pressure[end] + Jt = (TEQUILA.FE(rho, [sign(j) == sign(Ip_target) ? j : 0.0 for j in eqt1d.j_tor]), :toroidal) + Pbnd = eqt1d.pressure[end] end Fbnd = eqt.global_quantities.vacuum_toroidal_field.b0 * eqt.global_quantities.vacuum_toroidal_field.r0 @@ -124,8 +124,8 @@ function _step(actor::ActorTEQUILA) catch e plot(eqt.boundary.outline.r, eqt.boundary.outline.z; marker=:dot, aspect_ratio=:equal) display(plot!(IMAS.MXH(actor.shot.surfaces[:, end]))) - display(plot(rho, eq1d.pressure; marker=:dot, xlabel="ρ", title="Pressure [Pa]")) - display(plot(rho, eq1d.j_tor; marker=:dot, xlabel="ρ", title="Jtor [A]")) + display(plot(rho, eqt1d.pressure; marker=:dot, xlabel="ρ", title="Pressure [Pa]")) + display(plot(rho, eqt1d.j_tor; marker=:dot, xlabel="ρ", title="Jtor [A]")) rethrow(e) end @@ -148,7 +148,7 @@ function tequila2imas(shot::TEQUILA.Shot, dd::IMAS.dd, par::FUSEparameters__Acto free_boundary = par.free_boundary eq = dd.equilibrium eqt = eq.time_slice[] - eq1d = eqt.profiles_1d + eqt1d = eqt.profiles_1d R0 = shot.R0fe(1.0) Z0 = shot.Z0fe(1.0) @@ -173,12 +173,12 @@ function tequila2imas(shot::TEQUILA.Shot, dd::IMAS.dd, par::FUSEparameters__Acto psit = shot.C[2:2:end, 1] psia = psit[1] psib = psit[end] - eq1d.psi = range(psia, psib, nψ_grid) - rhoi = TEQUILA.ρ.(Ref(shot), eq1d.psi) - eq1d.pressure = MXHEquilibrium.pressure.(Ref(shot), eq1d.psi) - eq1d.dpressure_dpsi = MXHEquilibrium.pressure_gradient.(Ref(shot), eq1d.psi) - eq1d.f = shot.F.(rhoi) - eq1d.f_df_dpsi = TEQUILA.Fpol_dFpol_dψ.(Ref(shot), rhoi; shot.invR, shot.invR2) + eqt1d.psi = range(psia, psib, nψ_grid) + rhoi = TEQUILA.ρ.(Ref(shot), eqt1d.psi) + eqt1d.pressure = MXHEquilibrium.pressure.(Ref(shot), eqt1d.psi) + eqt1d.dpressure_dpsi = MXHEquilibrium.pressure_gradient.(Ref(shot), eqt1d.psi) + eqt1d.f = shot.F.(rhoi) + eqt1d.f_df_dpsi = TEQUILA.Fpol_dFpol_dψ.(Ref(shot), rhoi; shot.invR, shot.invR2) resize!(eqt.profiles_2d, 2) @@ -254,7 +254,7 @@ function tequila2imas(shot::TEQUILA.Shot, dd::IMAS.dd, par::FUSEparameters__Acto eq2d.psi[i, j] = shot(r, z; extrapolate=true) + ψbound end end - eq1d.psi .+= ψbound + eqt1d.psi .+= ψbound end end From 0be21dd1b2001fa132d84f0ff77d2f519e70dd86 Mon Sep 17 00:00:00 2001 From: Orso Meneghini Date: Fri, 13 Sep 2024 15:22:28 -0700 Subject: [PATCH 04/19] bug fix --- src/actors/equilibrium/equilibrium_actor.jl | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/src/actors/equilibrium/equilibrium_actor.jl b/src/actors/equilibrium/equilibrium_actor.jl index 7cc0e3c83..719535fda 100644 --- a/src/actors/equilibrium/equilibrium_actor.jl +++ b/src/actors/equilibrium/equilibrium_actor.jl @@ -77,6 +77,7 @@ function _step(actor::ActorEquilibrium) eqt = dd.equilibrium.time_slice[] idxeqt = IMAS.index(eqt) if idxeqt != 1 && !ismissing(dd.equilibrium.time_slice[idxeqt-1].global_quantities, :ip) + eqt1d = eqt.profiles_1d eqt1d__1 = dd.equilibrium.time_slice[idxeqt-1].profiles_1d psi__1 = IMAS.interp1d(eqt1d__1.psi_norm, eqt1d__1.psi).(eqt1d.psi_norm) R__1 = IMAS.interp1d(eqt1d__1.psi_norm, eqt1d__1.gm8).(eqt1d.psi_norm) @@ -156,7 +157,7 @@ Prepare `dd.equilibrium` to run equilibrium actors - Clear equilibrium__time_slice - Set Ip, Bt from core_profiles, equilibrium, or pulse_schedule - Use position control from pulse_schedule - - Use j_tor,pressure from core_profiles or equilibrium + - Use j_tor,pressure from core_profiles (for self-consistent iterations) or equilibrium (to re-solve equilibrium with different solver) """ function prepare(actor::ActorEquilibrium) dd = actor.dd @@ -172,21 +173,21 @@ function prepare(actor::ActorEquilibrium) index = cp1d.grid.psi_norm .> 0.05 psi0 = cp1d.grid.psi rho_tor_norm0 = cp1d.grid.rho_tor_norm - rho_pol_norm0 = vcat(-reverse(sqrt.(cp1d.grid.psi_norm[index])), sqrt.(cp1d.grid.psi_norm[index])) + rho_pol_norm_sqrt0 = vcat(-reverse(sqrt.(cp1d.grid.psi_norm[index])), sqrt.(cp1d.grid.psi_norm[index])) j_tor0 = vcat(reverse(cp1d.j_tor[index]), cp1d.j_tor[index]) pressure0 = vcat(reverse(cp1d.pressure[index]), cp1d.pressure[index]) - j_itp = IMAS.interp1d(rho_pol_norm0, j_tor0, :cubic) - p_itp = IMAS.interp1d(rho_pol_norm0, pressure0, :cubic) + j_itp = IMAS.interp1d(rho_pol_norm_sqrt0, j_tor0, :cubic) + p_itp = IMAS.interp1d(rho_pol_norm_sqrt0, pressure0, :cubic) elseif par.j_p_from == :equilibrium @assert !isempty(dd.equilibrium.time) eqt1d = dd.equilibrium.time_slice[].profiles_1d psi0 = eqt1d.psi rho_tor_norm0 = eqt1d.rho_tor_norm - rho_pol_norm0 = sqrt.(eqt1d.psi_norm) + rho_pol_norm_sqrt0 = sqrt.(eqt1d.psi_norm) j_tor0 = eqt1d.j_tor pressure0 = eqt1d.pressure - j_itp = IMAS.interp1d(rho_pol_norm0, j_tor0, :cubic) - p_itp = IMAS.interp1d(rho_pol_norm0, pressure0, :cubic) + j_itp = IMAS.interp1d(rho_pol_norm_sqrt0, j_tor0, :cubic) + p_itp = IMAS.interp1d(rho_pol_norm_sqrt0, pressure0, :cubic) else @assert par.j_p_from in (:core_profiles, :equilibrium) end From 816c91378df3ad49a9a3834616bd43173db6cbad Mon Sep 17 00:00:00 2001 From: Orso Meneghini Date: Sat, 14 Sep 2024 21:21:31 -0700 Subject: [PATCH 05/19] remove unecessary ip_from --- src/actors/equilibrium/chease_actor.jl | 2 -- src/actors/equilibrium/equilibrium_actor.jl | 6 +++--- src/actors/equilibrium/fresco_actor.jl | 1 - src/actors/equilibrium/tequila_actor.jl | 2 -- 4 files changed, 3 insertions(+), 8 deletions(-) diff --git a/src/actors/equilibrium/chease_actor.jl b/src/actors/equilibrium/chease_actor.jl index 99f3ca395..0ded00d7d 100644 --- a/src/actors/equilibrium/chease_actor.jl +++ b/src/actors/equilibrium/chease_actor.jl @@ -11,8 +11,6 @@ Base.@kwdef mutable struct FUSEparameters__ActorCHEASE{T<:Real} <: ParametersAct free_boundary::Entry{Bool} = Entry{Bool}("-", "Convert fixed boundary equilibrium to free boundary one"; default=true) clear_workdir::Entry{Bool} = Entry{Bool}("-", "Clean the temporary workdir for CHEASE"; default=true) rescale_eq_to_ip::Entry{Bool} = Entry{Bool}("-", "Scale equilibrium to match Ip"; default=true) - #== data flow parameters ==# - ip_from::Switch{Symbol} = switch_get_from(:ip) end mutable struct ActorCHEASE{D,P} <: SingleAbstractActor{D,P} diff --git a/src/actors/equilibrium/equilibrium_actor.jl b/src/actors/equilibrium/equilibrium_actor.jl index 719535fda..e9a952201 100644 --- a/src/actors/equilibrium/equilibrium_actor.jl +++ b/src/actors/equilibrium/equilibrium_actor.jl @@ -39,11 +39,11 @@ function ActorEquilibrium(dd::IMAS.dd, par::FUSEparameters__ActorEquilibrium, ac logging_actor_init(ActorEquilibrium) par = par(kw...) if par.model == :FRESCO - eq_actor = ActorFRESCO(dd, act.ActorFRESCO; par.ip_from) + eq_actor = ActorFRESCO(dd, act.ActorFRESCO) elseif par.model == :CHEASE - eq_actor = ActorCHEASE(dd, act.ActorCHEASE; par.ip_from) + eq_actor = ActorCHEASE(dd, act.ActorCHEASE) elseif par.model == :TEQUILA - eq_actor = ActorTEQUILA(dd, act.ActorTEQUILA; par.ip_from) + eq_actor = ActorTEQUILA(dd, act.ActorTEQUILA) else error("ActorEquilibrium: model = `$(par.model)` can only be `:TEQUILA`, `:FRESCO`, `:CHEASE`") end diff --git a/src/actors/equilibrium/fresco_actor.jl b/src/actors/equilibrium/fresco_actor.jl index f12ecf211..ef2cd4f5f 100644 --- a/src/actors/equilibrium/fresco_actor.jl +++ b/src/actors/equilibrium/fresco_actor.jl @@ -12,7 +12,6 @@ Base.@kwdef mutable struct FUSEparameters__ActorFRESCO{T<:Real} <: ParametersAct relax::Entry{Float64} = Entry{Float64}("-", "Relaxation on the Picard iterations"; default=0.5) #tolerance::Entry{Float64} = Entry{Float64}("-", "Tolerance for terminating iterations"; default=1e-4) #== data flow parameters ==# - ip_from::Switch{Symbol} = switch_get_from(:ip) # fixed_grid::Switch{Symbol} = Switch{Symbol}([:poloidal, :toroidal], "-", "Fix P and Jt on this rho grid"; default=:toroidal) #== display and debugging parameters ==# do_plot::Entry{Bool} = act_common_parameters(; do_plot=false) diff --git a/src/actors/equilibrium/tequila_actor.jl b/src/actors/equilibrium/tequila_actor.jl index 90c1175c6..0af1d0094 100644 --- a/src/actors/equilibrium/tequila_actor.jl +++ b/src/actors/equilibrium/tequila_actor.jl @@ -16,7 +16,6 @@ Base.@kwdef mutable struct FUSEparameters__ActorTEQUILA{T<:Real} <: ParametersAc relax::Entry{Float64} = Entry{Float64}("-", "Relaxation on the Picard iterations"; default=0.25) tolerance::Entry{Float64} = Entry{Float64}("-", "Tolerance for terminating iterations"; default=1e-4) #== data flow parameters ==# - ip_from::Switch{Symbol} = switch_get_from(:ip) fixed_grid::Switch{Symbol} = Switch{Symbol}([:poloidal, :toroidal], "-", "Fix P and Jt on this rho grid"; default=:toroidal) #== display and debugging parameters ==# do_plot::Entry{Bool} = act_common_parameters(; do_plot=false) @@ -69,7 +68,6 @@ function _step(actor::ActorTEQUILA) actor.ψbound = 0.0 Ip_target = eqt.global_quantities.ip - if par.fixed_grid === :poloidal rho = sqrt.(eqt1d.psi_norm) rho[1] = 0.0 From 98afcb89aa70788369a6a0c9eade18f60045a6ab Mon Sep 17 00:00:00 2001 From: Orso Meneghini Date: Sat, 14 Sep 2024 21:22:19 -0700 Subject: [PATCH 06/19] fix show of pf_active --- src/actors/pf/pf_active_utils.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/actors/pf/pf_active_utils.jl b/src/actors/pf/pf_active_utils.jl index edcad4c9e..29124f157 100644 --- a/src/actors/pf/pf_active_utils.jl +++ b/src/actors/pf/pf_active_utils.jl @@ -329,7 +329,7 @@ function DataFrames.DataFrame(coils::IMAS.IDSvector{<:IMAS.pf_active__coil}) return df end -function Base.show(io::IO, mime, ::MIME"text/plain", coils::IMAS.IDSvector{<:Union{IMAS.pf_active__coil,IMAS.pf_active__supply}}) +function Base.show(io::IO, mime::MIME"text/plain", coils::IMAS.IDSvector{<:Union{IMAS.pf_active__coil,IMAS.pf_active__supply}}) old_lines = get(ENV, "LINES", missing) old_columns = get(ENV, "COLUMNS", missing) df = DataFrames.DataFrame(coils) From 5056f170b94bef15c86b73276f828949af147440 Mon Sep 17 00:00:00 2001 From: Orso Meneghini Date: Sat, 14 Sep 2024 21:22:56 -0700 Subject: [PATCH 07/19] expose more FRESCO parameters --- src/actors/equilibrium/equilibrium_actor.jl | 8 +++++--- src/actors/equilibrium/fresco_actor.jl | 12 +++++++----- 2 files changed, 12 insertions(+), 8 deletions(-) diff --git a/src/actors/equilibrium/equilibrium_actor.jl b/src/actors/equilibrium/equilibrium_actor.jl index e9a952201..6e72ea584 100644 --- a/src/actors/equilibrium/equilibrium_actor.jl +++ b/src/actors/equilibrium/equilibrium_actor.jl @@ -6,7 +6,7 @@ Base.@kwdef mutable struct FUSEparameters__ActorEquilibrium{T<:Real} <: Paramete _name::Symbol = :not_set _time::Float64 = NaN #== actor parameters ==# - model::Switch{Symbol} = Switch{Symbol}([:FRESCO, :CHEASE, :TEQUILA], "-", "Equilibrium actor to run"; default=:TEQUILA) + model::Switch{Symbol} = Switch{Symbol}([:TEQUILA, :FRESCO, :CHEASE], "-", "Equilibrium actor to run"; default=:TEQUILA) symmetrize::Entry{Bool} = Entry{Bool}("-", "Force equilibrium up-down symmetry with respect to magnetic axis"; default=false) #== data flow parameters ==# j_p_from::Switch{Symbol} = Switch{Symbol}([:equilibrium, :core_profiles], "-", "Take j_tor and pressure profiles from this IDS"; default=:core_profiles) @@ -76,7 +76,7 @@ function _step(actor::ActorEquilibrium) if par.model == :FRESCO eqt = dd.equilibrium.time_slice[] idxeqt = IMAS.index(eqt) - if idxeqt != 1 && !ismissing(dd.equilibrium.time_slice[idxeqt-1].global_quantities, :ip) + if false && idxeqt != 1 && !ismissing(dd.equilibrium.time_slice[idxeqt-1].global_quantities, :ip) eqt1d = eqt.profiles_1d eqt1d__1 = dd.equilibrium.time_slice[idxeqt-1].profiles_1d psi__1 = IMAS.interp1d(eqt1d__1.psi_norm, eqt1d__1.psi).(eqt1d.psi_norm) @@ -90,7 +90,9 @@ function _step(actor::ActorEquilibrium) eqt1d.f_df_dpsi = f_df_dpsi eqt1d.f = f else - ActorTEQUILA(dd, actor.act; par.ip_from) + # NOTE: using free_boundary=true will also update pf_active + # which is then used as initial guess by FRESCO + ActorTEQUILA(dd, actor.act; free_boundary=true) fw = IMAS.first_wall(dd.wall) IMAS.flux_surfaces(eqt, fw.r, fw.z) end diff --git a/src/actors/equilibrium/fresco_actor.jl b/src/actors/equilibrium/fresco_actor.jl index ef2cd4f5f..c761c2213 100644 --- a/src/actors/equilibrium/fresco_actor.jl +++ b/src/actors/equilibrium/fresco_actor.jl @@ -9,13 +9,14 @@ Base.@kwdef mutable struct FUSEparameters__ActorFRESCO{T<:Real} <: ParametersAct _time::Float64 = NaN #== actor parameters ==# control::Switch{Symbol} = Switch{Symbol}([:vertical, :shape], "-", ""; default=:shape) + number_of_iterations::Entry{Tuple{Int,Int}} = Entry{Tuple{Int,Int}}("-", "Number of outer and inner iterations"; default=(100, 3)) relax::Entry{Float64} = Entry{Float64}("-", "Relaxation on the Picard iterations"; default=0.5) - #tolerance::Entry{Float64} = Entry{Float64}("-", "Tolerance for terminating iterations"; default=1e-4) + tolerance::Entry{Float64} = Entry{Float64}("-", "Tolerance for terminating iterations"; default=1e-4) #== data flow parameters ==# - # fixed_grid::Switch{Symbol} = Switch{Symbol}([:poloidal, :toroidal], "-", "Fix P and Jt on this rho grid"; default=:toroidal) + fixed_grid::Switch{Symbol} = Switch{Symbol}([:poloidal, :toroidal], "-", "Fix P and Jt on this rho grid"; default=:toroidal) #== display and debugging parameters ==# do_plot::Entry{Bool} = act_common_parameters(; do_plot=false) - debug::Entry{Bool} = Entry{Bool}("-", "Print debug information withing FRESCO solve"; default=false) + debug::Entry{Bool} = Entry{Bool}("-", "Print debug information withing FRESCO solve"; default=true) #== IMAS psi grid settings ==# nR::Entry{Int} = Entry{Int}("-", "Grid resolution along R"; default=129) nZ::Entry{Int} = Entry{Int}("-", "Grid resolution along Z"; default=129) @@ -57,7 +58,6 @@ function _step(actor::ActorFRESCO) eqt1d = eqt.profiles_1d # Ip_target = eqt.global_quantities.ip - # # if par.fixed_grid === :poloidal # rho = sqrt.(eqt1d.psi_norm) # rho[1] = 0.0 @@ -79,7 +79,7 @@ function _step(actor::ActorFRESCO) actor.canvas = FRESCO.Canvas(dd, par.nR, par.nZ) - FRESCO.solve!(actor.canvas, profile, 100, 3; relax=par.relax, par.debug, par.control) + FRESCO.solve!(actor.canvas, profile, par.number_of_iterations...; par.relax, par.debug, par.control, par.tolerance) # using Plots # p1 = Plots.heatmap(Rs, Zs, psi', aspect_ratio=:equal, xrange=(0,12.5), yrange=(-9,9), size=(400,500)) @@ -102,6 +102,8 @@ function _finalize(actor::ActorFRESCO) Raxis, Zaxis, _ = FRESCO.find_axis(canvas) eqt.global_quantities.magnetic_axis.r = Raxis eqt.global_quantities.magnetic_axis.z = Zaxis + eqt.global_quantities.psi_boundary = canvas.Ψbnd + eqt.global_quantities.psi_axis = canvas.Ψaxis eqt1d.psi = range(canvas.Ψaxis, canvas.Ψbnd, length(eqt1d.psi)) # p, p', f, ff' don't change From 929db2705b723413339312113d5252d3c274fc7e Mon Sep 17 00:00:00 2001 From: Brendan Lyons Date: Mon, 23 Sep 2024 17:44:13 -0700 Subject: [PATCH 08/19] Use properly computed Ip_non_inductive in ActorQED --- src/actors/current/qed_actor.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/actors/current/qed_actor.jl b/src/actors/current/qed_actor.jl index 4dd136c7b..a76afaf3f 100644 --- a/src/actors/current/qed_actor.jl +++ b/src/actors/current/qed_actor.jl @@ -60,7 +60,6 @@ function _step(actor::ActorQED) cp1d = dd.core_profiles.profiles_1d[] # non_inductive contribution - Ip_non_inductive = trapz(cp1d.grid.area, cp1d.j_non_inductive) B0 = eqt.global_quantities.vacuum_toroidal_field.b0 JBni = QED.FE(cp1d.grid.rho_tor_norm, cp1d.j_non_inductive .* B0) @@ -93,6 +92,7 @@ function _step(actor::ActorQED) for time0 in range(t0 + δt / 2.0, t1 + δt / 2.0, No + 1)[1:end-1] if par.solve_for == :ip Ip = IMAS.get_from(dd, Val{:ip}, par.ip_from; time0) + Ip_non_inductive = dd.core_profiles.global_quantities.current_non_inductive[time0] Vedge = nothing if par.allow_floating_plasma_current && abs(Ip) < abs(Ip_non_inductive) Ip = nothing @@ -109,6 +109,7 @@ function _step(actor::ActorQED) # steady state solution if par.solve_for == :ip Ip = IMAS.get_from(dd, Val{:ip}, par.ip_from) + Ip_non_inductive = dd.core_profiles.global_quantities.current_non_inductive[] Vedge = nothing if par.allow_floating_plasma_current && abs(Ip) < abs(Ip_non_inductive) Ip = nothing From adafa3983ccde8380ebb5dc3b1dc7ed0b8c97740 Mon Sep 17 00:00:00 2001 From: Brendan Lyons Date: Wed, 25 Sep 2024 14:12:41 -0700 Subject: [PATCH 09/19] Require QED v1.4 for coil evolution --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 76f64288c..44c6e09f5 100644 --- a/Project.toml +++ b/Project.toml @@ -97,7 +97,7 @@ OrderedCollections = "1" Plots = "1" PolygonOps = "0.1" ProgressMeter = "1" -QED = "1" +QED = "1.4" QuadGK = "2" RABBIT = "1" SimulationParameters = "1" From 8ae2b5a26455dd75216050cae1df481554af633e Mon Sep 17 00:00:00 2001 From: Brendan Lyons Date: Wed, 25 Sep 2024 14:36:14 -0700 Subject: [PATCH 10/19] Add FRESCO back to Makefile --- Makefile | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/Makefile b/Makefile index e572619a7..7e83e80a4 100644 --- a/Makefile +++ b/Makefile @@ -15,7 +15,7 @@ else FUSE_LOCAL_BRANCH=$(shell echo $(GITHUB_REF) | sed 's/refs\/heads\///') endif -FUSE_PACKAGES_MAKEFILE := ADAS BoundaryPlasmaModels CHEASE CoordinateConventions EPEDNN FiniteElementHermite Fortran90Namelists FuseUtils FusionMaterials FuseExchangeProtocol IMAS IMASdd MXHEquilibrium MeshTools MillerExtendedHarmonic NEO NNeutronics QED RABBIT SimulationParameters TEQUILA TGLFNN TJLF VacuumFields XSteam ThermalSystemModels +FUSE_PACKAGES_MAKEFILE := ADAS BoundaryPlasmaModels CHEASE CoordinateConventions EPEDNN FiniteElementHermite Fortran90Namelists FRESCO FuseUtils FusionMaterials FuseExchangeProtocol IMAS IMASdd MXHEquilibrium MeshTools MillerExtendedHarmonic NEO NNeutronics QED RABBIT SimulationParameters TEQUILA TGLFNN TJLF VacuumFields XSteam ThermalSystemModels FUSE_PACKAGES_MAKEFILE := $(sort $(FUSE_PACKAGES_MAKEFILE)) FUSE_PACKAGES := $(shell echo '$(FUSE_PACKAGES_MAKEFILE)' | awk '{printf("[\"%s\"", $$1); for (i=2; i<=NF; i++) printf(", \"%s\"", $$i); print "]"}') DEV_PACKAGES_MAKEFILE := $(shell find ../*/.git/config -exec grep ProjectTorreyPines \{\} /dev/null \; | cut -d'/' -f 2) @@ -115,6 +115,9 @@ FiniteElementHermite: Fortran90Namelists: $(call clone_pull_repo,$@) +FRESCO: + $(call clone_pull_repo,$@) + CHEASE: $(call clone_pull_repo,$@) From 86c38bbeb2ed9d7aeb9115c4bc347444add89a2a Mon Sep 17 00:00:00 2001 From: Brendan Lyons Date: Wed, 25 Sep 2024 22:28:18 -0700 Subject: [PATCH 11/19] WIP: Create ActorPlasmaBuild --- src/actors/current/qed-build_actor.jl | 207 ++++++++++++++++++++++++++ 1 file changed, 207 insertions(+) create mode 100644 src/actors/current/qed-build_actor.jl diff --git a/src/actors/current/qed-build_actor.jl b/src/actors/current/qed-build_actor.jl new file mode 100644 index 000000000..cdb349305 --- /dev/null +++ b/src/actors/current/qed-build_actor.jl @@ -0,0 +1,207 @@ +import QED + +#= ======== =# +# ActorPlasmaBuild # +#= ======== =# +Base.@kwdef mutable struct FUSEparameters__ActorPlasmaBuild{T<:Real} <: ParametersActor{T} + _parent::WeakRef = WeakRef(nothing) + _name::Symbol = :not_set + _time::Float64 = NaN + Δt::Entry{Float64} = Entry{Float64}("s", "Evolve for Δt") + Nt::Entry{Int} = Entry{Int}("-", "Number of time steps during evolution"; default=100) +end + +mutable struct ActorPlasmaBuild{D,P} <: SingleAbstractActor{D,P} + dd::IMAS.dd{D} + par::FUSEparameters__ActorPlasmaBuild{P} + QO::Union{Nothing,QED.QED_state} + build::Union{Nothing,QED.QED_build} +end + +""" + ActorPlasmaBuild(dd::IMAS.dd, act::ParametersAllActors; kw...) + +Evolves the plasma and coil/vessel currents using the QED current diffusion solver. + +!!! note + + This actor operates at "dd.global_time", any time advance must be done outside of the actor + + IMAS.new_timeslice!(dd, dd.global_time + Δt) + dd.global_time += Δt + ActorPlasmaBuild(dd, act) + +!!! note + + Stores data in `dd.core_profiles.profiles_1d[].j_ohmic` +""" +function ActorPlasmaBuild(dd::IMAS.dd, act::ParametersAllActors; kw...) + actor = ActorPlasmaBuild(dd, act.ActorPlasmaBuild; kw...) + step(actor) + finalize(actor) + return actor +end + +function ActorPlasmaBuild(dd::IMAS.dd, par::FUSEparameters__ActorPlasmaBuild; kw...) + logging_actor_init(ActorPlasmaBuild) + par = par(kw...) + return ActorPlasmaBuild(dd, par, nothing, nothing) +end + +function _step(actor::ActorPlasmaBuild) + dd = actor.dd + par = actor.par + + eqt = dd.equilibrium.time_slice[] + cp1d = dd.core_profiles.profiles_1d[] + + # non_inductive contribution + B0 = eqt.global_quantities.vacuum_toroidal_field.b0 + JBni = QED.FE(cp1d.grid.rho_tor_norm, cp1d.j_non_inductive .* B0) + + # initialize QED + if actor.QO === nothing + actor.QO = qed_init_from_imas(eqt, cp1d; uniform_rho = 501) + else + actor.QO.JBni = JBni + end + + actor.build = qed_build_from_imas(dd) + + # current diffusion + actor.QO = QED.evolve(actor.QO, η_imas(dd.core_profiles.profiles_1d[time0]), actor.build, par.Δt, par.Nt) + + return actor +end + +function _finalize(actor::ActorPlasmaBuild) + dd = actor.dd + + eqt = dd.equilibrium.time_slice[] + B0 = eqt.global_quantities.vacuum_toroidal_field.b0 + + cp1d = dd.core_profiles.profiles_1d[] + j_total = QED.JB(actor.QO; ρ=cp1d.grid.rho_tor_norm) ./ B0 + + if ismissing(cp1d, :j_non_inductive) + cp1d.j_ohmic = j_total + else + cp1d.j_ohmic = j_total .- cp1d.j_non_inductive + end + + # NEED TO POPULATE COIL CURRENTS BACK IN dd + build = actor.build + current_per_turn = build.Ic + + active_coils = dd.pf_active.coil + passive_loops = dd.pf_passive.loop + for (k, coil) in enumerate(active_coils) + IMAS.@ddtime(coil.current.data = current_per_turn[k]) + end + Nactive = length(active_coils) + for (k, loop) in enumerate(passive_loops) + IMAS.@ddtime(loop.current = current_per_turn[k + Nactive]) + end + + # MAYBE POPULATE RESISTANCE SO IT'S NOT RECALCULATED? + + return actor +end + +# utils +""" + qed_init_from_imas(dd::IMAS.dd) + +Setup QED build from data in IMAS `dd` +""" + +# function loopelement2coil(loop, element) +# outline = element.geometry.outline +# @assert length(outline.r) == 4 "For the time being passive structures must be composed of quadrilateral elements" +# passive_coil = VacuumFields.QuadCoil(outline.r, outline.z) +# element_turns = ismissing(element, :turns_with_sign) ? 1 : abs(element.turns_with_sign) +# loop_turns = sum(ismissing(elm, :turns_with_sign) ? 1 : abs(elm.turns_with_sign) for elm in loop.element) +# Ic = ismissing(loop, :current) ? 0.0 : loop.current / Nturns) +# VacuumFields.set_current!(passive_coil, Ic ) +# passive_coil.resistance = VacuumFields.resistance(passive_coil, loop.resistivity) +# return passive_coil +# end + +# elements(), turns(), QuadCoil(), and loop2multi() should all wind up in VacuumFields +elements(loop) = Iterators.filter(!isempty, loop.element) +turns(element::IMAS.pf_passive__loop___element) = isempty(element) ? 0.0 : (ismissing(element, :turns_with_sign) ? 1.0 : abs(element.turns_with_sign)) + +function VacuumFields.QuadCoil(elm::IMAS.pf_passive__loop___element, current_per_turn::Real = 0.0, resistance_per_turn::Real = 0.0) + R, Z = IMAS.outline(elm) + @assert length(R) == length(Z) == 4 + Nt = turns(elm) + return VacuumFields.QuadCoil(R, Z, current_per_turn * Nt, resistance_per_turn * Nt, Nt) +end + +function loop2multi(loop) + total_turns = sum(abs(ismissing(element, :turns_with_sign) ? 1.0 : element.turns_with_sign) for element in elements(loop)) + if !ismissing(loop, :resistance) + resistance_per_turn = loop.resistance / total_turns + else + resistance_per_turn = 0.0 + end + + # I'm assuming that pf_passive is like pf_passive and loop.current is current/turn like coil.current is + current_per_turn = ismissing(loop, :current) ? 0.0 : loop.current + + coils = [VacuumFields.QuadCoil(elm, current_per_turn, resistance_per_turn) for elm in elements(loop)] + if resistance_per_turn == 0.0 && !ismissing(loop, :resistivity) + for coil in coils + coil.resistance = VacuumFields.resistance(coil, loop.resistivity) + end + end + + return VacuumFields.MultiCoil(coils) +end + +function qed_build_from_imas(dd::IMAS.dd{D}) where {D <: Real} + + active_coils = VacuumFields.MultiCoils(dd); + passive_coils = [loopelement2coil(loop, element) for loop in dd.pf_passive.loop for element in loop.element] + coils = vcat(active_coils, passive_coils) + + # Coil-only quantities + Mcc = [VacuumFields.mutual(c1, c2) for c1 in coils, c2 in coils] + Ic = [VacuumFields.current(c) / VacuumFields.turns(c) for c in coils] #c urrent per turn + Rc = [VacuumFields.resistance(c) for c in coils]; + Vc = zero(Ic); + + eqt = dd.equilibrium.time_slice[] + cp1d = dd.core_profiles.profiles_1d[] + Ip = eqt.global_quantities.ip + + # plasma-coil mutuals + image = VacuumFields.Image(eqt) + Mpc = [VacuumFields.mutual(image, coil, Ip) for coil in coils] + @warn "ActorPlasmaBuild doesn't contain dMpc_dt yet" + dMpc_dt = zero(Mpc) # How Mpc changes in time (like shape)... to test later + + # self inductance + It = IMAS.cumtrapz(cp1d.grid.area, cp1d.j_tor) + Wp = 0.5 * IMAS.trapz(cp1d.grid.psi, It) + Li = 2 * Wp / Ip^2 # internal + ψb = eqt.profiles_1d.psi[end] + ψc = sum(Mpc[k] * Ic[k] for k in eachindex(coils)) + Le = (ψb - ψc) / Ip # external + Lp = Li + Le # total self + + # plasma resistance + # BCL 9/25/24: from Pohm, which may be wrong + Pohm = dd.core_sources.source[:ohmic].profiles_1d[].electrons.power_inside[end] + Ini = dd.core_profiles.global_quantities.current_non_inductive[] + Iohm = Ip - Ini + Rp = Pohm / (Ip * Iohm) + + # Waveforms + # These should come from pulse_schedule + W0 = QED.Waveform{Float64}(t -> 0.0) + V_waveforms = fill(W0, length(coils)); + + return QED.QED_build(Ic, Vc, Rc, Mcc, Vni, Rp, Lp, Mpc, dMpc_dt, V_waveforms) +end + From 90e52144ca91d1cc7d8bf7af44439441b6b006f4 Mon Sep 17 00:00:00 2001 From: Brendan Lyons Date: Thu, 26 Sep 2024 11:23:04 -0700 Subject: [PATCH 12/19] Rename to ActorQEDcoupled --- ...ed-build_actor.jl => qed-coupled_actor.jl} | 28 +++++++++---------- 1 file changed, 14 insertions(+), 14 deletions(-) rename src/actors/current/{qed-build_actor.jl => qed-coupled_actor.jl} (88%) diff --git a/src/actors/current/qed-build_actor.jl b/src/actors/current/qed-coupled_actor.jl similarity index 88% rename from src/actors/current/qed-build_actor.jl rename to src/actors/current/qed-coupled_actor.jl index cdb349305..2325e7d78 100644 --- a/src/actors/current/qed-build_actor.jl +++ b/src/actors/current/qed-coupled_actor.jl @@ -1,9 +1,9 @@ import QED #= ======== =# -# ActorPlasmaBuild # +# ActorQEDcoupled # #= ======== =# -Base.@kwdef mutable struct FUSEparameters__ActorPlasmaBuild{T<:Real} <: ParametersActor{T} +Base.@kwdef mutable struct FUSEparameters__ActorQEDcoupled{T<:Real} <: ParametersActor{T} _parent::WeakRef = WeakRef(nothing) _name::Symbol = :not_set _time::Float64 = NaN @@ -11,15 +11,15 @@ Base.@kwdef mutable struct FUSEparameters__ActorPlasmaBuild{T<:Real} <: Paramete Nt::Entry{Int} = Entry{Int}("-", "Number of time steps during evolution"; default=100) end -mutable struct ActorPlasmaBuild{D,P} <: SingleAbstractActor{D,P} +mutable struct ActorQEDcoupled{D,P} <: SingleAbstractActor{D,P} dd::IMAS.dd{D} - par::FUSEparameters__ActorPlasmaBuild{P} + par::FUSEparameters__ActorQEDcoupled{P} QO::Union{Nothing,QED.QED_state} build::Union{Nothing,QED.QED_build} end """ - ActorPlasmaBuild(dd::IMAS.dd, act::ParametersAllActors; kw...) + ActorQEDcoupled(dd::IMAS.dd, act::ParametersAllActors; kw...) Evolves the plasma and coil/vessel currents using the QED current diffusion solver. @@ -29,26 +29,26 @@ Evolves the plasma and coil/vessel currents using the QED current diffusion solv IMAS.new_timeslice!(dd, dd.global_time + Δt) dd.global_time += Δt - ActorPlasmaBuild(dd, act) + ActorQEDcoupled(dd, act) !!! note Stores data in `dd.core_profiles.profiles_1d[].j_ohmic` """ -function ActorPlasmaBuild(dd::IMAS.dd, act::ParametersAllActors; kw...) - actor = ActorPlasmaBuild(dd, act.ActorPlasmaBuild; kw...) +function ActorQEDcoupled(dd::IMAS.dd, act::ParametersAllActors; kw...) + actor = ActorQEDcoupled(dd, act.ActorQEDcoupled; kw...) step(actor) finalize(actor) return actor end -function ActorPlasmaBuild(dd::IMAS.dd, par::FUSEparameters__ActorPlasmaBuild; kw...) - logging_actor_init(ActorPlasmaBuild) +function ActorQEDcoupled(dd::IMAS.dd, par::FUSEparameters__ActorQEDcoupled; kw...) + logging_actor_init(ActorQEDcoupled) par = par(kw...) - return ActorPlasmaBuild(dd, par, nothing, nothing) + return ActorQEDcoupled(dd, par, nothing, nothing) end -function _step(actor::ActorPlasmaBuild) +function _step(actor::ActorQEDcoupled) dd = actor.dd par = actor.par @@ -74,7 +74,7 @@ function _step(actor::ActorPlasmaBuild) return actor end -function _finalize(actor::ActorPlasmaBuild) +function _finalize(actor::ActorQEDcoupled) dd = actor.dd eqt = dd.equilibrium.time_slice[] @@ -178,7 +178,7 @@ function qed_build_from_imas(dd::IMAS.dd{D}) where {D <: Real} # plasma-coil mutuals image = VacuumFields.Image(eqt) Mpc = [VacuumFields.mutual(image, coil, Ip) for coil in coils] - @warn "ActorPlasmaBuild doesn't contain dMpc_dt yet" + @warn "ActorQEDcoupled doesn't contain dMpc_dt yet" dMpc_dt = zero(Mpc) # How Mpc changes in time (like shape)... to test later # self inductance From 86609c7f36f24a76045e174b5102003e57601041 Mon Sep 17 00:00:00 2001 From: Brendan Lyons Date: Thu, 26 Sep 2024 13:57:40 -0700 Subject: [PATCH 13/19] Bug fixes --- src/actors/current/qed-coupled_actor.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/actors/current/qed-coupled_actor.jl b/src/actors/current/qed-coupled_actor.jl index 2325e7d78..bfebdcc57 100644 --- a/src/actors/current/qed-coupled_actor.jl +++ b/src/actors/current/qed-coupled_actor.jl @@ -69,6 +69,7 @@ function _step(actor::ActorQEDcoupled) actor.build = qed_build_from_imas(dd) # current diffusion + time0 = dd.global_time + 0.5 * par.Δt actor.QO = QED.evolve(actor.QO, η_imas(dd.core_profiles.profiles_1d[time0]), actor.build, par.Δt, par.Nt) return actor @@ -162,7 +163,7 @@ end function qed_build_from_imas(dd::IMAS.dd{D}) where {D <: Real} active_coils = VacuumFields.MultiCoils(dd); - passive_coils = [loopelement2coil(loop, element) for loop in dd.pf_passive.loop for element in loop.element] + passive_coils = [loop2multi(loop) for loop in dd.pf_passive.loop] coils = vcat(active_coils, passive_coils) # Coil-only quantities @@ -197,6 +198,9 @@ function qed_build_from_imas(dd::IMAS.dd{D}) where {D <: Real} Iohm = Ip - Ini Rp = Pohm / (Ip * Iohm) + # non-inductive voltage + Vni = Rp * Ini + # Waveforms # These should come from pulse_schedule W0 = QED.Waveform{Float64}(t -> 0.0) From 6084934c6c356ab1e42d7ed13b0a19117401929b Mon Sep 17 00:00:00 2001 From: Brendan Lyons Date: Thu, 26 Sep 2024 13:57:54 -0700 Subject: [PATCH 14/19] Expose ActorQEDcoupled in ActorCurrent --- src/FUSE.jl | 1 + src/actors/current/current_actor.jl | 6 ++++-- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/src/FUSE.jl b/src/FUSE.jl index b515d654b..b82e0d5ff 100644 --- a/src/FUSE.jl +++ b/src/FUSE.jl @@ -94,6 +94,7 @@ include(joinpath("actors", "nuclear", "blanket_actor.jl")) include(joinpath("actors", "nuclear", "neutronics_actor.jl")) include(joinpath("actors", "current", "qed_actor.jl")) +include(joinpath("actors", "current", "qed-coupled_actor.jl")) include(joinpath("actors", "current", "steadycurrent_actor.jl")) include(joinpath("actors", "current", "current_actor.jl")) diff --git a/src/actors/current/current_actor.jl b/src/actors/current/current_actor.jl index 45ac1839d..e0c01211e 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, :QEDcoupled, :none], "-", "Current actor to run"; default=:SteadyStateCurrent) 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) @@ -15,7 +15,7 @@ end mutable struct ActorCurrent{D,P} <: CompoundAbstractActor{D,P} dd::IMAS.dd{D} par::FUSEparameters__ActorCurrent{P} - jt_actor::Union{ActorSteadyStateCurrent{D,P},ActorQED{D,P},ActorNoOperation{D,P}} + jt_actor::Union{ActorSteadyStateCurrent{D,P},ActorQED{D,P},ActorQEDcoupled{D,P},ActorNoOperation{D,P}} end """ @@ -43,6 +43,8 @@ function ActorCurrent(dd::IMAS.dd, par::FUSEparameters__ActorCurrent, act::Param jt_actor = ActorSteadyStateCurrent(dd, act.ActorSteadyStateCurrent; par.ip_from, par.allow_floating_plasma_current) elseif par.model == :QED jt_actor = ActorQED(dd, act.ActorQED; par.ip_from, par.vloop_from, par.allow_floating_plasma_current) + elseif par.model == :QEDcoupled + jt_actor = ActorQEDcoupled(dd, act.ActorQEDcoupled) end return ActorCurrent(dd, par, jt_actor) end From 338b9c0ef2b4834c35599038d3974a73fe6a9b0e Mon Sep 17 00:00:00 2001 From: Brendan Lyons Date: Sun, 29 Sep 2024 15:28:21 -0700 Subject: [PATCH 15/19] Fix conventions in from_TokSys --- src/cases/_toksys.jl | 50 ++++++++++++++++++++++++++------------------ 1 file changed, 30 insertions(+), 20 deletions(-) diff --git a/src/cases/_toksys.jl b/src/cases/_toksys.jl index c0ec76864..27dc3f935 100644 --- a/src/cases/_toksys.jl +++ b/src/cases/_toksys.jl @@ -1,7 +1,7 @@ -function from_TokSys(dd::IMAS.dd, mat::Dict) - +function from_TokSys(dd::IMAS.dd, mat::Dict) + # ========== - + eqs = mat["eqs"] nr = Int(eqs["nr"][1,1]) @@ -14,9 +14,9 @@ function from_TokSys(dd::IMAS.dd, mat::Dict) b0 = Float64.(eqs["bzero"][1,:]) r0 = Float64.(eqs["rzero"][1,:][1]) - + rmaxis = Float64.(eqs["rmaxis"][1,:]) - zmaxis = Float64.(eqs["zmaxis"][1,:]) + zmaxis = Float64.(eqs["zmaxis"][1,:]) ffprim = Float64.(vcat(eqs["ffprim"][1, :]...)') pprime = Float64.(vcat(eqs["pprime"][1, :]...)') @@ -36,16 +36,22 @@ function from_TokSys(dd::IMAS.dd, mat::Dict) psirz = [Float64.(x) for x in eqs["psizr"][1, :]] # =========== - + pf_time = Float64.(mat["Vct"])[:,1] pf_voltages = Float64.(mat["Vcd"]) pf_names = mat["Tokamak"]["c"]["ccnames"] + iy_ic = dropdims(Int.(mat["Tokamak"]["c"]["iy"]["ic"]); dims=1) + pf_currents = Float64.(mat["yd"][:, iy_ic]) + # =========== dd.equilibrium.time = time dd.equilibrium.vacuum_toroidal_field.r0 = r0 dd.equilibrium.vacuum_toroidal_field.b0 = b0 + for k in eachindex(dd.equilibrium.time_slice) + IMAS.retime!(dd.equilibrium.time_slice[k], time[k]) + end resize!(dd.equilibrium.time_slice, length(time)) for (k,eqt) in enumerate(dd.equilibrium.time_slice) eq1d = eqt.profiles_1d @@ -60,32 +66,34 @@ function from_TokSys(dd::IMAS.dd, mat::Dict) eqt.global_quantities.ip = ip[k] + + # N.B. psi is poloidal flux, but pprime and ffprim are wrt poloidal flux per radian eq1d.psi = range(psimag[k], psibry[k], length(psi_norm)) eq1d.q = q[:,k] - eq1d.pressure = pres[:,k] - eq1d.dpressure_dpsi = pprime[:,k] + eq1d.pressure = pres[:,k] .+ 5e-4 .* pres[1,k] + eq1d.dpressure_dpsi = pprime[:,k] ./ 2π # make psi consistent eq1d.f = fpol[:,k] - eq1d.f_df_dpsi = ffprim[:,k] + eq1d.f_df_dpsi = ffprim[:,k] ./ 2π # make psi consistent eq2d.grid.dim1 = rg eq2d.grid.dim2 = zg eq2d.psi = collect(psirz[k]') end - + resize!(dd.wall.description_2d, 1) resize!(dd.wall.description_2d[1].limiter.unit, 1) dd.wall.description_2d[1].limiter.unit[1].outline.r = rlim dd.wall.description_2d[1].limiter.unit[1].outline.z = zlim - + dd.pulse_schedule.profiles_control.time = time dd.pulse_schedule.profiles_control.psi_norm = psi_norm - dd.pulse_schedule.profiles_control.dpressure_dpsi.reference = pprime - dd.pulse_schedule.profiles_control.f_df_dpsi.reference = ffprim - + dd.pulse_schedule.profiles_control.dpressure_dpsi.reference = pprime ./ 2π + dd.pulse_schedule.profiles_control.f_df_dpsi.reference = ffprim ./ 2π + dd.pulse_schedule.flux_control.time = time dd.pulse_schedule.flux_control.i_plasma.reference = ip - + dd.pulse_schedule.position_control.time = time dd.pulse_schedule.position_control.magnetic_axis.r.reference = rmaxis dd.pulse_schedule.position_control.magnetic_axis.z.reference = zmaxis @@ -94,14 +102,16 @@ function from_TokSys(dd::IMAS.dd, mat::Dict) dd.pulse_schedule.position_control.boundary_outline[k].r.reference = rbbbs[k,:] dd.pulse_schedule.position_control.boundary_outline[k].z.reference = zbbbs[k,:] end - + dd.pulse_schedule.pf_active.time = pf_time - resize!(dd.pulse_schedule.pf_active.coil, size(pf_voltages)[2]) + resize!(dd.pulse_schedule.pf_active.supply, size(pf_voltages)[2]) + # TokSys plasma current was the wrong sign, so coils are all the wrong sign. Flip voltage and current for (k,voltage) in enumerate(eachcol(pf_voltages)) - dd.pulse_schedule.pf_active.coil[k].identifier = pf_names[k] - dd.pulse_schedule.pf_active.coil[k].voltage.reference = voltage + dd.pulse_schedule.pf_active.supply[k].identifier = "S$k" + dd.pulse_schedule.pf_active.supply[k].voltage.reference = -voltage + dd.pulse_schedule.pf_active.supply[k].current.reference = -pf_currents[:,k] end - + IMAS.flux_surfaces(dd.equilibrium, IMAS.first_wall(dd.wall)...) end From 570c9cbc62fe67030716df199f9abddede30ffcb Mon Sep 17 00:00:00 2001 From: Brendan Lyons Date: Sun, 29 Sep 2024 15:31:14 -0700 Subject: [PATCH 16/19] Write FRESCO ff', f, and p to dd in finalize --- src/actors/equilibrium/fresco_actor.jl | 22 +++++++++++++++++----- 1 file changed, 17 insertions(+), 5 deletions(-) diff --git a/src/actors/equilibrium/fresco_actor.jl b/src/actors/equilibrium/fresco_actor.jl index c761c2213..05bee0b71 100644 --- a/src/actors/equilibrium/fresco_actor.jl +++ b/src/actors/equilibrium/fresco_actor.jl @@ -26,6 +26,7 @@ mutable struct ActorFRESCO{D,P} <: SingleAbstractActor{D,P} dd::IMAS.dd{D} par::FUSEparameters__ActorFRESCO{P} canvas::Union{Nothing,FRESCO.Canvas} + profile::Union{Nothing,FRESCO.CurrentProfile} end """ @@ -43,7 +44,7 @@ end function ActorFRESCO(dd::IMAS.dd{D}, par::FUSEparameters__ActorFRESCO{P}; kw...) where {D<:Real,P<:Real} logging_actor_init(ActorFRESCO) par = par(kw...) - return ActorFRESCO(dd, par, nothing) + return ActorFRESCO(dd, par, nothing, nothing) end """ @@ -75,11 +76,11 @@ function _step(actor::ActorFRESCO) gpp = IMAS.interp1d(eqt1d.psi_norm, eqt1d.dpressure_dpsi, :cubic) gffp = IMAS.interp1d(eqt1d.psi_norm, eqt1d.f_df_dpsi, :cubic) - profile = FRESCO.PprimeFFprime(x -> gpp(x), x -> gffp(x)) + actor.profile = FRESCO.PprimeFFprime(x -> gpp(x), x -> gffp(x)) - actor.canvas = FRESCO.Canvas(dd, par.nR, par.nZ) + actor.canvas = FRESCO.Canvas(dd, par.nR, par.nZ; load_pf_active=true, load_pf_passive=true) - FRESCO.solve!(actor.canvas, profile, par.number_of_iterations...; par.relax, par.debug, par.control, par.tolerance) + FRESCO.solve!(actor.canvas, actor.profile, par.number_of_iterations...; par.relax, par.debug, par.control, par.tolerance) # using Plots # p1 = Plots.heatmap(Rs, Zs, psi', aspect_ratio=:equal, xrange=(0,12.5), yrange=(-9,9), size=(400,500)) @@ -93,6 +94,7 @@ end # finalize by converting FRESCO canvas to dd.equilibrium function _finalize(actor::ActorFRESCO) canvas = actor.canvas + profile = actor.profile dd = actor.dd eq = dd.equilibrium eqt = eq.time_slice[] @@ -105,7 +107,17 @@ function _finalize(actor::ActorFRESCO) eqt.global_quantities.psi_boundary = canvas.Ψbnd eqt.global_quantities.psi_axis = canvas.Ψaxis eqt1d.psi = range(canvas.Ψaxis, canvas.Ψbnd, length(eqt1d.psi)) - # p, p', f, ff' don't change + # p' doesn't change + eqt1d.f_df_dpsi .*= profile.ffp_scale + + fend = dd.equilibrium.vacuum_toroidal_field.b0[] * dd.equilibrium.vacuum_toroidal_field.r0[] + f2 = 2 * IMAS.cumtrapz(eqt1d.psi, eqt1d.f_df_dpsi) + f2 .= f2 .- f2[end] .+ fend ^ 2 + eqt1d.f = sign(fend) .* sqrt.(f2) + + pend = eqt1d.pressure[end] + eqt1d.pressure = IMAS.cumtrapz(eqt1d.psi, eqt1d.dpressure_dpsi) + eqt1d.pressure .+= pend .- eqt1d.pressure[end] eq2d.grid_type.index = 1 eq2d.grid.dim1 = canvas.Rs From 435a061a466df9f377e6ea4d6fa6c30e2dce22fa Mon Sep 17 00:00:00 2001 From: Brendan Lyons Date: Sun, 29 Sep 2024 15:35:54 -0700 Subject: [PATCH 17/19] Equilibrium hack: initialize from pp' and FF' in pulse_schedule --- src/actors/equilibrium/equilibrium_actor.jl | 74 ++++++++++++++------- 1 file changed, 50 insertions(+), 24 deletions(-) diff --git a/src/actors/equilibrium/equilibrium_actor.jl b/src/actors/equilibrium/equilibrium_actor.jl index 5afadf265..3cdc4f5a2 100644 --- a/src/actors/equilibrium/equilibrium_actor.jl +++ b/src/actors/equilibrium/equilibrium_actor.jl @@ -9,7 +9,7 @@ Base.@kwdef mutable struct FUSEparameters__ActorEquilibrium{T<:Real} <: Paramete model::Switch{Symbol} = Switch{Symbol}([:TEQUILA, :FRESCO, :CHEASE], "-", "Equilibrium actor to run"; default=:TEQUILA) symmetrize::Entry{Bool} = Entry{Bool}("-", "Force equilibrium up-down symmetry with respect to magnetic axis"; default=false) #== data flow parameters ==# - j_p_from::Switch{Symbol} = Switch{Symbol}([:equilibrium, :core_profiles], "-", "Take j_tor and pressure profiles from this IDS"; default=:core_profiles) + j_p_from::Switch{Symbol} = Switch{Symbol}([:equilibrium, :core_profiles, :pulse_schedule_pp_ffp], "-", "Take j_tor and pressure profiles from this IDS"; default=:core_profiles) ip_from::Switch{Symbol} = switch_get_from(:ip) vacuum_r0_b0_from::Switch{Symbol} = switch_get_from(:vacuum_r0_b0; default=:pulse_schedule) #== display and debugging parameters ==# @@ -76,26 +76,29 @@ function _step(actor::ActorEquilibrium) if par.model == :FRESCO eqt = dd.equilibrium.time_slice[] idxeqt = IMAS.index(eqt) - if false && idxeqt != 1 && !ismissing(dd.equilibrium.time_slice[idxeqt-1].global_quantities, :ip) - eqt1d = eqt.profiles_1d - eqt1d__1 = dd.equilibrium.time_slice[idxeqt-1].profiles_1d - psi__1 = IMAS.interp1d(eqt1d__1.psi_norm, eqt1d__1.psi).(eqt1d.psi_norm) - R__1 = IMAS.interp1d(eqt1d__1.psi_norm, eqt1d__1.gm8).(eqt1d.psi_norm) - one_R__1 = IMAS.interp1d(eqt1d__1.psi_norm, eqt1d__1.gm9).(eqt1d.psi_norm) - one_R2__1 = IMAS.interp1d(eqt1d__1.psi_norm, eqt1d__1.gm1).(eqt1d.psi_norm) - b0 = eqt.global_quantities.vacuum_toroidal_field.b0 - r0 = eqt.global_quantities.vacuum_toroidal_field.r0 - dpressure_dpsi, f_df_dpsi, f = IMAS.calc_pprime_ffprim_f(psi__1, R__1, one_R__1, one_R2__1, r0, b0; eqt1d.pressure, eqt1d.j_tor) - eqt1d.dpressure_dpsi = dpressure_dpsi - eqt1d.f_df_dpsi = f_df_dpsi - eqt1d.f = f - else - # NOTE: using free_boundary=true will also update pf_active - # which is then used as initial guess by FRESCO - ActorTEQUILA(dd, actor.act; free_boundary=true) - fw = IMAS.first_wall(dd.wall) - IMAS.flux_surfaces(eqt, fw.r, fw.z) + if false + if false && idxeqt != 1 && !ismissing(dd.equilibrium.time_slice[idxeqt-1].global_quantities, :ip) + eqt1d = eqt.profiles_1d + eqt1d__1 = dd.equilibrium.time_slice[idxeqt-1].profiles_1d + psi__1 = IMAS.interp1d(eqt1d__1.psi_norm, eqt1d__1.psi).(eqt1d.psi_norm) + R__1 = IMAS.interp1d(eqt1d__1.psi_norm, eqt1d__1.gm8).(eqt1d.psi_norm) + one_R__1 = IMAS.interp1d(eqt1d__1.psi_norm, eqt1d__1.gm9).(eqt1d.psi_norm) + one_R2__1 = IMAS.interp1d(eqt1d__1.psi_norm, eqt1d__1.gm1).(eqt1d.psi_norm) + b0 = eqt.global_quantities.vacuum_toroidal_field.b0 + r0 = eqt.global_quantities.vacuum_toroidal_field.r0 + dpressure_dpsi, f_df_dpsi, f = IMAS.calc_pprime_ffprim_f(psi__1, R__1, one_R__1, one_R2__1, r0, b0; eqt1d.pressure, eqt1d.j_tor) + eqt1d.dpressure_dpsi = dpressure_dpsi + eqt1d.f_df_dpsi = f_df_dpsi + eqt1d.f = f + else + # NOTE: using free_boundary=true will also update pf_active + # which is then used as initial guess by FRESCO + ActorTEQUILA(dd, actor.act; free_boundary=true) + fw = IMAS.first_wall(dd.wall) + IMAS.flux_surfaces(eqt, fw.r, fw.z) + end end + end # step selected equilibrium actor @@ -190,8 +193,15 @@ function prepare(actor::ActorEquilibrium) pressure0 = eqt1d.pressure j_itp = IMAS.interp1d(rho_pol_norm_sqrt0, j_tor0, :cubic) p_itp = IMAS.interp1d(rho_pol_norm_sqrt0, pressure0, :cubic) + elseif par.j_p_from == :pulse_schedule_pp_ffp + @assert !isempty(dd.equilibrium.time) + eqt1d = dd.equilibrium.time_slice[].profiles_1d + psi0 = eqt1d.psi + rho_tor_norm0 = eqt1d.rho_tor_norm + rho_pol_norm_sqrt0 = sqrt.(eqt1d.psi_norm) + pend = eqt1d.pressure[end] else - @assert par.j_p_from in (:core_profiles, :equilibrium) + @assert par.j_p_from in (:core_profiles, :equilibrium, :pulse_schedule_pp_ffp) end # get ip and b0 before wiping eqt in case ip_from=:equilibrium @@ -252,12 +262,28 @@ function prepare(actor::ActorEquilibrium) end end - # set j_tor and pressure, forcing zero derivative on axis eqt1d = dd.equilibrium.time_slice[].profiles_1d eqt1d.psi = psi0 eqt1d.rho_tor_norm = rho_tor_norm0 - eqt1d.j_tor = j_itp.(sqrt.(eqt1d.psi_norm)) - eqt1d.pressure = p_itp.(sqrt.(eqt1d.psi_norm)) + if par.j_p_from === :pulse_schedule_pp_ffp + time0 = dd.global_time + profcon = dd.pulse_schedule.profiles_control + psin_ps = profcon.psi_norm + + pprime = IMAS.interp1d(psin_ps, IMAS.get_time_array(profcon.dpressure_dpsi, :reference, [time0])[:,1], :cubic) + ffprime = IMAS.interp1d(psin_ps, IMAS.get_time_array(profcon.f_df_dpsi, :reference, [time0])[:,1], :cubic) + #fpol = IMAS.interp1d(psin_ps, IMAS.get_time_array(profcon.f, :reference, [time0])[:,1], :cubic) + + psin_eq = eqt1d.psi_norm + eqt1d.dpressure_dpsi = pprime.(psin_eq) + eqt1d.pressure = IMAS.cumtrapz(eqt1d.psi, eqt1d.dpressure_dpsi) + eqt1d.pressure .+= pend .- eqt1d.pressure[end] + eqt1d.f_df_dpsi = ffprime.(psin_eq) + else + # set j_tor and pressure, forcing zero derivative on axis + eqt1d.j_tor = j_itp.(sqrt.(eqt1d.psi_norm)) + eqt1d.pressure = p_itp.(sqrt.(eqt1d.psi_norm)) + end return dd end From a891abb6492eafcd65331a9db9bf60ec8e233cc0 Mon Sep 17 00:00:00 2001 From: Brendan Lyons Date: Sun, 29 Sep 2024 15:39:45 -0700 Subject: [PATCH 18/19] QEDcoupled: use coil voltages from pulse_schedule --- Project.toml | 2 +- src/actors/current/qed-coupled_actor.jl | 107 ++++++++++++++---------- 2 files changed, 64 insertions(+), 45 deletions(-) diff --git a/Project.toml b/Project.toml index c372d5cf0..fd24d5032 100644 --- a/Project.toml +++ b/Project.toml @@ -110,6 +110,6 @@ TGLFNN = "1" TJLF = "1" ThermalSystemModels = "1" TimerOutputs = "0.5" -VacuumFields = "1.2, 2" +VacuumFields = "2.3" Weave = "0.10" YAML = "0.4" diff --git a/src/actors/current/qed-coupled_actor.jl b/src/actors/current/qed-coupled_actor.jl index bfebdcc57..73ecff907 100644 --- a/src/actors/current/qed-coupled_actor.jl +++ b/src/actors/current/qed-coupled_actor.jl @@ -9,6 +9,8 @@ Base.@kwdef mutable struct FUSEparameters__ActorQEDcoupled{T<:Real} <: Parameter _time::Float64 = NaN Δt::Entry{Float64} = Entry{Float64}("s", "Evolve for Δt") Nt::Entry{Int} = Entry{Int}("-", "Number of time steps during evolution"; default=100) + #== display and debugging parameters ==# + debug::Entry{Bool} = Entry{Bool}("-", "Turn on QED debugging/plotting"; default=false) end mutable struct ActorQEDcoupled{D,P} <: SingleAbstractActor{D,P} @@ -61,16 +63,16 @@ function _step(actor::ActorQEDcoupled) # initialize QED if actor.QO === nothing - actor.QO = qed_init_from_imas(eqt, cp1d; uniform_rho = 501) + actor.QO = qed_init_from_imas(eqt, cp1d; uniform_rho = 33) else actor.QO.JBni = JBni end - actor.build = qed_build_from_imas(dd) + actor.build = qed_build_from_imas(dd, dd.global_time - par.Δt) # current diffusion - time0 = dd.global_time + 0.5 * par.Δt - actor.QO = QED.evolve(actor.QO, η_imas(dd.core_profiles.profiles_1d[time0]), actor.build, par.Δt, par.Nt) + time0 = dd.global_time - 0.5 * par.Δt + actor.QO = QED.evolve(actor.QO, η_imas(dd.core_profiles.profiles_1d[time0]), actor.build, par.Δt, par.Nt; par.debug) return actor end @@ -129,42 +131,10 @@ Setup QED build from data in IMAS `dd` # end # elements(), turns(), QuadCoil(), and loop2multi() should all wind up in VacuumFields -elements(loop) = Iterators.filter(!isempty, loop.element) -turns(element::IMAS.pf_passive__loop___element) = isempty(element) ? 0.0 : (ismissing(element, :turns_with_sign) ? 1.0 : abs(element.turns_with_sign)) - -function VacuumFields.QuadCoil(elm::IMAS.pf_passive__loop___element, current_per_turn::Real = 0.0, resistance_per_turn::Real = 0.0) - R, Z = IMAS.outline(elm) - @assert length(R) == length(Z) == 4 - Nt = turns(elm) - return VacuumFields.QuadCoil(R, Z, current_per_turn * Nt, resistance_per_turn * Nt, Nt) -end - -function loop2multi(loop) - total_turns = sum(abs(ismissing(element, :turns_with_sign) ? 1.0 : element.turns_with_sign) for element in elements(loop)) - if !ismissing(loop, :resistance) - resistance_per_turn = loop.resistance / total_turns - else - resistance_per_turn = 0.0 - end - - # I'm assuming that pf_passive is like pf_passive and loop.current is current/turn like coil.current is - current_per_turn = ismissing(loop, :current) ? 0.0 : loop.current - - coils = [VacuumFields.QuadCoil(elm, current_per_turn, resistance_per_turn) for elm in elements(loop)] - if resistance_per_turn == 0.0 && !ismissing(loop, :resistivity) - for coil in coils - coil.resistance = VacuumFields.resistance(coil, loop.resistivity) - end - end - return VacuumFields.MultiCoil(coils) -end - -function qed_build_from_imas(dd::IMAS.dd{D}) where {D <: Real} +function qed_build_from_imas(dd::IMAS.dd{D}, time0::D) where {D <: Real} - active_coils = VacuumFields.MultiCoils(dd); - passive_coils = [loop2multi(loop) for loop in dd.pf_passive.loop] - coils = vcat(active_coils, passive_coils) + coils = VacuumFields.MultiCoils(dd; load_pf_active=true, load_pf_passive=true); # Coil-only quantities Mcc = [VacuumFields.mutual(c1, c2) for c1 in coils, c2 in coils] @@ -172,8 +142,8 @@ function qed_build_from_imas(dd::IMAS.dd{D}) where {D <: Real} Rc = [VacuumFields.resistance(c) for c in coils]; Vc = zero(Ic); - eqt = dd.equilibrium.time_slice[] - cp1d = dd.core_profiles.profiles_1d[] + eqt = dd.equilibrium.time_slice[time0] + cp1d = dd.core_profiles.profiles_1d[time0] Ip = eqt.global_quantities.ip # plasma-coil mutuals @@ -194,7 +164,7 @@ function qed_build_from_imas(dd::IMAS.dd{D}) where {D <: Real} # plasma resistance # BCL 9/25/24: from Pohm, which may be wrong Pohm = dd.core_sources.source[:ohmic].profiles_1d[].electrons.power_inside[end] - Ini = dd.core_profiles.global_quantities.current_non_inductive[] + Ini = IMAS.get_time_array(dd.core_profiles.global_quantities, :current_non_inductive, time0, :linear) Iohm = Ip - Ini Rp = Pohm / (Ip * Iohm) @@ -203,9 +173,58 @@ function qed_build_from_imas(dd::IMAS.dd{D}) where {D <: Real} # Waveforms # These should come from pulse_schedule - W0 = QED.Waveform{Float64}(t -> 0.0) - V_waveforms = fill(W0, length(coils)); - + #V_waveforms = fill(W0, length(coils)); + #W = QED.Waveform{D}(t -> dd.pulse_schedule.voltage[1](t)) + + Nc = length(coils) + V_waveforms = Vector{QED.Waveform{D}}(undef, length(coils)) + Vactive = Vwaveforms_from_pulse_schedule(dd, time0) + @assert length(Vactive) == length(dd.pf_active.coil) + Nactive = length(Vactive) + V_waveforms[1:Nactive] .= Vactive + V_waveforms[Nactive+1:end] .= fill(QED.Waveform{D}(t -> 0.0), Nc - Nactive) return QED.QED_build(Ic, Vc, Rc, Mcc, Vni, Rp, Lp, Mpc, dMpc_dt, V_waveforms) end +function Vwaveforms_from_pulse_schedule(dd::IMAS.dd{D}, t_start::D=IMAS.global_time(dd)) where {D<:Real} + Nc = length(dd.pf_active.coil) + time = dd.pulse_schedule.pf_active.time .- t_start + WFs = Vector{QED.Waveform{D}}(undef, Nc) + for (k, supply) in enumerate(dd.pulse_schedule.pf_active.supply) + if k in (1, 2) + WFs[k] = QED.Waveform(time, supply.voltage.reference) + elseif k == 3 + Cu = VacuumFields.MultiCoil(dd.pf_active.coil[3]) + Cl = VacuumFields.MultiCoil(dd.pf_active.coil[4]) + Lu = VacuumFields.mutual(Cu, Cu) + Ll = VacuumFields.mutual(Cl, Cl) + M = VacuumFields.mutual(Cu, Cl) + Lt = (Lu + 2M + Ll) + + facu = (Lu + M) / Lt + WFs[3] = QED.Waveform(time , facu .* supply.voltage.reference) + + facl = (Ll + M) / Lt + WFs[4] = QED.Waveform(time , facl .* supply.voltage.reference) + elseif k < 12 + coil = dd.pf_active.coil[k+1] + WFs[k+1] = QED.Waveform(time, supply.voltage.reference) + else + Cu = VacuumFields.MultiCoil(dd.pf_active.coil[13]) + Cl = VacuumFields.MultiCoil(dd.pf_active.coil[14]) + Lu = VacuumFields.mutual(Cu, Cu) + Ll = VacuumFields.mutual(Cl, Cl) + M = -VacuumFields.mutual(Cu, Cl) # oppositely connected + Lt = (Lu + 2M + Ll) + + # positive voltage gives negative current in VS3U + facu = -(Lu + M) / Lt + WFs[13] = QED.Waveform(time , facu .* supply.voltage.reference) + + # positive voltage gives positive current in VS3U + facl = (Ll + M) / Lt + WFs[14] = QED.Waveform(time , facl .* supply.voltage.reference) + end + end + return WFs +end From 6f04b5739b81e8fa4a379da8324270e37403877b Mon Sep 17 00:00:00 2001 From: Orso Meneghini Date: Mon, 30 Sep 2024 00:13:24 -0700 Subject: [PATCH 19/19] from_TokSys without retime! --- src/cases/_toksys.jl | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/src/cases/_toksys.jl b/src/cases/_toksys.jl index 866810e45..4e21bd9dd 100644 --- a/src/cases/_toksys.jl +++ b/src/cases/_toksys.jl @@ -61,11 +61,8 @@ function from_TokSys(dd::IMAS.dd, mat::Dict; what...) dd.equilibrium.time = eq_time dd.equilibrium.vacuum_toroidal_field.r0 = r0 dd.equilibrium.vacuum_toroidal_field.b0 = b0 - for k in eachindex(dd.equilibrium.time_slice) - IMAS.retime!(dd.equilibrium.time_slice[k], eq_time[k]) - end - resize!(dd.equilibrium.time_slice, length(eq_time)) - for (k, eqt) in enumerate(dd.equilibrium.time_slice) + for (k, time) in enumerate(collect(eq_time)) + eqt = resize!(dd.equilibrium.time_slice, time) eq1d = eqt.profiles_1d eq2d = resize!(eqt.profiles_2d, 1)[1] eq2d.grid_type.index = 1