Cantera/Surface Coverage

In Cantera, the surface coverage is equivalent to a mole fraction of a surface species.

Getting Coverages

Reactor::getInitialConditions()

Starting from the Reactor class, we see in the getInitialConditions method that we're looping over each Wall object associated with the Reactor, and we set the surface coverages for the wall.

    // set the first component to the total internal
    // energy
    y[0] = m_thermo->intEnergy_mass() * mass;

    // set the second component to the total volume
    y[1] = m_vol;

    // set the remaining components to the surface species
    // coverages on the walls
    size_t loc = m_nsp + 2;
    SurfPhase* surf;
    for (size_t m = 0; m < m_nwalls; m++) {
        surf = m_wall[m]->surface(m_lr[m]);
        if (surf) {
            m_wall[m]->getCoverages(m_lr[m], y + loc);
            loc += surf->nSpecies();
        }
    }

This piece of code is going to the Wall class and calling getCoverages. The arguments it passes are m_lr[m], which indicates whether the wall is left-sided or right-sided, and y+loc, which is a pointer to an array. The code y+loc means "the array y indexed at loc," which is set to the number of species plus two (the two being mass and volume).

The coverages, stored internally by the Wall object, are copied to the vector passed in the getCoverages function call, namely, y+loc.

Now, looking around more at the y array and the role that it plays in the updateState method, we can see that y is the solution vector - it holds all the solutions. Which means that when we call the getCoverages() method on the Wall object, we're passing it a pointer to our solution vector, and telling the Wall object to populate our solution vector with the coverages - that is, with the mole fractions of the surface species associated with the Wall object.

Wall::getCoverages()

The Wall::getCoverages method is not particularly interesting: it just copies data from one array into another.

void Wall::getCoverages(int leftright, doublereal* cov)
{
    if (leftright == 0) {
        copy(m_leftcov.begin(), m_leftcov.end(), cov);
    } else {
        copy(m_rightcov.begin(), m_rightcov.end(), cov);
    }
}

The m_leftcov and m_rightcov are C++ float vectors that hold the coverages of surface species on the left and right sides of the Wall.

SurfPhase

The Cantera Wall object basically wraps a Cantera Kinetics object and a Cantera SurfPhase object. It uses these to handle the underlying surface chemistry. The Wall class handles the interaction between the surface kinetics and the reactor. This happens via the Wall class's m_chem member (an array with two Kinetics objects, one for each side of the wall) and its m_surf member (an array with two SurfPhase objects, one for each side of the wall).

Setting Coverages

More interesting than getting the wall coverages is the setting the wall coverages.

Reactor::updateState()

The Reactor class's updateState method updates the Wall's surface coverages. In the file Reactor.cpp, the Reactor sets the coverages in the following code block:

    size_t loc = m_nsp + 2;
    SurfPhase* surf;
    for (size_t m = 0; m < m_nwalls; m++) {
        surf = m_wall[m]->surface(m_lr[m]);
        if (surf) {
            m_wall[m]->setCoverages(m_lr[m], y+loc);
            loc += surf->nSpecies();
        }
    }

Again, we're passing a flag to indicate whether the wall is left-handed or right-handed, and then we're passing our solution vector, indexed to the surface species (that's y+loc).

Wall::setCoverages()

In the Wall class, the setCoverages method is not doing anything too fancy. It's basically setting the contents of the Wall's coverage arrays (m_leftcov or m_rightcov) equal to the contents of another array (cov):

void Wall::setCoverages(int leftright, const doublereal* cov)
{
    if (leftright == 0) {
        copy(cov, cov + m_nsp[0], m_leftcov.begin());
    } else {
        copy(cov, cov + m_nsp[1], m_rightcov.begin());
    }
}

Calculating Coverages

More interesting than getting or setting coverages, is calculating coverages.

The money shot is in the Reactor::evalEqs method, which computes the wall terms for the Reactor:

    // compute wall terms
    size_t loc = m_nsp+2;
    fill(m_sdot.begin(), m_sdot.end(), 0.0);
    for (size_t i = 0; i < m_nwalls; i++) {
        int lr = 1 - 2*m_lr[i];
        double vdot = lr*m_wall[i]->vdot(time);
        m_vdot += vdot;
        m_Q += lr*m_wall[i]->Q(time);
        Kinetics* kin = m_wall[i]->kinetics(m_lr[i]);
        SurfPhase* surf = m_wall[i]->surface(m_lr[i]);
        if (surf && kin) {
            double rs0 = 1.0/surf->siteDensity();
            size_t nk = surf->nSpecies();
            double sum = 0.0;
            surf->setTemperature(m_state[0]);
            m_wall[i]->syncCoverages(m_lr[i]);
            kin->getNetProductionRates(DATA_PTR(m_work));
            size_t ns = kin->surfacePhaseIndex();
            size_t surfloc = kin->kineticsSpeciesIndex(0,ns);
            for (size_t k = 1; k < nk; k++) {
                ydot[loc + k] = m_work[surfloc+k]*rs0*surf->size(k);
                sum -= ydot[loc + k];
            }
            ydot[loc] = sum;
            loc += nk;

            double wallarea = m_wall[i]->area();
            for (size_t k = 0; k < m_nsp; k++) {
                m_sdot[k] += m_work[k]*wallarea;
            }
        }
    }

Now let's break down this block of code.

First, we check if we have a surface, and whether there are reactions on that surface - that's the kinetics() and surface() methods of Wall.

    // compute wall terms
    size_t loc = m_nsp+2;
    fill(m_sdot.begin(), m_sdot.end(), 0.0);
    for (size_t i = 0; i < m_nwalls; i++) {
        int lr = 1 - 2*m_lr[i];
        double vdot = lr*m_wall[i]->vdot(time);
        m_vdot += vdot;
        m_Q += lr*m_wall[i]->Q(time);
        Kinetics* kin = m_wall[i]->kinetics(m_lr[i]);
        SurfPhase* surf = m_wall[i]->surface(m_lr[i]);

Now, if we do have a surface, and we do have kinetics, we do two things: 1. Advance surface coverages 2. Compute the right-hand side of our species coverage equation, to advance the coverages

First, we compute some stuff that will be used in finding the RHS of our species coverage equation:

        if (surf && kin) {
            double rs0 = 1.0/surf->siteDensity();
            size_t nk = surf->nSpecies();
            double sum = 0.0;

Next, we set the surface temperature, sync our coverages, and get our net production rates from the surface kinetics:

            surf->setTemperature(m_state[0]);
            m_wall[i]->syncCoverages(m_lr[i]);
            kin->getNetProductionRates(DATA_PTR(m_work));

Finally, we compute the RHS for each of our surface species. This is a net reaction rate, and is found by multiplying production rate, inverse site density, and the number of surface sites occupied by species k (that's what surf->size(k) computes; see documentation for the SurfPhase::setCoverages() method):

            size_t ns = kin->surfacePhaseIndex();
            size_t surfloc = kin->kineticsSpeciesIndex(0,ns);
            for (size_t k = 1; k < nk; k++) {
                ydot[loc + k] = m_work[surfloc+k]*rs0*surf->size(k);
                sum -= ydot[loc + k];
            }
            ydot[loc] = sum;
            loc += nk;

This allows us to compute RHS for surface coverage ODEs.