add jacobianp! for ConstantPDomain and ParametrizedPDomain

Author: hwpang

src/Domain.jl

@@ -2232,6 +2232,46 @@ function jacobianp!(jacp::Q,y::U,p::W,t::Z,domain::D,interfaces::Q3,colorvec::Q2
     end
 end
 
+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}}
+    jacp.=0.0
+    ns,cs,T,P,V,C,N,mu,kfs,krevs,Hs,Us,Gs,diffs,Cvave,cpdivR = calcthermo(domain,y,t,p)
+    
+    Nspcs = length(cs)
+    Nrxns = size(domain.rxnarray)[2]
+    
+    dydt = zeros(size(y))
+    addreactionratecontributions!(dydt,domain.rxnarray,cs,kfs,krevs)
+    dydt .*= V
+    
+    jacobianpnsderiv!(jacp,y,p,t,domain,domain.rxnarray,cs,T,V,kfs,krevs,Nspcs,Nrxns)
+    
+    @fastmath Cpave = Cvave+R
+    @simd for i in 1:Nspcs
+        @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 ω
+        @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
+    end
+
+    @simd for i in Nspcs+1:Nspcs+Nrxns
+        @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 ω
+        @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
+    end
+    
+    @simd for ind in domain.constantspeciesinds
+        @inbounds jacp[ind,:] .= 0.
+    end
+    
+    @simd for inter in interfaces
+        if isa(inter,Inlet) && domain == inter.domain
+            flow = inter.F(t)
+            for i in 1:Nspcs
+                ddGidTdt = flow*(- ns[i]/N)/(N*Cpave)
+                jacp[domain.indexes[3],i] += ddGidTdt
+                jacp[domain.indexes[4],i] += ddGidTdt*V/T 
+            end
+        end
+    end
+end
+
 function getreactionindices(ig::Q) where {Q<:AbstractPhase}
     arr = zeros(Int64,(6,length(ig.reactions)))
     for (i,rxn) in enumerate(ig.reactions)