Add fast density calculation to IMAS.sources

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,37 +593,37 @@ 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:
     # bkeflog = logy*(3 + emzrat) +
     #           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
 

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

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
 

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)