Add test for ConstantTPDomain reactor with interfaces

Author: hwpang

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