diff --git a/.github/workflows/runtests.yml b/.github/workflows/runtests.yml index a55c1ebcd..c64190c2a 100644 --- a/.github/workflows/runtests.yml +++ b/.github/workflows/runtests.yml @@ -35,7 +35,7 @@ jobs: using Pkg Pkg.activate(".") dependencies = PackageSpec[] - for package in ["IMAS", "IMASDD", "CoordinateConventions", "MillerExtendedHarmonic", "FusionMaterials", "VacuumFields", "Equilibrium", "TAUENN", "EPEDNN", "TGLFNN", "QED", "FiniteElementHermite", "CHEASE", "Fortran90Namelists", "EFIT"] + for package in ["IMAS", "IMASDD", "CoordinateConventions", "MillerExtendedHarmonic", "FusionMaterials", "VacuumFields", "Equilibrium", "TAUENN", "EPEDNN", "TGLFNN", "QED", "FiniteElementHermite", "Fortran90Namelists", "CHEASE", "EFIT", "NNeutronics"] push!(dependencies, PackageSpec(url="https://project-torrey-pines:${{secrets.PTP_READ_TOKEN}}@github.com/ProjectTorreyPines/$package.jl.git")) end Pkg.add(dependencies) diff --git a/Makefile b/Makefile index 4300c2eea..d1d04d136 100644 --- a/Makefile +++ b/Makefile @@ -25,9 +25,9 @@ install_no_registry: julia -e '\ using Pkg;\ Pkg.activate(".");\ -Pkg.develop(["IMAS", "IMASDD", "CoordinateConventions", "MillerExtendedHarmonic", "FusionMaterials", "VacuumFields", "Equilibrium", "TAUENN", "EPEDNN", "TGLFNN", "QED", "FiniteElementHermite", "Fortran90Namelists", "CHEASE", "EFIT"]);\ +Pkg.develop(["IMAS", "IMASDD", "CoordinateConventions", "MillerExtendedHarmonic", "FusionMaterials", "VacuumFields", "Equilibrium", "TAUENN", "EPEDNN", "TGLFNN", "QED", "FiniteElementHermite", "Fortran90Namelists", "CHEASE", "EFIT", "NNeutronics"]);\ Pkg.activate();\ -Pkg.develop(["FUSE", "IMAS", "IMASDD", "CoordinateConventions", "MillerExtendedHarmonic", "FusionMaterials", "VacuumFields", "Equilibrium", "TAUENN", "EPEDNN", "TGLFNN", "QED", "FiniteElementHermite", "Fortran90Namelists", "CHEASE", "EFIT"]);\ +Pkg.develop(["FUSE", "IMAS", "IMASDD", "CoordinateConventions", "MillerExtendedHarmonic", "FusionMaterials", "VacuumFields", "Equilibrium", "TAUENN", "EPEDNN", "TGLFNN", "QED", "FiniteElementHermite", "Fortran90Namelists", "CHEASE", "EFIT", "NNeutronics"]);\ ' install: clone_update_all install_no_registry precompile @@ -61,7 +61,7 @@ precompile: julia -e 'using Pkg; Pkg.activate("."); Pkg.precompile()' clone_update_all: - make -j 100 FUSE IMAS IMASDD CoordinateConventions MillerExtendedHarmonic FusionMaterials VacuumFields Equilibrium TAUENN EPEDNN TGLFNN QED FiniteElementHermite Fortran90Namelists CHEASE EFIT + make -j 100 FUSE IMAS IMASDD CoordinateConventions MillerExtendedHarmonic FusionMaterials VacuumFields Equilibrium TAUENN EPEDNN TGLFNN QED FiniteElementHermite Fortran90Namelists CHEASE EFIT NNeutronics update: install clone_update_all precompile @@ -101,16 +101,19 @@ EPEDNN: QED: $(call clone_update_repo,$@) -CHEASE: +FiniteElementHermite: $(call clone_update_repo,$@) -EFIT: +Fortran90Namelists: $(call clone_update_repo,$@) -Fortran90Namelists: +CHEASE: $(call clone_update_repo,$@) -FiniteElementHermite: +EFIT: + $(call clone_update_repo,$@) + +NNeutronics: $(call clone_update_repo,$@) docker_image: diff --git a/Project.toml b/Project.toml index 1912ca37f..3d36c27d2 100644 --- a/Project.toml +++ b/Project.toml @@ -21,6 +21,7 @@ FiniteElementHermite = "6ce10167-d368-4c0e-af34-9787cdd175e5" Fortran90Namelists = "8fb689aa-71ff-4044-8071-0cffc910b57d" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" FusionMaterials = "4c86da02-02c8-4634-8460-96566129f8e0" +GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f" IMAS = "13ead8c1-b7d1-41bb-a6d0-5b8b65ed587a" IMASDD = "06b86afa-9f21-11ec-2ef8-e51b8960cfc5" Images = "916415d5-f1e6-5110-898d-aaa5f9f070e0" @@ -33,6 +34,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LoggingExtras = "e6f89c97-d47a-5376-807f-9c37f3926c36" Metaheuristics = "bcdb8e00-2c21-11e9-3065-2b553b22f898" MillerExtendedHarmonic = "c82744c2-dc08-461a-8c37-87ab04d0f9b8" +NNeutronics = "a9424c20-d414-11ec-167b-9106c24d956c" NumericalIntegration = "e7bfaba1-d571-5449-8927-abc22e82249b" Optim = "429524aa-4258-5aef-a3af-852621145aeb" OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" diff --git a/cases/ARC.jl b/cases/ARC.jl index a2e401e39..f08dc86e4 100644 --- a/cases/ARC.jl +++ b/cases/ARC.jl @@ -3,7 +3,7 @@ CFS/MIT ARC design """ -function case_parameters(::Type{Val{:ARC}})::Tuple{ParametersAllInits, ParametersAllActors} +function case_parameters(::Type{Val{:ARC}})::Tuple{ParametersAllInits,ParametersAllActors} ini = ParametersAllInits() act = ParametersAllActors() ini.general.casename = "ARC" @@ -16,7 +16,7 @@ function case_parameters(::Type{Val{:ARC}})::Tuple{ParametersAllInits, Parameter ini.equilibrium.B0 = -11.5 ini.equilibrium.Z0 = 0.0 ini.equilibrium.ip = 9.9e6 - ini.equilibrium.βn = 1.4 + ini.equilibrium.pressure_core = 1.45e6 ini.equilibrium.x_point = (3.1, 1.85) ini.equilibrium.symmetric = true act.ActorCXbuild.rebuild_wall = false @@ -50,9 +50,12 @@ function case_parameters(::Type{Val{:ARC}})::Tuple{ParametersAllInits, Parameter ini.tf.n_coils = 18 ini.tf.technology = coil_technology(:HTS) ini.oh.technology = coil_technology(:HTS) - ini.oh.flattop_duration = 1800 - ini.core_profiles.ne_ped = 7e19 #estimate (from ITER) + #ini.target.power_electric_net = 50E6 ? + ini.target.flattop_duration = 1800 + #ini.target.tritium_breeding_ratio = 1.0 ? + + ini.core_profiles.ne_ped = 1.0e20 ini.core_profiles.greenwald_fraction = 0.49 ini.core_profiles.helium_fraction = 0.10 #estimate ini.core_profiles.T_shaping = 1.8 #estimate (from ITER) @@ -65,6 +68,7 @@ function case_parameters(::Type{Val{:ARC}})::Tuple{ParametersAllInits, Parameter ini.ic_antennas.power_launched = 4 * 1e6 #rf power coupled act.ActorPFcoilsOpt.symmetric = true #note: symmetric, but not evenly spaced + # act.ActorEquilibrium.model = :CHEASE return set_new_base!(ini), set_new_base!(act) end diff --git a/cases/CAT.jl b/cases/CAT.jl index 8d0f1d525..1841e2ff5 100644 --- a/cases/CAT.jl +++ b/cases/CAT.jl @@ -3,7 +3,7 @@ GA Compact Advanced Tokamak design """ -function case_parameters(::Type{Val{:CAT}})::Tuple{ParametersAllInits, ParametersAllActors} +function case_parameters(::Type{Val{:CAT}})::Tuple{ParametersAllInits,ParametersAllActors} ini = ParametersAllInits() act = ParametersAllActors() @@ -29,7 +29,6 @@ function case_parameters(::Type{Val{:CAT}})::Tuple{ParametersAllInits, Parameter ini.tf.technology = coil_technology(:ITER, :TF) ini.oh.technology = coil_technology(:ITER, :OH) - ini.oh.flattop_duration = 1000 ini.core_profiles.ne_ped = 7E19 ini.core_profiles.greenwald_fraction = 0.8 ini.core_profiles.helium_fraction = 0.01 @@ -40,6 +39,8 @@ function case_parameters(::Type{Val{:CAT}})::Tuple{ParametersAllInits, Parameter ini.core_profiles.bulk = :DT ini.core_profiles.impurity = :Ne + ini.target.flattop_duration = 1000 + ini.nbi.power_launched = 20E6 ini.nbi.beam_energy = 200e3 ini.nbi.beam_mass = 2 diff --git a/cases/D3D.jl b/cases/D3D.jl index 781210175..3a41081e5 100644 --- a/cases/D3D.jl +++ b/cases/D3D.jl @@ -3,7 +3,7 @@ DIII-D """ -function case_parameters(::Type{Val{:D3D}})::Tuple{ParametersAllInits, ParametersAllActors} +function case_parameters(::Type{Val{:D3D}})::Tuple{ParametersAllInits,ParametersAllActors} ini = ParametersAllInits() act = ParametersAllActors() @@ -27,7 +27,6 @@ function case_parameters(::Type{Val{:D3D}})::Tuple{ParametersAllInits, Parameter ini.tf.technology = coil_technology(:copper) ini.oh.technology = coil_technology(:copper) - ini.oh.flattop_duration = 5 ini.core_profiles.ne_ped = 5E19 ini.core_profiles.greenwald_fraction = 0.7 @@ -44,6 +43,8 @@ function case_parameters(::Type{Val{:D3D}})::Tuple{ParametersAllInits, Parameter ini.nbi.beam_mass = 2 ini.nbi.toroidal_angle = 20.0 / 180 * pi + ini.target.flattop_duration = 5 + act.ActorPFcoilsOpt.symmetric = true return set_new_base!(ini), set_new_base!(act) diff --git a/cases/FPP.jl b/cases/FPP.jl index 111c854f9..1f2a15f25 100644 --- a/cases/FPP.jl +++ b/cases/FPP.jl @@ -32,6 +32,9 @@ function case_parameters(::Type{Val{:FPP}}; version::Symbol, init_from::Symbol): act.ActorHFSsizing.fixed_aspect_ratio = true end + ini.target.tritium_breeding_ratio = 1.1 + ini.target.cost = 5E9 + ini.core_profiles.bulk = :DT ini.core_profiles.rot_core = 0.0 ini.tf.shape = :princeton_D_scaled @@ -42,7 +45,7 @@ function case_parameters(::Type{Val{:FPP}}; version::Symbol, init_from::Symbol): if init_from == :ods ini.pf_active.n_pf_coils_outside = 8 else - ini.pf_active.n_pf_coils_outside = 6 + ini.pf_active.n_pf_coils_outside = 5 end ini.material.shield = "Tungsten" @@ -66,10 +69,26 @@ function case_parameters(::Type{Val{:FPP}}; version::Symbol, init_from::Symbol): # greenwald_fraction is a powerful knob ini.core_profiles.greenwald_fraction = 0.9 - # ini.equilibrium.δ *= -1 # negative triangularity + # negative triangularity + # ini.equilibrium.δ *= -1 - # ini.equilibrium.ζ = 0.1 # squareness + # squareness + # ini.equilibrium.ζ = 0.1 # act.ActorEquilibrium.model = :CHEASE + # add wall layer + if true + gasc_add_wall_layers!(ini.build.layers; thickness=0.02) + if version != :v1 + ini.build.n_first_wall_conformal_layers = 2 + end + end + + # bucking + ini.center_stack.bucked = false + if ini.center_stack.bucked + gasc_buck_OH_TF!(ini.build.layers) + end + return ini, act end diff --git a/cases/HDB5.jl b/cases/HDB5.jl index 04a3ce30b..963b5335e 100644 --- a/cases/HDB5.jl +++ b/cases/HDB5.jl @@ -6,7 +6,7 @@ import CSV For description of cases/variables see https://osf.io/593q6/ """ -function case_parameters(::Type{Val{:HDB5}}; tokamak::Union{String,Symbol}=:any, case=missing, database_case=missing)::Tuple{ParametersAllInits, ParametersAllActors} +function case_parameters(::Type{Val{:HDB5}}; tokamak::Union{String,Symbol}=:any, case=missing, database_case=missing)::Tuple{ParametersAllInits,ParametersAllActors} if !ismissing(database_case) data_row = load_hdb5(database_case=database_case) elseif !ismissing(case) @@ -30,9 +30,8 @@ function case_parameters(data_row::DataFrames.DataFrameRow) ini.equilibrium.ϵ = data_row[:AMIN] / data_row[:RGEO] ini.equilibrium.κ = data_row[:KAPPA] ini.equilibrium.δ = data_row[:DELTA] - ini.equilibrium.βn = 1.0 ini.equilibrium.ip = data_row[:IP] - + ini.equilibrium.pressure_core = IMAS.pressure_avg_from_beta_n(1.0, data_row[:AMIN], data_row[:BT], data_row[:IP]) * 3.0 act.ActorSolovev.area = data_row[:AREA] act.ActorSolovev.volume = data_row[:VOL] @@ -62,8 +61,8 @@ function case_parameters(data_row::DataFrames.DataFrameRow) # Core_profiles parameters ini.core_profiles.ne_ped = data_row[:NEL] / 1.3 - ini.core_profiles.greenwald_fraction = data_row[:NEL]*1e-20 / (data_row[:IP]/1e6 / (pi * data_row[:AMIN]^2 )) - ini.core_profiles.helium_fraction = 0. + ini.core_profiles.greenwald_fraction = data_row[:NEL] * 1e-20 / (data_row[:IP] / 1e6 / (pi * data_row[:AMIN]^2)) + ini.core_profiles.helium_fraction = 0.0 ini.core_profiles.T_shaping = 1.8 ini.core_profiles.w_ped = 0.03 ini.core_profiles.zeff = data_row[:ZEFF] @@ -98,10 +97,10 @@ end function load_hdb5(tokamak::Union{String,Symbol}=:all; maximum_ohmic_fraction::Float64=0.25, database_case::Union{Int,Missing}=missing, extra_signal_names=Union{String,Symbol}[]) # For description of variables see https://osf.io/593q6/ run_df = CSV.read(joinpath(@__DIR__, "..", "sample", "HDB5_compressed.csv"), DataFrames.DataFrame) - run_df[:,"database_case"] = collect(1:length(run_df[:,"TOK"])) + run_df[:, "database_case"] = collect(1:length(run_df[:, "TOK"])) if !ismissing(database_case) - return run_df[run_df.database_case .== database_case, :] + return run_df[run_df.database_case.==database_case, :] end signal_names = ["TOK", "SHOT", "AMIN", "KAPPA", "DELTA", "NEL", "ZEFF", "TAUTH", "RGEO", "BT", "IP", "PNBI", "ENBI", "PICRH", "PECRH", "POHM", "MEFF", "VOL", "AREA", "WTH", "CONFIG"] @@ -113,8 +112,8 @@ function load_hdb5(tokamak::Union{String,Symbol}=:all; maximum_ohmic_fraction::F # some basic filters run_df = run_df[(run_df.TOK.!="T10").&(run_df.TOK.!="TDEV").&(run_df.KAPPA.>1.0).&(run_df.DELTA.<0.79).&(1.6 .< run_df.MEFF .< 2.2).&(1.1 .< run_df.ZEFF .< 5.9), :] # Filter cases where the ohmic power is dominating - run_df[:,"Paux"] = run_df[:,"PNBI"] .+ run_df[:,"PECRH"] .+ run_df[:,"PICRH"] .+ run_df[:,"POHM"] - run_df = run_df[run_df[:,"POHM"] .< maximum_ohmic_fraction .* (run_df[:,"Paux"] .- run_df[:,"POHM"]),:] + run_df[:, "Paux"] = run_df[:, "PNBI"] .+ run_df[:, "PECRH"] .+ run_df[:, "PICRH"] .+ run_df[:, "POHM"] + run_df = run_df[run_df[:, "POHM"]. "BOP_gasc", "index" => 2) sys.power = 0.07 .* bop_thermal.power_electric_generated - else + + elseif model == :EU_DEMO # More realistic DEMO numbers bop_systems = [:cryostat, :tritium_handling, :pumping, :pf_active] # index 2 : 5 for (idx, system) in enumerate(bop_systems) sys = resize!(bop_electric.system, "name" => string(system), "index" => (idx + 1)) sys.power = electricity(system, bop.time) end + else + error("ActorBalanceOfPlant: model = $(model) not recognized") end return actor end @@ -109,7 +119,7 @@ function electricity(nbi::IMAS.nbi, time_array::Vector{<:Real}) end function electricity(ec_launchers::IMAS.ec_launchers, time_array::Vector{<:Real}) - return heating_and_current_drive_calc(ec_launchers.launcher, time_array) + return heating_and_current_drive_calc(ec_launchers.beam, time_array) end function electricity(ic_antennas::IMAS.ic_antennas, time_array::Vector{<:Real}) diff --git a/src/actors/blanket_actors.jl b/src/actors/blanket_actors.jl index 52444644d..dfefae2fb 100644 --- a/src/actors/blanket_actors.jl +++ b/src/actors/blanket_actors.jl @@ -2,8 +2,11 @@ # ActorBlanket # #= ============ =# -Base.@kwdef mutable struct ActorBlanket <: ReactorAbstractActor +import NNeutronics + +mutable struct ActorBlanket <: ReactorAbstractActor dd::IMAS.dd + par::ParametersActor blanket_multiplier::Real thermal_power_extraction_efficiency::Real end @@ -15,7 +18,7 @@ function ParametersActor(::Type{Val{:ActorBlanket}}) Real, "", "Fraction of thermal power that is carried out by the coolant at the blanket interface, rather than being lost in the surrounding strutures."; - default=1.0, + default=1.0 ) return par end @@ -30,16 +33,20 @@ Blanket actor """ function ActorBlanket(dd::IMAS.dd, act::ParametersAllActors; kw...) par = act.ActorBlanket(kw...) - actor = ActorBlanket(dd, par.blanket_multiplier, par.thermal_power_extraction_efficiency) + actor = ActorBlanket(dd, par) step(actor) finalize(actor) return actor end +function ActorBlanket(dd::IMAS.dd, par::ParametersActor; kw...) + par = par(kw...) + return ActorBlanket(dd, par, par.blanket_multiplier, par.thermal_power_extraction_efficiency) +end + function step(actor::ActorBlanket) dd = actor.dd - empty!(dd.blanket) eqt = dd.equilibrium.time_slice[] nnt = dd.neutronics.time_slice[] @@ -49,15 +56,14 @@ function step(actor::ActorBlanket) total_power_radiated = 0.0 # IMAS.radiative_power(dd.core_profiles.profiles_1d[]) blankets = [structure for structure in dd.build.structure if structure.type == Int(_blanket_)] - - tritium_breeding_ratio = 0.0 resize!(dd.blanket.module, length(blankets)) + for (k, structure) in enumerate(blankets) bm = dd.blanket.module[k] bm.name = structure.name # evaluate neutron_capture_fraction - tmp = convex_hull(collect(zip(vcat(eqt.boundary.outline.r, structure.outline.r), vcat(eqt.boundary.outline.z, structure.outline.z))); closed_polygon=true) + tmp = convex_hull(vcat(eqt.boundary.outline.r, structure.outline.r), vcat(eqt.boundary.outline.z, structure.outline.z); closed_polygon=true) index = findall(x -> x == 1, [IMAS.PolygonOps.inpolygon((r, z), tmp) for (r, z) in zip(dd.neutronics.first_wall.r, dd.neutronics.first_wall.z)]) neutron_capture_fraction = sum(nnt.wall_loading.power[index]) / total_power_neutrons @@ -75,9 +81,60 @@ function step(actor::ActorBlanket) bmt.power_thermal_extracted = actor.thermal_power_extraction_efficiency * (bmt.power_thermal_neutrons + bmt.power_thermal_radiated) - bmt.tritium_breeding_ratio = 1.0 # some function - tritium_breeding_ratio += bmt.tritium_breeding_ratio * bmt.power_incident_neutrons + # blanket layer structure (designed to handle missing wall and/or missing shield) + resize!(bm.layer, 3) + if sum(structure.outline.r) / length(structure.outline.r) < eqt.boundary.geometric_axis.r # HFS + d1 = IMAS.get_build(dd.build, type=_wall_, fs=_hfs_, raise_error_on_missing=false) + d1 = missing + d2 = IMAS.get_build(dd.build, type=_blanket_, fs=_hfs_) + d3 = IMAS.get_build(dd.build, type=_shield_, fs=_hfs_, return_only_one=false, raise_error_on_missing=false) + if d3 !== missing + d3 = d3[end] + end + d3 = missing + else # LFS + d1 = IMAS.get_build(dd.build, type=_wall_, fs=_lfs_) + d2 = IMAS.get_build(dd.build, type=_blanket_, fs=_lfs_) + d3 = IMAS.get_build(dd.build, type=_shield_, fs=_lfs_, return_only_one=false, raise_error_on_missing=false) + if d3 !== missing + d3 = d3[1] + end + end + for (kl, dl) in enumerate(bm.layer) + dl.name = "dummy layer $kl" + dl.thickness = 0.0 + dl.material = "vacuum" + end + for (kl, dl) in enumerate([d1, d2, d3]) + if dl !== missing + for field in [:name, :thickness, :material] + setproperty!(bm.layer[kl], field, getproperty(dl, field)) + end + end + end end - @ddtime(dd.blanket.tritium_breeding_ratio = tritium_breeding_ratio / total_power_neutrons) + # Optimize Li6/Li7 ratio to obtain target TBR + function target_TBR(blanket_model, Li6, dd, target=nothing) + total_tritium_breeding_ratio = 0.0 + for bm in dd.blanket.module + bmt = bm.time_slice[] + bm.layer[2].material = @sprintf("lithium-lead: Li6/7=%3.3f", Li6) + bmt.tritium_breeding_ratio = NNeutronics.TBR(blanket_model, [dl.thickness for dl in bm.layer]..., Li6) + total_tritium_breeding_ratio += bmt.tritium_breeding_ratio * bmt.power_incident_neutrons / total_power_neutrons + end + if target === nothing + return total_tritium_breeding_ratio + else + return (total_tritium_breeding_ratio - target)^2 + end + end + + blanket_model_1d = NNeutronics.Blanket() + res = Optim.optimize(Li6 -> target_TBR(blanket_model_1d, Li6, dd, dd.target.tritium_breeding_ratio), 0.0, 100.0, Optim.GoldenSection(), rel_tol=1E-6) + total_tritium_breeding_ratio = target_TBR(blanket_model_1d, res.minimizer, dd) + + @ddtime(dd.blanket.tritium_breeding_ratio = total_tritium_breeding_ratio) + + return actor end diff --git a/src/actors/build_actors.jl b/src/actors/build_actors.jl index c9f78873c..12ca307d0 100644 --- a/src/actors/build_actors.jl +++ b/src/actors/build_actors.jl @@ -44,7 +44,7 @@ function oh_required_J_B!(bd::IMAS.build; double_swing::Bool=true) innerSolenoidRadius = OH.start_radius outerSolenoidRadius = OH.end_radius - totalOhFluxReq = bd.flux_swing_estimates.rampup + bd.flux_swing_estimates.flattop + bd.flux_swing_estimates.pf + totalOhFluxReq = bd.flux_swing.rampup + bd.flux_swing.flattop + bd.flux_swing.pf # Calculate magnetic field at solenoid bore required to match flux swing request RiRo_factor = innerSolenoidRadius / outerSolenoidRadius @@ -58,11 +58,11 @@ function oh_required_J_B!(bd::IMAS.build; double_swing::Bool=true) end """ - flattop_estimate!(bd::IMAS.build, eqt::IMAS.equilibrium__time_slice, cp1d::IMAS.core_profiles__profiles_1d; double_swing::Bool=true) + flattop_duration!(bd::IMAS.build, eqt::IMAS.equilibrium__time_slice, cp1d::IMAS.core_profiles__profiles_1d; double_swing::Bool=true) Estimate OH flux requirement during flattop (if j_ohmic profile is missing then steady state ohmic profile is assumed) """ -function flattop_estimate!(bd::IMAS.build, eqt::IMAS.equilibrium__time_slice, cp1d::IMAS.core_profiles__profiles_1d; double_swing::Bool=true) +function flattop_duration!(bd::IMAS.build, eqt::IMAS.equilibrium__time_slice, cp1d::IMAS.core_profiles__profiles_1d; double_swing::Bool=true) OH = IMAS.get_build(bd, type=_oh_) innerSolenoidRadius = OH.start_radius outerSolenoidRadius = OH.end_radius @@ -70,8 +70,8 @@ function flattop_estimate!(bd::IMAS.build, eqt::IMAS.equilibrium__time_slice, cp # estimate oh flattop flux and duration RiRo_factor = innerSolenoidRadius / outerSolenoidRadius totalOhFlux = bd.oh.max_b_field * (pi * outerSolenoidRadius^2 * (RiRo_factor^2 + RiRo_factor + 1.0) * (double_swing ? 2 : 1)) / 3.0 - bd.flux_swing_estimates.flattop = totalOhFlux - bd.flux_swing_estimates.rampup - bd.flux_swing_estimates.pf - bd.oh.flattop_estimate = bd.flux_swing_estimates.flattop / abs(integrate(cp1d.grid.area, cp1d.j_ohmic ./ cp1d.conductivity_parallel)) + bd.flux_swing.flattop = totalOhFlux - bd.flux_swing.rampup - bd.flux_swing.pf + bd.oh.flattop_duration = bd.flux_swing.flattop / abs(integrate(cp1d.grid.area, cp1d.j_ohmic ./ cp1d.conductivity_parallel)) end #= ========= =# @@ -127,8 +127,9 @@ end #= ========== =# # flux-swing # #= ========== =# -Base.@kwdef mutable struct ActorFluxSwing <: ReactorAbstractActor +mutable struct ActorFluxSwing <: ReactorAbstractActor dd::IMAS.dd + par::ParametersActor operate_at_j_crit::Bool j_tolerance::Real end @@ -162,17 +163,22 @@ OH flux consumption based on: * vertical field from PF coils !!! note - Stores data in `dd.build.flux_swing_estimates`, `dd.build.tf`, and `dd.build.oh` + Stores data in `dd.build.flux_swing`, `dd.build.tf`, and `dd.build.oh` """ function ActorFluxSwing(dd::IMAS.dd, act::ParametersAllActors; kw...) par = act.ActorFluxSwing(kw...) - actor = ActorFluxSwing(; dd, par...) + actor = ActorFluxSwing(dd, par) step(actor) finalize(actor) return actor end +function ActorFluxSwing(dd::IMAS.dd, par::ParametersActor; kw...) + par = par(kw...) + return ActorFluxSwing(dd, par, par.operate_at_j_crit, par.j_tolerance) +end + """ step(actor::ActorFluxSwing; operate_at_j_crit::Bool=actor.operate_at_j_crit, j_tolerance::Real=actor.j_tolerance, only=:all) @@ -186,25 +192,26 @@ The `only` parameter controls if :tf, :oh, or :all (both) should be calculated function step(actor::ActorFluxSwing; operate_at_j_crit::Bool=actor.operate_at_j_crit, j_tolerance::Real=actor.j_tolerance, only=:all) bd = actor.dd.build + target = actor.dd.target eq = actor.dd.equilibrium eqt = eq.time_slice[] cp = actor.dd.core_profiles cp1d = cp.profiles_1d[] if only ∈ [:all, :oh] - bd.flux_swing_estimates.rampup = rampup_flux_estimates(eqt, cp) - bd.flux_swing_estimates.pf = pf_flux_estimates(eqt) + bd.flux_swing.rampup = rampup_flux_estimates(eqt, cp) + bd.flux_swing.pf = pf_flux_estimates(eqt) if operate_at_j_crit oh_maximum_J_B!(bd; j_tolerance) - bd.flux_swing_estimates.flattop = flattop_flux_estimates(bd) # target flattop flux based on available current + bd.flux_swing.flattop = flattop_flux_estimates(bd) # flattop flux based on available current else - bd.flux_swing_estimates.flattop = flattop_flux_estimates(bd, cp1d) # target flattop flux based on target duration + bd.flux_swing.flattop = flattop_flux_estimates(target, cp1d) # flattop flux based on target duration oh_required_J_B!(bd) end - # estimate flattop duration - flattop_estimate!(bd, eqt, cp1d) + # flattop duration + flattop_duration!(bd, eqt, cp1d) end if only ∈ [:all, :tf] @@ -244,12 +251,12 @@ function rampup_flux_estimates(eqt::IMAS.equilibrium__time_slice, cp::IMAS.core_ end """ - flattop_flux_estimates(bd::IMAS.build, cp1d::IMAS.core_profiles__profiles_1d) + flattop_flux_estimates(target::IMAS.target, cp1d::IMAS.core_profiles__profiles_1d) Estimate OH flux requirement during flattop (if j_ohmic profile is missing then steady state ohmic profile is assumed) """ -function flattop_flux_estimates(bd::IMAS.build, cp1d::IMAS.core_profiles__profiles_1d) - return abs(integrate(cp1d.grid.area, cp1d.j_ohmic ./ cp1d.conductivity_parallel)) * bd.oh.flattop_duration # V*s +function flattop_flux_estimates(target::IMAS.target, cp1d::IMAS.core_profiles__profiles_1d) + return abs(integrate(cp1d.grid.area, cp1d.j_ohmic ./ cp1d.conductivity_parallel)) * target.flattop_duration # V*s end """ @@ -264,7 +271,7 @@ function flattop_flux_estimates(bd::IMAS.build; double_swing::Bool=true) magneticFieldSolenoidBore = bd.oh.max_b_field RiRo_factor = innerSolenoidRadius / outerSolenoidRadius totalOhFluxReq = magneticFieldSolenoidBore / 3.0 * pi * outerSolenoidRadius^2 * (RiRo_factor^2 + RiRo_factor + 1.0) * (double_swing ? 2 : 1) - bd.flux_swing_estimates.flattop = totalOhFluxReq - bd.flux_swing_estimates.rampup - bd.flux_swing_estimates.pf + bd.flux_swing.flattop = totalOhFluxReq - bd.flux_swing.rampup - bd.flux_swing.pf end """ @@ -292,8 +299,13 @@ end #= ========== =# # LFS sizing # #= ========== =# -Base.@kwdef mutable struct ActorLFSsizing <: ReactorAbstractActor +mutable struct ActorLFSsizing <: ReactorAbstractActor dd::IMAS.dd + par::ParametersActor + function ActorLFSsizing(dd::IMAS.dd, par::ParametersActor; kw...) + par = par(kw...) + return new(dd, par) + end end function ParametersActor(::Type{Val{:ActorLFSsizing}}) @@ -315,10 +327,10 @@ Actor that resizes the Low Field Side of the build. """ function ActorLFSsizing(dd::IMAS.dd, act::ParametersAllActors; kw...) par = act.ActorLFSsizing(kw...) + actor = ActorLFSsizing(dd, par) if par.do_plot plot(dd.build) end - actor = ActorLFSsizing(dd) step(actor; par.verbose) finalize(actor) if par.do_plot @@ -360,7 +372,9 @@ end #= ========== =# # HFS sizing # #= ========== =# -Base.@kwdef mutable struct ActorHFSsizing <: ReactorAbstractActor +mutable struct ActorHFSsizing <: ReactorAbstractActor + dd::IMAS.dd + par::ParametersActor stresses_actor::ActorStresses fluxswing_actor::ActorFluxSwing end @@ -386,15 +400,13 @@ Actor that resizes the High Field Side of the build. !!! note Manipulates radial build information in `dd.build.layer` """ -function ActorHFSsizing(dd::IMAS.dd, act::ParametersAllActors; kw...) +function ActorHFSsizing(dd::IMAS.dd, act::ParametersAllActors; kw_ActorFluxSwing=Dict(), kw_ActorStresses=Dict(), kw...) par = act.ActorHFSsizing(kw...) + actor = ActorHFSsizing(dd, par, act; kw_ActorFluxSwing, kw_ActorStresses) if par.do_plot p = plot(dd.build) end - fluxswing_actor = ActorFluxSwing(; dd, operate_at_j_crit=par.unconstrained_flattop_duration, par.j_tolerance) - stresses_actor = ActorStresses(dd) - actor = ActorHFSsizing(stresses_actor, fluxswing_actor) - step(actor; par.verbose, par.j_tolerance, par.stress_tolerance, par.fixed_aspect_ratio, par.unconstrained_flattop_duration) + step(actor) finalize(actor) if par.do_plot display(plot!(p, dd.build; cx=false)) @@ -402,14 +414,21 @@ function ActorHFSsizing(dd::IMAS.dd, act::ParametersAllActors; kw...) return actor end +function ActorHFSsizing(dd::IMAS.dd, par::ParametersActor, act::ParametersAllActors; kw_ActorFluxSwing=Dict(), kw_ActorStresses=Dict(), kw...) + par = act.ActorHFSsizing(kw...) + fluxswing_actor = ActorFluxSwing(dd, act.ActorFluxSwing; kw_ActorFluxSwing...) + stresses_actor = ActorStresses(dd, act.ActorStresses; kw_ActorStresses...) + return ActorHFSsizing(dd, par, stresses_actor, fluxswing_actor) +end + function step( actor::ActorHFSsizing; - j_tolerance::Real=0.4, - stress_tolerance::Real=0.2, - fixed_aspect_ratio::Bool=true, - unconstrained_flattop_duration::Bool=true, - verbose::Bool=false, - do_plot=false + j_tolerance::Real=actor.par.j_tolerance, + stress_tolerance::Real=actor.par.stress_tolerance, + fixed_aspect_ratio::Bool=actor.par.fixed_aspect_ratio, + unconstrained_flattop_duration::Bool=actor.par.unconstrained_flattop_duration, + verbose::Bool=actor.par.verbose, + debug_plot=false ) function target_value(value, target, tolerance) # relative error with tolerance @@ -478,7 +497,7 @@ function step( c_stf = target_value(maximum(dd.solid_mechanics.center_stack.stress.vonmises.tf), stainless_steel.yield_strength, stress_tolerance) end - if do_plot + if debug_plot push!(C_JOH, c_joh) push!(C_SOH, c_soh) push!(C_JTF, c_jtf) @@ -506,7 +525,7 @@ function step( old_a = plasma.thickness / 2.0 old_ϵ = old_R0 / old_a - if do_plot + if debug_plot C_JOH = [] C_SOH = [] C_JTF = [] @@ -575,7 +594,7 @@ function step( a = plasma.thickness / 2.0 ϵ = R0 / a - if do_plot + if debug_plot p = plot(yscale=:log10) plot!(p, C_JOH ./ (C_JOH .> 0.0), label="Jcrit OH") plot!(p, C_SOH ./ (C_SOH .> 0.0), label="Stresses OH") @@ -589,8 +608,8 @@ function step( @show target_B0 @show dd.build.tf.max_b_field * TFhfs.end_radius / R0 - @show dd.build.oh.flattop_estimate @show dd.build.oh.flattop_duration + @show dd.target.flattop_duration @show dd.build.oh.max_j @show dd.build.oh.critical_j @@ -620,7 +639,7 @@ function step( @assert maximum(dd.solid_mechanics.center_stack.stress.vonmises.oh) < stainless_steel.yield_strength @assert maximum(dd.solid_mechanics.center_stack.stress.vonmises.tf) < stainless_steel.yield_strength if !unconstrained_flattop_duration - @assert rel_error(dd.build.oh.flattop_estimate, dd.build.oh.flattop_duration) < 0.1 "Relative error on flattop duration is more than 10% ($(dd.build.oh.flattop_estimate) --> $(dd.build.oh.flattop_duration))" + @assert rel_error(dd.build.oh.flattop_duration, dd.target.flattop_duration) < 0.1 "Relative error on flattop duration is more than 10% ($(dd.build.oh.flattop_duration) --> $(dd.target.flattop_duration))" end if fixed_aspect_ratio @assert rel_error(ϵ, old_ϵ) < 0.1 "ActorHFSsizing: plasma aspect ratio changed more than 10% ($old_ϵ --> $ϵ)" @@ -632,8 +651,13 @@ end #= ============= =# # cross-section # #= ============= =# -Base.@kwdef mutable struct ActorCXbuild <: ReactorAbstractActor +mutable struct ActorCXbuild <: ReactorAbstractActor dd::IMAS.dd + par::ParametersActor + function ActorCXbuild(dd::IMAS.dd, par::ParametersActor; kw...) + par = par(kw...) + return new(dd, par) + end end function ParametersActor(::Type{Val{:ActorCXbuild}}) @@ -653,8 +677,8 @@ Actor that builds the 2D cross section of the build. """ function ActorCXbuild(dd::IMAS.dd, act::ParametersAllActors; kw...) par = act.ActorCXbuild(kw...) - actor = ActorCXbuild(dd) - step(actor; par.rebuild_wall) + actor = ActorCXbuild(dd, par) + step(actor) finalize(actor) if par.do_plot plot(dd.build) @@ -663,7 +687,7 @@ function ActorCXbuild(dd::IMAS.dd, act::ParametersAllActors; kw...) return actor end -function step(actor::ActorCXbuild; rebuild_wall::Bool=true) +function step(actor::ActorCXbuild; rebuild_wall::Bool=actor.par.rebuild_wall) build_cx!(actor.dd; rebuild_wall) end @@ -689,8 +713,8 @@ function wall_from_eq(bd::IMAS.build, eqt::IMAS.equilibrium__time_slice; diverto wall_poly = LibGEOS.buffer(plasma_poly, a) # hfs and lfs spacing may not be symmetric - R = [v[1] for v in LibGEOS.coordinates(wall_poly)[1]] - Z = [v[2] for v in LibGEOS.coordinates(wall_poly)[1]] + R = [v[1] for v in GeoInterface.coordinates(wall_poly)[1]] + Z = [v[2] for v in GeoInterface.coordinates(wall_poly)[1]] R .+= ((R_lfs_plasma + R_hfs_plasma) - (maximum(rlcfs) + minimum(rlcfs))) / 2.0 R[R.R_lfs_plasma] .= R_lfs_plasma @@ -700,7 +724,8 @@ function wall_from_eq(bd::IMAS.build, eqt::IMAS.equilibrium__time_slice; diverto t = LinRange(0, 2pi, 31) # divertor lengths - max_divertor_length = (maximum(zlcfs) - minimum(zlcfs)) * divertor_length_fraction + linear_plasma_size = sqrt((maximum(zlcfs) - minimum(zlcfs)) * (maximum(rlcfs) - minimum(rlcfs))) + max_divertor_length = linear_plasma_size * divertor_length_fraction # private flux regions private = IMAS.flux_surface(eqt, ψb, false) @@ -708,7 +733,7 @@ function wall_from_eq(bd::IMAS.build, eqt::IMAS.equilibrium__time_slice; diverto if sign(pz[1] - Z0) != sign(pz[end] - Z0) # open flux surface does not encicle the plasma continue - elseif IMAS.minimum_distance_two_shapes(pr, pz, rlcfs, zlcfs) > (maximum(zlcfs) - minimum(zlcfs)) / 20 + elseif IMAS.minimum_distance_two_shapes(pr, pz, rlcfs, zlcfs) > linear_plasma_size / 5 # secondary Xpoint far away continue elseif (sum(pz) - Z0) < 0 @@ -754,8 +779,8 @@ function wall_from_eq(bd::IMAS.build, eqt::IMAS.equilibrium__time_slice; diverto wall_poly = LibGEOS.buffer(wall_poly, -a / 4) wall_poly = LibGEOS.buffer(wall_poly, a / 4) - pr = [v[1] for v in LibGEOS.coordinates(wall_poly)[1]] - pz = [v[2] for v in LibGEOS.coordinates(wall_poly)[1]] + pr = [v[1] for v in GeoInterface.coordinates(wall_poly)[1]] + pz = [v[2] for v in GeoInterface.coordinates(wall_poly)[1]] pr, pz = IMAS.resample_2d_line(pr, pz; step=0.1) @@ -769,6 +794,15 @@ function divertor_regions!(bd::IMAS.build, eqt::IMAS.equilibrium__time_slice) ipl = IMAS.get_build(bd, type=_plasma_, return_index=true) plasma_poly = xy_polygon(bd.layer[ipl]) + wall_poly = xy_polygon(bd.layer[ipl-1]) + for ltype in [_blanket_, _shield_, _wall_,] + iwl = IMAS.get_build(bd, type=ltype, fs=_hfs_, return_index=true, raise_error_on_missing=false) + if iwl !== missing + wall_poly = xy_polygon(bd.layer[iwl]) + break + end + end + divertors = IMAS.IDSvectorElement[] for x_point in eqt.boundary.x_point Zx = x_point.z @@ -780,12 +814,11 @@ function divertor_regions!(bd::IMAS.build, eqt::IMAS.equilibrium__time_slice) pz = vcat(yy, [Zx * 5, Zx * 5], yy[1]) domain = xy_polygon(pr, pz) - wall_poly = xy_polygon(bd.layer[ipl-1]) divertor_poly = LibGEOS.intersection(wall_poly, domain) divertor_poly = LibGEOS.difference(divertor_poly, plasma_poly) - pr = [v[1] for v in LibGEOS.coordinates(divertor_poly)[1]] - pz = [v[2] for v in LibGEOS.coordinates(divertor_poly)[1]] + pr = [v[1] for v in GeoInterface.coordinates(divertor_poly)[1]] + pz = [v[2] for v in GeoInterface.coordinates(divertor_poly)[1]] # assign to build structure if Zx > Z0 @@ -810,7 +843,7 @@ function blanket_regions!(bd::IMAS.build, eqt::IMAS.equilibrium__time_slice) layers = bd.layer iblanket = IMAS.get_build(bd; type=_blanket_, fs=_lfs_, return_index=true, raise_error_on_missing=false) - if iblanket === nothing + if iblanket === missing return IMAS.IDSvectorElement[] end layer = layers[iblanket] @@ -827,7 +860,7 @@ function blanket_regions!(bd::IMAS.build, eqt::IMAS.equilibrium__time_slice) geometries = LibGEOS.getGeometries(ring_poly) blankets = IMAS.IDSvectorElement[] for poly in geometries - coords = LibGEOS.coordinates(poly) + coords = GeoInterface.coordinates(poly) pr = [v[1] for v in coords[1]] pz = [v[2] for v in coords[1]] @@ -872,6 +905,7 @@ function build_cx!(dd::IMAS.dd; rebuild_wall::Bool=true) build_cx!(dd.build, pr, pz) divertor_regions!(dd.build, dd.equilibrium.time_slice[]) + blanket_regions!(dd.build, dd.equilibrium.time_slice[]) if wall === missing || rebuild_wall @@ -911,18 +945,20 @@ function build_cx!(bd::IMAS.build, pr::Vector{Float64}, pz::Vector{Float64}) # k-1 means the layer outside (ie. towards the tf) # k is the current layer # k+1 means the layer inside (ie. towards the plasma) + # # forward pass: from plasma to TF _convex_hull_ and then desired TF shape tf_to_plasma = IMAS.get_build(bd, fs=_hfs_, return_only_one=false, return_index=true) plasma_to_tf = reverse(tf_to_plasma) for k in plasma_to_tf - layer_shape = BuildLayerShape(mod(mod(bd.layer[k].shape, 1000), 100)) + original_shape = bd.layer[k].shape if k == itf + 1 # layer that is inside of the TF sets TF shape + layer_shape = BuildLayerShape(mod(mod(bd.layer[k].shape, 1000), 100)) optimize_shape(bd, k + 1, k, layer_shape; tight=!coils_inside) else # everything else is conformal convex hull optimize_shape(bd, k + 1, k, _convex_hull_) - bd.layer[k].shape = Int(layer_shape) + bd.layer[k].shape = original_shape end end # reverse pass: from TF to plasma only with negative offset @@ -931,6 +967,15 @@ function build_cx!(bd::IMAS.build, pr::Vector{Float64}, pz::Vector{Float64}) optimize_shape(bd, k, k + 1, _negative_offset_) end end + # forward pass: from plasma to TF with desired shapes + for k in plasma_to_tf[1:end-1] + if bd.layer[k].shape == Int(_negative_offset_) + break + else + layer_shape = BuildLayerShape(mod(mod(bd.layer[k].shape, 1000), 100)) + optimize_shape(bd, k + 1, k, layer_shape) + end + end # _in_ D = minimum(IMAS.get_build(bd, type=_tf_, fs=_hfs_).outline.z) @@ -1013,11 +1058,10 @@ function optimize_shape(bd::IMAS.build, obstr_index::Int, layer_index::Int, shap # handle offset, negative offset, negative offset, and convex-hull if layer.shape in [Int(_offset_), Int(_negative_offset_), Int(_convex_hull_)] poly = LibGEOS.buffer(xy_polygon(oR, oZ), (hfs_thickness + lfs_thickness) / 2.0) - R = [v[1] .+ r_offset for v in LibGEOS.coordinates(poly)[1]] - Z = [v[2] for v in LibGEOS.coordinates(poly)[1]] + R = [v[1] .+ r_offset for v in GeoInterface.coordinates(poly)[1]] + Z = [v[2] for v in GeoInterface.coordinates(poly)[1]] if layer.shape == Int(_convex_hull_) - h = [[r, z] for (r, z) in collect(zip(R, Z))] - hull = convex_hull(h; closed_polygon=true) + hull = convex_hull(R, Z; closed_polygon=true) R = [r for (r, z) in hull] Z = [z for (r, z) in hull] # resample disabled because this can lead to outlines of different layers to be crossing @@ -1075,6 +1119,8 @@ function optimize_shape(bd::IMAS.build, obstr_index::Int, layer_index::Int, shap layer.shape_parameters = optimize_shape(oR, oZ, target_minimum_distance, func, l_start, l_end, layer.shape_parameters) layer.outline.r, layer.outline.z = func(l_start, l_end, layer.shape_parameters...; resample=false) end + + IMAS.reorder_flux_surface!(layer.outline.r, layer.outline.z) # display(plot!(layer.outline.r, layer.outline.z)) end diff --git a/src/actors/chease_actor.jl b/src/actors/chease_actor.jl new file mode 100644 index 000000000..cb85c2c15 --- /dev/null +++ b/src/actors/chease_actor.jl @@ -0,0 +1,157 @@ +import CHEASE + +#= =========== =# +# ActorCHEASE # +#= =========== =# +mutable struct ActorCHEASE <: PlasmaAbstractActor + dd::IMAS.dd + par::ParametersActor + chease::Union{Nothing,CHEASE.Chease} +end + +function ParametersActor(::Type{Val{:ActorCHEASE}}) + par = ParametersActor(nothing) + par.free_boundary = Entry(Bool, "", "Convert fixed boundary equilibrium to free boundary one"; default=true) + par.clear_workdir = Entry(Bool, "", "Clean the temporary workdir for CHEASE"; default=true) + par.rescale_eq_to_ip = Entry(Bool, "", "Scale equilibrium to match Ip"; default=false) + return par +end + +""" + ActorCHEASE(dd::IMAS.dd, act::ParametersAllActors; kw...) + +This actor runs the Fixed boundary equilibrium solver CHEASE +""" +function ActorCHEASE(dd::IMAS.dd, act::ParametersAllActors; kw...) + par = act.ActorCHEASE(kw...) + actor = ActorCHEASE(dd, par) + step(actor) + finalize(actor) + return actor +end + +function ActorCHEASE(dd::IMAS.dd, par::ParametersActor; kw...) + par = par(kw...) + ActorCHEASE(dd, par, nothing) +end + +""" + prepare(dd::IMAS.dd, :ActorCHEASE, act::ParametersAllActors; kw...) + +Prepare dd to run ActorCHEASE +* Copy pressure from core_profiles to equilibrium +* Copy j_parallel from core_profiles to equilibrium +""" +function prepare(dd::IMAS.dd, ::Type{Val{:ActorCHEASE}}, act::ParametersAllActors; kw...) + eq1d = dd.equilibrium.time_slice[].profiles_1d + cp1d = dd.core_profiles.profiles_1d[] + eq1d.j_tor = IMAS.interp1d(cp1d.grid.psi_norm, cp1d.j_tor).(eq1d.psi_norm) + eq1d.pressure = IMAS.interp1d(cp1d.grid.psi_norm, cp1d.pressure).(eq1d.psi_norm) +end + +""" + step(actor::ActorCHEASE) + +Runs CHEASE on the r_z boundary, equilibrium pressure and equilibrium j_tor +""" +function step(actor::ActorCHEASE) + dd = actor.dd + eqt = dd.equilibrium.time_slice[] + eq1d = eqt.profiles_1d + + # remove points at high curvature points (ie. X-points) + r_bound = eqt.boundary.outline.r + z_bound = eqt.boundary.outline.z + r_bound, z_bound = IMAS.resample_2d_line(r_bound, z_bound; n_points=201) + index = abs.(IMAS.curvature(r_bound, z_bound)) .< 0.9 + r_bound = r_bound[index] + z_bound = z_bound[index] + + # scalars + Ip = eqt.global_quantities.ip + Bt_center = @ddtime(dd.equilibrium.vacuum_toroidal_field.b0) + r_center = dd.equilibrium.vacuum_toroidal_field.r0 + r_geo = eqt.boundary.geometric_axis.r + z_geo = eqt.boundary.geometric_axis.z + Bt_geo = Bt_center * r_center / r_geo + ϵ = eqt.boundary.minor_radius / r_geo + + # pressure and j_tor + psin = eq1d.psi_norm + j_tor = eq1d.j_tor + pressure = eq1d.pressure + rho_pol = sqrt.(psin) + pressure_sep = pressure[end] + + # run and handle errors + try + actor.chease = CHEASE.run_chease( + ϵ, z_geo, pressure_sep, Bt_geo, + r_geo, Ip, r_bound, z_bound, 82, + rho_pol, pressure, j_tor, + rescale_eq_to_ip=actor.par.rescale_eq_to_ip, + clear_workdir=actor.par.clear_workdir) + catch + display(plot(r_bound, z_bound; marker=:dot, aspect_ratio=:equal)) + display(plot(psin, pressure)) + display(plot(psin, abs.(j_tor))) + rethrow() + end + + # convert from fixed to free boundary equilibrium + if actor.par.free_boundary + EQ = Equilibrium.efit(actor.chease.gfile, 1) + psi_free_rz = VacuumFields.fixed2free(EQ, actor.chease.gfile.nbbbs) + actor.chease.gfile.psirz = psi_free_rz + EQ = Equilibrium.efit(actor.chease.gfile, 1) + psi_b = Equilibrium.psi_boundary(EQ; r=EQ.r, z=EQ.z) + psi_a = EQ.psi_rz(EQ.axis...) + actor.chease.gfile.psirz = (psi_free_rz .- psi_a) * ((actor.chease.gfile.psi[end] - actor.chease.gfile.psi[1]) / (psi_b - psi_a)) .+ actor.chease.gfile.psi[1] + end + + return actor +end + +# define `finalize` function for this actor +function finalize(actor::ActorCHEASE) + gEQDSK2IMAS(actor.chease.gfile, actor.dd.equilibrium) + return actor +end + +""" + gEQDSK2IMAS(GEQDSKFile::GEQDSKFile,eq::IMAS.equilibrium) + +Convert IMAS.equilibrium__time_slice to Equilibrium.jl EFIT structure +""" +function gEQDSK2IMAS(g::EFIT.GEQDSKFile, eq::IMAS.equilibrium) + + tc = transform_cocos(1, 11) # chease output is cocos 1 , dd is cocos 11 + + eqt = eq.time_slice[] + eq1d = eqt.profiles_1d + resize!(eqt.profiles_2d, 1) + eq2d = eqt.profiles_2d[1] + + @ddtime(eq.vacuum_toroidal_field.b0 = g.bcentr) + eq.vacuum_toroidal_field.r0 = g.rcentr + + eqt.global_quantities.magnetic_axis.r = g.rmaxis + eqt.boundary.geometric_axis.r = g.rcentr + eqt.boundary.geometric_axis.z = g.zmid + 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"] + + IMAS.flux_surfaces(eqt) +end \ No newline at end of file diff --git a/src/actors/compound_actors.jl b/src/actors/compound_actors.jl index e8a7f3acc..97b5a156c 100644 --- a/src/actors/compound_actors.jl +++ b/src/actors/compound_actors.jl @@ -1,11 +1,13 @@ #= ========================= =# # ActorEquilibriumTransport # #= ========================= =# -Base.@kwdef mutable struct ActorEquilibriumTransport <: PlasmaAbstractActor +mutable struct ActorEquilibriumTransport <: PlasmaAbstractActor dd::IMAS.dd - actor_jt::Union{Nothing,ActorSteadyStateCurrent} - actor_eq::Union{Nothing,ActorEquilibrium} - actor_tr::Union{Nothing,ActorTauenn} + par::ParametersActor + act::ParametersAllActors + actor_jt::ActorSteadyStateCurrent + actor_eq::ActorEquilibrium + actor_tr::ActorTauenn end function ParametersActor(::Type{Val{:ActorEquilibriumTransport}}) @@ -19,53 +21,58 @@ end ActorEquilibriumTransport(dd::IMAS.dd, act::ParametersAllActors; kw...) Compound actor that runs the following actors in succesion: -```julia -ActorSteadyStateCurrent(dd, act) # Current evolution to steady-state -ActorTauenn(dd, act) # For transport -ActorEquilibrium(dd, act) # Equilibrium -``` +* ActorSteadyStateCurrent +* ActorTauenn +* ActorEquilibrium !!! note Stores data in `dd.equilibrium, dd.core_profiles, dd.core_sources` """ function ActorEquilibriumTransport(dd::IMAS.dd, act::ParametersAllActors; kw...) par = act.ActorEquilibriumTransport(kw...) - actor = ActorEquilibriumTransport(dd, nothing, nothing, nothing) - step(actor; act, iterations=par.iterations, do_plot=par.do_plot) + actor = ActorEquilibriumTransport(dd, par, act) + step(actor) finalize(actor) return actor end -function step(actor::ActorEquilibriumTransport; act::Union{Missing,ParametersAllActors}=missing, iterations::Int=1, do_plot::Bool=false) +function ActorEquilibriumTransport(dd::IMAS.dd, par::ParametersActor, act::ParametersAllActors; kw...) + par = par(kw...) + actor_jt = ActorSteadyStateCurrent(dd, act.ActorSteadyStateCurrent) + actor_eq = ActorEquilibrium(dd, act.ActorEquilibrium, act) + actor_tr = ActorTauenn(dd, act.ActorTauenn) + return ActorEquilibriumTransport(dd, par, act, actor_jt, actor_eq, actor_tr) +end + +function step(actor::ActorEquilibriumTransport) dd = actor.dd - if act === missing - act = ParametersAllActors() - end + par = actor.par + act = actor.act - if do_plot + if par.do_plot pe = plot(dd.equilibrium; color=:gray, label="old", coordinate=:rho_tor_norm) pp = plot(dd.core_profiles; color=:gray, label=" (old)") ps = plot(dd.core_sources; color=:gray, label=" (old)") end # Set j_ohmic to steady state - ActorSteadyStateCurrent(dd, act) + finalize(step(actor.actor_jt)) - for iteration in 1:iterations + for iteration in 1:par.iterations # run transport actor - actor.actor_tr = ActorTauenn(dd, act) + finalize(step(actor.actor_tr)) # prepare equilibrium input based on transport core_profiles output prepare(dd, :ActorEquilibrium, act) # run equilibrium actor with the updated beta - actor.actor_eq = ActorEquilibrium(dd, act) + finalize(step(actor.actor_eq)) # Set j_ohmic to steady state - actor.actor_jt = ActorSteadyStateCurrent(dd, act) + finalize(step(actor.actor_jt)) end - if do_plot + if par.do_plot display(plot!(pe, dd.equilibrium, coordinate=:rho_tor_norm)) display(plot!(pp, dd.core_profiles)) display(plot!(ps, dd.core_sources)) @@ -78,8 +85,20 @@ end #= ================== =# # ActorWholeFacility # #= ================== =# -Base.@kwdef mutable struct ActorWholeFacility <: FacilityAbstractActor +mutable struct ActorWholeFacility <: FacilityAbstractActor dd::IMAS.dd + par::ParametersActor + act::ParametersAllActors + EquilibriumTransport::Union{Nothing,ActorEquilibriumTransport} + HFSsizing::Union{Nothing,ActorHFSsizing} + LFSsizing::Union{Nothing,ActorLFSsizing} + CXbuild::Union{Nothing,ActorCXbuild} + PFcoilsOpt::Union{Nothing,ActorPFcoilsOpt} + Neutronics::Union{Nothing,ActorNeutronics} + Blanket::Union{Nothing,ActorBlanket} + Divertors::Union{Nothing,ActorDivertors} + BalanceOfPlant::Union{Nothing,ActorBalanceOfPlant} + Costing::Union{Nothing,ActorCosting} end function ParametersActor(::Type{Val{:ActorWholeFacility}}) @@ -90,30 +109,57 @@ end """ ActorWholeFacility(dd::IMAS.dd, act::ParametersAllActors; kw...) -Compound actor that runs all the actors needed to model the whole plant +Compound actor that runs all the actors needed to model the whole plant: +* ActorEquilibriumTransport +* ActorHFSsizing +* ActorLFSsizing +* ActorCXbuild +* ActorPFcoilsOpt +* ActorNeutronics +* ActorBlanket +* ActorDivertors +* ActorBalanceOfPlant +* ActorCosting !!! note Stores data in `dd` """ function ActorWholeFacility(dd::IMAS.dd, act::ParametersAllActors; kw...) par = act.ActorWholeFacility(kw...) - actor = ActorWholeFacility(dd) - step(actor; act) + actor = ActorWholeFacility(dd, par, act) + step(actor) finalize(actor) return actor end -function step(actor::ActorWholeFacility; act::Union{Missing,ParametersAllActors}=missing, iterations::Int=1, do_plot::Bool=false) +function ActorWholeFacility(dd::IMAS.dd, par::ParametersActor, act::ParametersAllActors; kw...) + + par = par(kw...) + ActorWholeFacility(dd, par, act, + nothing, + nothing, + nothing, + nothing, + nothing, + nothing, + nothing, + nothing, + nothing, + nothing) +end + +function step(actor::ActorWholeFacility) dd = actor.dd - ActorEquilibriumTransport(dd, act) - ActorHFSsizing(dd, act) - ActorLFSsizing(dd, act) - ActorCXbuild(dd, act) - ActorPFcoilsOpt(dd, act) - ActorNeutronics(dd, act) - ActorBlanket(dd, act) - ActorDivertors(dd, act) - ActorBalanceOfPlant(dd,act) - ActorCosting(dd, act) + act = actor.act + actor.EquilibriumTransport = ActorEquilibriumTransport(dd, act) + actor.HFSsizing = ActorHFSsizing(dd, act) + actor.LFSsizing = ActorLFSsizing(dd, act) + actor.CXbuild = ActorCXbuild(dd, act) + actor.PFcoilsOpt = ActorPFcoilsOpt(dd, act) + actor.Neutronics = ActorNeutronics(dd, act) + actor.Blanket = ActorBlanket(dd, act) + actor.Divertors = ActorDivertors(dd, act) + actor.BalanceOfPlant = ActorBalanceOfPlant(dd, act) + actor.Costing = ActorCosting(dd, act) return actor end \ No newline at end of file diff --git a/src/actors/costing_actors.jl b/src/actors/costing_actors.jl index 19d4e5720..b166ce75e 100644 --- a/src/actors/costing_actors.jl +++ b/src/actors/costing_actors.jl @@ -1,9 +1,9 @@ #= ============== =# -# materials cost # +# materials cost # #= ============== =# #NOTE: material should be priced by Kg #NOTE: if something is priced by m^3 then it is for a specific part already -function unit_cost(material::String) +function unit_cost(material::AbstractString) if material == "Vacuum" return 0.0 # $M/m^3 elseif material == "ReBCO" @@ -29,13 +29,49 @@ function unit_cost(material::String) end end -function cost(layer::IMAS.build__layer) +function unit_cost(coil_tech::Union{IMAS.build__tf__technology,IMAS.build__oh__technology,IMAS.build__pf_active__technology}) + if coil_tech.material == "Copper" + return unit_cost("Copper") + else + fraction_cable = 1 - coil_tech.fraction_stainless - coil_tech.fraction_void + fraction_SC = fraction_cable * coil_tech.ratio_SC_to_copper + fraction_copper = fraction_cable - fraction_SC + return (coil_tech.fraction_stainless * unit_cost("Steel, Stainless 316") + fraction_copper * unit_cost("Copper") + fraction_SC * unit_cost(coil_tech.material)) + end +end + +#= ================== =# +# Dispatch on symbol # +#= ================== =# + +function cost_direct_capital(item::Symbol, args...; kw...) + return cost_direct_capital(Val{item}, args...; kw...) +end + +function cost_operations(item::Symbol, args...; kw...) + return cost_operations(Val{item}, args...; kw...) +end + +function cost_decomissioning(item::Symbol, args...; kw...) + return cost_decomissioning(Val{item}, args...; kw...) +end + +#= =================== =# +# direct capital cost # +#= =================== =# + +""" + cost_direct_capital(layer::IMAS.build__layer) + +Capital cost for each layer in the build +""" +function cost_direct_capital(layer::IMAS.build__layer) if layer.type == Int(_oh_) build = IMAS.parent(IMAS.parent(layer)) - return unit_cost(build.oh.technology) * layer.volume + return layer.volume * unit_cost(build.oh.technology) elseif layer.type == Int(_tf_) build = IMAS.parent(IMAS.parent(layer)) - return unit_cost(build.tf.technology) * layer.volume + return layer.volume * unit_cost(build.tf.technology) elseif layer.type == Int(_shield_) return layer.volume * 0.29 # $M/m^3 elseif layer.type == Int(_blanket_) @@ -43,64 +79,248 @@ function cost(layer::IMAS.build__layer) elseif layer.type ∈ [Int(_wall_), Int(_vessel_), Int(_cryostat_)] return layer.volume * 0.36 # $M/m^3 else - return unit_cost(layer.material) * layer.volume + return layer.volume * unit_cost(layer.material) end end -function cost(ecl::IMAS.ec_launchers__launcher) - ecl.available_launch_power / 1E6 * 3.0 # $/W #ARIES -end +""" + cost_direct_capital(ecl::IMAS.ec_launchers__beam) -function cost(ica::IMAS.ic_antennas__antenna) - ica.available_launch_power / 1E6 * 1.64 #$/W ARIES +Capital cost for each EC launcer is proportional to its power + +NOTE: ARIES https://cer.ucsd.edu/_files/publications/UCSD-CER-13-01.pdf +""" +function cost_direct_capital(ecl::IMAS.ec_launchers__beam) + unit_cost = 3.0 # $/W + return ecl.available_launch_power / 1E6 * unit_cost end -function cost(lha::IMAS.lh_antennas__antenna) - lha.available_launch_power / 1E6 * 2.13 #$/W ARIES +""" + cost_direct_capital(ica::IMAS.ic_antennas__antenna) + +Capital cost for each IC antenna is proportional to its power + +NOTE: ARIES https://cer.ucsd.edu/_files/publications/UCSD-CER-13-01.pdf +""" +function cost_direct_capital(ica::IMAS.ic_antennas__antenna) + unit_cost = 1.64 # $/W + return ica.available_launch_power / 1E6 * unit_cost end -function cost(nbu::IMAS.nbi__unit) - nbu.available_launch_power / 1E6 * 4.93 #$/W ARIES +""" + cost_direct_capital(lha::IMAS.lh_antennas__antenna) + +Capital cost for each LH antenna is proportional to its power + +NOTE: ARIES https://cer.ucsd.edu/_files/publications/UCSD-CER-13-01.pdf +""" +function cost_direct_capital(lha::IMAS.lh_antennas__antenna) + unit_cost = 2.13 # $/W + return lha.available_launch_power / 1E6 * unit_cost end -function unit_cost(coil_tech::Union{IMAS.build__tf__technology,IMAS.build__oh__technology,IMAS.build__pf_active__technology}) - if coil_tech.material == "Copper" - return unit_cost("Copper") - else - fraction_cable = 1 - coil_tech.fraction_stainless - coil_tech.fraction_void - fraction_SC = fraction_cable * coil_tech.ratio_SC_to_copper - fraction_copper = fraction_cable - fraction_SC - return (coil_tech.fraction_stainless * unit_cost("Steel, Stainless 316") + fraction_copper * unit_cost("Copper") + fraction_SC * unit_cost(coil_tech.material)) - end +""" + cost_direct_capital(nbu::IMAS.nbi__unit) + +Capital cost for each NBI unit is proportional to its power + +NOTE: ARIES https://cer.ucsd.edu/_files/publications/UCSD-CER-13-01.pdf +""" +function cost_direct_capital(nbu::IMAS.nbi__unit) + unit_cost = 4.93 # $/W + return nbu.available_launch_power / 1E6 * unit_cost end -function cost(pf_active::IMAS.pf_active) +""" + cost_direct_capital(pf_active::IMAS.pf_active) +""" +function cost_direct_capital(pf_active::IMAS.pf_active) dd = IMAS.top_dd(pf_active) c = Dict("OH" => 0.0, "PF" => 0.0) for coil in pf_active.coil if coil.name == "OH" - c["OH"] += cost(coil, dd.build.oh.technology) + c["OH"] += cost_direct_capital(coil, dd.build.oh.technology) else - c["PF"] += cost(coil, dd.build.pf_active.technology) + c["PF"] += cost_direct_capital(coil, dd.build.pf_active.technology) end end return c end -function cost(coil::IMAS.pf_active__coil, technology::Union{IMAS.build__tf__technology,IMAS.build__oh__technology,IMAS.build__pf_active__technology}) +""" + cost_direct_capital(coil::IMAS.pf_active__coil, technology::Union{IMAS.build__tf__technology,IMAS.build__oh__technology,IMAS.build__pf_active__technology}) + +""" +function cost_direct_capital(coil::IMAS.pf_active__coil, technology::Union{IMAS.build__tf__technology,IMAS.build__oh__technology,IMAS.build__pf_active__technology}) return IMAS.volume(coil) * unit_cost(technology) end +""" + cost_direct_capital(::Type{Val{:land}}, land::Real, power_electric_net::Real) + +NOTE: ARIES https://cer.ucsd.edu/_files/publications/UCSD-CER-13-01.pdf +""" +function cost_direct_capital(::Type{Val{:land}}, land::Real, power_electric_net::Real) + power_electric_net = power_electric_net / 1E6 + return 1.2 * land * 20.0e-3 * (power_electric_net / 1000.0)^0.3 +end + +""" + cost_direct_capital(::Type{Val{:buildings}}, land::Real, building_volume::Real, power_electric_generated::Real, power_thermal::Real, power_electric_net::Real) + +NOTE: ARIES https://cer.ucsd.edu/_files/publications/UCSD-CER-13-01.pdf +""" +function cost_direct_capital(::Type{Val{:buildings}}, land::Real, building_volume::Real, power_electric_generated::Real, power_thermal::Real, power_electric_net::Real) + power_electric_generated = power_electric_generated / 1E6 + power_thermal = power_thermal / 1E6 + power_electric_net = power_electric_net / 1E6 + cost = 27.0 * (land / 1000.0)^0.2 # site + cost += 111.661 * (building_volume / 80.0e3)^0.62 # tokamak building + cost += 78.9 * (power_electric_generated / 1246)^0.5 # turbine/generator + cost += 16.804 * ((power_thermal - power_electric_net) / 1860.0)^0.5 # heat rejection + cost += 22.878 * (power_electric_net / 1000.0)^0.3 # electrical equipment + cost += 21.96 * (power_electric_generated / 1246)^0.3 # Plant auxiliary systems + cost += 4.309 * (power_electric_generated / 1000.0)^0.3 # power core service building + cost += 1.513 * (power_electric_generated / 1000.0)^0.3 # service water + cost += 25.0 * (power_thermal / 1759.0)^0.3 # tritium plant and systems + cost += 7.11 # control room + cost += 4.7 * (power_electric_net / 1000.0)^0.3 # On-site AC power supply building + cost += 2.0 # administrative + cost += 2.0 # site service + cost += 2.09 # cyrogenic and inert gas storage + cost += 0.71 # security + cost += 4.15 # ventilation stack + return cost +end + +""" + cost_direct_capital(::Type{Val{:hot_cell}}, building_volume) + +NOTE: https://www.iter.org/mach/HotCell +""" +function cost_direct_capital(::Type{Val{:hot_cell}}, building_volume::Real) + return 0.34 * 111.661 * (building_volume / 80.0e3)^0.62 +end + +""" + cost_direct_capital(::Type{Val{:heat_transfer_loop_materials}}, power_thermal::Real) + +NOTE: ARIES https://cer.ucsd.edu/_files/publications/UCSD-CER-13-01.pdf (warning uses LiPb for blanket) +""" +function cost_direct_capital(::Type{Val{:heat_transfer_loop_materials}}, power_thermal::Real) + power_thermal = power_thermal / 1E6 + cost = 50.0 * (power_thermal / 2000.0)^0.55 # water + cost += 125.0 * (power_thermal / 2000.0)^0.55 # LiPb + cost += 110.0 * (power_thermal / 2000.0)^0.55 # He + cost += 0.01 * power_thermal # NbIHX + cost += 50.0 * (power_thermal / 2000.0)^0.55 # Na + return cost +end + +""" + cost_direct_capital(::Type{Val{:balance_of_plant_equipment}}, power_thermal::Real, power_electric_generated::Real) + +NOTE: ARIES https://cer.ucsd.edu/_files/publications/UCSD-CER-13-01.pdf +""" +function cost_direct_capital(::Type{Val{:balance_of_plant_equipment}}, power_thermal::Real, power_electric_generated::Real) + power_thermal = power_thermal / 1E6 + power_electric_generated = power_electric_generated / 1E6 + cost = 350.0 * (power_thermal / 2620.0)^0.7 # Turbine equipment + cost += 182.98 * (power_electric_generated / 1200.0)^0.5 # Electrical plant equipment + cost += 87.52 * ((power_thermal - power_electric_generated) / 2300.0) # Heat rejection equipment + cost += 88.89 * (power_electric_generated / 1200.0)^0.6 # Miscellenous equipment + return cost +end + +""" + cost_direct_capital(::Type{Val{:fuel_cycle_rad_handling}}, power_thermal::Real, power_electric_net::Real) + +NOTE: ARIES https://cer.ucsd.edu/_files/publications/UCSD-CER-13-01.pdf (warning uses LiPb for blanket) +""" +function cost_direct_capital(::Type{Val{:fuel_cycle_rad_handling}}, power_thermal::Real, power_electric_net::Real) + power_thermal = power_thermal / 1E6 + power_electric_net = power_electric_net / 1E6 + cost = 15.0 * (power_thermal / 1758.0)^0.85 # radioactive material treatment and management + cost += 70.0 * (power_thermal / 1758.0)^0.8 # Fuel handling and storage + cost += 100.0 * (power_electric_net / 2000.0)^0.55 # Hot cell maintanance + cost += 60.0 # Instrumentation and Control + cost += 8.0 * (power_thermal / 1000.0)^0.8 # Misc power core equipment + return cost +end + +#= ====================== =# +# yearly operations cost # +#= ====================== =# + +""" + cost_operations(::Type{Val{:operation_maintanance}}, power_electric_generated::Real) + +Yearly cost for maintenance [\$M/year] +NOTE: ARIES https://cer.ucsd.edu/_files/publications/UCSD-CER-13-01.pdf +""" +function cost_operations(::Type{Val{:operation_maintanance}}, power_electric_generated::Real) + power_electric_generated = power_electric_generated / 1E6 + return 80.0 * (power_electric_generated / 1200.0)^0.5 +end + +""" + cost_operations(::Type{Val{:tritium_handling}}) + +Yearly cost for tritium_handling [\$M/year] +!!!!WRONG!!!! Needs estiamte +""" +function cost_operations(::Type{Val{:tritium_handling}}) + return 1.0 +end + +""" + cost_operations(::Type{Val{:blanket_replacement}}, cost_blanket::Real, blanket_lifetime::Real) + +Yearly cost for blanket replacement [\$M/year] +""" +function cost_operations(::Type{Val{:blanket_replacement}}, cost_blanket::Real, blanket_lifetime::Real) + return cost_blanket / blanket_lifetime +end + +#= =================== =# +# Decomissioning cost # +#= =================== =# + +""" + cost_decomissioning(::Type{Val{:decom_wild_guess}}, lifetime::Real) + +Cost to decommission the plant [\$M] + +LIKELY NEEDS FIXING +""" +function cost_decomissioning(::Type{Val{:decom_wild_guess}}, lifetime::Real) + unit_cost = 2.76 # [$M/year]From GASC + return unit_cost * lifetime +end + #= ============ =# # ActorCosting # #= ============ =# -Base.@kwdef mutable struct ActorCosting <: FacilityAbstractActor +mutable struct ActorCosting <: FacilityAbstractActor dd::IMAS.dd + par::ParametersActor + function ActorCosting(dd::IMAS.dd, par::ParametersActor; kw...) + par = par(kw...) + return new(dd, par) + end end function ParametersActor(::Type{Val{:ActorCosting}}) par = ParametersActor(nothing) + par.land_space = Entry(Real, "acres", "Plant site space required in acres"; default=1000.0) + par.building_volume = Entry(Real, "m^3", "Volume of the tokmak building"; default=140.0e3) + par.interest_rate = Entry(Real, "", "Anual interest rate fraction of direct capital cost"; default=0.05) + par.indirect_cost_rate = Entry(Real, "", "Indirect cost associated with construction, equipment, services, energineering construction management and owners cost"; default=0.4) + par.lifetime = Entry(Integer, "years", "lifetime of the plant"; default=40) + par.availability = Entry(Real, "", "availability fraction of the plant"; default=0.803) + par.escalation_fraction = Entry(Real, "", "yearly escalation fraction based on risk assessment"; default=0.05) + par.blanket_lifetime = Entry(Real, "years", "lifetime of the blanket"; default=6.8) return par end @@ -114,53 +334,134 @@ This actor estimates the cost of the fusion power plant. """ function ActorCosting(dd::IMAS.dd, act::ParametersAllActors; kw...) par = act.ActorCosting(kw...) - actor = ActorCosting(dd) + actor = ActorCosting(dd, par) step(actor) finalize(actor) return actor end function step(actor::ActorCosting) + par = actor.par dd = actor.dd cst = dd.costing - empty!(cst) + cost_direct = cst.cost_direct_capital + cost_ops = cst.cost_operations + cost_decom = cst.cost_decommissioning + + ###### Direct Capital ###### + + empty!(cost_direct) - # TOKAMAK - sys = resize!(cst.system, "name" => "tokamak") + ### Tokamak + + # build layers + sys = resize!(cost_direct.system, "name" => "tokamak") for layer in dd.build.layer if layer.fs == Int(_lfs_) continue # avoid double counting of hfs and lfs layers elseif layer.type == Int(_oh_) continue # avoid double counting of oh end - c = cost(layer) + c = cost_direct_capital(layer) if c > 0 sub = resize!(sys.subsystem, "name" => replace(layer.name, r"^hfs " => "")) sub.cost = c end end - for (name, c) in cost(dd.pf_active) + + # PF coils + for (name, c) in cost_direct_capital(dd.pf_active) sub = resize!(sys.subsystem, "name" => name) sub.cost = c end - # HCD - sys = resize!(cst.system, "name" => "hcd") - for hcd in vcat(dd.ec_launchers.launcher, dd.ic_antennas.antenna, dd.lh_antennas.antenna, dd.nbi.unit) - c = cost(hcd) + # Heating and current drive + for hcd in vcat(dd.ec_launchers.beam, dd.ic_antennas.antenna, dd.lh_antennas.antenna, dd.nbi.unit) + c = cost_direct_capital(hcd) if c > 0 sub = resize!(sys.subsystem, "name" => hcd.name) sub.cost = c end end + ### Facility + sys = resize!(cost_direct.system, "name" => "Facility structures, buildings and site") + + if ismissing(dd.balance_of_plant.thermal_cycle, :power_electric_generated) || @ddtime(dd.balance_of_plant.thermal_cycle.power_electric_generated) < 0 + @warn("The plant doesn't generate net electricity therefore costing excludes facility estimates") + power_electric_net = 0.0 + power_thermal = 0.0 + power_electric_generated = 0.0 + else + power_electric_net = @ddtime(dd.balance_of_plant.power_electric_net) # should be pulse average + power_thermal = @ddtime(dd.balance_of_plant.thermal_cycle.power_thermal_convertable_total) + power_electric_generated = @ddtime(dd.balance_of_plant.thermal_cycle.power_electric_generated) + + for item in vcat(:land, :buildings, :hot_cell, :heat_transfer_loop_materials, :balance_of_plant_equipment, :fuel_cycle_rad_handling) + sub = resize!(sys.subsystem, "name" => string(item)) + if item == :land + sub.cost = cost_direct_capital(item, par.land_space, power_electric_generated) + elseif item == :buildings + sub.cost = cost_direct_capital(item, par.building_volume, + par.land_space, power_electric_generated, + power_thermal, power_electric_net) + elseif item == :hot_cell + sub.cost = cost_direct_capital(item, par.building_volume) + elseif item == :heat_transfer_loop_materials + sub.cost = cost_direct_capital(item, power_thermal) + elseif item == :balance_of_plant_equipment + sub.cost = cost_direct_capital(item, power_thermal, power_electric_generated) + elseif item == :fuel_cycle_rad_handling + sub.cost = cost_direct_capital(item, power_thermal, power_electric_net) + else + sub.cost = cost_direct_capital(item) + end + end + end + + ###### Operations ###### + empty!(cost_ops) + + sys = resize!(cost_ops.system, "name" => "tritium handling") + sys.yearly_cost = cost_operations(:tritium_handling) + + sys = resize!(cost_ops.system, "name" => "maintanance and operators") + sys.yearly_cost = cost_operations(:operation_maintanance, power_electric_generated) + + sys = resize!(cost_ops.system, "name" => "replacements") + for item in [:blanket_replacement] + sub = resize!(sys.subsystem, "name" => string(item)) + if item == :blanket_replacement + tokamak = cost_direct.system[findfirst(system -> system.name == "tokamak", cost_direct.system)] + blanket_cost = sum([item.cost for item in tokamak.subsystem if item.name == "blanket"]) + sub.yearly_cost = cost_operations(:blanket_replacement, blanket_cost, par.blanket_lifetime) + else + sub.yearly_cost = cost_operations(item) + end + end + + ###### Decomissioning ###### + empty!(cost_decom) + + sys = resize!(cost_decom.system, "name" => "decommissioning") + sys.cost = cost_decomissioning(:decom_wild_guess, par.lifetime) + + ###### Levelized Cost Of Electricity ###### + capital_cost_rate = par.interest_rate / (1 - (1 + par.interest_rate)^(-1.0 * par.lifetime)) + lifetime_cost = 0.0 + for year in 1:par.lifetime + yearly_cost = (capital_cost_rate * cost_direct.cost + cost_ops.yearly_cost + cost_decom.cost / par.lifetime) + lifetime_cost += (1.0 + par.escalation_fraction) * (1.0 + par.indirect_cost_rate) * yearly_cost + end + dd.costing.cost_lifetime = lifetime_cost + dd.costing.levelized_CoE = (dd.costing.cost_lifetime * 1E6) / (par.lifetime * 24 * 365 * power_electric_net / 1e3 * par.availability) return actor end function finalize(actor::ActorCosting) - # sort system/subsystem costs - sort!(actor.dd.costing.system, by=x -> x.cost, rev=true) - for sys in actor.dd.costing.system + # sort system/subsystem by their costs + sort!(actor.dd.costing.cost_direct_capital.system, by=x -> x.cost, rev=true) + for sys in actor.dd.costing.cost_direct_capital.system sort!(sys.subsystem, by=x -> x.cost, rev=true) end -end \ No newline at end of file +end diff --git a/src/actors/current_actors.jl b/src/actors/current_actors.jl index 26ece8394..7663faaf1 100644 --- a/src/actors/current_actors.jl +++ b/src/actors/current_actors.jl @@ -3,8 +3,9 @@ import QED #= =============== =# # ActorQEDcurrent # #= =============== =# -Base.@kwdef mutable struct ActorQEDcurrent <: PlasmaAbstractActor +mutable struct ActorQEDcurrent <: PlasmaAbstractActor dd::IMAS.dd + par::ParametersActor QI::QED.QED_state η#::Base.Callable QO::Union{QED.QED_state,Missing} @@ -27,14 +28,15 @@ This actor evolves the current using QED. """ function ActorQEDcurrent(dd::IMAS.dd, act::ParametersAllActors; kw...) par = act.ActorQEDcurrent(kw...) - actor = ActorQEDcurrent(dd) + actor = ActorQEDcurrent(dd, par) step(actor) finalize(actor) return actor end -function ActorQEDcurrent(dd::IMAS.dd) - ActorQEDcurrent(dd, from_imas(dd), η_imas(dd), missing, @ddtime(dd.equilibrium.time), 0.0) +function ActorQEDcurrent(dd::IMAS.dd, par::ParametersActor; kw...) + par = par(kw...) + return ActorQEDcurrent(dd, par, from_imas(dd), η_imas(dd), missing, @ddtime(dd.equilibrium.time), dd.global_time) end function step(actor::ActorQEDcurrent, tmax, Nt, Vedge=nothing, Ip=nothing; resume=false) @@ -45,7 +47,6 @@ function step(actor::ActorQEDcurrent, tmax, Nt, Vedge=nothing, Ip=nothing; resum end actor.tmax += tmax else - actor.initial_time = @ddtime(dd.equilibrium.time) actor.tmax = tmax end actor.QO = QED.diffuse(actor.QI, actor.η, tmax, Nt; Vedge, Ip) @@ -113,8 +114,13 @@ end #= ======================= =# # ActorSteadyStateCurrent # #= ======================= =# -Base.@kwdef mutable struct ActorSteadyStateCurrent <: PlasmaAbstractActor +mutable struct ActorSteadyStateCurrent <: PlasmaAbstractActor dd::IMAS.dd + par::ParametersActor + function ActorSteadyStateCurrent(dd::IMAS.dd, par::ParametersActor; kw...) + par = par(kw...) + return new(dd, par) + end end function ParametersActor(::Type{Val{:ActorSteadyStateCurrent}}) @@ -134,7 +140,7 @@ Also sets the ohmic, bootstrap and non-inductive current profiles in `dd.core_pr """ function ActorSteadyStateCurrent(dd::IMAS.dd, act::ParametersAllActors; kw...) par = act.ActorSteadyStateCurrent(kw...) - actor = ActorSteadyStateCurrent(dd, par...) + actor = ActorSteadyStateCurrent(dd, par) step(actor) finalize(actor) return actor @@ -146,4 +152,5 @@ function step(actor::ActorSteadyStateCurrent) # update core_sources related to current IMAS.bootstrap_source!(dd) IMAS.ohmic_source!(dd) + return actor end diff --git a/src/actors/divertors_actors.jl b/src/actors/divertors_actors.jl index 8e6200b52..d0e1634de 100644 --- a/src/actors/divertors_actors.jl +++ b/src/actors/divertors_actors.jl @@ -2,8 +2,9 @@ # ActorDivertors # #= ============== =# -Base.@kwdef mutable struct ActorDivertors <: ReactorAbstractActor +mutable struct ActorDivertors <: ReactorAbstractActor dd::IMAS.dd + par::ParametersActor thermal_power_extraction_efficiency::Real end @@ -24,12 +25,17 @@ Simple Divertor actor """ function ActorDivertors(dd::IMAS.dd, act::ParametersAllActors; kw...) par = act.ActorDivertors(kw...) - actor = ActorDivertors(dd, par.thermal_power_extraction_efficiency) + actor = ActorDivertors(dd, par) step(actor) finalize(actor) return actor end +function ActorDivertors(dd::IMAS.dd, par::ParametersActor; kw...) + par = par(kw...) + return ActorDivertors(dd, par, par.thermal_power_extraction_efficiency) +end + function step(actor::ActorDivertors) dd = actor.dd diff --git a/src/actors/equilibrium_actors.jl b/src/actors/equilibrium_actors.jl index 90c0287fa..d851631a7 100644 --- a/src/actors/equilibrium_actors.jl +++ b/src/actors/equilibrium_actors.jl @@ -1,13 +1,7 @@ -import Equilibrium -import EFIT -import CHEASE -import ForwardDiff -import Optim - #= ================ =# # ActorEquilibrium # #= ================ =# -Base.@kwdef mutable struct ActorEquilibrium <: PlasmaAbstractActor +mutable struct ActorEquilibrium <: PlasmaAbstractActor dd::IMAS.dd par::ParametersActor eq_actor::PlasmaAbstractActor @@ -32,7 +26,8 @@ function ActorEquilibrium(dd::IMAS.dd, act::ParametersAllActors; kw...) return actor end -function ActorEquilibrium(dd::IMAS.dd, par::ParametersActor, act::ParametersAllActors) +function ActorEquilibrium(dd::IMAS.dd, par::ParametersActor, act::ParametersAllActors; kw...) + par = par(kw...) if par.model == :Solovev eq_actor = ActorSolovev(dd, act.ActorSolovev) elseif par.model == :CHEASE @@ -40,7 +35,7 @@ function ActorEquilibrium(dd::IMAS.dd, par::ParametersActor, act::ParametersAllA else error("ActorEquilibrium: model = $(par.model) is unknown") end - return ActorEquilibrium(dd, deepcopy(par), eq_actor) + return ActorEquilibrium(dd, par, eq_actor) end """ @@ -81,343 +76,9 @@ end #= ============ =# # ActorSolovev # #= ============ =# -Base.@kwdef mutable struct ActorSolovev <: PlasmaAbstractActor - eq::IMAS.equilibrium - par::ParametersActor - S::Equilibrium.SolovevEquilibrium -end - -function ParametersActor(::Type{Val{:ActorSolovev}}) - par = ParametersActor(nothing) - par.ngrid = Entry(Integer, "", "Grid size (for R, Z follows proportionally to plasma elongation)"; default=129) - par.qstar = Entry(Real, "", "Initial guess of kink safety factor"; default=1.5) - par.alpha = Entry(Real, "", "Initial guess of constant relating to beta regime"; default=0.0) - par.volume = Entry(Real, "m³", "Scalar volume to match (optional)"; default=missing) - par.area = Entry(Real, "m²", "Scalar area to match (optional)"; default=missing) - par.verbose = Entry(Bool, "", "verbose"; default=false) - return par -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...) - par = act.ActorSolovev(kw...) - actor = ActorSolovev(dd, par) - step(actor) - finalize(actor) - # record optimized values of qstar and alpha in `act` for subsequent ActorSolovev calls - par.qstar = actor.S.qstar - par.alpha = actor.S.alpha - return actor -end - -function ActorSolovev(dd::IMAS.dd, par::ParametersActor) - # extract info from dd - eq = dd.equilibrium - eqt = eq.time_slice[] - a = eqt.boundary.minor_radius - R0 = eqt.boundary.geometric_axis.r - κ = eqt.boundary.elongation - δ = eqt.boundary.triangularity - ϵ = a / R0 - B0 = @ddtime eq.vacuum_toroidal_field.b0 - - # 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)) - else - x_point = nothing - end - - # run Solovev - S = Equilibrium.solovev(abs(B0), R0, ϵ, δ, κ, par.alpha, par.qstar, B0_dir=Int64(sign(B0)), Ip_dir=1, x_point=x_point, symmetric=symmetric) - - return ActorSolovev(dd.equilibrium, deepcopy(par), S) -end - -""" - prepare(dd::IMAS.dd, :ActorSolovev, act::ParametersAllActors; kw...) - -Prepare dd to run ActorSolovev -* Copy beta_n from summary to equilibrium -""" -function prepare(dd::IMAS.dd, ::Type{Val{:ActorSolovev}}, act::ParametersAllActors; kw...) - dd.equilibrium.time_slice[].global_quantities.beta_normal = @ddtime(dd.summary.global_quantities.beta_tor_thermal_norm.value) -end - -""" - step(actor::ActorSolovev) - -Non-linear optimization to obtain a target `ip` and `beta_normal` -""" -function step(actor::ActorSolovev) - S0 = actor.S - par = actor.par - - eqt = actor.eq.time_slice[] - target_ip = abs(eqt.global_quantities.ip) - target_beta = eqt.global_quantities.beta_normal - - B0, R0, epsilon, delta, kappa, alpha, qstar, target_ip, target_beta = promote(S0.B0, S0.R0, S0.epsilon, S0.delta, S0.kappa, S0.alpha, S0.qstar, target_ip, target_beta) - - function cost(x) - # NOTE: Ip/Beta calculation is very much off in Equilibrium.jl for diverted plasmas because boundary calculation is wrong - S = Equilibrium.solovev(abs(B0), R0, epsilon, delta, kappa, x[1], x[2], B0_dir=sign(B0), Ip_dir=1, symmetric=true, x_point=nothing) - beta_cost = (Equilibrium.beta_n(S) - target_beta) / target_beta - ip_cost = (Equilibrium.plasma_current(S) - target_ip) / target_ip - c = beta_cost^2 + ip_cost^2 - return c - end - - res = Optim.optimize(cost, [alpha, qstar], Optim.NelderMead(), Optim.Options(g_tol=1E-3)) - - if par.verbose - println(res) - end - - actor.S = Equilibrium.solovev(abs(B0), R0, epsilon, delta, kappa, res.minimizer[1], res.minimizer[2], B0_dir=sign(B0), Ip_dir=1, symmetric=S0.symmetric, x_point=S0.x_point) - - return res -end - -""" - finalize( - actor::ActorSolovev; - rlims::NTuple{2,<:Real}=(maximum([actor.S.R0 * (1 - actor.S.epsilon * 2), 0.0]), actor.S.R0 * (1 + actor.S.epsilon * 2)), - zlims::NTuple{2,<:Real}=(-actor.S.R0 * actor.S.epsilon * actor.S.kappa * 1.7, actor.S.R0 * actor.S.epsilon * actor.S.kappa * 1.7) - )::IMAS.equilibrium__time_slice - -Store ActorSolovev data in IMAS.equilibrium format -""" -function finalize( - actor::ActorSolovev; - rlims::NTuple{2,<:Real}=(maximum([actor.S.R0 * (1 - actor.S.epsilon * 2), 0.0]), actor.S.R0 * (1 + actor.S.epsilon * 2)), - zlims::NTuple{2,<:Real}=(-actor.S.R0 * actor.S.epsilon * actor.S.kappa * 1.7, actor.S.R0 * actor.S.epsilon * actor.S.kappa * 1.7) -)::IMAS.equilibrium__time_slice - - ngrid = actor.par.ngrid - tc = transform_cocos(3, 11) - - eq = actor.eq - 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 = actor.S.R0 - @ddtime eq.vacuum_toroidal_field.b0 = actor.S.B0 * sign_Bt - - empty!(eqt) - - eqt.global_quantities.ip = ip - eqt.boundary.geometric_axis.r = actor.S.R0 - eqt.boundary.geometric_axis.z = Z0 - orig_psi = collect(range(Equilibrium.psi_limits(actor.S)..., length=ngrid)) - eqt.profiles_1d.psi = orig_psi * (tc["PSI"] * sign_Ip) - - eqt.profiles_1d.pressure = Equilibrium.pressure(actor.S, orig_psi) - eqt.profiles_1d.dpressure_dpsi = Equilibrium.pressure_gradient(actor.S, orig_psi) / (tc["PSI"] * sign_Ip) - - eqt.profiles_1d.f = Equilibrium.poloidal_current(actor.S, orig_psi) * (tc["F"] * sign_Bt) - eqt.profiles_1d.f_df_dpsi = Equilibrium.poloidal_current(actor.S, orig_psi) .* Equilibrium.poloidal_current_gradient(actor.S, orig_psi) * (tc["F_FPRIME"] * sign_Bt * sign_Ip) - - resize!(eqt.profiles_2d, 1) - eqt.profiles_2d[1].grid_type.index = 1 - eqt.profiles_2d[1].grid.dim1 = range(rlims[1], rlims[2], length=ngrid) - eqt.profiles_2d[1].grid.dim2 = range(zlims[1] + Z0, zlims[2] + Z0, length=Int(ceil(ngrid * actor.S.kappa))) - - eqt.profiles_2d[1].psi = [actor.S(rr, flip_z * (zz - Z0)) * (tc["PSI"] * sign_Ip) for rr in eqt.profiles_2d[1].grid.dim1, zz in eqt.profiles_2d[1].grid.dim2] - IMAS.flux_surfaces(eqt) - - # correct equilibrium volume and area - if !ismissing(actor.par, :volume) - eqt.profiles_1d.volume .*= actor.par.volume / eqt.profiles_1d.volume[end] - end - if !ismissing(actor.par, :area) - eqt.profiles_1d.area .*= actor.par.area / eqt.profiles_1d.area[end] - end - - return eqt -end - -""" - IMAS2Equilibrium(eqt::IMAS.equilibrium__time_slice) - -Convert IMAS.equilibrium__time_slice to Equilibrium.jl EFIT structure -""" -function IMAS2Equilibrium(eqt::IMAS.equilibrium__time_slice) - dim1 = range(eqt.profiles_2d[1].grid.dim1[1], eqt.profiles_2d[1].grid.dim1[end], length=length(eqt.profiles_2d[1].grid.dim1)) - @assert collect(dim1) ≈ eqt.profiles_2d[1].grid.dim1 - dim2 = range(eqt.profiles_2d[1].grid.dim2[1], eqt.profiles_2d[1].grid.dim2[end], length=length(eqt.profiles_2d[1].grid.dim2)) - @assert collect(dim2) ≈ eqt.profiles_2d[1].grid.dim2 - psi = range(eqt.profiles_1d.psi[1], eqt.profiles_1d.psi[end], length=length(eqt.profiles_1d.psi)) - @assert collect(psi) ≈ eqt.profiles_1d.psi - - Equilibrium.efit(Equilibrium.cocos(11), # COCOS - dim1, # Radius/R range - dim2, # Elevation/Z range - psi, # Polodial Flux range (polodial flux from magnetic axis) - eqt.profiles_2d[1].psi, # Polodial Flux on RZ grid (polodial flux from magnetic axis) - eqt.profiles_1d.f, # Polodial Current - eqt.profiles_1d.pressure, # Plasma pressure - eqt.profiles_1d.q, # Q profile - eqt.profiles_1d.psi .* 0, # Electric Potential - (eqt.global_quantities.magnetic_axis.r, eqt.global_quantities.magnetic_axis.z), # Magnetic Axis (raxis,zaxis) - Int(sign(eqt.profiles_1d.f[end]) * sign(eqt.global_quantities.ip)) # sign(dot(J,B)) - ) -end - -""" - gEQDSK2IMAS(GEQDSKFile::GEQDSKFile,eq::IMAS.equilibrium) - -Convert IMAS.equilibrium__time_slice to Equilibrium.jl EFIT structure -""" -function gEQDSK2IMAS(g::EFIT.GEQDSKFile, eq::IMAS.equilibrium) - - tc = transform_cocos(1, 11) # chease output is cocos 1 , dd is cocos 11 - - eqt = eq.time_slice[] - eq1d = eqt.profiles_1d - resize!(eqt.profiles_2d, 1) - eq2d = eqt.profiles_2d[1] - - @ddtime(eq.vacuum_toroidal_field.b0 = g.bcentr) - eq.vacuum_toroidal_field.r0 = g.rcentr - - eqt.global_quantities.magnetic_axis.r = g.rmaxis - eqt.boundary.geometric_axis.r = g.rcentr - eqt.boundary.geometric_axis.z = g.zmid - 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"] - - IMAS.flux_surfaces(eqt) -end +include("solovev_actor.jl") #= =========== =# # ActorCHEASE # #= =========== =# -Base.@kwdef mutable struct ActorCHEASE <: PlasmaAbstractActor - dd::IMAS.dd - par::ParametersActor - chease::Union{Nothing,CHEASE.Chease} -end - -function ParametersActor(::Type{Val{:ActorCHEASE}}) - par = ParametersActor(nothing) - par.free_boundary = Entry(Bool, "", "Convert fixed boundary equilibrium to free boundary one"; default=true) - par.clear_workdir = Entry(Bool, "", "Clean the temporary workdir for CHEASE"; default=true) - return par -end - -""" - ActorCHEASE(dd::IMAS.dd, act::ParametersAllActors; kw...) - -This actor runs the Fixed boundary equilibrium solver CHEASE""" -function ActorCHEASE(dd::IMAS.dd, act::ParametersAllActors; kw...) - par = act.ActorCHEASE(kw...) - actor = ActorCHEASE(dd, par) - step(actor) - finalize(actor) - return actor -end - -function ActorCHEASE(dd::IMAS.dd, par::ParametersActor) - ActorCHEASE(dd, deepcopy(par), nothing) -end - -""" - prepare(dd::IMAS.dd, :ActorCHEASE, act::ParametersAllActors; kw...) - -Prepare dd to run ActorCHEASE -* Copy pressure from core_profiles to equilibrium -* Copy j_parallel from core_profiles to equilibrium -""" -function prepare(dd::IMAS.dd, ::Type{Val{:ActorCHEASE}}, act::ParametersAllActors; kw...) - eq1d = dd.equilibrium.time_slice[].profiles_1d - cp1d = dd.core_profiles.profiles_1d[] - eq1d.j_tor = IMAS.interp1d(cp1d.grid.psi_norm, cp1d.j_tor).(eq1d.psi_norm) - eq1d.pressure = IMAS.interp1d(cp1d.grid.psi_norm, cp1d.pressure).(eq1d.psi_norm) -end - -""" - step(actor::ActorCHEASE) - -Runs CHEASE on the r_z boundary, equilibrium pressure and equilibrium j_tor -""" -function step(actor::ActorCHEASE) - dd = actor.dd - eqt = dd.equilibrium.time_slice[] - eq1d = eqt.profiles_1d - - r_bound = eqt.boundary.outline.r - z_bound = eqt.boundary.outline.z - index = abs.(IMAS.curvature(r_bound,z_bound)) .< 1.0 - r_bound = r_bound[index] - z_bound = z_bound[index] - - Ip = eqt.global_quantities.ip - Bt_center = @ddtime(dd.equilibrium.vacuum_toroidal_field.b0) - r_center = dd.equilibrium.vacuum_toroidal_field.r0 - - r_geo = eqt.boundary.geometric_axis.r - z_geo = eqt.boundary.geometric_axis.z - Bt_geo = Bt_center * r_center / r_geo - - ϵ = eqt.boundary.minor_radius / r_geo - - j_tor = eq1d.j_tor - pressure = eq1d.pressure - rho_pol = sqrt.(eq1d.psi_norm) - pressure_sep = pressure[end] - - actor.chease = CHEASE.run_chease(ϵ, z_geo, pressure_sep, Bt_geo, r_geo, Ip, r_bound, z_bound, 82, rho_pol, pressure, j_tor, clear_workdir=actor.par.clear_workdir) - - if actor.par.free_boundary - # convert from fixed to free boundary equilibrium - EQ = Equilibrium.efit(actor.chease.gfile, 1) - psi_free_rz = VacuumFields.fixed2free(EQ, actor.chease.gfile.nbbbs) - actor.chease.gfile.psirz = psi_free_rz - EQ = Equilibrium.efit(actor.chease.gfile, 1) - psi_b = Equilibrium.psi_boundary(EQ; r=EQ.r, z=EQ.z) - psi_a = EQ.psi_rz(EQ.axis...) - actor.chease.gfile.psirz = (psi_free_rz .- psi_a) * ((actor.chease.gfile.psi[end] - actor.chease.gfile.psi[1]) / (psi_b - psi_a)) .+ actor.chease.gfile.psi[1] - end - - return actor -end - -# define `finalize` function for this actor -function finalize(actor::ActorCHEASE) - gEQDSK2IMAS(actor.chease.gfile, actor.dd.equilibrium) - return actor -end \ No newline at end of file +include("chease_actor.jl") \ No newline at end of file diff --git a/src/actors/neutronics_actors.jl b/src/actors/neutronics_actors.jl index 8d77ea596..dbc18d683 100644 --- a/src/actors/neutronics_actors.jl +++ b/src/actors/neutronics_actors.jl @@ -1,3 +1,9 @@ +#= =============== =# +# ActorNeutronics # +#= =============== =# + +using MillerExtendedHarmonic + mutable struct neutron_particle x y @@ -15,12 +21,13 @@ function Zcoord(n::neutron_particle) n.z end -#= =============== =# -# ActorNeutronics # -#= =============== =# - -Base.@kwdef mutable struct ActorNeutronics <: PlasmaAbstractActor +mutable struct ActorNeutronics <: PlasmaAbstractActor dd::IMAS.dd + par::ParametersActor + function ActorNeutronics(dd::IMAS.dd, par::ParametersActor; kw...) + par = par(kw...) + return new(dd, par) + end end function ParametersActor(::Type{Val{:ActorNeutronics}}) @@ -41,13 +48,13 @@ This actor estimates the neutron loading on the wall using the fusion source fro """ function ActorNeutronics(dd::IMAS.dd, act::ParametersAllActors; kw...) par = act.ActorNeutronics(kw...) - actor = ActorNeutronics(dd) - step(actor; N=par.N, step=par.step, do_plot=par.do_plot) + actor = ActorNeutronics(dd, par) + step(actor) finalize(actor) return actor end -function step(actor::ActorNeutronics; N::Integer=100000, step=0.05, do_plot::Bool=false) +function step(actor::ActorNeutronics; N::Integer=actor.par.N, step=actor.par.step, do_plot::Bool=actor.par.do_plot) dd = actor.dd cp1d = dd.core_profiles.profiles_1d[] eqt = dd.equilibrium.time_slice[] @@ -81,14 +88,9 @@ function step(actor::ActorNeutronics; N::Integer=100000, step=0.05, do_plot::Boo # resample wall and make sure it's clockwise (for COCOS = 11) wall = IMAS.first_wall(dd.wall) wall_r, wall_z = IMAS.resample_2d_line(wall.r, wall.z; step=0.1) - #wall_r, wall_z = deepcopy(wall.r), deepcopy(wall.z) R0 = eqt.global_quantities.magnetic_axis.r Z0 = eqt.global_quantities.magnetic_axis.z - θ = unwrap(atan.(wall_z .- Z0, wall_r .- R0)) - if θ[1] < θ[end] - reverse!(wall_r) - reverse!(wall_z) - end + IMAS.reorder_flux_surface!(wall_r, wall_z, R0, Z0) # advance neutrons until they hit the wall rz_wall = collect(zip(wall_r, wall_z)) diff --git a/src/actors/pedestal_actors.jl b/src/actors/pedestal_actors.jl new file mode 100644 index 000000000..43d61960b --- /dev/null +++ b/src/actors/pedestal_actors.jl @@ -0,0 +1,117 @@ +import EPEDNN + +#= ============= =# +# ActorPedestal # +#= ============= =# +mutable struct ActorPedestal <: PlasmaAbstractActor + dd::IMAS.dd + par::ParametersActor + epedmod::EPEDNN.EPEDmodel + inputs::EPEDNN.InputEPED + wped::Union{Missing,Real} + pped::Union{Missing,Real} +end + +function ParametersActor(::Type{Val{:ActorPedestal}}) + par = ParametersActor(nothing) + par.temp_pedestal_ratio = Entry(Real, "", "Ratio of ion to electron temperatures"; default=1.0) + par.eped_factor = Entry(Real, "", "Pedestal height multiplier (affects width by the squareroot)"; default=1.0) + par.warn_nn_train_bounds = Entry(Bool, "", "Raise warnings if querying cases that are certainly outside of the training range"; default=false) + par.only_powerlaw = Entry(Bool, "", "Use power-law pedestal fit (without NN correction)"; default=false) + return par +end + +""" + ActorPedestal(dd::IMAS.dd, act::ParametersAllActors; kw...) + +The ActorPedestal evaluates the pedestal boundary condition (height and width) +""" +function ActorPedestal(dd::IMAS.dd, act::ParametersAllActors; kw...) + par = act.ActorPedestal(kw...) + actor = ActorPedestal(dd, par) + step(actor) + finalize(actor) + return actor +end + +function ActorPedestal(dd::IMAS.dd, par::ParametersActor; kw...) + par = par(kw...) + + epedmod = EPEDNN.loadmodelonce("EPED1NNmodel.bson") + + eq = dd.equilibrium + eqt = eq.time_slice[] + + m = [ion.element[1].a for ion in dd.core_profiles.profiles_1d[].ion if Int(floor(ion.element[1].z_n)) == 1] + m = sum(m) / length(m) + if m < 2 + m = 2 + elseif m > 2 + m = 2.5 + end + + neped = @ddtime dd.summary.local.pedestal.n_e.value + zeffped = @ddtime dd.summary.local.pedestal.zeff.value + Bt = abs(@ddtime(eq.vacuum_toroidal_field.b0)) * eq.vacuum_toroidal_field.r0 / eqt.boundary.geometric_axis.r + βn = @ddtime(dd.summary.global_quantities.beta_tor_thermal_norm.value) + + inputs = EPEDNN.InputEPED( + eqt.boundary.minor_radius, + βn, + Bt, + eqt.boundary.triangularity, + abs(eqt.global_quantities.ip / 1e6), + eqt.boundary.elongation, + m, + neped / 1e19, + eqt.boundary.geometric_axis.r, + zeffped) + + return ActorPedestal(dd, par, epedmod, inputs, missing, missing) +end + +""" + step(actor::ActorPedestal; + warn_nn_train_bounds::Bool=actor.par.warn_nn_train_bounds, + only_powerlaw::Bool=false) + +Runs pedestal actor to evaluate pedestal width and height +""" +function step(actor::ActorPedestal; + warn_nn_train_bounds::Bool=actor.par.warn_nn_train_bounds, + only_powerlaw::Bool=false) + + sol = actor.epedmod(actor.inputs; only_powerlaw, warn_nn_train_bounds) + + actor.wped = sol.width.GH.H + actor.pped = sol.pressure.GH.H + + return actor +end + +""" + finalize(actor::ActorPedestal; + temp_pedestal_ratio::Real=actor.par.temp_pedestal_ratio, + eped_factor::Real=actor.par.eped_factor) + +Writes results to dd.summary.local.pedestal +""" +function finalize(actor::ActorPedestal; + temp_pedestal_ratio::Real=actor.par.temp_pedestal_ratio, + eped_factor::Real=actor.par.eped_factor) + + dd = actor.dd + + impurity = [ion.element[1].z_n for ion in dd.core_profiles.profiles_1d[].ion if Int(floor(ion.element[1].z_n)) != 1][1] + zi = sum(impurity) / length(impurity) + + nival = actor.inputs.neped * (actor.inputs.zeffped - 1) / (zi^2 - zi) + nval = actor.inputs.neped - zi * nival + nsum = actor.inputs.neped + nval + nival + tped = actor.pped * 1e6 / nsum / 1.60218e-19 + + @ddtime dd.summary.local.pedestal.t_e.value = 2.0 * tped / (1.0 + temp_pedestal_ratio) * eped_factor + @ddtime dd.summary.local.pedestal.position.rho_tor_norm = 1 - actor.wped * sqrt(eped_factor) + + return actor +end diff --git a/src/actors/pf_active_actors.jl b/src/actors/pf_active_actors.jl index 8b5667892..30bd8784c 100644 --- a/src/actors/pf_active_actors.jl +++ b/src/actors/pf_active_actors.jl @@ -13,7 +13,8 @@ Base.@kwdef mutable struct PFcoilsOptTrace cost_total::Vector{Real} = Real[] end -Base.@kwdef mutable struct ActorPFcoilsOpt <: ReactorAbstractActor +mutable struct ActorPFcoilsOpt <: ReactorAbstractActor + par::ParametersActor eq_in::IMAS.equilibrium eq_out::IMAS.equilibrium pf_active::IMAS.pf_active @@ -34,10 +35,11 @@ function ParametersActor(::Type{Val{:ActorPFcoilsOpt}}) ] par.green_model = Switch(options, "", "Model used for the coils Green function calculations"; default=:simple) par.symmetric = Entry(Bool, "", "Force PF coils location to be up-down symmetric"; default=true) - par.λ_currents = Entry(Real, "", "Weight of current limit constraint"; default=0.5) - par.λ_strike = Entry(Real, "", "Weight given to matching the strike-points"; default=0.0) - par.λ_lcfs = Entry(Real, "", "Weight given to matching last closed flux surface"; default=1.0) - par.λ_null = Entry(Real, "", "Weight given to get field null for plasma breakdown"; default=1E-3) + par.weight_currents = Entry(Real, "", "Weight of current limit constraint"; default=0.5) + par.weight_strike = Entry(Real, "", "Weight given to matching the strike-points"; default=0.0) + par.weight_lcfs = Entry(Real, "", "Weight given to matching last closed flux surface"; default=1.0) + par.weight_null = Entry(Real, "", "Weight given to get field null for plasma breakdown"; default=1E-3) + par.maxiter = Entry(Integer, "", "Maximum number of optimizer iterations"; default=1000) options = [ :none => "Do not optimize", :currents => "Find optimial coil currents but do not change coil positions", @@ -61,7 +63,7 @@ to match the equilibrium boundary shape and obtain a field-null region at plasma """ function ActorPFcoilsOpt(dd::IMAS.dd, act::ParametersAllActors; kw...) par = act.ActorPFcoilsOpt(kw...) - actor = ActorPFcoilsOpt(dd; par.green_model, par.symmetric) + actor = ActorPFcoilsOpt(dd, par) if par.optimization_scheme == :none if par.do_plot @@ -73,13 +75,13 @@ function ActorPFcoilsOpt(dd::IMAS.dd, act::ParametersAllActors; kw...) else if par.optimization_scheme == :currents # find coil currents - step(actor; par.λ_lcfs, par.λ_null, par.λ_currents, par.λ_strike, par.verbose, maxiter=1000, optimization_scheme=:currents) - finalize(actor) + step(actor) + finalize(actor; update_equilibrium=false) elseif par.optimization_scheme == :rail # optimize coil location and currents - step(actor; par.λ_lcfs, par.λ_null, par.λ_currents, par.λ_strike, par.verbose, maxiter=1000, optimization_scheme=:rail) - finalize(actor) + step(actor) + finalize(actor; update_equilibrium=false) if par.do_plot display(plot(actor.trace, :cost)) @@ -96,106 +98,106 @@ function ActorPFcoilsOpt(dd::IMAS.dd, act::ParametersAllActors; kw...) display(plot(actor, equilibrium=true, time_index=length(dd.equilibrium.time))) end - finalize(actor; update_eq_in=par.update_equilibrium) + finalize(actor) end return actor end -function ActorPFcoilsOpt(dd::IMAS.dd; kw...) - return ActorPFcoilsOpt(dd.equilibrium, dd.build, dd.pf_active; kw...) -end - -function ActorPFcoilsOpt( - eq_in::IMAS.equilibrium, - bd::IMAS.build, - pf::IMAS.pf_active; - λ_regularize=1E-3, - green_model=:simple, - symmetric=false) +function ActorPFcoilsOpt(dd::IMAS.dd, par::ParametersActor; kw...) + par = par(kw...) + eq_in = dd.equilibrium + eq_out = deepcopy(eq_in) + pf = dd.pf_active + bd = dd.build + par.symmetric + λ_regularize=1E-3 + trace=PFcoilsOptTrace() + par.green_model # reset pf coil rails n_coils = [rail.coils_number for rail in bd.pf_active.rail] init_pf_active(pf, bd, n_coils) - # basic constructors - eq_out = deepcopy(eq_in) - - # constructor - pfactor = ActorPFcoilsOpt(eq_in, eq_out, pf, bd, symmetric, λ_regularize, PFcoilsOptTrace(), green_model) - - return pfactor + return ActorPFcoilsOpt(par, eq_in, eq_out, pf, bd, par.symmetric, λ_regularize, trace, par.green_model) end """ - step(pfactor::ActorPFcoilsOpt; - symmetric=pfactor.symmetric, - λ_regularize=pfactor.λ_regularize, - λ_lcfs=1.0, - λ_null=1E-3, - λ_currents=0.5, - λ_strike=0.0, - maxiter=10000, - optimization_scheme=:rail, - verbose=false) + step(actor::ActorPFcoilsOpt; + symmetric::Bool=actor.symmetric, + λ_regularize::Real=actor.λ_regularize, + weight_lcfs::Real=actor.par.weight_lcfs, + weight_null::Real=actor.par.weight_null, + weight_currents::Real=actor.par.weight_currents, + weight_strike::Real=actor.par.weight_strike, + maxiter::Real=actor.par.maxiter, + optimization_scheme::Symbol=actor.par.optimization_scheme, + verbose::Bool=actor.par.verbose) Optimize coil currents and positions to produce sets of equilibria while minimizing coil currents """ -function step(pfactor::ActorPFcoilsOpt; - symmetric=pfactor.symmetric, - λ_regularize=pfactor.λ_regularize, - λ_lcfs=1.0, - λ_null=1E-3, - λ_currents=0.5, - λ_strike=0.0, - maxiter=10000, - optimization_scheme=:rail, - verbose=false) - - fixed_coils, pinned_coils, optim_coils = fixed_pinned_optim_coils(pfactor, optimization_scheme) +function step(actor::ActorPFcoilsOpt; + symmetric::Bool=actor.symmetric, + λ_regularize::Real=actor.λ_regularize, + weight_lcfs::Real=actor.par.weight_lcfs, + weight_null::Real=actor.par.weight_null, + weight_currents::Real=actor.par.weight_currents, + weight_strike::Real=actor.par.weight_strike, + maxiter::Real=actor.par.maxiter, + optimization_scheme::Symbol=actor.par.optimization_scheme, + verbose::Bool=actor.par.verbose) + + fixed_coils, pinned_coils, optim_coils = fixed_pinned_optim_coils(actor, optimization_scheme) coils = vcat(pinned_coils, optim_coils, fixed_coils) for coil in coils - coil.current_time = pfactor.eq_in.time - coil.current_data = zeros(size(pfactor.eq_in.time)) + coil.current_time = actor.eq_in.time + coil.current_data = zeros(size(actor.eq_in.time)) end - bd = pfactor.bd + bd = actor.bd # run rail type optimizer if optimization_scheme in [:rail, :currents] - (λ_regularize, trace) = optimize_coils_rail(pfactor.eq_in; pinned_coils, optim_coils, fixed_coils, symmetric, λ_regularize, λ_lcfs, λ_null, λ_currents, λ_strike, bd, maxiter, verbose) + (λ_regularize, trace) = optimize_coils_rail(actor.eq_in; pinned_coils, optim_coils, fixed_coils, symmetric, λ_regularize, weight_lcfs, weight_null, weight_currents, weight_strike, bd, maxiter, verbose) else error("Supported ActorPFcoilsOpt optimization_scheme are `:currents` or `:rail`") end - pfactor.λ_regularize = λ_regularize - pfactor.trace = trace + actor.λ_regularize = λ_regularize + actor.trace = trace # transfer the results to IMAS.pf_active for coil in coils transfer_info_GS_coil_to_IMAS(bd, coil) end - return pfactor + return actor end """ - finalize(pfactor::ActorPFcoilsOpt; scale_eq_domain_size = 1.0) + finalize( + actor::ActorPFcoilsOpt; + update_equilibrium::Bool=actor.par.update_equilibrium, + scale_eq_domain_size::Real=1.0) -Update pfactor.eq_out 2D equilibrium PSI based on coils positions and currents +Update actor.eq_out 2D equilibrium PSI based on coils positions and currents """ -function finalize(pfactor::ActorPFcoilsOpt; scale_eq_domain_size=1.0, update_eq_in=false) +function finalize( + actor::ActorPFcoilsOpt; + update_equilibrium::Bool=actor.par.update_equilibrium, + scale_eq_domain_size::Real=1.0) + coils = GS_IMAS_pf_active__coil[] - for (k, coil) in enumerate(pfactor.pf_active.coil) - if k <= pfactor.bd.pf_active.rail[1].coils_number - coil_tech = pfactor.bd.oh.technology + for (k, coil) in enumerate(actor.pf_active.coil) + if k <= actor.bd.pf_active.rail[1].coils_number + coil_tech = actor.bd.oh.technology else - coil_tech = pfactor.bd.pf_active.technology + coil_tech = actor.bd.pf_active.technology end - push!(coils, GS_IMAS_pf_active__coil(coil, coil_tech, pfactor.green_model)) + push!(coils, GS_IMAS_pf_active__coil(coil, coil_tech, actor.green_model)) end # update equilibrium - for time_index in 1:length(pfactor.eq_in.time_slice) - if ismissing(pfactor.eq_in.time_slice[time_index].global_quantities, :ip) + for time_index in 1:length(actor.eq_in.time_slice) + if ismissing(actor.eq_in.time_slice[time_index].global_quantities, :ip) continue end for coil in coils @@ -203,24 +205,28 @@ function finalize(pfactor::ActorPFcoilsOpt; scale_eq_domain_size=1.0, update_eq_ end # convert equilibrium to Equilibrium.jl format, since this is what VacuumFields uses - EQfixed = IMAS2Equilibrium(pfactor.eq_in.time_slice[time_index]) + EQfixed = IMAS2Equilibrium(actor.eq_in.time_slice[time_index]) # # update ψ map R = range(EQfixed.r[1] / scale_eq_domain_size, EQfixed.r[end] * scale_eq_domain_size, length=length(EQfixed.r)) Z = range(EQfixed.z[1] * scale_eq_domain_size, EQfixed.z[end] * scale_eq_domain_size, length=length(EQfixed.z)) - ψ_f2f = VacuumFields.fixed2free(EQfixed, coils, R, Z) - pfactor.eq_out.time_slice[time_index].profiles_2d[1].grid.dim1 = R - pfactor.eq_out.time_slice[time_index].profiles_2d[1].grid.dim2 = Z - pfactor.eq_out.time_slice[time_index].profiles_2d[1].psi = transpose(ψ_f2f) + ψ_f2f = transpose(VacuumFields.fixed2free(EQfixed, coils, R, Z)) + actor.eq_out.time_slice[time_index].profiles_2d[1].grid.dim1 = R + actor.eq_out.time_slice[time_index].profiles_2d[1].grid.dim2 = Z + if false # hack to force up-down symmetric equilibrium + actor.eq_out.time_slice[time_index].profiles_2d[1].psi = (ψ_f2f .+ ψ_f2f[1:end, end:-1:1]) ./ 2.0 + else + actor.eq_out.time_slice[time_index].profiles_2d[1].psi = ψ_f2f + end end # update psi - if update_eq_in - for time_index in 1:length(pfactor.eq_out.time_slice) - if !ismissing(pfactor.eq_out.time_slice[time_index].global_quantities, :ip) - psi1 = pfactor.eq_out.time_slice[time_index].profiles_2d[1].psi - pfactor.eq_in.time_slice[time_index].profiles_2d[1].psi = psi1 - IMAS.flux_surfaces(pfactor.eq_in.time_slice[time_index]) + if update_equilibrium + for time_index in 1:length(actor.eq_out.time_slice) + if !ismissing(actor.eq_out.time_slice[time_index].global_quantities, :ip) + psi1 = actor.eq_out.time_slice[time_index].profiles_2d[1].psi + actor.eq_in.time_slice[time_index].profiles_2d[1].psi = psi1 + IMAS.flux_surfaces(actor.eq_in.time_slice[time_index]) end end end @@ -496,12 +502,12 @@ function optimize_coils_rail( fixed_coils::Vector{GS_IMAS_pf_active__coil}, symmetric::Bool, λ_regularize::Real, - λ_lcfs::Real, - λ_null::Real, - λ_currents::Real, - λ_strike::Real, + weight_lcfs::Real, + weight_null::Real, + weight_currents::Real, + weight_strike::Real, bd::IMAS.build, - maxiter::Int, + maxiter::Integer, verbose::Bool) R0 = eq.vacuum_toroidal_field.r0 @@ -524,7 +530,7 @@ function optimize_coils_rail( # private flux regions Rx = Float64[] Zx = Float64[] - if λ_strike > 0.0 + if weight_strike > 0.0 private = IMAS.flux_surface(eqt, eqt.profiles_1d.psi[end], false) vessel = IMAS.get_build(bd, type=_plasma_) for (pr, pz) in private @@ -533,15 +539,18 @@ function optimize_coils_rail( append!(Zx, pvy) end if isempty(Rx) - @warn "λ_strike>0 but no strike point found" + @warn "weight_strike>0 but no strike point found" end end # find ψp Bp_fac, ψp, Rp, Zp = VacuumFields.ψp_on_fixed_eq_boundary(fixed_eq, fixed_coils; Rx, Zx) push!(fixed_eqs, (Bp_fac, ψp, Rp, Zp)) + # weight more near the x-points + h = (maximum(Zp) - minimum(Zp)) / 2.0 + o = (maximum(Zp) + minimum(Zp)) / 2.0 + weight = sqrt.(((Zp .- o) ./ h) .^ 2 .+ h) / h # give each strike point the same weight as the lcfs - weight = Rp .* 0.0 .+ 1.0 - weight[end-length(Rx)+1:end] .= length(Rp) / (1 + length(Rx)) * λ_strike + weight[end-length(Rx)+1:end] .= length(Rp) / (1 + length(Rx)) * weight_strike if all(weight .== 1.0) weight = Float64[] end @@ -560,7 +569,7 @@ function optimize_coils_rail( index = findall(.>(1.0), abs.(packed[1:end-1])) if length(index) > 0 - cost_1to1 = sum(abs.(packed[index]) .- 1.0) * 10 + cost_1to1 = sum(abs.(packed[index]) .- 1.0) else cost_1to1 = 0.0 end @@ -582,17 +591,22 @@ function optimize_coils_rail( max_current_densities = vcat(oh_max_current_densities, pf_max_current_densities) fraction_max_current_densities = abs.(current_densities ./ max_current_densities) #currents cost - push!(all_cost_currents, norm(exp.(fraction_max_current_densities / λ_currents) / exp(1)) / length(currents)) - # boundary cost + push!(all_cost_currents, norm(exp.(fraction_max_current_densities * weight_currents) / exp(1)) / length(currents)) + # boundary and oh costs if ismissing(eq.time_slice[time_index].global_quantities, :ip) - push!(all_cost_lcfs, cost_lcfs0 / λ_null) + push!(all_cost_lcfs, cost_lcfs0 * weight_null) push!(all_cost_oh, 0.0) else - #OH cost - oh_current_densities = current_densities[oh_indexes] - avg_oh = sum(oh_current_densities) / length(oh_current_densities) - cost_oh = norm(oh_current_densities .- avg_oh) / avg_oh - push!(all_cost_lcfs, cost_lcfs0 / λ_lcfs) + push!(all_cost_lcfs, cost_lcfs0 * weight_lcfs) + # OH cost (to avoid OH shrinking too much) + index = findfirst(rail -> rail.name === "OH", bd.pf_active.rail) + if index !== nothing + oh_rail_length = diff(collect(extrema(bd.pf_active.rail[index].outline.z)))[1] + total_oh_coils_length = sum([coil.height for coil in coils[oh_indexes.==true]]) + cost_oh = oh_rail_length / total_oh_coils_length + else + cost_oh = 0.0 + end push!(all_cost_oh, cost_oh) end end @@ -608,13 +622,20 @@ function optimize_coils_rail( end end end - cost_spacing = cost_spacing / R0 + cost_spacing = cost_spacing / length(optim_coils)^2 / R0 + + cost_lcfs_2 = cost_lcfs^2 * 10000.0 + cost_currents_2 = cost_currents^2 + cost_oh_2 = cost_oh^2 + cost_1to1_2 = cost_1to1^2 + cost_spacing_2 = cost_spacing^2 + # total cost - cost = sqrt(cost_lcfs^2 + cost_currents^2 + 0.1 * cost_oh^2 + cost_1to1^2 + 10 * cost_spacing^2) + cost = sqrt(cost_lcfs_2 + cost_currents_2 + cost_oh_2 + cost_1to1_2 + cost_spacing_2) if do_trace push!(trace.params, packed) - push!(trace.cost_lcfs, cost_lcfs) - push!(trace.cost_currents, cost_currents) + push!(trace.cost_lcfs, sqrt(cost_lcfs_2)) + push!(trace.cost_currents, sqrt(cost_currents_2)) push!(trace.cost_total, cost) end if isnan(cost) @@ -660,31 +681,31 @@ function optimize_coils_rail( end """ - fixed_pinned_optim_coils(pfactor, optimization_scheme) + fixed_pinned_optim_coils(actor, optimization_scheme) Returns tuple of GS_IMAS_pf_active__coil coils organized by their function: - fixed: fixed position and current - pinned: coils with fixed position but current is optimized - optim: coils that have theri position and current optimized """ -function fixed_pinned_optim_coils(pfactor, optimization_scheme) +function fixed_pinned_optim_coils(actor, optimization_scheme) fixed_coils = GS_IMAS_pf_active__coil[] pinned_coils = GS_IMAS_pf_active__coil[] optim_coils = GS_IMAS_pf_active__coil[] - for (k, coil) in enumerate(pfactor.pf_active.coil) - if k <= pfactor.bd.pf_active.rail[1].coils_number - coil_tech = pfactor.bd.oh.technology + for (k, coil) in enumerate(actor.pf_active.coil) + if k <= actor.bd.pf_active.rail[1].coils_number + coil_tech = actor.bd.oh.technology else - coil_tech = pfactor.bd.pf_active.technology + coil_tech = actor.bd.pf_active.technology end if coil.identifier == "pinned" - push!(pinned_coils, GS_IMAS_pf_active__coil(coil, coil_tech, pfactor.green_model)) + push!(pinned_coils, GS_IMAS_pf_active__coil(coil, coil_tech, actor.green_model)) elseif (coil.identifier == "optim") && (optimization_scheme == :currents) - push!(pinned_coils, GS_IMAS_pf_active__coil(coil, coil_tech, pfactor.green_model)) + push!(pinned_coils, GS_IMAS_pf_active__coil(coil, coil_tech, actor.green_model)) elseif coil.identifier == "optim" - push!(optim_coils, GS_IMAS_pf_active__coil(coil, coil_tech, pfactor.green_model)) + push!(optim_coils, GS_IMAS_pf_active__coil(coil, coil_tech, actor.green_model)) elseif coil.identifier == "fixed" - push!(fixed_coils, GS_IMAS_pf_active__coil(coil, coil_tech, pfactor.green_model)) + push!(fixed_coils, GS_IMAS_pf_active__coil(coil, coil_tech, actor.green_model)) else error("Accepted type of coil.identifier are only \"optim\", \"pinned\", or \"fixed\"") end @@ -696,12 +717,12 @@ end # plotting # #= ======== =# """ - plot_pfcoilsactor_cx(pfactor::ActorPFcoilsOpt; time_index=1, equilibrium=true, rail=true) + plot_pfcoilsactor_cx(actor::ActorPFcoilsOpt; time_index=1, equilibrium=true, rail=true) Plot ActorPFcoilsOpt optimization cross-section """ @recipe function plot_ActorPFcoilsOpt_cx( - pfactor::ActorPFcoilsOpt; + actor::ActorPFcoilsOpt; time_index=nothing, equilibrium=true, build=true, @@ -717,13 +738,13 @@ Plot ActorPFcoilsOpt optimization cross-section @assert typeof(plot_r_buffer) <: Real if time_index === nothing - time_index = length(pfactor.eq_out.time_slice) + time_index = length(actor.eq_out.time_slice) end - time = pfactor.eq_out.time_slice[time_index].time + time = actor.eq_out.time_slice[time_index].time # if there is no equilibrium then treat this as a field_null plot field_null = false - if length(pfactor.eq_out.time_slice[time_index].profiles_2d) == 0 || ismissing(pfactor.eq_out.time_slice[time_index].profiles_2d[1], :psi) + if length(actor.eq_out.time_slice[time_index].profiles_2d) == 0 || ismissing(actor.eq_out.time_slice[time_index].profiles_2d[1], :psi) coils_flux = equilibrium field_null = true end @@ -734,8 +755,8 @@ Plot ActorPFcoilsOpt optimization cross-section end # setup plotting area - xlim = [0.0, maximum(pfactor.bd.layer[end].outline.r)] - ylim = [minimum(pfactor.bd.layer[end].outline.z), maximum(pfactor.bd.layer[end].outline.z)] + xlim = [0.0, maximum(actor.bd.layer[end].outline.r)] + ylim = [minimum(actor.bd.layer[end].outline.z), maximum(actor.bd.layer[end].outline.z)] xlim --> xlim * plot_r_buffer ylim --> ylim aspect_ratio --> :equal @@ -744,7 +765,7 @@ Plot ActorPFcoilsOpt optimization cross-section if build @series begin exclude_layers --> [:oh] - pfactor.bd + actor.bd end end @@ -755,19 +776,19 @@ Plot ActorPFcoilsOpt optimization cross-section Z = range(ylim[1], ylim[2], length=Int(ceil(ngrid * (ylim[2] - ylim[1]) / (xlim[2] - xlim[1])))) coils = GS_IMAS_pf_active__coil[] - for (k, coil) in enumerate(pfactor.pf_active.coil) - if k <= pfactor.bd.pf_active.rail[1].coils_number - coil_tech = pfactor.bd.oh.technology + for (k, coil) in enumerate(actor.pf_active.coil) + if k <= actor.bd.pf_active.rail[1].coils_number + coil_tech = actor.bd.oh.technology else - coil_tech = pfactor.bd.pf_active.technology + coil_tech = actor.bd.pf_active.technology end - coil = GS_IMAS_pf_active__coil(coil, coil_tech, pfactor.green_model) + coil = GS_IMAS_pf_active__coil(coil, coil_tech, actor.green_model) coil.time_index = time_index push!(coils, coil) end # ψ coil currents - ψbound = pfactor.eq_out.time_slice[time_index].global_quantities.psi_boundary + ψbound = actor.eq_out.time_slice[time_index].global_quantities.psi_boundary ψ = VacuumFields.coils_flux(2 * pi, coils, R, Z) ψmin = minimum(x -> isnan(x) ? Inf : x, ψ) @@ -802,7 +823,7 @@ Plot ActorPFcoilsOpt optimization cross-section @series begin wireframe --> true exclude_layers --> [:oh] - pfactor.bd + actor.bd end end @@ -813,14 +834,14 @@ Plot ActorPFcoilsOpt optimization cross-section cx := true label --> "Field null region" seriescolor --> :red - pfactor.eq_out.time_slice[time_index] + actor.eq_out.time_slice[time_index] end else @series begin cx := true label --> "Final" seriescolor --> :red - pfactor.eq_out.time_slice[time_index] + actor.eq_out.time_slice[time_index] end @series begin cx := true @@ -828,7 +849,7 @@ Plot ActorPFcoilsOpt optimization cross-section seriescolor --> :blue lcfs --> true linestyle --> :dash - pfactor.eq_in.time_slice[time_index] + actor.eq_in.time_slice[time_index] end end end @@ -836,14 +857,14 @@ Plot ActorPFcoilsOpt optimization cross-section # plot pf_active coils @series begin time --> time - pfactor.pf_active + actor.pf_active end # plot optimization rails if rail @series begin label --> (build ? "Coil opt. rail" : "") - pfactor.bd.pf_active.rail + actor.bd.pf_active.rail end end diff --git a/src/actors/pf_passive_actors.jl b/src/actors/pf_passive_actors.jl new file mode 100644 index 000000000..1374e083b --- /dev/null +++ b/src/actors/pf_passive_actors.jl @@ -0,0 +1,77 @@ +#= ========== =# +# PF passive # +#= ========== =# + +mutable struct ActorPassiveStructures <: ReactorAbstractActor + dd::IMAS.dd + par::ParametersActor + function ActorPassiveStructures(dd::IMAS.dd, par::ParametersActor; kw...) + par = par(kw...) + return new(dd, par) + end +end + +function ParametersActor(::Type{Val{:ActorPassiveStructures}}) + par = ParametersActor(nothing) + par.do_plot = Entry(Bool, "", "plot"; default=false) + return par +end + +""" + ActorPassiveStructures(dd::IMAS.dd, act::ParametersAllActors; kw...) + +Populates `pf_passive` structures +""" +function ActorPassiveStructures(dd::IMAS.dd, act::ParametersAllActors; kw...) + par = act.ActorPassiveStructures(kw...) + actor = ActorPassiveStructures(dd, par) + step(actor) + finalize(actor) + if par.do_plot + display(plot(dd.pf_passive)) + end + return actor +end + +function step(actor::ActorPassiveStructures) + dd = actor.dd + + # all LFS + ilayers = IMAS.get_build(dd.build, fs=IMAS._lfs_, return_only_one=false, return_index=true) + ilayers = vcat(ilayers[1] - 1, ilayers) + for k in ilayers + l = dd.build.layer[k] + l1 = dd.build.layer[k+1] + if l1.material == "Vacuum" + continue + end + poly = IMAS.join_outlines(l.outline.r, l.outline.z, l1.outline.r, l1.outline.z) + add_pf_passive_loop(dd.pf_passive, l1, poly[1], poly[2]) + end + + # OH and plug + add_pf_passive_loop(dd.pf_passive, dd.build.layer[2], dd.build.layer[2].outline.r, dd.build.layer[2].outline.z) + if dd.build.layer[1].material != "Vacuum" + add_pf_passive_loop(dd.pf_passive, dd.build.layer[1], dd.build.layer[1].outline.r, dd.build.layer[1].outline.z) + end + + # cryostat + layer = dd.build.layer[end] + if layer.type == Int(IMAS._cryostat_) + i = findfirst(layer.outline.r .== 0) + r = vcat(layer.outline.r[i+1:end-1], layer.outline.r[1:i]) + z = vcat(layer.outline.z[i+1:end-1], layer.outline.z[1:i]) + layer = dd.build.layer[end-1] + i = findfirst(layer.outline.r .== 0) + r = vcat(layer.outline.r[i+1:end-1], layer.outline.r[1:i], reverse(r), layer.outline.r[i+1:i+1]) + z = vcat(layer.outline.z[i+1:end-1], layer.outline.z[1:i], reverse(z), layer.outline.z[i+1:i+1]) + add_pf_passive_loop(dd.pf_passive, dd.build.layer[end], r, z) + end +end + +function add_pf_passive_loop(pf_passive::IMAS.pf_passive, layer::IMAS.build__layer, r::AbstractVector{T}, z::AbstractVector{T}) where {T<:Real} + resize!(resize!(pf_passive.loop, length(pf_passive.loop) + 1)[end].element, 1) + pf_passive.loop[end].name = replace(layer.name, r"[lh]fs " => "") + pf_passive.loop[end].element[end].geometry.outline.r = r + pf_passive.loop[end].element[end].geometry.outline.z = z +end \ No newline at end of file diff --git a/src/actors/solovev_actor.jl b/src/actors/solovev_actor.jl new file mode 100644 index 000000000..142c09463 --- /dev/null +++ b/src/actors/solovev_actor.jl @@ -0,0 +1,216 @@ +import Equilibrium +import EFIT +import ForwardDiff +import Optim + +#= ============ =# +# ActorSolovev # +#= ============ =# +mutable struct ActorSolovev <: PlasmaAbstractActor + eq::IMAS.equilibrium + par::ParametersActor + S::Equilibrium.SolovevEquilibrium +end + +function ParametersActor(::Type{Val{:ActorSolovev}}) + par = ParametersActor(nothing) + par.ngrid = Entry(Integer, "", "Grid size (for R, Z follows proportionally to plasma elongation)"; default=129) + par.qstar = Entry(Real, "", "Initial guess of kink safety factor"; default=1.5) + par.alpha = Entry(Real, "", "Initial guess of constant relating to pressure"; default=0.0) + par.volume = Entry(Real, "m³", "Scalar volume to match (optional)"; default=missing) + par.area = Entry(Real, "m²", "Scalar area to match (optional)"; default=missing) + par.verbose = Entry(Bool, "", "verbose"; default=false) + return par +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...) + par = act.ActorSolovev(kw...) + actor = ActorSolovev(dd, par) + step(actor) + finalize(actor) + # record optimized values of qstar and alpha in `act` for subsequent ActorSolovev calls + par.qstar = actor.S.qstar + par.alpha = actor.S.alpha + return actor +end + +function ActorSolovev(dd::IMAS.dd, par::ParametersActor; kw...) + par = par(kw...) + + # extract info from dd + eq = dd.equilibrium + eqt = eq.time_slice[] + a = eqt.boundary.minor_radius + R0 = eqt.boundary.geometric_axis.r + κ = eqt.boundary.elongation + δ = eqt.boundary.triangularity + ϵ = a / R0 + B0 = @ddtime eq.vacuum_toroidal_field.b0 + + # 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)) + else + x_point = nothing + end + + # run Solovev + S = Equilibrium.solovev(abs(B0), R0, ϵ, δ, κ, par.alpha, par.qstar, B0_dir=Int64(sign(B0)), Ip_dir=1, x_point=x_point, symmetric=symmetric) + + return ActorSolovev(dd.equilibrium, par, S) +end + +""" + prepare(dd::IMAS.dd, :ActorSolovev, act::ParametersAllActors; kw...) + +Prepare dd to run ActorSolovev +* Copy pressure from core_profiles to equilibrium +""" +function prepare(dd::IMAS.dd, ::Type{Val{:ActorSolovev}}, act::ParametersAllActors; kw...) + eq1d = dd.equilibrium.time_slice[].profiles_1d + cp1d = dd.core_profiles.profiles_1d[] + eq1d.pressure = IMAS.interp1d(cp1d.grid.psi_norm, cp1d.pressure).(eq1d.psi_norm) +end + +""" + step(actor::ActorSolovev) + +Non-linear optimization to obtain a target `ip` and `pressure_core` +""" +function step(actor::ActorSolovev) + S0 = actor.S + par = actor.par + + eqt = actor.eq.time_slice[] + target_ip = abs(eqt.global_quantities.ip) + target_pressure_core = eqt.profiles_1d.pressure[1] + + B0, R0, epsilon, delta, kappa, alpha, qstar, target_ip, target_pressure_core = promote(S0.B0, S0.R0, S0.epsilon, S0.delta, S0.kappa, S0.alpha, S0.qstar, target_ip, target_pressure_core) + + function cost(x) + # NOTE: Ip/pressure calculation is very much off in Equilibrium.jl for diverted plasmas because boundary calculation is wrong + S = Equilibrium.solovev(abs(B0), R0, epsilon, delta, kappa, x[1], x[2], B0_dir=sign(B0), Ip_dir=1, symmetric=true, x_point=nothing) + psimag, psibry = Equilibrium.psi_limits(S) + pressure_cost = (Equilibrium.pressure(S, psimag) - target_pressure_core) / target_pressure_core + ip_cost = (Equilibrium.plasma_current(S) - target_ip) / target_ip + c = pressure_cost^2 + ip_cost^2 + return c + end + + res = Optim.optimize(cost, [alpha, qstar], Optim.NelderMead(), Optim.Options(g_tol=1E-3)) + + if par.verbose + println(res) + end + + actor.S = Equilibrium.solovev(abs(B0), R0, epsilon, delta, kappa, res.minimizer[1], res.minimizer[2], B0_dir=sign(B0), Ip_dir=1, symmetric=S0.symmetric, x_point=S0.x_point) + + return actor +end + +""" + finalize( + actor::ActorSolovev; + rlims::NTuple{2,<:Real}=(maximum([actor.S.R0 * (1 - actor.S.epsilon * 2), 0.0]), actor.S.R0 * (1 + actor.S.epsilon * 2)), + zlims::NTuple{2,<:Real}=(-actor.S.R0 * actor.S.epsilon * actor.S.kappa * 1.7, actor.S.R0 * actor.S.epsilon * actor.S.kappa * 1.7) + )::IMAS.equilibrium__time_slice + +Store ActorSolovev data in IMAS.equilibrium format +""" +function finalize( + actor::ActorSolovev; + rlims::NTuple{2,<:Real}=(maximum([actor.S.R0 * (1 - actor.S.epsilon * 2), 0.0]), actor.S.R0 * (1 + actor.S.epsilon * 2)), + zlims::NTuple{2,<:Real}=(-actor.S.R0 * actor.S.epsilon * actor.S.kappa * 1.7, actor.S.R0 * actor.S.epsilon * actor.S.kappa * 1.7) +)::IMAS.equilibrium__time_slice + + ngrid = actor.par.ngrid + tc = transform_cocos(3, 11) + + eq = actor.eq + 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 = actor.S.R0 + @ddtime eq.vacuum_toroidal_field.b0 = actor.S.B0 * sign_Bt + + empty!(eqt) + + eqt.global_quantities.ip = ip + eqt.boundary.geometric_axis.r = actor.S.R0 + eqt.boundary.geometric_axis.z = Z0 + orig_psi = collect(range(Equilibrium.psi_limits(actor.S)..., length=ngrid)) + eqt.profiles_1d.psi = orig_psi * (tc["PSI"] * sign_Ip) + + eqt.profiles_1d.pressure = Equilibrium.pressure(actor.S, orig_psi) + eqt.profiles_1d.dpressure_dpsi = Equilibrium.pressure_gradient(actor.S, orig_psi) / (tc["PSI"] * sign_Ip) + + eqt.profiles_1d.f = Equilibrium.poloidal_current(actor.S, orig_psi) * (tc["F"] * sign_Bt) + eqt.profiles_1d.f_df_dpsi = Equilibrium.poloidal_current(actor.S, orig_psi) .* Equilibrium.poloidal_current_gradient(actor.S, orig_psi) * (tc["F_FPRIME"] * sign_Bt * sign_Ip) + + resize!(eqt.profiles_2d, 1) + eqt.profiles_2d[1].grid_type.index = 1 + eqt.profiles_2d[1].grid.dim1 = range(rlims[1], rlims[2], length=ngrid) + eqt.profiles_2d[1].grid.dim2 = range(zlims[1] + Z0, zlims[2] + Z0, length=Int(ceil(ngrid * actor.S.kappa))) + + eqt.profiles_2d[1].psi = [actor.S(rr, flip_z * (zz - Z0)) * (tc["PSI"] * sign_Ip) for rr in eqt.profiles_2d[1].grid.dim1, zz in eqt.profiles_2d[1].grid.dim2] + IMAS.flux_surfaces(eqt) + + # correct equilibrium volume and area + if !ismissing(actor.par, :volume) + eqt.profiles_1d.volume .*= actor.par.volume / eqt.profiles_1d.volume[end] + end + if !ismissing(actor.par, :area) + eqt.profiles_1d.area .*= actor.par.area / eqt.profiles_1d.area[end] + end + + return eqt +end + +""" + IMAS2Equilibrium(eqt::IMAS.equilibrium__time_slice) + +Convert IMAS.equilibrium__time_slice to Equilibrium.jl EFIT structure +""" +function IMAS2Equilibrium(eqt::IMAS.equilibrium__time_slice) + dim1 = range(eqt.profiles_2d[1].grid.dim1[1], eqt.profiles_2d[1].grid.dim1[end], length=length(eqt.profiles_2d[1].grid.dim1)) + @assert collect(dim1) ≈ eqt.profiles_2d[1].grid.dim1 + dim2 = range(eqt.profiles_2d[1].grid.dim2[1], eqt.profiles_2d[1].grid.dim2[end], length=length(eqt.profiles_2d[1].grid.dim2)) + @assert collect(dim2) ≈ eqt.profiles_2d[1].grid.dim2 + psi = range(eqt.profiles_1d.psi[1], eqt.profiles_1d.psi[end], length=length(eqt.profiles_1d.psi)) + @assert collect(psi) ≈ eqt.profiles_1d.psi + + Equilibrium.efit(Equilibrium.cocos(11), # COCOS + dim1, # Radius/R range + dim2, # Elevation/Z range + psi, # Polodial Flux range (polodial flux from magnetic axis) + eqt.profiles_2d[1].psi, # Polodial Flux on RZ grid (polodial flux from magnetic axis) + eqt.profiles_1d.f, # Polodial Current + eqt.profiles_1d.pressure, # Plasma pressure + eqt.profiles_1d.q, # Q profile + eqt.profiles_1d.psi .* 0, # Electric Potential + (eqt.global_quantities.magnetic_axis.r, eqt.global_quantities.magnetic_axis.z), # Magnetic Axis (raxis,zaxis) + Int(sign(eqt.profiles_1d.f[end]) * sign(eqt.global_quantities.ip)) # sign(dot(J,B)) + ) +end diff --git a/src/actors/sources_actors.jl b/src/actors/sources_actors.jl index 996060f04..7150ac3be 100644 --- a/src/actors/sources_actors.jl +++ b/src/actors/sources_actors.jl @@ -3,18 +3,19 @@ import NumericalIntegration: integrate #= === =# # NBI # #= === =# -Base.@kwdef mutable struct ActorNBIsimple <: HCDAbstractActor +mutable struct ActorNBIsimple <: HCDAbstractActor dd::IMAS.dd - width::Vector{Real} - rho_0::Vector{Real} - current_efficiency::Vector{Real} + par::ParametersActor + width::AbstractVector{<:Real} + rho_0::AbstractVector{<:Real} + current_efficiency::AbstractVector{<:Real} end function ParametersActor(::Type{Val{:ActorNBIsimple}}) par = ParametersActor(nothing) - par.width = Entry(Real, "", "Width of the deposition profile"; default=0.3) - par.rho_0 = Entry(Real, "", "Radial location of the deposition profile"; default=0.0) - par.current_efficiency = Entry(Real, "A/W", "Current drive efficiency"; default=0.3) + par.width = Entry(Union{Real,AbstractVector{<:Real}}, "", "Width of the deposition profile"; default=0.3) + par.rho_0 = Entry(Union{Real,AbstractVector{<:Real}}, "", "Radial location of the deposition profile"; default=0.0) + par.current_efficiency = Entry(Union{Real,AbstractVector{<:Real}}, "A/W", "Current drive efficiency"; default=0.3) return par end @@ -28,15 +29,17 @@ This actor estimates the NBI ion/electron energy deposition, particle source, ro """ function ActorNBIsimple(dd::IMAS.dd, act::ParametersAllActors; kw...) par = act.ActorNBIsimple(kw...) - actor = ActorNBIsimple(dd; par.width, par.rho_0, par.current_efficiency) + actor = ActorNBIsimple(dd, par) step(actor) finalize(actor) return actor end -function ActorNBIsimple(dd::IMAS.dd; width::Real=0.3, rho_0::Real=0.0, current_efficiency::Real=0.3) - nbeam = ones(length(dd.nbi.unit)) - return ActorNBIsimple(dd, nbeam .* width, nbeam .* rho_0, nbeam .* current_efficiency) +function ActorNBIsimple(dd::IMAS.dd, par::ParametersActor; kw...) + par = par(kw...) + n_beams = length(dd.nbi.unit) + _, width, rho_0, current_efficiency = same_length_vectors(1:n_beams, par.width, par.rho_0, par.current_efficiency) + return ActorNBIsimple(dd, par, width, rho_0, current_efficiency) end function step(actor::ActorNBIsimple) @@ -87,18 +90,19 @@ end #= == =# # EC # #= == =# -Base.@kwdef mutable struct ActorECsimple <: HCDAbstractActor +mutable struct ActorECsimple <: HCDAbstractActor dd::IMAS.dd - width::Vector{Real} - rho_0::Vector{Real} - current_efficiency::Vector{Real} + par::ParametersActor + width::AbstractVector{<:Real} + rho_0::AbstractVector{<:Real} + current_efficiency::AbstractVector{<:Real} end function ParametersActor(::Type{Val{:ActorECsimple}}) par = ParametersActor(nothing) - par.width = Entry(Real, "", "Width of the deposition profile"; default=0.1) - par.rho_0 = Entry(Real, "", "Radial location of the deposition profile"; default=0.0) - par.current_efficiency = Entry(Real, "A/W", "Current drive efficiency"; default=0.2) + par.width = Entry(Union{Real,AbstractVector{<:Real}}, "", "Width of the deposition profile"; default=0.1) + par.rho_0 = Entry(Union{Real,AbstractVector{<:Real}}, "", "Radial location of the deposition profile"; default=0.0) + par.current_efficiency = Entry(Union{Real,AbstractVector{<:Real}}, "A/W", "Current drive efficiency"; default=0.2) return par end @@ -112,19 +116,21 @@ This actor estimates the EC electron energy deposition and current drive as a ga """ function ActorECsimple(dd::IMAS.dd, act::ParametersAllActors; kw...) par = act.ActorECsimple(kw...) - actor = ActorECsimple(dd; par.width, par.rho_0, par.current_efficiency) + actor = ActorECsimple(dd, par) step(actor) finalize(actor) return actor end -function ActorECsimple(dd::IMAS.dd; width::Real=0.1, rho_0::Real=0.0, current_efficiency::Real=0.2) - n_launchers = ones(length(dd.ec_launchers.launcher)) - return ActorECsimple(dd, n_launchers .* width, n_launchers .* rho_0, n_launchers .* current_efficiency) +function ActorECsimple(dd::IMAS.dd, par::ParametersActor; kw...) + par = par(kw...) + n_launchers = length(dd.ec_launchers.beam) + _, width, rho_0, current_efficiency = same_length_vectors(1:n_launchers, par.width, par.rho_0, par.current_efficiency) + return ActorECsimple(dd, par, width, rho_0, current_efficiency) end function step(actor::ActorECsimple) - for (idx, ecl) in enumerate(actor.dd.ec_launchers.launcher) + for (idx, ecl) in enumerate(actor.dd.ec_launchers.beam) eqt = actor.dd.equilibrium.time_slice[] cp1d = actor.dd.core_profiles.profiles_1d[] cs = actor.dd.core_sources @@ -163,18 +169,19 @@ end #= == =# # IC # #= == =# -Base.@kwdef mutable struct ActorICsimple <: HCDAbstractActor +mutable struct ActorICsimple <: HCDAbstractActor dd::IMAS.dd - width::Vector{Real} - rho_0::Vector{Real} - current_efficiency::Vector{Real} + par::ParametersActor + width::AbstractVector{<:Real} + rho_0::AbstractVector{<:Real} + current_efficiency::AbstractVector{<:Real} end function ParametersActor(::Type{Val{:ActorICsimple}}) par = ParametersActor(nothing) - par.width = Entry(Real, "", "Width of the deposition profile"; default=0.1) - par.rho_0 = Entry(Real, "", "Radial location of the deposition profile"; default=0.0) - par.current_efficiency = Entry(Real, "A/W", "Current drive efficiency"; default=0.125) + par.width = Entry(Union{Real,AbstractVector{<:Real}}, "", "Width of the deposition profile"; default=0.1) + par.rho_0 = Entry(Union{Real,AbstractVector{<:Real}}, "", "Radial location of the deposition profile"; default=0.0) + par.current_efficiency = Entry(Union{Real,AbstractVector{<:Real}}, "A/W", "Current drive efficiency"; default=0.125) return par end @@ -188,15 +195,17 @@ This actor estimates the ion-cyclotron electron/ion energy deposition and curren """ function ActorICsimple(dd::IMAS.dd, act::ParametersAllActors; kw...) par = act.ActorICsimple(kw...) - actor = ActorICsimple(dd; par.width, par.rho_0, par.current_efficiency) + actor = ActorICsimple(dd, par) step(actor) finalize(actor) return actor end -function ActorICsimple(dd::IMAS.dd; width::Real=0.1, rho_0::Real=0.0, current_efficiency::Real=0.125) - n_antennas = ones(length(dd.ic_antennas.antenna)) - return ActorICsimple(dd, n_antennas .* width, n_antennas .* rho_0, n_antennas .* current_efficiency) +function ActorICsimple(dd::IMAS.dd, par::ParametersActor; kw...) + par = par(kw...) + n_antennas = length(dd.ic_antennas.antenna) + _, width, rho_0, current_efficiency = same_length_vectors(1:n_antennas, par.width, par.rho_0, par.current_efficiency) + return ActorICsimple(dd, par, width, rho_0, current_efficiency) end function step(actor::ActorICsimple) @@ -239,18 +248,19 @@ end #= == =# # LH # #= == =# -Base.@kwdef mutable struct ActorLHsimple <: HCDAbstractActor +mutable struct ActorLHsimple <: HCDAbstractActor dd::IMAS.dd - width::Vector{Real} - rho_0::Vector{Real} - current_efficiency::Vector{Real} + par::ParametersActor + width::AbstractVector{<:Real} + rho_0::AbstractVector{<:Real} + current_efficiency::AbstractVector{<:Real} end function ParametersActor(::Type{Val{:ActorLHsimple}}) par = ParametersActor(nothing) - par.width = Entry(Real, "", "Width of the deposition profile"; default=0.15) - par.rho_0 = Entry(Real, "", "Radial location of the deposition profile"; default=0.6) - par.current_efficiency = Entry(Real, "A/W", "Current drive efficiency"; default=0.4) + par.width = Entry(Union{Real,AbstractVector{<:Real}}, "", "Width of the deposition profile"; default=0.15) + par.rho_0 = Entry(Union{Real,AbstractVector{<:Real}}, "", "Radial location of the deposition profile"; default=0.6) + par.current_efficiency = Entry(Union{Real,AbstractVector{<:Real}}, "A/W", "Current drive efficiency"; default=0.4) return par end @@ -264,15 +274,17 @@ This actor estimates the Lower-hybrid electron energy deposition and current dri """ function ActorLHsimple(dd::IMAS.dd, act::ParametersAllActors; kw...) par = act.ActorLHsimple(kw...) - actor = ActorLHsimple(dd; par.width, par.rho_0, par.current_efficiency) + actor = ActorLHsimple(dd, par) step(actor) finalize(actor) return actor end -function ActorLHsimple(dd::IMAS.dd; width::Real=0.15, rho_0::Real=0.6, current_efficiency::Real=0.4) - n_antennas = ones(length(dd.lh_antennas.antenna)) - return ActorICsimple(dd, n_antennas .* width, n_antennas .* rho_0, n_antennas .* current_efficiency) +function ActorLHsimple(dd::IMAS.dd, par::ParametersActor; kw...) + par = par(kw...) + n_antennas = length(dd.lh_antennas.antenna) + _, width, rho_0, current_efficiency = same_length_vectors(1:n_antennas, par.width, par.rho_0, par.current_efficiency) + return ActorLHsimple(dd, par, width, rho_0, current_efficiency) end function step(actor::ActorLHsimple) diff --git a/src/actors/stresses_actors.jl b/src/actors/stresses_actors.jl index 36cbc8f58..f2d343f18 100644 --- a/src/actors/stresses_actors.jl +++ b/src/actors/stresses_actors.jl @@ -1,8 +1,13 @@ #= ============== =# # OH TF stresses # #= ============== =# -Base.@kwdef mutable struct ActorStresses <: ReactorAbstractActor +mutable struct ActorStresses <: ReactorAbstractActor dd::IMAS.dd + par::ParametersActor + function ActorStresses(dd::IMAS.dd, par::ParametersActor; kw...) + par = par(kw...) + return new(dd, par) + end end function ParametersActor(::Type{Val{:ActorStresses}}) @@ -15,15 +20,14 @@ end """ ActorStresses(dd::IMAS.dd, act::ParametersAllActors; kw...) -This actor estimates vertical field from PF coils and its contribution to flux swing, where -`eqt` is supposed to be the equilibrium right at the end of the rampup phase, beginning of flattop +This actor estimates vertical field from PF coils and its contribution to flux swing !!! note Stores data in `dd.solid_mechanics` """ function ActorStresses(dd::IMAS.dd, act::ParametersAllActors; kw...) par = act.ActorStresses(kw...) - actor = ActorStresses(dd) + actor = ActorStresses(dd, par) step(actor; par.n_points) finalize(actor) if par.do_plot @@ -69,7 +73,7 @@ function step(actor::ActorStresses; n_points::Integer=5) f_struct_oh=f_struct_oh, f_struct_pl=1.0, n_points=n_points, - verbose=false, + verbose=false ) end diff --git a/src/actors/transport_actors.jl b/src/actors/transport_actors.jl index 1af7cee82..aab48ecab 100644 --- a/src/actors/transport_actors.jl +++ b/src/actors/transport_actors.jl @@ -4,8 +4,9 @@ import TAUENN # TAUENN actor # #= ================ =# -Base.@kwdef mutable struct ActorTauenn <: PlasmaAbstractActor +mutable struct ActorTauenn <: PlasmaAbstractActor dd::IMAS.dd + par::ParametersActor tauenn_parameters::TAUENN.TauennParameters tauenn_outputs::TAUENN.TauennOutputs end @@ -28,7 +29,7 @@ end """ ActorTauenn(dd::IMAS.dd, act::ParametersAllActors; kw...) -This actor estimates the core-transport using Tauenn and evolves the kinetic profiles according to heat and particle flux matching. +This actor estimates the core-transport using Tauenn, which evolves the kinetic profiles according to heat and particle flux matching. The pedestal in this actor is evolved using EPED-NN. @@ -37,11 +38,15 @@ The pedestal in this actor is evolved using EPED-NN. """ function ActorTauenn(dd::IMAS.dd, act::ParametersAllActors; kw...) par = act.ActorTauenn(kw...) - if par.do_plot - ps = plot(dd.core_sources; color=:gray) - pp = plot(dd.core_profiles; color=:gray, label="") - end - actor = ActorTauenn(dd; + actor = ActorTauenn(dd, par) + step(actor) + finalize(actor) + return actor +end + +function ActorTauenn(dd::IMAS.dd, par::ParametersActor; kw...) + par = par(kw...) + tauenn_parameters = TAUENN.TauennParameters(; par.error, par.eped_factor, par.rho_fluxmatch, @@ -50,8 +55,17 @@ function ActorTauenn(dd::IMAS.dd, act::ParametersAllActors; kw...) par.transport_model, par.confinement_factor, par.warn_nn_train_bounds) - step(actor; verbose=par.verbose) - finalize(actor) + return ActorTauenn(dd, par, tauenn_parameters, TAUENN.TauennOutputs()) +end + +function step(actor::ActorTauenn) + dd = actor.dd + par = actor.par + if par.do_plot + ps = plot(dd.core_sources; color=:gray) + pp = plot(dd.core_profiles; color=:gray, label="") + end + actor.tauenn_outputs = TAUENN.tau_enn(dd, actor.tauenn_parameters; par.verbose) if par.do_plot display(plot!(ps, dd.core_sources)) display(plot!(pp, dd.core_profiles)) @@ -60,17 +74,4 @@ function ActorTauenn(dd::IMAS.dd, act::ParametersAllActors; kw...) display(actor.tauenn_parameters) end return actor -end - -function ActorTauenn(dd::IMAS.dd; kw...) - tauenn_parameters = TAUENN.TauennParameters() - for key in keys(kw) - setfield!(tauenn_parameters, key, kw[key]) - end - return ActorTauenn(dd, tauenn_parameters, TAUENN.TauennOutputs()) -end - -function step(actor::ActorTauenn; verbose=false) - actor.tauenn_outputs = TAUENN.tau_enn(actor.dd, actor.tauenn_parameters; verbose) - return actor end \ No newline at end of file diff --git a/src/gasc.jl b/src/ddinit/gasc.jl similarity index 59% rename from src/gasc.jl rename to src/ddinit/gasc.jl index f124766c4..7e6509f90 100644 --- a/src/gasc.jl +++ b/src/ddinit/gasc.jl @@ -1,13 +1,27 @@ import JSON import PeriodicTable: elements -struct GASC +mutable struct GASC filename::String + data::Dict case::Int - solution::Dict version::Int end +function Base.getproperty(gasc::GASC, field::Symbol) + if field == :solution + return getfield(gasc, :data)["SOLUTIONS"][gasc.case+1] + elseif field == :outputs + return getfield(gasc, :data)["SOLUTIONS"][gasc.case+1]["OUTPUTS"] + elseif field == :inputs + return getfield(gasc, :data)["SOLUTIONS"][gasc.case+1]["INPUTS"] + elseif field == :constraints + return getfield(gasc, :data)["SETUP"]["SETTINGS"]["constraints"] + else + return getfield(gasc, field) + end +end + """ GASC(filename::String, case::Int) @@ -24,12 +38,11 @@ function GASC(filename::String, case::Int) version = 1 end # GASC struct - gasc = GASC(filename, case, data["SOLUTIONS"][case_j], version) - # convert list of floats to arrays - for item in keys(gasc.solution["OUTPUTS"]["numerical profiles"]) - if item in ["boundaryInnerTF", "boundaryOuterTF"] - else - gasc.solution["OUTPUTS"]["numerical profiles"][item] = Vector{Float64}(gasc.solution["OUTPUTS"]["numerical profiles"][item]) + gasc = GASC(filename, data, case, version) + # convert list of floats to arrays + for item in keys(gasc.outputs["numerical profiles"]) + if item ∉ ["boundaryInnerTF", "boundaryOuterTF"] + gasc.outputs["numerical profiles"][item] = Vector{Float64}(gasc.outputs["numerical profiles"][item]) end end return gasc @@ -46,11 +59,13 @@ function case_parameters(gasc::GASC) ini.gasc.filename = gasc.filename ini.gasc.case = gasc.case - ini.general.casename = "GASC" + ini.general.casename = gasc.data["LABEL"] ini.general.init_from = :scalars gasc_2_build(gasc, ini, act) + gasc_2_target(gasc, ini, act) + gasc_2_equilibrium(gasc, ini, act) gasc_2_sources(gasc, ini, act) @@ -66,19 +81,17 @@ end Convert core_profiles information in GASC solution to FUSE `ini` and `act` parameters """ function gasc_2_core_profiles(gasc::GASC, ini::ParametersAllInits, act::ParametersAllActors) - gascsol = gasc.solution - - ini.core_profiles.ne_ped = gascsol["OUTPUTS"]["plasma parameters"]["neped"] * 1e20 - ini.core_profiles.greenwald_fraction = gascsol["OUTPUTS"]["plasma parameters"]["greenwaldFraction"] + ini.core_profiles.ne_ped = gasc.outputs["plasma parameters"]["neped"] * 1e20 + ini.core_profiles.greenwald_fraction = gasc.outputs["plasma parameters"]["greenwaldFraction"] ini.core_profiles.T_shaping = 1.8 - i_ped = argmin(abs.(gascsol["OUTPUTS"]["numerical profiles"]["neProf"] .- gascsol["OUTPUTS"]["plasma parameters"]["neped"] / gascsol["OUTPUTS"]["plasma parameters"]["ne0"])) - ini.core_profiles.w_ped = 1 - gascsol["OUTPUTS"]["numerical profiles"]["rProf"][i_ped] - ini.core_profiles.zeff = gascsol["OUTPUTS"]["impurities"]["effectiveZ"] + i_ped = argmin(abs.(gasc.outputs["numerical profiles"]["neProf"] .- gasc.outputs["plasma parameters"]["neped"] / gasc.outputs["plasma parameters"]["ne0"])) + ini.core_profiles.w_ped = 1 - gasc.outputs["numerical profiles"]["rProf"][i_ped] + ini.core_profiles.zeff = gasc.outputs["impurities"]["effectiveZ"] ini.core_profiles.rot_core = 0.0 # Not in GASC ini.core_profiles.bulk = :DT - ini.core_profiles.helium_fraction = gascsol["INPUTS"]["impurities"]["heliumFraction"] - ini.core_profiles.impurity = Symbol(elements[Int(gascsol["INPUTS"]["impurities"]["impurityZ"])].symbol) - ini.core_profiles.ejima = gascsol["INPUTS"]["plasma parameters"]["ejimaCoeff"] + ini.core_profiles.helium_fraction = gasc.inputs["impurities"]["heliumFraction"] + ini.core_profiles.impurity = Symbol(elements[Int(gasc.inputs["impurities"]["impurityZ"])].symbol) + ini.core_profiles.ejima = gasc.inputs["plasma parameters"]["ejimaCoeff"] return ini end @@ -88,17 +101,22 @@ end Convert equilibrium information in GASC solution to FUSE `ini` and `act` parameters """ function gasc_2_equilibrium(gasc::GASC, ini::ParametersAllInits, act::ParametersAllActors) - gascsol = gasc.solution - ini.equilibrium.B0 = gascsol["INPUTS"]["conductors"]["magneticFieldOnAxis"] - ini.equilibrium.R0 = gascsol["INPUTS"]["radial build"]["majorRadius"] + ini.equilibrium.B0 = gasc.inputs["conductors"]["magneticFieldOnAxis"] + ini.equilibrium.R0 = gasc.inputs["radial build"]["majorRadius"] ini.equilibrium.Z0 = 0.0 - ini.equilibrium.ϵ = 1 / gascsol["INPUTS"]["radial build"]["aspectRatio"] - ini.equilibrium.κ = gascsol["OUTPUTS"]["plasma parameters"]["elongation"] - ini.equilibrium.δ = gascsol["INPUTS"]["plasma parameters"]["triangularity"] - ini.equilibrium.βn = gascsol["OUTPUTS"]["plasma parameters"]["betaN"] - ini.equilibrium.ip = gascsol["INPUTS"]["plasma parameters"]["plasmaCurrent"] * 1E6 - ini.equilibrium.x_point = gascsol["INPUTS"]["divertor metrics"]["numberDivertors"] > 0 - ini.equilibrium.symmetric = (mod(gascsol["INPUTS"]["divertor metrics"]["numberDivertors"], 2) == 0) + ini.equilibrium.ϵ = 1 / gasc.inputs["radial build"]["aspectRatio"] + ini.equilibrium.κ = gasc.outputs["plasma parameters"]["elongation"] + ini.equilibrium.δ = gasc.inputs["plasma parameters"]["triangularity"] + + Pavg = gasc.outputs["plasma parameters"]["pressureVolAvg"] + V = gasc.outputs["plasma parameters"]["plasmaVolume"] + vol = gasc.outputs["numerical profiles"]["volumeProf"] .* V + P1 = sum(IMAS.gradient(vol) .* LinRange(1.0, 0.0, length(vol))) / V + ini.equilibrium.pressure_core = Pavg / P1 + + ini.equilibrium.ip = gasc.inputs["plasma parameters"]["plasmaCurrent"] * 1E6 + ini.equilibrium.x_point = gasc.inputs["divertor metrics"]["numberDivertors"] > 0 + ini.equilibrium.symmetric = (mod(gasc.inputs["divertor metrics"]["numberDivertors"], 2) == 0) return ini end @@ -108,10 +126,8 @@ end Convert sources (NBI, EC, IC, LH) information in GASC solution to FUSE `ini` and `act` parameters """ function gasc_2_sources(gasc::GASC, ini::ParametersAllInits, act::ParametersAllActors) - gascsol = gasc.solution - - inputs = gascsol["INPUTS"]["current drive"] - outputs = gascsol["OUTPUTS"]["current drive"] + inputs = gasc.inputs["current drive"] + outputs = gasc.outputs["current drive"] cd_powers = Float64[] ini.nbi.power_launched = Float64[] @@ -185,29 +201,37 @@ end Convert radial build information in GASC solution to FUSE `ini` and `act` parameters """ function gasc_2_build(gasc::GASC, ini::ParametersAllInits, act::ParametersAllActors) - gascsol = gasc.solution - ini.build.layers = gasc_2_layers(gascsol) - ini.build.symmetric = (mod(gascsol["INPUTS"]["divertor metrics"]["numberDivertors"], 2) == 0) + ini.build.layers = gasc_2_layers(gasc) + ini.build.symmetric = (mod(gasc.inputs["divertor metrics"]["numberDivertors"], 2) == 0) ini.tf.technology = gasc_2_coil_technology(gasc, :TF) ini.oh.technology = gasc_2_coil_technology(gasc, :OH) ini.pf_active.technology = gasc_2_coil_technology(gasc, :PF) - ini.center_stack.bucked = gascsol["INPUTS"]["radial build"]["isBucked"] - ini.center_stack.noslip = gascsol["INPUTS"]["radial build"]["nonSlip"] - ini.center_stack.plug = gascsol["INPUTS"]["radial build"]["hasPlug"] + ini.center_stack.bucked = gasc.inputs["radial build"]["isBucked"] + ini.center_stack.noslip = gasc.inputs["radial build"]["nonSlip"] + ini.center_stack.plug = gasc.inputs["radial build"]["hasPlug"] + return ini +end + +""" + gasc_2_target(gasc::GASC, ini::ParametersAllInits, act::ParametersAllActors) - ini.oh.flattop_duration = gascsol["INPUTS"]["plasma parameters"]["flattopDuration"] +Convert nominal target design information in GASC solution to FUSE `ini` and `act` parameters +""" +function gasc_2_target(gasc::GASC, ini::ParametersAllInits, act::ParametersAllActors) + ini.target.flattop_duration = gasc.inputs["plasma parameters"]["flattopDuration"] + ini.target.power_electric_net = gasc.constraints["lowerOutputConstraints"]["powerNet"] * 1E6 return ini end """ - function gasc_2_layers(gascsol::Dict) + function gasc_2_layers(gasc::GASC) Convert GASC ["OUTPUTS"]["radial build"] to FUSE build layers dictionary """ -function gasc_2_layers(gascsol::Dict) - gascrb = gascsol["OUTPUTS"]["radial build"] +function gasc_2_layers(gasc::GASC) + gascrb = gasc.outputs["radial build"] layers = DataStructures.OrderedDict() mapper = Dict( @@ -252,7 +276,7 @@ function gasc_2_layers(gascsol::Dict) d = gascrb[replace(g2, "InnerTF" => "TF")] - gascrb[replace(g1, "InnerTF" => "TF")] f1 = mapper[replace(replace(replace(replace(g1, "Ri" => ""), "Ro" => ""), "Inner" => ""), "Outer" => "")] f2 = mapper[replace(replace(replace(replace(g2, "Ri" => ""), "Ro" => ""), "Inner" => ""), "Outer" => "")] - if startswith(g1, "Ri") + if startswith(g1, "Ri") # hfs if contains(g1, "Inner") f1 = "hfs_" * f1 elseif contains(g1, "Outer") @@ -261,7 +285,8 @@ function gasc_2_layers(gascsol::Dict) f = f1 f = replace(f, r".fs_gap_TF_OH" => "gap_TF_OH") layers[f] = d - elseif startswith(g1, "Ro") && d > 0 + + elseif startswith(g1, "Ro") && (d > 0) # lfs if contains(g2, "Outer") f = "lfs_gap_$(f1)_$(f2)" else @@ -281,15 +306,46 @@ function gasc_2_layers(gascsol::Dict) end end - if false # to remove gap that prevents bucking - for k in collect(keys(layers)) - if k == "gap_TF_OH" - layers["OH"] += layers["gap_TF_OH"] - delete!(layers, k) - end + return layers +end + +""" + gasc_buck_OH_TF!(layers::DataStructures.OrderedDict) + +Remove gap between OH and TF to allow bucking (gap gets added to OH thickness) +""" +function gasc_buck_OH_TF!(layers::DataStructures.OrderedDict) + for k in collect(keys(layers)) + if k == "gap_TF_OH" + layers["OH"] += layers["gap_TF_OH"] + delete!(layers, k) end end + return layers +end + +""" + gasc_add_wall_layers!(layers::DataStructures.OrderedDict, thickness::Float64) +Add wall layer of given thickness expressed [meters] (gets subtracted from blanket layer) +""" +function gasc_add_wall_layers!(layers::DataStructures.OrderedDict; thickness::Float64) + tmp = DataStructures.OrderedDict() + for layer in keys(layers) + if layer == "hfs_blanket" + tmp[layer] = layers[layer] - thickness + tmp["hfs_first_wall"] = thickness + elseif layer == "lfs_blanket" + tmp["lfs_first_wall"] = thickness + tmp[layer] = layers[layer] - thickness + else + tmp[layer] = layers[layer] + end + end + empty!(layers) + for layer in keys(tmp) + layers[layer] = tmp[layer] + end return layers end @@ -299,35 +355,34 @@ end Convert coil technology information in GASC solution to FUSE `coil_technology` data structure """ function gasc_2_coil_technology(gasc::GASC, coil_type::Symbol) - gascsol = gasc.solution if coil_type ∉ [:OH, :TF, :PF] error("Supported coil type are [:OH, :TF, :PF]") end - if gascsol["INPUTS"]["conductors"]["superConducting"] == "copper" + if gasc.inputs["conductors"]["superConducting"] == "copper" coil_tech = coil_technology(:copper) else - if gascsol["INPUTS"]["conductors"]["superConducting"] == "LTS" + if gasc.inputs["conductors"]["superConducting"] == "LTS" coil_tech = coil_technology(:LTS) - elseif gascsol["INPUTS"]["conductors"]["superConducting"] == "HTS" + elseif gasc.inputs["conductors"]["superConducting"] == "HTS" coil_tech = coil_technology(:HTS) end if coil_type == :PF # assume PF coils are always LTS coil_tech = coil_technology(:LTS) end - if "thermalStrain$coil_type" ∉ keys(gascsol["INPUTS"]["conductors"]) + if "thermalStrain$coil_type" ∉ keys(gasc.inputs["conductors"]) coil_tech.thermal_strain = 0.0 coil_tech.JxB_strain = 0.0 else - coil_tech.thermal_strain = gascsol["INPUTS"]["conductors"]["thermalStrain$coil_type"] - coil_tech.JxB_strain = gascsol["INPUTS"]["conductors"]["structuralStrain$coil_type"] + coil_tech.thermal_strain = gasc.inputs["conductors"]["thermalStrain$coil_type"] + coil_tech.JxB_strain = gasc.inputs["conductors"]["structuralStrain$coil_type"] end end - if "fractionVoid$coil_type" ∉ keys(gascsol["INPUTS"]["conductors"]) - coil_tech.fraction_void = gascsol["INPUTS"]["conductors"]["fractionVoidOH"] - coil_tech.fraction_stainless = gascsol["INPUTS"]["conductors"]["fractionStainlessOH"] + if "fractionVoid$coil_type" ∉ keys(gasc.inputs["conductors"]) + coil_tech.fraction_void = gasc.inputs["conductors"]["fractionVoidOH"] + coil_tech.fraction_stainless = gasc.inputs["conductors"]["fractionStainlessOH"] else - coil_tech.fraction_void = gascsol["INPUTS"]["conductors"]["fractionVoid$coil_type"] - coil_tech.fraction_stainless = gascsol["INPUTS"]["conductors"]["fractionStainless$coil_type"] + coil_tech.fraction_void = gasc.inputs["conductors"]["fractionVoid$coil_type"] + coil_tech.fraction_stainless = gasc.inputs["conductors"]["fractionStainless$coil_type"] end return set_new_base!(coil_tech) end @@ -337,12 +392,12 @@ function compare(dd::IMAS.dd, gasc::GASC) # collisionless bootstrap coefficient FUSE = IMAS.collisionless_bootstrap_coefficient(dd) - GASC = gasc.solution["INPUTS"]["plasma parameters"]["user_bootstrapCoefficient"] + GASC = gasc.inputs["plasma parameters"]["user_bootstrapCoefficient"] df.Cbs = [FUSE, GASC] # fusion power [MW] FUSE = IMAS.fusion_power(dd.core_profiles.profiles_1d[]) / 1E6 - GASC = gasc.solution["OUTPUTS"]["power balance"]["powerFusion"] + GASC = gasc.outputs["power balance"]["powerFusion"] df.Pfusion = [FUSE, GASC] df diff --git a/src/ddinit/init.jl b/src/ddinit/init.jl index fff36a85a..6a37e6bb9 100644 --- a/src/ddinit/init.jl +++ b/src/ddinit/init.jl @@ -13,6 +13,7 @@ function init(dd::IMAS.dd, ini::ParametersAllInits, act::ParametersAllActors; do if ini.general.init_from == :ods ods_items = keys(IMAS.json2imas(ini.ods.filename)) end + # initialize equilibrium if !ismissing(ini.equilibrium, :B0) || :equilibrium ∈ ods_items init_equilibrium(dd, ini, act) diff --git a/src/ddinit/init_build.jl b/src/ddinit/init_build.jl index 26eaf4bd1..4125c9919 100644 --- a/src/ddinit/init_build.jl +++ b/src/ddinit/init_build.jl @@ -1,4 +1,5 @@ import LibGEOS +import GeoInterface import Interpolations import DataStructures @@ -57,7 +58,7 @@ Initialize `dd.build` starting from `ini` and `act` parameters """ function init_build(dd::IMAS.dd, ini::ParametersAllInits, act::ParametersAllActors) init_from = ini.general.init_from - + if init_from == :ods dd1 = IMAS.json2imas(ini.ods.filename) if length(keys(dd1.wall)) > 0 @@ -70,39 +71,67 @@ function init_build(dd::IMAS.dd, ini::ParametersAllInits, act::ParametersAllActo end end - if init_from == :scalars - if ismissing(ini.build, :layers) - init_build( - dd.build, - dd.equilibrium.time_slice[], - IMAS.first_wall(dd.wall); - shield=ini.build.shield, - blanket=ini.build.blanket, - vessel=ini.build.vessel, - pf_inside_tf=(ini.pf_active.n_pf_coils_inside > 0), - pf_outside_tf=(ini.pf_active.n_pf_coils_outside > 0)) + # if layers are not filled explicitly, then generate them from fractions in ini.build. + if ismissing(ini.build, :layers) + layers = layers_meters_from_fractions( + dd.equilibrium.time_slice[], + IMAS.first_wall(dd.wall); + shield=ini.build.shield, + blanket=ini.build.blanket, + vessel=ini.build.vessel, + plasma_gap=ini.build.plasma_gap, + pf_inside_tf=(ini.pf_active.n_pf_coils_inside > 0), + pf_outside_tf=(ini.pf_active.n_pf_coils_outside > 0)) + else + layers = deepcopy(ini.build.layers) + end + + # scale radial build layers based on equilibrium R0 and ϵ + iplasma = findfirst(key -> key in ["plasma", :plasma], collect(keys(layers))) + R_hfs_build = sum([d for (k, d) in enumerate(values(layers)) if k < iplasma]) + if ismissing(ini.equilibrium, :R0) + R0 = R_hfs_build + collect(values(layers))[iplasma] / 2.0 + else + R0 = ini.equilibrium.R0 + end + if ismissing(ini.equilibrium, :ϵ) + a_gap = collect(values(layers))[iplasma] / 2.0 + else + a_gap = ini.equilibrium.ϵ * R0 * (1.0 + ini.build.plasma_gap) + end + R_hfs = R0 - a_gap + for key in keys(layers) + if key in ["plasma", :plasma] + layers[key] = a_gap * 2.0 else - init_build(dd.build, ini.build.layers) + layers[key] = layers[key] * R_hfs / R_hfs_build end end + # populate dd.build with radial build layers + init_build(dd.build, layers) + # set the TF shape tf_to_plasma = IMAS.get_build(dd.build, fs=_hfs_, return_only_one=false, return_index=true) dd.build.layer[tf_to_plasma[1]].shape = Int(_offset_) dd.build.layer[tf_to_plasma[2]].shape = Int(to_enum(ini.tf.shape)) + # set all other shapes for k in tf_to_plasma[2:end] dd.build.layer[k+1].shape = Int(_convex_hull_) end - for k in tf_to_plasma[2:end-ini.build.n_first_wall_conformal_layers] - dd.build.layer[k+1].shape = Int(_negative_offset_) + k = tf_to_plasma[end] + if (dd.build.layer[k].type == Int(_wall_)) && ((dd.build.layer[k-1].type == Int(_blanket_)) || (dd.build.layer[k-1].type == Int(_shield_))) + dd.build.layer[k].shape = Int(_offset_) + end + if ini.build.n_first_wall_conformal_layers >= 0 + for k in tf_to_plasma[2:end-ini.build.n_first_wall_conformal_layers] + dd.build.layer[k+1].shape = Int(_negative_offset_) + end end # 2D build cross-section ActorCXbuild(dd, act) - # flattop duration - dd.build.oh.flattop_duration = ini.oh.flattop_duration - # TF coils dd.build.tf.coils_n = ini.tf.n_coils # set the toroidal thickness of the TF coils based on the innermost radius and the number of coils @@ -214,13 +243,13 @@ end Initialization of build IDS based on equilibrium time_slice """ -function init_build( - bd::IMAS.build, +function layers_meters_from_fractions( eqt::IMAS.equilibrium__time_slice, wall::T where {T<:Union{IMAS.wall__description_2d___limiter__unit___outline,Missing}}; blanket::Float64=1.0, shield::Float64=0.5, vessel::Float64=0.125, + plasma_gap::Float64=0.1, pf_inside_tf::Bool=false, pf_outside_tf::Bool=true) @@ -230,8 +259,7 @@ function init_build( else rmin = eqt.boundary.geometric_axis.r - eqt.boundary.minor_radius rmax = eqt.boundary.geometric_axis.r + eqt.boundary.minor_radius - gap = (rmax - rmin) / 20.0 # plasma-wall gap - gap = (rmax - rmin) / 20.0 # plasma-wall gap + gap = (rmax - rmin) / 2.0 * plasma_gap rmin -= gap rmax += gap end @@ -287,8 +315,5 @@ function init_build( end end - # radial build - init_build(bd; layers...) - - return bd + return layers end diff --git a/src/ddinit/init_core_profiles.jl b/src/ddinit/init_core_profiles.jl index 197aec906..6d54bfb62 100644 --- a/src/ddinit/init_core_profiles.jl +++ b/src/ddinit/init_core_profiles.jl @@ -24,6 +24,7 @@ function init_core_profiles(dd::IMAS.dd, ini::ParametersAllInits, act::Parameter dd.equilibrium, dd.summary; ne_ped=ini.core_profiles.ne_ped, + pressure_core=dd.equilibrium.time_slice[].profiles_1d.pressure[1], greenwald_fraction=ini.core_profiles.greenwald_fraction, helium_fraction=ini.core_profiles.helium_fraction, T_shaping=ini.core_profiles.T_shaping, @@ -48,6 +49,7 @@ function init_core_profiles( eq::IMAS.equilibrium, summary::IMAS.summary; ne_ped::Real, + pressure_core::Real, greenwald_fraction::Real, helium_fraction::Real, w_ped::Real, @@ -68,10 +70,10 @@ function init_core_profiles( cp1d.zeff = ones(ngrid) .* zeff cp1d.rotation_frequency_tor_sonic = rot_core .* (1.0 .- cp1d.grid.rho_tor_norm) - # Set ions - # DT == 1 - # Imp == 2 - # He == 3 + # Set ions: + # 1. DT + # 2. Imp + # 3. He ion = resize!(cp1d.ion, "label" => String(bulk)) fill!(ion, IMAS.ion_element(ion_symbol=bulk)) @assert ion.element[1].z_n == 1 "Bulk ion must be a Hydrogen isotope [:H, :D, :DT, :T]" @@ -107,18 +109,16 @@ function init_core_profiles( niFraction[1] = (zimp - zeff + 4 * niFraction[3] - 2 * zimp * niFraction[3]) / (zimp - 1) niFraction[2] = (zeff - niFraction[1] - 4 * niFraction[3]) / zimp^2 @assert !any(niFraction .< 0.0) "zeff impossible to match for given helium fraction [$helium_fraction] and zeff [$zeff]" + ni_core = 0.0 for i = 1:length(cp1d.ion) cp1d.ion[i].density_thermal = cp1d.electrons.density_thermal .* niFraction[i] + ni_core += cp1d.electrons.density_thermal[1] * niFraction[i] end # Set temperatures - eqt = eq.time_slice[] - betaN = eqt.global_quantities.beta_normal - Bt = @ddtime eq.vacuum_toroidal_field.b0 - Ip = eqt.global_quantities.ip - a = eqt.boundary.minor_radius - Te_core = 10.0 * betaN * abs(Bt * (Ip / 1e6)) / a / (ne_core / 1e20) / (2.0 * 1.6e1 * 4.0 * pi * 1.0e-4) - Te_ped = Te_core / 4 + Te_core = pressure_core / (ni_core + ne_core) / IMAS.constants.e + Te_ped = sqrt(Te_core / 1000.0 / 3.0) * 1000.0 + cp1d.electrons.temperature = Hmode_profiles(80.0, Te_ped, Te_core, ngrid, T_shaping, T_shaping, w_ped) for i = 1:length(cp1d.ion) cp1d.ion[i].temperature = cp1d.electrons.temperature ./ T_ratio @@ -126,7 +126,7 @@ function init_core_profiles( # remove He if not present if sum(niFraction[3]) == 0.0 - deleteat!(cp1d.ion,3) + deleteat!(cp1d.ion, 3) end # ejima diff --git a/src/ddinit/init_core_sources.jl b/src/ddinit/init_core_sources.jl index c036486b6..4f501435f 100644 --- a/src/ddinit/init_core_sources.jl +++ b/src/ddinit/init_core_sources.jl @@ -18,12 +18,12 @@ end function init_nbi( dd::IMAS.dd, - power_launched::Union{Real,Vector}, - beam_energy::Union{Real,Vector}, - beam_mass::Union{Real,Vector}, - toroidal_angle::Union{Real,Vector}, - efficiency_conversion::Union{Missing,Real,Vector}, - efficiency_transmission::Union{Missing,Real,Vector}) + power_launched::Union{Real,AbstractVector{<:Real}}, + beam_energy::Union{Real,AbstractVector{<:Real}}, + beam_mass::Union{Real,AbstractVector{<:Real}}, + toroidal_angle::Union{Real,AbstractVector{<:Real}}, + efficiency_conversion::Union{Missing,Real,AbstractVector{<:Real}}, + efficiency_transmission::Union{Missing,Real,AbstractVector{<:Real}}) power_launched, beam_energy, beam_mass, toroidal_angle, efficiency_conversion, efficiency_transmission = same_length_vectors(power_launched, beam_energy, beam_mass, toroidal_angle, efficiency_conversion, efficiency_transmission) @@ -60,13 +60,13 @@ end function init_ec_launchers( dd::IMAS.dd, - power_launched::Union{Real,Vector}, - efficiency_conversion::Union{Real,Vector,Missing}, - efficiency_transmission::Union{Real,Vector,Missing}) + power_launched::Union{Real,AbstractVector{<:Real}}, + efficiency_conversion::Union{Real,AbstractVector{<:Real},Missing}, + efficiency_transmission::Union{Real,AbstractVector{<:Real},Missing}) (power_launched, efficiency_conversion, efficiency_transmission) = same_length_vectors(power_launched, efficiency_conversion, efficiency_transmission) for idx in 1:length(power_launched) - ecl = resize!(dd.ec_launchers.launcher, idx)[idx] + ecl = resize!(dd.ec_launchers.beam, idx)[idx] ecl.name = length(power_launched) > 1 ? "ec_$idx" : "ec" @ddtime(ecl.power_launched.data = power_launched[idx]) ecl.available_launch_power = power_launched[idx] @@ -91,10 +91,10 @@ end function init_ic_antennas( dd::IMAS.dd, - power_launched::Union{Real,Vector}, - efficiency_conversion::Union{Real,Vector,Missing}, - efficiency_transmission::Union{Real,Vector,Missing}, - efficiency_coupling::Union{Real,Vector,Missing}) + power_launched::Union{Real,AbstractVector{<:Real}}, + efficiency_conversion::Union{Real,AbstractVector{<:Real},Missing}, + efficiency_transmission::Union{Real,AbstractVector{<:Real},Missing}, + efficiency_coupling::Union{Real,AbstractVector{<:Real},Missing}) (power_launched, efficiency_conversion, efficiency_transmission, efficiency_coupling) = same_length_vectors(power_launched, efficiency_conversion, efficiency_transmission, efficiency_coupling) for idx in 1:length(power_launched) @@ -124,10 +124,10 @@ end function init_lh_antennas( dd::IMAS.dd, - power_launched::Union{Real,Vector}, - efficiency_conversion::Union{Real,Vector,Missing}, - efficiency_transmission::Union{Real,Vector,Missing}, - efficiency_coupling::Union{Real,Vector,Missing}) + power_launched::Union{Real,AbstractVector{<:Real}}, + efficiency_conversion::Union{Real,AbstractVector{<:Real},Missing}, + efficiency_transmission::Union{Real,AbstractVector{<:Real},Missing}, + efficiency_coupling::Union{Real,AbstractVector{<:Real},Missing}) (power_launched, efficiency_conversion, efficiency_transmission, efficiency_coupling) = same_length_vectors(power_launched, efficiency_conversion, efficiency_transmission, efficiency_coupling) for idx in 1:length(power_launched) lha = resize!(dd.lh_antennas.antenna, idx)[idx] diff --git a/src/ddinit/init_equilibrium.jl b/src/ddinit/init_equilibrium.jl index 5f5f6dfed..c8677bb0e 100644 --- a/src/ddinit/init_equilibrium.jl +++ b/src/ddinit/init_equilibrium.jl @@ -31,7 +31,7 @@ function init_equilibrium(dd::IMAS.dd, ini::ParametersAllInits, act::ParametersA κ=ini.equilibrium.κ, δ=ini.equilibrium.δ, ζ=ini.equilibrium.ζ, - βn=ini.equilibrium.βn, + pressure_core=ini.equilibrium.pressure_core, ip=ini.equilibrium.ip, boundary_switch=ini.equilibrium.boundary_from, MXH_params=getproperty(ini.equilibrium, :MXH_params, missing), @@ -39,7 +39,10 @@ function init_equilibrium(dd::IMAS.dd, ini::ParametersAllInits, act::ParametersA symmetric=ini.equilibrium.symmetric) # solve equilibrium - ActorEquilibrium(dd, act) + act_copy = deepcopy(act) + act_copy.ActorCHEASE.rescale_eq_to_ip = true + ActorEquilibrium(dd, act_copy) + end # field null surface @@ -54,7 +57,7 @@ function init_equilibrium(dd::IMAS.dd, ini::ParametersAllInits, act::ParametersA end """ - function init_equilibrium( + init_equilibrium( eq::IMAS.equilibrium; B0::Real, R0::Real, @@ -63,9 +66,12 @@ end κ::Real, δ::Real, ζ::Real, - βn::Real, + pressure_core::Real, ip::Real, - x_point::Union{Vector,NTuple{2},Bool} = false, + boundary_switch::Symbol, + rz_points::Union{Missing,Vector{Vector{<:Real}}}=missing, + MXH_params::Union{Missing,Vector{<:Real}}=missing, + x_point::Union{AbstractVector,NTuple{2},Bool}=false, symmetric::Bool=true) Initialize equilibrium IDS based on some basic Miller geometry parameters @@ -79,7 +85,7 @@ function init_equilibrium( κ::Real, δ::Real, ζ::Real, - βn::Real, + pressure_core::Real, ip::Real, boundary_switch::Symbol, rz_points::Union{Missing,Vector{Vector{<:Real}}}=missing, @@ -115,34 +121,32 @@ function init_equilibrium( # scalar quantities eqt.global_quantities.ip = ip - eqt.global_quantities.beta_normal = βn eq.vacuum_toroidal_field.r0 = R0 @ddtime eq.vacuum_toroidal_field.b0 = B0 # initial guesses for pressure and j_tor eq1d = eqt.profiles_1d psin = eq1d.psi = LinRange(0, 1, 129) - p_core_estimate = 6.0 * IMAS.pressure_avg_from_beta_n(βn, minor_radius, B0, ip) eq1d.j_tor = eqt.global_quantities.ip .* (1.0 .- psin .^ 2) ./ eqt.boundary.geometric_axis.r - eq1d.pressure = p_core_estimate .- p_core_estimate .* psin + eq1d.pressure = pressure_core .* (1.0 .- psin) - # R,Z boundary from: points if boundary_switch == :rz_points + # R,Z boundary from: points if ismissing(rz_points) error("ini.equilibrium.boundary_from is set as $boundary_switch but rz_points wasn't set") end eqt.boundary.outline.r, eqt.boundary.outline.z = rz_points[1], rz_points[2] - # R,Z boundary from: MXH elseif boundary_switch == :MXH_params + # R,Z boundary from: MXH if ismissing(MXH_params) error("ini.equilibrium.boundary_from is set as $boundary_switch but MXH_params wasn't set") end mxh = IMAS.MXH(MXH_params)() eqt.boundary.outline.r, eqt.boundary.outline.z = mxh[1], mxh[2] - # R,Z boundary from: scalars elseif boundary_switch == :scalars + # R,Z boundary from: scalars eqt.boundary.outline.r, eqt.boundary.outline.z = square_miller(R0, ϵ, κ, δ, ζ; exact=true, x_points=x_point !== false) eqt.boundary.outline.z .+= Z0 end diff --git a/src/ddinit/init_others.jl b/src/ddinit/init_others.jl index 182818e17..4f2ff35c2 100644 --- a/src/ddinit/init_others.jl +++ b/src/ddinit/init_others.jl @@ -19,5 +19,13 @@ function init_missing_from_ods(dd::IMAS.dd, ini::ParametersAllInits, act::Parame end end + # target + for field in [:power_electric_net, :flattop_duration, :tritium_breeding_ratio, :cost] + value = getproperty(ini.target, field, missing) + if value !== missing + setproperty!(dd.target, field, value) + end + end + return dd end \ No newline at end of file diff --git a/src/ddinit/init_pf_active.jl b/src/ddinit/init_pf_active.jl index db63b0e02..050dbb127 100644 --- a/src/ddinit/init_pf_active.jl +++ b/src/ddinit/init_pf_active.jl @@ -125,7 +125,7 @@ function init_pf_active( for k in lfs_out_indexes layer = bd.layer[k] - if (k == gap_cryostat_index) && (length(n_coils) >= krail+1) && (n_coils[krail+1] > 0) + if (k == gap_cryostat_index) && (length(n_coils) >= krail + 1) && (n_coils[krail+1] > 0) #pass elseif !contains(lowercase(layer.name), "coils") continue @@ -149,9 +149,9 @@ function init_pf_active( dcoil = (coil_size + coils_cleareance[krail]) / 2 * sqrt(2) inner_layer = IMAS.get_build(bd, identifier=bd.layer[k-1].identifier, fs=_hfs_) poly = LibGEOS.buffer(xy_polygon(inner_layer.outline.r, inner_layer.outline.z), dcoil) - rail_r = [v[1] for v in LibGEOS.coordinates(poly)[1]] - rail_z = [v[2] for v in LibGEOS.coordinates(poly)[1]] - rail_r, rail_z = IMAS.resample_2d_line(rail_r, rail_z, step=dr/3) + rail_r = [v[1] for v in GeoInterface.coordinates(poly)[1]] + rail_z = [v[2] for v in GeoInterface.coordinates(poly)[1]] + rail_r, rail_z = IMAS.resample_2d_line(rail_r, rail_z, step=dr / 3) # mark what regions on that rail do not intersect solid structures and can hold coils valid_k = [] diff --git a/src/logging.jl b/src/logging.jl index 3f9c94475..1e5fc12f0 100644 --- a/src/logging.jl +++ b/src/logging.jl @@ -4,7 +4,7 @@ function uwanted_warnings_filter(log_args) return !any([ contains(log_args.message, "both ImageMetadata and ImageAxes export"), startswith(log_args.message, "Keyword argument letter not supported with Plots.GRBackend") - ]) + ]) end logger = ActiveFilteredLogger(uwanted_warnings_filter, global_logger()) diff --git a/src/optimization.jl b/src/optimization.jl index c4a02cc81..269684486 100644 --- a/src/optimization.jl +++ b/src/optimization.jl @@ -48,8 +48,8 @@ const ObjectivesFunctionsLibrary = Dict{Symbol,ObjectiveFunction}() ObjectiveFunction(:min_cost, "\$M", dd -> dd.costing.cost, -Inf) ObjectiveFunction(:max_fusion, "MW", dd -> IMAS.fusion_power(dd.core_profiles.profiles_1d[]) / 1E6, Inf) ObjectiveFunction(:max_power_electric_net, "MW", dd -> @ddtime(dd.balance_of_plant.power_electric_net) / 1E6, Inf) -ObjectiveFunction(:max_flattop, "hours", dd -> dd.build.oh.flattop_estimate / 3600, Inf) -ObjectiveFunction(:max_log10_flattop, "log₁₀(hours)", dd -> log10(dd.build.oh.flattop_estimate / 3600), Inf) +ObjectiveFunction(:max_flattop, "hours", dd -> dd.build.oh.flattop_duration / 3600, Inf) +ObjectiveFunction(:max_log10_flattop, "log₁₀(hours)", dd -> log10(dd.build.oh.flattop_duration / 3600), Inf) function Base.show(io::IO, f::ObjectiveFunction) printstyled(io, f.name; bold=true, color=:blue) diff --git a/src/parameters.jl b/src/parameters.jl index fa842f0e5..92473780d 100644 --- a/src/parameters.jl +++ b/src/parameters.jl @@ -1,9 +1,7 @@ -using InteractiveUtils: subtypes import AbstractTrees abstract type AbstractParameter end abstract type AbstractParameters end -abstract type AbstractParametersActor end #= ===== =# # Entry # @@ -199,7 +197,7 @@ Generates all initalization parameters """ function ParametersAllInits() ini = ParametersAllInits(missing, WeakRef(missing), Dict{Symbol,Union{AbstractParameter,ParametersInit}}()) - for item in [:general, :equilibrium, :core_profiles, :pf_active, :oh, :tf, :center_stack, :nbi, :ec_launchers, :ic_antennas, :lh_antennas, :build, :gasc, :ods, :material] + for item in [:general, :equilibrium, :core_profiles, :pf_active, :oh, :tf, :center_stack, :nbi, :ec_launchers, :ic_antennas, :lh_antennas, :build, :gasc, :ods, :material, :target] setproperty!(ini, item, ParametersInit(item)) end ini._name = :ini @@ -254,7 +252,7 @@ function ParametersAllActors() setproperty!(act, par, ParametersActor(par)) catch e if typeof(e) <: InexistentParameterException - @warn e + @warn sprint(showerror, e) else rethrow() end @@ -288,6 +286,10 @@ function Base.keys(p::AbstractParameters) return keys(getfield(p, :_parameters)) end +function Base.values(p::AbstractParameters) + return values(getfield(p, :_parameters)) +end + function Base.getindex(p::AbstractParameters, field::Symbol) return getfield(p, :_parameters)[field] end @@ -398,7 +400,7 @@ function Base.iterate(par::AbstractParameters, state) return nothing end key = popfirst!(state) - data = par[key].value + data = par[key]#.value return key => data, state end @@ -406,24 +408,6 @@ function Base.show(io::IO, ::MIME"text/plain", pars::AbstractParameters, depth:: return AbstractTrees.print_tree(io, pars) end -function AbstractTrees.children(obj::AbstractDict) - return [obj[k] for k in sort(collect(keys(obj)))] -end - -function AbstractTrees.children(obj::Pair) - return [] -end - -function AbstractTrees.printnode(io::IO, obj::Pair) - printstyled(io, obj.first) - printstyled(io, " ➡ ") - if typeof(obj.second) <: AbstractFloat - printstyled(io, @sprintf("%3.3f", obj.second)) - else - printstyled(io, "$(repr(obj.second))") - end -end - function AbstractTrees.children(pars::AbstractParameters) return [pars[k] for k in sort(collect(keys(pars)))] end @@ -448,7 +432,7 @@ function AbstractTrees.printnode(io::IO, par::AbstractParameter) printstyled(io, par._name) printstyled(io, " ➡ ") printstyled(io, "$(repr(par.value))"; color=color) - if length(par.units) > 0 && par.value !== missing + if length(replace(par.units, "-" => "")) > 0 && par.value !== missing printstyled(io, " [$(par.units)]"; color=color) end end @@ -476,8 +460,8 @@ end This functor is used to override the parameters at function call """ function (par::AbstractParameters)(kw...) + par = deepcopy(par) if !isempty(kw) - par = deepcopy(par) for (key, value) in kw setproperty!(par, key, value) end @@ -518,19 +502,19 @@ Convert FUSE parameters to dictionary """ function par2dict(par::AbstractParameters) ret = Dict() - return par2dict(par, ret) + return par2dict!(par, ret) end -function par2dict(par::AbstractParameters, ret::AbstractDict) +function par2dict!(par::AbstractParameters, ret::AbstractDict) data = getfield(par, :_parameters) - return par2dict(data, ret) + return par2dict!(data, ret) end -function par2dict(data::AbstractDict, ret::AbstractDict) +function par2dict!(data::AbstractDict, ret::AbstractDict) for item in keys(data) if typeof(data[item]) <: AbstractParameters ret[item] = Dict() - par2dict(data[item], ret[item]) + par2dict!(data[item], ret[item]) elseif typeof(data[item]) <: AbstractParameter ret[item] = Dict() ret[item][:value] = data[item].value @@ -542,17 +526,76 @@ function par2dict(data::AbstractDict, ret::AbstractDict) end """ - par2json(@nospecialize(par::AbstractParameters), filename::String; kw...) + ini2json(ini::ParametersAllInits, filename::String; kw...) + +Save the FUSE parameters to a JSON file with give `filename` +`kw` arguments are passed to the JSON.print function +""" +function ini2json(ini::ParametersAllInits, filename::String; kw...) + return par2json(ini, filename; kw...) +end + +""" + act2json(act::ParametersAllActors, filename::String; kw...) Save the FUSE parameters to a JSON file with give `filename` `kw` arguments are passed to the JSON.print function """ +function act2json(act::ParametersAllActors, filename::String; kw...) + return par2json(act, filename; kw...) +end + function par2json(@nospecialize(par::AbstractParameters), filename::String; kw...) open(filename, "w") do io - JSON.print(io, par2dict(par); kw...) + JSON.print(io, par2dict(par), 1; kw...) end end +function dict2par!(dct::AbstractDict, par::AbstractParameters) + for (key, val) in par + if key ∈ keys(dct) + # this is if dct was par2dict function + dkey = key + dvalue = :value + else + # this is if dct was generated from json + dkey = string(key) + dvalue = "value" + end + if typeof(val) <: AbstractParameters + dict2par!(dct[dkey], val) + elseif dct[dkey][dvalue] === nothing + setproperty!(par, key, missing) + elseif typeof(dct[dkey][dvalue]) <: AbstractVector # this could be done more generally + setproperty!(par, key, Real[k for k in dct[dkey][dvalue]]) + else + try + setproperty!(par, key, Symbol(dct[dkey][dvalue])) + catch e + try + setproperty!(par, key, dct[dkey][dvalue]) + catch e + display((key, e)) + end + end + end + end + return par +end + +function json2par(filename::AbstractString, par_data::AbstractParameters) + json_data = JSON.parsefile(filename) + return dict2par!(json_data, par_data) +end + +function json2ini(filename::AbstractString) + return json2par(filename, ParametersAllInits()) +end + +function json2act(filename::AbstractString) + return json2par(filename, ParametersAllActors()) +end + #= ================= =# # Parameters errors # #= ================= =# @@ -560,7 +603,7 @@ struct InexistentParameterException <: Exception parameter_type::DataType path::Vector{Symbol} end -Base.showerror(io::IO, e::InexistentParameterException) = print(io, "ERROR: $(e.parameter_type) $(join(e.path,".")) does not exist") +Base.showerror(io::IO, e::InexistentParameterException) = print(io, "$(e.parameter_type).$(join(e.path,".")) does not exist") struct NotsetParameterException <: Exception path::Vector{Symbol} @@ -569,9 +612,9 @@ end NotsetParameterException(path::Vector{Symbol}) = NotsetParameterException(path, []) function Base.showerror(io::IO, e::NotsetParameterException) if length(e.options) > 0 - print(io, "ERROR: Parameter $(join(e.path,".")) is not set. Valid options are: $(join(map(repr,e.options),", "))") + print(io, "Parameter $(join(e.path,".")) is not set. Valid options are: $(join(map(repr,e.options),", "))") else - print(io, "ERROR: Parameter $(join(e.path,".")) is not set") + print(io, "Parameter $(join(e.path,".")) is not set") end end @@ -581,7 +624,7 @@ struct BadParameterException <: Exception options::Vector{Any} end Base.showerror(io::IO, e::BadParameterException) = - print(io, "ERROR: Parameter $(join(e.path,".")) = $(repr(e.value)) is not one of the valid options: $(join(map(repr,e.options),", "))") + print(io, "Parameter $(join(e.path,".")) = $(repr(e.value)) is not one of the valid options: $(join(map(repr,e.options),", "))") #= ============ =# # case studies # diff --git a/src/parameters_init.jl b/src/parameters_init.jl index 47c51b34e..8daf94603 100644 --- a/src/parameters_init.jl +++ b/src/parameters_init.jl @@ -22,21 +22,21 @@ end function ParametersInit(::Type{Val{:equilibrium}}) equilibrium = ParametersInit(nothing) equilibrium.B0 = Entry(Real, IMAS.equilibrium__vacuum_toroidal_field, :b0) - equilibrium.R0 = Entry(Real, IMAS.equilibrium__vacuum_toroidal_field, :r0) + equilibrium.R0 = Entry(Real, "m", "Geometric genter of the plasma. NOTE: This also scales the radial build layers.") equilibrium.Z0 = Entry(Real, "m", "Z offset of the machine midplane"; default=0.0) - equilibrium.ϵ = Entry(Real, "", "Plasma aspect ratio") + equilibrium.ϵ = Entry(Real, "", "Plasma inverse aspect ratio. NOTE: This also scales the radial build layers.") equilibrium.κ = Entry(Real, IMAS.equilibrium__time_slice___boundary, :elongation) equilibrium.δ = Entry(Real, IMAS.equilibrium__time_slice___boundary, :triangularity) equilibrium.ζ = Entry(Real, IMAS.equilibrium__time_slice___boundary, :squareness; default=0.0) - equilibrium.βn = Entry(Real, IMAS.equilibrium__time_slice___global_quantities, :beta_normal) + equilibrium.pressure_core = Entry(Real, "Pa", "On axis pressure") equilibrium.ip = Entry(Real, IMAS.equilibrium__time_slice___global_quantities, :ip) equilibrium.x_point = Entry(Union{NTuple{2},Bool}, IMAS.equilibrium__time_slice___boundary, :x_point) equilibrium.symmetric = Entry(Bool, "", "Is plasma up-down symmetric") equilibrium.ngrid = Entry(Int, "", "Resolution of the equilibrium grid"; default=129) - equilibrium.field_null_surface = Entry(Real, "", "ψn value of the field_null_surface. Disable with 0.0"; default=0.25)#, min=0.0, max=1.0) - equilibrium.boundary_from = Switch([:scalars, :MXH_params, :rz_points], "" ,"The starting r, z boundary taken from"; default=:scalars) + equilibrium.field_null_surface = Entry(Real, "", "ψn value of the field_null_surface. Disable with 0.0"; default=0.25) + equilibrium.boundary_from = Switch([:scalars, :MXH_params, :rz_points], "", "The starting r, z boundary taken from"; default=:scalars) equilibrium.MXH_params = Entry(Union{Nothing,Vector{<:Real}}, "", "Vector of MXH flats", default=missing) - equilibrium.rz_points = Entry(Union{Nothing, Vector{Vector{<:Real}}}, "m", "R_Z boundary as Vector{Vector{<:Real}}} : r = rz_points[1], z = rz_points[2]", default=missing) + equilibrium.rz_points = Entry(Union{Nothing,Vector{Vector{<:Real}}}, "m", "R_Z boundary as Vector{Vector{<:Real}}} : r = rz_points[1], z = rz_points[2]", default=missing) return equilibrium end @@ -78,7 +78,6 @@ end function ParametersInit(::Type{Val{:oh}}) oh = ParametersInit(nothing) oh.technology = ParametersInit(:coil_technology) - oh.flattop_duration = Entry(Real, "s", "Duration of the flattop (use Inf for steady-state)") return oh end @@ -90,7 +89,6 @@ function ParametersInit(::Type{Val{:center_stack}}) return center_stack end - function ParametersInit(::Type{Val{:nbi}}) nbi = ParametersInit(nothing) nbi.power_launched = Entry(Union{X,Vector{X}} where {X<:Real}, "W", "Beam power") @@ -105,8 +103,8 @@ end function ParametersInit(::Type{Val{:ec_launchers}}) ec_launchers = ParametersInit(nothing) ec_launchers.power_launched = Entry(Union{X,Vector{X}} where {X<:Real}, "W", "EC launched power") - ec_launchers.efficiency_conversion = Entry(Union{X,Vector{X}} where {X<:Real}, IMAS.ec_launchers__launcher___efficiency, :conversion) - ec_launchers.efficiency_transmission = Entry(Union{X,Vector{X}} where {X<:Real}, IMAS.ec_launchers__launcher___efficiency, :transmission) + ec_launchers.efficiency_conversion = Entry(Union{X,Vector{X}} where {X<:Real}, IMAS.ec_launchers__beam___efficiency, :conversion) + ec_launchers.efficiency_transmission = Entry(Union{X,Vector{X}} where {X<:Real}, IMAS.ec_launchers__beam___efficiency, :transmission) return ec_launchers end @@ -134,6 +132,7 @@ function ParametersInit(::Type{Val{:build}}) build.blanket = Entry(Float64, "", "Fraction of blanket in radial build") build.shield = Entry(Float64, "", "Fraction of shield in radial build") build.vessel = Entry(Float64, "", "Fraction of vessel in radial build") + build.plasma_gap = Entry(Real, "", "Fraction of vacuum gap between first wall and plasma separatrix in radial build"; default=0.1) build.symmetric = Entry(Bool, "", "Is the build up-down symmetric") build.n_first_wall_conformal_layers = Entry(Integer, "", "Number of layers that are conformal to the first wall"; default=1) return build @@ -163,3 +162,12 @@ function ParametersInit(::Type{Val{:coil_technology}}) coil_tech.ratio_SC_to_copper = Entry(Real, "", "Fraction of superconductor to copper cross-sectional areas") return coil_tech end + +function ParametersInit(::Type{Val{:target}}) + target = ParametersInit(nothing) + target.power_electric_net = Entry(Real, "W", "Target net electric power generated by the fusion power plant") + target.flattop_duration = Entry(Real, "s", "Target duration of the flattop (use Inf for steady-state)") + target.tritium_breeding_ratio = Entry(Real, "", "Target tritium breeding ratio of the whole plant") + target.cost = Entry(Real, "\$M", "Target total FPP cost") + return target +end diff --git a/src/physics.jl b/src/physics.jl index 96298bf5c..fdb2a87cb 100644 --- a/src/physics.jl +++ b/src/physics.jl @@ -176,7 +176,7 @@ function optimize_shape(r_obstruction, z_obstruction, target_clearance, func, r_ return shape_parameters end - function cost_TF_shape(r_obstruction, z_obstruction, rz_obstruction, target_clearance, func, r_start, r_end, shape_parameters; verbose=false) + function cost_shape(r_obstruction, z_obstruction, rz_obstruction, target_clearance, func, r_start, r_end, shape_parameters; verbose=false) R, Z = func(r_start, r_end, shape_parameters...) # disregard near r_start and r_end where optimizer has no control and shape is allowed to go over obstruction @@ -195,7 +195,7 @@ function optimize_shape(r_obstruction, z_obstruction, target_clearance, func, r_ cost_mean_distance = mean_distance_error / target_clearance # favor up/down symmetric solutions - cost_up_down_symmetry = abs(maximum(Z) + minimum(Z)) / maximum(abs.(Z)) + cost_up_down_symmetry = abs(maximum(Z) + minimum(Z)) / (maximum(Z) - minimum(Z)) if verbose @show minimum_distance @@ -213,10 +213,10 @@ function optimize_shape(r_obstruction, z_obstruction, target_clearance, func, r_ rz_obstruction = collect(zip(r_obstruction, z_obstruction)) initial_guess = copy(shape_parameters) - # res = optimize(shape_parameters-> cost_TF_shape(r_obstruction, z_obstruction, rz_obstruction, target_clearance, func, r_start, r_end, shape_parameters), + # res = optimize(shape_parameters-> cost_shape(r_obstruction, z_obstruction, rz_obstruction, target_clearance, func, r_start, r_end, shape_parameters), # initial_guess, Newton(), Optim.Options(time_limit=time_limit); autodiff=:forward) - res = Optim.optimize(shape_parameters -> cost_TF_shape(r_obstruction, z_obstruction, rz_obstruction, target_clearance, func, r_start, r_end, shape_parameters), - initial_guess, length(shape_parameters) == 1 ? Optim.BFGS() : Optim.NelderMead(), Optim.Options(time_limit=time_limit); autodiff=:forward) + res = Optim.optimize(shape_parameters -> cost_shape(r_obstruction, z_obstruction, rz_obstruction, target_clearance, func, r_start, r_end, shape_parameters), + initial_guess, length(shape_parameters) == 1 ? Optim.BFGS() : Optim.NelderMead(), Optim.Options(time_limit=time_limit)) if verbose println(res) end @@ -225,7 +225,7 @@ function optimize_shape(r_obstruction, z_obstruction, target_clearance, func, r_ # plot(func(r_start, r_end, initial_guess...); markershape=:x) # plot!(r_obstruction, z_obstruction, ; markershape=:x) # display(plot!(R, Z; markershape=:x, aspect_ratio=:equal)) - # cost_TF_shape(r_obstruction, z_obstruction, rz_obstruction, obstruction_area, target_clearance, func, r_start, r_end, shape_parameters; verbose=true) + # cost_shape(r_obstruction, z_obstruction, rz_obstruction, obstruction_area, target_clearance, func, r_start, r_end, shape_parameters; verbose=true) return shape_parameters end @@ -549,7 +549,7 @@ function layer_structure_intersect_volume(layer::IMAS.build__layer, structure::I end vol = 0.0 - for poly in LibGEOS.coordinates(LibGEOS.intersection(ring_poly, structure_poly)) + for poly in GeoInterface.coordinates(LibGEOS.intersection(ring_poly, structure_poly)) pr = [v[1] for v in poly] pz = [v[2] for v in poly] vol += IMAS.area(pr, pz) * structure.toroidal_extent * length(toroidal_angles) @@ -612,9 +612,9 @@ end function add_xpoint(mr::AbstractVector{T}, mz::AbstractVector{T}, i::Integer, R0::T, Z0::T, α::T) where {T<:Real} RX = mr[i] .* α .+ R0 .* (1.0 .- α) ZX = mz[i] .* α .+ Z0 .* (1.0 .- α) - RZ = FUSE.convex_hull(collect(zip(vcat(mr, RX), vcat(mz, ZX))); closed_polygon=true) - R = [r for (r, z) in RZ] - Z = [z for (r, z) in RZ] + RZ = FUSE.convex_hull(vcat(mr, RX), vcat(mz, ZX); closed_polygon=true) + R = T[r for (r, z) in RZ] + Z = T[z for (r, z) in RZ] return RX, ZX, R, Z end @@ -631,9 +631,13 @@ Control of the X-point location can be achieved by modifying R0, Z0 """ function add_xpoint(mr::AbstractVector{T}, mz::AbstractVector{T}, R0::Union{Nothing,T}=nothing, Z0::Union{Nothing,T}=nothing; upper::Bool) where {T<:Real} - function cost(mr, mz, i, R0, Z0, α) - RX, ZX, R, Z = add_xpoint(mr, mz, i, R0, Z0, α[1]) - return (1.0 - maximum(abs.(IMAS.curvature(R, Z) .* (upper ? Z .> Z0 : Z .< Z0))))^2 + function cost(mr::AbstractVector{T}, mz::AbstractVector{T}, i::Integer, R0::T, Z0::T, α::Float64) + RX, ZX, R, Z = add_xpoint(mr, mz, i, R0, Z0, α) + if upper + return (1.0 - maximum(abs.(IMAS.curvature(R, Z) .* (Z .> Z0))))^2.0 + else + return (1.0 - maximum(abs.(IMAS.curvature(R, Z) .* (Z .< Z0))))^2.0 + end end if Z0 === nothing @@ -674,7 +678,7 @@ function MXH_boundary(mxh::IMAS.MXH; upper_x_point::Bool, lower_x_point::Bool, n push!(ZX, ZXL) end - RZ = FUSE.convex_hull(collect(zip(vcat(mr, RX), vcat(mz, ZX))); closed_polygon=true) + RZ = FUSE.convex_hull(vcat(mr, RX), vcat(mz, ZX); closed_polygon=true) R = [r for (r, z) in RZ] Z = [z for (r, z) in RZ] diff --git a/src/technology.jl b/src/technology.jl index de6df39af..354c7f071 100644 --- a/src/technology.jl +++ b/src/technology.jl @@ -274,7 +274,7 @@ function solve_1D_solid_mechanics!( f_tf_sash=0.873, # : (float), conversion factor from hoop stress to axial stress for TF coil (nominally 0.873) f_oh_sash=0.37337, # : (float), conversion factor from hoop stress to axial stress for OH coil (nominally 0.37337) n_points::Integer=21, # : (int), number of radial points - verbose::Bool=false, # : (bool), flag for verbose output to terminal + verbose::Bool=false # : (bool), flag for verbose output to terminal ) tp = typeof(promote(R0, B0, R_tf_in, R_tf_out, Bz_oh, R_oh_in, R_oh_out)[1]) @@ -535,11 +535,11 @@ function solve_1D_solid_mechanics!( hoop_stress_tf = sh(r_tf, em_tf, gam_tf, displacement_tf, ddiplacementdr_tf) if axial_stress_tf_avg === nothing - hoop_stress_tf_avg = sum(hoop_stress_tf)/length(hoop_stress_tf) + hoop_stress_tf_avg = sum(hoop_stress_tf) / length(hoop_stress_tf) axial_stress_tf_avg = -f_tf_sash * hoop_stress_tf_avg end if axial_stress_oh_avg === nothing - hoop_stress_oh_avg = sum(hoop_stress_oh)/length(hoop_stress_oh) + hoop_stress_oh_avg = sum(hoop_stress_oh) / length(hoop_stress_oh) axial_stress_oh_avg = -f_oh_sash * hoop_stress_oh_avg end diff --git a/src/utils.jl b/src/utils.jl index ebe3833bb..47d1bd192 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1,3 +1,4 @@ +using InteractiveUtils: subtypes import ForwardDiff function unwrap(v, inplace=false) @@ -81,45 +82,7 @@ end # ****************************************** # Convex Hull # ****************************************** -struct Point - x::Float64 - y::Float64 -end - -function Base.isless(p::Point, q::Point) - p.x < q.x || (p.x == q.x && p.y < q.y) -end - -function isrightturn(p::Point, q::Point, r::Point) - (q.x - p.x) * (r.y - p.y) - (q.y - p.y) * (r.x - p.x) < 0 -end - -function halfhull(points::Vector{Point}) - halfhull = points[1:2] - for p in points[3:end] - push!(halfhull, p) - while length(halfhull) > 2 && !isrightturn(halfhull[end-2:end]...) - deleteat!(halfhull, length(halfhull) - 1) - end - end - halfhull -end - -function grahamscan(points::Vector{Point}) - sorted = sort(points) - upperhull = halfhull(sorted) - lowerhull = halfhull(reverse(sorted)) - [upperhull..., lowerhull[2:end-1]...] -end - -function convex_hull(xy::Vector; closed_polygon::Bool) - tmp = [(k.x, k.y) for k in grahamscan([Point(xx, yx) for (xx, yx) in xy])] - if closed_polygon - return push!(tmp, tmp[1]) - else - return tmp - end -end +import VacuumFields: convex_hull # ****************************************** # TraceCAD @@ -204,3 +167,16 @@ function fuse() ╚═╝ ╚═════╝ ╚══════╝╚══════╝ """ end + +""" + warmup() + +Function to execute to precompile the majority of FUSE +""" +function warmup() + ini, act = FUSE.case_parameters(:FPP; version=:v1_demount, init_from=:scalars) + dd = IMAS.dd() + FUSE.init(dd, ini, act) + FUSE.ActorWholeFacility(dd,act) + IMAS.freeze(dd) +end \ No newline at end of file diff --git a/test/runtests_parameters.jl b/test/runtests_parameters.jl index f71169299..33f6cc2c2 100644 --- a/test/runtests_parameters.jl +++ b/test/runtests_parameters.jl @@ -1,7 +1,6 @@ using Revise using FUSE using Test -using InteractiveUtils: subtypes @testset "ParametersInit" begin par = FUSE.ParametersAllInits() @@ -54,4 +53,8 @@ using InteractiveUtils: subtypes @test typeof(FUSE.ParametersActor(par)) <: FUSE.ParametersActor end + # save load + FUSE.dict2par!(FUSE.par2dict(ini), FUSE.ParametersAllInits()) + FUSE.dict2par!(FUSE.par2dict(act), FUSE.ParametersAllActors()) + end