Domain1D.h

Go to the documentation of this file.

 //! @file Domain1D.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_DOMAIN1D_H
#define CT_DOMAIN1D_H
 
#include "cantera/base/ctexceptions.h"
#include "cantera/base/global.h"
 
namespace Cantera
{
 
class MultiJac;
class OneDim;
class Refiner;
class AnyMap;
class Kinetics;
class Transport;
class Solution;
class SolutionArray;
class SystemJacobian;
 
/**
 * Base class for one-dimensional domains.
 * @ingroup flowGroup
 */
class Domain1D
{
protected:
    /**
     * Constructor.
     * @param nv      Number of variables at each grid point.
     * @param points  Number of grid points.
     * @param time    (unused)
     */
    explicit Domain1D(size_t nv=1, size_t points=1, double time=0.0);
 
public:
    virtual ~Domain1D();
    Domain1D(const Domain1D&) = delete;
    Domain1D& operator=(const Domain1D&) = delete;
 
    //! Domain type flag.
    //! @since Starting in %Cantera 3.1, the return type is a `string`.
    virtual string domainType() const { return "domain"; }
 
    //! The left-to-right location of this domain.
    size_t domainIndex() {
        return m_index;
    }
 
    //! True if the domain is a connector domain.
    virtual bool isConnector() {
        return false;
    }
 
    //! Set transport model by name.
    //! @param model  String specifying model name.
    //! @since New in %Cantera 3.2.
    virtual void setTransportModel(const string& model) {
        throw NotImplementedError("Domain1D::setTransportModel");
    }
 
protected:
    //! Update transport model to existing instance
    //! @since New in %Cantera 3.2.
    virtual void _setKinetics(shared_ptr<Kinetics> kin) {
        throw NotImplementedError("Domain1D::_setKinetics");
    }
 
    //! Update transport model to existing instance
    //! @since New in %Cantera 3.2.
    virtual void _setTransport(shared_ptr<Transport> trans) {
        throw NotImplementedError("Domain1D::_setTransport");
    }
 
public:
    //! The container holding this domain.
    const OneDim& container() const {
        return *m_container;
    }
 
    //! Specify the container object for this domain, and the position of this
    //! domain in the list.
    void setContainer(OneDim* c, size_t index) {
        m_container = c;
        m_index = index;
    }
 
    //! Set the Jacobian bandwidth. See the discussion of method bandwidth().
    void setBandwidth(int bw = -1) {
        m_bw = bw;
    }
 
    //! Set the Jacobian bandwidth for this domain.
    /**
     * When class OneDim computes the bandwidth of the overall multi-domain
     * problem (in OneDim::resize()), it calls this method for the bandwidth
     * of each domain. If setBandwidth has not been called, then a negative
     * bandwidth is returned, in which case OneDim assumes that this domain is
     * dense -- that is, at each point, all components depend on the value of
     * all other components at that point. In this case, the bandwidth is bw =
     * 2*nComponents() - 1. However, if this domain contains some components
     * that are uncoupled from other components at the same point, then this
     * default bandwidth may greatly overestimate the true bandwidth, with a
     * substantial penalty in performance. For such domains, use method
     * setBandwidth to specify the bandwidth before passing this domain to the
     * Sim1D or OneDim constructor.
     */
    size_t bandwidth() {
        return m_bw;
    }
 
    /**
     * Initialize. This method is called by OneDim::init() for each domain once
     * at the beginning of a simulation. Base class method does nothing, but may
     * be overloaded.
     */
    virtual void init() {  }
 
    /**
     * When called, this function should reset "bad" values in the state vector
     * such as negative species concentrations. This function may be called
     * after a failed solution attempt.
     */
    virtual void resetBadValues(span<double> x) {}
 
    /**
     * Resize the domain to have nv components and np grid points. This method
     * is virtual so that subclasses can perform other actions required to
     * resize the domain.
     */
    virtual void resize(size_t nv, size_t np);
 
    //! Return a reference to the grid refiner.
    Refiner& refiner() {
        return *m_refiner;
    }
 
    //! Number of components at each grid point.
    size_t nComponents() const {
        return m_nv;
    }
 
    //! Number of grid points in this domain.
    size_t nPoints() const {
        return m_points;
    }
 
    //! Name of component `n`. May be overloaded.
    virtual string componentName(size_t n) const;
 
    //! Set the name of the component `n` to `name`.
    void setComponentName(size_t n, const string& name) {
        m_name[n] = name;
    }
 
    /**
     * Index of component with name `name`.
     * @param name  name of component
     * @param checkAlias  if `true` (default), check alias mapping
     */
    virtual size_t componentIndex(const string& name, bool checkAlias=true) const;
 
    /**
     * Check whether the Domain contains a component.
     * @param name  name of component
     * @param checkAlias  if `true` (default), check alias mapping
     *
     * @since New in %Cantera 3.2.
     */
    virtual bool hasComponent(const string& name, bool checkAlias=true) const;
 
    //! Index of component `name` at grid point `j` within the global solution
    //! vector of the containing OneDim/Sim1D object.
    //! @since New in %Cantera 4.0.
    size_t globalComponentIndex(const string& name, size_t j) const {
        return loc() + index(componentIndex(name), j);
    }
 
    /**
     * Update state at given location to state of associated Solution object.
     */
    virtual void updateState(size_t loc) {
        throw NotImplementedError("Domain1D::updateState",
            "Not implemented for domain type '{}'.", domainType());
    }
 
    /**
     * Set the upper and lower bounds for a solution component, n.
     *
     * @param n  solution component index
     * @param lower  lower bound on component n
     * @param upper  upper bound on component n
     */
    void setBounds(size_t n, double lower, double upper) {
        m_min[n] = lower;
        m_max[n] = upper;
    }
 
    //! Set tolerances for time-stepping mode
    /*!
     * @param rtol  Relative tolerance
     * @param atol  Absolute tolerance
     * @param n  component index these tolerances apply to. If set to -1 (the
     *      default), these tolerances will be applied to all solution
     *      components.
     */
    void setTransientTolerances(double rtol, double atol, size_t n=npos);
 
    //! Set tolerances for steady-state mode
    /*!
     * @param rtol  Relative tolerance
     * @param atol  Absolute tolerance
     * @param n  component index these tolerances apply to. If set to -1 (the
     *     default), these tolerances will be applied to all solution
     *     components.
     */
    void setSteadyTolerances(double rtol, double atol, size_t n=npos);
 
    //! Relative tolerance of the nth component.
    double rtol(size_t n) {
        return (m_rdt == 0.0 ? m_rtol_ss[n] : m_rtol_ts[n]);
    }
 
    //! Absolute tolerance of the nth component.
    double atol(size_t n) {
        return (m_rdt == 0.0 ? m_atol_ss[n] : m_atol_ts[n]);
    }
 
    //! Steady relative tolerance of the nth component
    double steady_rtol(size_t n) {
        return m_rtol_ss[n];
    }
 
    //! Steady absolute tolerance of the nth component
    double steady_atol(size_t n) {
        return m_atol_ss[n];
    }
 
    //! Transient relative tolerance of the nth component
    double transient_rtol(size_t n) {
        return m_rtol_ts[n];
    }
 
    //! Transient absolute tolerance of the nth component
    double transient_atol(size_t n) {
        return m_atol_ts[n];
    }
 
    //! Upper bound on the nth component.
    double upperBound(size_t n) const {
        return m_max[n];
    }
 
    //! Lower bound on the nth component
    double lowerBound(size_t n) const {
        return m_min[n];
    }
 
    //! Set the method used to evaluate this domain's Jacobian columns.
    //! @param mode  One of `"auto"` (default), `"analytic"`, or
    //!     `"finite-difference"`. In `"auto"` mode the domain uses analytic
    //!     Jacobian columns where supported (see hasAnalyticJacobian()) and
    //!     silently falls back to finite differences otherwise. `"analytic"`
    //!     additionally raises a CanteraError if analytic evaluation was
    //!     requested but cannot be used because the kinetics object lacks the
    //!     required composition derivatives or multicomponent transport is
    //!     active. `"finite-difference"` always uses finite differences.
    //! @since New in %Cantera 4.0.
    void setJacobianMode(const string& mode) {
        if (mode != "finite-difference" && mode != "analytic" && mode != "auto") {
            throw CanteraError("Domain1D::setJacobianMode",
                "Unknown Jacobian mode '{}'", mode);
        }
        m_jacobianMode = mode;
    }
 
    //! Get the method used to evaluate this domain's Jacobian columns.
    //! @since New in %Cantera 4.0.
    const string& jacobianMode() const {
        return m_jacobianMode;
    }
 
    //! Returns `true` if this domain computes the Jacobian column for
    //! component `n` at (domain-local) grid point `j` analytically, in which
    //! case the finite-difference evaluation of that column is skipped and the
    //! domain must provide all entries of the column in evalJacobianAnalytic().
    //! @since New in %Cantera 4.0.
    virtual bool hasAnalyticJacobian(size_t j, size_t n) const {
        return false;
    }
 
    //! Add this domain's analytic Jacobian entries (for columns claimed by
    //! hasAnalyticJacobian()) to `jac` via SystemJacobian::setValue(), using
    //! global row/column indices (domain-local index + loc()). Entries must be
    //! the steady-state Jacobian; transient diagonal terms are handled
    //! separately. Base class implementation does nothing.
    //! @since New in %Cantera 4.0.
    virtual void evalJacobianAnalytic(span<const double> x, SystemJacobian& jac) {}
 
    //! Validate that an explicitly requested analytic Jacobian
    //! (`jacobian_mode == "analytic"`) can be used for the current
    //! configuration, throwing CanteraError otherwise. Called once per
    //! Jacobian evaluation before the finite-difference column loop. The base
    //! implementation does nothing; the `"auto"` and `"finite-difference"`
    //! modes never raise.
    //! @since New in %Cantera 4.0.
    virtual void checkAnalyticJacobian() const {}
 
    /**
     * Set grid refinement criteria. @see Refiner::setCriteria.
     * @since New in %Cantera 3.2
     */
    void setRefineCriteria(double ratio = 10.0,
                           double slope = 0.8, double curve = 0.8,
                           double prune = -0.1);
 
    /**
     * Get the grid refinement criteria. @see Refiner::getCriteria
     * @since New in %Cantera 3.2
     */
    vector<double> getRefineCriteria();
 
    /**
     * Performs the setup required before starting a time-stepping solution.
     * Stores the solution provided in `x0` to the internal storage, and sets
     * the reciprocal of the time step to `1/dt`.
     *
     * @param[in] dt  Time step
     * @param[in] x0  Array to store the solution at the last time step
     */
    void initTimeInteg(double dt, span<const double> x0) {
        auto local = x0.subspan(loc(), size());
        std::copy(local.begin(), local.end(), m_slast.begin());
        m_rdt = 1.0/dt;
    }
 
    /**
     * Set the internally-stored reciprocal of the time step to 0.0, which is
     * used to indicate that the problem is in steady-state mode.
     */
    void setSteadyMode() {
        m_rdt = 0.0;
    }
 
    //! True if in steady-state mode
    bool steady() {
        return (m_rdt == 0.0);
    }
 
    //! True if not in steady-state mode
    bool transient() {
        return (m_rdt != 0.0);
    }
 
    /**
     * Set this if something has changed in the governing
     * equations (for example, the value of a constant has been changed,
     * so that the last-computed Jacobian is no longer valid.
     */
    void needJacUpdate();
 
    //! Evaluate the residual function at point j. If j == npos,
    //! evaluate the residual function at all points.
    /*!
     *  This function must be implemented in classes derived from Domain1D.
     *
     *  @param[in] j  Grid point at which to update the residual
     *  @param[in] x  State vector
     *  @param[out] r  residual vector
     *  @param[out] mask  Boolean mask indicating whether each solution
     *      component has a time derivative (1) or not (0).
     *  @param[in] rdt  Reciprocal of the timestep (`rdt=0` implies steady-state.)
     */
    virtual void eval(size_t j, span<const double> x, span<double> r, span<int> mask,
                      double rdt=0.0) {
        throw NotImplementedError("Domain1D::eval");
    }
 
    /**
     * Returns the index of the solution vector, which corresponds to component
     * n at grid point j.
     *
     * @param n  component index
     * @param j  grid point index
     */
    size_t index(size_t n, size_t j) const {
        return m_nv*j + n;
    }
 
    /**
     * Set a single component value at a boundary.
     * @param component  Name of the component.
     *
     * @since New in %Cantera 3.2.
     */
    virtual double value(const string& component) const {
        throw NotImplementedError("Domain1D::value",
            "Not implemented for domain type '{}'.", domainType());
    }
 
    /**
     * Set a single component value in a flow domain or at a boundary.
     * @param component  Name of the component.
     * @param value  Value of the component.
     *
     * @since New in %Cantera 3.2.
     */
    virtual void setValue(const string& component, double value) {
        throw NotImplementedError("Domain1D::setValue",
            "Not implemented for domain type '{}'.", domainType());
    }
 
    /**
     * Retrieve component values.
     * @param component  Name of the component.
     * @returns  Vector of length nPoints() containing values at grid points.
     *
     * @since New in %Cantera 3.2.
     */
    vector<double> values(const string& component) const {
        vector<double> data(nPoints());
        getValues(component, data);
        return data;
    }
 
    /**
     * Retrieve component values.
     * @param component  Name of the component.
     * @param[out] values  Vector of length nPoints() containing values at grid points.
     *
     * @since New in %Cantera 3.2.
     */
    virtual void getValues(const string& component, span<double> values) const {
        throw NotImplementedError("Domain1D::getValues",
            "Not implemented for domain type '{}'.", domainType());
    }
 
    /**
     * Specify component values.
     * @param component  Name of the component.
     * @param[in] values  Vector of length nPoints() containing values at grid points.
     *
     * @since New in %Cantera 3.2.
     */
    virtual void setValues(const string& component, span<const double> values) {
        throw NotImplementedError("Domain1D::setValues",
            "Not implemented for domain type '{}'.", domainType());
    }
 
    /**
     * Retrieve internal work array values for a component.
     * After calling Sim1D::eval(), this array contains the values of the residual
     * function.
     * @param component  Name of the component.
     * @returns  Vector of length nPoints() containing residuals at grid  points.
     *
     * @since New in %Cantera 3.2.
     */
    vector<double> residuals(const string& component) const {
        vector<double> data(nPoints());
        getResiduals(component, data);
        return data;
    }
 
    /**
     * Retrieve internal work array values for a component.
     * After calling Sim1D::eval(), this array contains the values of the residual
     * function.
     * @param component  Name of the component.
     * @param[out] values  Vector of length nPoints() containing residuals at grid
     *      points.
     *
     * @since New in %Cantera 3.2.
     */
    virtual void getResiduals(const string& component, span<double> values) const {
        throw NotImplementedError("Domain1D::getResiduals",
            "Not applicable or not implemented for domain type '{}'.", domainType());
    }
 
    /**
     * Specify a profile for a component.
     * @param component  Name of the component.
     * @param[in] pos  A vector of relative positions, beginning with 0.0 at the
     *     left of the domain, and ending with 1.0 at the right of the domain.
     * @param[in] values  A vector of values corresponding to the relative position
     *     locations.
     *
     * Note that the vector pos and values can have lengths different than the
     * number of grid points, but their lengths must be equal. The values at
     * the grid points will be linearly interpolated based on the (pos,
     * values) specification.
     *
     * @since New in %Cantera 3.2.
     */
    virtual void setProfile(const string& component,
                            span<const double> pos, span<const double> values) {
        throw NotImplementedError("Domain1D::setProfile",
            "Not implemented for domain type '{}'.", domainType());
    }
 
    /**
     * Specify a flat profile for a component.
     * @param component  Name of the component.
     * @param value  Constant value.
     *
     * @since New in %Cantera 3.2.
     */
    virtual void setFlatProfile(const string& component, double value) {
        throw NotImplementedError("Domain1D::setFlatProfile",
            "Not implemented for domain type '{}'.", domainType());
    }
 
    //! Save the state of this domain to a SolutionArray.
    /*!
     * This method serves as an external interface for high-level API's; it does not
     * provide direct access to memory.
     * @param normalize  If true, normalize concentrations (default=false)
     *
     * @since New in %Cantera 3.0.
     */
    virtual shared_ptr<SolutionArray> toArray(bool normalize=false) {
        throw NotImplementedError("Domain1D::toArray", "Needs to be overloaded.");
    }
 
    //! Restore the solution for this domain from a SolutionArray.
    /*!
     * This method serves as an external interface for high-level API's.
     * @param arr  SolutionArray defining the state of this domain
     * @since New in %Cantera 3.0.
     */
    virtual void fromArray(const shared_ptr<SolutionArray>& arr) {
        throw NotImplementedError("Domain1D::fromArray", "Needs to be overloaded.");
    }
 
    /**
     * Return a concise summary of a Domain.
     * @see SolutionArray.info()
     * @param keys  List of components to be displayed; if empty, all components are
     *      considered.
     * @param rows  Maximum number of rendered rows.
     * @param width  Maximum width of rendered output.
     * @since New in %Cantera 3.2
     */
    string info(const vector<string>& keys, int rows=10, int width=80);
 
    /**
     * Return a concise summary of a Domain.
     * Skips keys input while `vector<string>` is not implemented in sourcegen.
     * @see SolutionArray.info()
     * @param rows  Maximum number of rendered rows.
     * @param width  Maximum width of rendered output.
     * @since New in %Cantera 3.2
     */
    string _info(int rows=10, int width=80) {
        return info({}, rows, width);
    }
 
    //! Return thermo/kinetics/transport manager used in the domain
    //! @since New in %Cantera 3.2.
    shared_ptr<Solution> phase() const {
        return m_solution;
    }
 
    //! Return the size of the solution vector (the product of #m_nv and #m_points).
    size_t size() const {
        return m_nv*m_points;
    }
 
    /**
     * Find the index of the first grid point in this domain, and
     * the start of its variables in the global solution vector.
     */
    void locate();
 
    //! Location of the start of the local solution vector in the global solution vector
    virtual size_t loc(size_t j = 0) const {
        return m_iloc;
    }
 
    //! The index of the first (that is, left-most) grid point belonging to this domain.
    size_t firstPoint() const {
        return m_jstart;
    }
 
    //! The index of the last (that is, right-most) grid point belonging to this domain.
    size_t lastPoint() const {
        return m_jstart + m_points - 1;
    }
 
    /**
     * Set the left neighbor to domain 'left.' Method 'locate' is called to
     * update the global positions of this domain and all those to its right.
     */
    void linkLeft(Domain1D* left) {
        m_left = left;
        locate();
    }
 
    //! Set the right neighbor to domain 'right.'
    void linkRight(Domain1D* right) {
        m_right = right;
    }
 
    //! Append domain 'right' to this one, and update all links.
    void append(Domain1D* right) {
        linkRight(right);
        right->linkLeft(this);
    }
 
    //! Return a pointer to the left neighbor.
    Domain1D* left() const {
        return m_left;
    }
 
    //! Return a pointer to the right neighbor.
    Domain1D* right() const {
        return m_right;
    }
 
    //! Value of component n at point j in the previous solution.
    double prevSoln(size_t n, size_t j) const {
        return m_slast[m_nv*j + n];
    }
 
    //! Specify an identifying tag for this domain.
    void setID(const string& s) {
        m_id = s;
    }
 
    //! Returns the identifying tag for this domain.
    string id() const {
        if (m_id != "") {
            return m_id;
        } else {
            return fmt::format("domain {}", m_index);
        }
    }
 
    //! Print the solution.
    //! @param x  Pointer to the local portion of the system state vector
    virtual void show(span<const double> x);
 
    //! Get the coordinate [m] of the point with local index `jlocal`
    double z(size_t jlocal) const {
        return m_z[jlocal];
    }
 
    //! Get the coordinate [m] of the first (leftmost) grid point in this domain
    double zmin() const {
        return m_z[0];
    }
 
    //! Get the coordinate [m] of the last (rightmost) grid point in this domain
    double zmax() const {
        return m_z[m_points - 1];
    }
 
    //! Access the array of grid coordinates [m]
    span<double> grid() {
        return m_z;
    }
 
    //! Access the array of grid coordinates [m]
    span<const double> grid() const {
        return m_z;
    }
 
    //! Set up initial grid.
    //! @since New in %Cantera 3.2.
    //! called to set up initial grid, and after grid refinement
    virtual void setupGrid(span<const double> z);
 
    //! Set up uniform grid.
    //! @param points  Number of grid points
    //! @param length  Length of domain
    //! @param start  Start position of domain (default=0.)
    //! @since New in %Cantera 3.2.
    void setupUniformGrid(size_t points, double length, double start=0.);
 
    /**
     * Writes some or all initial solution values into the global solution
     * array, beginning at the location pointed to by x. This method is called
     * by the Sim1D constructor, and allows default values or ones that have
     * been set locally prior to installing this domain into the container to be
     * written to the global solution vector.
     */
    virtual void _getInitialSoln(span<double> x);
 
    //! Initial value of solution component @e n at grid point @e j.
    virtual double initialValue(size_t n, size_t j);
 
    /**
     * In some cases, a domain may need to set parameters that depend on the
     * initial solution estimate. In such cases, the parameters may be set in
     * method _finalize. This method is called just before the Newton solver is
     * called, and the x array is guaranteed to be the local solution vector for
     * this domain that will be used as the initial guess. If no such parameters
     * need to be set, then method _finalize does not need to be overloaded.
     */
    virtual void _finalize(span<const double> x) {}
 
    /**
     * In some cases, for computational efficiency some properties (such as
     * transport coefficients) may not be updated during Jacobian evaluations.
     * Set this to `true` to force these properties to be updated even while
     * calculating Jacobian elements.
     */
    void forceFullUpdate(bool update) {
        m_force_full_update = update;
    }
 
    //! Set shared data pointer
    void setData(shared_ptr<vector<double>>& data) {
        m_state = data;
    }
 
protected:
    //! Retrieve meta data
    virtual AnyMap getMeta() const;
 
    //! Retrieve meta data
    virtual void setMeta(const AnyMap& meta);
 
    double m_press = -1.0; //!< pressure [Pa]
 
    shared_ptr<vector<double>> m_state; //!< data pointer shared from OneDim
 
    double m_rdt = 0.0; //!< Reciprocal of the time step
    size_t m_nv = 0; //!< Number of solution components
    size_t m_points; //!< Number of grid points
    vector<double> m_slast; //!< Solution vector at the last time step
    vector<double> m_max; //!< Upper bounds on solution components
    vector<double> m_min; //!< Lower bounds on solution components
    vector<double> m_rtol_ss; //!< Relative tolerances for steady mode
    vector<double> m_rtol_ts; //!< Relative tolerances for transient mode
    vector<double> m_atol_ss; //!< Absolute tolerances for steady mode
    vector<double> m_atol_ts; //!< Absolute tolerances for transient mode
    vector<double> m_z; //!< 1D spatial grid coordinates
 
    //! Parent OneDim simulation containing this and adjacent domains
    OneDim* m_container = nullptr;
 
    size_t m_index; //!< Left-to-right location of this domain
 
    //! Starting location within the solution vector for unknowns that
    //! correspond to this domain
    /*!
     * Remember there may be multiple domains associated with this problem
     */
    size_t m_iloc = 0;
 
    //! Index of the first point in this domain in the global point list.
    //! @see firstPoint(), lastPoint()
    size_t m_jstart = 0;
 
    Domain1D* m_left = nullptr; //!< Pointer to the domain to the left
    Domain1D* m_right = nullptr; //!< Pointer to the domain to the right
 
    //! Identity tag for the domain
    string m_id;
    unique_ptr<Refiner> m_refiner; //!< Refiner object used for placing grid points
    vector<string> m_name; //!< Names of solution components
    int m_bw = -1; //!< See bandwidth()
    bool m_force_full_update = false; //!< see forceFullUpdate()
    string m_jacobianMode = "auto"; //!< see setJacobianMode()
 
    //! Composite thermo/kinetics/transport handler
    shared_ptr<Solution> m_solution;
};
}
 
#endif