Skip to content

Commit

Permalink
API improvements
Browse files Browse the repository at this point in the history
  • Loading branch information
@moataz-sabry
moataz-sabry committed Feb 26, 2023
1 parent 9e1f518 commit d79fe15
Show file tree
Hide file tree
Showing 11 changed files with 194 additions and 136 deletions.
4 changes: 3 additions & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
@@ -1,16 +1,18 @@
name = "Apophis"
uuid = "cf72b054-1e85-4718-a3f4-02cfdc10a054"
authors = ["moataz-sabry <[email protected]>", "Ahmed Hassan <[email protected]>"]
version = "1.3.4"
version = "1.3.5"

[deps]
PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
SplitApplyCombine = "03a91e81-4c3e-53e1-a0a4-9c0c8f19dd66"
Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
YAML = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6"
Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"

[compat]
PrettyTables = "2.2.2"
PyCall = "1.95.1"
SplitApplyCombine = "1.2.2"
Unitful = "1.12.3"
Expand Down
29 changes: 20 additions & 9 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,35 +3,46 @@ Welcome to Apophis! This package is designed to provide sensitivity analysis of

## Installation
```julia
] add Apophis
pkg> add Apophis
```
## Quick Start
### Reading a Mechanism
```julia
using Apophis
julia> using Apophis

gas = Gas(:GRI3; T = 1000, P = 101325, Y = "CH4: 0.05, O2: 0.20, N2: 0.75")
julia> gas = Gas(:GRI3; T = 1000, P = 101325, Y = "CH4: 0.05, O2: 0.20, N2: 0.75")
# Creates a Gas based on the given mechanism data files

temperature(gas) ## T [K]
T: 1000 K // P: 101325 Pa // ρ: 0.3372 kg/m³
--------- ------ ------ --------- ------------ -------------- ------------
Species Y X C cₚ h s
– – kmol/m³ J/(kmol⋅K) J/kmol J/(kmol⋅K)
--------- ------ ------ --------- ------------ -------------- ------------
N2 0.75 0.74 0.0090 32762 2.14699e+07 228089
O2 0.20 0.17 0.0021 34883 2.27068e+07 243586
CH4 0.05 0.09 0.0011 73616.7 -3.59484e+07 248279
⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮
--------- ------ ------ --------- ------------ -------------- ------------
50 rows omitted
julia> temperature(gas) # T [K]
1000.0

density(gas) ## ρ [kg/m³]
julia> density(gas) # ρ [kg/m³]
0.337205
```
### Changing the State
```julia
TρY!(gas, 1250.0, 0.35, mass_fractions(gas))
julia> set!(gas; T = 1250.0, ρ = 0.35, Y = mass_fractions(gas)); # T && (P || ρ) && (Y || X || C)

pressure(gas) ## P [Pa]
julia> pressure(gas) # P [Pa]
131461.83
```
### Chemical Kinetics
```julia
update(gas)
julia> update(gas);
# Updates internal variables based on the current gas state

production_rates(gas) # ω̇ [kmol/(m³⋅s)]
julia> production_rates(gas) # ω̇ [kmol/(m³⋅s)]
53-element Vector{Float64}:
0.0
3.7072899145013144e-11
Expand Down
7 changes: 4 additions & 3 deletions src/Apophis.jl
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,8 @@
module Apophis

using Base.Iterators: filter, flatten, reverse, take, partition
using Base.Iterators: filter, flatten, partition, reverse, take
using Base: Fix1, Fix2, OneTo, rest
using PrettyTables
using SparseArrays
using SplitApplyCombine: combinedimsview, mapview
using Unitful
Expand All @@ -14,8 +15,8 @@ abstract type AbstractReaction{N<:Number} end
const R = 8.31446261815324e3 # J K⁻¹ kmol⁻¹
const Rc = 1.987261815324 # cal K⁻¹ mol⁻¹
const kB = 1.380649e-23 # J K⁻¹
const Pa = 101325.0 # Pa
const Tᵣ = 300.0 # K
const Pa = 101325 # Pa
const Tᵣ = 300 # K
const d = 0.14

const Maybe{T} = Union{T, Nothing}
Expand Down
4 changes: 2 additions & 2 deletions src/calc.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,6 @@ 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}) = ((update_thermodynamics, update_reaction_rates, update_production_rates)(gas); return gas)
update(gas::Gas{<:Number}, d::Symbol) = ((update_thermodynamics, update_reaction_rates, update_production_rates)(Val(d), gas); return gas) ## Val(d) allocates 1
update(gas::Gas{<:Number}, ds::Vararg{Symbol}) = foreach(d -> update(gas, d), ds)
46 changes: 23 additions & 23 deletions src/gas.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,33 +14,33 @@ end

State(mechanism::Mechanism{N}) where {N<:Number} = State(N(Tᵣ), N(Pa), zero(N), (zeros(N, length(mechanism.species)) for _ in OneTo(3))...)

XW2W̅(X::Vector{N}, W::AbstractVector{N}) where {N<:Number} = sum(x * w for (x, w) in zip(X, W))
YW2W̅(Y::Vector{N}, W::AbstractVector{N}) where {N<:Number} = sum(y / w for (y, w) in zip(Y, W)) |> inv
CW2W̅(C::Vector{N}, W::AbstractVector{N}) where {N<:Number} = sum(c * w for (c, w) in zip(C, W)) / sum(C)
XW2W̅(X::AbstractVector{N}, W::AbstractVector{N}) where {N<:Number} = sum(x * w for (x, w) in zip(X, W))
YW2W̅(Y::AbstractVector{N}, W::AbstractVector{N}) where {N<:Number} = sum(y / w for (y, w) in zip(Y, W)) |> inv
CW2W̅(C::AbstractVector{N}, W::AbstractVector{N}) where {N<:Number} = sum(c * w for (c, w) in zip(C, W)) / sum(C)

W̅YW2Y!(Y::Vector{N}, X::Vector{N}, W::AbstractVector{N}, W̅::N) where {N<:Number} = map!((x, w) -> x * w / W̅, Y, X, W)
W̅YW2X!(X::Vector{N}, Y::Vector{N}, W::AbstractVector{N}, W̅::N) where {N<:Number} = map!((y, w) -> y * W̅ / w, X, Y, W)
W̅YW2Y!(Y::AbstractVector{N}, X::AbstractVector{N}, W::AbstractVector{N}, W̅::N) where {N<:Number} = map!((x, w) -> x * w / W̅, Y, X, W)
W̅YW2X!(X::AbstractVector{N}, Y::AbstractVector{N}, W::AbstractVector{N}, W̅::N) where {N<:Number} = map!((y, w) -> y * W̅ / w, X, Y, W)

function C2X!(X::Vector{N}, C::Vector{N}) where {N<:Number}
function C2X!(X::AbstractVector{N}, C::AbstractVector{N}) where {N<:Number}
∑C = sum(C)
map!(c -> c / ∑C, X, C)
return nothing
end

function CW2Y!(Y::Vector{N}, C::Vector{N}, W::AbstractVector{N}) where {N<:Number}
function CW2Y!(Y::AbstractVector{N}, C::AbstractVector{N}, W::AbstractVector{N}) where {N<:Number}
∑CW = sum(c * w for (c, w) in zip(C, W))
map!((c, w) -> c * w / ∑CW, Y, C, W)
return nothing
end

Cρ2Y!(Y::Vector{N}, C::Vector{N}, W::AbstractVector{N}, ρ::N) where {N<:Number} = map!((c, w) -> c * w / ρ, Y, C, W)
YWρ2C!(C::Vector{N}, Y::Vector{N}, W::AbstractVector{N}, ρ::N) where {N<:Number} = map!((y, w) -> y * ρ / w, C, Y, W)
XW̅ρ2C!(C::Vector{N}, X::Vector{N}, W̅::N, ρ::N) where {N<:Number} = map!(x -> x * ρ / W̅, C, X)
Cρ2Y!(Y::AbstractVector{N}, C::AbstractVector{N}, W::AbstractVector{N}, ρ::N) where {N<:Number} = map!((c, w) -> c * w / ρ, Y, C, W)
YWρ2C!(C::AbstractVector{N}, Y::AbstractVector{N}, W::AbstractVector{N}, ρ::N) where {N<:Number} = map!((y, w) -> y * ρ / w, C, Y, W)
XW̅ρ2C!(C::AbstractVector{N}, X::AbstractVector{N}, W̅::N, ρ::N) where {N<:Number} = map!(x -> x * ρ / W̅, C, X)

CT2P(C::Vector{N}, T::N) where {N<:Number} = sum(C) * R * T
CT2P(C::AbstractVector{N}, T::N) where {N<:Number} = sum(C) * R * T
ρW̅T2P(ρ::N, W̅::N, T::N) where {N<:Number} = ρ * R * T / W̅

CW2ρ(C::Vector{N}, W::AbstractVector{N}) where {N<:Number} = sum(c * w for (c, w) in zip(C, W))
CW2ρ(C::AbstractVector{N}, W::AbstractVector{N}) where {N<:Number} = sum(c * w for (c, w) in zip(C, W))
PW̅T2ρ(P::N, W̅::N, T::N) where {N<:Number} = P * W̅ / (R * T)

struct Gas{N<:Number, S<:AbstractSpecies{N}, R<:AbstractReaction{N}}
Expand All @@ -66,7 +66,7 @@ end

Gas(mechanism::Mechanism{N}) where {N<:Number} = Gas(mechanism, State(mechanism))

function TPY!(gas::Gas{N}, T::N, P::N, Y::AbstractVector{N}) where {N<:Number}
function TPY!(gas::Gas{N}; T::N=temperature(gas), P::N=pressure(gas), Y::AbstractVector{N}=mass_fractions(gas)) where {N<:Number}
state = gas.state
W = molecular_weights(gas)
state.T = T
Expand All @@ -75,11 +75,11 @@ function TPY!(gas::Gas{N}, T::N, P::N, Y::AbstractVector{N}) where {N<:Number}
W̅ = YW2W̅(state.Y, W)
W̅YW2X!(state.X, state.Y, W, W̅)
state.ρ = PW̅T2ρ(P, W̅, T)
XW̅ρ2C!(state.C, state.X, W̅, state.ρ)
XW̅ρ2C!(state.C, state.X, W̅, state.ρ)
return gas
end

function TPX!(gas::Gas{N}, T::N, P::N, X::AbstractVector{N}) where {N<:Number}
function TPX!(gas::Gas{N}; T::N=temperature(gas), P::N=pressure(gas), X::AbstractVector{N}=molar_fractions(gas)) where {N<:Number}
state = gas.state
W = molecular_weights(gas)
state.T = T
Expand All @@ -88,11 +88,11 @@ function TPX!(gas::Gas{N}, T::N, P::N, X::AbstractVector{N}) where {N<:Number}
W̅ = XW2W̅(state.X, W)
W̅YW2Y!(state.Y, state.X, W, W̅)
state.ρ = PW̅T2ρ(P, W̅, T)
XW̅ρ2C!(state.C, X, W̅, state.ρ)
XW̅ρ2C!(state.C, X, W̅, state.ρ)
return gas
end

function TPC!(gas::Gas{N}, T::N, P::N, C::AbstractVector{N}) where {N<:Number}
function TPC!(gas::Gas{N}; T::N=temperature(gas), P::N=pressure(gas), C::AbstractVector{N}=molar_concentrations(gas)) where {N<:Number}
state = gas.state
W = molecular_weights(gas)
state.T = T
Expand All @@ -104,7 +104,7 @@ function TPC!(gas::Gas{N}, T::N, P::N, C::AbstractVector{N}) where {N<:Number}
return gas
end

function TρY!(gas::Gas{N}, T::N, ρ::N, Y::AbstractVector{N}) where {N<:Number}
function TρY!(gas::Gas{N}; T::N=temperature(gas), ρ::N=density(gas), Y::AbstractVector{N}=mass_fractions(gas)) where {N<:Number}
state = gas.state
W = molecular_weights(gas)
state.T = T
Expand All @@ -113,11 +113,11 @@ function TρY!(gas::Gas{N}, T::N, ρ::N, Y::AbstractVector{N}) where {N<:Number}
W̅ = YW2W̅(state.Y, W)
W̅YW2X!(state.X, state.Y, W, W̅)
state.P = ρW̅T2P(ρ, W̅, T)
XW̅ρ2C!(state.C, state.X, W̅, ρ)
XW̅ρ2C!(state.C, state.X, W̅, ρ)
return gas
end

function TρX!(gas::Gas{N}, T::N, ρ::N, X::AbstractVector{N}) where {N<:Number}
function TρX!(gas::Gas{N}; T::N=temperature(gas), ρ::N=density(gas), X::AbstractVector{N}=molar_fractions(gas)) where {N<:Number}
state = gas.state
W = molecular_weights(gas)
state.T = T
Expand All @@ -126,18 +126,18 @@ function TρX!(gas::Gas{N}, T::N, ρ::N, X::AbstractVector{N}) where {N<:Number}
W̅ = XW2W̅(state.X, W)
W̅YW2Y!(state.Y, state.X, W, W̅)
state.P = ρW̅T2P(ρ, W̅, T)
XW̅ρ2C!(state.C, X, W̅, state.ρ)
XW̅ρ2C!(state.C, X, W̅, state.ρ)
return gas
end

function TρC!(gas::Gas{N}, T::N, ρ::N, C::AbstractVector{N}) where {N<:Number}
function TρC!(gas::Gas{N}; T::N=temperature(gas), ρ::N=density(gas), C::AbstractVector{N}=molar_concentrations(gas)) where {N<:Number}
state = gas.state
W = molecular_weights(gas)
state.T = T
state.ρ = ρ
map!(identity, state.C, C)
state.P = CT2P(state.C, state.T)
C2X!(state.X, state.C)
Cρ2Y!(state.Y, state.C, W, state.ρ)
Cρ2Y!(state.Y, state.C, W, state.ρ)
return gas
end
30 changes: 17 additions & 13 deletions src/reactions.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@ struct Arrhenius{N<:Number} ## [A] := M^(∑νᵣ - 1) ⋅ s^-1; where M = cm^3
E::N
end

function Arrhenius(itr; order=1, mechunits = chemkin_default_units)
function Arrhenius(itr; order=1, mechunits=chemkin_default_units)
activation_energy_unit, length_unit, quantity_unit = (mechunits[k] for k in ("activation-energy", "length", "quantity"))
isempty(itr) && return nothing
A, β, E = itr
Expand All @@ -31,7 +31,7 @@ struct Plog{N<:Number}
arrhenius::Vector{Arrhenius{N}}
end

Plog(itr...; order = 1) = isempty(first(itr)) ? nothing : Plog(first(itr), Arrhenius{Float64}[Arrhenius(p; order) for p in zip(rest(itr, 2)...)])
Plog(itr...; order=1) = isempty(first(itr)) ? nothing : Plog(first(itr), Arrhenius{Float64}[Arrhenius(p; order) for p in zip(rest(itr, 2)...)])

function isin(p::N, x::Vector{N}) where {N<:Number}
for i in eachindex(x)
Expand All @@ -41,7 +41,7 @@ end

isout(p::N, x::Vector{N}) where {N<:Number} = p > last(x) ? lastindex(x) : p < first(x) ? firstindex(x) : nothing ## on the assumption that plogs are sorted, To-Do: sort in advance
interpolate(kᵢ::N, kᵢ₊₁::N, Pᵢ::N, Pᵢ₊₁::N, P::N) where {N<:Number} = log(kᵢ) + (log(kᵢ₊₁) - log(kᵢ)) * (log(P) - log(Pᵢ)) / (log(Pᵢ₊₁) - log(Pᵢ))
interpolate(::Val{:dP}, kᵢ::N, kᵢ₊₁::N, Pᵢ::N, Pᵢ₊₁::N, P::N) where {N<:Number} = (log(kᵢ₊₁) - log(kᵢ)) / (log(Pᵢ₊₁) - log(Pᵢ))P
interpolate(::Val{:dP}, kᵢ::N, kᵢ₊₁::N, Pᵢ::N, Pᵢ₊₁::N, P::N) where {N<:Number} = (log(kᵢ₊₁) - log(kᵢ)) / (log(Pᵢ₊₁) - log(Pᵢ))P

function ((; pressures, arrhenius)::Plog{N})(T::N, P::N) where {N<:Number}
i = isout(P, pressures)
Expand Down Expand Up @@ -145,7 +145,7 @@ function troe_function(::Val{:dC}, Fc::N, Pᵣ::N) where {N<:Number}
t₃ = t₂^2
t₄ = one(N) + (t₁ / t₃)

F = 10.0^(log10Fc / t₄)
F = 10.0^(log10Fc / t₄)
t₅ = log(10.0)F
t₆ = t₄^2
t₇ = 2.0log10Fc * t₁ * t₅
Expand Down Expand Up @@ -247,7 +247,7 @@ function forward_rate(v::Val{:dT}, reaction::FallOffReaction{N}, (; T, C)::State
else
Fc = troe_parameters(T)
dFcdT = troe_parameters(v, T)

F, dFdFc, dFdPᵣ = troe_function(v, Fc, Pᵣ)
dFdT = dFdFc * dFcdT + dFdPᵣ * dPᵣdT

Expand Down Expand Up @@ -287,10 +287,10 @@ function forward_rate(v::Val{:dP}, (; rates, forward_rate_parameters, plog_param
return nothing
end

change_enthalpy((; reactants, products)::AbstractReaction{<:Number}, i::Int = 1) =
change_enthalpy((; reactants, products)::AbstractReaction{<:Number}, i::Int=1) =
@inbounds sum(getfield(s.thermo.h, i)[] * ν for (s, ν) in flatten((reactants, products)))

change_entropy((; reactants, products)::AbstractReaction{<:Number}, i::Int = 1) =
change_entropy((; reactants, products)::AbstractReaction{<:Number}, i::Int=1) =
@inbounds sum(getfield(s.thermo.s, i)[] * ν for (s, ν) in flatten((reactants, products)))

function equilibrium_constants(reaction::AbstractReaction{N}, T::N) where {N<:Number}
Expand All @@ -307,20 +307,20 @@ end
function equilibrium_constants(::Val{:dT}, reaction::AbstractReaction{N}, T::N) where {N<:Number}
∆H, d∆HdT = (change_enthalpy(reaction, f) for f in OneTo(2))
∆S, d∆SdT = (change_entropy(reaction, f) for f in OneTo(2))

∑ν = reaction.reaction_order
t₀ = inv(R * T)
t₁ = ∆S * T - ∆H

Kp = exp(t₁ * t₀)
∂Kp∂∆S = Kp / R
∂Kp∂∆H = -Kp * t₀
∂Kp∂T = (∆S - t₁ / T) * Kp * t₀
∂Kp∂T = (∆S - t₁ / T) * Kp * t₀

Kc = Kp * (Pa * t₀)^∑ν
∂Kc∂Kp = (Pa * t₀)^∑ν
∂Kc∂T = -Kp * ∑ν * (Pa * t₀)^∑ν / T

dKcdT = ∂Kc∂T + (∂Kp∂T + ∂Kp∂∆H * d∆HdT + ∂Kp∂∆S * d∆SdT) * ∂Kc∂Kp
return Kc, dKcdT
end
Expand All @@ -331,6 +331,7 @@ function reverse_rate(reaction::AbstractReaction{N}, (; T)::State{N}) where {N<:
isnothing(reverse_rate_parameters) || (@inbounds rates.kr.val[] = reverse_rate_parameters(T);
return nothing
)

Kc = equilibrium_constants(reaction, T)
@inbounds rates.kr.val[] = rates.kf.val[] / Kc
return nothing
Expand All @@ -344,6 +345,7 @@ function reverse_rate(v::Val{:dT}, reaction::AbstractReaction{N}, (; T)::State{N
isnothing(reverse_rate_parameters) || ((kr.val[], kr.dT[]) = reverse_rate_parameters(v, T);
return nothing
)

Kc, dKcdT = equilibrium_constants(v, reaction, T)
t = inv(Kc)

Expand All @@ -363,6 +365,7 @@ function reverse_rate(::Val{:dP}, reaction::AbstractReaction{N}, (; T)::State{N}
isnothing(reverse_rate_parameters) || (@inbounds rates.kr.val[] = reverse_rate_parameters(T);
return nothing
)

Kc = equilibrium_constants(reaction, T)
@inbounds rates.kr.dP[] = rates.kf.dP[] / Kc
return nothing
Expand All @@ -376,6 +379,7 @@ function reverse_rate(::Val{:dC}, reaction::FallOffReaction{N}, (; T)::State{N})
isnothing(reverse_rate_parameters) || (@inbounds kr.val[] = reverse_rate_parameters(T);
return nothing
)

Kc = equilibrium_constants(reaction, T)
t = inv(Kc)

Expand Down Expand Up @@ -405,7 +409,7 @@ function progress_rate(::Val{:dT}, reaction::AbstractReaction{N}, (; C)::State{N
(; kf, kr, q) = reaction.rates
∏ᴵ = step(reaction.reactants, C)
∏ᴵᴵ = reaction.isreversible ? step(reaction.products, C) : zero(N)

M = reaction isa ThreeBodyReaction ? total_molar_concentration(C, reaction.enhancement_factors) : one(N)
@inbounds q.val[] = M * (kf.val[] * ∏ᴵ - kr.val[] * ∏ᴵᴵ)
@inbounds q.dT[] = M * (kf.dT[] * ∏ᴵ - kr.dT[] * ∏ᴵᴵ)
Expand Down Expand Up @@ -436,7 +440,7 @@ function progress_rate(::Val{:dC}, reaction::ThreeBodyReaction{N}, (; C)::State{
d∏ᴵᴵdCₖ = step(reaction.products, C, k)
@inbounds q.dC[k] = M * (kf.val[] * d∏ᴵdCₖ - kr.val[] * d∏ᴵᴵdCₖ)
end

for ((; k), dMdCₖ) in reaction.enhancement_factors
@inbounds q.dC[k] += dMdCₖ * (kf.val[] * ∏ᴵ - kr.val[] * ∏ᴵᴵ)
end
Expand All @@ -462,7 +466,7 @@ function progress_rate(::Val{:dP}, reaction::AbstractReaction{N}, (; C)::State{N
(; kf, kr, q) = reaction.rates
∏ᴵ = step(reaction.reactants, C)
∏ᴵᴵ = reaction.isreversible ? step(reaction.products, C) : zero(N)

M = reaction isa ThreeBodyReaction ? total_molar_concentration(C, reaction.enhancement_factors) : one(N)
@inbounds q.dP[] = M * (kf.dP[] * ∏ᴵ - kr.dP[] * ∏ᴵᴵ)
return nothing
Expand Down
Loading

0 comments on commit d79fe15

Please sign in to comment.