BoundaryValueProblem.h (Source)

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

// 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 BVP_H
#define BVP_H

#include "cantera/onedim.h"
#include <fstream>

/// 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
    std::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) {}
};

/**
 * 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 Cantera::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, double zmin, double zmax) :
        m_left(0), m_right(0), m_sim(0)
    {
        // Create the initial uniform grid
        Cantera::vector_fp z(np);
        int iz;
        for (iz = 0; iz < np; iz++) {
            z[iz] = zmin + iz*(zmax - zmin)/(np-1);
        }
        setupGrid(np, z.data());
        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, double* 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(size_t n, Component& c) {
        if (m_sim == 0) {
            start();
        }
        if (n >= m_nv) {
            throw Cantera::CanteraError("BoundaryValueProblem::setComponent",
                                        "Illegal solution component number");
        }
        // set the upper and lower bounds for this component
        setBounds(n, c.lower, c.upper);
        // set the error tolerances
        setSteadyTolerances(c.rtol, c.atol, n);
        setTransientTolerances(c.rtol, c.atol, n);
        // 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 = fmt::format("Component {}", 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(std::string filename = "output.csv",
                  bool dotitles = true, std::string ztitle = "z") const {
        std::ofstream f(filename);
        int np = nPoints();
        int nc = nComponents();
        int n, m;
        if (dotitles) {
            f << ztitle << ", ";
            for (m = 0; m < nc; m++) {
                f << componentName(m);
                if (m != nc - 1) {
                    f << ", ";
                }
            }
            f << std::endl;
        }
        for (n = 0; n < np; n++) {
            f << z(n) << ", ";
            for (m = 0; m < nc; m++) {
                f << m_sim->value(1, m, n);
                if (m != nc - 1) {
                    f << ", ";
                }
            }
            f << std::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 double initialValue(size_t n, size_t j) {
        return 0.0;
    }

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

    /**
     * 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 Cantera::Empty1D;
        m_right = new Cantera::Empty1D;
        std::vector<Cantera::Domain1D*> domains { m_left, this, m_right };

        // create the Sim1D instance that will control the
        // solution process
        m_sim = new Cantera::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.
     */
    double cdif2(const double* x, int n, int j) const {
        double c1 = value(x,n,j) - value(x,n,j-1);
        double 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.
    */
    double firstDeriv(const double* 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.
     */
    double rightFirstDeriv(const double* 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.
     */

    double leftFirstDeriv(const double* 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.
     */
    double centralFirstDeriv(const double* x, int n, int j) const {
        double c1 = value(x,n,j+1) - value(x,n,j-1);
        return c1/(z(j+1) - z(j-1));
    }

    //@}
};
}
#endif

blasius.cpp (Source)

/*!
 * @file blasius.cpp
 *
 * Blasius BVP Solver
 *
 * This example solves the Blasius boundary value problem for the velocity
 * profile of a laminar boundary layer over a flat plate. It defines class
 * BoundaryValueProblem, which provides a simplified interface to the boundary
 * value problem capabilities of Cantera.
 */

// 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.

#include "BoundaryValueProblem.h"

using Cantera::npos;

/**
 * This class solves the Blasius boundary value problem on the domain (0,L):
 * \f[
 *             \frac{d\zeta}{dz} = u.
 * \f]
 * \f[
 *             \frac{d^2u}{dz^2} + 0.5\zeta \frac{du}{dz} = 0.
 * \f]
 * with boundary conditions
 * \f[
 * \zeta(0) = 0, u(0) = 0, u(L) = 1.
 * \f]
 * Note that this is formulated as a system of two equations, with maximum
 * order of 2, rather than as a single third-order boundary value problem.
 * For reasons having to do with the band structure of the Jacobian, no
 * equation in the system should have order greater than 2.
 */
class Blasius : public BVP::BoundaryValueProblem
{
public:
    // This problem has two components (zeta and u)
    Blasius(int np, double L) : BVP::BoundaryValueProblem(2, np, 0.0, L) {
        // specify the component bounds, error tolerances, and names.
        BVP::Component A;
        A.lower = -200.0;
        A.upper = 200.0;
        A.rtol = 1.0e-12;
        A.atol = 1.0e-15;
        A.name = "zeta";
        setComponent(0, A); // zeta will be component 0

        BVP::Component B;
        B.lower = -200.0;
        B.upper = 200.0;
        B.rtol = 1.0e-12;
        B.atol = 1.0e-15;
        B.name = "u";
        setComponent(1, B); // u will be component 1
    }

    // destructor
    virtual ~Blasius() {}

    // specify guesses for the initial values. These can be anything
    // that leads to a converged solution.
    virtual double initialValue(size_t n, size_t j) {
        switch (n) {
        case 0:
            return 0.1*z(j);
        case 1:
            return 0.5*z(j);
        default:
            return 0.0;
        }
    }

    // Specify the residual function. This is where the ODE system and boundary
    // conditions are specified. The solver will attempt to find a solution
    // x so that rsd is zero.
    void eval(size_t jg, double* x, double* rsd, int* diag, double rdt) {
        size_t jpt = jg - firstPoint();
        size_t jmin, jmax;
        if (jg == npos) {  // evaluate all points
            jmin = 0;
            jmax = m_points - 1;
        } else {  // evaluate points for Jacobian
            jmin = std::max<size_t>(jpt, 1) - 1;
            jmax = std::min(jpt+1,m_points-1);
        }

        for (size_t j = jmin; j <= jmax; j++) {
            if (j == 0) {
                rsd[index(0,j)] = zeta(x,j);
                rsd[index(1,j)] = u(x,j);
            } else if (j == m_points - 1) {
                rsd[index(0,j)] = leftFirstDeriv(x,0,j) - u(x,j);
                rsd[index(1,j)] = u(x,j) - 1.0;
            } else {
                rsd[index(0,j)] = leftFirstDeriv(x,0,j) - u(x,j);
                rsd[index(1,j)] = cdif2(x,1,j) + 0.5*zeta(x,j)*centralFirstDeriv(x,1,j)
                                  - rdt*(value(x,1,j) - prevSoln(1,j));
                diag[index(1,j)] = 1;
            }
        }
    }

private:
    // for convenience only. Note that the compiler will inline these.
    double zeta(double* x, int j) {
        return value(x,0,j);
    }
    double u(double* x, int j) {
        return value(x,1,j);
    }
};

int main()
{
    try {
        // Specify a problem on (0,10), with an initial uniform grid of
        // 6 points.
        Blasius eqs(6, 10.0);
        // Solve the equations, refining the grid as needed
        eqs.solve(1);
        // write the solution to a CSV file.
        eqs.writeCSV("blasius.csv");
        return 0;
    } catch (Cantera::CanteraError& err) {
        std::cerr << err.what() << std::endl;
        return -1;
    }
}