OneDim.h

Go to the documentation of this file.

/**
 * @file OneDim.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_ONEDIM_H
#define CT_ONEDIM_H
 
#include "Domain1D.h"
#include "cantera/base/global.h"
#include "cantera/numerics/SystemJacobian.h"
#include "cantera/numerics/SteadyStateSystem.h"
 
namespace Cantera
{
 
/**
 * Container class for multiple-domain 1D problems. Each domain is
 * represented by an instance of Domain1D.
 * @ingroup onedGroup
 */
class OneDim : public SteadyStateSystem
{
public:
    //! Default constructor
    OneDim() = default;
 
    //! Construct a OneDim container for the domains in the list *domains*.
    OneDim(span<const shared_ptr<Domain1D>> domains);
 
    //! Add a domain. Domains are added left-to-right.
    void addDomain(shared_ptr<Domain1D> d);
 
    //! Compute the weighted norm of a step vector
    //!
    //! The weighted norm of a step vector @f$ \Delta x @f$ is defined as
    //! @f[
    //!   ||\Delta x|| = \sqrt{ \frac{1}{N}
    //!       \sum_{d,n,j} \left( \frac{\Delta x_{d,n,j}}{w_{d,n}} \right)^2 }
    //! @f]
    //! where the error weight for solution component @f$ n @f$ in domain @f$ d @f$ is
    //! given by
    //! @f[
    //!   w_{d,n} = \frac{\epsilon_{{\rm r};d,n}}{J_d} \sum_j |x_{d,n,j}|
    //!           + \epsilon_{{\rm a};d,n}
    //! @f]
    //! Here, @f$ \epsilon_{{\rm r};d,n} @f$ is the relative error tolerance for
    //! component @f$ n @f$ in domain @f$ d @f$, and multiplies the average magnitude of
    //! solution component @f$ n @f$ in the domain. The second term,
    //! @f$ \epsilon_{{\rm a};d,n} @f$, is the absolute error tolerance for component
    //! @f$ n @f$ in domain @f$ d @f$. @f$ N @f$ is the total number of state variables
    //! across all domains and @f$ J_d @f$ is the number of grid points in domain
    //! @f$ d @f$.
    double weightedNorm(span<const double> step) const override;
 
    //! Number of domains.
    size_t nDomains() const {
        return m_dom.size();
    }
 
    //! Return a reference to domain i.
    Domain1D& domain(size_t i) const {
        return *m_dom[i];
    }
 
    //! Get the index of the domain named `name`.
    size_t domainIndex(const string& name) const;
 
    //! Check that the specified domain index is in range.
    //! Throws an exception if n is greater than nDomains()-1
    void checkDomainIndex(size_t n) const {
        if (n >= m_dom.size()) {
            throw IndexError("OneDim::checkDomainIndex", "domains", n, m_dom.size());
        }
    }
 
    //! Check that an array size is at least nDomains().
    //! Throws an exception if nn is less than nDomains(). Used before calls
    //! which take an array pointer.
    void checkDomainArraySize(size_t nn) const {
        if (m_dom.size() > nn) {
            throw ArraySizeError("OneDim::checkDomainArraySize", nn,
                                 m_dom.size());
        }
    }
 
    //! The index of the start of domain i in the solution vector.
    size_t start(size_t i) const {
        if (m_dom[i]->nComponents()) {
            return m_dom[i]->loc();
        } else {
            // Special case for domains with no solution components to avoid
            // spurious out-of-bounds memory access
            return 0;
        }
    }
 
    //! Pointer to left-most domain (first added).
    Domain1D* left() {
        return m_dom[0].get();
    }
 
    //! Pointer to right-most domain (last added).
    Domain1D* right() {
        return m_dom.back().get();
    }
 
    //! Number of solution components at global point `jg`.
    size_t nVars(size_t jg) {
        return m_nvars[jg];
    }
 
    //! Location in the solution vector of the first component of global point `jg`.
    size_t loc(size_t jg) {
        return m_loc[jg];
    }
 
    //! Return the domain, local point index, and component name for the i-th
    //! component of the global solution vector
    std::tuple<string, size_t, string> component(size_t i) const;
 
    string componentName(size_t i) const override;
    pair<string, string> componentTableHeader() const override;
    string componentTableLabel(size_t i) const override;
    double upperBound(size_t i) const override;
    double lowerBound(size_t i) const override;
 
    /**
     * Initialize all domains. On the first call, this methods calls the init
     * method of each domain, proceeding from left to right. Subsequent calls
     * do nothing.
     */
    void init();
 
    //! Total number of points.
    size_t points() {
        return m_pts;
    }
 
    void initTimeInteg(double dt, span<const double> x) override;
    void setSteadyMode() override;
 
    /**
     * Evaluate the multi-domain residual function
     *
     * @param j       if j != npos, only evaluate residual for points j-1, j,
     *                and j + 1; otherwise, evaluate at all grid points.
     * @param x       solution vector
     * @param r       on return, contains the residual vector
     * @param rdt     Reciprocal of the time step. if omitted, then
     *                  the default value is used.
     * @param count   Set to zero to omit this call from the statistics
     */
    void eval(size_t j, span<const double> x, span<double> r,
              double rdt=-1.0, int count=1);
 
    void eval(span<const double> x, span<double> r,
              double rdt=-1.0, int count=1) override
    {
        return eval(npos, x, r, rdt, count);
    }
 
    void evalJacobian(span<const double> x0) override;
 
    //! Return a pointer to the domain global point *i* belongs to.
    /*!
     * The domains are scanned right-to-left, and the first one with starting
     * location less or equal to i is returned.
     */
    Domain1D* pointDomain(size_t i);
 
    void resize() override;
    void resetBadValues(span<double> x) override;
 
    //! Write statistics about the number of iterations and Jacobians at each
    //! grid level
    /*!
     *  @param printTime  Boolean that indicates whether time should be printed
     *                    out The default is true. It's turned off for test
     *                    problems where we don't want to print any times
     */
    void writeStats(int printTime = 1);
 
    //! Return total grid size in each call to solve()
    //! @deprecated To be removed after %Cantera 4.0. Use solverStats() instead.
    vector<size_t> gridSizeStats() {
        warn_deprecated("OneDim::gridSizeStats",
            "To be removed after Cantera 4.0. Use solverStats() instead.");
        saveStats();
        auto& v = m_stats["grid_points"].asVector<long int>();
        return vector<size_t>(v.begin(), v.end());
    }
 
    //! Return wall-clock time spent evaluating Jacobians in each call to solve()
    //! @deprecated To be removed after %Cantera 4.0. Use solverStats() instead.
    vector<double> jacobianTimeStats() {
        warn_deprecated("OneDim::jacobianTimeStats",
            "To be removed after Cantera 4.0. Use solverStats() instead.");
        saveStats();
        return m_stats["jacobian_time"].asVector<double>();
    }
 
    //! Return wall-clock time spent on non-Jacobian function evaluations in each
    //! call to solve()
    //! @deprecated To be removed after %Cantera 4.0. Use solverStats() instead.
    vector<double> evalTimeStats() {
        warn_deprecated("OneDim::evalTimeStats",
            "To be removed after Cantera 4.0. Use solverStats() instead.");
        saveStats();
        return m_stats["residual_time"].asVector<double>();
    }
 
    //! Return number of Jacobian evaluations made in each call to solve()
    //! @deprecated To be removed after %Cantera 4.0. Use solverStats() instead.
    vector<int> jacobianCountStats() {
        warn_deprecated("OneDim::jacobianCountStats",
            "To be removed after Cantera 4.0. Use solverStats() instead.");
        saveStats();
        auto& v = m_stats["jacobian_evals"].asVector<long int>();
        return vector<int>(v.begin(), v.end());
    }
 
    //! Return number of non-Jacobian function evaluations made in each call to
    //! solve()
    //! @deprecated To be removed after %Cantera 4.0. Use solverStats() instead.
    vector<int> evalCountStats() {
        warn_deprecated("OneDim::evalCountStats",
            "To be removed after Cantera 4.0. Use solverStats() instead.");
        saveStats();
        auto& v = m_stats["residual_evals"].asVector<long int>();
        return vector<int>(v.begin(), v.end());
    }
 
    //! Return number of time steps taken in each call to solve()
    //! @deprecated To be removed after %Cantera 4.0. Use solverStats() instead.
    vector<int> timeStepStats() {
        warn_deprecated("OneDim::timeStepStats",
            "To be removed after Cantera 4.0. Use solverStats() instead.");
        saveStats();
        auto& v = m_stats["steps"].asVector<long int>();
        return vector<int>(v.begin(), v.end());
    }
 
    //! Access internal work array.
    const vector<double>& _workVector() const {
        return m_xnew;
    }
 
protected:
    size_t gridSize() const override { return m_pts; }
 
    //! All domains comprising the system
    vector<shared_ptr<Domain1D>> m_dom;
 
    //! All connector and boundary domains
    vector<shared_ptr<Domain1D>> m_connect;
 
    //! All bulk/flow domains
    vector<shared_ptr<Domain1D>> m_bulk;
 
    //! Indicates whether one-time initialization for each domain has been completed.
    bool m_init = false;
 
    //! Number of variables at each point, across all domains. Length points().
    //! Accessed with nVars().
    vector<size_t> m_nvars;
 
    //! Location in the state vector of the first component of each point, across all
    //! domains. Accessed with loc().
    vector<size_t> m_loc;
 
    //! Domain, grid point, and component indices for each element of the global state
    //! vector. Length size()
    vector<std::tuple<size_t, size_t, size_t>>m_componentInfo;
 
    //! Total number of points.
    size_t m_pts = 0;
 
};
 
}
 
#endif