ImplicitSurfChem.h

Go to the documentation of this file.

/**
 *  @file ImplicitSurfChem.h
 * Declarations for the implicit integration of surface site density equations
 *  (see \ref  kineticsmgr and class
 *  \link Cantera::ImplicitSurfChem ImplicitSurfChem\endlink).
 */
// Copyright 2001  California Institute of Technology

#ifndef CT_IMPSURFCHEM_H
#define CT_IMPSURFCHEM_H

#include "cantera/numerics/FuncEval.h"
#include "cantera/numerics/Integrator.h"
#include "cantera/kinetics/InterfaceKinetics.h"
#include "cantera/thermo/SurfPhase.h"
#include "solveSP.h"

namespace Cantera
{

class solveSP;


//! Advances the surface coverages of the associated set of SurfacePhase
//! objects in time
/*!
 *  This function advances a set of SurfacePhase objects, each
 *  associated with one InterfaceKinetics object, in time.
 *  The following equation is used for each surface phase, <I>i</I>.
 *
 *   \f[
 *        \dot \theta_k = \dot s_k (\sigma_k / s_0)
 *   \f]
 *
 *  In this equation,
 *    \f$ \theta_k \f$ is the site coverage for the kth species.
 *    \f$ \dot s_k \f$ is the source term for the kth species
 *   \f$ \sigma_k \f$ is the number of surface sites covered by
 *  each species k.
 *   \f$ s_0 \f$ is the total site density of the interfacial phase.
 *
 *  Additionally, the 0'th equation in the set is discarded. Instead the
 *  alternate equation is solved for
 *
 *   \f[
 *        \sum_{k=0}^{N-1}  \dot \theta_k = 0
 *   \f]
 *
 *  This last equation serves to ensure that sum of the \f$ \theta_k \f$
 *  values stays constant.
 *
 *  The object uses the CVODE software to advance the surface equations.
 *
 *  The solution vector used by this object is as follows.
 *   For each surface phase with \f$ N_s \f$ surface sites,
 *   it consists of the surface coverages
 *       \f$ \theta_k \f$ for \f$ k = 0, N_s - 1 \f$
 *
 * @ingroup  kineticsmgr
 *
 */
class ImplicitSurfChem : public FuncEval
{

public:


    //! Constructor for multiple surfaces.
    /*!
     * @param k  Vector of pointers to InterfaceKinetics objects
     *           Each object consists of a surface or an edge containing
     *           internal degrees of freedom representing the concentration
     *           of surface adsorbates.
     */
    ImplicitSurfChem(std::vector<InterfaceKinetics*> k);

    /**
     * Destructor. Deletes the integrator.
     */
    virtual ~ImplicitSurfChem();

    /**
     * Overloads the virtual function
     * declared in FuncEval.
     */
    virtual void initialize(doublereal t0 = 0.0);


    //! Integrate from t0 to t1. The integrator is reinitialized first.
    /*!
     *   This routine does a time accurate solve from t = t0 to t = t1.
     *   of the surface problem.
     *
     *  @param t0  Initial Time -> this is an input
     *  @param t1  Final Time -> This is an input
     */
    void integrate(doublereal t0, doublereal t1);


    //! Integrate from t0 to t1 without reinitializing the integrator.
    /*!
     *  Use when the coverages have not changed from
     *  their values on return from the last call to integrate or
     *  integrate0.
     *
     *  @param t0  Initial Time -> this is an input
     *  @param t1  Final Time -> This is an input
     */
    void integrate0(doublereal t0, doublereal t1);


    //! Solve for the pseudo steady-state of the surface problem
    /*!
     * Solve for the steady state of the surface problem.
     * This is the same thing as the advanceCoverages() function,
     * but at infinite times.
     *
     * Note, a direct solve is carried out under the hood here,
     * to reduce the computational time.
     *
     * @param ifuncOverride 4 values are possible
     *                    1  SFLUX_INITIALIZE
     *                    2  SFLUX_RESIDUAL
     *                    3  SFLUX_JACOBIAN
     *                    4  SFLUX_TRANSIENT
     *         The default is -1, which means that the program
     *         will decide.
     * @param timeScaleOverride When a pseudo transient is
     *             selected this value can be used to override
     *             the default time scale for integration which
     *             is one.
     *             When SFLUX_TRANSIENT is used, this is equal to the
     *             time over which the equations are integrated.
     *             When SFLUX_INITIALIZE is used, this is equal to the
     *             time used in the initial transient algorithm,
     *             before the equation system is solved directly.
     */
    void solvePseudoSteadyStateProblem(int ifuncOverride = -1,
                                       doublereal timeScaleOverride = 1.0);


    // overloaded methods of class FuncEval

    //! Return the number of equations
    virtual size_t neq() {
        return m_nv;
    }

    //! Evaluate the value of ydot[k] at the current conditions
    /*!
     *  @param t   Time (seconds)
     *  @param y   Vector containing the current solution vector
     *  @param ydot   Output vector containing the value of the
     *                derivative of the surface coverages.
     *  @param p   Unused parameter pass-through parameter vector
     */
    virtual void eval(doublereal t, doublereal* y, doublereal* ydot,
                      doublereal* p);

    //! Set the initial conditions for the solution vector
    /*!
     *  @param t0  Initial time
     *  @param leny  Length of the solution vector
     *  @param y   Value of the solution vector to be used.
     *            On output, this contains the initial value
     *           of the solution.
     */
    virtual void getInitialConditions(doublereal t0,
                                      size_t leny, doublereal* y);

    /*
     * Get the specifications for the problem from the values
     * in the ThermoPhase objects for all phases.
     *
     *  1) concentrations of all species in all phases, m_concSpecies[]
     *  2) Temperature and pressure
     *
     *
     *  @param vecConcSpecies Vector of concentrations. The
     *                  phase concentration vectors are contiguous
     *                  within the object, in the same order as the
     *                  unknown vector.
     */
    void getConcSpecies(doublereal* const vecConcSpecies) const;

    //! Sets the concentrations within phases that are unknowns in
    //! the surface problem
    /*!
     * Fills the local concentration vector for all of the
     * species in all of the phases that are unknowns in the surface
     * problem.
     *
     *  @param vecConcSpecies Vector of concentrations. The
     *                  phase concentration vectors are contiguous
     *                  within the object, in the same order as the
     *                  unknown vector.
     */
    void setConcSpecies(const doublereal* const vecConcSpecies);

    //! Sets the state variable in all thermodynamic phases (surface and
    //! surrounding bulk phases) to the input temperature and pressure
    /*!
     *  @param TKelvin input temperature (kelvin)
     *  @param PresPa   input pressure in pascal.
     */
    void setCommonState_TP(doublereal TKelvin, doublereal PresPa);


    //! Returns a reference to the vector of pointers to the
    //! InterfaceKinetics objects
    /*!
     * This should probably go away in the future, as it opens up the
     * class.
     */
    std::vector<InterfaceKinetics*> & getObjects() {
        return m_vecKinPtrs;
    }

    int checkMatch(std::vector<ThermoPhase*> m_vec, ThermoPhase* thPtr);

    void setIOFlag(int ioFlag) {
        m_ioFlag = ioFlag;
    }

protected:


    //! Set the mixture to a state consistent with solution
    //! vector y.
    /*!
     *  This function will set the surface site factions
     *  in the underlying %SurfPhase objects to the current
     *  value of the solution vector.
     *
     * @param y Current value of the solution vector.
     *          The lenth is equal to the sum of the number of surface
     *          sites in all the surface phases
     */
    void updateState(doublereal* y);

    //! vector of pointers to surface phases.
    std::vector<SurfPhase*> m_surf;

    //! Vector of pointers to bulk phases
    std::vector<ThermoPhase*> m_bulkPhases;

    //! vector of pointers to InterfaceKinetics objects
    std::vector<InterfaceKinetics*> m_vecKinPtrs;

    //! Vector of number of species in each Surface Phase
    std::vector<size_t>  m_nsp;

    //! index of the surface phase in each InterfaceKinetics object
    std::vector<size_t> m_surfindex;

    std::vector<size_t> m_specStartIndex;

    //! Total number of surface phases.
    /*!
     * This is also equal to the number of InterfaceKinetics objects
     * as there is a 1-1 correspondence between InterfaceKinetics objects
     * and surface phases.
     */
    size_t m_nsurf;

    //! Total number of surface species in all surface phases
    /*!
     * This is the total number of unknowns in m_mode 0 problem
     */
    size_t m_nv;

    size_t m_numBulkPhases;
    std::vector<size_t> m_nspBulkPhases;
    size_t  m_numTotalBulkSpecies;
    size_t m_numTotalSpecies;

    std::vector<vector_int> pLocVec;
    //! Pointer to the cvode integrator
    Integrator* m_integ;
    doublereal m_atol, m_rtol;   // tolerances
    doublereal m_maxstep;        // max step size
    vector_fp m_work;



    /**
     * Temporary vector - length num species in the Kinetics object.
     *                    This is the sum of the number of species
     *                    in each phase included in the kinetics object.
     */
    vector_fp m_concSpecies;
    vector_fp m_concSpeciesSave;

    //std::vector<vector_fp> m_vectorConcKinSpecies;
    //std::vector<vector_fp> m_vectorNetSpeciesProdRate;
    /**
     * Index into the species vector of the kinetics manager,
     * pointing to the first species from the surrounding medium.
     */
    int m_mediumSpeciesStart;
    /**
     * Index into the species vector of the kinetics manager,
     * pointing to the first species from the condensed phase
     * of the particles.
     */
    int m_bulkSpeciesStart;
    /**
     * Index into the species vector of the kinetics manager,
     * pointing to the first species from the surface
     * of the particles
     */
    int m_surfSpeciesStart;
    /**
     * Pointer to the helper method, Placid, which solves the
     * surface problem.
     */
    solveSP* m_surfSolver;

    //! If true, a common temperature and pressure for all
    //! surface and bulk phases associated with the surface problem
    //! is imposed
    bool m_commonTempPressForPhases;

    //! We make the solveSS class a friend because we need
    //! to access all of the above information directly.
    //! Adding the members into the class is also a possibility.
    friend class solveSS;

private:

    //! Controls the amount of printing from this routine
    //! and underlying routines.
    int m_ioFlag;
};
}

#endif