ResidEval.h

Go to the documentation of this file.

/**
 *  @file ResidEval.h
 */

// Copyright 2006  California Institute of Technology

#ifndef CT_RESIDEVAL_H
#define CT_RESIDEVAL_H

#include "cantera/base/ct_defs.h"
#include "cantera/base/ctexceptions.h"

#include <cstdio>

namespace Cantera
{

const int c_NONE = 0;
const int c_GE_ZERO = 1;
const int c_GT_ZERO = 2;
const int c_LE_ZERO = -1;
const int c_LT_ZERO = -2;

/**
 *  Virtual base class for DAE residual function evaluators.
 *  Classes derived from ResidEval evaluate the residual function
 * \f[
 *             \vec{F}(t,\vec{y}, \vec{y^\prime})
 * \f]
 * The DAE solver attempts to find a solution y(t) such that F = 0.
 *  @ingroup DAE_Group
 */
class ResidEval
{
public:
    ResidEval() {}
    virtual ~ResidEval() {}

    /**
     * Constrain solution component k. Possible values for
     * 'flag' are:
     *   - c_NONE       no constraint
     *   - c_GE_ZERO    >= 0
     *   - c_GT_ZERO    >  0
     *   - c_LE_ZERO    <= 0
     *   - c_LT_ZERO    <  0
     */
    virtual void constrain(const int k, const int flag) {
        m_constrain[k] = flag;
    }
    int constraint(const int k) const {
        std::map<int,int>::const_iterator i = m_constrain.find(k);
        if (i != m_constrain.end()) {
            return i->second;
        }
        return c_NONE;
    }

    //! Initialization function
    virtual void initSizes() {
        int neq = nEquations();
        m_alg.resize(neq, 0);
    }

    /**
     * Specify that solution component k is purely algebraic -
     * that is, the derivative of this component does not appear
     * in the residual function.
     */
    virtual void setAlgebraic(const int k) {
        if ((int) m_alg.size() < (k+1)) {
            initSizes();
        }
        m_alg[k] = 1;
    }

    virtual bool isAlgebraic(const int k) {
        return (m_alg[k] == 1);
    }


    /**
     * Evaluate the residual function. Called by the
     * integrator.
     * @param t time. (input)
     * @param y solution vector. (input)
     * @param ydot rate of change of solution vector. (input)
     * @param r residual vector (output)
     */
    virtual int eval(const doublereal t, const doublereal* const y,
                     const doublereal* const ydot,
                     doublereal* const r) {
        throw CanteraError("ResidEval::eval()", "base class called");
    }

    virtual int evalSS(const doublereal t, const doublereal* const y,
                       doublereal* const r) {
        return eval(t, y, 0, r);
    }

    virtual int evalSimpleTD(const doublereal t, const doublereal* const y,
                             const doublereal* const yold, doublereal deltaT,
                             doublereal* const r) {
        int nn = nEquations();
        vector_fp ydot(nn);
        for (int i = 0; i < nn; i++) {
            ydot[i] = (y[i] - yold[i]) / deltaT;
        }
        return eval(t, y, DATA_PTR(ydot), r);
    }

    /**
     * Fill the solution and derivative vectors with the initial
     * conditions at initial time t0.
     * @return 1   Everything is fine
     *         0 or neg Something went wrong
     */
    virtual int getInitialConditions(const doublereal t0, doublereal* const y,
                                     doublereal* const ydot) {
        initSizes();
        throw CanteraError("ResidEval::GetInitialConditions()", "base class called");
        return 1;
    }

    //! Return the number of equations in the equation system
    virtual int nEquations() const = 0;


    //!      Write out to a file or to standard output the current solution
    /*!
     *      ievent  is a description of the event that caused this
     *      function to be called.
     */
    virtual void writeSolution(int ievent, const double time,
                               const double deltaT,
                               const int time_step_num,
                               const double* y, const double* ydot) {
        int k;
        printf("ResidEval::writeSolution\n");
        printf("     Time = %g, ievent = %d, deltaT = %g\n", time, ievent, deltaT);
        if (ydot) {
            printf(" k    y[]  ydot[]\n");
            for (k = 0; k < nEquations(); k++) {
                printf("%d %g %g\n", k, y[k], ydot[k]);
            }
        } else {
            printf(" k    y[]\n");
            for (k = 0; k < nEquations(); k++) {
                printf("%d %g \n", k, y[k]);
            }
        }
    }

    //!  Return the number of parameters in the calculation
    /*!
     *  This is the number of parameters in the sensitivity calculation. We have
     *  set this to zero and have included it for later expansion
     */
    int nparams() const {
        return 0;
    }

protected:

    //! Mapping vector that stores whether a degree of freedom is a DAE or not
    /*!
     *   The first index is the equation number. The second index is 1 if it is a DAE,
     *   and zero if it is not.
     */
    std::vector<int> m_alg;
    std::map<int, int> m_constrain;

private:

};

}

#endif