diff --git a/src/Domain.jl b/src/Domain.jl index 0fd269c3..b1868bc7 100644 --- a/src/Domain.jl +++ b/src/Domain.jl @@ -2652,32 +2652,40 @@ export jacobiany! @fastmath @inbounds fderiv = kfs[rxnind]*cs[rxnarray[1,rxnind]] elseif @inbounds rxnarray[3,rxnind] == 0 @fastmath @inbounds fderiv = kfs[rxnind]*cs[rxnarray[1,rxnind]]*cs[rxnarray[2,rxnind]] - else + elseif @inbounds rxnarray[4,rxnind] == 0 @fastmath @inbounds fderiv = kfs[rxnind]*cs[rxnarray[1,rxnind]]*cs[rxnarray[2,rxnind]]*cs[rxnarray[3,rxnind]] + else + @fastmath @inbounds fderiv = kfs[rxnind]*cs[rxnarray[1,rxnind]]*cs[rxnarray[2,rxnind]]*cs[rxnarray[3,rxnind]]*cs[rxnarray[4,rxnind]] end - if @inbounds rxnarray[5,rxnind] == 0 - @fastmath @inbounds rderiv = krevs[rxnind]*cs[rxnarray[4,rxnind]] - elseif @inbounds rxnarray[6,rxnind] == 0 - @fastmath @inbounds rderiv = krevs[rxnind]*cs[rxnarray[4,rxnind]]*cs[rxnarray[5,rxnind]] + if @inbounds rxnarray[6,rxnind] == 0 + @fastmath @inbounds rderiv = krevs[rxnind]*cs[rxnarray[5,rxnind]] + elseif @inbounds rxnarray[7,rxnind] == 0 + @fastmath @inbounds rderiv = krevs[rxnind]*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]] + elseif @inbounds rxnarray[8,rxnind] == 0 + @fastmath @inbounds rderiv = krevs[rxnind]*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]]*cs[rxnarray[7,rxnind]] else - @fastmath @inbounds rderiv = krevs[rxnind]*cs[rxnarray[4,rxnind]]*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]] + @fastmath @inbounds rderiv = krevs[rxnind]*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]]*cs[rxnarray[7,rxnind]]*cs[rxnarray[8,rxnind]] end else if @inbounds rxnarray[2,rxnind] == 0 @inbounds fderiv = cs[rxnarray[1,rxnind]] elseif @inbounds rxnarray[3,rxnind] == 0 @fastmath @inbounds fderiv = cs[rxnarray[1,rxnind]]*cs[rxnarray[2,rxnind]] - else + elseif @inbounds rxnarray[4,rxnind] == 0 @fastmath @inbounds fderiv = cs[rxnarray[1,rxnind]]*cs[rxnarray[2,rxnind]]*cs[rxnarray[3,rxnind]] + else + @fastmath @inbounds fderiv = cs[rxnarray[1,rxnind]]*cs[rxnarray[2,rxnind]]*cs[rxnarray[3,rxnind]]*cs[rxnarray[4,rxnind]] end - if @inbounds rxnarray[5,rxnind] == 0 - @fastmath @inbounds rderiv = krevs[rxnind]/kfs[rxnind]*cs[rxnarray[4,rxnind]] - elseif @inbounds rxnarray[6,rxnind] == 0 - @fastmath @inbounds rderiv = krevs[rxnind]/kfs[rxnind]*cs[rxnarray[4,rxnind]]*cs[rxnarray[5,rxnind]] + if @inbounds rxnarray[6,rxnind] == 0 + @fastmath @inbounds rderiv = krevs[rxnind]/kfs[rxnind]*cs[rxnarray[5,rxnind]] + elseif @inbounds rxnarray[7,rxnind] == 0 + @fastmath @inbounds rderiv = krevs[rxnind]/kfs[rxnind]*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]] + elseif @inbounds rxnarray[8,rxnind] == 0 + @fastmath @inbounds rderiv = krevs[rxnind]/kfs[rxnind]*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]]*cs[rxnarray[7,rxnind]] else - @fastmath @inbounds rderiv = krevs[rxnind]/kfs[rxnind]*cs[rxnarray[4,rxnind]]*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]] + @fastmath @inbounds rderiv = krevs[rxnind]/kfs[rxnind]*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]]*cs[rxnarray[7,rxnind]]*cs[rxnarray[8,rxnind]] end end @@ -2701,23 +2709,43 @@ export jacobiany! @inbounds jacp[rxnarray[2,rxnind],rxnarray[3,rxnind]] -= gderiv @inbounds jacp[rxnarray[3,rxnind],rxnarray[3,rxnind]] -= gderiv @inbounds _spreadreactantpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[3,rxnind]) + if @inbounds rxnarray[4,rxnind] !== 0 + @inbounds jacp[rxnarray[4,rxnind],rxnarray[1,rxnind]] -= gderiv + @inbounds jacp[rxnarray[4,rxnind],rxnarray[2,rxnind]] -= gderiv + @inbounds jacp[rxnarray[4,rxnind],rxnarray[3,rxnind]] -= gderiv + @inbounds jacp[rxnarray[1,rxnind],rxnarray[4,rxnind]] -= gderiv + @inbounds jacp[rxnarray[2,rxnind],rxnarray[4,rxnind]] -= gderiv + @inbounds jacp[rxnarray[3,rxnind],rxnarray[4,rxnind]] -= gderiv + @inbounds jacp[rxnarray[4,rxnind],rxnarray[4,rxnind]] -= gderiv + @inbounds _spreadreactantpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[4,rxnind]) + end end end - @inbounds jacp[rxnarray[4,rxnind],rxnarray[4,rxnind]] -= gderiv - @inbounds _spreadproductpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[4,rxnind]) - if @inbounds rxnarray[5,rxnind] !== 0 - @inbounds jacp[rxnarray[5,rxnind],rxnarray[4,rxnind]] -= gderiv - @inbounds jacp[rxnarray[4,rxnind],rxnarray[5,rxnind]] -= gderiv - @inbounds jacp[rxnarray[5,rxnind],rxnarray[5,rxnind]] -= gderiv - @inbounds _spreadproductpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[5,rxnind]) - if @inbounds rxnarray[6,rxnind] !== 0 - @inbounds jacp[rxnarray[6,rxnind],rxnarray[4,rxnind]] -= gderiv - @inbounds jacp[rxnarray[6,rxnind],rxnarray[5,rxnind]] -= gderiv - @inbounds jacp[rxnarray[4,rxnind],rxnarray[6,rxnind]] -= gderiv - @inbounds jacp[rxnarray[5,rxnind],rxnarray[6,rxnind]] -= gderiv - @inbounds jacp[rxnarray[6,rxnind],rxnarray[6,rxnind]] -= gderiv - @inbounds _spreadproductpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[6,rxnind]) + @inbounds jacp[rxnarray[5,rxnind],rxnarray[5,rxnind]] -= gderiv + @inbounds _spreadproductpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[5,rxnind]) + if @inbounds rxnarray[6,rxnind] !== 0 + @inbounds jacp[rxnarray[6,rxnind],rxnarray[5,rxnind]] -= gderiv + @inbounds jacp[rxnarray[5,rxnind],rxnarray[6,rxnind]] -= gderiv + @inbounds jacp[rxnarray[6,rxnind],rxnarray[6,rxnind]] -= gderiv + @inbounds _spreadproductpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[6,rxnind]) + if @inbounds rxnarray[7,rxnind] !== 0 + @inbounds jacp[rxnarray[7,rxnind],rxnarray[5,rxnind]] -= gderiv + @inbounds jacp[rxnarray[7,rxnind],rxnarray[6,rxnind]] -= gderiv + @inbounds jacp[rxnarray[5,rxnind],rxnarray[7,rxnind]] -= gderiv + @inbounds jacp[rxnarray[6,rxnind],rxnarray[7,rxnind]] -= gderiv + @inbounds jacp[rxnarray[7,rxnind],rxnarray[7,rxnind]] -= gderiv + @inbounds _spreadproductpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[7,rxnind]) + if @inbounds rxnarray[8,rxnind] !== 0 + @inbounds jacp[rxnarray[8,rxnind],rxnarray[5,rxnind]] -= gderiv + @inbounds jacp[rxnarray[8,rxnind],rxnarray[6,rxnind]] -= gderiv + @inbounds jacp[rxnarray[8,rxnind],rxnarray[7,rxnind]] -= gderiv + @inbounds jacp[rxnarray[5,rxnind],rxnarray[8,rxnind]] -= gderiv + @inbounds jacp[rxnarray[6,rxnind],rxnarray[8,rxnind]] -= gderiv + @inbounds jacp[rxnarray[7,rxnind],rxnarray[8,rxnind]] -= gderiv + @inbounds jacp[rxnarray[8,rxnind],rxnarray[8,rxnind]] -= gderiv + @inbounds _spreadproductpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[8,rxnind]) + end end end @@ -2732,16 +2760,20 @@ end @fastmath @inbounds fderiv = kfs[rxnind]*cs[rxnarray[1,rxnind]] elseif @inbounds rxnarray[3,rxnind] == 0 @fastmath @inbounds fderiv = kfs[rxnind]*cs[rxnarray[1,rxnind]]*cs[rxnarray[2,rxnind]] - else + elseif @inbounds rxnarray[4,rxnind] == 0 @fastmath @inbounds fderiv = kfs[rxnind]*cs[rxnarray[1,rxnind]]*cs[rxnarray[2,rxnind]]*cs[rxnarray[3,rxnind]] + else + @fastmath @inbounds fderiv = kfs[rxnind]*cs[rxnarray[1,rxnind]]*cs[rxnarray[2,rxnind]]*cs[rxnarray[3,rxnind]]*cs[rxnarray[4,rxnind]] end - if @inbounds rxnarray[5,rxnind] == 0 - @fastmath @inbounds rderiv = krevs[rxnind]*cs[rxnarray[4,rxnind]] - elseif @inbounds rxnarray[6,rxnind] == 0 - @fastmath @inbounds rderiv = krevs[rxnind]*cs[rxnarray[4,rxnind]]*cs[rxnarray[5,rxnind]] + if @inbounds rxnarray[6,rxnind] == 0 + @fastmath @inbounds rderiv = krevs[rxnind]*cs[rxnarray[5,rxnind]] + elseif @inbounds rxnarray[7,rxnind] == 0 + @fastmath @inbounds rderiv = krevs[rxnind]*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]] + elseif @inbounds rxnarray[8,rxnind] == 0 + @fastmath @inbounds rderiv = krevs[rxnind]*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]]*cs[rxnarray[7,rxnind]] else - @fastmath @inbounds rderiv = krevs[rxnind]*cs[rxnarray[4,rxnind]]*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]] + @fastmath @inbounds rderiv = krevs[rxnind]*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]]*cs[rxnarray[7,rxnind]]*cs[rxnarray[8,rxnind]] end @fastmath flux = fderiv-rderiv @@ -2764,23 +2796,43 @@ end @inbounds jacp[rxnarray[2,rxnind],rxnarray[3,rxnind]] -= gderiv @inbounds jacp[rxnarray[3,rxnind],rxnarray[3,rxnind]] -= gderiv @inbounds _spreadreactantpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[3,rxnind]) + if @inbounds rxnarray[4,rxnind] !== 0 + @inbounds jacp[rxnarray[4,rxnind],rxnarray[1,rxnind]] -= gderiv + @inbounds jacp[rxnarray[4,rxnind],rxnarray[2,rxnind]] -= gderiv + @inbounds jacp[rxnarray[4,rxnind],rxnarray[3,rxnind]] -= gderiv + @inbounds jacp[rxnarray[1,rxnind],rxnarray[4,rxnind]] -= gderiv + @inbounds jacp[rxnarray[2,rxnind],rxnarray[4,rxnind]] -= gderiv + @inbounds jacp[rxnarray[3,rxnind],rxnarray[4,rxnind]] -= gderiv + @inbounds jacp[rxnarray[4,rxnind],rxnarray[4,rxnind]] -= gderiv + @inbounds _spreadreactantpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[4,rxnind]) + end end end - @inbounds jacp[rxnarray[4,rxnind],rxnarray[4,rxnind]] -= gderiv - @inbounds _spreadproductpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[4,rxnind]) - if @inbounds rxnarray[5,rxnind] !== 0 - @inbounds jacp[rxnarray[5,rxnind],rxnarray[4,rxnind]] -= gderiv - @inbounds jacp[rxnarray[4,rxnind],rxnarray[5,rxnind]] -= gderiv - @inbounds jacp[rxnarray[5,rxnind],rxnarray[5,rxnind]] -= gderiv - @inbounds _spreadproductpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[5,rxnind]) - if @inbounds rxnarray[6,rxnind] !== 0 - @inbounds jacp[rxnarray[6,rxnind],rxnarray[4,rxnind]] -= gderiv - @inbounds jacp[rxnarray[6,rxnind],rxnarray[5,rxnind]] -= gderiv - @inbounds jacp[rxnarray[4,rxnind],rxnarray[6,rxnind]] -= gderiv - @inbounds jacp[rxnarray[5,rxnind],rxnarray[6,rxnind]] -= gderiv - @inbounds jacp[rxnarray[6,rxnind],rxnarray[6,rxnind]] -= gderiv - @inbounds _spreadproductpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[6,rxnind]) + @inbounds jacp[rxnarray[5,rxnind],rxnarray[5,rxnind]] -= gderiv + @inbounds _spreadproductpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[5,rxnind]) + if @inbounds rxnarray[6,rxnind] !== 0 + @inbounds jacp[rxnarray[6,rxnind],rxnarray[5,rxnind]] -= gderiv + @inbounds jacp[rxnarray[5,rxnind],rxnarray[6,rxnind]] -= gderiv + @inbounds jacp[rxnarray[6,rxnind],rxnarray[6,rxnind]] -= gderiv + @inbounds _spreadproductpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[6,rxnind]) + if @inbounds rxnarray[7,rxnind] !== 0 + @inbounds jacp[rxnarray[7,rxnind],rxnarray[5,rxnind]] -= gderiv + @inbounds jacp[rxnarray[7,rxnind],rxnarray[6,rxnind]] -= gderiv + @inbounds jacp[rxnarray[5,rxnind],rxnarray[7,rxnind]] -= gderiv + @inbounds jacp[rxnarray[6,rxnind],rxnarray[7,rxnind]] -= gderiv + @inbounds jacp[rxnarray[7,rxnind],rxnarray[7,rxnind]] -= gderiv + @inbounds _spreadproductpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[7,rxnind]) + if @inbounds rxnarray[8,rxnind] !== 0 + @inbounds jacp[rxnarray[8,rxnind],rxnarray[5,rxnind]] -= gderiv + @inbounds jacp[rxnarray[8,rxnind],rxnarray[6,rxnind]] -= gderiv + @inbounds jacp[rxnarray[8,rxnind],rxnarray[7,rxnind]] -= gderiv + @inbounds jacp[rxnarray[5,rxnind],rxnarray[8,rxnind]] -= gderiv + @inbounds jacp[rxnarray[6,rxnind],rxnarray[8,rxnind]] -= gderiv + @inbounds jacp[rxnarray[7,rxnind],rxnarray[8,rxnind]] -= gderiv + @inbounds jacp[rxnarray[8,rxnind],rxnarray[8,rxnind]] -= gderiv + @inbounds _spreadproductpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[8,rxnind]) + end end end end @@ -2794,16 +2846,20 @@ end @fastmath @inbounds fderiv = kfs[rxnind]*kfs[rxnind]/(p[Nspcs+rxnind]*p[Nspcs+rxnind])*cs[rxnarray[1,rxnind]] elseif @inbounds rxnarray[3,rxnind] == 0 @fastmath @inbounds fderiv = kfs[rxnind]*kfs[rxnind]/(p[Nspcs+rxnind]*p[Nspcs+rxnind])*cs[rxnarray[1,rxnind]]*cs[rxnarray[2,rxnind]] - else + elseif rxnarray[4,rxnind] == 0 @fastmath @inbounds fderiv = kfs[rxnind]*kfs[rxnind]/(p[Nspcs+rxnind]*p[Nspcs+rxnind])*cs[rxnarray[1,rxnind]]*cs[rxnarray[2,rxnind]]*cs[rxnarray[3,rxnind]] + else + @fastmath @inbounds fderiv = kfs[rxnind]*kfs[rxnind]/(p[Nspcs+rxnind]*p[Nspcs+rxnind])*cs[rxnarray[1,rxnind]]*cs[rxnarray[2,rxnind]]*cs[rxnarray[3,rxnind]]*cs[rxnarray[4,rxnind]] end - if @inbounds rxnarray[5,rxnind] == 0 - @fastmath @inbounds rderiv = kfs[rxnind]*krevs[rxnind]/(p[Nspcs+rxnind]*p[Nspcs+rxnind])*cs[rxnarray[4,rxnind]] - elseif @inbounds rxnarray[6,rxnind] == 0 - @fastmath @inbounds rderiv = kfs[rxnind]*krevs[rxnind]/(p[Nspcs+rxnind]*p[Nspcs+rxnind])*cs[rxnarray[4,rxnind]]*cs[rxnarray[5,rxnind]] + if @inbounds rxnarray[6,rxnind] == 0 + @fastmath @inbounds rderiv = kfs[rxnind]*krevs[rxnind]/(p[Nspcs+rxnind]*p[Nspcs+rxnind])*cs[rxnarray[5,rxnind]] + elseif @inbounds rxnarray[7,rxnind] == 0 + @fastmath @inbounds rderiv = kfs[rxnind]*krevs[rxnind]/(p[Nspcs+rxnind]*p[Nspcs+rxnind])*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]] + elseif @inbounds rxnarray[8,rxnind] == 0 + @fastmath @inbounds rderiv = kfs[rxnind]*krevs[rxnind]/(p[Nspcs+rxnind]*p[Nspcs+rxnind])*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]]*cs[rxnarray[7,rxnind]] else - @fastmath @inbounds rderiv = kfs[rxnind]*krevs[rxnind]/(p[Nspcs+rxnind]*p[Nspcs+rxnind])*cs[rxnarray[4,rxnind]]*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]] + @fastmath @inbounds rderiv = kfs[rxnind]*krevs[rxnind]/(p[Nspcs+rxnind]*p[Nspcs+rxnind])*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]]*cs[rxnarray[7,rxnind]]*cs[rxnarray[8,rxnind]] end @fastmath flux = fderiv-rderiv @@ -2826,23 +2882,43 @@ end @inbounds jacp[rxnarray[2,rxnind],rxnarray[3,rxnind]] -= gderiv @inbounds jacp[rxnarray[3,rxnind],rxnarray[3,rxnind]] -= gderiv @inbounds _spreadreactantpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[3,rxnind]) + if @inbounds rxnarray[4,rxnind] !== 0 + @inbounds jacp[rxnarray[4,rxnind],rxnarray[1,rxnind]] -= gderiv + @inbounds jacp[rxnarray[4,rxnind],rxnarray[2,rxnind]] -= gderiv + @inbounds jacp[rxnarray[4,rxnind],rxnarray[3,rxnind]] -= gderiv + @inbounds jacp[rxnarray[1,rxnind],rxnarray[4,rxnind]] -= gderiv + @inbounds jacp[rxnarray[2,rxnind],rxnarray[4,rxnind]] -= gderiv + @inbounds jacp[rxnarray[3,rxnind],rxnarray[4,rxnind]] -= gderiv + @inbounds jacp[rxnarray[4,rxnind],rxnarray[4,rxnind]] -= gderiv + @inbounds _spreadreactantpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[4,rxnind]) + end end end - @inbounds jacp[rxnarray[4,rxnind],rxnarray[4,rxnind]] -= gderiv - @inbounds _spreadproductpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[4,rxnind]) - if @inbounds rxnarray[5,rxnind] !== 0 - @inbounds jacp[rxnarray[5,rxnind],rxnarray[4,rxnind]] -= gderiv - @inbounds jacp[rxnarray[4,rxnind],rxnarray[5,rxnind]] -= gderiv - @inbounds jacp[rxnarray[5,rxnind],rxnarray[5,rxnind]] -= gderiv - @inbounds _spreadproductpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[5,rxnind]) - if @inbounds rxnarray[6,rxnind] !== 0 - @inbounds jacp[rxnarray[6,rxnind],rxnarray[4,rxnind]] -= gderiv - @inbounds jacp[rxnarray[6,rxnind],rxnarray[5,rxnind]] -= gderiv - @inbounds jacp[rxnarray[4,rxnind],rxnarray[6,rxnind]] -= gderiv - @inbounds jacp[rxnarray[5,rxnind],rxnarray[6,rxnind]] -= gderiv - @inbounds jacp[rxnarray[6,rxnind],rxnarray[6,rxnind]] -= gderiv - @inbounds _spreadproductpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[6,rxnind]) + @inbounds jacp[rxnarray[5,rxnind],rxnarray[5,rxnind]] -= gderiv + @inbounds _spreadproductpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[5,rxnind]) + if @inbounds rxnarray[6,rxnind] !== 0 + @inbounds jacp[rxnarray[6,rxnind],rxnarray[5,rxnind]] -= gderiv + @inbounds jacp[rxnarray[5,rxnind],rxnarray[6,rxnind]] -= gderiv + @inbounds jacp[rxnarray[6,rxnind],rxnarray[6,rxnind]] -= gderiv + @inbounds _spreadproductpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[6,rxnind]) + if @inbounds rxnarray[7,rxnind] !== 0 + @inbounds jacp[rxnarray[7,rxnind],rxnarray[5,rxnind]] -= gderiv + @inbounds jacp[rxnarray[7,rxnind],rxnarray[6,rxnind]] -= gderiv + @inbounds jacp[rxnarray[5,rxnind],rxnarray[7,rxnind]] -= gderiv + @inbounds jacp[rxnarray[6,rxnind],rxnarray[7,rxnind]] -= gderiv + @inbounds jacp[rxnarray[7,rxnind],rxnarray[7,rxnind]] -= gderiv + @inbounds _spreadproductpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[7,rxnind]) + if @inbounds rxnarray[8,rxnind] !== 0 + @inbounds jacp[rxnarray[8,rxnind],rxnarray[5,rxnind]] -= gderiv + @inbounds jacp[rxnarray[8,rxnind],rxnarray[6,rxnind]] -= gderiv + @inbounds jacp[rxnarray[8,rxnind],rxnarray[7,rxnind]] -= gderiv + @inbounds jacp[rxnarray[5,rxnind],rxnarray[8,rxnind]] -= gderiv + @inbounds jacp[rxnarray[6,rxnind],rxnarray[8,rxnind]] -= gderiv + @inbounds jacp[rxnarray[7,rxnind],rxnarray[8,rxnind]] -= gderiv + @inbounds jacp[rxnarray[8,rxnind],rxnarray[8,rxnind]] -= gderiv + @inbounds _spreadproductpartials!(jacp,gderiv,rxnarray,rxnind,rxnarray[8,rxnind]) + end end end end diff --git a/src/EdgeAnalysis.jl b/src/EdgeAnalysis.jl index 0e6da777..fc1d9f4f 100644 --- a/src/EdgeAnalysis.jl +++ b/src/EdgeAnalysis.jl @@ -208,8 +208,8 @@ function getkeyselectioninds(coreedgedomains,coreedgeinters,domains,inters) rxnindexedge += length(coreedgedomains[i].phase.reactions) indend = length(domains[i].phase.reactions) - reactantindices[:,ind:ind+indend-1] = domains[i].rxnarray[1:3,:] - productindices[:,ind:ind+indend-1] = domains[i].rxnarray[4:6,:] + reactantindices[:,ind:ind+indend-1] = domains[i].rxnarray[1:4,:] + productindices[:,ind:ind+indend-1] = domains[i].rxnarray[5:8,:] ind += indend end @@ -220,8 +220,8 @@ function getkeyselectioninds(coreedgedomains,coreedgeinters,domains,inters) index += length(coreedgeinters[i].phase.reactions) indend = length(inters[i].reactions) - reactantindices[:,ind:ind+indend] = inters[i].rxnarray[1:3,:] - productindices[:,ind:ind+indend] = inters[i].rxnarray[4:6,:] + reactantindices[:,ind:ind+indend] = inters[i].rxnarray[1:4,:] + productindices[:,ind:ind+indend] = inters[i].rxnarray[5:8,:] ind += indend end end @@ -251,8 +251,8 @@ function getkeyselectioninds(coreedgedomain::AbstractDomain,coreedgeinters,domai end end - @views @inbounds reactantindices = coreedgedomain.rxnarray[1:3,:] - @views @inbounds productindices = coreedgedomain.rxnarray[4:6,:] + @views @inbounds reactantindices = coreedgedomain.rxnarray[1:4,:] + @views @inbounds productindices = coreedgedomain.rxnarray[5:8,:] coretoedgespcmap = Dict([i=>findfirst(isequal(spc),coreedgedomain.phase.species) for (i,spc) in enumerate(domain.phase.species)]) @simd for j = 3:length(domain.indexes) coretoedgespcmap[domain.indexes[j]] = coreedgedomain.indexes[j] @@ -286,7 +286,7 @@ function processfluxes(sim::SystemSimulation, if any(d.rxnarray[:,i].>length(corespeciesconcentrations)) continue end - for j = 1:3 + for j = 1:4 if d.rxnarray[j,i] != 0 corespeciesconsumptionrates[d.rxnarray[j,i]] += frts[i+index] corespeciesproductionrates[d.rxnarray[j,i]] += rrts[i+index] @@ -294,7 +294,7 @@ function processfluxes(sim::SystemSimulation, break end end - for j = 4:6 + for j = 5:8 if d.rxnarray[j,i] != 0 corespeciesproductionrates[d.rxnarray[j,i]] += frts[i+index] corespeciesconsumptionrates[d.rxnarray[j,i]] += rrts[i+index] @@ -311,7 +311,7 @@ function processfluxes(sim::SystemSimulation, if any(d.rxnarray[:,i].>length(corespeciesconcentrations)) continue end - for j = 1:3 + for j = 1:4 if d.rxnarray[j,i] != 0 corespeciesconsumptionrates[d.rxnarray[j,i]] += frts[i+index] corespeciesproductionrates[d.rxnarray[j,i]] += rrts[i+index] @@ -319,7 +319,7 @@ function processfluxes(sim::SystemSimulation, break end end - for j = 4:6 + for j = 5:8 if d.rxnarray[j,i] != 0 corespeciesproductionrates[d.rxnarray[j,i]] += frts[i+index] corespeciesconsumptionrates[d.rxnarray[j,i]] += rrts[i+index] @@ -359,7 +359,7 @@ function processfluxes(sim::Simulation,corespcsinds,corerxninds,edgespcsinds,edg if any(d.rxnarray[:,i].>length(corespeciesconcentrations)) continue end - for j = 1:3 + for j = 1:4 if d.rxnarray[j,i] != 0 corespeciesconsumptionrates[d.rxnarray[j,i]] += frts[i] corespeciesproductionrates[d.rxnarray[j,i]] += rrts[i] @@ -367,7 +367,7 @@ function processfluxes(sim::Simulation,corespcsinds,corerxninds,edgespcsinds,edg break end end - for j = 4:6 + for j = 5:8 if d.rxnarray[j,i] != 0 corespeciesproductionrates[d.rxnarray[j,i]] += frts[i] corespeciesconsumptionrates[d.rxnarray[j,i]] += rrts[i] diff --git a/src/Phase.jl b/src/Phase.jl index 510fb98c..45e14a51 100644 --- a/src/Phase.jl +++ b/src/Phase.jl @@ -257,7 +257,7 @@ Broadcast.broadcastable(p::T) where {T<:AbstractPhase} = Ref(p) export broadcastable function getreactionindices(spcs,rxns) where {Q<:AbstractPhase} - arr = zeros(Int64,(6,length(rxns))) + arr = zeros(Int64,(8,length(rxns))) names = [spc.name for spc in spcs] for (i,rxn) in enumerate(rxns) inds = [findfirst(isequal(spc),spcs) for spc in rxn.reactants] @@ -268,7 +268,7 @@ function getreactionindices(spcs,rxns) where {Q<:AbstractPhase} end for (j,spc) in enumerate(rxn.products) ind = findfirst(isequal(spc),spcs) - arr[j+3,i] = ind + arr[j+4,i] = ind rxn.productinds[j] = ind end if hasproperty(rxn.kinetics,:efficiencies) && length(rxn.kinetics.nameefficiencies) > 0 diff --git a/src/PhaseState.jl b/src/PhaseState.jl index 31a5fdf0..d4beb4cb 100644 --- a/src/PhaseState.jl +++ b/src/PhaseState.jl @@ -111,6 +111,8 @@ export getDiffusiveRate @fastmath @inbounds dGrxn -= Gs[rxn.reactantinds[1]]+Gs[rxn.reactantinds[2]] elseif Nreact == 3 @fastmath @inbounds dGrxn -= Gs[rxn.reactantinds[1]]+Gs[rxn.reactantinds[2]]+Gs[rxn.reactantinds[3]] + elseif Nreact == 4 + @fastmath @inbounds dGrxn -= Gs[rxn.reactantinds[1]]+Gs[rxn.reactantinds[2]]+Gs[rxn.reactantinds[3]]+Gs[rxn.reactantinds[4]] end if Nprod == 1 @fastmath @inbounds dGrxn += Gs[rxn.productinds[1]] @@ -118,6 +120,8 @@ export getDiffusiveRate @fastmath @inbounds dGrxn += Gs[rxn.productinds[1]]+Gs[rxn.productinds[2]] elseif Nprod == 3 @fastmath @inbounds dGrxn += Gs[rxn.productinds[1]]+Gs[rxn.productinds[2]]+Gs[rxn.productinds[3]] + elseif Nprod == 4 + @fastmath @inbounds dGrxn += Gs[rxn.productinds[1]]+Gs[rxn.productinds[2]]+Gs[rxn.productinds[3]]+Gs[rxn.productinds[4]] end return @inbounds @fastmath exp(-(dGrxn+rxn.electronchange*phi)/(R*T))*(getC0(ph,T))^(Nprod-Nreact) end diff --git a/src/Reactor.jl b/src/Reactor.jl index ee66640d..1ca0042d 100644 --- a/src/Reactor.jl +++ b/src/Reactor.jl @@ -94,7 +94,7 @@ function Reactor(domains::T,y0s::W1,tspan::W2,interfaces::Z=Tuple(),ps::X=DiffEq for (j,domain) in enumerate(domains) Nspcs = length(domain.phase.species) Ntherm = length(domain.indexes) - 2 - for i = 1:6, j = 1:size(domain.rxnarray)[2] + for i = 1:8, j = 1:size(domain.rxnarray)[2] if domain.rxnarray[i,j] != 0 domain.rxnarray[i,j] += k-1 end @@ -293,6 +293,8 @@ end @fastmath @inbounds R += kfs[rxn.index]*cs[rxn.reactantinds[1]]*cs[rxn.reactantinds[2]] elseif Nreact == 3 @fastmath @inbounds R += kfs[rxn.index]*cs[rxn.reactantinds[1]]*cs[rxn.reactantinds[2]]*cs[rxn.reactantinds[3]] + elseif Nreact == 4 + @fastmath @inbounds R += kfs[rxn.index]*cs[rxn.reactantinds[1]]*cs[rxn.reactantinds[2]]*cs[rxn.reactantinds[3]]*cs[rxn.reactantinds[4]] end if Nprod == 1 @@ -301,6 +303,8 @@ end @fastmath @inbounds R -= krevs[rxn.index]*cs[rxn.productinds[1]]*cs[rxn.productinds[2]] elseif Nprod == 3 @fastmath @inbounds R -= krevs[rxn.index]*cs[rxn.productinds[1]]*cs[rxn.productinds[2]]*cs[rxn.productinds[3]] + elseif Nprod == 4 + @fastmath @inbounds R -= krevs[rxn.index]*cs[rxn.productinds[1]]*cs[rxn.productinds[2]]*cs[rxn.productinds[3]]*cs[rxn.productinds[4]] end return R @@ -338,15 +342,19 @@ export getrates @inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]] elseif @inbounds rarray[3,i] == 0 @inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]*cs[rarray[2,i]] - else + elseif @inbounds rarray[4,i] == 0 @inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]*cs[rarray[2,i]]*cs[rarray[3,i]] + else + @inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]*cs[rarray[2,i]]*cs[rarray[3,i]]*cs[rarray[4,i]] end - if @inbounds rarray[5,i] == 0 - @inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]] - elseif @inbounds rarray[6,i] == 0 - @inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]]*cs[rarray[5,i]] + if @inbounds rarray[6,i] == 0 + @inbounds @fastmath rR = krevs[i]*cs[rarray[5,i]] + elseif @inbounds rarray[7,i] == 0 + @inbounds @fastmath rR = krevs[i]*cs[rarray[5,i]]*cs[rarray[6,i]] + elseif @inbounds rarray[8,i] == 0 + @inbounds @fastmath rR = krevs[i]*cs[rarray[5,i]]*cs[rarray[6,i]]*cs[rarray[7,i]] else - @inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]]*cs[rarray[5,i]]*cs[rarray[6,i]] + @inbounds @fastmath rR = krevs[i]*cs[rarray[5,i]]*cs[rarray[6,i]]*cs[rarray[7,i]]*cs[rarray[8,i]] end @fastmath R = fR - rR @inbounds @fastmath dydt[rarray[1,i]] -= R @@ -354,13 +362,19 @@ export getrates @inbounds @fastmath dydt[rarray[2,i]] -= R if @inbounds rarray[3,i] != 0 @inbounds @fastmath dydt[rarray[3,i]] -= R + if @inbounds rarray[4,i] != 0 + @inbounds @fastmath dydt[rarray[4,i]] -= R + end end end - @inbounds @fastmath dydt[rarray[4,i]] += R - if @inbounds rarray[5,i] != 0 - @inbounds @fastmath dydt[rarray[5,i]] += R - if @inbounds rarray[6,i] != 0 - @inbounds @fastmath dydt[rarray[6,i]] += R + @inbounds @fastmath dydt[rarray[5,i]] += R + if @inbounds rarray[6,i] != 0 + @inbounds @fastmath dydt[rarray[6,i]] += R + if @inbounds rarray[7,i] != 0 + @inbounds @fastmath dydt[rarray[7,i]] += R + if @inbounds rarray[8,i] != 0 + @inbounds @fastmath dydt[rarray[8 ,i]] += R + end end end end @@ -410,15 +424,19 @@ export addreactionratecontributions! @inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]] elseif @inbounds rarray[3,i] == 0 @inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]*cs[rarray[2,i]] - else + elseif @inbounds rarray[4,i] == 0 @inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]*cs[rarray[2,i]]*cs[rarray[3,i]] + else + @inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]*cs[rarray[2,i]]*cs[rarray[3,i]]*cs[rarray[4,i]] end - if @inbounds rarray[5,i] == 0 - @inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]] - elseif @inbounds rarray[6,i] == 0 - @inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]]*cs[rarray[5,i]] + if @inbounds rarray[6,i] == 0 + @inbounds @fastmath rR = krevs[i]*cs[rarray[5,i]] + elseif @inbounds rarray[7,i] == 0 + @inbounds @fastmath rR = krevs[i]*cs[rarray[5,i]]*cs[rarray[6,i]] + elseif @inbounds rarray[8,i] == 0 + @inbounds @fastmath rR = krevs[i]*cs[rarray[5,i]]*cs[rarray[6,i]]*cs[rarray[7,i]] else - @inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]]*cs[rarray[5,i]]*cs[rarray[6,i]] + @inbounds @fastmath rR = krevs[i]*cs[rarray[5,i]]*cs[rarray[6,i]]*cs[rarray[7,i]]*cs[rarray[8,i]] end @inbounds @fastmath frts[i] = fR*V @inbounds @fastmath rrts[i] = rR*V @@ -429,13 +447,19 @@ export addreactionratecontributions! @inbounds @fastmath dydt[rarray[2,i]] -= R if @inbounds rarray[3,i] != 0 @inbounds @fastmath dydt[rarray[3,i]] -= R + if @inbounds rarray[4,i] != 0 + @inbounds @fastmath dydt[rarray[4,i]] -= R + end end end - @inbounds @fastmath dydt[rarray[4,i]] += R - if @inbounds rarray[5,i] != 0 - @inbounds @fastmath dydt[rarray[5,i]] += R - if @inbounds rarray[6,i] != 0 - @inbounds @fastmath dydt[rarray[6,i]] += R + @inbounds @fastmath dydt[rarray[5,i]] += R + if @inbounds rarray[6,i] != 0 + @inbounds @fastmath dydt[rarray[6,i]] += R + if @inbounds rarray[7,i] != 0 + @inbounds @fastmath dydt[rarray[7,i]] += R + if @inbounds rarray[8,i] != 0 + @inbounds @fastmath dydt[rarray[8,i]] += R + end end end end @@ -603,11 +627,14 @@ end export jacobianp @inline function _spreadreactantpartials!(jac::S,deriv::Float64,rxnarray::Array{Int64,2},rxnind::Int64,ind::Int64) where {S<:AbstractArray} - @inbounds jac[rxnarray[4,rxnind],ind] += deriv - if @inbounds rxnarray[5,rxnind] !== 0 - @inbounds jac[rxnarray[5,rxnind],ind] += deriv - if @inbounds rxnarray[6,rxnind] !== 0 - @inbounds jac[rxnarray[6,rxnind],ind] += deriv + @inbounds jac[rxnarray[5,rxnind],ind] += deriv + if @inbounds rxnarray[6,rxnind] !== 0 + @inbounds jac[rxnarray[6,rxnind],ind] += deriv + if @inbounds rxnarray[7,rxnind] !== 0 + @inbounds jac[rxnarray[7,rxnind],ind] += deriv + if @inbounds rxnarray[8,rxnind] !== 0 + @inbounds jac[rxnarray[8,rxnind],ind] += deriv + end end end end @@ -617,6 +644,9 @@ end @inbounds jac[rxnarray[2,rxnind],ind] += deriv if @inbounds rxnarray[3,rxnind] !== 0 @inbounds jac[rxnarray[3,rxnind],ind] += deriv + if @inbounds rxnarray[4,rxnind] !== 0 + @inbounds jac[rxnarray[4,rxnind],ind] += deriv + end end end end @@ -642,7 +672,7 @@ end @inbounds jac[rxnarray[2,rxnind],rxnarray[2,rxnind]] -= deriv @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[2,rxnind]) end - else + elseif rxnarray[4,rxnind] == 0 if rxnarray[1,rxnind]==rxnarray[2,rxnind] && rxnarray[1,rxnind]==rxnarray[3,rxnind] @inbounds @fastmath deriv = 3.0*k*cs[rxnarray[1,rxnind]]*cs[rxnarray[1,rxnind]] @inbounds jac[rxnarray[1,rxnind],rxnarray[1,rxnind]] -= 3.0*deriv @@ -691,75 +721,401 @@ end @inbounds jac[rxnarray[3,rxnind],rxnarray[3,rxnind]] -= deriv @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[3,rxnind]) end + else + if rxnarray[1,rxnind]==rxnarray[2,rxnind] && rxnarray[1,rxnind]==rxnarray[3,rxnind] && rxnarray[1,rxnind]==rxnarray[4,rxnind] + @inbounds @fastmath deriv = 4.0*k*cs[rxnarray[1,rxnind]]*cs[rxnarray[1,rxnind]]*cs[rxnarray[1,rxnind]] + @inbounds jac[rxnarray[1,rxnind],rxnarray[1,rxnind]] -= 4.0*deriv + @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[1,rxnind]) + elseif rxnarray[1,rxnind]==rxnarray[2,rxnind] && rxnarray[1,rxnind]==rxnarray[3,rxnind] + @inbounds @fastmath deriv = 3.0*k*cs[rxnarray[1,rxnind]]*cs[rxnarray[1,rxnind]]*cs[rxnarray[4,rxnind]] + @inbounds jac[rxnarray[1,rxnind],rxnarray[1,rxnind]] -= 3.0*deriv + @inbounds jac[rxnarray[4,rxnind],rxnarray[1,rxnind]] -= deriv + @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[1,rxnind]) + @inbounds @fastmath deriv = k*cs[rxnarray[1,rxnind]]*cs[rxnarray[1,rxnind]]*cs[rxnarray[1,rxnind]] + @inbounds jac[rxnarray[1,rxnind],rxnarray[4,rxnind]] -= 3.0*deriv + @inbounds jac[rxnarray[4,rxnind],rxnarray[4,rxnind]] -= deriv + @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[4,rxnind]) + elseif rxnarray[1,rxnind]==rxnarray[3,rxnind] && rxnarray[1,rxnind]==rxnarray[4,rxnind] + @inbounds @fastmath deriv = 3.0*k*cs[rxnarray[1,rxnind]]*cs[rxnarray[1,rxnind]]*cs[rxnarray[2,rxnind]] + @inbounds jac[rxnarray[1,rxnind],rxnarray[1,rxnind]] -= 3.0*deriv + @inbounds jac[rxnarray[2,rxnind],rxnarray[1,rxnind]] -= deriv + @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[1,rxnind]) + @inbounds @fastmath deriv = k*cs[rxnarray[1,rxnind]]*cs[rxnarray[1,rxnind]]*cs[rxnarray[1,rxnind]] + @inbounds jac[rxnarray[1,rxnind],rxnarray[2,rxnind]] -= 3.0*deriv + @inbounds jac[rxnarray[2,rxnind],rxnarray[2,rxnind]] -= deriv + @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[2,rxnind]) + elseif rxnarray[1,rxnind]==rxnarray[2,rxnind] && rxnarray[1,rxnind]==rxnarray[4,rxnind] + @inbounds @fastmath deriv = 3.0*k*cs[rxnarray[1,rxnind]]*cs[rxnarray[1,rxnind]]*cs[rxnarray[3,rxnind]] + @inbounds jac[rxnarray[1,rxnind],rxnarray[1,rxnind]] -= 3.0*deriv + @inbounds jac[rxnarray[3,rxnind],rxnarray[1,rxnind]] -= deriv + @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[1,rxnind]) + @inbounds @fastmath deriv = k*cs[rxnarray[1,rxnind]]*cs[rxnarray[1,rxnind]]*cs[rxnarray[1,rxnind]] + @inbounds jac[rxnarray[1,rxnind],rxnarray[3,rxnind]] -= 3.0*deriv + @inbounds jac[rxnarray[3,rxnind],rxnarray[3,rxnind]] -= deriv + @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[3,rxnind]) + elseif rxnarray[2,rxnind]==rxnarray[3,rxnind] && rxnarray[2,rxnind]==rxnarray[4,rxnind] + @inbounds @fastmath deriv = k*cs[rxnarray[2,rxnind]]*cs[rxnarray[2,rxnind]]*cs[rxnarray[2,rxnind]] + @inbounds jac[rxnarray[1,rxnind],rxnarray[1,rxnind]] -= deriv + @inbounds jac[rxnarray[2,rxnind],rxnarray[1,rxnind]] -= 3.0*deriv + @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[1,rxnind]) + @inbounds @fastmath deriv = 3.0*k*cs[rxnarray[1,rxnind]]*cs[rxnarray[2,rxnind]]*cs[rxnarray[2,rxnind]] + @inbounds jac[rxnarray[1,rxnind],rxnarray[2,rxnind]] -= deriv + @inbounds jac[rxnarray[2,rxnind],rxnarray[2,rxnind]] -= 3.0*deriv + @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[2,rxnind]) + elseif rxnarray[1,rxnind]==rxnarray[2,rxnind] + @inbounds @fastmath deriv = 2.0*k*cs[rxnarray[1,rxnind]]*cs[rxnarray[3,rxnind]]*cs[rxnarray[4,rxnind]] + @inbounds jac[rxnarray[1,rxnind],rxnarray[1,rxnind]] -= 2.0*deriv + @inbounds jac[rxnarray[3,rxnind],rxnarray[1,rxnind]] -= deriv + @inbounds jac[rxnarray[4,rxnind],rxnarray[1,rxnind]] -= deriv + @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[1,rxnind]) + @inbounds @fastmath deriv = k*cs[rxnarray[1,rxnind]]*cs[rxnarray[1,rxnind]]*cs[rxnarray[4,rxnind]] + @inbounds jac[rxnarray[1,rxnind],rxnarray[3,rxnind]] -= 2.0*deriv + @inbounds jac[rxnarray[3,rxnind],rxnarray[3,rxnind]] -= deriv + @inbounds jac[rxnarray[4,rxnind],rxnarray[3,rxnind]] -= deriv + @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[3,rxnind]) + @inbounds @fastmath deriv = k*cs[rxnarray[1,rxnind]]*cs[rxnarray[1,rxnind]]*cs[rxnarray[3,rxnind]] + @inbounds jac[rxnarray[1,rxnind],rxnarray[4,rxnind]] -= 2.0*deriv + @inbounds jac[rxnarray[3,rxnind],rxnarray[4,rxnind]] -= deriv + @inbounds jac[rxnarray[4,rxnind],rxnarray[4,rxnind]] -= deriv + @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[4,rxnind]) + elseif rxnarray[1,rxnind]==rxnarray[3,rxnind] + @inbounds @fastmath deriv = 2.0*k*cs[rxnarray[1,rxnind]]*cs[rxnarray[2,rxnind]]*cs[rxnarray[4,rxnind]] + @inbounds jac[rxnarray[1,rxnind],rxnarray[1,rxnind]] -= 2.0*deriv + @inbounds jac[rxnarray[2,rxnind],rxnarray[1,rxnind]] -= deriv + @inbounds jac[rxnarray[4,rxnind],rxnarray[1,rxnind]] -= deriv + @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[1,rxnind]) + @inbounds @fastmath deriv = k*cs[rxnarray[1,rxnind]]*cs[rxnarray[1,rxnind]]*cs[rxnarray[4,rxnind]] + @inbounds jac[rxnarray[1,rxnind],rxnarray[2,rxnind]] -= 2.0*deriv + @inbounds jac[rxnarray[2,rxnind],rxnarray[2,rxnind]] -= deriv + @inbounds jac[rxnarray[4,rxnind],rxnarray[2,rxnind]] -= deriv + @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[2,rxnind]) + @inbounds @fastmath deriv = k*cs[rxnarray[1,rxnind]]*cs[rxnarray[1,rxnind]]*cs[rxnarray[2,rxnind]] + @inbounds jac[rxnarray[1,rxnind],rxnarray[4,rxnind]] -= 2.0*deriv + @inbounds jac[rxnarray[2,rxnind],rxnarray[4,rxnind]] -= deriv + @inbounds jac[rxnarray[4,rxnind],rxnarray[4,rxnind]] -= deriv + @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[4,rxnind]) + elseif rxnarray[1,rxnind]==rxnarray[4,rxnind] + @inbounds @fastmath deriv = 2.0*k*cs[rxnarray[1,rxnind]]*cs[rxnarray[2,rxnind]]*cs[rxnarray[3,rxnind]] + @inbounds jac[rxnarray[1,rxnind],rxnarray[1,rxnind]] -= 2.0*deriv + @inbounds jac[rxnarray[2,rxnind],rxnarray[1,rxnind]] -= deriv + @inbounds jac[rxnarray[3,rxnind],rxnarray[1,rxnind]] -= deriv + @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[1,rxnind]) + @inbounds @fastmath deriv = k*cs[rxnarray[1,rxnind]]*cs[rxnarray[1,rxnind]]*cs[rxnarray[3,rxnind]] + @inbounds jac[rxnarray[1,rxnind],rxnarray[2,rxnind]] -= 2.0*deriv + @inbounds jac[rxnarray[2,rxnind],rxnarray[2,rxnind]] -= deriv + @inbounds jac[rxnarray[3,rxnind],rxnarray[2,rxnind]] -= deriv + @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[2,rxnind]) + @inbounds @fastmath deriv = k*cs[rxnarray[1,rxnind]]*cs[rxnarray[1,rxnind]]*cs[rxnarray[2,rxnind]] + @inbounds jac[rxnarray[1,rxnind],rxnarray[3,rxnind]] -= 2.0*deriv + @inbounds jac[rxnarray[2,rxnind],rxnarray[3,rxnind]] -= deriv + @inbounds jac[rxnarray[3,rxnind],rxnarray[3,rxnind]] -= deriv + @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[3,rxnind]) + elseif rxnarray[2,rxnind]==rxnarray[3,rxnind] + @inbounds @fastmath deriv = k*cs[rxnarray[2,rxnind]]*cs[rxnarray[2,rxnind]]*cs[rxnarray[4,rxnind]] + @inbounds jac[rxnarray[1,rxnind],rxnarray[1,rxnind]] -= deriv + @inbounds jac[rxnarray[2,rxnind],rxnarray[1,rxnind]] -= 2.0*dderiv + @inbounds jac[rxnarray[4,rxnind],rxnarray[1,rxnind]] -= deriv + @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[1,rxnind]) + @inbounds @fastmath deriv = 2.0*k*cs[rxnarray[1,rxnind]]*cs[rxnarray[2,rxnind]]*cs[rxnarray[4,rxnind]] + @inbounds jac[rxnarray[1,rxnind],rxnarray[2,rxnind]] -= deriv + @inbounds jac[rxnarray[2,rxnind],rxnarray[2,rxnind]] -= 2.0*deriv + @inbounds jac[rxnarray[4,rxnind],rxnarray[2,rxnind]] -= deriv + @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[2,rxnind]) + @inbounds @fastmath deriv = k*cs[rxnarray[1,rxnind]]*cs[rxnarray[2,rxnind]]*cs[rxnarray[2,rxnind]] + @inbounds jac[rxnarray[1,rxnind],rxnarray[4,rxnind]] -= deriv + @inbounds jac[rxnarray[2,rxnind],rxnarray[4,rxnind]] -= 2.0*deriv + @inbounds jac[rxnarray[4,rxnind],rxnarray[4,rxnind]] -= deriv + @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[4,rxnind]) + elseif rxnarray[2,rxnind]==rxnarray[4,rxnind] + @inbounds @fastmath deriv = k*cs[rxnarray[2,rxnind]]*cs[rxnarray[2,rxnind]]*cs[rxnarray[3,rxnind]] + @inbounds jac[rxnarray[1,rxnind],rxnarray[1,rxnind]] -= deriv + @inbounds jac[rxnarray[2,rxnind],rxnarray[1,rxnind]] -= 2.0*dderiv + @inbounds jac[rxnarray[3,rxnind],rxnarray[1,rxnind]] -= deriv + @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[1,rxnind]) + @inbounds @fastmath deriv = 2.0*k*cs[rxnarray[1,rxnind]]*cs[rxnarray[2,rxnind]]*cs[rxnarray[3,rxnind]] + @inbounds jac[rxnarray[1,rxnind],rxnarray[2,rxnind]] -= deriv + @inbounds jac[rxnarray[2,rxnind],rxnarray[2,rxnind]] -= 2.0*deriv + @inbounds jac[rxnarray[3,rxnind],rxnarray[2,rxnind]] -= deriv + @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[2,rxnind]) + @inbounds @fastmath deriv = k*cs[rxnarray[1,rxnind]]*cs[rxnarray[2,rxnind]]*cs[rxnarray[2,rxnind]] + @inbounds jac[rxnarray[1,rxnind],rxnarray[3,rxnind]] -= deriv + @inbounds jac[rxnarray[2,rxnind],rxnarray[3,rxnind]] -= 2.0*deriv + @inbounds jac[rxnarray[3,rxnind],rxnarray[3,rxnind]] -= deriv + @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[3,rxnind]) + elseif rxnarray[3,rxnind]==rxnarray[4,rxnind] + @inbounds @fastmath deriv = k*cs[rxnarray[2,rxnind]]*cs[rxnarray[3,rxnind]]*cs[rxnarray[3,rxnind]] + @inbounds jac[rxnarray[1,rxnind],rxnarray[1,rxnind]] -= deriv + @inbounds jac[rxnarray[2,rxnind],rxnarray[1,rxnind]] -= dderiv + @inbounds jac[rxnarray[3,rxnind],rxnarray[1,rxnind]] -= 2.0*deriv + @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[1,rxnind]) + @inbounds @fastmath deriv = k*cs[rxnarray[1,rxnind]]*cs[rxnarray[3,rxnind]]*cs[rxnarray[3,rxnind]] + @inbounds jac[rxnarray[1,rxnind],rxnarray[2,rxnind]] -= deriv + @inbounds jac[rxnarray[2,rxnind],rxnarray[2,rxnind]] -= deriv + @inbounds jac[rxnarray[3,rxnind],rxnarray[2,rxnind]] -= 2.0*deriv + @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[2,rxnind]) + @inbounds @fastmath deriv = 2.0*k*cs[rxnarray[1,rxnind]]*cs[rxnarray[2,rxnind]]*cs[rxnarray[3,rxnind]] + @inbounds jac[rxnarray[1,rxnind],rxnarray[3,rxnind]] -= deriv + @inbounds jac[rxnarray[2,rxnind],rxnarray[3,rxnind]] -= deriv + @inbounds jac[rxnarray[3,rxnind],rxnarray[3,rxnind]] -= 2.0*deriv + @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[3,rxnind]) + else + @inbounds @fastmath deriv = k*cs[rxnarray[2,rxnind]]*cs[rxnarray[3,rxnind]]*cs[rxnarray[4,rxnind]] + @inbounds jac[rxnarray[1,rxnind],rxnarray[1,rxnind]] -= deriv + @inbounds jac[rxnarray[2,rxnind],rxnarray[1,rxnind]] -= deriv + @inbounds jac[rxnarray[3,rxnind],rxnarray[1,rxnind]] -= deriv + @inbounds jac[rxnarray[4,rxnind],rxnarray[1,rxnind]] -= deriv + @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[1,rxnind]) + @inbounds @fastmath deriv = k*cs[rxnarray[1,rxnind]]*cs[rxnarray[3,rxnind]]*cs[rxnarray[4,rxnind]] + @inbounds jac[rxnarray[1,rxnind],rxnarray[2,rxnind]] -= deriv + @inbounds jac[rxnarray[2,rxnind],rxnarray[2,rxnind]] -= deriv + @inbounds jac[rxnarray[3,rxnind],rxnarray[2,rxnind]] -= deriv + @inbounds jac[rxnarray[4,rxnind],rxnarray[2,rxnind]] -= deriv + @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[2,rxnind]) + @inbounds @fastmath deriv = k*cs[rxnarray[1,rxnind]]*cs[rxnarray[2,rxnind]]*cs[rxnarray[4,rxnind]] + @inbounds jac[rxnarray[1,rxnind],rxnarray[3,rxnind]] -= deriv + @inbounds jac[rxnarray[2,rxnind],rxnarray[3,rxnind]] -= deriv + @inbounds jac[rxnarray[3,rxnind],rxnarray[3,rxnind]] -= deriv + @inbounds jac[rxnarray[4,rxnind],rxnarray[3,rxnind]] -= deriv + @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[3,rxnind]) + @inbounds @fastmath deriv = k*cs[rxnarray[1,rxnind]]*cs[rxnarray[2,rxnind]]*cs[rxnarray[3,rxnind]] + @inbounds jac[rxnarray[1,rxnind],rxnarray[4,rxnind]] -= deriv + @inbounds jac[rxnarray[2,rxnind],rxnarray[4,rxnind]] -= deriv + @inbounds jac[rxnarray[3,rxnind],rxnarray[4,rxnind]] -= deriv + @inbounds jac[rxnarray[4,rxnind],rxnarray[4,rxnind]] -= deriv + @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[4,rxnind]) + end end k=krev - if rxnarray[5,rxnind] == 0 + if rxnarray[6,rxnind] == 0 deriv = k - @inbounds jac[rxnarray[4,rxnind],rxnarray[4,rxnind]] -= deriv - @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[4,rxnind]) - elseif rxnarray[6,rxnind] == 0 - if rxnarray[4,rxnind] == rxnarray[5,rxnind] - @inbounds @fastmath deriv = 2.0*k*cs[rxnarray[4,rxnind]] - @inbounds jac[rxnarray[4,rxnind],rxnarray[4,rxnind]] -= 2.0*deriv - @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[4,rxnind]) + @inbounds jac[rxnarray[5,rxnind],rxnarray[5,rxnind]] -= deriv + @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[5,rxnind]) + elseif rxnarray[7,rxnind] == 0 + if rxnarray[5,rxnind] == rxnarray[6,rxnind] + @inbounds @fastmath deriv = 2.0*k*cs[rxnarray[5,rxnind]] + @inbounds jac[rxnarray[5,rxnind],rxnarray[5,rxnind]] -= 2.0*deriv + @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[5,rxnind]) else + @inbounds @fastmath deriv = k*cs[rxnarray[6,rxnind]] + @inbounds jac[rxnarray[5,rxnind],rxnarray[5,rxnind]] -= deriv + @inbounds jac[rxnarray[6,rxnind],rxnarray[5,rxnind]] -= deriv + @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[5,rxnind]) @inbounds @fastmath deriv = k*cs[rxnarray[5,rxnind]] - @inbounds jac[rxnarray[4,rxnind],rxnarray[4,rxnind]] -= deriv - @inbounds jac[rxnarray[5,rxnind],rxnarray[4,rxnind]] -= deriv - @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[4,rxnind]) - @inbounds @fastmath deriv = k*cs[rxnarray[4,rxnind]] - @inbounds jac[rxnarray[4,rxnind],rxnarray[5,rxnind]] -= deriv + @inbounds jac[rxnarray[5,rxnind],rxnarray[6,rxnind]] -= deriv + @inbounds jac[rxnarray[6,rxnind],rxnarray[6,rxnind]] -= deriv + @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[6,rxnind]) + end + elseif rxnarray[8,rxnind] == 0 + if rxnarray[5,rxnind]==rxnarray[6,rxnind] && rxnarray[5,rxnind]==rxnarray[7,rxnind] + @inbounds @fastmath deriv = 3.0*k*cs[rxnarray[5,rxnind]]*cs[rxnarray[5,rxnind]] + @inbounds jac[rxnarray[5,rxnind],rxnarray[5,rxnind]] -= 3.0*deriv + @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[5,rxnind]) + elseif rxnarray[5,rxnind]==rxnarray[6,rxnind] + @inbounds @fastmath deriv = 2.0*k*cs[rxnarray[5,rxnind]]*cs[rxnarray[7,rxnind]] + @inbounds jac[rxnarray[5,rxnind],rxnarray[5,rxnind]] -= 2.0*deriv + @inbounds jac[rxnarray[7,rxnind],rxnarray[5,rxnind]] -= deriv + @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[5,rxnind]) + @inbounds @fastmath deriv = k*cs[rxnarray[5,rxnind]]*cs[rxnarray[5,rxnind]] + @inbounds jac[rxnarray[5,rxnind],rxnarray[7,rxnind]] -= 2.0*deriv + @inbounds jac[rxnarray[7,rxnind],rxnarray[7,rxnind]] -= deriv + @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[7,rxnind]) + elseif rxnarray[6,rxnind]==rxnarray[7,rxnind] + @inbounds @fastmath deriv = k*cs[rxnarray[6,rxnind]]*cs[rxnarray[6,rxnind]] @inbounds jac[rxnarray[5,rxnind],rxnarray[5,rxnind]] -= deriv + @inbounds jac[rxnarray[6,rxnind],rxnarray[5,rxnind]] -= 2.0*deriv @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[5,rxnind]) + @inbounds @fastmath deriv = 2.0*k*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]] + @inbounds jac[rxnarray[5,rxnind],rxnarray[6,rxnind]] -= deriv + @inbounds jac[rxnarray[6,rxnind],rxnarray[6,rxnind]] -= 2.0*deriv + @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[6,rxnind]) + elseif rxnarray[5,rxnind]==rxnarray[7,rxnind] + @inbounds @fastmath deriv = 2.0*k*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]] + @inbounds jac[rxnarray[5,rxnind],rxnarray[5,rxnind]] -= 2.0*deriv + @inbounds jac[rxnarray[6,rxnind],rxnarray[5,rxnind]] -= deriv + @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[5,rxnind]) + @inbounds @fastmath deriv = k*cs[rxnarray[5,rxnind]]*cs[rxnarray[5,rxnind]] + @inbounds jac[rxnarray[5,rxnind],rxnarray[6,rxnind]] -= 2.0*deriv + @inbounds jac[rxnarray[6,rxnind],rxnarray[6,rxnind]] -= deriv + @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[6,rxnind]) + else + @inbounds @fastmath deriv = k*cs[rxnarray[6,rxnind]]*cs[rxnarray[7,rxnind]] + @inbounds jac[rxnarray[5,rxnind],rxnarray[5,rxnind]] -= deriv + @inbounds jac[rxnarray[6,rxnind],rxnarray[5,rxnind]] -= deriv + @inbounds jac[rxnarray[7,rxnind],rxnarray[5,rxnind]] -= deriv + @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[5,rxnind]) + @inbounds @fastmath deriv = k*cs[rxnarray[5,rxnind]]*cs[rxnarray[7,rxnind]] + @inbounds jac[rxnarray[5,rxnind],rxnarray[6,rxnind]] -= deriv + @inbounds jac[rxnarray[6,rxnind],rxnarray[6,rxnind]] -= deriv + @inbounds jac[rxnarray[7,rxnind],rxnarray[6,rxnind]] -= deriv + @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[6,rxnind]) + @inbounds @fastmath deriv = k*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]] + @inbounds jac[rxnarray[5,rxnind],rxnarray[7,rxnind]] -= deriv + @inbounds jac[rxnarray[6,rxnind],rxnarray[7,rxnind]] -= deriv + @inbounds jac[rxnarray[7,rxnind],rxnarray[7,rxnind]] -= deriv + @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[7,rxnind]) end else - if rxnarray[4,rxnind]==rxnarray[5,rxnind] && rxnarray[4,rxnind]==rxnarray[6,rxnind] - @inbounds @fastmath deriv = 3.0*k*cs[rxnarray[4,rxnind]]*cs[rxnarray[4,rxnind]] - @inbounds jac[rxnarray[4,rxnind],rxnarray[4,rxnind]] -= 3.0*deriv - @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[4,rxnind]) - elseif rxnarray[4,rxnind]==rxnarray[5,rxnind] - @inbounds @fastmath deriv = 2.0*k*cs[rxnarray[4,rxnind]]*cs[rxnarray[6,rxnind]] - @inbounds jac[rxnarray[4,rxnind],rxnarray[4,rxnind]] -= 2.0*deriv - @inbounds jac[rxnarray[6,rxnind],rxnarray[4,rxnind]] -= deriv - @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[4,rxnind]) - @inbounds @fastmath deriv = k*cs[rxnarray[4,rxnind]]*cs[rxnarray[4,rxnind]] - @inbounds jac[rxnarray[4,rxnind],rxnarray[6,rxnind]] -= 2.0*deriv + if rxnarray[5,rxnind]==rxnarray[6,rxnind] && rxnarray[5,rxnind]==rxnarray[7,rxnind] && rxnarray[5,rxnind]==rxnarray[8,rxnind] + @inbounds @fastmath deriv = 4.0*k*cs[rxnarray[5,rxnind]]*cs[rxnarray[5,rxnind]]*cs[rxnarray[5,rxnind]] + @inbounds jac[rxnarray[5,rxnind],rxnarray[5,rxnind]] -= 4.0*deriv + @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[5,rxnind]) + elseif rxnarray[5,rxnind]==rxnarray[6,rxnind] && rxnarray[5,rxnind]==rxnarray[7,rxnind] + @inbounds @fastmath deriv = 3.0*k*cs[rxnarray[5,rxnind]]*cs[rxnarray[5,rxnind]]*cs[rxnarray[8,rxnind]] + @inbounds jac[rxnarray[5,rxnind],rxnarray[5,rxnind]] -= 3.0*deriv + @inbounds jac[rxnarray[8,rxnind],rxnarray[5,rxnind]] -= deriv + @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[5,rxnind]) + @inbounds @fastmath deriv = k*cs[rxnarray[5,rxnind]]*cs[rxnarray[5,rxnind]]*cs[rxnarray[5,rxnind]] + @inbounds jac[rxnarray[5,rxnind],rxnarray[8,rxnind]] -= 3.0*deriv + @inbounds jac[rxnarray[8,rxnind],rxnarray[8,rxnind]] -= deriv + @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[8,rxnind]) + elseif rxnarray[5,rxnind]==rxnarray[7,rxnind] && rxnarray[5,rxnind]==rxnarray[8,rxnind] + @inbounds @fastmath deriv = 3.0*k*cs[rxnarray[5,rxnind]]*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]] + @inbounds jac[rxnarray[5,rxnind],rxnarray[5,rxnind]] -= 3.0*deriv + @inbounds jac[rxnarray[6,rxnind],rxnarray[5,rxnind]] -= deriv + @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[5,rxnind]) + @inbounds @fastmath deriv = k*cs[rxnarray[5,rxnind]]*cs[rxnarray[5,rxnind]]*cs[rxnarray[5,rxnind]] + @inbounds jac[rxnarray[5,rxnind],rxnarray[6,rxnind]] -= 3.0*deriv @inbounds jac[rxnarray[6,rxnind],rxnarray[6,rxnind]] -= deriv @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[6,rxnind]) + elseif rxnarray[5,rxnind]==rxnarray[6,rxnind] && rxnarray[5,rxnind]==rxnarray[8,rxnind] + @inbounds @fastmath deriv = 3.0*k*cs[rxnarray[5,rxnind]]*cs[rxnarray[5,rxnind]]*cs[rxnarray[7,rxnind]] + @inbounds jac[rxnarray[5,rxnind],rxnarray[5,rxnind]] -= 3.0*deriv + @inbounds jac[rxnarray[7,rxnind],rxnarray[5,rxnind]] -= deriv + @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[5,rxnind]) + @inbounds @fastmath deriv = k*cs[rxnarray[5,rxnind]]*cs[rxnarray[5,rxnind]]*cs[rxnarray[5,rxnind]] + @inbounds jac[rxnarray[5,rxnind],rxnarray[7,rxnind]] -= 3.0*deriv + @inbounds jac[rxnarray[7,rxnind],rxnarray[7,rxnind]] -= deriv + @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[7,rxnind]) + elseif rxnarray[6,rxnind]==rxnarray[7,rxnind] && rxnarray[6,rxnind]==rxnarray[8,rxnind] + @inbounds @fastmath deriv = k*cs[rxnarray[6,rxnind]]*cs[rxnarray[6,rxnind]]*cs[rxnarray[6,rxnind]] + @inbounds jac[rxnarray[5,rxnind],rxnarray[5,rxnind]] -= deriv + @inbounds jac[rxnarray[6,rxnind],rxnarray[5,rxnind]] -= 3.0*deriv + @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[5,rxnind]) + @inbounds @fastmath deriv = 3.0*k*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]]*cs[rxnarray[6,rxnind]] + @inbounds jac[rxnarray[5,rxnind],rxnarray[6,rxnind]] -= deriv + @inbounds jac[rxnarray[6,rxnind],rxnarray[6,rxnind]] -= 3.0*deriv + @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[6,rxnind]) elseif rxnarray[5,rxnind]==rxnarray[6,rxnind] - @inbounds @fastmath deriv = k*cs[rxnarray[5,rxnind]]*cs[rxnarray[5,rxnind]] - @inbounds jac[rxnarray[4,rxnind],rxnarray[4,rxnind]] -= deriv - @inbounds jac[rxnarray[5,rxnind],rxnarray[4,rxnind]] -= 2.0*deriv - @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[4,rxnind]) - @inbounds @fastmath deriv = 2.0*k*cs[rxnarray[4,rxnind]]*cs[rxnarray[5,rxnind]] - @inbounds jac[rxnarray[4,rxnind],rxnarray[5,rxnind]] -= deriv + @inbounds @fastmath deriv = 2.0*k*cs[rxnarray[5,rxnind]]*cs[rxnarray[7,rxnind]]*cs[rxnarray[8,rxnind]] + @inbounds jac[rxnarray[5,rxnind],rxnarray[5,rxnind]] -= 2.0*deriv + @inbounds jac[rxnarray[7,rxnind],rxnarray[5,rxnind]] -= deriv + @inbounds jac[rxnarray[8,rxnind],rxnarray[5,rxnind]] -= deriv + @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[5,rxnind]) + @inbounds @fastmath deriv = k*cs[rxnarray[5,rxnind]]*cs[rxnarray[5,rxnind]]*cs[rxnarray[8,rxnind]] + @inbounds jac[rxnarray[5,rxnind],rxnarray[7,rxnind]] -= 2.0*deriv + @inbounds jac[rxnarray[7,rxnind],rxnarray[7,rxnind]] -= deriv + @inbounds jac[rxnarray[8,rxnind],rxnarray[7,rxnind]] -= deriv + @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[7,rxnind]) + @inbounds @fastmath deriv = k*cs[rxnarray[5,rxnind]]*cs[rxnarray[5,rxnind]]*cs[rxnarray[7,rxnind]] + @inbounds jac[rxnarray[5,rxnind],rxnarray[8,rxnind]] -= 2.0*deriv + @inbounds jac[rxnarray[7,rxnind],rxnarray[8,rxnind]] -= deriv + @inbounds jac[rxnarray[8,rxnind],rxnarray[8,rxnind]] -= deriv + @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[8,rxnind]) + elseif rxnarray[5,rxnind]==rxnarray[7,rxnind] + @inbounds @fastmath deriv = 2.0*k*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]]*cs[rxnarray[8,rxnind]] + @inbounds jac[rxnarray[5,rxnind],rxnarray[5,rxnind]] -= 2.0*deriv + @inbounds jac[rxnarray[6,rxnind],rxnarray[5,rxnind]] -= deriv + @inbounds jac[rxnarray[8,rxnind],rxnarray[5,rxnind]] -= deriv + @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[5,rxnind]) + @inbounds @fastmath deriv = k*cs[rxnarray[5,rxnind]]*cs[rxnarray[5,rxnind]]*cs[rxnarray[8,rxnind]] + @inbounds jac[rxnarray[5,rxnind],rxnarray[6,rxnind]] -= 2.0*deriv + @inbounds jac[rxnarray[6,rxnind],rxnarray[6,rxnind]] -= deriv + @inbounds jac[rxnarray[8,rxnind],rxnarray[6,rxnind]] -= deriv + @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[6,rxnind]) + @inbounds @fastmath deriv = k*cs[rxnarray[5,rxnind]]*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]] + @inbounds jac[rxnarray[5,rxnind],rxnarray[8,rxnind]] -= 2.0*deriv + @inbounds jac[rxnarray[6,rxnind],rxnarray[8,rxnind]] -= deriv + @inbounds jac[rxnarray[8,rxnind],rxnarray[8,rxnind]] -= deriv + @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[8,rxnind]) + elseif rxnarray[5,rxnind]==rxnarray[8,rxnind] + @inbounds @fastmath deriv = 2.0*k*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]]*cs[rxnarray[7,rxnind]] @inbounds jac[rxnarray[5,rxnind],rxnarray[5,rxnind]] -= 2.0*deriv + @inbounds jac[rxnarray[6,rxnind],rxnarray[5,rxnind]] -= deriv + @inbounds jac[rxnarray[7,rxnind],rxnarray[5,rxnind]] -= deriv @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[5,rxnind]) - elseif rxnarray[4,rxnind]==rxnarray[6,rxnind] - @inbounds @fastmath deriv = 2.0*k*cs[rxnarray[4,rxnind]]*cs[rxnarray[5,rxnind]] - @inbounds jac[rxnarray[4,rxnind],rxnarray[4,rxnind]] -= 2.0*deriv - @inbounds jac[rxnarray[5,rxnind],rxnarray[4,rxnind]] -= deriv - @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[4,rxnind]) - @inbounds @fastmath deriv = k*cs[rxnarray[4,rxnind]]*cs[rxnarray[4,rxnind]] - @inbounds jac[rxnarray[4,rxnind],rxnarray[5,rxnind]] -= 2.0*deriv + @inbounds @fastmath deriv = k*cs[rxnarray[5,rxnind]]*cs[rxnarray[5,rxnind]]*cs[rxnarray[7,rxnind]] + @inbounds jac[rxnarray[5,rxnind],rxnarray[6,rxnind]] -= 2.0*deriv + @inbounds jac[rxnarray[6,rxnind],rxnarray[6,rxnind]] -= deriv + @inbounds jac[rxnarray[7,rxnind],rxnarray[6,rxnind]] -= deriv + @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[6,rxnind]) + @inbounds @fastmath deriv = k*cs[rxnarray[5,rxnind]]*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]] + @inbounds jac[rxnarray[5,rxnind],rxnarray[7,rxnind]] -= 2.0*deriv + @inbounds jac[rxnarray[6,rxnind],rxnarray[7,rxnind]] -= deriv + @inbounds jac[rxnarray[7,rxnind],rxnarray[7,rxnind]] -= deriv + @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[7,rxnind]) + elseif rxnarray[6,rxnind]==rxnarray[7,rxnind] + @inbounds @fastmath deriv = k*cs[rxnarray[6,rxnind]]*cs[rxnarray[6,rxnind]]*cs[rxnarray[8,rxnind]] @inbounds jac[rxnarray[5,rxnind],rxnarray[5,rxnind]] -= deriv + @inbounds jac[rxnarray[6,rxnind],rxnarray[5,rxnind]] -= 2.0*deriv + @inbounds jac[rxnarray[8,rxnind],rxnarray[5,rxnind]] -= deriv @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[5,rxnind]) + @inbounds @fastmath deriv = 2.0*k*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]]*cs[rxnarray[8,rxnind]] + @inbounds jac[rxnarray[5,rxnind],rxnarray[6,rxnind]] -= deriv + @inbounds jac[rxnarray[6,rxnind],rxnarray[6,rxnind]] -= 2.0*deriv + @inbounds jac[rxnarray[8,rxnind],rxnarray[6,rxnind]] -= deriv + @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[6,rxnind]) + @inbounds @fastmath deriv = k*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]]*cs[rxnarray[6,rxnind]] + @inbounds jac[rxnarray[5,rxnind],rxnarray[8,rxnind]] -= deriv + @inbounds jac[rxnarray[6,rxnind],rxnarray[8,rxnind]] -= 2.0*deriv + @inbounds jac[rxnarray[8,rxnind],rxnarray[8,rxnind]] -= deriv + @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[8,rxnind]) + elseif rxnarray[6,rxnind]==rxnarray[8,rxnind] + @inbounds @fastmath deriv = k*cs[rxnarray[6,rxnind]]*cs[rxnarray[6,rxnind]]*cs[rxnarray[7,rxnind]] + @inbounds jac[rxnarray[5,rxnind],rxnarray[5,rxnind]] -= deriv + @inbounds jac[rxnarray[6,rxnind],rxnarray[5,rxnind]] -= 2.0*deriv + @inbounds jac[rxnarray[7,rxnind],rxnarray[5,rxnind]] -= deriv + @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[5,rxnind]) + @inbounds @fastmath deriv = 2.0*k*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]]*cs[rxnarray[7,rxnind]] + @inbounds jac[rxnarray[5,rxnind],rxnarray[6,rxnind]] -= deriv + @inbounds jac[rxnarray[6,rxnind],rxnarray[6,rxnind]] -= 2.0*deriv + @inbounds jac[rxnarray[7,rxnind],rxnarray[6,rxnind]] -= deriv + @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[6,rxnind]) + @inbounds @fastmath deriv = k*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]]*cs[rxnarray[6,rxnind]] + @inbounds jac[rxnarray[5,rxnind],rxnarray[7,rxnind]] -= deriv + @inbounds jac[rxnarray[6,rxnind],rxnarray[7,rxnind]] -= 2.0*deriv + @inbounds jac[rxnarray[7,rxnind],rxnarray[7,rxnind]] -= deriv + @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[7,rxnind]) + elseif rxnarray[7,rxnind]==rxnarray[8,rxnind] + @inbounds @fastmath deriv = k*cs[rxnarray[6,rxnind]]*cs[rxnarray[7,rxnind]]*cs[rxnarray[7,rxnind]] + @inbounds jac[rxnarray[5,rxnind],rxnarray[5,rxnind]] -= deriv + @inbounds jac[rxnarray[6,rxnind],rxnarray[5,rxnind]] -= deriv + @inbounds jac[rxnarray[7,rxnind],rxnarray[5,rxnind]] -= 2.0*deriv + @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[5,rxnind]) + @inbounds @fastmath deriv = k*cs[rxnarray[5,rxnind]]*cs[rxnarray[7,rxnind]]*cs[rxnarray[7,rxnind]] + @inbounds jac[rxnarray[5,rxnind],rxnarray[6,rxnind]] -= deriv + @inbounds jac[rxnarray[6,rxnind],rxnarray[6,rxnind]] -= deriv + @inbounds jac[rxnarray[7,rxnind],rxnarray[6,rxnind]] -= 2.0*deriv + @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[6,rxnind]) + @inbounds @fastmath deriv = 2.0*k*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]]*cs[rxnarray[7,rxnind]] + @inbounds jac[rxnarray[5,rxnind],rxnarray[7,rxnind]] -= deriv + @inbounds jac[rxnarray[6,rxnind],rxnarray[7,rxnind]] -= deriv + @inbounds jac[rxnarray[7,rxnind],rxnarray[7,rxnind]] -= 2.0*deriv + @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[7,rxnind]) else - @inbounds @fastmath deriv = k*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]] - @inbounds jac[rxnarray[4,rxnind],rxnarray[4,rxnind]] -= deriv - @inbounds jac[rxnarray[5,rxnind],rxnarray[4,rxnind]] -= deriv - @inbounds jac[rxnarray[6,rxnind],rxnarray[4,rxnind]] -= deriv - @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[4,rxnind]) - @inbounds @fastmath deriv = k*cs[rxnarray[4,rxnind]]*cs[rxnarray[6,rxnind]] - @inbounds jac[rxnarray[4,rxnind],rxnarray[5,rxnind]] -= deriv + @inbounds @fastmath deriv = k*cs[rxnarray[6,rxnind]]*cs[rxnarray[7,rxnind]]*cs[rxnarray[8,rxnind]] @inbounds jac[rxnarray[5,rxnind],rxnarray[5,rxnind]] -= deriv @inbounds jac[rxnarray[6,rxnind],rxnarray[5,rxnind]] -= deriv + @inbounds jac[rxnarray[7,rxnind],rxnarray[5,rxnind]] -= deriv + @inbounds jac[rxnarray[8,rxnind],rxnarray[5,rxnind]] -= deriv @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[5,rxnind]) - @inbounds @fastmath deriv = k*cs[rxnarray[4,rxnind]]*cs[rxnarray[5,rxnind]] - @inbounds jac[rxnarray[4,rxnind],rxnarray[6,rxnind]] -= deriv + @inbounds @fastmath deriv = k*cs[rxnarray[5,rxnind]]*cs[rxnarray[7,rxnind]]*cs[rxnarray[8,rxnind]] @inbounds jac[rxnarray[5,rxnind],rxnarray[6,rxnind]] -= deriv @inbounds jac[rxnarray[6,rxnind],rxnarray[6,rxnind]] -= deriv + @inbounds jac[rxnarray[7,rxnind],rxnarray[6,rxnind]] -= deriv + @inbounds jac[rxnarray[8,rxnind],rxnarray[6,rxnind]] -= deriv @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[6,rxnind]) + @inbounds @fastmath deriv = k*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]]*cs[rxnarray[8,rxnind]] + @inbounds jac[rxnarray[5,rxnind],rxnarray[7,rxnind]] -= deriv + @inbounds jac[rxnarray[6,rxnind],rxnarray[7,rxnind]] -= deriv + @inbounds jac[rxnarray[7,rxnind],rxnarray[7,rxnind]] -= deriv + @inbounds jac[rxnarray[8,rxnind],rxnarray[7,rxnind]] -= deriv + @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[7,rxnind]) + @inbounds @fastmath deriv = k*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]]*cs[rxnarray[7,rxnind]] + @inbounds jac[rxnarray[5,rxnind],rxnarray[8,rxnind]] -= deriv + @inbounds jac[rxnarray[6,rxnind],rxnarray[8,rxnind]] -= deriv + @inbounds jac[rxnarray[7,rxnind],rxnarray[8,rxnind]] -= deriv + @inbounds jac[rxnarray[8,rxnind],rxnarray[8,rxnind]] -= deriv + @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[8,rxnind]) end end end @@ -773,27 +1129,41 @@ end @inbounds jac[rxnarray[1,rxnind],Vind] -= deriv @inbounds jac[rxnarray[2,rxnind],Vind] -= deriv _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,Vind) - else + elseif rxnarray[4,rxnind]== 0 @inbounds @fastmath deriv = -2.0*k*cs[rxnarray[1,rxnind]]*cs[rxnarray[2,rxnind]]*cs[rxnarray[3,rxnind]] @inbounds jac[rxnarray[1,rxnind],Vind] -= deriv @inbounds jac[rxnarray[2,rxnind],Vind] -= deriv @inbounds jac[rxnarray[3,rxnind],Vind] -= deriv _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,Vind) + else + @inbounds @fastmath deriv = -3.0*k*cs[rxnarray[1,rxnind]]*cs[rxnarray[2,rxnind]]*cs[rxnarray[3,rxnind]]*cs[rxnarray[4,rxnind]] + @inbounds jac[rxnarray[1,rxnind],Vind] -= deriv + @inbounds jac[rxnarray[2,rxnind],Vind] -= deriv + @inbounds jac[rxnarray[3,rxnind],Vind] -= deriv + @inbounds jac[rxnarray[4,rxnind],Vind] -= deriv + _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,Vind) end k=krev - if rxnarray[5,rxnind]==0 + if rxnarray[6,rxnind]==0 nothing - elseif rxnarray[6,rxnind] == 0 - @inbounds @fastmath deriv = -k*cs[rxnarray[4,rxnind]]*cs[rxnarray[5,rxnind]] - @inbounds jac[rxnarray[4,rxnind],Vind] -= deriv + elseif rxnarray[7,rxnind] == 0 + @inbounds @fastmath deriv = -k*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]] @inbounds jac[rxnarray[5,rxnind],Vind] -= deriv + @inbounds jac[rxnarray[6,rxnind],Vind] -= deriv _spreadproductpartials!(jac,deriv,rxnarray,rxnind,Vind) - else - @inbounds @fastmath deriv = -2.0*k*cs[rxnarray[4,rxnind]]*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]] - @inbounds jac[rxnarray[4,rxnind],Vind] -= deriv + elseif rxnarray[8,rxnind] == 0 + @inbounds @fastmath deriv = -2.0*k*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]]*cs[rxnarray[7,rxnind]] @inbounds jac[rxnarray[5,rxnind],Vind] -= deriv @inbounds jac[rxnarray[6,rxnind],Vind] -= deriv + @inbounds jac[rxnarray[7,rxnind],Vind] -= deriv _spreadproductpartials!(jac,deriv,rxnarray,rxnind,Vind) + else + @inbounds @fastmath deriv = -3.0*k*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]]*cs[rxnarray[7,rxnind]]*cs[rxnarray[8,rxnind]] + @inbounds jac[rxnarray[5,rxnind],Vind] -= deriv + @inbounds jac[rxnarray[6,rxnind],Vind] -= deriv + @inbounds jac[rxnarray[7,rxnind],Vind] -= deriv + @inbounds jac[rxnarray[8,rxnind],Vind] -= deriv + _spreadproductpartials!(jac,deriv,rxnarray,rxnind,Vind) end end @@ -811,28 +1181,42 @@ This function calculates the ns partials in jacobiany involving k derivatives. d @inbounds jac[rxnarray[1,rxnind],xind] -= deriv @inbounds jac[rxnarray[2,rxnind],xind] -= deriv _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,xind) - else + elseif rxnarray[4,rxnind] == 0 @inbounds @fastmath deriv = dkdx*cs[rxnarray[1,rxnind]]*cs[rxnarray[2,rxnind]]*cs[rxnarray[3,rxnind]]*V @inbounds jac[rxnarray[1,rxnind],xind] -= deriv @inbounds jac[rxnarray[2,rxnind],xind] -= deriv @inbounds jac[rxnarray[3,rxnind],xind] -= deriv _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,xind) + else + @inbounds @fastmath deriv = dkdx*cs[rxnarray[1,rxnind]]*cs[rxnarray[2,rxnind]]*cs[rxnarray[3,rxnind]]*cs[rxnarray[4,rxnind]]*V + @inbounds jac[rxnarray[1,rxnind],xind] -= deriv + @inbounds jac[rxnarray[2,rxnind],xind] -= deriv + @inbounds jac[rxnarray[3,rxnind],xind] -= deriv + @inbounds jac[rxnarray[4,rxnind],xind] -= deriv + _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,xind) end dkdx = dkrevdx - if rxnarray[5,rxnind] == 0 - @inbounds @fastmath deriv = dkdx*cs[rxnarray[4,rxnind]]*V - @inbounds jac[rxnarray[4,rxnind],xind] -= deriv + if rxnarray[6,rxnind] == 0 + @inbounds @fastmath deriv = dkdx*cs[rxnarray[5,rxnind]]*V + @inbounds jac[rxnarray[5,rxnind],xind] -= deriv _spreadproductpartials!(jac,deriv,rxnarray,rxnind,xind) - elseif rxnarray[6,rxnind] == 0 - @inbounds @fastmath deriv = dkdx*cs[rxnarray[4,rxnind]]*cs[rxnarray[5,rxnind]]*V - @inbounds jac[rxnarray[4,rxnind],xind] -= deriv + elseif rxnarray[7,rxnind] == 0 + @inbounds @fastmath deriv = dkdx*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]]*V + @inbounds jac[rxnarray[5,rxnind],xind] -= deriv + @inbounds jac[rxnarray[6,rxnind],xind] -= deriv + _spreadproductpartials!(jac,deriv,rxnarray,rxnind,xind) + elseif rxnarray[8,rxnind] == 0 + @inbounds @fastmath deriv = dkdx*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]]*cs[rxnarray[7,rxnind]]*V @inbounds jac[rxnarray[5,rxnind],xind] -= deriv + @inbounds jac[rxnarray[6,rxnind],xind] -= deriv + @inbounds jac[rxnarray[7,rxnind],xind] -= deriv _spreadproductpartials!(jac,deriv,rxnarray,rxnind,xind) else - @inbounds @fastmath deriv = dkdx*cs[rxnarray[4,rxnind]]*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]]*V - @inbounds jac[rxnarray[4,rxnind],xind] -= deriv + @inbounds @fastmath deriv = dkdx*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]]*cs[rxnarray[7,rxnind]]*cs[rxnarray[8,rxnind]]*V @inbounds jac[rxnarray[5,rxnind],xind] -= deriv @inbounds jac[rxnarray[6,rxnind],xind] -= deriv + @inbounds jac[rxnarray[7,rxnind],xind] -= deriv + @inbounds jac[rxnarray[8,rxnind],xind] -= deriv _spreadproductpartials!(jac,deriv,rxnarray,rxnind,xind) end end diff --git a/src/Simulation.jl b/src/Simulation.jl index eeb87e36..48fc44e3 100644 --- a/src/Simulation.jl +++ b/src/Simulation.jl @@ -304,15 +304,19 @@ function rops!(ropmat,rarray,cs,kfs,krevs,V,start) @inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]] elseif @inbounds rarray[3,i] == 0 @inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]*cs[rarray[2,i]] - else + elseif @inbounds rarray[4,i] == 0 @inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]*cs[rarray[2,i]]*cs[rarray[3,i]] + else + @inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]*cs[rarray[2,i]]*cs[rarray[3,i]]*cs[rarray[4,i]] end - if @inbounds rarray[5,i] == 0 - @inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]] - elseif @inbounds rarray[6,i] == 0 - @inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]]*cs[rarray[5,i]] + if @inbounds rarray[6,i] == 0 + @inbounds @fastmath rR = krevs[i]*cs[rarray[5,i]] + elseif @inbounds rarray[7,i] == 0 + @inbounds @fastmath rR = krevs[i]*cs[rarray[5,i]]*cs[rarray[6,i]] + elseif rarray[8,i] == 0 + @inbounds @fastmath rR = krevs[i]*cs[rarray[5,i]]*cs[rarray[6,i]]*cs[rarray[7,i]] else - @inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]]*cs[rarray[5,i]]*cs[rarray[6,i]] + @inbounds @fastmath rR = krevs[i]*cs[rarray[5,i]]*cs[rarray[6,i]]*cs[rarray[7,i]]*cs[rarray[8,i]] end @fastmath R = (fR - rR)*V @@ -321,13 +325,19 @@ function rops!(ropmat,rarray,cs,kfs,krevs,V,start) @inbounds @fastmath ropmat[i+start,rarray[2,i]] -= R if @inbounds rarray[3,i] != 0 @inbounds @fastmath ropmat[i+start,rarray[3,i]] -= R + if @inbounds rarray[4,i] != 0 + @inbounds @fastmath ropmat[i+start,rarray[4,i]] -= R + end end end - @inbounds @fastmath ropmat[i+start,rarray[4,i]] += R - if @inbounds rarray[5,i] != 0 - @inbounds @fastmath ropmat[i+start,rarray[5,i]] += R - if @inbounds rarray[6,i] != 0 - @inbounds @fastmath ropmat[i+start,rarray[6,i]] += R + @inbounds @fastmath ropmat[i+start,rarray[5,i]] += R + if @inbounds rarray[6,i] != 0 + @inbounds @fastmath ropmat[i+start,rarray[6,i]] += R + if @inbounds rarray[7,i] != 0 + @inbounds @fastmath ropmat[i+start,rarray[7,i]] += R + if @inbounds rarray[8,i] != 0 + @inbounds @fastmath ropmat[i+start,rarray[8,i]] += R + end end end end @@ -335,21 +345,25 @@ end function rops!(ropvec,rarray,cs,kfs,krevs,V,start,ind) for i = 1:length(kfs) - c = count(isequal(ind),rarray[4:6,i])-count(isequal(ind),rarray[1:3,i]) + c = count(isequal(ind),rarray[5:8,i])-count(isequal(ind),rarray[1:4,i]) if c != 0.0 if @inbounds rarray[2,i] == 0 @inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]] elseif @inbounds rarray[3,i] == 0 @inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]*cs[rarray[2,i]] - else + elseif @inbounds rarray[4,i] == 0 @inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]*cs[rarray[2,i]]*cs[rarray[3,i]] + else + @inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]*cs[rarray[2,i]]*cs[rarray[3,i]]*cs[rarray[4,i]] end - if @inbounds rarray[5,i] == 0 - @inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]] - elseif @inbounds rarray[6,i] == 0 - @inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]]*cs[rarray[5,i]] + if @inbounds rarray[6,i] == 0 + @inbounds @fastmath rR = krevs[i]*cs[rarray[5,i]] + elseif @inbounds rarray[7,i] == 0 + @inbounds @fastmath rR = krevs[i]*cs[rarray[5,i]]*cs[rarray[6,i]] + elseif rarray[8,i] == 0 + @inbounds @fastmath rR = krevs[i]*cs[rarray[5,i]]*cs[rarray[6,i]]*cs[rarray[7,i]] else - @inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]]*cs[rarray[5,i]]*cs[rarray[6,i]] + @inbounds @fastmath rR = krevs[i]*cs[rarray[5,i]]*cs[rarray[6,i]]*cs[rarray[7,i]]*cs[rarray[8,i]] end @fastmath R = (fR - rR)*V @fastmath @inbounds ropvec[i+start] = c*R