research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

initial import

Browse source
This commit is contained in:
Dave Goodwin 2006-10-10 16:02:57 +00:00
parent 601abe3feb
commit 10add09884
3 changed files with 518 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

383
apps/bvp/BoundaryValueProblem.h Normal file
View file

@ -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 <cantera/Cantera.h>
#include <cantera/onedim.h>
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<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

15
apps/bvp/README Normal file
View file

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

120
apps/bvp/blasius.cpp Normal file
View file

@ -0,0 +1,120 @@
/// @file blasius.cpp
/// The Blasius boundary layer
#include <cantera/Cantera.h>
#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;
}
}