Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Threaded expressions #799

Merged
merged 3 commits into from
Jan 17, 2025
Merged

Threaded expressions #799

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
df72dbe
empty commit
orso82 Jan 16, 2025
f25f1f4
init D3D with experimental data
orso82 Jan 17, 2025
2e5e37b
init_expressions
orso82 Jan 17, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
91 changes: 80 additions & 11 deletions src/cases/D3D.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,12 +16,11 @@
`scenario_sources` keywods says whether core_sources will be taken from scenario or not
"""
function case_parameters(::Type{Val{:D3D}};
scenario::Union{Symbol,AbstractString}=:default,
scenario::Union{Symbol,AbstractString,Tuple{Int,Float64}}=:default,
scenario_sources::Bool=true,
scenario_core_profiles::Bool=true,
realistic_pf_active::Bool=true)::Tuple{ParametersAllInits,ParametersAllActors}

@assert scenario in (:default, :H_mode, :L_mode) || isfile(scenario)
realistic_pf_active::Bool=true
)::Tuple{ParametersAllInits,ParametersAllActors}

ini = ParametersInits()
act = ParametersActors()
Expand All @@ -30,9 +29,73 @@
ini.general.init_from = :ods
ini.equilibrium.boundary_from = :ods

ini.time.simulation_start = 0.0
shot_ods_dir = ""
shot = 0
if typeof(scenario) <: Symbol
@assert scenario in (:default, :H_mode, :L_mode)
scenario_selector = scenario
elseif typeof(scenario) <: AbstractString
@assert isfile(scenario)
scenario_selector = :file

Check warning on line 40 in src/cases/D3D.jl

View check run for this annotation

Codecov / codecov/patch

src/cases/D3D.jl#L38-L40 

Added lines #L38 - L40 were not covered by tests
else
scenario_selector = :experiment
shot, time0 = scenario
ini.time.simulation_start = time0
shot_ods_dir = tempdir()
omas_py = """

Check warning on line 46 in src/cases/D3D.jl

View check run for this annotation

Codecov / codecov/patch

src/cases/D3D.jl#L42-L46 

Added lines #L42 - L46 were not covered by tests
print("Importing packages")
import time
import omas
from numpy import *

tic = time.time()
ods = omas.ODS()

print("Fetching ec_launcher data")
omas.omas_machine.d3d.ec_launcher_active_hardware(ods, $shot)

print("Fetching core_profiles data")
omas.omas_machine.d3d.core_profiles_profile_1d(ods, $shot, PROFILES_tree="ZIPFIT01")

print("Fetching equilibrium data")
with ods.open('d3d', $shot, options={'EFIT_tree': 'EFIT02'}):
ods["equilibrium.time"]
k = argmin(abs(ods["equilibrium.time"] - $time0))
ods["equilibrium.time"]
for k in range(len(ods["equilibrium.time"])):
ods["equilibrium.time_slice"][k]["time"]
ods["equilibrium.time_slice"][k]["global_quantities.ip"]
ods["equilibrium.time_slice"][k]["profiles_1d.psi"]
ods["equilibrium.time_slice"][k]["profiles_1d.f"]
ods["equilibrium.time_slice"][k]["profiles_1d.pressure"]
ods["equilibrium.time_slice"][k]["profiles_2d[0].psi"]
ods["equilibrium.time_slice"][k]["profiles_2d[0].grid.dim1"]
ods["equilibrium.time_slice"][k]["profiles_2d[0].grid.dim2"]
ods["equilibrium.time_slice"][k]["profiles_2d[0].grid_type.index"] = 1
ods["equilibrium.vacuum_toroidal_field.r0"]
ods["equilibrium.vacuum_toroidal_field.b0"]

toc = time.time()
print(f"Data fetched in {toc-tic} seconds")

print("Saving ODS to json")
tic = time.time()
ods.save("$shot_ods_dir/D3D_$shot.json")
toc = time.time()
print(f"Saved in {toc-tic} seconds")

"""
println(omas_py)
open(joinpath(shot_ods_dir, "omas_data_fetch.py"), "w") do io

Check warning on line 90 in src/cases/D3D.jl

View check run for this annotation

Codecov / codecov/patch

src/cases/D3D.jl#L89-L90 

Added lines #L89 - L90 were not covered by tests
return write(io, omas_py)
end
Base.run(`python -u $(joinpath(shot_ods_dir,"omas_data_fetch.py"))`)

Check warning on line 93 in src/cases/D3D.jl

View check run for this annotation

Codecov / codecov/patch

src/cases/D3D.jl#L93 

Added line #L93 was not covered by tests
end

machine_description = joinpath("__FUSE__", "sample", "D3D_machine.json")
shot_mappings = Dict(
scenario => Dict(
:file => Dict(
:nbi_power => 0.0,
:filename => "$(machine_description),$(scenario)"
),
Expand All @@ -46,15 +109,22 @@
),
:default => Dict(
:nbi_power => 5.0e6,
:filename => "$(machine_description),$(joinpath("__FUSE__", "sample", "D3D_eq_ods.json"))")
:filename => "$(machine_description),$(joinpath("__FUSE__", "sample", "D3D_eq_ods.json"))"),
:experiment => Dict(
:nbi_power => 5.0e6,
:filename => "$(machine_description),$(joinpath(shot_ods_dir, "D3D_$shot.json"))")
)

ini.time.simulation_start = 0.0
ini.ods.filename = shot_mappings[scenario][:filename]
ini.ods.filename = shot_mappings[scenario_selector][:filename]
ini.general.dd = load_ods(ini; error_on_missing_coordinates=false)
if !realistic_pf_active
empty!(ini.general.dd.pf_active)
end
if scenario_selector == :experiment
tt = round.(ini.general.dd.equilibrium.time, digits=3)
ini.time.pulse_shedule_time_basis = range(tt[1], tt[end], Int(round((tt[end] - tt[1]) / minimum(diff(tt)))))
IMAS.flux_surfaces(ini.general.dd.equilibrium, IMAS.first_wall(ini.general.dd.wall)...)

Check warning on line 126 in src/cases/D3D.jl

View check run for this annotation

Codecov / codecov/patch

src/cases/D3D.jl#L124-L126 

Added lines #L124 - L126 were not covered by tests
end

ini.build.layers = OrderedCollections.OrderedDict(
:gap_plug => 1.2,
Expand All @@ -80,7 +150,6 @@
ini.build.layers[:hfs_gap_coils].coils_inside = [7:10; 16:19]
ini.build.layers[:lfs_gap_coils].coils_inside = [11:15; 20:24]


ini.oh.technology = :copper
ini.pf_active.technology = :copper
ini.tf.technology = :copper
Expand All @@ -101,10 +170,10 @@
ini.core_profiles.rot_core = 5E3
end

if isempty(ini.general.dd.core_profiles) || !scenario_sources
if isempty(ini.general.dd.core_sources) || !scenario_sources
empty!(ini.general.dd.core_sources)
resize!(ini.nb_unit, 1)
ini.nb_unit[1].power_launched = shot_mappings[scenario][:nbi_power]
ini.nb_unit[1].power_launched = shot_mappings[scenario_selector][:nbi_power]
ini.nb_unit[1].beam_energy = 80e3
ini.nb_unit[1].beam_mass = 2.0
ini.nb_unit[1].toroidal_angle = 18.0 * deg
Expand Down
27 changes: 14 additions & 13 deletions src/ddinit/init.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -206,20 +206,21 @@ function init(

# add strike point information to pulse_schedule
if ps_was_set
eqt = dd.equilibrium.time_slice[]
fw = IMAS.first_wall(dd.wall)
psi_first_open = IMAS.find_psi_boundary(eqt, fw.r, fw.z; raise_error_on_not_open=true).first_open
Rxx, Zxx, _ = IMAS.find_strike_points(eqt, fw.r, fw.z, psi_first_open, dd.divertors)
RXX = []
ZXX = []
for eqt in dd.equilibrium.time_slice
fw = IMAS.first_wall(dd.wall)
psi_first_open = IMAS.find_psi_boundary(eqt, fw.r, fw.z; raise_error_on_not_open=true).first_open
Rxx, Zxx, _ = IMAS.find_strike_points(eqt, fw.r, fw.z, psi_first_open, dd.divertors)
push!(RXX, Rxx)
push!(ZXX, Zxx)
end
N = maximum(collect(map(length, RXX)))
pc = dd.pulse_schedule.position_control
resize!(pc.strike_point, 4)
for k in 1:4
if k <= length(Rxx)
pc.strike_point[k].r.reference = fill(Rxx[k], size(pc.time))
pc.strike_point[k].z.reference = fill(Zxx[k], size(pc.time))
else
pc.strike_point[k].r.reference = zeros(size(pc.time))
pc.strike_point[k].z.reference = zeros(size(pc.time))
end
resize!(pc.strike_point, N)
for k in 1:N
pc.strike_point[k].r.reference = IMAS.interp1d(dd.equilibrium.time, [k<=length(Rxx) ? Rxx[k] : 0.0 for Rxx in RXX]).(pc.time)
pc.strike_point[k].z.reference = IMAS.interp1d(dd.equilibrium.time, [k<=length(Zxx) ? Zxx[k] : 0.0 for Zxx in ZXX]).(pc.time)
end
end

Expand Down
1 change: 1 addition & 0 deletions src/ddinit/init_core_profiles.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -205,6 +205,7 @@ function init_core_profiles!(
end

# temperatures
@assert Te_core !== missing || pressure_core !== missing "Must specify either `ini.core_profiles.Te_core` or `ini.equilibrium.pressure_core`"
if Te_core === missing
Te_core = pressure_core / (Ti_Te_ratio * ni_core + ne_core) / IMAS.mks.e
end
Expand Down
3 changes: 1 addition & 2 deletions src/ddinit/init_equilibrium.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -92,10 +92,9 @@
pc.boundary_outline[k].z.reference[1] = z
end
else
pushfirst!(pc.time, -Inf)

Check warning on line 95 in src/ddinit/init_equilibrium.jl

View check run for this annotation

Codecov / codecov/patch

src/ddinit/init_equilibrium.jl#L95 

Added line #L95 was not covered by tests
for (k, (r, z)) in enumerate(zip(pr, pz))
pushfirst!(pc.time, -Inf)
pushfirst!(pc.boundary_outline[k].r.reference, r)
pushfirst!(pc.time, -Inf)
pushfirst!(pc.boundary_outline[k].z.reference, z)
end
end
Expand Down
36 changes: 2 additions & 34 deletions src/ddinit/init_expressions.json
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,8 +1,4 @@
[
"balance_of_plant.Q_plant",
"balance_of_plant.power_electric_net",
"balance_of_plant.power_electric_plant_operation.total_power",
"build.tf.wedge_thickness",
"core_profiles.global_quantities.beta_tor_norm",
"core_profiles.global_quantities.current_bootstrap",
"core_profiles.global_quantities.current_non_inductive",
Expand All @@ -28,9 +24,6 @@
"core_sources.source[:].profiles_1d[:].grid.psi",
"core_sources.source[:].profiles_1d[:].grid.psi_norm",
"core_sources.source[:].profiles_1d[:].grid.surface",
"costing.cost_decommissioning.cost",
"costing.cost_direct_capital.cost",
"costing.cost_operations.yearly_cost",
"equilibrium.time_slice[:].boundary.elongation",
"equilibrium.time_slice[:].boundary.elongation_lower",
"equilibrium.time_slice[:].boundary.elongation_upper",
Expand Down Expand Up @@ -60,6 +53,7 @@
"equilibrium.time_slice[:].profiles_1d.dpsi_drho_tor",
"equilibrium.time_slice[:].profiles_1d.dvolume_dpsi",
"equilibrium.time_slice[:].profiles_1d.dvolume_drho_tor",
"equilibrium.time_slice[:].profiles_1d.f_df_dpsi",
"equilibrium.time_slice[:].profiles_1d.geometric_axis.r",
"equilibrium.time_slice[:].profiles_1d.geometric_axis.z",
"equilibrium.time_slice[:].profiles_1d.j_parallel",
Expand All @@ -69,38 +63,12 @@
"equilibrium.time_slice[:].profiles_2d[:].r",
"equilibrium.time_slice[:].profiles_2d[:].z",
"pulse_schedule.tf.b_field_tor_vacuum_r.reference",
"pulse_schedule.time",
"limits.all_cleared",
"summary.fusion.power.value",
"summary.global_quantities.beta_pol_mhd.value",
"summary.global_quantities.beta_tor.value",
"summary.global_quantities.beta_tor_mhd.value",
"summary.global_quantities.beta_tor_norm.value",
"summary.global_quantities.beta_tor_norm_mhd.value",
"summary.global_quantities.beta_tor_thermal_norm.value",
"summary.global_quantities.current_bootstrap.value",
"summary.global_quantities.current_non_inductive.value",
"summary.global_quantities.current_ohm.value",
"summary.global_quantities.energy_thermal.value",
"summary.global_quantities.h_98.value",
"summary.global_quantities.ip.value",
"summary.global_quantities.tau_energy.value",
"summary.global_quantities.tau_energy_98.value",
"summary.heating_current_drive.power_launched_ec.value",
"summary.heating_current_drive.power_launched_ic.value",
"summary.heating_current_drive.power_launched_lh.value",
"summary.heating_current_drive.power_launched_nbi.value",
"summary.heating_current_drive.power_launched_total.value",
"summary.local.magnetic_axis.n_e.value",
"summary.local.magnetic_axis.t_e.value",
"summary.local.magnetic_axis.t_i_average.value",
"summary.local.magnetic_axis.zeff.value",
"summary.local.separatrix.n_e.value",
"summary.local.separatrix.t_e.value",
"summary.local.separatrix.t_i_average.value",
"summary.local.separatrix.zeff.value",
"summary.volume_average.n_e.value",
"summary.volume_average.t_e.value",
"summary.volume_average.t_i_average.value",
"summary.volume_average.zeff.value"
"summary.local.separatrix.zeff.value"
]
Loading
Loading