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<double> 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 @e 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 @e n at initial grid
* point @e j. The default is zero for all components at all
* grid points. Overload in derived classes to specify other
* choices for initial values.
*/
double initialValue(size_t n, size_t j) override {
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 @e n at point @e 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 @e n at point @e 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
/*!
* @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.
*
* Keywords: user-defined model
*/
// 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
}
// specify guesses for the initial values. These can be anything
// that leads to a converged solution.
double initialValue(size_t n, size_t j) override {
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) override {
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;
}
}