diff --git a/src/Interface.jl b/src/Interface.jl index 99ab65d2..38283e2e 100644 --- a/src/Interface.jl +++ b/src/Interface.jl @@ -89,7 +89,7 @@ struct ReactiveInternalInterfaceConstantTPhi{J,N,B,B2,B3,C,C2,Q<:AbstractReactio forwardability::Array{Bool,1} end function ReactiveInternalInterfaceConstantTPhi(domain1,domain2,reactions,T,A,phi=0.0) - @assert domain1.T == domain2.T + @assert domain1.T == domain2.T reactions = upgradekinetics(reactions,domain1,domain2) rxnarray = getinterfacereactioninds(domain1,domain2,reactions) M,Nrp1,Nrp2 = getstoichmatrix(domain1,domain2,reactions) @@ -301,7 +301,7 @@ export ConstantReservoirDiffusion kLAkHCondensationEvaporationWithReservoir adds evaporation and condensation to (1) a liquid phase domain with a constant composition vapor resevoir, where number of moles, P, and T need to be specified, or (2) a gas phase domain with a constant composition liquid resevoir, where number of moles, V, and T need to be specified. -kLA and kH are used to model cond/evap. +kLA and kH are used to model cond/evap. kLA is liquid volumetric mass transfer coefficient with unit 1/s , and kH is Henry's law constant defined as gas phase partial pressure of solute over liquid phase concentration of solute. """ @@ -389,7 +389,7 @@ end export VolumetricFlowRateInlet """ -VolumeMaintainingOutlet is designed for gas phase domain such that the flow rate of this outlet will adjust to maintain the volume of the +VolumeMaintainingOutlet is designed for gas phase domain such that the flow rate of this outlet will adjust to maintain the volume of the domain to be constant. This is particularly useful to simulate any vapor-liquid phase system where the gas phase outlet is determined by the amount of evaporation. """ @@ -451,6 +451,58 @@ function getinterfaceignoremasstransferspcinds(domaingas, domainliq, ignoremasst return indices end +""" +FickDiffusionInternalInterface models Fickian diffusive mass transfer across a stagnant +film of finite thickness between two finite volume domains. Unlike +ConstantReservoirDiffusion, both sides carry their own state, so species +accumulate/deplete over time and mass is conserved across the interface. + +The molar flux into domain1 (and equal-and-opposite out of domain2) for species i is + J_i = A * D_i * (c2_i - c1_i) / layer_thickness [mol/s] +where c1, c2 are the per-species concentrations of domain1, domain2, A is the +interfacial area, layer_thickness is the film spacing, and D_i is the per-species +diffusivity across the film. +""" +struct FickDiffusionInternalInterface{D1,D2} <: AbstractInternalInterface + domain1::D1 + domain2::D2 + A::Float64 + layer_thickness::Float64 + diffusivity::Array{Float64,1} + parameterindexes::Array{Int64,1} + domaininds::Array{Int64,1} + p::Array{Float64,1} +end + +function FickDiffusionInternalInterface(domain1, domain2, A, layer_thickness; + diffusivity=nothing) + @assert getfield.(domain1.phase.species, :name) == getfield.(domain2.phase.species, :name) "FickDiffusionInternalInterface requires both domains to share the same ordered species list" + nspc = length(domain1.phase.species) + D = diffusivity === nothing ? convert(Array{Float64,1}, domain1.diffusivity) : convert(Array{Float64,1}, diffusivity) + @assert length(D) == nspc "diffusivity must have one entry per species ($nspc); got $(length(D)). If domain1.diffusivity is empty, pass diffusivity explicitly." + return FickDiffusionInternalInterface(domain1, domain2, Float64(A), Float64(layer_thickness), + D, [1, nspc], [0, 0], ones(nspc)), ones(nspc) +end +export FickDiffusionInternalInterface + +@inline function fickdiffusionflux(fd::FickDiffusionInternalInterface, cstot) + @views @inbounds c1 = cstot[fd.domain1.indexes[1]:fd.domain1.indexes[2]] + @views @inbounds c2 = cstot[fd.domain2.indexes[1]:fd.domain2.indexes[2]] + @fastmath return fd.A .* fd.diffusivity .* (c2 .- c1) ./ fd.layer_thickness +end + +function evaluate(fd::FickDiffusionInternalInterface, dydt, cstot, p::W) where {W<:SciMLBase.NullParameters} + flux = fickdiffusionflux(fd, cstot) + @views @inbounds @fastmath dydt[fd.domain1.indexes[1]:fd.domain1.indexes[2]] .+= flux + @views @inbounds @fastmath dydt[fd.domain2.indexes[1]:fd.domain2.indexes[2]] .-= flux +end + +function evaluate(fd::FickDiffusionInternalInterface, dydt, cstot, p) + flux = fickdiffusionflux(fd, cstot) .* p[fd.parameterindexes[1]:fd.parameterindexes[2]] + @views @inbounds @fastmath dydt[fd.domain1.indexes[1]:fd.domain1.indexes[2]] .+= flux + @views @inbounds @fastmath dydt[fd.domain2.indexes[1]:fd.domain2.indexes[2]] .-= flux +end + struct FragmentBasedReactiveFilmGrowthInterfaceConstantT{D1,D2,Q<:AbstractReaction,M1} <: AbstractReactiveInternalInterface domainfilm::D1 domain2::D2 @@ -567,4 +619,4 @@ function evaluate(ri::FragmentBasedReactiveFilmGrowthInterfaceConstantT, dydt, V end end -export evaluate \ No newline at end of file +export evaluate diff --git a/src/Reactor.jl b/src/Reactor.jl index bb4460b6..c89b475c 100644 --- a/src/Reactor.jl +++ b/src/Reactor.jl @@ -196,6 +196,12 @@ function Reactor(domains::T, y0s::W1, tspan::W2, interfaces::Z=Tuple(), ps::X=Sc inter.parameterindexes[2] = length(p) + length(ps[i+length(domains)]) inter.ignoremasstransferspcinds .= getinterfaceignoremasstransferspcinds(inter.domaingas, inter.domainliq, inter.ignoremasstransferspcnames) p = vcat(p, ps[i+length(domains)]) + elseif isa(inter, FickDiffusionInternalInterface) + inter.domaininds[1] = findfirst(isequal(inter.domain1), domains) + inter.domaininds[2] = findfirst(isequal(inter.domain2), domains) + inter.parameterindexes[1] = length(p) + 1 + inter.parameterindexes[2] = length(p) + length(ps[i+length(domains)]) + p = vcat(p, ps[i+length(domains)]) end end @@ -390,7 +396,7 @@ end u: the current ODE state p: the ODE parameters t: the current ODE time - newW: a Bool which specifies whether the W matrix has been updated since the last call to precs. + newW: a Bool which specifies whether the W matrix has been updated since the last call to precs. It is recommended that this is checked to only update the preconditioner when newW == true. Plprev: the previous Pl. Prprev: the previous Pr. @@ -716,6 +722,8 @@ end evaluate(inter, dydt, domains, vT[inter.domaininds[1]], vT[inter.domaininds[2]], vphi[inter.domaininds[1]], vphi[inter.domaininds[2]], vGs[inter.domaininds[1]], vGs[inter.domaininds[2]], cstot, p) elseif isa(inter, VaporLiquidMassTransferInternalInterfaceConstantT) evaluate(inter, dydt, vV[inter.domaininds[1]], vV[inter.domaininds[2]], vT[inter.domaininds[1]], vT[inter.domaininds[2]], vN[inter.domaininds[1]], vN[inter.domaininds[2]], vP[inter.domaininds[1]], vP[inter.domaininds[2]], vCvave[inter.domaininds[1]], vCvave[inter.domaininds[2]], vns[inter.domaininds[1]], vns[inter.domaininds[2]], vUs[inter.domaininds[1]], vUs[inter.domaininds[2]], cstot, p) + elseif isa(inter, FickDiffusionInternalInterface) + evaluate(inter, dydt, cstot, p) end end for (i, domain) in enumerate(domains) @@ -1482,11 +1490,11 @@ end # Nrxns = length(rxns) # RTinv = 1.0/(R*T) # ratederiv .= 0.0 -# +# # for (j,rxn) in enumerate(rxns) # Nreact = length(rxn.reactantinds) # Nprod = length(rxn.productinds) -# +# # if Nreact == 1 # rind1 = rxn.reactantinds[1] # fderiv = cs[rind1] @@ -1497,7 +1505,7 @@ end # rind1,rind2,rind3 = rxn.reactantinds # fderiv = cs[rind1]*cs[rind2]*cs[rind3] # end -# +# # if Nprod == 1 # pind1 = rxn.productinds[1] # rderiv = krevs[j]/kfs[j]*cs[pind1] @@ -1508,12 +1516,12 @@ end # pind1,pind2,pind3 = rxn.productinds # rderiv = krevs[j]/kfs[j]*cs[pind1]*cs[pind2]*cs[pind3] # end -# +# # flux = fderiv-rderiv # gderiv = rderiv*kfs[j]*RTinv -# +# # deriv = zeros(Nspcs) -# +# # deriv[rind1] += gderiv # if Nreact > 1 # deriv[rind2] += gderiv @@ -1521,7 +1529,7 @@ end # deriv[rind3] == gderiv # end # end -# +# # deriv[pind1] -= gderiv # if Nprod > 1 # deriv[pind2] -= gderiv @@ -1529,7 +1537,7 @@ end # deriv[pind3] -= gderiv # end # end -# +# # ratederiv[rind1,j] -= flux # ratederiv[rind1,Nrxns+1:Nrxns+Nspcs] .-= deriv # if Nreact > 1 @@ -1540,7 +1548,7 @@ end # ratederiv[rind3,Nrxns+1:Nrxns+Nspcs] .-= deriv # end # end -# +# # ratederiv[pind1,j] += flux # ratederiv[pind1,Nrxns+1:Nrxns+Nspcs] .+= deriv # if Nprod > 1 @@ -1554,18 +1562,18 @@ end # end # return V*ratederiv # end -# +# # function jacobianp!(d::W; cs::Q,V::Y,T::Y2,Us::Z3,Cvave::Y3,N::Y2,kfs::Z,krevs::X,wV::Q2,ratederiv::Q3) where {W<:Union{ConstantVDomain,ParametrizedVDomain},Q3,Z3<:AbstractArray,Q<:AbstractArray,Q2<:AbstractArray,Y3<:Real,Y2<:Real,Y<:Real,Z<:AbstractArray,X<:AbstractArray} # Nspcs = length(cs) # rxns = d.phase.reactions # Nrxns = length(rxns) # RTinv = 1.0/(R*T) # ratederiv .= 0.0 -# +# # for (j,rxn) in enumerate(rxns) # Nreact = length(rxn.reactantinds) # Nprod = length(rxn.productinds) -# +# # if Nreact == 1 # rind1 = rxn.reactantinds[1] # fderiv = cs[rind1] @@ -1576,7 +1584,7 @@ end # rind1,rind2,rind3 = rxn.reactantinds # fderiv = cs[rind1]*cs[rind2]*cs[rind3] # end -# +# # if Nprod == 1 # pind1 = rxn.productinds[1] # rderiv = krevs[j]/kfs[j]*cs[pind1] @@ -1587,12 +1595,12 @@ end # pind1,pind2,pind3 = rxn.productinds # rderiv = krevs[j]/kfs[j]*cs[pind1]*cs[pind2]*cs[pind3] # end -# +# # flux = fderiv-rderiv # gderiv = rderiv*kfs[j]*RTinv -# +# # deriv = zeros(Nspcs) -# +# # deriv[rind1] += gderiv # if Nreact > 1 # deriv[rind2] += gderiv @@ -1600,7 +1608,7 @@ end # deriv[rind3] == gderiv # end # end -# +# # deriv[pind1] -= gderiv # if Nprod > 1 # deriv[pind2] -= gderiv @@ -1608,7 +1616,7 @@ end # deriv[pind3] -= gderiv # end # end -# +# # ratederiv[rind1,j] -= flux # ratederiv[rind1,Nrxns+1:Nrxns+Nspcs] .-= deriv # if Nreact > 1 @@ -1619,7 +1627,7 @@ end # ratederiv[rind3,Nrxns+1:Nrxns+Nspcs] .-= deriv # end # end -# +# # ratederiv[pind1,j] += flux # ratederiv[pind1,Nrxns+1:Nrxns+Nspcs] .+= deriv # if Nprod > 1 @@ -1651,7 +1659,7 @@ end # jac[inds[3],ind] += deriv # end # end -# +# # @inline function spreadpartials!(jac::S,deriv::T,inds::V,ind::Q,N::Q) where {S<:AbstractArray, T<:Real, V<:AbstractArray, Q<:Integer} # if N == 1 # jac[inds[1],ind] += deriv @@ -1664,7 +1672,7 @@ end # jac[inds[3],ind] += deriv # end # end -# +# # function jacobiany!(y::Array{T,1},t::T,domain::ConstantTPDomain,kfs::Array{T,1},krevs::Array{T,1},jac::P;zero::Bool=true) where {P<:AbstractArray,T<:Real,J<:Integer} # if zero # jac .= 0 @@ -1681,7 +1689,7 @@ end # krev = krevs[i] # if rxnarray[2,i] == 0 # jac[rxnarray[1,i],rxnarray[1,i]] -= kf -# if rxnarray[5,i] == 0 +# if rxnarray[5,i] == 0 # jac[rxnarray[4,i],rxnarray[1,i]] += kf # elseif rxnarray[6,i] == 0 # jac[rxnarray[4,i],rxnarray[1,i]] += kf @@ -1697,101 +1705,101 @@ end # deriv = 2*kf*cs[rxnarray[1,i]] # jac[rxnarray[1,i],rxnarray[1,i]] -= 2.0*deriv # for j in 1:Nspcs -# jac[rxnarray[1,i],j] -= 2.0*corr +# jac[rxnarray[1,i],j] -= 2.0*corr # end # jac[rxnarray[4,i],rxnarray[1,i]] += deriv -# for j in 1:Nspcs -# jac[rxnarray[4,i],j] += corr -# end +# for j in 1:Nspcs +# jac[rxnarray[4,i],j] += corr +# end # if rxnarray[5,i] != 0 -# jac[rxnarray[5,i],rxnarray[1,i]] += deriv -# for j = 1:Nspcs -# jac[rxnarray[5,i],j] += corr -# end -# if rxnarray[6,i] != 0 -# jac[rxnarray[6,i],rxnarray[1,i]] += deriv -# for j = 1:Nspcs -# jac[rxnarray[6,i],j] += corr -# end -# end -# end -# else +# jac[rxnarray[5,i],rxnarray[1,i]] += deriv +# for j = 1:Nspcs +# jac[rxnarray[5,i],j] += corr +# end +# if rxnarray[6,i] != 0 +# jac[rxnarray[6,i],rxnarray[1,i]] += deriv +# for j = 1:Nspcs +# jac[rxnarray[6,i],j] += corr +# end +# end +# end +# else # #derivative with respect to reactant 1 # deriv = kf*cs[rxnarray[2,i]] # jac[rxnarray[1,i],rxnarray[1,i]] -= deriv # jac[rxnarray[2,i],rxnarray[1,i]] -= deriv -# -# jac[rxnarray[4,i],rxnarray[1,i]] += deriv -# if rxnarray[5,i] != 0 +# +# jac[rxnarray[4,i],rxnarray[1,i]] += deriv +# if rxnarray[5,i] != 0 # jac[rxnarray[5,i],rxnarray[1,i]] += deriv -# if rxnarray[6,i] != 0 -# jac[rxnarray[6,i],rxnarray[1,i]] += deriv -# end -# end -# +# if rxnarray[6,i] != 0 +# jac[rxnarray[6,i],rxnarray[1,i]] += deriv +# end +# end +# # #derivative with respect to reactant 2 # deriv = kf*cs[rxnarray[1,i]] -# jac[rxnarray[1,i],rxnarray[2,i]] -= deriv -# jac[rxnarray[2,i],rxnarray[2,i]] -= deriv -# for j = 1:Nspcs -# jac[rxnarray[1,i],j] -= corr -# jac[rxnarray[2,i],j] -= corr +# jac[rxnarray[1,i],rxnarray[2,i]] -= deriv +# jac[rxnarray[2,i],rxnarray[2,i]] -= deriv +# for j = 1:Nspcs +# jac[rxnarray[1,i],j] -= corr +# jac[rxnarray[2,i],j] -= corr +# end +# jac[rxnarray[4,i],rxnarray[2,j]] += deriv +# if rxnarray[5,i] != 0 +# jac[rxnarray[5,i],rxnarray[2,i]] += deriv +# for j = 1:Nspcs +# jac[rxnarray[5,i],j] += corr +# end +# if rxnarray[6,i] != 0 +# jac[rxnarray[6,i],rxnarray[2,i]] += deriv +# for j = 1:Nspcs +# jac[rxnarray[6,i],j] += corr +# end +# end # end -# jac[rxnarray[4,i],rxnarray[2,j]] += deriv -# if rxnarray[5,i] != 0 -# jac[rxnarray[5,i],rxnarray[2,i]] += deriv -# for j = 1:Nspcs -# jac[rxnarray[5,i],j] += corr -# end -# if rxnarray[6,i] != 0 -# jac[rxnarray[6,i],rxnarray[2,i]] += deriv -# for j = 1:Nspcs -# jac[rxnarray[6,i],j] += corr -# end -# end -# end # end # else # corr = -2.0*kf*cs[rxnarray[1,i]]*cs[rxnarray[2,i]]*cs[rxnarray[3,i]]/C # if (rxnarray[1,i] == rxnarray[2,i] && rxnarray[1,i] == rxnarray[3,i]) # deriv = 3.0*kf*cs[rxnarray[1,i]]*cs[rxnarray[1,i]] -# jac[rxnarray[1,i],rxnarray[1,i]] -= 3.0*deriv -# for j = 1:Nspcs -# jac[rxnarray[1,i],j] -= 3.0*corr -# end -# jac[rxnarray[4,i],rxnarray[1,i]] += deriv -# for j = 1:Nspcs -# jac[rxnarray[4,i],j] += corr -# end -# if rxnarray[5,i] != 0 -# jac[rxnarray[5,i],rxnarray[1,i]] += deriv -# for j = 1:Nspcs -# jac[rxnarray[5,i],j] += corr -# end -# if rxnarray[6,i] != 0 -# jac[rxnarray[6,i],rxnarray[1,i]] += deriv -# for j = 1:Nspcs -# jac[rxnarray[6,i],j] += corr -# end -# end -# end +# jac[rxnarray[1,i],rxnarray[1,i]] -= 3.0*deriv +# for j = 1:Nspcs +# jac[rxnarray[1,i],j] -= 3.0*corr +# end +# jac[rxnarray[4,i],rxnarray[1,i]] += deriv +# for j = 1:Nspcs +# jac[rxnarray[4,i],j] += corr +# end +# if rxnarray[5,i] != 0 +# jac[rxnarray[5,i],rxnarray[1,i]] += deriv +# for j = 1:Nspcs +# jac[rxnarray[5,i],j] += corr +# end +# if rxnarray[6,i] != 0 +# jac[rxnarray[6,i],rxnarray[1,i]] += deriv +# for j = 1:Nspcs +# jac[rxnarray[6,i],j] += corr +# end +# end +# end # elseif rxnarray[1,i] == rxnarray[2,i] -# #derivative with respect to reactant 1 +# #derivative with respect to reactant 1 # deriv = 2.0*kf*cs[rxnarray[1,i]]*cs[rxnarray[3,i]] -# jac[rxnarray[1,i],rxnarray[1,i]] -= 2.0*deriv +# jac[rxnarray[1,i],rxnarray[1,i]] -= 2.0*deriv # jac[rxnarray[3,i],rxnarray[1,i]] -= deriv -# -# jac[rxnarray[4,i],rxnarray[1,i]] += deriv -# if rxnarray[5,i] != 0 -# jac[rxnarray[5,i],rxnarray[1,i]] += deriv -# if rxnarray[6,i] != 0 -# jac[rxnarray[6,i],rxnarray[1,i]] += deriv -# end -# end -# +# +# jac[rxnarray[4,i],rxnarray[1,i]] += deriv +# if rxnarray[5,i] != 0 +# jac[rxnarray[5,i],rxnarray[1,i]] += deriv +# if rxnarray[6,i] != 0 +# jac[rxnarray[6,i],rxnarray[1,i]] += deriv +# end +# end +# # #derivative with respect to reactant 3 # deriv = kf*cs[rxnarray] -# +# # ind1,ind2,ind3 = rxn.reactantinds # corr = -2.0*kf*cs[ind1]*cs[ind2]*cs[ind3]/C # deriv = kf*cs[ind1]*cs[ind2] @@ -1877,7 +1885,7 @@ end # end # return jac # end -# +# # function jacobiany!(y::Array{T,1},t::T,domain::ConstantTPDomain,kfs::Array{T,1},krevs::Array{T,1},jac::P;zero::Bool=true) where {P<:AbstractArray,T<:Real,J<:Integer} # if zero # jac .= 0 diff --git a/src/TestReactors.jl b/src/TestReactors.jl index 0dfbe1c0..bc4e616c 100644 --- a/src/TestReactors.jl +++ b/src/TestReactors.jl @@ -141,10 +141,71 @@ using SciMLSensitivity name = "oxygen" ind = findfirst(x -> x == name, liqspcnames) - @test sol(sol.t[end])[ind] ≈ 0.11758959354431776 rtol = 1e-4 #test there are oxygen dissolved into the liquid + @test sol(sol.t[end])[ind] ≈ 0.11758959354431776 rtol = 1e-4 #test there are oxygen dissolved into the liquid end + @testset "Test cross-domain diffusion problem with FickDiffusionInternalInterface" begin + # Unit test of the Fickian diffusion interface: an inert, dilute O2 tracer + # diffusing between two liquid n-octane volumes at 298 K. + # Two well-mixed volumes connected by conductance g = A*D/delta obey the + # following governing equation: + # d(c1-c2)/dt = -k*(c1-c2), where k = g*(1/V1 + 1/V2) + # => (c1-c2)(t) = (c1-c2)(0)*exp(-k*t), equilibrating to c_eq = n_tot/(V1+V2). + phaseDict = readinput("../src/testing/liquid_phase.rms") + spcs = phaseDict["phase"]["Species"] + solv = phaseDict["Solvents"][1] + liq = IdealDiluteSolution(spcs, [], solv; name="phase", diffusionlimited=true) + + spcnames = getfield.(liq.species, :name) + oxygenind = findfirst(isequal("oxygen"), spcnames) + + T = 298.15 # K + V1 = 1.0e-3 # m^3 (1 L) + V2 = 3.0e-3 # m^3 (3 L) + A = 1.0e-2 # m^2, interfacial area + delta = 1.0e-4 # m, layer thickness (0.1 mm) + + c_octane = 6.154e3 # mol/m^3, liquid n-octane molar density + c_O2 = 5.0 # mol/m^3, dilute dissolved O2 + n_O2 = c_O2 * V1 # total O2 moles + + initialconds1 = Dict(["T" => T, "V" => V1, "octane" => c_octane * V1, "oxygen" => n_O2]) + domain1, y01, p1 = ConstantTVDomain(phase=liq, initialconds=initialconds1) + initialconds2 = Dict(["T" => T, "V" => V2, "octane" => c_octane * V2, "oxygen" => 0.0]) + domain2, y02, p2 = ConstantTVDomain(phase=liq, initialconds=initialconds2) + + fick, pfick = FickDiffusionInternalInterface(domain1, domain2, A, delta) + + D_O2 = domain1.diffusivity[oxygenind] + g = A * D_O2 / delta # effective mass transfer coefficient + k = g * (1.0 / V1 + 1.0 / V2) # time constant for exponential relaxation of the concentration difference across the boundary + c_eq = n_O2 / (V1 + V2) # equilibrium concentration of O2 in both domains + delta0 = n_O2 / V1 # initial concentration difference across the boundary + tf = 10.0 / k # arbitrary long time to reach equilibrium (10 times the time constant 1/k) + + domains = (domain1, domain2) + interfaces = [fick] + react, y0, p = Reactor(domains, (y01, y02), (0.0, tf), interfaces, (p1, p2, pfick)) + sol = solve(react.ode, react.recommendedsolver, abstol=1e-16, reltol=1e-10) + + i1 = domain1.indexes[1] - 1 + oxygenind + i2 = domain2.indexes[1] - 1 + oxygenind + + # mass conservation + y_end = sol.u[end] + @test y_end[i1] + y_end[i2] ≈ n_O2 rtol = 1e-6 + + # concentration difference at t = 1/k + y_tau = sol(1.0 / k) + delta_tau = y_tau[i1] / V1 - y_tau[i2] / V2 + @test delta_tau ≈ delta0 * exp(-1.0) rtol = 1e-3 + + # test whether equilibrium is reached at tf + @test y_end[i1] / V1 ≈ c_eq rtol = 1e-2 + @test y_end[i2] / V2 ≈ c_eq rtol = 1e-2 + end + @testset "Test liquid phase Parametrized T Constant V reactor jacobian" begin #Parametrized T constant V Ideal Dilute Liquid initialconds = Dict(["ts" => [0.0, 600.0, 1200.0], "T" => [450.0, 490.0, 500.0], "V" => 1.0e-6 * 1e6, "octane" => 6.154e-3 * 1e6, "oxygen" => 4.953e-6 * 1e6]) @@ -212,10 +273,10 @@ using SciMLSensitivity @test y[o2ind] / N ≈ 0.200419093 rtol = 1e-4 @test y[h2oind] / N ≈ 0.386618602 rtol = 1e-4 - #plotting + #plotting plotrops(sim, "H2", 20.4402454, N=10) getfluxdiagram(sim, 20.4402454) - + #analytic jacobian vs. ForwardDiff jacobian t = 20.44002454 y = sol(t) @@ -645,4 +706,4 @@ using SciMLSensitivity @test length(dmech.spcs) == 4 @test dmech.rxns[1].index == 11 end -end; \ No newline at end of file +end;