FuncEval.h

Go to the documentation of this file.

/**
 *  @file FuncEval.h
 */
 
// This file is part of Cantera. See License.txt in the top-level directory or
// at https://cantera.org/license.txt for license and copyright information.
 
#ifndef CT_FUNCEVAL_H
#define CT_FUNCEVAL_H
 
#include "cantera/base/ct_defs.h"
#include "cantera/base/ctexceptions.h"
#include "cantera/base/global.h"
 
namespace Cantera
{
 
/**
 * @defgroup odeGroup  ODE and DAE Function Utilities
*/
 
/**
 *  Virtual base class for ODE/DAE right-hand-side function evaluators.
 *  Classes derived from FuncEval evaluate the right-hand-side function
 * @f$ \vec{F}(t,\vec{y}) @f$ in
 * @f[
 *  \dot{\vec{y}} = \vec{F}(t,\vec{y}).
 * @f]
 *  @ingroup odeGroup
 */
class FuncEval
{
public:
    FuncEval() = default;
    virtual ~FuncEval() = default;
 
    /**
     * Evaluate the right-hand-side ODE function. Called by the integrator.
     * @param[in] t time.
     * @param[in] y solution vector, length neq()
     * @param[out] ydot rate of change of solution vector, length neq()
     * @param[in] p sensitivity parameter vector, length nparams()
     */
    virtual void eval(double t, span<const double> y, span<double> ydot,
                      span<const double> p) {
        throw NotImplementedError("FuncEval::eval");
    }
 
    /**
     * Evaluate the right-hand-side DAE function. Called by the integrator.
     * @param[in] t time.
     * @param[in] y solution vector, length neq()
     * @param[in] ydot rate of change of solution vector, length neq()
     * @param[in] p sensitivity parameter vector, length nparams()
     * @param[out] residual the DAE residuals, length nparams()
     */
    virtual void evalDae(double t, span<const double> y, span<const double> ydot,
                         span<const double> p, span<double> residual) {
        throw NotImplementedError("FuncEval::evalDae");
    }
 
    //! Given a vector of length neq(), mark which variables should be
    //! considered algebraic constraints
    virtual void getConstraints(span<double> constraints) {
        throw NotImplementedError("FuncEval::getConstraints");
    }
 
    //! Evaluate the right-hand side using return code to indicate status.
    /*!
     * Errors are indicated using the return value, rather than by throwing
     *  exceptions. This method is used when calling from a C-based integrator
     *  such as CVODES. Exceptions may either be stored or printed, based on the
     *  setting of suppressErrors().
     *  @returns 0 for a successful evaluation; 1 after a potentially-
     *      recoverable error; -1 after an unrecoverable error.
     */
    int evalNoThrow(double t, span<const double> y, span<double> ydot);
 
    //! Evaluate the right-hand side using return code to indicate status.
    /*!
     * Errors are indicated using the return value, rather than by throwing
     *  exceptions. This method is used when calling from a C-based integrator
     *  such as CVODES. Exceptions may either be stored or printed, based on the
     *  setting of suppressErrors().
     *  @returns 0 for a successful evaluation; 1 after a potentially-
     *      recoverable error; -1 after an unrecoverable error.
     */
    int evalDaeNoThrow(double t, span<const double> y, span<const double> ydot,
                       span<double> residual);
 
    /**
     * Evaluate the setup processes for the Jacobian preconditioner.
     * @param[in] t time.
     * @param[in] y solution vector, length neq()
     * @param gamma the gamma in M=I-gamma*J
     * @warning This function is an experimental part of the %Cantera API and may be
     * changed or removed without notice.
     */
    virtual void preconditionerSetup(double t, span<const double> y, double gamma) {
        throw NotImplementedError("FuncEval::preconditionerSetup");
    }
 
    /**
     * Evaluate the linear system Ax=b where A is the preconditioner.
     * @param[in] rhs right hand side vector used in linear system
     * @param[out] output output vector for solution
     * @warning This function is an experimental part of the %Cantera API and may be
     * changed or removed without notice.
     */
    virtual void preconditionerSolve(span<const double> rhs, span<double> output) {
        throw NotImplementedError("FuncEval::preconditionerSolve");
    }
 
    //! Update the preconditioner based on already computed jacobian values
    virtual void updatePreconditioner(double gamma) {
        throw NotImplementedError("FuncEval::updatePreconditioner");
    }
 
    /**
     * Preconditioner setup that doesn't throw an error but returns a
     * CVODES flag. It also helps as a first level of polymorphism
     * which identifies the specific FuncEval, e.g., ReactorNet.
     * @param[in] t time.
     * @param[in] y solution vector, length neq()
     * @param gamma the gamma in M=I-gamma*J
     * @warning This function is an experimental part of the %Cantera API and may be
     * changed or removed without notice.
     */
    int preconditioner_setup_nothrow(double t, span<const double> y, double gamma);
 
    /**
     * Preconditioner solve that doesn't throw an error but returns a
     * CVODES flag. It also helps as a first level of polymorphism
     * which identifies the specific FuncEval, e.g., ReactorNet.
     * @param[in] rhs right hand side vector used in linear system
     * @param[out] output output vector for solution
     * @warning This function is an experimental part of the %Cantera API and may be
     * changed or removed without notice.
     */
    int preconditioner_solve_nothrow(span<const double> rhs, span<double> output);
 
    //! Number of event/root functions exposed to the integrator.
    //! 0 indicates root finding is disabled.
    virtual size_t nRootFunctions() const {
        return 0;
    }
 
    /**
     * Evaluate the event/root functions currently in play.
     * Integrators invoke this whenever root finding is enabled; implementations
     * should fill `gout` with the function values.
     * @param[in] t Time at which to evaluate the root functions
     * @param[in] y Current solution vector at time *t* of length neq()
     * @param[out] gout Array of length nRootFunctions() to be filled with the
     *      values of the root functions
     */
    virtual void evalRootFunctions(double t, span<const double> y, span<double> gout) {}
 
    //! Wrapper for evalRootFunctions that converts exceptions to return codes.
    //! @returns 0 for a successful evaluation, 1 after a potentially-
    //!     recoverable error, or -1 after an unrecoverable error.
    int evalRootFunctionsNoThrow(double t, span<const double> y, span<double> gout);
 
    //! Fill in the vector *y* with the current state of the system.
    //! Used for getting the initial state for ODE systems.
    virtual void getState(span<double> y) {
        throw NotImplementedError("FuncEval::getState");
    }
 
    //! Fill in the vectors *y* and *ydot* with the current state of the system.
    //! Used for getting the initial state for DAE systems.
    virtual void getStateDae(span<double> y, span<double> ydot) {
        throw NotImplementedError("FuncEval::getStateDae");
    }
 
    //! Number of equations.
    virtual size_t neq() const = 0;
 
    //! Number of sensitivity parameters.
    virtual size_t nparams() const {
        return m_sens_params.size();
    }
 
    //! Enable or disable suppression of errors when calling eval()
    void suppressErrors(bool suppress) {
        m_suppress_errors = suppress;
    }
 
    //! Get current state of error suppression
    bool suppressErrors() const {
        return m_suppress_errors;
    };
 
    //! Return a string containing the text of any suppressed errors
    string getErrors() const;
 
    //! Clear any previously-stored suppressed errors
    void clearErrors() {
        m_errors.clear();
    };
 
    //! Values for the problem parameters for which sensitivities are computed
    //! This is the array which is perturbed and passed back as the fourth
    //! argument to eval().
    vector<double> m_sens_params;
 
    //! Scaling factors for each sensitivity parameter
    vector<double> m_paramScales;
 
protected:
    // If true, errors are accumulated in m_errors. Otherwise, they are printed
    bool m_suppress_errors = false;
 
    //! Errors occurring during function evaluations
    vector<string> m_errors;
};
 
}
 
#endif