/// @file BoundaryValueProblem.h
/// Simplified interface to the capabilities provided by Cantera to
/// solve boundary value problems.

#ifndef BVP_H
#define BVP_H

#include <cantera/Cantera.h>
#include <cantera/onedim.h>
using namespace Cantera;
using namespace std;

/// Namespace for the boundary value problem package.
namespace BVP {

    // default grid refinement parameters
    const double max_grid_ratio  = 4.0;   ///< max ratio of neighboring grid intervals
    const double max_delta       = 0.01;  ///< max difference in function values
    const double max_delta_slope = 0.02;  ///< max difference in slopes
    const double prune           = 0.000; ///< don't remove grid points

    /**
     * Used to specify component-specific options for method
     * setComponent of method BoundaryValueProblem. An instance of
     * class Component should be created for each solution component,
     * and its values set appropriately.
     */
    class Component {
    public:
        double lower;  ///< lower bound
        double upper;  ///< upper bound
        double rtol;   ///< relative error tolerance
        double atol;   ///< absolute error tolerance
        bool refine;   ///< make this component active for grid refinement
        string name;   ///< component name

        /**
         * Constructor. Sets default values.
         */
        Component() : lower(0.0), upper(1.0), rtol(1.0e-9), atol(1.0e-12), 
                      refine(true), name("") {}
    };


    /**
     * Exception thrown for illegal parameter values when setting up
     * the problem.
     */
    class BVP_Error {
    public:
        /**
         * Constructor. Write an error message.
         */
        BVP_Error(string msg) {writelog("BVP Error: "+msg+"\n");}
    };


    /**
     * Base class for boundary value problems. This class is designed
     * to provide a simplified interface to the capabilities Cantera
     * provides to solve boundary value problems. Classes for specific 
     * boundary value problems should be derived from this one.
     *
     * Class BoundaryValueProblem derives from Cantera's Domain1D
     * class.
     */
    class BoundaryValueProblem : public Domain1D {

    public:

        /**
         * Constructor. This constructor begins with a uniform grid of
         * np points starting at zmin, and ending at zmax.
         *
         * @param nv Number of solution components
         * @param np Number of grid points in initial grid
         * @param zmin Location of left-hand side of domain
         * @param zmax Location of right-hand side of domain
         */
        BoundaryValueProblem(int nv, int np, 
            doublereal zmin, doublereal zmax) : 
            m_left(0), m_right(0), m_sim(0) {

            // Create the initial uniform grid
            vector_fp z(np);
            int iz;
            for (iz = 0; iz < np; iz++) z[iz] = zmin + iz*(zmax - zmin)/(np-1);
            setupGrid(np, DATA_PTR(z));
            resize(nv, np);
        }

        /** 
         * Constructor. This alternate constructor starts with a
         * specified grid, unlike the above that uses a uniform grid
         * to start. The array z must contain the z coordinates of np 
         * grid points. 
         */
        BoundaryValueProblem(int nv, int np, 
            doublereal* z) : 
            m_left(0), m_right(0), m_sim(0) {

            setupGrid(np, z);
            resize(nv, np);
        }


        /**
         * Destructor. Deletes the dummy terminator domains, and the
         * solver.
         */ 
        virtual ~BoundaryValueProblem() {
            delete m_left;
            delete m_right;
            delete m_sim;
        }
        

        /** 
         *  Set parameters and options for solution component \a n.
         *  This method should be invoked for each solution component
         *  before calling 'solve'. The parameter values should first
         *  be set by creating an instance of class Component, and 
         *  setting its member data appropriately.
         *
         *  @param n Component number. 
         *  @param c Component parameter values
         */
        void setComponent(int n, Component& c) {
            if (m_sim == 0) start();
            if (n < 0 || n >= m_nv) 
                throw BVP_Error("Illegal solution component number");
            // set the upper and lower bounds for this component
            setBounds(n, c.lower, c.upper);
            // set the error tolerances
            setTolerances(n, c.rtol, c.atol);
            // specify whether this component should be considered in 
            // refining the grid
            m_refiner->setActive(n, c.refine);
            // set a default name if one has not been entered
            if (c.name == "") c.name = "Component "+int2str(n);
            setComponentName(n, c.name);
        }


        /**
         * Solve the boundary value problem.
         * @param loglevel controls amount of diagnostic output.
         */
        void solve(int loglevel=0) {
            if (m_sim == 0) start();
            bool refine = true;
            m_sim->solve(loglevel, refine);
        }

            
        /**
         * Write the solution to a CSV file.
         * @param filename CSV file name.
         * @param ztitle Title for 'z' column.
         * @param dotitles If true, begin with a row of column titles.
         */
        void writeCSV(string filename = "output.csv", 
            bool dotitles = true, string ztitle = "z") const {
            ofstream f(filename.c_str());
            int np = nPoints();
            int nc = nComponents();
            int n, m;
            if (dotitles) {
                f << ztitle << ", ";
                for (m = 0; m < nc; m++) f << componentName(m) << ", ";
                f << endl;
            }
            for (n = 0; n < np; n++) {
                f << z(n) << ", ";
                for (m = 0; m < nc; m++) {
                    f << m_sim->value(1, m, n) << ", ";
                }
                f << endl;
            }
        }

        /**
         * Initial value of solution component \a n at initial grid
         * point \a j. The default is zero for all components at all
         * grid points. Overload in derived classes to specify other
         * choices for initial values.
         */
        virtual doublereal initialValue(int n, int j) { return 0.0;}


        /**
         * Value of component \a m at point \a j. This method is used
         * to access solution values once a converged solution has been 
         * attained. 
         */
        double v(int m, int j) const { return m_sim->value(1,m,j); }


    protected:

        Domain1D* m_left;  ///< dummy terminator
        Domain1D* m_right; ///< dummy terminator
        Sim1D* m_sim;      ///< controller for solution


        /**
         * True if n is the index of the left-most grid point (zero),
         * false otherwise.
         */
        bool isLeft(int n) const { return (n == 0); }

        /**
         * True if \a n is the index of the right-most grid point, false
         * otherwise.
         */
        bool isRight(int n) const { return (n == nPoints() - 1); }

        /**
         * Set up the problem. Creates the solver instance, and sets
         * default grid refinement parameters. This method is called
         * internally, and does not need to be invoked explicitly in
         * derived classes.
         */ 
        void start() {

            // Add dummy terminator domains on either side of this one. 
            m_left = new Empty1D;  
            m_right = new Empty1D;
            vector<Domain1D*> domains;
            domains.push_back(m_left);
            domains.push_back(this);
            domains.push_back(m_right);

            // create the Sim1D instance that will control the
            // solution process
            m_sim = new Sim1D(domains);

            // set default grid refinement parameters
            m_sim->setRefineCriteria(1, max_grid_ratio, max_delta, 
                max_delta_slope, prune);
        }


        /**
         * @name Trial Solution Derivatives
         * These methods return
         * derivatives of individual components at specified grid
         * points, using a given trial solution.  They are designed
         * for use in writing overloaded versions of method 'residual'
         * in derived classes.
         */

        //@{

        /**
         * This method is provided for use in method residual when 
         * central-differenced second derivatives are needed.
         * @param x The current trial solution vector.
         * @param n Component index.
         * @param j Grid point number.
         */
        doublereal cdif2(const doublereal* x, int n, int j) const {
            doublereal c1 = value(x,n,j) - value(x,n,j-1);
            doublereal c2 = value(x,n,j+1) - value(x,n,j);
            return 2.0*(c2/(z(j+1) - z(j)) - c1/(z(j) - z(j-1)))/
                (z(j+1) - z(j-1));
        }

        /**
         * The first derivative of solution component n at point j.
         * If type is -1, the first derivative is computed using the
         * value to the left of point j, if it is +1 then the
         * value to the right is used, and if it is zero (default) a
         * central-differenced first derivative is computed.
        */
        doublereal firstDeriv(const doublereal*x, int n, int j, 
            int type = 0) const {
            switch (type) {
            case -1:
                return leftFirstDeriv(x, n, j);
            case 1:
                return rightFirstDeriv(x, n, j);
            default:
                return centralFirstDeriv(x, n, j);
            }
        }

        /**
         * First derivative of component \a n at point \a j. The derivative
         * is formed to the right of point j, using values at point j
         * and point j + 1.
         */
        doublereal rightFirstDeriv(const doublereal* x, int n, int j) const {
            return (value(x,n,j+1) - value(x,n,j))/(z(j+1) - z(j));
        }

        /**
         * First derivative of component \a n at point \a j. The derivative
         * is formed to the left of point j, using values at point j
         * and point j - 1.
         */

        doublereal leftFirstDeriv(const doublereal* x, int n, int j) const {
            return (value(x,n,j) - value(x,n,j-1))/(z(j) - z(j-1));
        }

        /**
         * This method is provided for use in method residual when 
         * central-differenced first derivatives are needed.
         * @param x The current trial solution vector.
         * @param n Component index.
         * @param j Grid point number.
         */
        doublereal centralFirstDeriv(const doublereal* x, int n, int j) const {
            doublereal c1 = value(x,n,j+1) - value(x,n,j-1);
            return c1/(z(j+1) - z(j-1));
        }


        /**
         * This method is provided for use in method residual when 
         * central-differenced second derivatives are needed.
         * @param x The current trial solution vector.
         * @param n Component index.
         * @param j Grid point number.
         */
        doublereal central_secondDeriv(const doublereal* x, 
            int n, int j) const {
            doublereal c1 = leftFirstDeriv(x, n, j);
            doublereal c2 = rightFirstDeriv(x, n, j);
            return 2.0*(c2 - c1)/(z(j+1) - z(j-1));
        }

        /**
         * This method is provided for use in method residual when 
         * central-differenced evaluation of terms like 
         * \f[
         * \frac{d}{dz}\left(g \frac{df}{dz}\right)
         * \f]
         * is required.
         * @param x The current trial solution vector.
         * @param g The array of g values at the grid points.
         * @param n Component index.
         * @param j Grid point number.
         */
        doublereal central_Deriv_G_Deriv(const doublereal* x, 
            const doublereal* g, int n, int j) const {
            doublereal c1 = 0.5*(g[j] + g[j-1])*leftFirstDeriv(x, n, j);
            doublereal c2 = 0.5*(g[j+1] + g[j])*rightFirstDeriv(x, n, j);
            return 2.0*(c2 - c1)/(z(j+1) - z(j-1));
        }

        /** 
         * Value of component m between points j and j + 1. This is
         * computed as the mean of the values at j and j + 1.
         */
        doublereal midpointSolution(const doublereal* x, int m, int j) const {
            return 0.5*(value(x,m,j) + value(x,m,j+1));
        }

        /**
         * This method is provided for use in method residual when 
         * central-differenced evaluation of terms like 
         * \f[
         * \frac{d}{dz}\left(f_m \frac{df_n}{dz}\right)
         * \f]
         * is required.
         * @param x The current trial solution vector.
         * @param n Solution component for \f$ f_n \f$
         * @param m Solution component for \f$ f_m \f$
         * @param j Grid point number.
         */
        doublereal central_Deriv_S_Deriv(const doublereal* x, 
            int n, int m, int j) const {
            doublereal c1 = midpointSolution(x,m,j-1)*leftFirstDeriv(x, n, j);
            doublereal c2 = midpointSolution(x,m,j)*rightFirstDeriv(x, n, j);
            return 2.0*(c2 - c1)/(z(j+1) - z(j-1));
        }
        //@}

            
    };
}
#endif