Merge pull request #170 from ReactionMechanismGenerator/multidomainsens

Multidomain and Interfaces Adjoint Sensitivity Analysis

Author: mjohnson541

src/Domain.jl

@@ -1630,7 +1630,7 @@ end
             @simd for i in domain.indexes[1]:domain.indexes[2]
                 @inbounds @fastmath jac[i,i] -= flow/N
             end
-            @views @inbounds @fastmath jac[domain.indexes[1]:domain.indexes[2],domain.indexes[3]] .-= flow.*ns.*(-C)
+            @views @inbounds @fastmath jac[domain.indexes[1]:domain.indexes[2],domain.indexes[3]] .-= -flow.*ns/(N*V)
         end
     end
 end

src/Interface.jl

@@ -113,7 +113,7 @@ function evaluate(ri::ReactiveInternalInterfaceConstantTPhi,dydt,domains,T1,T2,p
 end
 
 function evaluate(ri::ReactiveInternalInterfaceConstantTPhi,dydt,domains,T1,T2,phi1,phi2,Gs1,Gs2,cstot,p)
-    if all(p[ri.parameterindexes[1]:ri.parameterindexes[2]] .== ri.kfs)
+    if p[ri.parameterindexes[1]:ri.parameterindexes[2]] == ri.kfs
         kfs = ri.kfs
     else 
         kfs = p[ri.parameterindexes[1]:ri.parameterindexes[2]]

src/Reactor.jl

@@ -69,6 +69,12 @@ function Reactor(domains::T,y0s::W,tspan::W2,interfaces::Z=Tuple(),ps::X=DiffEqB
                 domain.rxnarray[i,j] += k-1
             end
         end
+        for i in 1:length(domain.constantspeciesinds)
+            domain.constantspeciesinds[i] += k-1
+        end
+        for (thermovar,ind) in domain.thermovariabledict
+            domain.thermovariabledict[thermovar] += k-1
+        end
         domain.indexes[1] = k
         k += Nspcs
         domain.indexes[2] = k-1
@@ -241,7 +247,7 @@ export getrates
 end
 
 @inline function addreactionratecontributions!(dydt::Q,rarray::Array{W2,2},cs::W,kfs::Z,krevs::Y,V) where {Q,Z,Y,T,W,W2}
-    @inbounds for i = 1:size(rarray)[2]
+    @inbounds @simd for i = 1:size(rarray)[2]
         if @inbounds rarray[2,i] == 0
             @inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]
         elseif @inbounds rarray[3,i] == 0
@@ -326,7 +332,8 @@ export addreactionratecontributionsforwardreverse!
     return dydt
 end
 @inline function dydtreactor!(dydt::RC,y::U,t::Z,domains::Q,interfaces::B;p::RV=DiffEqBase.NullParameters(),sensitivity::Bool=true) where {RC,RV,B,Z,U,Q<:Tuple}    
-    cstot = zeros(typeof(y).parameters[1],length(y))
+    cstot = similar(y)
+    cstot .= 0.0
     dydt .= 0.0
     domain = domains[1]
     ns,cs,T,P,V,C,N,mu,kfs,krevs,Hs,Us,Gs,diffs,Cvave,cpdivR,phi = calcthermo(domain,y,t,p)

src/Simulation.jl

@@ -354,6 +354,9 @@ based alternative algorithm is slower, but avoids this concern.
 """
 function getadjointsensitivities(bsol::Q,target::String,solver::W;sensalg::W2=InterpolatingAdjoint(autojacvec=ReverseDiffVJP(false)),abstol::Float64=1e-6,reltol::Float64=1e-3,kwargs...) where {Q,W,W2}
     @assert target in bsol.names || target in ["T","V","P"]
+
+    pethane = 160
+
     if target in ["T","V","P"]
         if haskey(bsol.domain.thermovariabledict, target)
             ind = bsol.domain.thermovariabledict[target]
@@ -362,14 +365,16 @@ function getadjointsensitivities(bsol::Q,target::String,solver::W;sensalg::W2=In
         end
     else
         ind = findfirst(isequal(target),bsol.names)
-        sensdomain,sensspcnames,senstooriginspcind,senstooriginrxnind = getsensdomain(bsol.domain,ind)
-        if :thermovariabledict in fieldnames(typeof(bsol.domain))
-            yinds = vcat(senstooriginspcind,collect(values(bsol.domain.thermovariabledict)))
-        else
-            yinds = vcat(senstooriginspcind)
+        if isempty(bsol.interfaces)
+            sensdomain,sensspcnames,senstooriginspcind,senstooriginrxnind = getsensdomain(bsol.domain,ind)
+            if :thermovariabledict in fieldnames(typeof(bsol.domain))
+                yinds = vcat(senstooriginspcind,collect(values(bsol.domain.thermovariabledict)))
+            else
+                yinds = vcat(senstooriginspcind)
+            end
+            pinds = vcat(senstooriginspcind,length(bsol.domain.phase.species).+senstooriginrxnind)
+            ind = findfirst(isequal(target),sensspcnames)
         end
-        pinds = vcat(senstooriginspcind,length(bsol.domain.phase.species).+senstooriginrxnind)
-        ind = findfirst(isequal(target),sensspcnames)
     end
 
     function sensg(y::X,p::Array{Y,1},t::Z) where {Q,V,X,Y<:Float64,Z}
@@ -399,19 +404,19 @@ function getadjointsensitivities(bsol::Q,target::String,solver::W;sensalg::W2=In
 
     function g(y::X,p::Array{Y,1},t::Z) where {Q,V,X,Y<:Float64,Z} 
         dy = similar(y,length(y))
-        return dydtreactor!(dy,y,t,bsol.domain,[],p=p)[ind]
+        return dydtreactor!(dy,y,t,bsol.domain,bsol.interfaces,p=p)[ind]
     end
     function g(y::Array{X,1},p::Y,t::Z) where {Q,V,X<:Float64,Y,Z} 
         dy = similar(p,length(y))
-        return dydtreactor!(dy,y,t,bsol.domain,[],p=p)[ind]
+        return dydtreactor!(dy,y,t,bsol.domain,bsol.interfaces,p=p)[ind]
     end
     function g(y::Array{X,1},p::Array{Y,1},t::Z) where {Q,V,X<:Float64,Y<:Float64,Z} 
         dy = zeros(length(y))
-        return dydtreactor!(dy,y,t,bsol.domain,[],p=p)[ind]
+        return dydtreactor!(dy,y,t,bsol.domain,bsol.interfaces,p=p)[ind]
     end
     function g(y::Array{X,1},p::Array{Y,1},t::Z) where {Q,V,X<:ForwardDiff.Dual,Y<:ForwardDiff.Dual,Z} 
         dy = similar(y,length(y))
-        return dydtreactor!(dy,y,t,bsol.domain,[],p=p)[ind]
+        return dydtreactor!(dy,y,t,bsol.domain,bsol.interfaces,p=p)[ind]
     end
 
     dsensgdu(out, y, p, t) = ForwardDiff.gradient!(out, y -> sensg(y, p, t), y)
@@ -423,23 +428,64 @@ function getadjointsensitivities(bsol::Q,target::String,solver::W;sensalg::W2=In
     dgdurevdiff(out, y, p, t) = ReverseDiff.gradient!(out, y -> g(y, p, t), y)
     dgdprevdiff(out, y, p, t) = ReverseDiff.gradient!(out, p -> g(y, p, t), p)
 
-    pethane = 160
     if length(bsol.domain.p)<= pethane
-        if target in ["T","V","P"]
+        if target in ["T","V","P"] || !isempty(bsol.interfaces)
             du0,dpadj = adjoint_sensitivities(bsol.sol,solver,g,nothing,(dgdu,dgdp);sensealg=sensalg,abstol=abstol,reltol=reltol,kwargs...)
         else
             du0,dpadj = adjoint_sensitivities(bsol.sol,solver,sensg,nothing,(dsensgdu,dsensgdp);sensealg=sensalg,abstol=abstol,reltol=reltol,kwargs...)
         end
     else
-        if target in ["T","V","P"]
+        if target in ["T","V","P"] || !isempty(bsol.interfaces)
             du0,dpadj = adjoint_sensitivities(bsol.sol,solver,g,nothing,(dgdurevdiff,dgdprevdiff);sensealg=sensalg,abstol=abstol,reltol=reltol,kwargs...)
         else
             du0,dpadj = adjoint_sensitivities(bsol.sol,solver,sensg,nothing,(dsensgdurevdiff,dsensgdprevdiff);sensealg=sensalg,abstol=abstol,reltol=reltol,kwargs...)
         end
     end
     dpadj[length(bsol.domain.phase.species)+1:end] .*= bsol.domain.p[length(bsol.domain.phase.species)+1:end]
-    if !(target in ["T","V","P"])
+    if !(target in ["T","V","P"]) && isempty(bsol.interfaces)
         dpadj ./= bsol.sol(bsol.sol.t[end])[senstooriginspcind[ind]]
+    elseif !(target in ["T","V","P"]) && !isempty(bsol.interfaces)
+        dpadj ./= bsol.sol(bsol.sol.t[end])[ind]
+    end
+    return dpadj
+end
+
+function getadjointsensitivities(syssim::Q,bsol::W3,target::String,solver::W;sensalg::W2=InterpolatingAdjoint(autojacvec=ReverseDiffVJP(false)),abstol::Float64=1e-6,reltol::Float64=1e-3,kwargs...) where {Q,W,W2,W3}
+    @assert target in bsol.names || target in ["T","V","P"]
+    if target in ["T","V","P"]
+        if haskey(bsol.domain.thermovariabledict, target)
+            ind = bsol.domain.thermovariabledict[target]
+        else
+            throw(error("$(bsol.domain) doesn't have $target in its thermovariables"))
+        end
+    else
+        ind = findfirst(isequal(target),bsol.names)+bsol.domain.indexes[1]-1
+    end
+    domains = Tuple([x.domain for x in syssim.sims])
+    function g(y::X,p::Array{Y,1},t::Z) where {Q,V,X,Y<:Float64,Z} 
+        dy = similar(y,length(y))
+        return dydtreactor!(dy,y,t,domains,syssim.interfaces,p=p)[ind]
+    end
+    function g(y::Array{X,1},p::Y,t::Z) where {Q,V,X<:Float64,Y,Z} 
+        dy = similar(p,length(y))
+        return dydtreactor!(dy,y,t,domains,syssim.interfaces,p=p)[ind]
+    end
+    function g(y::Array{Float64,1},p::Array{Float64,1},t::Z) where {Q,V,Z} 
+        dy = similar(p,length(y))
+        return dydtreactor!(dy,y,t,domains,syssim.interfaces,p=p)[ind]
+    end
+    function g(y::Array{X,1},p::Array{Y,1},t::Z) where {Q,V,X<:ForwardDiff.Dual,Y<:ForwardDiff.Dual,Z} 
+        dy = similar(y,length(y))
+        return dydtreactor!(dy,y,t,domains,syssim.interfaces,p=p)[ind]
+    end
+    dgdu(out, y, p, t) = ForwardDiff.gradient!(out, y -> g(y, p, t), y)
+    dgdp(out, y, p, t) = ForwardDiff.gradient!(out, p -> g(y, p, t), p)
+    du0,dpadj = adjoint_sensitivities(syssim.sol,solver,g,nothing,(dgdu,dgdp);sensealg=sensalg,abstol=abstol,reltol=reltol,kwargs...)
+    for domain in domains
+       dpadj[domain.parameterindexes[1]+length(domain.phase.species):domain.parameterindexes[2]] .*= syssim.p[domain.parameterindexes[1]+length(domain.phase.species):domain.parameterindexes[2]]
+    end
+    if !(target in ["T","V","P"])
+       dpadj ./= bsol.sol(bsol.sol.t[end])[ind]
     end
     return dpadj
 end

src/TestReactors.jl

@@ -124,6 +124,46 @@ rerr = [isinf(x) ? 0.0 : x for x in rerr]
 @test all((abs.(rerr) .> 1e-1).==false)
 end;
 
+#Constant T and P Ideal Gas
+@testset "Test constant T and P reactor with interfaces simulation" begin
+    #Define the phase (how species thermodynamic and kinetic properties calculated)
+   initialconds = Dict(["T"=>1000.0,"P"=>1e5,"H2"=>0.67,"O2"=>0.33]) #Set simulation Initial Temp and Pressure
+   domain,y0,p = ConstantTPDomain(phase=ig,initialconds=initialconds) #Define the domain (encodes how system thermodynamic properties calculated)
+
+   interfaces = [Inlet(domain,Dict{String,Float64}("H2"=>0.67,"O2"=>0.33,"T"=>1000.0,"P"=>1e5),x->0.001),
+                Outlet(domain,x->0.001)]
+   
+   react = Reactor(domain,y0,(0.0,150.11094),interfaces;p=p) #Create the reactor object
+   sol = solve(react.ode,CVODE_BDF(),abstol=1e-20,reltol=1e-12); #solve the ode associated with the reactor
+   sim = Simulation(sol,domain,interfaces);
+   
+   #analytic jacobian vs. ForwardDiff jacobian
+   t = 20.44002454;
+   y = sol(t)
+   ja = jacobiany(y,p,t,domain,interfaces,nothing);
+   j = jacobianyforwarddiff(y,p,t,domain,interfaces,nothing);
+   @test all((abs.(ja.-j) .> 1e-4.*abs.(j).+1e-16).==false)
+   
+   jpa = jacobianp(y,p,t,domain,interfaces,nothing);
+   jp = jacobianpforwarddiff(y,p,t,domain,interfaces,nothing);
+   @test all((abs.(jpa.-jp) .> 1e-4.*abs.(jp).+1e-16).==false)
+   
+   #sensitivities
+   dps = getadjointsensitivities(sim,"H2",CVODE_BDF(linear_solver=:GMRES);sensealg=InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true)),abstol=1e-16,reltol=1e-6)
+   react2 = Reactor(domain,y0,(0.0,150.11094),interfaces;p=p,forwardsensitivities=true)
+   sol2 = solve(react2.ode,CVODE_BDF(linear_solver=:GMRES),abstol=1e-21,reltol=1e-7); #solve the ode associated with the reactor
+   sim2 = Simulation(sol2,domain,interfaces)
+   
+   x,dp = extract_local_sensitivities(sol2,150.11094);
+   ind = findfirst(isequal("H2"),sim2.names)
+   dpvs = [v[ind] for v in dp]
+   dpvs[length(domain.phase.species)+1:end] .*= domain.p[length(domain.phase.species)+1:end]
+   dpvs ./= sol2(150.11094)[ind]
+   rerr = (dpvs .- dps')./dpvs
+   rerr = [isinf(x) ? 0.0 : x for x in rerr]
+   @test all((abs.(rerr) .> 1e-1).==false)
+   end;
+
 #Constant V adiabatic Ideal Gas
 #uses superminimal.yml mechanism
 @testset "Constant volume adiabatic reactor simulation" begin
@@ -266,6 +306,37 @@ end;
 
 
 
+@testset "Multi-domain ConstantV sensitivity analysis" begin
+    phaseDict = readinput("../src/testing/superminimal.rms")
+    spcs = phaseDict["phase"]["Species"]
+    rxns = phaseDict["phase"]["Reactions"]
+    ig = IdealGas(spcs,rxns,name="phase")
+    
+    initialcondsV = Dict(["T"=>1000.0,"P"=>10.0e5,"H2"=>0.67,"O2"=>0.33]) 
+    domainV,y0V,pV = ConstantVDomain(phase=ig,initialconds=initialcondsV) #Define the domain (encodes how system thermodynamic properties calculated)
+    
+    reactV = Reactor(domainV,y0V,(0.0,0.037);p=pV) #Create the reactor object
+    solV = solve(reactV.ode,CVODE_BDF(),abstol=1e-16,reltol=1e-6); #solve the ode associated with the reactor
+    simV = Simulation(solV,domainV)
+    
+    initialcondsV1 = Dict(["T"=>1000.0,"P"=>10.0e5,"H2"=>0.67,"O2"=>0.33]) 
+    domainV1,y0V1,pV1 = ConstantVDomain(phase=ig,initialconds=initialcondsV1) #Define the domain (encodes how system thermodynamic properties calculated)
+    initialcondsV2 = Dict(["T"=>1000.0,"P"=>10.0e5,"H2"=>0.67,"O2"=>0.33]) 
+    domainV2,y0V2,pV2 = ConstantVDomain(phase=ig,initialconds=initialcondsV2) #Define the domain (encodes how system thermodynamic properties calculated)
+    
+    react,y0,p = Reactor((domainV1,domainV2),(y0V1,y0V2),(0.0,0.037),[],(pV1,pV2));
+    sol = solve(react.ode,CVODE_BDF(),abstol=1e-16,reltol=1e-6);
+    sysim = SystemSimulation(sol,(domainV1,domainV2),[],p);
+    
+    t = 0.03
+    @test sol(t)[1:length(spcs)] ≈ solV(t)[1:end-2] rtol=1e-5
+    @test sol(t)[length(spcs)+1:end-4] ≈ solV(t)[1:end-2] rtol=1e-5
+    
+    dpsV = getadjointsensitivities(simV,"H2",CVODE_BDF();sensealg=InterpolatingAdjoint(autojacvec=ReverseDiffVJP(false)),abstol=1e-16,reltol=1e-6)
+    dps = getadjointsensitivities(sysim,sysim.sims[1],"H2",CVODE_BDF();sensealg=InterpolatingAdjoint(autojacvec=ReverseDiffVJP(false)),abstol=1e-16,reltol=1e-6)
+    @test dpsV ≈ dps rtol=1e-4
+end;
+
 @testset "Multi-domain ConstantV and ConstantTP simulation" begin
     phaseDict = readinput("../src/testing/superminimal.rms")
     spcs = phaseDict["phase"]["Species"]