From 10add098847b7746c6b92d6cb659e1c8a202bb50 Mon Sep 17 00:00:00 2001 From: Dave Goodwin Date: Tue, 10 Oct 2006 16:02:57 +0000 Subject: [PATCH] initial import --- apps/bvp/BoundaryValueProblem.h | 383 ++++++++++++++++++++++++++++++++ apps/bvp/README | 15 ++ apps/bvp/blasius.cpp | 120 ++++++++++ 3 files changed, 518 insertions(+) create mode 100644 apps/bvp/BoundaryValueProblem.h create mode 100644 apps/bvp/README create mode 100644 apps/bvp/blasius.cpp diff --git a/apps/bvp/BoundaryValueProblem.h b/apps/bvp/BoundaryValueProblem.h new file mode 100644 index 000000000..6083fd42f --- /dev/null +++ b/apps/bvp/BoundaryValueProblem.h @@ -0,0 +1,383 @@ +/// @file BoundaryValueProblem.h +/// Simplified interface to the capabilities provided by Cantera to +/// solve boundary value problems. + +#ifndef BVP_H +#define BVP_H + +#include +#include +using namespace Cantera; + +/// 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 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 diff --git a/apps/bvp/README b/apps/bvp/README new file mode 100644 index 000000000..fb2bab387 --- /dev/null +++ b/apps/bvp/README @@ -0,0 +1,15 @@ +This example program solves the Blasius boundary value problem for the +velocity profile of a laminar boundary layer over a flat plate. It +uses class BoundaryValueProblem, which provides a simplified interface +to the boundary value problem capabilities of Cantera. + +To build this example, type "ctnew" to generate a demo c++ program and +a makefile (demo.mak) that is correctly configured for your Cantera +installation. It it is not on your path, you can find the ctnew script +in the "bin" directory of your Cantera installation directory. + +First make sure the Cantera demo works by typing "make -f demo.mak", +then "./demo". Assuming this works, now edit demo.mak and change the +line "OBJS = demo.o" to "OBJS = blasius.o". You can optionally change +the program name too. Now when you rebuild the executable and run it, +it will solve the blasius boundary layer problem. diff --git a/apps/bvp/blasius.cpp b/apps/bvp/blasius.cpp new file mode 100644 index 000000000..dd9215c38 --- /dev/null +++ b/apps/bvp/blasius.cpp @@ -0,0 +1,120 @@ +/// @file blasius.cpp +/// The Blasius boundary layer + +#include +#include "BoundaryValueProblem.h" + + +/** + * 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 doublereal initialValue(int n, int j) { + switch (n) { + case 0: + return 0.1*z(j); + case 1: + return 0.5*z(j); + default: + return 0.0; + } + } + + // Specify the residual. This is where the ODE system and boundary + // conditions are specified. The solver will attempt to find a solution + // x so that this function returns 0 for all n and j. + virtual doublereal residual(doublereal* x, int n, int j) { + + // if n = 0, return the residual for the first ODE + if (n == 0) { + if (isLeft(j)) // here we specify zeta(0) = 0 + return zeta(x,j); + else + // this implements d(zeta)/dz = u + return (zeta(x,j) - zeta(x,j-1))/(z(j)-z(j-1)) - u(x,j); + } + // if n = 1, then return the residual for the second ODE + else { + if (isLeft(j)) // here we specify u(0) = 0 + return u(x,j); + else if (isRight(j)) // and here we specify u(L) = 1 + return u(x,j) - 1.0; + else + // this implements the 2nd ODE + return cdif2(x,1,j) + 0.5*zeta(x,j)*centralFirstDeriv(x,1,j); + } + } + + +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, and print lots of diagnostic output (loglevel = 4) + eqs.solve(4); + // write the solution to a CSV file. + eqs.writeCSV(); + return 0; + } + catch (CanteraError) { + showErrors(cerr); + return -1; + } +} + +