Implemented SplitApplyCombine

LICENSE

@@ -1,6 +1,6 @@
 MIT License
 
-Copyright (c) 2022 moataz-sabry
+Copyright (c) 2023 moataz-sabry
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal

Manifest.toml

@@ -2,7 +2,7 @@
 
 julia_version = "1.8.3"
 manifest_format = "2.0"
-project_hash = "02b8adcbf7ffa58e50b1cd6b9fa5d6b318486a61"
+project_hash = "3cb89865eaa8c104f3136e377a234921cf7a5336"
 
 [[deps.AbstractFFTs]]
 deps = ["ChainRulesCore", "LinearAlgebra"]
@@ -117,6 +117,12 @@ git-tree-sha1 = "0fba8b706d0178b4dc7fd44a96a92382c9065c2c"
 uuid = "244e2a9f-e319-4986-a169-4d1fe445cd52"
 version = "0.1.2"
 
+[[deps.Dictionaries]]
+deps = ["Indexing", "Random", "Serialization"]
+git-tree-sha1 = "e82c3c97b5b4ec111f3c1b55228cebc7510525a2"
+uuid = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4"
+version = "0.3.25"
+
 [[deps.DiffResults]]
 deps = ["StaticArraysCore"]
 git-tree-sha1 = "782dd5f4561f5d267313f23853baaaa4c52ea621"
@@ -181,6 +187,11 @@ git-tree-sha1 = "2e99184fca5eb6f075944b04c22edec29beb4778"
 uuid = "7869d1d1-7146-5819-86e3-90919afe41df"
 version = "0.4.7"
 
+[[deps.Indexing]]
+git-tree-sha1 = "ce1566720fd6b19ff3411404d4b977acd4814f9f"
+uuid = "313cdc1a-70c2-5d6a-ae34-0150d3930a38"
+version = "1.1.1"
+
 [[deps.InitialValues]]
 git-tree-sha1 = "4da0f88e9a39111c2fa3add390ab15f3a44f3ca3"
 uuid = "22cec73e-a1b8-11e9-2c92-598750a2cf9c"
@@ -391,6 +402,12 @@ git-tree-sha1 = "d75bda01f8c31ebb72df80a46c88b25d1c79c56d"
 uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
 version = "2.1.7"
 
+[[deps.SplitApplyCombine]]
+deps = ["Dictionaries", "Indexing"]
+git-tree-sha1 = "48f393b0231516850e39f6c756970e7ca8b77045"
+uuid = "03a91e81-4c3e-53e1-a0a4-9c0c8f19dd66"
+version = "1.2.2"
+
 [[deps.SplittablesBase]]
 deps = ["Setfield", "Test"]
 git-tree-sha1 = "e08a62abc517eb79667d0a29dc08a3b589516bb5"

Project.toml

@@ -6,6 +6,7 @@ version = "1.0.0"
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+SplitApplyCombine = "03a91e81-4c3e-53e1-a0a4-9c0c8f19dd66"
 ThreadsX = "ac1d9e8a-700a-412c-b207-f0111f4b6c0d"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"

src/Apophis.jl

@@ -4,6 +4,7 @@ using Base.Iterators: filter, flatten, partition, take
 using Base: OneTo, rest
 using Combinatorics
 using SparseArrays
+using SplitApplyCombine: combinedimsview, mapview
 using Unitful
 using Zygote
 # using ThreadsX

src/calc.jl

@@ -9,9 +9,9 @@ update_reaction_rates(v::Val, (; mechanism, state)::Gas{<:Number}) = foreach(rea
 update_production_rates((; mechanism)::Gas{<:Number}) = foreach(species -> _update_production_rates(species), mechanism.species)
 update_production_rates(v::Val, (; mechanism)::Gas{<:Number}) = foreach(species -> _update_production_rates(v, species), mechanism.species)
 
-update(gas::Gas{<:Number}) = ((update_thermodynamics, update_reaction_rates, update_production_rates)(gas); return nothing)
-update(gas::Gas{<:Number}, d::Symbol) = ((update_thermodynamics, update_reaction_rates, update_production_rates)(Val(d), gas); return nothing) ## Val(d) allocates 1
-update(gas::Gas{<:Number}, ds::Vararg{Symbol}) = foreach(d -> update(gas, d), ds)
+update!(gas::Gas{<:Number}) = ((update_thermodynamics, update_reaction_rates, update_production_rates)(gas); return nothing)
+update!(gas::Gas{<:Number}, d::Symbol) = ((update_thermodynamics, update_reaction_rates, update_production_rates)(Val(d), gas); return nothing) ## Val(d) allocates 1
+update!(gas::Gas{<:Number}, ds::Vararg{Symbol}) = foreach(d -> update!(gas, d), ds)
 
 _kinetics_sensitivity((; mechanism, state)::Gas{<:Number}, d::Int = 1) = @inbounds map(reaction -> _kinetics_sensitivity(reaction, state)[d], mechanism.reactions)
 kinetics_sensitivity(gas::Gas{N}, d::Int = 1) where {N<:Number} = @inbounds _kinetics_sensitivity(gas, d) |> dqdk -> map(s -> sum(dqdk[i] * ν for ((; i), ν) in s.inreactions; init=zero(N)), species(gas)) |> Diagonal

src/utils.jl

@@ -1,5 +1,7 @@
 export 
-    Gas, State, Species, Reaction, TPC!, TρC!, TPY!, TPX!, TρY!, TρX!
+    Gas, Species, species, Reaction, reactions, TPC!, TρC!, TPY!, TPX!, TρY!, TρX!, update!, temperature, pressure, density, mass_fractions, molar_fractions, mass_fractions,
+    molecular_weights, production_rates, heat_capacities_pressure, average_heat_capacity_pressure, heat_capacities_volume, average_heat_capacity_volume,
+    enthalpies, internal_energies, entropies, forward_rates, reverse_rates, progress_rates, stoichiometry_matrix
 
 @generated (t::Union{NTuple{N, Function}, NTuple{N, Arrhenius}})(x...) where N = :($((:(t[$i](x...)) for i in OneTo(N))...),)
 
@@ -49,66 +51,63 @@ temperature((; state)::Gas) = state.T
 pressure((; state)::Gas; in=nothing) = state.P * (isnothing(in) || ustrip(in, 1u"dyn/cm^2"))
 density((; state)::Gas; in=nothing) = state.ρ * (isnothing(in) || ustrip(in, 1u"g/cm^3"))
 mass_fractions((; state)::Gas) = state.Y
+mass_fraction(gas::Gas, s::Union{String, Symbol}) = species(gas, s) |> s -> mass_fractions(gas)[s.k]
 molar_fractions((; state)::Gas) = state.X
-molar_concentrations((; state)::Gas) = state.C
+molar_concentrations((; state)::Gas; in=nothing) = isnothing(in) ? state.C : mapview(c -> c * ustrip(in, 1u"mol/cm^3"), state.C)
 
 species((; mechanism)::Gas) = mechanism.species
 species(gas::Gas, i::Int) = species(gas)[i]
 species(gas::Gas, s::Union{String, Symbol}) = assign(species(gas), s)
 
-molecular_weights(gas::Gas) = map(s -> s.weight, species(gas))
+molecular_weights(gas::Gas{<:Number}, f::Function=identity) = mapview(s -> f(s.weight), species(gas))
 stoichiometry_matrix(gas::Gas) = [(s.k, r.i, ν) for s in species(gas) for (r, ν) in s.inreactions] |> v -> sparse((getfield.(v, i) for i in OneTo(3))...)
 
-heat_capacity_pressure(species::Species; in=nothing) = heat_capacity_pressure(species, Val(:val); in)
-heat_capacity_pressure(species::Species, ::Val{:val}; in=nothing) = species.thermo.cₚ.val[] * (isnothing(in) || ustrip(in, 1u"erg/(mol*K)"))
-heat_capacity_pressure(species::Species, ::Val{:dT}; in=nothing) = species.thermo.cₚ.dT[] * (isnothing(in) || ustrip(in, 1u"erg/(mol*K^2)"))
-heat_capacities_pressure(gas::Gas, v::Val = Val(:val); in=nothing) = map(s -> heat_capacity_pressure(s, v; in), species(gas))
+heat_capacity_pressure(species::Species{<:Number}, ::Val{:val} = Val(:val); in=nothing) = species.thermo.cₚ.val[] * (isnothing(in) || ustrip(in, 1u"erg/(mol*K)"))
+heat_capacity_pressure(species::Species{<:Number}, ::Val{:dT}; in=nothing) = species.thermo.cₚ.dT[] * (isnothing(in) || ustrip(in, 1u"erg/(mol*K^2)"))
+heat_capacities_pressure(gas::Gas{<:Number}, v::Val = Val(:val); in=nothing) = mapview(s -> heat_capacity_pressure(s, v; in), species(gas))
+average_heat_capacity_pressure(gas::Gas{<:Number}, v::Val = Val(:val); in=nothing) = sum(cₚ * y / w for (cₚ, y, w) in zip(heat_capacities_pressure(gas, v; in = in * u"g/mol"), mass_fractions(gas), molecular_weights(gas)))
 
-average_heat_capacity_pressure(gas::Gas{<:Number}, f::Symbol = :val; in=nothing) = sum(heat_capacity_pressure(s, f) * y * inv(s.weight) for (s, y) in zip(species(gas), mass_fractions(gas))) * (isnothing(in) || ustrip(in, 1u"erg/(g*K)"))
+heat_capacities_volume(gas::Gas{<:Number}, ::Val{:val} = Val(:val); in=nothing) = mapview(cₚ -> cₚ - R, heat_capacities_pressure(gas; in))
+heat_capacities_volume(gas::Gas{<:Number}, v::Val{:dT}; in=nothing) = heat_capacities_pressure(gas, v; in)
+average_heat_capacity_volume(gas::Gas{<:Number}, v::Val = Val(:val); in=nothing) = sum(cᵥ * y / w for (cᵥ, y, w) in zip(heat_capacities_volume(gas, v; in = in * u"g/mol"), mass_fractions(gas), molecular_weights(gas)))
 
-heat_capacity_volume(species::Species, v::Val = Val(:val)) = heat_capacity_pressure(species, v) - R
-average_heat_capacity_volume(gas::Gas{<:Number}) = sum(heat_capacity_volume(s) * y * inv(s.weight) for (s, y) in zip(species(gas), mass_fractions(gas)))
+enthalpy(species::Species{<:Number}; in=nothing) = enthalpy(species, Val(:val); in)
+enthalpy(species::Species{<:Number}, ::Val{:val}; in=nothing) = species.thermo.h.val[] * (isnothing(in)  || ustrip(in, 1u"erg/mol"))
+enthalpy(species::Species{<:Number}, ::Val{:dT}; in=nothing) = species.thermo.h.dT[] * (isnothing(in) || ustrip(in, 1u"erg/(mol*K)"))
+enthalpies(gas::Gas{<:Number}, v::Val = Val(:val); in=nothing) = mapview(s -> enthalpy(s, v; in), species(gas)) 
 
-enthalpy(species::Species; in=nothing) = enthalpy(species, Val(:val); in)
-enthalpy(species::Species, ::Val{:val}; in=nothing) = species.thermo.h.val[] * (isnothing(in) || ustrip(in, 1u"erg/mol"))
-enthalpy(species::Species, ::Val{:dT}; in=nothing) = species.thermo.h.dT[] * (isnothing(in) || ustrip(in, 1u"erg/(mol*K)"))
-enthalpies(gas::Gas, v::Val = Val(:val); in=nothing) = map(s -> enthalpy(s, v; in), species(gas))
+internal_energies(gas::Gas{<:Number}, ::Val{:val} = Val(:val); in=nothing) = mapview(h -> h - R * temperature(gas), enthalpies(gas; in))
+internal_energies(gas::Gas{<:Number}, v::Val{:dT}; in=nothing) = mapview(h -> h - R, enthalpies(gas, v; in))
 
-internal_energy(gas::Gas, species::Species, v::Val = Val(:val)) = enthalpy(species, v) - R * temperature(gas)
+entropy(species::Species{<:Number}; in=nothing) = entropy(species, Val(:val); in)
+entropy(species::Species{<:Number}, ::Val{:val}; in=nothing) = species.thermo.s.val[] * (isnothing(in) || ustrip(in, 1u"erg/(mol*K)"))
+entropy(species::Species{<:Number}, ::Val{:dT}; in=nothing) = species.thermo.s.dT[] * (isnothing(in) || ustrip(in, 1u"erg/(mol*K^2)"))
+entropies(gas::Gas{<:Number}, v::Val = Val(:val); in=nothing) = mapview(s -> entropy(s, v; in), species(gas))
 
-entropy(species::Species; in=nothing) = entropy(species, Val(:val); in)
-entropy(species::Species, ::Val{:val}; in=nothing) = species.thermo.s.val[] * (isnothing(in) || ustrip(in, 1u"erg/(mol*K)"))
-entropy(species::Species, ::Val{:dT}; in=nothing) = species.thermo.s.dT[] * (isnothing(in) || ustrip(in, 1u"erg/(mol*K^2)"))
-entropies(gas::Gas, v::Val = Val(:val); in=nothing) = map(s -> entropy(s, v; in), species(gas))
-
-production_rate(species::Species; in=nothing) = production_rate(species, Val(:val); in)
-production_rate(species::Species, ::Val{:val}; in=nothing) = species.rates.ω̇.val[] * (isnothing(in) || ustrip(in, 1u"mol/(cm^3*s)"))
-production_rate(species::Species, ::Val{:dT}; in=nothing) = species.rates.ω̇.dT[] * (isnothing(in) || ustrip(in, 1u"mol/(cm^3*K*s)"))
-production_rate(species::Species, ::Val{:dC}; in=nothing) = species.rates.ω̇.dC * (isnothing(in) || ustrip(in, 1u"s^-1"))
 production_rate(gas::Gas, s::Union{String, Symbol}, v::Val = Val(:val); in=nothing) = species(gas, s) |> s -> production_rate(s, v, in)
-production_rates(gas::Gas{N}, v::Val = Val(:val); in=nothing) where {N<:Number} = map(s -> production_rate(s, v; in), species(gas)) |> s -> s isa Vector{Vector{N}} ? mapreduce(permutedims, vcat, s) : s
-production_rates!(dest::AbstractArray, gas::Gas{N}, ::Val{S}) where {N<:Number, S} = for s in species(gas), (j, d) in enumerate(getfield(s.rates.ω̇, S)) dest[s.k, j] = d end
+production_rate(species::Species{<:Number}, ::Val{:val} = Val(:val); in=nothing) = species.rates.ω̇.val[] * (isnothing(in) || ustrip(in, 1u"mol/(cm^3*s)"))
+production_rate(species::Species{<:Number}, ::Val{:dT}; in=nothing) = species.rates.ω̇.dT[] * (isnothing(in) || ustrip(in, 1u"mol/(cm^3*K*s)"))
+production_rate(species::Species{<:Number}, ::Val{:dC}; in=nothing) = species.rates.ω̇.dC
+production_rates(gas::Gas{<:Number}, v::Val = Val(:val); in=nothing) = mapview(s -> production_rate(s, v; in), species(gas)) |> x -> v isa Val{:dC} ? combinedimsview(x, 1) : x
+#production_rates(gas::Gas{<:Number}, v::Val{:dC}; in=nothing) =
 
 reactions((; mechanism)::Gas) = mechanism.reactions
 reaction(gas::Gas, i::Int) = reactions(gas)[i]
 reaction(gas::Gas, r::Union{String, Symbol}) = assign(reactions(gas), r)
 
 order_unit(reaction::Reaction{N}, part::Vector{Pair{Species{N}, N}}) where {N<:Number} = sum(abs ∘ last, part) + (reaction isa ElementaryReaction ? zero(N) : one(N))
 
-forward_rate(reaction::Reaction, f::Symbol = :val) = forward_rate(reaction, Val(f))
-forward_rate(reaction::Reaction, ::Val{:val}) = reaction.rates.kf.val[]
-forward_rate(reaction::Reaction, ::Val{:dT}) = reaction.rates.kf.dT[]
-forward_rate(reaction::Reaction, ::Val{:dC}) = reaction.rates.kf.dC
-forward_rates(gas::Gas{N}, f::Symbol = :val) where {N<:Number} = map(r -> forward_rate(r, f), reactions(gas)) |> r -> r isa Vector{Vector{N}} ? mapreduce(permutedims, vcat, r) : r
-
-reverse_rate(reaction::Reaction, f::Symbol = :val) = reverse_rate(reaction, Val(f))
-reverse_rate(reaction::Reaction, ::Val{:val}) = reaction.rates.kr.val[]
-reverse_rate(reaction::Reaction, ::Val{:dT}) = reaction.rates.kr.dT[]
-reverse_rate(reaction::Reaction, ::Val{:dC}) = reaction.rates.kr.dC
-reverse_rates(gas::Gas{N}, f::Symbol = :val) where {N<:Number} = map(r -> reverse_rate(r, f), reactions(gas)) |> r -> r isa Vector{Vector{N}} ? mapreduce(permutedims, vcat, r) : r
-
-progress_rate(reaction::Reaction, f::Symbol = :val; in=nothing) = progress_rate(reaction, Val(f); in)
-progress_rate(reaction::Reaction, ::Val{:val}; in=nothing) = reaction.rates.q.val[] * (isnothing(in) || ustrip(in, 1u"mol/(cm^3*s)"))
-progress_rate(reaction::Reaction, ::Val{:dT}; in=nothing) = reaction.rates.q.dT[] * (isnothing(in) || ustrip(in, 1u"mol/(cm^3*K*s)"))
-progress_rate(reaction::Reaction, ::Val{:dC}; in=nothing) = reaction.rates.q.dC * (isnothing(in) || ustrip(in, 1u"s^-1"))
-progress_rates(gas::Gas{N}, f::Symbol = :val; in=nothing) where {N<:Number} = map(r -> progress_rate(r, f; in), reactions(gas)) |> r -> r isa Vector{Vector{N}} ? mapreduce(permutedims, vcat, r) : r
+forward_rate(reaction::Reaction{<:Number}, ::Val{:val}=Val(:val)) = reaction.rates.kf.val[]
+forward_rate(reaction::Reaction{<:Number}, ::Val{:dT}) = reaction.rates.kf.dT[]
+forward_rate(reaction::Reaction{<:Number}, ::Val{:dC}) = reaction.rates.kf.dC
+forward_rates(gas::Gas{<:Number}, v::Val = Val(:val)) = mapview(r -> forward_rate(r, v), reactions(gas))
+
+reverse_rate(reaction::Reaction{<:Number}, ::Val{:val}=Val(:val)) = reaction.rates.kr.val[]
+reverse_rate(reaction::Reaction{<:Number}, ::Val{:dT}) = reaction.rates.kr.dT[]
+reverse_rate(reaction::Reaction{<:Number}, ::Val{:dC}) = reaction.rates.kr.dC
+reverse_rates(gas::Gas{<:Number}, v::Val = Val(:val)) = mapview(r -> forward_rate(r, v), reactions(gas))
+
+progress_rate(reaction::Reaction{<:Number}, ::Val{:val}=Val(:val); in=nothing) = reaction.rates.q.val[] * (isnothing(in) || ustrip(in, 1u"mol/(cm^3*s)"))
+progress_rate(reaction::Reaction{<:Number}, ::Val{:dT}; in=nothing) = reaction.rates.q.dT[] * (isnothing(in) || ustrip(in, 1u"mol/(cm^3*K*s)"))
+progress_rate(reaction::Reaction{<:Number}, ::Val{:dC}; in=nothing) = reaction.rates.q.dC
+progress_rates(gas::Gas{<:Number}, v::Val = Val(:val)) = mapview(r -> forward_rate(r, v), reactions(gas))

test/Manifest.toml

@@ -2,7 +2,7 @@
 
 julia_version = "1.8.3"
 manifest_format = "2.0"
-project_hash = "d9c53e86e2a09f376fee7818180905d4b58ec3aa"
+project_hash = "3eb3c26938ca79971810ff77468fc0f0d57806f3"
 
 [[deps.AbstractFFTs]]
 deps = ["ChainRulesCore", "LinearAlgebra"]

test/Project.toml

@@ -4,3 +4,4 @@ BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 Conda = "8f4d0f93-b110-5947-807f-2305c1781a2d"
 PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"