diff --git a/Project.toml b/Project.toml index 282aa67..c4e5e8a 100644 --- a/Project.toml +++ b/Project.toml @@ -1,23 +1,18 @@ name = "Apophis" uuid = "cf72b054-1e85-4718-a3f4-02cfdc10a054" authors = ["moataz-sabry ", "Ahmed Hassan "] -version = "1.3.3" +version = "1.3.4" [deps] -Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" -LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SplitApplyCombine = "03a91e81-4c3e-53e1-a0a4-9c0c8f19dd66" -ThreadsX = "ac1d9e8a-700a-412c-b207-f0111f4b6c0d" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" YAML = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [compat] -Combinatorics = "1.0.2" PyCall = "1.95.1" SplitApplyCombine = "1.2.2" -ThreadsX = "0.1.11" Unitful = "1.12.3" YAML = "0.4.8" Zygote = "0.6.55" diff --git a/src/Apophis.jl b/src/Apophis.jl index 6aea08f..cfbe089 100644 --- a/src/Apophis.jl +++ b/src/Apophis.jl @@ -2,14 +2,11 @@ module Apophis using Base.Iterators: filter, flatten, reverse, take, partition using Base: Fix1, Fix2, OneTo, rest -using Combinatorics -using LinearAlgebra using SparseArrays using SplitApplyCombine: combinedimsview, mapview using Unitful using YAML: load_file using Zygote -# using ThreadsX abstract type AbstractSpecies{N<:Number} end abstract type AbstractReaction{N<:Number} end @@ -22,7 +19,6 @@ const Tᵣ = 300.0 # K const d = 0.14 const Maybe{T} = Union{T, Nothing} -const Diffusions = (:binary_diffusions, :mean_diffusions) const Energies{N<:Number} = NamedTuple{(:val, :dT), Tuple{Vector{N}, Vector{N}}} const Thermodynamics{N<:Number} = NamedTuple{(:cₚ, :h, :s), NTuple{3, Energies{N}}} const Rates{N<:Number} = NamedTuple{(:val, :dT, :dC), NTuple{3, Vector{N}}} diff --git a/src/gas.jl b/src/gas.jl index 32a14f2..b7df334 100644 --- a/src/gas.jl +++ b/src/gas.jl @@ -51,15 +51,15 @@ end """ Main function to create a `Gas` object. """ -function Gas(name::Union{Nothing, Symbol,String} = nothing; - kinetics_path = "", - thermo_path = "", - transport_path = "", - as::Type{<:Number} = Float64, init...) - - mechanism = read_mechanism(name; kinetics_path = kinetics_path, thermo_path = thermo_path, transport_path = transport_path, as = as) +function Gas(name::Union{Nothing, Symbol, String} = nothing; + kinetics_path::Maybe{String} = nothing, + thermo_path::Maybe{String} = nothing, + transport_path::Maybe{String} = nothing, + as::Type{<:Number} = Float64, + init...) + + mechanism = read_mechanism(name; kinetics_path, thermo_path, transport_path, as) gas = Gas(mechanism, State(mechanism)) - set!(gas; init...) return gas end diff --git a/src/reactions.jl b/src/reactions.jl index 6dc776e..a1ffc4b 100644 --- a/src/reactions.jl +++ b/src/reactions.jl @@ -208,8 +208,9 @@ function forward_rate(reaction::FallOffReaction{N}, (; T, C)::State{N}) where {N k∞, kₒ = (high_pressure_parameters, low_pressure_parameters)(T) M = total_molar_concentration(C, enhancement_factors) - Pᵣ = reduced_pressure(kₒ, M, k∞) + iszero(M) && (@inbounds rates.kf.val[] = zero(N); return nothing) + Pᵣ = reduced_pressure(kₒ, M, k∞) F = isnothing(troe_parameters) ? one(N) : troe_parameters(T) |> Fix2(troe_function, Pᵣ) @inbounds rates.kf.val[] = F * k∞ * Pᵣ * inv(one(N) + Pᵣ) return nothing @@ -230,6 +231,7 @@ function forward_rate(v::Val{:dT}, reaction::FallOffReaction{N}, (; T, C)::State (k∞, dk∞dT), (kₒ, dkₒdT) = (high_pressure_parameters, low_pressure_parameters)(v, T) M = total_molar_concentration(C, enhancement_factors) + iszero(M) && (@inbounds rates.kf.val[] = rates.kf.dT[] = zero(N); return nothing) Pᵣ = reduced_pressure(kₒ, M, k∞) t = inv(one(N) + Pᵣ) @@ -250,7 +252,7 @@ function forward_rate(v::Val{:dT}, reaction::FallOffReaction{N}, (; T, C)::State dFdT = dFdFc * dFcdT + dFdPᵣ * dPᵣdT @inbounds rates.kf.val[] = F * k∞ * Pᵣ * t - @inbounds rates.kf.dT[] = t * ((dk∞dT * Pᵣ + dPᵣdT * k∞ * t)F + dFdT * k∞ * Pᵣ) + @inbounds rates.kf.dT[] = t * (F * (dk∞dT * Pᵣ + dPᵣdT * k∞ * t) + dFdT * k∞ * Pᵣ) end return nothing end @@ -263,6 +265,7 @@ function forward_rate(v::Val{:dC}, reaction::FallOffReaction{N}, (; T, C)::State k∞, kₒ = (high_pressure_parameters, low_pressure_parameters)(T) M = total_molar_concentration(C, enhancement_factors) + iszero(M) && (@inbounds rates.kf.dC .= rates.kf.val[] = zero(N); return nothing) Pᵣ = reduced_pressure(kₒ, M, k∞) t = inv(one(N) + Pᵣ) @@ -272,8 +275,8 @@ function forward_rate(v::Val{:dC}, reaction::FallOffReaction{N}, (; T, C)::State @inbounds rates.kf.val[] = F * k∞ * Pᵣ * t dkfdPᵣ = F * k∞ * t^2 - dkfdF = k∞ * Pᵣ * t - for ((; k), dMdCₖ) in enhancement_factors ## SparseArrays? + dkfdF = k∞ * Pᵣ * t + for ((; k), dMdCₖ) in enhancement_factors @inbounds rates.kf.dC[k] = (dkfdPᵣ + dkfdF * dFdPᵣ) * dPᵣdM * dMdCₖ end return nothing @@ -424,9 +427,9 @@ end function progress_rate(::Val{:dC}, reaction::ThreeBodyReaction{N}, (; C)::State{N}) where {N<:Number} (; kf, kr, q) = reaction.rates + M = total_molar_concentration(C, reaction.enhancement_factors) ∏ᴵ = step(reaction.reactants, C) ∏ᴵᴵ = reaction.isreversible ? step(reaction.products, C) : zero(N) - M = total_molar_concentration(C, reaction.enhancement_factors) for k in eachindex(q.dC) d∏ᴵdCₖ = step(reaction.reactants, C, k) diff --git a/src/reader.jl b/src/reader.jl index 16d46be..42dd624 100644 --- a/src/reader.jl +++ b/src/reader.jl @@ -41,7 +41,7 @@ const chemkin_default_units = Dict{String, String}( "quantity" => "mol" ) -without_spaces(text::AbstractString) = replace(text, r"\s+" => "") +remove_spaces(text::AbstractString) = replace(text, r"\s+" => "") #= Reads a file data located at `file_path` as String. Filters out any exclamation marks `!` and the text following it. =# @@ -66,7 +66,7 @@ find_composition(composition::Dict, N::Type{<:Number}) = Pair{Symbol, N}[element #= Finds molecular weight of the given species. =# -find_weight(composition::Vector{Pair{Symbol, N}}) where {N<:Number} = sum(weights[element] * count for (element, count) in composition) +find_weight(composition::Vector{Pair{Symbol, N}}) where {N<:Number} = sum(N(weights[element]) * count for (element, count) in composition) #= Finds temperature ranges of the given species. @@ -104,7 +104,7 @@ end function _find_species(k::Int, species::Symbol, thermo::T, transport::Maybe{T}, N::Type{<:Number}) where {T<:String} transport_parameters = isnothing(transport) ? nothing : find_transport_parameters(species, transport, N) thermodynamics_parameters = find_thermodynamic_parameters(species, thermo, N) - return Species(k, species, Pair{Reaction{N, Species{N}}, N}[], thermodynamics_parameters..., transport_parameters, (; zip(Diffusions, (zeros(N, 9), zero(N)))...), + return Species(k, species, Pair{Reaction{N, Species{N}}, N}[], thermodynamics_parameters..., transport_parameters, (; zip((:cₚ, :h, :s), Tuple((; zip((:val, :dT), (zeros(N, 1) for _ in OneTo(2)))...) for _ in OneTo(3)))...), (; zip((:ω̇,), tuple((; zip((:val, :dT, :dC), (zeros(N, 1), zeros(N, 1), zeros(N, 0)))...)))...)) end @@ -115,7 +115,7 @@ function _find_species(k::Int, species::Dict, N::Type{<:Number}) weight = find_weight(composition) thermodynamic_parameters = find_thermodynamic_parameters(species["thermo"], N) - return Species(k, name, Pair{Reaction{N, Species{N}}, N}[], composition, weight, thermodynamic_parameters, nothing, (; zip(Diffusions, (zeros(N, 9), zero(N)))...), + return Species(k, name, Pair{Reaction{N, Species{N}}, N}[], composition, weight, thermodynamic_parameters, nothing, (; zip((:cₚ, :h, :s), Tuple((; zip((:val, :dT), (zeros(N, 1) for _ in OneTo(2)))...) for _ in OneTo(3)))...), (; zip((:ω̇,), tuple((; zip((:val, :dT, :dC), (zeros(N, 1), zeros(N, 1), zeros(N, 0)))...)))...)) end @@ -135,7 +135,11 @@ end #= Finds index of the given `species` in a species list. =# -assign(list::Vector, item::Union{AbstractString, Symbol}) = list[findfirst(i -> isequal(getfield(i, 2), Symbol(item)), list)] +function assign(list::Vector, item::Union{AbstractString, Symbol}) + index = findfirst(i -> isequal(getfield(i, 2), Symbol(item)), list) + @assert !isnothing(index) "$item is not present in the specified mechanism. Please check the mechanism files and ensure that $item is included" + return list[index] +end #= Checks if the given `reaction` is reversible. Returns true if reversible, false otherwise. @@ -179,7 +183,7 @@ Decomposes the given `reaction` into pairs of species and number of moles, for b =# function decompose(reaction::AbstractString, species_list::Vector{Species{N}}) where {N<:Number} sides = split(reaction, r"\s*?\s*", keepempty=false) ## left and right side of the equation - sides_components = Tuple(without_spaces(side) |> side -> split(side, r"\+|\(\+|(?<=M)\)"i; keepempty=false) for side in sides) ## products or reactants + sides_components = Tuple(remove_spaces(side) |> side -> split(side, r"\+|\(\+|(? without_spaces |> uparse |> Fix1(ustrip, u"Pa") for p in eachindex(plogs)], [plogs[p] |> Fix2(get_arrhenius, order) for p in eachindex(plogs)]) +find_plog(plogs::Vector{<:Dict}, order::Int, N::Type{<:Number}) = Plog{N}([plogs[p]["P"] |> remove_spaces |> uparse |> Fix1(ustrip, u"Pa") for p in eachindex(plogs)], [plogs[p] |> Fix2(get_arrhenius, order) for p in eachindex(plogs)]) #= Finds the enhancement factors in the given reaction `data`. Returns vector of pairs with species index as keys and enhancement factors as values. @@ -343,10 +347,10 @@ function get_stoichiometry_matrix(species::Vector{<:Species}) ## fail in case a return sparse(K, I, V) end -function _read_mechanism(kinetics::T, thermo::T, transport::T, N::Type{<:Number}) where {T<:String} +function _read_mechanism(kinetics::T, thermo::T, transport::Maybe{T}, N::Type{<:Number}) where {T<:String} kinetics_data = read_data(kinetics) thermo_data = read_data(thermo) - trans_data = isfile(transport) ? read_data(transport) : nothing + trans_data = isnothing(transport) ? nothing : isfile(transport) ? read_data(transport) : nothing species = find_species(kinetics_data, thermo_data, trans_data, N) K = length(species) @@ -367,18 +371,25 @@ function _read_mechanism(mechanism::Dict, N::Type{<:Number}) return (; species, reactions, stoichiometry_matrix) end -function read_mechanism(name::Union{Nothing,Symbol,String} = nothing; kinetics_path = "", thermo_path = "", transport_path = "", as::Type{<:Number} = Float64) - name isa String && last(name, 5) == ".yaml" && return _read_mechanism(load_file(name), as) +function read_mechanism(name::Union{Nothing, Symbol, String} = nothing; + kinetics_path::Maybe{String} = nothing, + thermo_path::Maybe{String} = nothing, + transport_path::Maybe{String} = nothing, + as::Type{<:Number} = Float64) - if name isa Symbol - mechanism_path = pkgdir(Apophis, "test/mechanisms/$(name)") - kinetics = mechanism_path * "/kinetics.dat" - thermo = mechanism_path * "/thermo.dat" - transport = mechanism_path * "/transport.dat" - else + name isa String && endswith(name, r"\.(:?yaml|yml)") && return _read_mechanism(load_file(name), as) + + if name isa Symbol || name isa String + mechanism_path = pkgdir(Apophis, "test/mechanisms/$name") + kinetics = joinpath(mechanism_path, "kinetics.dat") + thermo = joinpath(mechanism_path, "thermo.dat") + transport = joinpath(mechanism_path, "transport.dat") + elseif !isnothing(kinetics_path) && !isnothing(thermo_path) kinetics = kinetics_path thermo = thermo_path - transport = transport_path + transport = isnothing(transport_path) ? nothing : transport_path + else + throw(ArgumentError("Invalid arguments")) end return _read_mechanism(kinetics, thermo, transport, as) end \ No newline at end of file diff --git a/src/sensitivity.jl b/src/sensitivity.jl index 8eb03c1..61fcc21 100644 --- a/src/sensitivity.jl +++ b/src/sensitivity.jl @@ -20,6 +20,8 @@ _dkfdA((; forward_rate_parameters)::Union{ElementaryReaction{N}, ThreeBodyReacti function _dkfdA((; high_pressure_parameters, low_pressure_parameters, enhancement_factors, troe_parameters)::FallOffReaction{N}, (; T, C)::State{N}) where {N<:Number} ## Add dkfdAₒ later k∞, kₒ = (high_pressure_parameters, low_pressure_parameters)(T) M = total_molar_concentration(C, enhancement_factors) + iszero(M) && return zero(N) + Pᵣ = reduced_pressure(kₒ, M, k∞) t = inv(one(N) + Pᵣ) diff --git a/src/species.jl b/src/species.jl index 5188b8a..51ac3ee 100644 --- a/src/species.jl +++ b/src/species.jl @@ -54,7 +54,6 @@ struct Species{N<:Number, R<:AbstractReaction{N}} <: AbstractSpecies{N} weight::N nasa_polynomial::NasaPolynomial{N} transport_parameters::Maybe{TransportParameters{N}} - diffusions::NamedTuple{Diffusions, Tuple{Vector{N}, N}} thermo::Thermodynamics{N} rates::SpeciesRates{N} end diff --git a/src/utils.jl b/src/utils.jl index eab87f2..d9b1df3 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -14,7 +14,7 @@ Base.show(io::IO, (; state)::Gas) = foreach(f -> println(io, "$f: $(getfield(sta function find_init((; mechanism, state)::Gas{N}, I) where {N<:Number} S = zero(state.X) - for (s, f) in eachmatch(regextionary[:inputs], I) + for (s, f) in eachmatch(:inputs, I) species = assign(mechanism.species, s) S[species.k] = parse(N, f) end @@ -37,7 +37,7 @@ function set!(gas::Gas{N}; init...) where {N<:Number} Xᵢ = init[:X] X = Xᵢ isa AbstractVector ? Xᵢ : find_init(gas, init[:X]) P = init[:P] |> (P -> P isa Quantity ? ustrip(u"Pa", P) : P) |> N - TPX!(gas, T, P, X) |> N + TPX!(gas, T, P, X) elseif haskey(init, :Y) && haskey(init, :ρ) Xᵢ = init[:X] X = Xᵢ isa AbstractVector ? Xᵢ : find_init(gas, init[:X])