Jacobian in Cantera

The Jacobian matrix is constructed in Cantera by the one-dimensional domain.

The hierarchy of domains/points/species is described in one of Dave Goodwin's presentations (available at file http://www.et.byu.edu/~tom/classes/641/Cantera/Workshop/flames.pdf )

Background Info

Jacobian Matrix

The Jacobian matrix of the chemical source terms looks like this:

$J=\left|{\begin{array}{ccc}{\frac {\partial {\dot {y}}_{1}}{\partial y_{1}}}&\cdots &{\frac {\partial {\dot {y}}_{n}}{\partial y_{n}}}\\\cdots &\cdots &\cdots \\{\frac {\partial {\dot {y}}_{n}}{\partial y_{1}}}&\cdots &{\frac {\partial {\dot {y}}_{n}}{\partial y_{n}}}\end{array}}\right|$

Cantera Classes

Entire Domain Objects

OneDim:

• located in /path/to/cantera/Cantera/src/oneD/OneDim.cpp
• contains multiple 1D domains
• used to solve multiple domain problems

Sim1D:

• located in /path/to/cantera/Cantera/src/oneD/Sim1D.h
• extended version of OneDim class (stores more information)

Single Domain Objects

Domain1D:

• located in /path/to/cantera/Cantera/src/oneD/Domain1D.cpp
• represents a single 1D domain
• actually creates the Jacobian for each point in the domain
• this class is a base class

Domain1D is a base class extended by the following classes:

• StFlow - 1D flow domain that satisfies 1D similarity solution for reacting,axisymmetric flows
  • located in /path/to/cantera/Cantera/src/oneD/StFlow.h
  • all derivative types also located in StFlow.h
  • derivative types:
    • AxiStagnFlow
    • FreeFlame
    • OneDFlow

• Bdry1D - base class for boundaries between 1D domains
  • located in /path/to/cantera/Cantera/src/oneD/Inlet1D.h
  • all derivative types also located in Inlet1D.h
  • derivative types:
    • Inlet1D
    • Symm1D
    • Outlet1D
    • OutletRes1D
    • Surf1D
    • ReactingSurf1D

• Solid1D - class for 1D reacting solids
  • located in /path/to/cantera/Cantera/src/oneD/Solid1D.h

• Surf1D - class for 0D surface
  • located in /path/to/cantera/Cantera/src/oneD/Surf1D.h

NOTE: given that StFlow is the most relevant type of 1D domain for a discussion of Jacobians, I will focus only on it.

Each 1D domain knows about neighboring domains (see the Domain1D::linkleft() and Domain1D::linkright() methods, as well as the Domain1D::locate() method).

Each 1D domain also knows what OneDim container contains it (stored in the m_container variable).

Jacobian Objects

MultiJac

• located in /path/to/cantera/Cantera/src/oneD/MultiJac.cpp
• Jacobian evaluator object

Jacobian Calculation

There is a MultiJac object contained in each 1D domain object (i.e. Domain1D base class), but there is also one associated with the collection of all 1D domains (i.e. OneDim base class). Only the Jacobian object in the OneDim class is actually utilized to solve the Jacobian (it appears to me that the Jacobian contained in all Domain1D classes are never used).

In the OneDim class, the Jacobian is stored in the variable m_jac.

OneDim::solve()

The OneDim::solve() method is the method from which the Jacobian is actually evaluated.

First, OneDim::eval() is called. This loops through all 1D domains (Domain1D objects) contained by the OneDim object and calls Domain1D::eval() (i.e. the eval() method for each individual 1D domain).

The Domain1D::eval() method, in turn, calculates the residual function at a given point (or at all points, if the first argument is less than 0, which it is when called by OneDim::eval()).

Once the OneDim::eval() method has been called, this information is passed to the MultiJac object's MultiJac::eval() method. (Too many eval() methods!)

MultiJac::eval()

This is the meat and potatoes of the Jacobian matrix solution. The method first loops over all points:

for (j = 0; j < m_points; j++) {

where MultiJac::m_points is equal to OneDim::m_size. Next, it loops over all species:

nv = m_resid->nVars(j);
for (n = 0; n < nv; n++) {

(where m_resid is a OneDim pointer).

Next there is some code to slightly perturb x0 at the given cell location:

// perturb x(n)
xsave = x0[ipt];
dx = m_atol + fabs(xsave)*m_rtol;
x0[ipt] = xsave + dx;                                  
dx = x0[ipt] - xsave;                                  
rdx = 1.0/dx;

And then it calculates the perturbed residual (it calls OneDim::eval(), which then calls Domain1D::eval(), but it is calling the variations of OneDim::eval() and Domain1D::eval() that evaluate the residual at one particular point in time.

After that, the Jacobian column corresponding to species n is computed:

                for (i = j - 1; i <= j+1; i++) {
                    if (i >= 0 && i < m_points) {                      
                        mv = m_resid->nVars(i);
                        iloc = m_resid->loc(i);
                        for (m = 0; m < mv; m++) {
                            value(m+iloc,ipt) = (m_r1[m+iloc]
                                - resid0[m+iloc])*rdx;
                        }
                    }
                }

this is assigned to the value() method, which saves it to an array named data[] (see /path/to/cantera/Cantera/src/numerics/BandMatrix.h, of which MultiJac is a derived type).

Printing the Jacobian

If you wish to print the Jacobian, you could print it at the end of the for loop over m_points:

for (j = 0; j < m_points; j++) {

    [...]

}

// print Jacobian here

This can be done by including the following two header files in MultiJac.cpp:

#include <fstream>
#include <iomanip>

and this block of code (which could be put inside a Jacobi::write() method) can be added to MultiJac.cpp:

char filename[28];
ofstream oStream;

// this will create a filename for the Jacobi matrix, indexed by m_age
// (it would be a better idea to create a separate counter, since m_age is not guaranteed to increment)
int sizeofit = sprintf( filename, "Jacobi_%.2d.mat", m_age );

oStream.open(filename);

for( int iRow = 0; iRow < nRows(); ++iRow ) {
    for( int iCol = 0; iCol < nColumns(); ++iCol ) {
        oStream << scientific << setw(8) << setprecision(8) << " " << value(iRow, iCol);
    }
    oStream << endl;
}
oStream.close();

Creating Symbolic Jacobian

It is possible in theory to create a symbolic Jacobian matrix.

For a given reaction mechanism, a CTML file is created. This is an XML file with metadata about the species and reactions involved in a mechanism.

Cantera provides a C++ interface for dealing with the XML tree in the CTML files. You could couple this interface with a C++ computer algebra system like Ginac to loop through all species, and create symbolic variables; you could then loop through all reactions, make functions for each k(T), and then construct the source terms and derivatives for each species with respect to each other species.

For the actual CTML/XML interface in Cantera, see:

/path/to/cantera/Cantera/src/base/ctml.h
/path/to/cantera/Cantera/src/base/xml.h

For a simple example of the CTML/XML interface in Cantera, see:

/path/to/cantera/test_problems/fracCoeff/fracCoeff.cc

and search for ctml::get_CTML_Tree.

An example CTML file is in

/path/to/cantera/test_problems/cxx_ex/gri30.xml

Pseudocode

Get base XML node
For each <phase>
  For each species in <speciesArray>
    Create new Ginac variable
  End
End

For each <reaction>
  For each species
    If reaction contains species i
      Construct source term expression for given reaction for given species
      Add source term expression to overall source term for given species
    End
  End
End