cantera/Cantera/src/oneD/Resid1D.h

/**
 *  @file Resid1D.h
 *
 *  $Author$
 *  $Date$
 *  $Revision$
 *
 *  Copyright 2002 California Institute of Technology
 *
 */

#ifndef CT_RESID1D_H
#define CT_RESID1D_H


#include "../ctexceptions.h"
#include "../xml.h"
#include "refine.h"


namespace Cantera {

    // domain types
    const int cFlowType         = 101;
    const int cSurfType         = 102;
    const int cConnectorType    = 103;
    const int cInletType        = 104;
    const int cSymmType         = 105;
    const int cOutletType       = 106;

    class MultiJac;
    class OneDim;


    /**
     * Base class for one-dimensional domains.
     */
    class Resid1D {
    public:

        /**
         * Constructor.
         * @param nv      Number of variables at each grid point.
         * @param points  Number of grid points.
         */
        Resid1D(int nv=1, int points=1, 
            doublereal time = 0.0) : 
            m_time(time),
            m_container(0), 
            m_index(-1),
            m_type(0),
            m_iloc(0),
            m_jstart(0),
            m_left(0), 
            m_right(0),
            m_refiner(0) {
            resize(nv, points);
        }

        /// Destructor. Does nothing
        virtual ~Resid1D(){ delete m_refiner; }

        /// Domain type flag.
        const int domainType() { return m_type; }

        /**
         * The left-to-right location of this domain.
         */
        const int domainIndex() { return m_index; }

        /**
         * 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, int index){
            m_container = c;
            m_index = index;
        }

        /** 
         * Initialize. Base class method does nothing, but may be
         * overloaded. 
         */
        virtual void init(){}

        virtual void setInitialState(doublereal* xlocal = 0){}
        virtual void setState(int point, const doublereal* state, doublereal* 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(int nv, int np) {
            if (nv != m_nv || !m_refiner) {
                m_nv = nv;
                delete m_refiner;
                m_refiner = new Refiner(*this);
                cout << "created refiner with nv = " << m_nv << endl;
            }
            m_nv = nv;
            m_max.resize(m_nv, 0.0);
            m_min.resize(m_nv, 0.0);
            m_rtol.resize(m_nv, 0.0);
            m_atol.resize(m_nv, 0.0);
            m_points = np;
            m_z.resize(np, 0.0);
            m_slast.resize(m_nv * m_points, 0.0);
            locate();
        }

        Refiner& refiner() { return *m_refiner; }

        /// Number of components at each grid point.
        int nComponents() const { return m_nv; }

        /// Number of grid points in this domain.
        int nPoints() const { return m_points; }

        /// Name of the nth component. May be overloaded.
        virtual string componentName(int n) const { 
            return "component " + int2str(n); }

        /**
         * Set the lower and upper bounds for each solution component.
         */
        void setBounds(int nl, const doublereal* lower, 
            int nu, const doublereal* upper) {
            if (nl < m_nv || nu < m_nv)
                throw CanteraError("Resid1D::setBounds",
                    "wrong array size for solution bounds. "
                    "Size should be at least "+int2str(m_nv));
            copy(upper, upper + m_nv, m_max.begin());
            copy(lower, lower + m_nv, m_min.begin());
        }

        void setTolerances(int nr, const doublereal* rtol, 
            int na, const doublereal* atol) {
            if (nr < m_nv || na < m_nv)
                throw CanteraError("Resid1D::setTolerances",
                    "wrong array size for solution error tolerances. "
                    "Size should be at least "+int2str(m_nv));
            copy(rtol, rtol + m_nv, m_rtol.begin());
            copy(atol, atol + m_nv, m_atol.begin());
        }

        /// Relative tolerance of the nth component.
        doublereal rtol(int n) { return m_rtol[n]; }

        /// Absolute tolerance of the nth component.
        doublereal atol(int n) { return m_atol[n]; }

        /// Upper bound on the nth component.
        doublereal upperBound(int n) const { return m_max[n]; }

        /// Lower bound on the nth component
        doublereal lowerBound(int n) const { return m_min[n]; }

        /**
         * Prepare to do time stepping with time step dt. Copy the
         * internally-stored solution at the last time step to array
         * x0.
         */
        void initTimeInteg(doublereal dt, const doublereal* x0) {
            copy(x0 + loc(), x0 + loc() + size(), m_slast.begin());
            m_rdt = 1.0/dt;
        } 
        
        /**
         * Prepare to solve the steady-state problem.
         * Set the internally-stored reciprocal of the time step to 0,0
         */
        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 (e.g. the value of a constant has been changed,
         * so that the last-computed Jacobian is no longer valid.
         */
        void needJacUpdate();

        /**
         * Evaluate the steady-state residual at all points, even if in 
         * transient mode. Used to print diagnostic output.
         */
        void evalss(doublereal* x, doublereal* r, integer* mask) {
            eval(-1,x,r,mask,0.0);
        }

        /**
         * Evaluate the residual function at point j. If j < 0,
         * evaluate the residual function at all points.
         */         
        virtual void eval(int j, doublereal* x, doublereal* r, 
            integer* mask, doublereal rdt=0.0) {
            throw CanteraError("Resid1D::eval",
                "residual function not defined.");
        }

        /**
         * Does nothing.
         */
        virtual void update(doublereal* x) {}

        doublereal time() const { return m_time;}
        void incrementTime(doublereal dt) { m_time += dt; }
        size_t index(int n, int j) const { return m_nv*j + n; }
        doublereal value(doublereal* x, int n, int j) const {
            return x[index(n,j)];
        }

        virtual void setJac(MultiJac* jac){}
        virtual void save(XML_Node& o, doublereal* sol) {
            throw CanteraError("Resid1D::save","base class method called");
        }

        int 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() {

            if (m_left) {
                // there is a domain on the left, so the first grid point
                // in this domain is one more than the last one on the left
                m_jstart = m_left->lastPoint() + 1;

                // the starting location in the solution vector
                m_iloc = m_left->loc() + m_left->size();
            }
            else {
                // this is the left-most domain
                m_jstart = 0;
                m_iloc = 0;
            }
            // if there is a domain to the right of this one, then 
            // repeat this for it
            if (m_right) m_right->locate();
        }

        /**
         * Location of the start of the local solution vector in the global
         * solution vector,
         */
        virtual int loc(int j = 0) const { return m_iloc; }

        /**
         * The index of the first (i.e., left-most) grid point
         * belonging to this domain.
         */
        int firstPoint() const { return m_jstart; }

        /**
         * The index of the last (i.e., right-most) grid point
         * belonging to this domain.
         */
        int 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(Resid1D* left) { 
            m_left = left;
            locate();
        }

        /**
         * Set the right neighbor to domain 'right.' 
         */
        void linkRight(Resid1D* right) { m_right = right; }

        /**
         * Append domain 'right' to this one, and update all links.
         */
        void append(Resid1D* right) {
            linkRight(right);
            right->linkLeft(this);
        }

        /**
         * Return a pointer to the left neighbor.
         */
        Resid1D* left() const { return m_left; }

        /**
         * Return a pointer to the right neighbor.
         */
        Resid1D* right() const { return m_right; }

        /**
         * Value of component n at point j in the previous solution.
         */
        double prevSoln(int n, int j) const {
            return m_slast[m_nv*j + n];
        }
        
        /**
         * Specify an identifying tag for this domain.
         */
        void setID(const string& s) {m_id = s;}

        /**
         * Specify descriptive text for this domain.
         */
        void setDesc(const string& s) {m_desc = s;}

        virtual void getTransientMask(integer* mask){}

        virtual void showSolution_s(ostream& s, const doublereal* x) {}

        doublereal z(int jlocal) const {
            return m_z[jlocal];
        }
        doublereal zmin() const { return m_z[0]; }
        doublereal zmax() const { return m_z[m_points - 1]; }


        void setProfile(string name, doublereal* values, doublereal* soln) {
            int n, j;
            for (n = 0; n < m_nv; n++) {
                if (name == componentName(n)) {
                    for (j = 0; j < m_points; j++) {
                        soln[index(n, j) + m_iloc] = values[j];
                    }
                    return;
                }
            }
            throw CanteraError("Resid1D::setProfile",
                "unknown component: "+name);
        }

        vector_fp& grid() { return m_z; }
        const vector_fp& grid() const { return m_z; }
        doublereal grid(int point) { return m_z[point]; }

        virtual void setupGrid(int n, const doublereal* z) {}

        /**
         * Writes some or all initial solution values into array x,
         * which is the solution vector for this domain. This allows
         * initial values that have been set prior to installing this
         * domain into the container to be written to the global
         * solution vector.
         */
        virtual void _getInitialSoln(doublereal* x) {cout << "base class method _getInitialSoln called!" << endl;}

        /**
         * Perform any necessary domain-specific initialization using 
         * local solution vector x.
         */ 
        virtual void _finalize(const doublereal* x) {cout << "base class method _finalize called!" << endl;}

    protected:

        doublereal m_rdt;
        int m_nv;
        int m_points;
        vector_fp m_slast;
        doublereal m_time;
        vector_fp m_max;
        vector_fp m_min;
        vector_fp m_rtol;
        vector_fp m_atol;
        vector_fp m_z;
        OneDim* m_container;
        int m_index;
        int m_type;
        int m_iloc;
        int m_jstart;
        Resid1D *m_left, *m_right;
        string m_id, m_desc;
        Refiner* m_refiner;

    private:

    };
}

#endif