diff --git a/src/physics/fast.jl b/src/physics/fast.jl index ce69855b..d4452dea 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,30 +337,21 @@ 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 end - - # re-evaluate all cp1d fast-ion expressions - IMAS.refreeze!(cp1d.electrons, :density) - IMAS.refreeze!(cp1d.electrons, :pressure) - for ion in cp1d.ion - IMAS.refreeze!(ion, :density) - IMAS.refreeze!(ion, :pressure) - end - IMAS.refreeze!(cp1d, :pressure_ion_total) - IMAS.refreeze!(cp1d, :pressure_parallel) - IMAS.refreeze!(cp1d, :pressure_perpendicular) - IMAS.refreeze!(cp1d, :pressure) - - # must decide how to handle quasineutrality end function fast_particles_profiles!(dd::IMAS.dd; verbose::Bool=false) @@ -601,22 +593,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 +616,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 +641,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