Jacobian in Cantera

The Jacobian matrix is constructed in Cantera by the one-dimensional domain. The hierarchy is described in the file http://files.charlesmartinreid.com/Cantera_Flames.pdf

Each domain creates a 1D domain.

Background Info

Jacobian Matrix

(insert stuff here)

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:

#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
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(20) << setprecision(20) << " " << value(iRow, iCol);
    }
    oStream << endl;
}
oStream.close();