Merge pull request #176 from ReactionMechanismGenerator/modelreduction

Enabling simulation of reduced model

Author: mjohnson541

.github/workflows/documentation.yml

@@ -14,7 +14,7 @@ jobs:
       - uses: actions/checkout@v2
       - uses: julia-actions/setup-julia@latest
         with:
-          version: 1.5
+          version: 1.6
       - name: Install dependencies
         run: |
             wget https://repo.anaconda.com/miniconda/Miniconda2-latest-Linux-x86_64.sh -O miniconda.sh;

Project.toml

@@ -4,11 +4,9 @@ authors = ["Matt Johnson <mjohnson541@gmail.com>"]
 version = "0.4.0"
 
 [deps]
-CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 Calculus = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9"
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
 Conda = "8f4d0f93-b110-5947-807f-2305c1781a2d"
-DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 DiffEqSensitivity = "41bf760c-e81c-5289-8e54-58b1f1f8abe2"
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
@@ -21,47 +19,51 @@ LsqFit = "2fda8390-95c7-5789-9bda-21331edee243"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
+PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
 PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
 PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
 QuartzImageIO = "dca85d43-d64c-5e67-8c65-017450d5d020"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
+SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SmoothingSplines = "102930c3-cf33-599f-b3b1-9a29a5acab30"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
+Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 YAML = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [compat]
-CSV = "0.7,0.8"
 Calculus = "0.4,0.5"
 Colors = "0.11,0.12"
 Conda = "1"
-DiffEqBase = "6"
 DiffEqSensitivity = "6"
 ForwardDiff = "0.10"
-Images = "0.23"
+Images = "0.24"
 IncompleteLU = "0.2.0"
 IterTools = "1.3.0"
 LsqFit = "0.12"
-ModelingToolkit = "3"
-OrdinaryDiffEq = "5"
+ModelingToolkit = "8"
+OrdinaryDiffEq = "6"
 Parameters = "0.12"
+PreallocationTools = "0"
 PyCall = "1"
 PyPlot = "2"
 QuartzImageIO = "0.7"
 RecursiveArrayTools = "2.17"
 ReverseDiff = "1.9"
-SmoothingSplines = "0.2.1"
-SpecialFunctions = "0.10"
-StaticArrays = "0.12"
+SciMLBase = "1"
+SmoothingSplines = "0.3"
+SpecialFunctions = "1"
+StaticArrays = "1"
 Sundials = "4"
+Symbolics = "4"
 Tracker = "0.2"
 Unitful = "1.3.0"
 YAML = "0.4"

iJulia/Quasi-Steady State Assumptions and Isomer Lumping.ipynb

@@ -5,7 +5,7 @@ Tools for model edge analysis for automatic mechanism generation
 using Logging
 using Sundials
 using SparseArrays
-using DiffEqBase: build_solution
+using SciMLBase: build_solution
 using Base.Iterators: flatten
 
 abstract type AbstractTerminationCriterion end
@@ -76,7 +76,7 @@ Calculate key flux and concentration related quantities for edge analysis
 @inline function calcfluxes(sim::Simulation)
     t = sim.sol.t[end]
     dydt = zeros(sim.domain.indexes[end])
-    ns,cs,T,P,V,C,N,mu,kfs,krevs,Hs,Us,Gs,diffs,Cvave,cpdivR,phi = calcthermo(sim.domain,sim.sol.u[end],sim.sol.t[end],DiffEqBase.NullParameters())
+    ns,cs,T,P,V,C,N,mu,kfs,krevs,Hs,Us,Gs,diffs,Cvave,cpdivR,phi = calcthermo(sim.domain,sim.sol.u[end],sim.sol.t[end],SciMLBase.NullParameters())
     rts,frts,rrts = addreactionratecontributionsforwardreverse!(dydt,sim.domain.rxnarray,cs,kfs,krevs,V)
     calcdomainderivatives!(sim.domain,dydt,[];t=t,T=T,P=P,Us=Us,Hs=Hs,V=V,C=C,ns=ns,N=N,Cvave=Cvave)
     return dydt,rts,frts,rrts,cs
@@ -94,7 +94,7 @@ Calculate key flux and concentration related quantities for edge analysis
     domains = [sim.domain for sim in ssys.sims]
     interfaces = ssys.interfaces
     domain = domains[1]
-    ns,cs,T,P,V,C,N,mu,kfs,krevs,Hs,Us,Gs,diffs,Cvave,cpdivR,phi = calcthermo(domain,y,t,DiffEqBase.NullParameters())
+    ns,cs,T,P,V,C,N,mu,kfs,krevs,Hs,Us,Gs,diffs,Cvave,cpdivR,phi = calcthermo(domain,y,t,SciMLBase.NullParameters())
     vns = Array{Any,1}(undef,length(domains))
     vns[1] = ns
     vcs = Array{Any,1}(undef,length(domains))
@@ -136,7 +136,7 @@ Calculate key flux and concentration related quantities for edge analysis
     rrtsall = [rrts]
     for (i,domain) in enumerate(@views domains[2:end])
         k = i + 1
-        vns[k],vcs[k],vT[k],vP[k],vV[k],vC[k],vN[k],vmu[k],vkfs[k],vkrevs[k],vHs[k],vUs[k],vGs[k],vdiffs[k],vCvave[k],vcpdivR[k],vphi[k] = calcthermo(domain,y,t,DiffEqBase.NullParameters())
+        vns[k],vcs[k],vT[k],vP[k],vV[k],vC[k],vN[k],vmu[k],vkfs[k],vkrevs[k],vHs[k],vUs[k],vGs[k],vdiffs[k],vCvave[k],vcpdivR[k],vphi[k] = calcthermo(domain,y,t,SciMLBase.NullParameters())
         cstot[domain.indexes[1]:domain.indexes[2]] .= vcs[k]
         rts,frts,rrts = addreactionratecontributionsforwardreverse!(dydt,domain.rxnarray,cstot,vkfs[k],vkrevs[k],vV[k])
         push!(rtsall,rts)

src/Domain.jl

@@ -1,5 +1,5 @@
 using LinearAlgebra
-using DiffEqBase
+using SciMLBase
 using LsqFit
 
 abstract type AbstractInterface end
@@ -54,7 +54,7 @@ function getkfskrevs(ri::ReactiveInternalInterface,T1,T2,phi1,phi2,Gs1,Gs2,cstot
     return kfs,krevs
 end
 
-function evaluate(ri::ReactiveInternalInterface,dydt,domains,T1,T2,phi1,phi2,Gs1,Gs2,cstot,p::W) where {W<:DiffEqBase.NullParameters}
+function evaluate(ri::ReactiveInternalInterface,dydt,domains,T1,T2,phi1,phi2,Gs1,Gs2,cstot,p::W) where {W<:SciMLBase.NullParameters}
     kfs,krevs = getkfskrevs(ri,T1,T2,phi1,phi2,Gs1,Gs2,cstot)
     addreactionratecontributions!(dydt,ri.rxnarray,cstot,kfs,krevs,ri.A)
 end
@@ -108,7 +108,7 @@ function getkfskrevs(ri::ReactiveInternalInterfaceConstantTPhi,T1,T2,phi1,phi2,G
     return ri.kfs,ri.krevs
 end
 
-function evaluate(ri::ReactiveInternalInterfaceConstantTPhi,dydt,domains,T1,T2,phi1,phi2,Gs1,Gs2,cstot,p::W) where {W<:DiffEqBase.NullParameters}
+function evaluate(ri::ReactiveInternalInterfaceConstantTPhi,dydt,domains,T1,T2,phi1,phi2,Gs1,Gs2,cstot,p::W) where {W<:SciMLBase.NullParameters}
     addreactionratecontributions!(dydt,ri.rxnarray,cstot,ri.kfs,ri.krevs,ri.A)
 end
 
@@ -237,11 +237,25 @@ struct Inlet{Q<:Real,S,V<:AbstractArray,U<:Real,X<:Real,FF<:Function} <: Abstrac
 end
 
 function Inlet(domain::V,conddict::Dict{X1,X},F::FF) where {V,X1,X,B<:Real,FF<:Function}
+    T = 0.0
+    P = 0.0
+
     y = makespcsvector(domain.phase,conddict)
-    T = conddict["T"]
-    P = conddict["P"]
     yout = y./sum(y)
-    H = dot(getEnthalpy.(getfield.(domain.phase.species,:thermo),T),yout)
+
+    if haskey(conddict,"T")
+        T = conddict["T"]
+    end
+    if haskey(conddict,"P")
+        P = conddict["P"]
+    end
+
+    if haskey(conddict,"Hin")
+        H = conddict["Hin"]
+    else
+        @assert T != 0.0 && P != 0.0 "T and P need to be provided if Hin is not provided for Inlet"
+        H = dot(getEnthalpy.(getfield.(domain.phase.species,:thermo),T),yout)
+    end
     return Inlet(domain,yout,F,T,P,H)
 end
 

src/EdgeAnalysis.jl

@@ -0,0 +1,88 @@
+using Symbolics
+using ModelingToolkit
+
+struct ReducedModelMappings{F<:Function}
+    qssindexes::Array{Int64,1}
+    lumpedindexes::Array{Int64,1}
+    reducedindexes::Array{Int64,1}
+    lumpedgroupmapping::Array{Dict{Int64,Float64},1}
+    qssc!::F
+end
+
+export ReducedModelMappings
+
+function getmapping(phase::T,qssnames::A1,isomergroups::A2) where {T<:AbstractPhase,A1<:AbstractArray,A2<:AbstractArray}
+    
+    spcnames = getfield.(phase.species,:name)
+    qssindexes = [index for (index, name) in enumerate(spcnames) if name in qssnames]
+    
+    lumpedgroupmapping = [Dict{Int64,Float64}() for group in isomergroups]
+    lumpedindexes = Array{Int64,1}()
+    for (i,group) in enumerate(isomergroups)
+        for (spcname, weight) in group
+            index = findfirst(x->x==spcname, spcnames)
+            lumpedgroupmapping[i][index] = weight
+            push!(lumpedindexes,index)
+        end
+    end
+    reducedindexes = [index for index in 1:length(spcnames) if !(index in qssindexes) && !(index in lumpedindexes)]
+    return qssindexes,lumpedindexes,reducedindexes,lumpedgroupmapping
+end
+export getmapping
+
+function generateqsscmapping(phase::T,qssnames::A1,isomergroups::A2;saveqssc::Bool=false,outputname::String="qssc.jl") where {T<:AbstractPhase,A1<:AbstractArray,A2<:AbstractArray}
+
+    qssindexes,lumpedindexes,reducedindexes,lumpedgroupmapping = getmapping(phase, qssnames, isomergroups)
+
+    # initialize symbolic variables
+    @variables symdc[1:length(phase.species)] symc[1:length(phase.species)] symkf[1:length(phase.reactions)] symkrev[1:length(phase.reactions)] symV
+    
+    # collect them so that we can use setindex
+    symdc = collect(symdc)
+    symc = collect(symc)
+
+    # setting the derivatives to 0
+    symdc .= Num(0);
+
+    # fill in the symbolic expression for derivatives
+    addreactionratecontributions!(symdc,phase.rxnarray,symc,symkf,symkrev)
+    
+    # sett the derivatives to 0 for QSS species
+    eqs = 0 .~ symdc[qssindexes];
+
+    # equate the quadratic nonlinear terms x^2 to 0
+    nonlinearterms = Dict([])
+    for i in symc[qssindexes]
+        nonlinearterms[i*i]=0 #x^2=>0
+    end
+
+    # ideally we would use the following to filter out both x^2 and x*y cases where x, y ∈ symc[qssindexes]
+    # But current SymbolicUtils doesn't substitute x*y in k*x*y 
+    # nonlinearterms = Dict([])
+    # for (i,j) in Iterators.product(symc[qssindexes], symc[qssindexes])
+    #     nonlinearterms[i*j]=0.0 #x^2=>0
+    # end
+
+    # a workaround to deal with x*y nonlinear terms, where x, y ∈ symc[qssindexes]
+    for ind in 1:length(eqs)
+        if !(Symbolics.linear_expansion([eqs[ind]],symc[qssindexes])[3])
+            eqs[ind] = 0~substitute(eqs[ind].rhs,nonlinearterms) #substitute x^2 nonlinear terms as 0
+            filter!(x->length(get_variables(x.first,symc[qssindexes])) < 2,eqs[ind].rhs.dict) #filter out all the terms in the format of k*a*b that doesn't contain more than 1 variables in symc[qssindexes]
+            eqs[ind].rhs.sorted_args_cache[] = nothing #remove cached arguments
+            @assert Symbolics.linear_expansion([eqs[ind]],symc[qssindexes])[3] == true #assert the equations are all linear with respect to the QSS species
+        end
+    end
+    
+    qsscsolved = Symbolics.expand(Symbolics.solve_for(eqs,symc[qssindexes];simplify=false)) #solve for the expression of QSS species
+    
+    symqssc! = build_function(qsscsolved, symc, symkf, symkrev)[2] #convert the symbolic expression to a function
+
+    if saveqssc
+        write(outputname,string(symqssc!)); #save the function to file for future use
+    end
+    
+    return ReducedModelMappings(qssindexes,lumpedindexes,reducedindexes,lumpedgroupmapping,eval(symqssc!))
+    
+end
+export generateqsscmapping
+

src/Interface.jl

@@ -37,6 +37,8 @@ function makespcsvector(phase,spcdict)
             continue
         elseif key == "V"
             continue
+        elseif key == "Hin"
+            continue
         else
             ind = findfirst(isequal(key),spnames)
             @assert typeof(ind)<: Integer  "$key not found in species list: $spnames"

src/ModelReduction.jl

@@ -53,6 +53,7 @@ module ReactionMechanismSimulator
     include("Domain.jl")
     include("yml.jl")
     include("Parse.jl")
+    include("ModelReduction.jl")
     include("Reactor.jl")
     include("Simulation.jl")
     include("TransitorySensitivities.jl")