diff --git a/LICENSE b/LICENSE index b0fa61f..75e45d5 100644 --- a/LICENSE +++ b/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 diff --git a/Manifest.toml b/Manifest.toml index 0f00b3d..78b7839 100644 --- a/Manifest.toml +++ b/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" diff --git a/Project.toml b/Project.toml index 5281617..b88719b 100644 --- a/Project.toml +++ b/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" diff --git a/src/Apophis.jl b/src/Apophis.jl index a41ce88..74ea96c 100644 --- a/src/Apophis.jl +++ b/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 diff --git a/src/calc.jl b/src/calc.jl index 579f83f..b72a5f6 100644 --- a/src/calc.jl +++ b/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 \ No newline at end of file diff --git a/src/utils.jl b/src/utils.jl index fbac383..b630986 100644 --- a/src/utils.jl +++ b/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,45 +51,45 @@ 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] @@ -95,20 +97,17 @@ 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 \ No newline at end of file +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)) \ No newline at end of file diff --git a/test/Manifest.toml b/test/Manifest.toml index 2877649..ba174a9 100644 --- a/test/Manifest.toml +++ b/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"] diff --git a/test/Project.toml b/test/Project.toml index 9f46076..147184a 100644 --- a/test/Project.toml +++ b/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"