Stoichiometric Reactions

[mjohnson541]

Currently RMS doesn't support reactions that can't be represented in elementary form or with more than 4 reactants or products. I got an email asking about support for specifying stochiometric coefficients.

I think I properly laid out what needs done for the derivative functions:

In order to simulate non-elementary reactions (non-integer stoichiometric coefficients) the simplest thing to do I think is to write a new class of AbstractReaction, probably something like StoichiometricReaction in Reaction.jl with the same attributes as ElementaryReaction plus the stoichiometric coefficients. When the phase objects are constructed these reactions need to be separated out before processing and stored in a separate attribute on each phase object in Phase.jl. Then add a function analogous to addreactionratecontributions! in Reactor.jl that takes the list of stoichiometric reactions instead of the rarray (which is a list of the indices of reactants (first 3) and products(last 3) for each elementary reaction), it should end up looking something like fR = cs[rxn.reactantinds[1]]^rxn.reactantstoich[1]*cs[rxn.reactantinds[2]]^rxn.reactantstoich[2] and so on. Then a call to that function needs to be added within the two dydtreactor! functions next to the addreactionratecontributions! functions in Reactor.jl you can get the phase object from domain.phase.

But changes also need made to the analytic jacobians and downstream processing (Simulation.jl). I think the downstream processing should be relatively straightforward, but @hwpang what do you think it would take for the analytic jacobians to support stochiometric reactions?

[hwpang]

I think we will need something analogous to the _jacobianynswrtns!, and put it inside jacobianynsderiv!. We will need a loop inside or outside this function to loop through the list of StoichiometricReaction, and I imagine we would need to loop though the reactants/products inside the function to calculate the derivative, something like

@inline function _jacobianynswrtns!(jac::S,rxn::StoichiometricReaction,rxnind::Int64,cs::Array{Float64,1},kf::Float64,krev::Float64) where {S<:AbstractArray}
    k=kf
    for i in 1:length(rxn.reactantinds)
		for j in 1:length(rxn.reactantinds)
			deriv = k * sum(ind == j ? rxn.reactantstoich[ind]*cs[ind]^(rxn.reactantstoich[ind]-1) : cs[ind]^rxn.reactantstoich[ind]))
			@inbounds jac[i,j] -= deriv
		end
	end

	k=krev
    for i in 1:length(rxn.productinds)
		for j in 1:length(rxn.productinds)
			deriv = k * sum(ind == j ? rxn.productstoich[ind]*cs[ind]^(rxn.productstoich[ind]-1) : cs[ind]^rxn.productstoich[ind]))
			@inbounds jac[i,j] += deriv
		end
	end
end