Skip to content

Commit b9fc0a0

Browse files
mjohnson541mjohnson541
authored
Merge pull request #92 from ReactionMechanismGenerator/analytic_Jp
Analytic parameter jacobian
2 parents 5a78065 + eb8b8dc commit b9fc0a0

File tree

3 files changed

+368
-19
lines changed

3 files changed

+368
-19
lines changed

‎src/Domain.jl‎

Lines changed: 300 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1988,6 +1988,306 @@ end
19881988

19891989
export jacobiany!
19901990

1991+
@inline function jacobianpnsderiv!(jacp::Q,y::U,p::W,t::Z,domain::D,rxnarray::Array{Int64,2},cs::Array{Float64,1},T::Float64,V::Float64,kfs::Array{Float64,1},krevs::Array{Float64,1},Nspcs::Int64,Nrxns::Int64) where {Q3<:AbstractArray,Q<:AbstractArray,U<:AbstractArray,W,Z<:Real,D<:Union{ConstantTPDomain,ConstantTADomain}}
1992+
@fastmath RTinv = 1.0/(R*T)
1993+
@simd for rxnind = 1:Nrxns
1994+
if @inbounds rxnarray[2,rxnind] == 0
1995+
@inbounds fderiv = cs[rxnarray[1,rxnind]]
1996+
elseif @inbounds rxnarray[3,rxnind] == 0
1997+
@fastmath @inbounds fderiv = cs[rxnarray[1,rxnind]]*cs[rxnarray[2,rxnind]]
1998+
else
1999+
@fastmath @inbounds fderiv = cs[rxnarray[1,rxnind]]*cs[rxnarray[2,rxnind]]*cs[rxnarray[3,rxnind]]
2000+
end
2001+
2002+
if @inbounds rxnarray[5,rxnind] == 0
2003+
@fastmath @inbounds rderiv = krevs[rxnind]/kfs[rxnind]*cs[rxnarray[4,rxnind]]
2004+
elseif @inbounds rxnarray[6,rxnind] == 0
2005+
@fastmath @inbounds rderiv = krevs[rxnind]/kfs[rxnind]*cs[rxnarray[4,rxnind]]*cs[rxnarray[5,rxnind]]
2006+
else
2007+
@fastmath @inbounds rderiv = krevs[rxnind]/kfs[rxnind]*cs[rxnarray[4,rxnind]]*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]]
2008+
end
2009+
2010+
flux = fderiv-rderiv
2011+
_spreadreactantpartials!(jacp,flux,rxnarray,rxnind,Nspcs+rxnind)
2012+
_spreadproductpartials!(jacp,-flux,rxnarray,rxnind,Nspcs+rxnind)
2013+
2014+
@fastmath @inbounds gderiv = rderiv*kfs[rxnind]*RTinv
2015+
2016+
@inbounds jacp[rxnarray[1,rxnind],rxnarray[1,rxnind]] -= gderiv
2017+
@inbounds _spreadreactantpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[1,rxnind])
2018+
if @inbounds rxnarray[2,rxnind] !== 0
2019+
@inbounds jacp[rxnarray[2,rxnind],rxnarray[1,rxnind]] -= gderiv
2020+
@inbounds jacp[rxnarray[1,rxnind],rxnarray[2,rxnind]] -= gderiv
2021+
@inbounds jacp[rxnarray[2,rxnind],rxnarray[2,rxnind]] -= gderiv
2022+
@inbounds _spreadreactantpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[2,rxnind])
2023+
if @inbounds rxnarray[3,rxnind] !== 0
2024+
@inbounds jacp[rxnarray[3,rxnind],rxnarray[1,rxnind]] -= gderiv
2025+
@inbounds jacp[rxnarray[3,rxnind],rxnarray[2,rxnind]] -= gderiv
2026+
@inbounds jacp[rxnarray[1,rxnind],rxnarray[3,rxnind]] -= gderiv
2027+
@inbounds jacp[rxnarray[2,rxnind],rxnarray[3,rxnind]] -= gderiv
2028+
@inbounds jacp[rxnarray[3,rxnind],rxnarray[3,rxnind]] -= gderiv
2029+
@inbounds _spreadreactantpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[3,rxnind])
2030+
end
2031+
end
2032+
2033+
@inbounds jacp[rxnarray[4,rxnind],rxnarray[4,rxnind]] -= gderiv
2034+
@inbounds _spreadproductpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[4,rxnind])
2035+
if @inbounds rxnarray[5,rxnind] !== 0
2036+
@inbounds jacp[rxnarray[5,rxnind],rxnarray[4,rxnind]] -= gderiv
2037+
@inbounds jacp[rxnarray[4,rxnind],rxnarray[5,rxnind]] -= gderiv
2038+
@inbounds jacp[rxnarray[5,rxnind],rxnarray[5,rxnind]] -= gderiv
2039+
@inbounds _spreadproductpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[5,rxnind])
2040+
if @inbounds rxnarray[6,rxnind] !== 0
2041+
@inbounds jacp[rxnarray[6,rxnind],rxnarray[4,rxnind]] -= gderiv
2042+
@inbounds jacp[rxnarray[6,rxnind],rxnarray[5,rxnind]] -= gderiv
2043+
@inbounds jacp[rxnarray[4,rxnind],rxnarray[6,rxnind]] -= gderiv
2044+
@inbounds jacp[rxnarray[5,rxnind],rxnarray[6,rxnind]] -= gderiv
2045+
@inbounds jacp[rxnarray[6,rxnind],rxnarray[6,rxnind]] -= gderiv
2046+
@inbounds _spreadproductpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[6,rxnind])
2047+
end
2048+
end
2049+
2050+
end
2051+
jacp .*= V
2052+
end
2053+
2054+
@inline function jacobianpnsderiv!(jacp::Q,y::U,p::W,t::Z,domain::D,rxnarray::Array{Int64,2},cs::Array{Float64,1},T::Float64,V::Float64,kfs::Array{Float64,1},krevs::Array{Float64,1},Nspcs::Int64,Nrxns::Int64) where {Q3<:AbstractArray,Q<:AbstractArray,U<:AbstractArray,W,Z<:Real,D<:Union{ConstantVDomain,ConstantPDomain,ParametrizedTPDomain,ParametrizedVDomain,ParametrizedPDomain,ParametrizedTConstantVDomain}}
2055+
@fastmath RTinv = 1.0/(R*T)
2056+
@simd for rxnind = 1:Nrxns
2057+
if @inbounds rxnarray[2,rxnind] == 0
2058+
@fastmath @inbounds fderiv = kfs[rxnind]*cs[rxnarray[1,rxnind]]
2059+
elseif @inbounds rxnarray[3,rxnind] == 0
2060+
@fastmath @inbounds fderiv = kfs[rxnind]*cs[rxnarray[1,rxnind]]*cs[rxnarray[2,rxnind]]
2061+
else
2062+
@fastmath @inbounds fderiv = kfs[rxnind]*cs[rxnarray[1,rxnind]]*cs[rxnarray[2,rxnind]]*cs[rxnarray[3,rxnind]]
2063+
end
2064+
2065+
if @inbounds rxnarray[5,rxnind] == 0
2066+
@fastmath @inbounds rderiv = krevs[rxnind]*cs[rxnarray[4,rxnind]]
2067+
elseif @inbounds rxnarray[6,rxnind] == 0
2068+
@fastmath @inbounds rderiv = krevs[rxnind]*cs[rxnarray[4,rxnind]]*cs[rxnarray[5,rxnind]]
2069+
else
2070+
@fastmath @inbounds rderiv = krevs[rxnind]*cs[rxnarray[4,rxnind]]*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]]
2071+
end
2072+
2073+
@fastmath flux = fderiv-rderiv
2074+
_spreadreactantpartials!(jacp,flux,rxnarray,rxnind,Nspcs+rxnind)
2075+
_spreadproductpartials!(jacp,-flux,rxnarray,rxnind,Nspcs+rxnind)
2076+
2077+
@fastmath gderiv = rderiv*RTinv
2078+
2079+
@inbounds jacp[rxnarray[1,rxnind],rxnarray[1,rxnind]] -= gderiv
2080+
@inbounds _spreadreactantpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[1,rxnind])
2081+
if @inbounds rxnarray[2,rxnind] !== 0
2082+
@inbounds jacp[rxnarray[2,rxnind],rxnarray[1,rxnind]] -= gderiv
2083+
@inbounds jacp[rxnarray[1,rxnind],rxnarray[2,rxnind]] -= gderiv
2084+
@inbounds jacp[rxnarray[2,rxnind],rxnarray[2,rxnind]] -= gderiv
2085+
@inbounds _spreadreactantpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[2,rxnind])
2086+
if @inbounds rxnarray[3,rxnind] !== 0
2087+
@inbounds jacp[rxnarray[3,rxnind],rxnarray[1,rxnind]] -= gderiv
2088+
@inbounds jacp[rxnarray[3,rxnind],rxnarray[2,rxnind]] -= gderiv
2089+
@inbounds jacp[rxnarray[1,rxnind],rxnarray[3,rxnind]] -= gderiv
2090+
@inbounds jacp[rxnarray[2,rxnind],rxnarray[3,rxnind]] -= gderiv
2091+
@inbounds jacp[rxnarray[3,rxnind],rxnarray[3,rxnind]] -= gderiv
2092+
@inbounds _spreadreactantpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[3,rxnind])
2093+
end
2094+
end
2095+
2096+
@inbounds jacp[rxnarray[4,rxnind],rxnarray[4,rxnind]] -= gderiv
2097+
@inbounds _spreadproductpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[4,rxnind])
2098+
if @inbounds rxnarray[5,rxnind] !== 0
2099+
@inbounds jacp[rxnarray[5,rxnind],rxnarray[4,rxnind]] -= gderiv
2100+
@inbounds jacp[rxnarray[4,rxnind],rxnarray[5,rxnind]] -= gderiv
2101+
@inbounds jacp[rxnarray[5,rxnind],rxnarray[5,rxnind]] -= gderiv
2102+
@inbounds _spreadproductpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[5,rxnind])
2103+
if @inbounds rxnarray[6,rxnind] !== 0
2104+
@inbounds jacp[rxnarray[6,rxnind],rxnarray[4,rxnind]] -= gderiv
2105+
@inbounds jacp[rxnarray[6,rxnind],rxnarray[5,rxnind]] -= gderiv
2106+
@inbounds jacp[rxnarray[4,rxnind],rxnarray[6,rxnind]] -= gderiv
2107+
@inbounds jacp[rxnarray[5,rxnind],rxnarray[6,rxnind]] -= gderiv
2108+
@inbounds jacp[rxnarray[6,rxnind],rxnarray[6,rxnind]] -= gderiv
2109+
@inbounds _spreadproductpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[6,rxnind])
2110+
end
2111+
end
2112+
end
2113+
jacp .*= V
2114+
end
2115+
2116+
@inline function jacobianpnsderiv!(jacp::Q,y::U,p::W,t::Z,domain::D,rxnarray::Array{Int64,2},cs::Array{Float64,1},T::Float64,V::Float64,kfs::Array{Float64,1},krevs::Array{Float64,1},Nspcs::Int64,Nrxns::Int64) where {Q3<:AbstractArray,Q<:AbstractArray,U<:AbstractArray,W,Z<:Real,D<:Union{ConstantTVDomain}}
2117+
@fastmath RTinv = 1.0/(R*T)
2118+
@simd for rxnind = 1:Nrxns
2119+
if @inbounds rxnarray[2,rxnind] == 0
2120+
@fastmath @inbounds fderiv = kfs[rxnind]*kfs[rxnind]/(p[Nspcs+rxnind]*p[Nspcs+rxnind])*cs[rxnarray[1,rxnind]]
2121+
elseif @inbounds rxnarray[3,rxnind] == 0
2122+
@fastmath @inbounds fderiv = kfs[rxnind]*kfs[rxnind]/(p[Nspcs+rxnind]*p[Nspcs+rxnind])*cs[rxnarray[1,rxnind]]*cs[rxnarray[2,rxnind]]
2123+
else
2124+
@fastmath @inbounds fderiv = kfs[rxnind]*kfs[rxnind]/(p[Nspcs+rxnind]*p[Nspcs+rxnind])*cs[rxnarray[1,rxnind]]*cs[rxnarray[2,rxnind]]*cs[rxnarray[3,rxnind]]
2125+
end
2126+
2127+
if @inbounds rxnarray[5,rxnind] == 0
2128+
@fastmath @inbounds rderiv = kfs[rxnind]*krevs[rxnind]/(p[Nspcs+rxnind]*p[Nspcs+rxnind])*cs[rxnarray[4,rxnind]]
2129+
elseif @inbounds rxnarray[6,rxnind] == 0
2130+
@fastmath @inbounds rderiv = kfs[rxnind]*krevs[rxnind]/(p[Nspcs+rxnind]*p[Nspcs+rxnind])*cs[rxnarray[4,rxnind]]*cs[rxnarray[5,rxnind]]
2131+
else
2132+
@fastmath @inbounds rderiv = kfs[rxnind]*krevs[rxnind]/(p[Nspcs+rxnind]*p[Nspcs+rxnind])*cs[rxnarray[4,rxnind]]*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]]
2133+
end
2134+
2135+
@fastmath flux = fderiv-rderiv
2136+
_spreadreactantpartials!(jacp,flux,rxnarray,rxnind,Nspcs+rxnind)
2137+
_spreadproductpartials!(jacp,-flux,rxnarray,rxnind,Nspcs+rxnind)
2138+
2139+
@fastmath @inbounds gderiv = rderiv*(p[Nspcs+rxnind]*p[Nspcs+rxnind])/kfs[rxnind]*RTinv
2140+
2141+
@inbounds jacp[rxnarray[1,rxnind],rxnarray[1,rxnind]] -= gderiv
2142+
@inbounds _spreadreactantpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[1,rxnind])
2143+
if @inbounds rxnarray[2,rxnind] !== 0
2144+
@inbounds jacp[rxnarray[2,rxnind],rxnarray[1,rxnind]] -= gderiv
2145+
@inbounds jacp[rxnarray[1,rxnind],rxnarray[2,rxnind]] -= gderiv
2146+
@inbounds jacp[rxnarray[2,rxnind],rxnarray[2,rxnind]] -= gderiv
2147+
@inbounds _spreadreactantpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[2,rxnind])
2148+
if @inbounds rxnarray[3,rxnind] !== 0
2149+
@inbounds jacp[rxnarray[3,rxnind],rxnarray[1,rxnind]] -= gderiv
2150+
@inbounds jacp[rxnarray[3,rxnind],rxnarray[2,rxnind]] -= gderiv
2151+
@inbounds jacp[rxnarray[1,rxnind],rxnarray[3,rxnind]] -= gderiv
2152+
@inbounds jacp[rxnarray[2,rxnind],rxnarray[3,rxnind]] -= gderiv
2153+
@inbounds jacp[rxnarray[3,rxnind],rxnarray[3,rxnind]] -= gderiv
2154+
@inbounds _spreadreactantpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[3,rxnind])
2155+
end
2156+
end
2157+
2158+
@inbounds jacp[rxnarray[4,rxnind],rxnarray[4,rxnind]] -= gderiv
2159+
@inbounds _spreadproductpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[4,rxnind])
2160+
if @inbounds rxnarray[5,rxnind] !== 0
2161+
@inbounds jacp[rxnarray[5,rxnind],rxnarray[4,rxnind]] -= gderiv
2162+
@inbounds jacp[rxnarray[4,rxnind],rxnarray[5,rxnind]] -= gderiv
2163+
@inbounds jacp[rxnarray[5,rxnind],rxnarray[5,rxnind]] -= gderiv
2164+
@inbounds _spreadproductpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[5,rxnind])
2165+
if @inbounds rxnarray[6,rxnind] !== 0
2166+
@inbounds jacp[rxnarray[6,rxnind],rxnarray[4,rxnind]] -= gderiv
2167+
@inbounds jacp[rxnarray[6,rxnind],rxnarray[5,rxnind]] -= gderiv
2168+
@inbounds jacp[rxnarray[4,rxnind],rxnarray[6,rxnind]] -= gderiv
2169+
@inbounds jacp[rxnarray[5,rxnind],rxnarray[6,rxnind]] -= gderiv
2170+
@inbounds jacp[rxnarray[6,rxnind],rxnarray[6,rxnind]] -= gderiv
2171+
@inbounds _spreadproductpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[6,rxnind])
2172+
end
2173+
end
2174+
end
2175+
jacp .*= V
2176+
end
2177+
2178+
function jacobianp!(jacp::Q,y::U,p::W,t::Z,domain::D,interfaces::Q3,colorvec::Q2=nothing) where {Q3<:AbstractArray,Q2,Q<:AbstractArray,U<:AbstractArray,W,Z<:Real,D<:Union{ConstantTPDomain,ParametrizedTPDomain}}
2179+
jacp.=0.0
2180+
ns,cs,T,P,V,C,N,mu,kfs,krevs,Hs,Us,Gs,diffs,Cvave,cpdivR = calcthermo(domain,y,t,p)
2181+
2182+
Nspcs=length(cs)
2183+
Nrxns = size(domain.rxnarray)[2]
2184+
2185+
jacobianpnsderiv!(jacp,y,p,t,domain,domain.rxnarray,cs,T,V,kfs,krevs,Nspcs,Nrxns)
2186+
2187+
@simd for i in 1:length(p)
2188+
@views @inbounds @fastmath jacp[domain.indexes[3],i] = sum(jacp[domain.indexes[1]:domain.indexes[2],i])*R*T/P
2189+
end
2190+
2191+
@simd for ind in domain.constantspeciesinds
2192+
@inbounds jacp[ind,:] .= 0.
2193+
end
2194+
end
2195+
2196+
function jacobianp!(jacp::Q,y::U,p::W,t::Z,domain::D,interfaces::Q3,colorvec::Q2=nothing) where {Q3<:AbstractArray,Q2,Q<:AbstractArray,U<:AbstractArray,W,Z<:Real,D<:Union{ConstantVDomain,ParametrizedVDomain}}
2197+
jacp.=0.0
2198+
ns,cs,T,P,V,C,N,mu,kfs,krevs,Hs,Us,Gs,diffs,Cvave,cpdivR = calcthermo(domain,y,t,p)
2199+
2200+
Nspcs = length(cs)
2201+
Nrxns = size(domain.rxnarray)[2]
2202+
2203+
dydt = zeros(size(y))
2204+
addreactionratecontributions!(dydt,domain.rxnarray,cs,kfs,krevs)
2205+
dydt .*= V
2206+
2207+
jacobianpnsderiv!(jacp,y,p,t,domain,domain.rxnarray,cs,T,V,kfs,krevs,Nspcs,Nrxns)
2208+
2209+
@simd for i in 1:Nspcs
2210+
@views @fastmath @inbounds jacp[domain.indexes[3],i] = -(dot(Us,jacp[domain.indexes[1]:domain.indexes[2],i])+dydt[i])/(N*Cvave)
2211+
@views @fastmath @inbounds jacp[domain.indexes[4],i] = sum(jacp[domain.indexes[1]:domain.indexes[2],i])*R*T/V + P/T*jacp[domain.indexes[3],i]
2212+
end
2213+
2214+
@simd for i in Nspcs+1:Nspcs+Nrxns
2215+
@views @fastmath @inbounds jacp[domain.indexes[3],i] = -dot(Us,jacp[domain.indexes[1]:domain.indexes[2],i])/(N*Cvave)
2216+
@views @fastmath @inbounds jacp[domain.indexes[4],i] = sum(jacp[domain.indexes[1]:domain.indexes[2],i])*R*T/V + P/T*jacp[domain.indexes[3],i]
2217+
end
2218+
2219+
@simd for ind in domain.constantspeciesinds
2220+
@inbounds jacp[ind,:] .= 0.
2221+
end
2222+
2223+
@simd for inter in interfaces
2224+
if isa(inter,Inlet) && domain == inter.domain
2225+
flow = inter.F(t)
2226+
@simd for i in 1:Nspcs
2227+
ddGidTdt = flow*(-ns[i]/N)/(N*Cvave)
2228+
jacp[domain.indexes[3],i] += ddGidTdt
2229+
jacp[domain.indexes[4],i] += P/T*ddGidTdt
2230+
end
2231+
end
2232+
end
2233+
end
2234+
2235+
function jacobianp!(jacp::Q,y::U,p::W,t::Z,domain::D,interfaces::Q3,colorvec::Q2=nothing) where {Q3<:AbstractArray,Q2,Q<:AbstractArray,U<:AbstractArray,W,Z<:Real,D<:Union{ConstantPDomain,ParametrizedPDomain}}
2236+
jacp.=0.0
2237+
ns,cs,T,P,V,C,N,mu,kfs,krevs,Hs,Us,Gs,diffs,Cvave,cpdivR = calcthermo(domain,y,t,p)
2238+
2239+
Nspcs = length(cs)
2240+
Nrxns = size(domain.rxnarray)[2]
2241+
2242+
dydt = zeros(size(y))
2243+
addreactionratecontributions!(dydt,domain.rxnarray,cs,kfs,krevs)
2244+
dydt .*= V
2245+
2246+
jacobianpnsderiv!(jacp,y,p,t,domain,domain.rxnarray,cs,T,V,kfs,krevs,Nspcs,Nrxns)
2247+
2248+
@fastmath Cpave = Cvave+R
2249+
@simd for i in 1:Nspcs
2250+
@views @fastmath @inbounds jacp[domain.indexes[3],i] = -(dot(Hs,jacp[domain.indexes[1]:domain.indexes[2],i])+dydt[i])/(N*Cpave) #divide by V to cancel ωV to ω
2251+
@views @fastmath @inbounds jacp[domain.indexes[4],i] = sum(jacp[domain.indexes[1]:domain.indexes[2],i])*R*T/P + jacp[domain.indexes[3],i]*V/T
2252+
end
2253+
2254+
@simd for i in Nspcs+1:Nspcs+Nrxns
2255+
@views @fastmath @inbounds jacp[domain.indexes[3],i] = -(dot(Hs,jacp[domain.indexes[1]:domain.indexes[2],i]))/(N*Cpave) #divide by V to cancel ωV to ω
2256+
@views @fastmath @inbounds jacp[domain.indexes[4],i] = sum(jacp[domain.indexes[1]:domain.indexes[2],i])*R*T/P + jacp[domain.indexes[3],i]*V/T
2257+
end
2258+
2259+
@simd for ind in domain.constantspeciesinds
2260+
@inbounds jacp[ind,:] .= 0.
2261+
end
2262+
2263+
@simd for inter in interfaces
2264+
if isa(inter,Inlet) && domain == inter.domain
2265+
flow = inter.F(t)
2266+
for i in 1:Nspcs
2267+
ddGidTdt = flow*(- ns[i]/N)/(N*Cpave)
2268+
jacp[domain.indexes[3],i] += ddGidTdt
2269+
jacp[domain.indexes[4],i] += ddGidTdt*V/T
2270+
end
2271+
end
2272+
end
2273+
end
2274+
2275+
function jacobianp!(jacp::Q,y::U,p::W,t::Z,domain::D,interfaces::Q3,colorvec::Q2=nothing) where {Q3<:AbstractArray,Q2,Q<:AbstractArray,U<:AbstractArray,W,Z<:Real,D<:Union{ConstantTVDomain,ParametrizedTConstantVDomain,ConstantTADomain}}
2276+
jacp.=0.0
2277+
ns,cs,T,P,V,C,N,mu,kfs,krevs,Hs,Us,Gs,diffs,Cvave,cpdivR = calcthermo(domain,y,t,p)
2278+
2279+
Nspcs = length(cs)
2280+
Nrxns = size(domain.rxnarray)[2]
2281+
2282+
jacobianpnsderiv!(jacp,y,p,t,domain,domain.rxnarray,cs,T,V,kfs,krevs,Nspcs,Nrxns)
2283+
2284+
@simd for ind in domain.constantspeciesinds
2285+
@inbounds jacp[ind,:] .= 0.
2286+
end
2287+
end
2288+
2289+
export jacobianp!
2290+
19912291
function getreactionindices(ig::Q) where {Q<:AbstractPhase}
19922292
arr = zeros(Int64,(6,length(ig.reactions)))
19932293
for (i,rxn) in enumerate(ig.reactions)

0 commit comments

Comments
 (0)