Error in the implementation of chemistryModel::jacobian?
 

Hi,

I had a look at the function calculating the jacobian for the ODE-solvers used for the chemistry in reactingFoam. It's called ODEChemistryModel::jacobian in 2.0.x (which I use) and chemistryModel::jacobian in 2.2.2 (both versions have identical implementation) and also defined in ODE.H in both versions. One of the outputs (in form of a reference to a scalarField) is supposed to be the dertivative of the RHS with respect to the independent variable (dfdx) of a general ODE system:

dydx = f(y(x),x)

This seems to be the case in the solving classes under /src/ODE/, from what I can see how they are using it in the ODE-solution methods. However, in it's implementation (in ODEChemistryModel.C or chemistryModel.C depending on version) it is put equal to the total reaction rate for each species:

template<class CompType, class ThermoType>
void Foam::ODEChemistryModel<CompType, ThermoType>::jacobian
(
const scalar t,
const scalarField& c,
scalarField& dcdt,
scalarSquareMatrix& dfdc
) const
{
....

scalarField c2(nSpecie_, 0.0);
forAll(c2, i)
{
c2[i] = max(c[i], 0.0);
}

......

dcdt = omega(c2, T, p);

......
}

In the case of the ODE-system describing the evolution of species concentrations with time, the total reaction rate for each species are analogous to rhe RHS (f) in the general ODE-system above. So, to me at least, it seems that ::jacobian calculates the RHS again instead when it is supposed to calculate the time derivative of the RHS (i.e. the time derivative of the total reaction rate for each species). There, is also already a function ::derivatives that calculates the RHS, which the solvers in /src/ODE/ use.

Have I misunderstood something?:confused: