diff --git a/src/extract/objectives.jl b/src/extract/objectives.jl index c6d862b1..3aa91d62 100644 --- a/src/extract/objectives.jl +++ b/src/extract/objectives.jl @@ -81,9 +81,35 @@ end # ========================= # const ObjectiveFunctionsLibrary = OrderedCollections.OrderedDict{Symbol,ObjectiveFunction}() +function calculate_greenwald_fraction(dd::IMAS.dd) + try + eqt = dd.equilibrium.time_slice[] + cp1d = dd.core_profiles.profiles_1d[] + ne_line = IMAS.ne_line(eqt, cp1d) + Ip_MA = eqt.global_quantities.ip / 1e6 # Convert to MA + a_minor = eqt.boundary.minor_radius # Get minor radius from equilibrium + n_Greenwald = (Ip_MA / (π * a_minor^2)) * 1e20 # Greenwald limit [m⁻³] + return ne_line / n_Greenwald + catch + return 0.0 + end +end + +function calculate_bootstrap_fraction(dd::IMAS.dd) + try + I_bootstrap = @ddtime(dd.summary.global_quantities.current_bootstrap.value) + I_p = dd.equilibrium.time_slice[].global_quantities.ip + return I_bootstrap / I_p + catch + return 0.0 + end +end + function update_ObjectiveFunctionsLibrary!() empty!(ObjectiveFunctionsLibrary) #! format: off + + # Original FUSE objectives ObjectiveFunction(:min_levelized_CoE, "\$/kWh", dd -> dd.costing.levelized_CoE, -Inf) ObjectiveFunction(:min_log10_levelized_CoE, "log₁₀(\$/kW)", dd -> log10(dd.costing.levelized_CoE), -Inf) ObjectiveFunction(:min_capital_cost, "\$B", dd -> dd.costing.cost_direct_capital.cost / 1E3, -Inf) @@ -97,6 +123,14 @@ function update_ObjectiveFunctionsLibrary!() ObjectiveFunction(:min_βn, "", dd -> dd.equilibrium.time_slice[].global_quantities.beta_normal, -Inf) ObjectiveFunction(:min_R0, "m", dd -> dd.equilibrium.time_slice[].boundary.geometric_axis.r, -Inf) ObjectiveFunction(:max_zeff, "", dd -> @ddtime(dd.summary.volume_average.zeff.value), Inf) + + # New physics objectives for my studies + ObjectiveFunction(:max_βp, "", dd -> dd.equilibrium.time_slice[].global_quantities.beta_pol, Inf) + ObjectiveFunction(:max_h98, "", dd -> @ddtime(dd.summary.global_quantities.h_98.value), Inf) + ObjectiveFunction(:max_tau_e, "s", dd -> @ddtime(dd.summary.global_quantities.tau_energy.value), Inf) + ObjectiveFunction(:max_fbs, "", dd -> calculate_bootstrap_fraction(dd), Inf) + ObjectiveFunction(:max_greenwald_fraction, "", dd -> calculate_greenwald_fraction(dd), Inf) + #! format: on return ObjectiveFunctionsLibrary end diff --git a/src/get_from.jl b/src/get_from.jl index 80ad7121..ce80bc3e 100644 --- a/src/get_from.jl +++ b/src/get_from.jl @@ -121,6 +121,20 @@ function get_from(dd::IMAS.dd{T}, what::Val{:zeff_ped}, from_where::Symbol, rho_ return error("`get_from(dd, $what, Val(:$from_where))` doesn't exist yet") end +# ne_sep [m^-3] +function get_from(dd::IMAS.dd{T}, what::Val{:ne_sep}, from_where::Symbol; time0::Float64=dd.global_time)::T where {T<:Real} + if from_where == :core_profiles + cp1d = dd.core_profiles.profiles_1d[time0] + return cp1d.electrons.density_thermal[end] + elseif from_where == :pulse_schedule + if !ismissing(dd.pulse_schedule.density_control.n_e_separatrix, :reference) + return get_time_array(dd.pulse_schedule.density_control.n_e_separatrix, :reference, time0, :linear) + end + error("`get_from(dd, $what, Val{:$from_where})` does not have data") + end + return error("`get_from(dd, $what, Val{:$from_where})` doesn't exist yet") +end + Base.Docs.@doc """ get_from(dd::IMAS.dd, what::Symbol, from_where::Symbol; time0::Float64=dd.global_time) @@ -141,6 +155,8 @@ Supported quantities for `what`: - Possible sources (`from_where`): `:core_profiles`, `:summary`, `:pulse_schedule` - `:zeff_ped` - Effective charge at the pedestal [-] - Possible sources (`from_where`): `:core_profiles`, `:summary`, `:pulse_schedule` +- `:ne_sep` - Electron density at the separatrix [m^-3] + - Possible sources (`from_where`): `:core_profiles`, `:pulse_schedule` `time0` defines the time point at which to retrieve the value, default is `dd.global_time`. diff --git a/src/physics/fast.jl b/src/physics/fast.jl index ce69855b..ec095998 100644 --- a/src/physics/fast.jl +++ b/src/physics/fast.jl @@ -179,7 +179,7 @@ function critical_energy(cp1d::IMAS.core_profiles__profiles_1d, mf::Real, Zf::Re ne = cp1d.electrons.density_thermal Te = cp1d.electrons.temperature - avg_cmr = sum(ion.density_thermal .* (avgZ(ion.element[1].z_n, sum(ion.temperature)/length(ion.temperature))^2) / ion.element[1].a for ion in cp1d.ion) ./ ne + avg_cmr = sum(ion.density_thermal .* (avgZ(ion.element[1].z_n, sum(ion.temperature) / length(ion.temperature))^2) / ion.element[1].a for ion in cp1d.ion) ./ ne Ec = 14.8 .* mf .* Te .* avg_cmr .^ 1.5 if !approximate Ec1 = critical_energy(cp1d, 1, mf, Zf; approximate) @@ -292,7 +292,7 @@ Fills the core_profiles fast ion densities and pressures that result from fast i This calculation is done based on the `slowing_down_time` and `thermalization_time` of the fast ion species. """ -function fast_particles_profiles!(cs::IMAS.core_sources, cp1d::IMAS.core_profiles__profiles_1d; verbose::Bool=false) +function fast_particles_profiles!(cs::IMAS.core_sources, cp1d::IMAS.core_profiles__profiles_1d; verbose::Bool=false, max_fast_frac::Float64=0.9) ne = cp1d.electrons.density_thermal Te = cp1d.electrons.temperature @@ -303,6 +303,7 @@ function fast_particles_profiles!(cs::IMAS.core_sources, cp1d::IMAS.core_profile for ion in cp1d.ion ion.pressure_fast_parallel = zeros(Npsi) ion.pressure_fast_perpendicular = zeros(Npsi) + ion.density_thermal .+= ion.density_fast ion.density_fast = zeros(Npsi) end @@ -336,12 +337,17 @@ function fast_particles_profiles!(cs::IMAS.core_sources, cp1d::IMAS.core_profile taus = slowing_down_time.(ne, Te, particle_mass, particle_charge) taut = thermalization_time(cp1d, particle_energy, particle_mass, particle_charge) _, i4 = estrada_I_integrals(cp1d, particle_energy, particle_mass, particle_charge) - + + density = cion.density sion_particles = interp1d(source1d.grid.rho_tor_norm, sion.particles).(cp1d.grid.rho_tor_norm) pressa = i4 .* taus .* 2.0 ./ 3.0 .* (sion_particles .* particle_energy .* mks.e) cion.pressure_fast_parallel .+= pressa ./ 3.0 cion.pressure_fast_perpendicular .+= pressa ./ 3.0 cion.density_fast .+= sion_particles .* taut + limit = max_fast_frac .* cion.density_thermal + mask = cion.density_fast .> limit + cion.density_fast[mask] .= limit[mask] + cion.density_thermal = density - cion.density_fast end end end @@ -601,22 +607,22 @@ end """ bkefun(y::T1, vcvo::T2, tstcx::T3, emzrat::T4) where {T1<:Real, T2<:Real, T3<:Real, T4<:Real} """ -function bkefun(y::T1, vcvo::T2, tstcx::T3, emzrat::T4) where {T1<:Real, T2<:Real, T3<:Real, T4<:Real} +function bkefun(y::T1, vcvo::T2, tstcx::T3, emzrat::T4) where {T1<:Real,T2<:Real,T3<:Real,T4<:Real} T = promote_type(T1, T2, T3, T4) y <= 0 && return zero(T) - y2 = y*y - y3 = y2*y - v2 = vcvo*vcvo - v3 = v2*vcvo + y2 = y * y + y3 = y2 * y + v2 = vcvo * vcvo + v3 = v2 * vcvo y3v3 = y3 + v3 inv_y3v3 = inv(y3v3) - arg = (1 + v3)*inv_y3v3 # (1+v³)/(y³+v³) + arg = (1 + v3) * inv_y3v3 # (1+v³)/(y³+v³) # log(arg), but with a cheap accuracy win when arg ≈ 1 - alogarg = log(arg) + alogarg = log(arg) - logy = log(y) + logy = log(y) logy3v3 = log(y3v3) # algebraic simplification of the original expression: @@ -624,14 +630,14 @@ function bkefun(y::T1, vcvo::T2, tstcx::T3, emzrat::T4) where {T1<:Real, T2<:Rea # alogarg*(emzrat - tstcx)/3 - # logy3v3 bkeflog = muladd(logy, 3 + emzrat, - muladd(alogarg, (emzrat - tstcx)/3, -logy3v3)) + muladd(alogarg, (emzrat - tstcx) / 3, -logy3v3)) return bkeflog < -30 ? zero(T) : exp(bkeflog) end """ ion_momentum_fraction(vpar::T1, tpar::T2, emzpar::T3; N::Int=100) where {T1<:Real, T2<:Real, T3<:Real} """ -function ion_momentum_fraction(vpar::T1, tpar::T2, emzpar::T3; N::Int=100) where {T1<:Real, T2<:Real, T3<:Real} +function ion_momentum_fraction(vpar::T1, tpar::T2, emzpar::T3; N::Int=100) where {T1<:Real,T2<:Real,T3<:Real} T = promote_type(T1, T2, T3) out = zero(T) @inbounds @simd for y1 in range(0.0, 1.0, N) @@ -649,7 +655,7 @@ function ion_momentum_slowingdown_time(cp1d::IMAS.core_profiles__profiles_1d, Ef Te = cp1d.electrons.temperature a = cp1d.ion[1].element[1].a Zeff = cp1d.zeff - tau_mom = @. slowing_down_time(ne, Te, mf, Zf) * ion_momentum_fraction(sqrt(E_c / Ef), 0.0, a * Zeff / (mf * Zf)) + tau_mom = @. slowing_down_time(ne, Te, mf, Zf) * ion_momentum_fraction(sqrt(E_c / Ef), 0.0, a * Zeff / (mf * Zf)) return tau_mom end diff --git a/src/physics/profiles.jl b/src/physics/profiles.jl index 49d6f809..957280d6 100644 --- a/src/physics/profiles.jl +++ b/src/physics/profiles.jl @@ -1054,7 +1054,7 @@ function enforce_quasi_neutrality!(cp1d::IMAS.core_profiles__profiles_1d, specie end q_density_difference = ne .- sum(ion.density .* avgZ(ion) for ion in cp1d.ion if ion !== ion0) if hasdata(ion0, :density_fast) - q_density_difference .-= .-ion0.density_fast .* avgZ(ion0) + q_density_difference .-= ion0.density_fast .* avgZ(ion0) end # positive difference is assigned to target ion density_thermal diff --git a/src/physics/sources.jl b/src/physics/sources.jl index 07661ea9..2403a452 100644 --- a/src/physics/sources.jl +++ b/src/physics/sources.jl @@ -130,7 +130,7 @@ push!(document[Symbol("Physics sources")], :radiation_source!) Calculates intrisic sources and sinks, and adds them to `dd.core_sources` """ -function sources!(dd::IMAS.dd; bootstrap::Bool=true, DD_fusion::Bool=false) +function sources!(dd::IMAS.dd; bootstrap::Bool=true, DD_fusion::Bool=false, fast_ion_densities::Bool=true) if bootstrap bootstrap_source!(dd) # current end @@ -143,6 +143,10 @@ function sources!(dd::IMAS.dd; bootstrap::Bool=true, DD_fusion::Bool=false) fusion_source!(dd; DD_fusion) # electron and ion energy, particles + if fast_ion_densities + fast_particles_profiles!(dd) # update fast ion densities + end + return nothing end diff --git a/src/physics/transport.jl b/src/physics/transport.jl index 8f69c289..124cbfab 100644 --- a/src/physics/transport.jl +++ b/src/physics/transport.jl @@ -71,7 +71,7 @@ function profile_from_rotation_shear_transport( transport_indices = [IMAS.argmin_abs(rho, rho_x) for rho_x in transport_grid] index_ped = IMAS.argmin_abs(rho, rho_ped) index_last = transport_indices[end] - + if index_ped > index_last # If rho_ped is beyond transport_grid[end], extend interpolation dw_dr_old = gradient(rho, profile_old; method=:backward) @@ -89,7 +89,7 @@ function profile_from_rotation_shear_transport( # Set up the profile profile_new = similar(profile_old) profile_new[index_last:end] .= @views profile_old[index_last:end] - + # Integrate rotation shear from boundary inward to axis # Integration: ω(rho) = ω(rho_boundary) - ∫_{rho}^{rho_boundary} (dω/dr') dr' profile_new[index_last] = profile_old[index_last] @@ -121,7 +121,7 @@ function total_fluxes( rho_total_fluxes::AbstractVector{<:Real}; time0::Float64) where {T<:Real} total_flux1d = core_transport__model___profiles_1d{T}() - total_fluxes!(total_flux1d, core_transport, cp1d, rho_total_fluxes; time0) + return total_fluxes!(total_flux1d, core_transport, cp1d, rho_total_fluxes; time0) end function total_fluxes!( @@ -150,7 +150,7 @@ function total_fluxes!( end # defines paths to fill - paths = Union{Tuple{Symbol, Symbol}, Tuple{Symbol, Int64, Symbol}, Tuple{Symbol}}[] + paths = Union{Tuple{Symbol,Symbol},Tuple{Symbol,Int64,Symbol},Tuple{Symbol}}[] push!(paths, (:electrons, :energy)) push!(paths, (:electrons, :particles)) for k in eachindex(total_flux1d.ion) @@ -222,5 +222,63 @@ function total_fluxes(dd::IMAS.dd, rho_total_fluxes::AbstractVector{<:Real}=dd.c return total_fluxes(dd.core_transport, dd.core_profiles.profiles_1d[], rho_total_fluxes; time0) end + +""" + calculate_diffusivities(dd; ne=:none, Te=:none, Ti=:none, omega=:none) + +Compute transport particle, heat, and rotation diffusivities from desired profiles. +""" +function calculate_diffusivities(dd::IMAS.dd{T}; ne::Vector{T}=T[], Te::Vector{T}=T[], Ti::Vector{T}=T[], ω::Vector{T}=T[]) where {T<:Real} + + cs1d = IMAS.total_sources(dd) + cp1d = dd.core_profiles.profiles_1d[] + eqt = dd.equilibrium.time_slice[] + eqt1d = eqt.profiles_1d + rho_cp = cs1d.grid.rho_tor_norm + rho_eq = eqt1d.rho_tor_norm + + # Volumes and surfaces + surf = IMAS.interp1d(rho_eq, eqt1d.surface).(rho_cp) + + # Default profiles if not provided + ne = isempty(ne) ? cp1d.electrons.density_thermal : ne + Te = isempty(Te) ? cp1d.electrons.temperature : Te + Ti = isempty(Ti) ? cp1d.ion[1].temperature : Ti + ω = isempty(ω) ? cp1d.rotation_frequency_tor_sonic : ω + + # omega currently not used + zeff = cp1d.zeff + + # Logarithmic gradients + dlntedr = .-IMAS.calc_z(rho_cp, Te, :backward) + dlnnedr = .-IMAS.calc_z(rho_cp, ne, :backward) + dlntidr = .-IMAS.calc_z(rho_cp, Ti, :backward) + + # Pressure gradient terms + dpe = @. ne * IMAS.mks.e * Te * (dlntedr + dlnnedr) + dpi = @. ne * IMAS.mks.e * Ti * (dlntidr + dlnnedr) / zeff + + # momentum pressure + mi = cp1d.ion[1].element[1].a * IMAS.mks.m_p + Rmaj = eqt.boundary.geometric_axis.r # major radius + dωdr = .-IMAS.calc_z(rho_cp, ω, :backward) .* ω + pφ = dωdr .* ne .* mi .* Rmaj .^ 2 + + # Diffusivities + Dn = @. cs1d.electrons.particles_inside / (surf * dlnnedr * ne) + χe = @. cs1d.electrons.power_inside / (surf * dpe) + χi = @. cs1d.total_ion_power_inside / (surf * dpi) + χφ = @. cs1d.torque_tor_inside / (surf * pφ) + + # Patch singular boundary values + Dn[1] = Dn[2] + χe[1] = χe[2] + χi[1] = χi[2] + χφ[1] = χφ[2] + + return Dn, χe, χi, χφ +end + + @compat public total_fluxes push!(document[Symbol("Physics transport")], :total_fluxes) \ No newline at end of file