Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix symbolic bug #252

Closed
wants to merge 20 commits into from
Closed

Fix symbolic bug #252

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
20 commits
Select commit Hold shift + click to select a range
864abb9
fix dispatch for symbolic types
mjohnson541 Mar 6, 2024
8e117d3
add propane mechanism from paper to repository
mjohnson541 Mar 6, 2024
f422946
properly handle jac_prototype base on sparsity
mjohnson541 Mar 9, 2024
9019ee0
Update src/Reactor.jl
mjohnson541 Mar 18, 2024
8953444
remove zygote
mjohnson541 Mar 18, 2024
54f95f3
use new rmgmolecule builds
mjohnson541 Mar 18, 2024
da8fb5a
Formatting
hwpang Mar 19, 2024
e00aa22
use rmgmolecule from mjohnson541
mjohnson541 Mar 21, 2024
e68d462
add install test workflow
mjohnson541 Mar 21, 2024
fafa1f9
avoid Optim 1.9.3
mjohnson541 Mar 21, 2024
94e9aef
fix
mjohnson541 Mar 21, 2024
8bca5a5
fix
mjohnson541 Mar 21, 2024
060df8d
fix
mjohnson541 Mar 21, 2024
100808e
try including Enzyme
mjohnson541 Mar 21, 2024
443c736
add Enzyme
mjohnson541 Mar 21, 2024
3281070
specify LinearSolve
mjohnson541 Mar 21, 2024
0209f43
tweak
mjohnson541 Mar 21, 2024
3ba3ace
use newer julia and python in builds
mjohnson541 Mar 21, 2024
e60d482
use julia 1.10
mjohnson541 Mar 21, 2024
38cc087
switch to julia 1.9
mjohnson541 Mar 21, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
44 changes: 40 additions & 4 deletions src/PhaseState.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@
using Tracker
using ReverseDiff
using RecursiveArrayTools
using ModelingToolkit

@inline function calcgibbs(ph::U,T::W) where {U<:IdealPhase,W<:Real}
return getGibbs.(getfield.(ph.species,:thermo),T)
Expand Down Expand Up @@ -181,12 +182,12 @@
Calculates the forward and reverse rate coefficients for a given reaction, phase and state
Maintains diffusion limitations if the phase has diffusionlimited=true
"""
@inline function getkfkrev(rxn::ElementaryReaction,ph::U,T::W1,P::W2,C::W3,N::W4,ns::Q1,Gs::Q2,diffs::Q3,V::W5,phi::W8;kf::W6=-1.0,f::W7=-1.0) where {U<:AbstractPhase,W8,W6,W7,W5,W4,W1,W2,W3<:Real,Q1,Q2,Q3<:AbstractArray}
@inline function getkfkrev(rxn::ElementaryReaction,ph::U,T::W1,P::W2,C::W3,N::W4,ns::Vector{Q1},Gs::Q2,diffs,V::W5,phi::W8;kf::W6=-1.0,f=-1.0) where {U<:AbstractPhase,W8,W6,W5,W4,W1,W2,W3,Q1<:Real,Q2}
if signbit(kf)
if signbit(f)
kf = getkf(rxn,ph,T,P,C,ns,V,phi)
else
if isa(f,Num) || !signbit(f)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you run a formatting script over your changes or this file?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm having difficulty getting the formatter extension to connect with the julia kernel. Was there anything in particular you had in mind?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could you add blank after comma?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I ran formatter on the files changed in this PR, so I'm good with the formatting

kf = getkf(rxn,ph,T,P,C,ns,V,phi)*f
else
kf = getkf(rxn,ph,T,P,C,ns,V,phi)
end
end
Kc = getKc(rxn,ph,T,Gs,phi)
Expand Down Expand Up @@ -221,6 +222,41 @@
krev *= rxn.reversible
return (kf,krev)
end

@inline function getkfkrev(rxn::ElementaryReaction,ph::U,T::W1,P::W2,C::W3,N::W4,ns::Vector{Q1},Gs::Q2,diffs,V::W5,phi::W8;kf::W6=-1.0,f=-1.0) where {U<:AbstractPhase,W8,W6,W5,W4,W1,W2,W3,Q1<:Num,Q2}
kf = getkf(rxn,ph,T,P,C,ns,V,phi)*f
Kc = getKc(rxn,ph,T,Gs,phi)
@fastmath krev = kf/Kc
if ph.diffusionlimited
if length(rxn.reactants) == 1
if length(rxn.products) > 1
krevdiff = getDiffusiveRate(rxn.products,diffs)
@fastmath krev = krev*krevdiff/(krev+krevdiff)
@fastmath kf = Kc*krev

Check warning on line 235 in src/PhaseState.jl

View check run for this annotation

Codecov / codecov/patch

src/PhaseState.jl#L226-L235 

Added lines #L226 - L235 were not covered by tests
end
elseif length(rxn.products) == 1
kfdiff = getDiffusiveRate(rxn.reactants,diffs)
@fastmath kf = kf*kfdiff/(kf+kfdiff)
@fastmath krev = kf/Kc
elseif length(rxn.products) == length(rxn.reactants)
kfdiff = getDiffusiveRate(rxn.reactants,diffs)
krevdiff = getDiffusiveRate(rxn.products,diffs)
@fastmath kff = kf*kfdiff/(kf+kfdiff)
@fastmath krevr = krev*krevdiff/(krev+krevdiff)
@fastmath kfr = Kc*krevr
if kff > kfr
kf = kfr
krev = krevr

Check warning on line 249 in src/PhaseState.jl

View check run for this annotation

Codecov / codecov/patch

src/PhaseState.jl#L237-L249 

Added lines #L237 - L249 were not covered by tests
else
kf = kff
@fastmath krev = kf/Kc

Check warning on line 252 in src/PhaseState.jl

View check run for this annotation

Codecov / codecov/patch

src/PhaseState.jl#L251-L252 

Added lines #L251 - L252 were not covered by tests
end
end
end
kf *= rxn.forwardable
krev *= rxn.reversible
return (kf,krev)

Check warning on line 258 in src/PhaseState.jl

View check run for this annotation

Codecov / codecov/patch

src/PhaseState.jl#L256-L258 

Added lines #L256 - L258 were not covered by tests
end
export getkfkrev

@inline function getkfkrevs(phase::U,T::W1,P::W2,C::W3,N::W4,ns::Q1,Gs::Q2,diffs::Q3,V::W5,phi::W7;kfs::W6=nothing) where {U,W7,W6,W5<:Real,W1<:Real,W2<:Real,W3,W4,Q1<:AbstractArray,Q2,Q3<:AbstractArray}
Expand Down
39 changes: 28 additions & 11 deletions src/Reactor.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -54,14 +54,20 @@
prectmp = ilu(W, τ=tau)
preccache = Ref(prectmp)

if sparsity > 0.8
jac_prototype = J

Check warning on line 58 in src/Reactor.jl

View check run for this annotation

Codecov / codecov/patch

src/Reactor.jl#L58 

Added line #L58 was not covered by tests
else
jac_prototype = nothing
end

if (forwardsensitivities || !forwarddiff) && domain isa Union{ConstantTPDomain,ConstantVDomain,ConstantPDomain,ParametrizedTPDomain,ParametrizedVDomain,ParametrizedPDomain,ConstantTVDomain,ParametrizedTConstantVDomain,ConstantTAPhiDomain}
if !forwardsensitivities
odefcn = ODEFunction(dydt; jac=jacy!, paramjac=jacp!)
else
odefcn = ODEFunction(dydt; paramjac=jacp!)
end
else
odefcn = ODEFunction(dydt; jac=jacyforwarddiff!, paramjac=jacpforwarddiff!, jac_prototype=float.(J)) #jac_prototype is not needed/used for Sundials solvers but maybe needed for Julia solvers
odefcn = ODEFunction(dydt; jac=jacyforwarddiff!, paramjac=jacpforwarddiff!, jac_prototype=jac_prototype) #jac_prototype is not needed/used for Sundials solvers but maybe needed for Julia solvers

Check warning on line 70 in src/Reactor.jl

View check run for this annotation

Codecov / codecov/patch

src/Reactor.jl#L70 

Added line #L70 was not covered by tests
end
if forwardsensitivities
ode = ODEForwardSensitivityProblem(odefcn, y0, tspan, p)
Expand All @@ -78,9 +84,9 @@
sys = modelingtoolkitize(ode)
jac = eval(ModelingToolkit.generate_jacobian(sys)[2])
if (forwardsensitivities || !forwarddiff) && domain isa Union{ConstantTPDomain,ConstantVDomain,ConstantPDomain,ParametrizedTPDomain,ParametrizedVDomain,ParametrizedPDomain,ConstantTVDomain,ParametrizedTConstantVDomain,ConstantTAPhiDomain}
odefcn = ODEFunction(dydt; jac=jac, paramjac=jacp!)
odefcn = ODEFunction(dydt; jac=jac, paramjac=jacp!, jac_prototype=jac_prototype)

Check warning on line 87 in src/Reactor.jl

View check run for this annotation

Codecov / codecov/patch

src/Reactor.jl#L87 

Added line #L87 was not covered by tests
else
odefcn = ODEFunction(dydt; jac=jac, paramjac=jacpforwarddiff!)
odefcn = ODEFunction(dydt; jac=jac, paramjac=jacpforwarddiff!, jac_prototype=jac_prototype)

Check warning on line 89 in src/Reactor.jl

View check run for this annotation

Codecov / codecov/patch

src/Reactor.jl#L89 

Added line #L89 was not covered by tests
end
if forwardsensitivities
ode = ODEForwardSensitivityProblem(odefcn, y0, tspan, p)
Expand All @@ -90,7 +96,7 @@
end
return Reactor(domain, interfaces, y0, tspan, p, ode, recsolver, forwardsensitivities, forwarddiff, modelingtoolkit, tau, precsundials, psetupsundials, precsjulia)
end
function Reactor(domains::T, y0s::W1, tspan::W2, interfaces::Z=Tuple(), ps::X=SciMLBase.NullParameters(); forwardsensitivities=false, modelingtoolkit=false, tau=1e-3) where {T<:Tuple,W1<:Tuple,Z,X,W2}
function Reactor(domains::T, y0s::W1, tspan::W2, interfaces::Z=Tuple(), ps::X=SciMLBase.NullParameters(); forwardsensitivities=false, forwarddiff=false, modelingtoolkit=false, tau=1e-3) where {T<:Tuple,W1<:Tuple,Z,X,W2}
#adjust indexing
y0 = zeros(sum(length(y) for y in y0s))
Nvars = 0
Expand Down Expand Up @@ -217,18 +223,24 @@
prectmp = ilu(W, τ=tau)
preccache = Ref(prectmp)

if sparsity > 0.8
jac_prototype = J

Check warning on line 227 in src/Reactor.jl

View check run for this annotation

Codecov / codecov/patch

src/Reactor.jl#L227 

Added line #L227 was not covered by tests
else
jac_prototype = nothing
end

if forwardsensitivities
odefcn = ODEFunction(dydt; paramjac=jacp!)
odefcn = ODEFunction(dydt; paramjac=jacp!, jac_prototype=jac_prototype)

Check warning on line 233 in src/Reactor.jl

View check run for this annotation

Codecov / codecov/patch

src/Reactor.jl#L233 

Added line #L233 was not covered by tests
ode = ODEForwardSensitivityProblem(odefcn, y0, tspan, p)
recsolver = Sundials.CVODE_BDF(linear_solver=:GMRES)
if modelingtoolkit
sys = modelingtoolkitize(ode)
jac = eval(ModelingToolkit.generate_jacobian(sys)[2])
odefcn = ODEFunction(dydt; jac=jac, paramjac=jacp!)
odefcn = ODEFunction(dydt; jac=jac, paramjac=jacp!, jac_prototype=jac_prototype)

Check warning on line 239 in src/Reactor.jl

View check run for this annotation

Codecov / codecov/patch

src/Reactor.jl#L239 

Added line #L239 was not covered by tests
ode = ODEForwardSensitivityProblem(odefcn, y0, tspan, p)
end
else
odefcn = ODEFunction(dydt; jac=jacy!, paramjac=jacp!, jac_prototype=float.(J))
odefcn = ODEFunction(dydt; jac=jacy!, paramjac=jacp!, jac_prototype=jac_prototype)
ode = ODEProblem(odefcn, y0, tspan, p)
if sparsity > 0.8 #empirical threshold to use preconditioner
recsolver = Sundials.CVODE_BDF(linear_solver=:GMRES, prec=precsundials, psetup=psetupsundials, prec_side=1)
Expand All @@ -238,7 +250,7 @@
if modelingtoolkit
sys = modelingtoolkitize(ode)
jac = eval(ModelingToolkit.generate_jacobian(sys)[2])
odefcn = ODEFunction(dydt; jac=jac, paramjac=jacp!)
odefcn = ODEFunction(dydt; jac=jac, paramjac=jacp!, jac_prototype=jac_prototype)

Check warning on line 253 in src/Reactor.jl

View check run for this annotation

Codecov / codecov/patch

src/Reactor.jl#L253 

Added line #L253 was not covered by tests
ode = ODEProblem(odefcn, y0, tspan, p)
end
end
Expand All @@ -256,7 +268,6 @@
dydt(dy::X, y::T, p::V, t::Q) where {X,T,Q,V} = dydtreactor!(dy, y, t, domain, interfaces, reducedmodelmappings, reducedmodelcache, p=p)
jacy!(J::Q2, y::T, p::V, t::Q) where {Q2,T,Q,V} = jacobianyforwarddiff!(J, y, p, t, domain, interfaces, reducedmodelmappings, reducedmodelcache)
jacp!(J::Q2, y::T, p::V, t::Q) where {Q2,T,Q,V} = jacobianpforwarddiff!(J, y, p, t, domain, interfaces, reducedmodelmappings, reducedmodelcache)

#y0 in Y space
y0 = zeros(length(reducedmodelmappings.reducedindexes) + length(reducedmodelmappings.lumpedgroupmapping) + length(domain.thermovariabledict))
@inbounds @views y0[1:end-length(domain.thermovariabledict)-length(reducedmodelmappings.lumpedgroupmapping)] .= y0unlumped[reducedmodelmappings.reducedindexes]
Expand Down Expand Up @@ -294,7 +305,13 @@
prectmp = ilu(W, τ=tau)
preccache = Ref(prectmp)

odefcn = ODEFunction(dydt; jac=jacy!, paramjac=jacp!, jac_prototype=float.(J)) #jac_prototype is not needed/used for Sundials solvers but maybe needed for Julia solvers
if sparsity > 0.8
jac_prototype = J

Check warning on line 309 in src/Reactor.jl

View check run for this annotation

Codecov / codecov/patch

src/Reactor.jl#L309 

Added line #L309 was not covered by tests
else
jac_prototype = nothing
end

odefcn = ODEFunction(dydt; jac=jacy!, paramjac=jacp!, jac_prototype=jac_prototype) #jac_prototype is not needed/used for Sundials solvers but maybe needed for Julia solvers

if forwardsensitivities
ode = ODEForwardSensitivityProblem(odefcn, y0, tspan, p)
Expand All @@ -310,7 +327,7 @@
if modelingtoolkit
sys = modelingtoolkitize(ode)
jac = eval(ModelingToolkit.generate_jacobian(sys)[2])
odefcn = ODEFunction(dydt; jac=jac, paramjac=jacp!)
odefcn = ODEFunction(dydt; jac=jac, paramjac=jacp!, jac_prototype=jac_prototype)

Check warning on line 330 in src/Reactor.jl

View check run for this annotation

Codecov / codecov/patch

src/Reactor.jl#L330 

Added line #L330 was not covered by tests
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is modelingtoolkitize tested in our test? Maybe we should add a test

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We should, but we need a very small mechanism to make it fast.

if forwardsensitivities
ode = ODEForwardSensitivityProblem(odefcn, y0, tspan, p)
else
Expand Down
Loading
Loading