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

Clean up Doxygen docs for 1D simulation classes

Browse source
This commit is contained in:
Ray Speth 2013-06-05 03:09:03 +00:00
parent 39770d6a48
commit f75ba112f9
16 changed files with 294 additions and 655 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

177
include/cantera/oneD/Domain1D.h
View file

@ -8,13 +8,11 @@
#ifndef CT_DOMAIN1D_H
#define CT_DOMAIN1D_H
#include "cantera/base/xml.h"
#include "cantera/base/stringUtils.h"
#include "cantera/base/ctexceptions.h"
#include "refine.h"
namespace Cantera
{
@ -32,18 +30,17 @@ const int cPorousType = 109;
class MultiJac;
class OneDim;
/**
* Base class for one-dimensional domains.
*/
class Domain1D
{
public:
/**
* Constructor.
* @param nv Number of variables at each grid point.
* @param points Number of grid points.
* @param time (unused)
*/
Domain1D(size_t nv=1, size_t points=1,
doublereal time = 0.0) :
@ -62,38 +59,31 @@ public:
resize(nv, points);
}
/// Destructor. Does nothing
virtual ~Domain1D() {
delete m_refiner;
}
/// Domain type flag.
//! Domain type flag.
int domainType() {
return m_type;
}
/**
* The left-to-right location of this domain.
*/
//! The left-to-right location of this domain.
size_t domainIndex() {
return m_index;
}
/**
* True if the domain is a connector domain.
*/
//! True if the domain is a connector domain.
bool isConnector() {
return (m_type >= cConnectorType);
}
/**
* The container holding this domain.
*/
//! The container holding this domain.
const OneDim& container() const {
return *m_container;
}
/**
/*!
* Specify the container object for this domain, and the
* position of this domain in the list.
*/
@ -102,35 +92,31 @@ public:
m_index = index;
}
/*
* Set the Jacobian bandwidth. See the discussion of method bandwidth.
*/
//! Set the Jacobian bandwidth. See the discussion of method bandwidth().
void setBandwidth(int bw = -1) {
m_bw = bw;
}
//! Set the Jacobian bandwidth for this domain.
/**
* Set the Jacobian bandwidth for this domain. When class
* OneDim computes the bandwidth of the overall multi-domain
* problem (in OneDim::resize()), it calls this method for the
* bandwidth of each domain. If setBandwidth has not been
* called, then a negative bandwidth is returned, in which
* case OneDim assumes that this domain is dense -- that is,
* at each point, all components depend on the value of all
* other components at that point. In this case, the bandwidth
* is bw = 2*nComponents() - 1. However, if this domain
* contains some components that are uncoupled from other
* components at the same point, then this default bandwidth
* may greatly overestimate the true bandwidth, with a
* substantial penalty in performance. For such domains, use
* method setBandwidth to specify the bandwidth before passing
* this domain to the Sim1D or OneDim constructor.
* When class OneDim computes the bandwidth of the overall multi-domain
* problem (in OneDim::resize()), it calls this method for the bandwidth
* of each domain. If setBandwidth has not been called, then a negative
* bandwidth is returned, in which case OneDim assumes that this domain is
* dense -- that is, at each point, all components depend on the value of
* all other components at that point. In this case, the bandwidth is bw =
* 2*nComponents() - 1. However, if this domain contains some components
* that are uncoupled from other components at the same point, then this
* default bandwidth may greatly overestimate the true bandwidth, with a
* substantial penalty in performance. For such domains, use method
* setBandwidth to specify the bandwidth before passing this domain to the
* Sim1D or OneDim constructor.
*/
size_t bandwidth() {
return m_bw;
}
/**
/*!
* Initialize. This method is called by OneDim::init() for
* each domain once at the beginning of a simulation. Base
* class method does nothing, but may be overloaded.
@ -140,7 +126,7 @@ public:
virtual void setInitialState(doublereal* xlocal = 0) {}
virtual void setState(size_t point, const doublereal* state, doublereal* x) {}
/**
/*!
* Resize the domain to have nv components and np grid points.
* This method is virtual so that subclasses can perform other
* actions required to resize the domain.
@ -168,12 +154,12 @@ public:
locate();
}
/// Return a reference to the grid refiner.
//! Return a reference to the grid refiner.
Refiner& refiner() {
return *m_refiner;
}
/// Number of components at each grid point.
//! Number of components at each grid point.
size_t nComponents() const {
return m_nv;
}
@ -195,7 +181,7 @@ public:
}
}
/// Number of grid points in this domain.
//! Number of grid points in this domain.
size_t nPoints() const {
return m_points;
}
@ -217,7 +203,7 @@ public:
}
}
/// Name of the nth component. May be overloaded.
//! Name of the nth component. May be overloaded.
virtual std::string componentName(size_t n) const {
if (m_name[n] != "") {
return m_name[n];
@ -236,7 +222,7 @@ public:
}
}
/// index of component with name \a name.
//! index of component with name \a name.
size_t componentIndex(const std::string& name) const {
size_t nc = nComponents();
for (size_t n = 0; n < nc; n++) {
@ -248,9 +234,7 @@ public:
"no component named "+name);
}
/**
* Set the lower and upper bounds for each solution component.
*/
//! Set the lower and upper bounds for each solution component.
void setBounds(size_t nl, const doublereal* lower,
size_t nu, const doublereal* upper) {
if (nl < m_nv || nu < m_nv)
@ -266,71 +250,69 @@ public:
m_max[n] = upper;
}
/// set the error tolerances for all solution components.
//! set the error tolerances for all solution components.
void setTolerances(size_t nr, const doublereal* rtol,
size_t na, const doublereal* atol, int ts = 0);
/// set the error tolerances for solution component \a n.
//! set the error tolerances for solution component \a n.
void setTolerances(size_t n, doublereal rtol, doublereal atol, int ts = 0);
/// set scalar error tolerances. All solution components will
/// have the same relative and absolute error tolerances.
//! set scalar error tolerances. All solution components will
//! have the same relative and absolute error tolerances.
void setTolerances(doublereal rtol, doublereal atol,int ts=0);
void setTolerancesTS(doublereal rtol, doublereal atol);
void setTolerancesSS(doublereal rtol, doublereal atol);
/// Relative tolerance of the nth component.
//! Relative tolerance of the nth component.
doublereal rtol(size_t n) {
return (m_rdt == 0.0 ? m_rtol_ss[n] : m_rtol_ts[n]);
}
/// Absolute tolerance of the nth component.
//! Absolute tolerance of the nth component.
doublereal atol(size_t n) {
return (m_rdt == 0.0 ? m_atol_ss[n] : m_atol_ts[n]);
}
/// Upper bound on the nth component.
//! Upper bound on the nth component.
doublereal upperBound(size_t n) const {
return m_max[n];
}
/// Lower bound on the nth component
//! Lower bound on the nth component
doublereal lowerBound(size_t n) const {
return m_min[n];
}
/**
* Prepare to do time stepping with time step dt. Copy the
* internally-stored solution at the last time step to array
* x0.
/*!
* Prepare to do time stepping with time step dt. Copy the internally-
* stored solution at the last time step to array x0.
*/
void initTimeInteg(doublereal dt, const doublereal* x0) {
std::copy(x0 + loc(), x0 + loc() + size(), m_slast.begin());
m_rdt = 1.0/dt;
}
/**
* Prepare to solve the steady-state problem.
* Set the internally-stored reciprocal of the time step to 0,0
/*!
* Prepare to solve the steady-state problem. Set the internally-stored
* reciprocal of the time step to 0,0
*/
void setSteadyMode() {
m_rdt = 0.0;
}
/// True if in steady-state mode
//! True if in steady-state mode
bool steady() {
return (m_rdt == 0.0);
}
/// True if not in steady-state mode
//! True if not in steady-state mode
bool transient() {
return (m_rdt != 0.0);
}
/**
/*!
* Set this if something has changed in the governing
* equations (e.g. the value of a constant has been changed,
* so that the last-computed Jacobian is no longer valid.
@ -338,7 +320,7 @@ public:
*/
void needJacUpdate();
/**
/*!
* Evaluate the steady-state residual at all points, even if in
* transient mode. Used only to print diagnostic output.
*/
@ -349,9 +331,13 @@ public:
//! Evaluate the residual function at point j. If j == npos,
//! evaluate the residual function at all points.
/*!
* @param j Grid point j
* @param x Soln vector. This is the input.
* @param r residual this is the output.
* @param j Grid point at which to update the residual
* @param[in] x State vector
* @param[out] r residual vector
* @param[out] mask Boolean mask indicating whether each solution
* component has a time derivative (1) or not (0).
* @param[in] rdt Reciprocal of the timestep (`rdt=0` implies steady-
* state.)
*/
virtual void eval(size_t j, doublereal* x, doublereal* r,
integer* mask, doublereal rdt=0.0);
@ -367,9 +353,6 @@ public:
m_td[n] = 0;
}
/**
* Does nothing.
*/
virtual void update(doublereal* x) {}
doublereal time() const {
@ -407,7 +390,7 @@ public:
* classes should call the base class method in addition to restoring
* their own data.
*
* @param o XML_Node for this domain
* @param dom XML_Node for this domain
* @param soln Current value of the solution vector, local to this object.
* @param loglevel 0 to suppress all output; 1 to show warnings; 2 for
* verbose output
@ -468,54 +451,41 @@ public:
}
/**
* Set the left neighbor to domain 'left.' Method 'locate' is
* called to update the global positions of this domain and
* all those to its right.
* Set the left neighbor to domain 'left.' Method 'locate' is called to
* update the global positions of this domain and all those to its right.
*/
void linkLeft(Domain1D* left) {
m_left = left;
locate();
}
/**
* Set the right neighbor to domain 'right.'
*/
//! Set the right neighbor to domain 'right.'
void linkRight(Domain1D* right) {
m_right = right;
}
/**
* Append domain 'right' to this one, and update all links.
*/
//! Append domain 'right' to this one, and update all links.
void append(Domain1D* right) {
linkRight(right);
right->linkLeft(this);
}
/**
* Return a pointer to the left neighbor.
*/
//! Return a pointer to the left neighbor.
Domain1D* left() const {
return m_left;
}
/**
* Return a pointer to the right neighbor.
*/
//! Return a pointer to the right neighbor.
Domain1D* right() const {
return m_right;
}
/**
* Value of component n at point j in the previous solution.
*/
//! Value of component n at point j in the previous solution.
double prevSoln(size_t n, size_t j) const {
return m_slast[m_nv*j + n];
}
/**
* Specify an identifying tag for this domain.
*/
//! Specify an identifying tag for this domain.
void setID(const std::string& s) {
m_id = s;
}
@ -528,9 +498,7 @@ public:
}
}
/**
* Specify descriptive text for this domain.
*/
//! Specify descriptive text for this domain.
void setDesc(const std::string& s) {
m_desc = s;
}
@ -541,6 +509,8 @@ public:
virtual void getTransientMask(integer* mask) {}
virtual void showSolution_s(std::ostream& s, const doublereal* x) {}
//! Print the solution.
virtual void showSolution(const doublereal* x);
doublereal z(size_t jlocal) const {
@ -553,7 +523,6 @@ public:
return m_z[m_points - 1];
}
void setProfile(const std::string& name, doublereal* values, doublereal* soln) {
for (size_t n = 0; n < m_nv; n++) {
if (name == componentName(n)) {
@ -577,6 +546,7 @@ public:
return m_z[point];
}
//! called to set up initial grid, and after grid refinement
virtual void setupGrid(size_t n, const doublereal* z);
void setGrid(size_t n, const doublereal* z);
@ -591,9 +561,7 @@ public:
*/
virtual void _getInitialSoln(doublereal* x);
/**
* Initial value of solution component \a n at grid point \a j.
*/
//! Initial value of solution component \a n at grid point \a j.
virtual doublereal initialValue(size_t n, size_t j);
/**
@ -614,7 +582,6 @@ public:
bool m_adiabatic;
protected:
doublereal m_rdt;
size_t m_nv;
size_t m_points;
@ -629,11 +596,10 @@ protected:
size_t m_index;
int m_type;
//! Starting location within the solution vector for unknowns
//! that correspond to this domain
//! Starting location within the solution vector for unknowns that
//! correspond to this domain
/*!
* Remember there may be multiple domains associated with
* this problem
* Remember there may be multiple domains associated with this problem
*/
size_t m_iloc;
@ -648,9 +614,6 @@ protected:
vector_int m_td;
std::vector<std::string> m_name;
int m_bw;
private:
};
}

45
include/cantera/oneD/Inlet1D.h
View file

@ -38,10 +38,8 @@ const int RightInlet = -1;
class Bdry1D : public Domain1D
{
public:
Bdry1D();
/// Initialize.
virtual void init() {
_init(1);
}
@ -93,7 +91,6 @@ public:
virtual void setupGrid(size_t n, const doublereal* z) {}
protected:
void _init(size_t n);
StFlow* m_flow_left, *m_flow_right;
@ -119,14 +116,7 @@ private:
*/
class Inlet1D : public Bdry1D
{
public:
/**
* Constructor. Create a new Inlet1D instance. If invoked
* without parameters, a left inlet (facing right) is
* constructed).
*/
Inlet1D() : Bdry1D(), m_V0(0.0), m_nsp(0), m_flow(0) {
m_type = cInletType;
m_xstr = "";
@ -143,7 +133,6 @@ public:
return m_V0;
}
virtual void showSolution(const doublereal* x) {
char buf[80];
sprintf(buf, " Mass Flux: %10.4g kg/m^2/s \n", m_mdot);
@ -187,7 +176,6 @@ public:
virtual void restore(const XML_Node& dom, doublereal* soln, int loglevel);
protected:
int m_ilr;
doublereal m_V0;
size_t m_nsp;
@ -196,15 +184,12 @@ protected:
StFlow* m_flow;
};
/**
* A terminator that does nothing.
*/
class Empty1D : public Domain1D
{
public:
Empty1D() : Domain1D() {
m_type = cEmptyType;
}
@ -223,9 +208,6 @@ public:
virtual void _getInitialSoln(doublereal* x) {
x[0] = 0.0;
}
protected:
};
/**
@ -234,7 +216,6 @@ protected:
*/
class Symm1D : public Bdry1D
{
public:
Symm1D() : Bdry1D() {
@ -256,19 +237,15 @@ public:
virtual void _getInitialSoln(doublereal* x) {
x[0] = m_temp;
}
protected:
};
/**
* An outlet.
*/
class Outlet1D : public Bdry1D
{
public:
Outlet1D() : Bdry1D() {
m_type = cOutletType;
}
@ -288,23 +265,15 @@ public:
virtual void _getInitialSoln(doublereal* x) {
x[0] = m_temp;
}
protected:
};
/**
* An outlet with specified composition.
*/
class OutletRes1D : public Bdry1D
{
public:
/**
* Constructor.
*/
OutletRes1D() : Bdry1D(), m_nsp(0), m_flow(0) {
m_type = cOutletResType;
m_xstr = "";
@ -337,14 +306,12 @@ public:
virtual void restore(const XML_Node& dom, doublereal* soln, int loglevel);
protected:
size_t m_nsp;
vector_fp m_yres;
std::string m_xstr;
StFlow* m_flow;
};
/**
* A non-reacting surface. The axial velocity is zero
* (impermeable), as is the transverse velocity (no slip). The
@ -353,9 +320,7 @@ protected:
*/
class Surf1D : public Bdry1D
{
public:
Surf1D() : Bdry1D() {
m_type = cSurfType;
}
@ -389,21 +354,14 @@ public:
writelog(buf);
writelog("\n");
}
protected:
};
/**
* A reacting surface.
*
*/
class ReactingSurf1D : public Bdry1D
{
public:
ReactingSurf1D() : Bdry1D(),
m_kin(0), m_surfindex(0), m_nsp(0) {
m_type = cSurfType;
@ -455,7 +413,6 @@ public:
}
protected:
InterfaceKinetics* m_kin;
SurfPhase* m_sphase;
size_t m_surfindex, m_nsp;

32
include/cantera/oneD/MultiJac.h
View file

@ -19,17 +19,12 @@ namespace Cantera
/**
* Class MultiJac evaluates the Jacobian of a system of equations
* defined by a residual function supplied by an instance of class
* 'OneDim.' The residual function may consist of several linked
* OneDim. The residual function may consist of several linked
* 1D domains, with different variables in each domain.
*/
class MultiJac : public BandMatrix
{
public:
/**
* Constructor.
*/
MultiJac(OneDim& r);
/**
@ -40,36 +35,30 @@ public:
*/
void eval(doublereal* x0, doublereal* resid0, double rdt);
/**
* Elapsed CPU time spent computing the Jacobian.
*/
//! Elapsed CPU time spent computing the Jacobian.
doublereal elapsedTime() const {
return m_elapsed;
}
/// Number of Jacobian evaluations.
//! Number of Jacobian evaluations.
int nEvals() const {
return m_nevals;
}
/**
* Number of times 'incrementAge' has been called since the
* last evaluation
*/
//! Number of times 'incrementAge' has been called since the last
//! evaluation
int age() const {
return m_age;
}
/**
* Increment the Jacobian age.
*/
//! Increment the Jacobian age.
void incrementAge() {
m_age++;
}
void updateTransient(doublereal rdt, integer* mask);
/// Set the age.
//! Set the Jacobian age.
void setAge(int age) {
m_age = age;
}
@ -81,11 +70,10 @@ public:
void incrementDiagonal(int j, doublereal d);
protected:
//! Residual evaluator for this jacobian
/*!
* This is a pointer to the residual evaluator. This
* object isn't owned by this jacobian object.
* This is a pointer to the residual evaluator. This object isn't owned
* by this jacobian object.
*/
OneDim* m_resid;
@ -102,5 +90,3 @@ protected:
}
#endif

35
include/cantera/oneD/MultiNewton.h
View file

@ -20,9 +20,7 @@ namespace Cantera
*/
class MultiNewton
{
public:
MultiNewton(int sz);
virtual ~MultiNewton();
@ -30,21 +28,40 @@ public:
return m_n;
}
/// Compute undamped step
//! Compute the undamped Newton step. The residual function is evaluated
//! at `x`, but the Jacobian is not recomputed.
void step(doublereal* x, doublereal* step,
OneDim& r, MultiJac& jac, int loglevel);
/// Compute factor to keep all components in bounds.
/**
* Return the factor by which the undamped Newton step 'step0'
* must be multiplied in order to keep all solution components in
* all domains between their specified lower and upper bounds.
*/
doublereal boundStep(const doublereal* x0, const doublereal* step0,
const OneDim& r, int loglevel);
/**
* On entry, step0 must contain an undamped Newton step for the
* solution x0. This method attempts to find a damping coefficient
* such that the next undamped step would have a norm smaller than
* that of step0. If successful, the new solution after taking the
* damped step is returned in x1, and the undamped step at x1 is
* returned in step1.
*/
int dampStep(const doublereal* x0, const doublereal* step0,
doublereal* x1, doublereal* step1, doublereal& s1,
OneDim& r, MultiJac& jac, int loglevel, bool writetitle);
//! Compute the weighted 2-norm of `step`.
doublereal norm2(const doublereal* x, const doublereal* step,
OneDim& r) const;
/**
* Find the solution to F(X) = 0 by damped Newton iteration. On
* entry, x0 contains an initial estimate of the solution. On
* successful return, x1 contains the converged solution.
*/
int solve(doublereal* x0, doublereal* x1, OneDim& r, MultiJac& jac,
int loglevel);
@ -56,22 +73,22 @@ public:
/// Change the problem size.
void resize(size_t points);
protected:
//! Get a pointer to an array of length m_n for temporary work space.
doublereal* getWorkArray();
//! Release a work array by pushing its pointer onto the stack of
//! available arrays.
void releaseWorkArray(doublereal* work);
std::vector<doublereal*> m_workarrays;
int m_maxAge;
size_t m_nv, m_np, m_n;
doublereal m_elapsed;
private:
char m_buf[100];
};
}
#endif

69
include/cantera/oneD/OneDim.h
View file

@ -20,19 +20,14 @@ class Func1;
*/
class OneDim
{
public:
// Default constructor.
OneDim();
// Constructor.
//! Construct a OneDim container for the domains in the list *domains*.
OneDim(std::vector<Domain1D*> domains);
/// Destructor.
virtual ~OneDim();
/// Add a domain.
/// Add a domain. Domains are added left-to-right.
void addDomain(Domain1D* d);
//! Return a reference to the Jacobian evaluator.
@ -116,7 +111,11 @@ public:
return m_bw;
}
/// Initialize.
/*!
* Initialize all domains. On the first call, this methods calls the init
* method of each domain, proceeding from left to right. Subsequent calls
* do nothing.
*/
void init();
/// Total number of points.
@ -125,8 +124,9 @@ public:
}
/**
* Steady-state max norm of the residual evaluated using solution x.
* On return, array r contains the steady-state residual values.
* Steady-state max norm (infinity norm) of the residual evaluated using
* solution x. On return, array r contains the steady-state residual
* values. Used only for diagnostic output.
*/
doublereal ssnorm(doublereal* x, doublereal* r);
@ -135,7 +135,7 @@ public:
return m_rdt;
}
/// Prepare for time stepping beginning with solution x.
//! Prepare for time stepping beginning with solution *x* and timestep *dt*.
void initTimeInteg(doublereal dt, doublereal* x);
/// True if transient mode.
@ -148,14 +148,14 @@ public:
return (m_rdt == 0.0);
}
/**
* Set steady mode. After invoking this method, subsequent
* calls to solve() will solve the steady-state problem.
/*!
* Prepare to solve the steady-state problem. After invoking this method,
* subsequent calls to solve() will solve the steady-state problem. Sets
* the reciprocal of the time step to zero, and, if it was previously non-
* zero, signals that a new Jacobian will be needed.
*/
void setSteadyMode();
/**
* Evaluate the multi-domain residual function
*
@ -170,18 +170,31 @@ public:
void eval(size_t j, double* x, double* r, doublereal rdt=-1.0,
int count = 1);
/// Pointer to the domain global point i belongs to.
//! Return a pointer to the domain global point *i* belongs to.
/*!
* The domains are scanned right-to-left, and the first one with starting
* location less or equal to i is returned.
*/
Domain1D* pointDomain(size_t i);
//! Call after one or more grids has been refined.
void resize();
//doublereal solveTime() { return m_solve_time; }
//void setTransientMask();
vector_int& transientMask() {
return m_mask;
}
/*!
* Take time steps using Backward Euler.
*
* @param nsteps number of steps
* @param dt initial step size
* @param x current solution vector
* @param r solution vector after time stepping
* @param loglevel controls amount of printed diagnostics
* @returns size of last timestep taken
*/
double timeStep(int nsteps, double dt, double* x,
double* r, int loglevel);
@ -214,6 +227,18 @@ public:
m_ts_jac_age = m_ss_jac_age;
}
}
/**
* Save statistics on function and Jacobian evaluation, and reset the
* counters. Statistics are saved only if the number of Jacobian
* evaluations is greater than zero. The statistics saved are:
*
* - number of grid points
* - number of Jacobian evaluations
* - CPU time spent evaluating Jacobians
* - number of non-Jacobian function evaluations
* - CPU time spent evaluating functions
*/
void saveStats();
//! Set a function that will be called every time #eval is called.
@ -224,7 +249,6 @@ public:
}
protected:
void evalSSJacobian(doublereal* x, doublereal* xnew);
doublereal m_tmin; // minimum timestep size
@ -259,7 +283,6 @@ protected:
Func1* m_interrupt;
private:
// statistics
int m_nevals;
doublereal m_evaltime;
@ -268,12 +291,8 @@ private:
vector_fp m_jacElapsed;
vector_int m_funcEvals;
vector_fp m_funcElapsed;
};
}
#endif

77
include/cantera/oneD/Sim1D.h
View file

@ -12,25 +12,21 @@ namespace Cantera
{
/**
* One-dimensional simulations. Class Sim1D extends class OneDim
* by storing the solution vector, and by adding a hybrid
* Newton/time-stepping solver.
* One-dimensional simulations. Class Sim1D extends class OneDim by storing
* the solution vector, and by adding a hybrid Newton/time-stepping solver.
*/
class Sim1D : public OneDim
{
public:
//! Default constructor.
/*!
* This constructor is provided to make
* the class default-constructible, but is not meant to be
* used in most applications. Use the next constructor
* This constructor is provided to make the class default-constructible,
* but is not meant to be used in most applications. Use the next
* constructor
*/
Sim1D();
/**
* Standard constructor.
* @param domains A vector of pointers to the domains to be linked together.
@ -52,15 +48,41 @@ public:
void setInitialGuess(const std::string& component, vector_fp& locs,
vector_fp& vals);
/// Set one entry in the solution vector.
/**
* Set a single value in the solution vector.
* @param dom domain number, beginning with 0 for the leftmost domain.
* @param comp component number
* @param localPoint grid point within the domain, beginning with 0 for
* the leftmost grid point in the domain.
* @param value the value.
*/
void setValue(size_t dom, size_t comp, size_t localPoint, doublereal value);
/// Get one entry in the solution vector.
/**
* Get one entry in the solution vector.
* @param dom domain number, beginning with 0 for the leftmost domain.
* @param comp component number
* @param localPoint grid point within the domain, beginning with 0 for
* the leftmost grid point in the domain.
*/
doublereal value(size_t dom, size_t comp, size_t localPoint) const;
doublereal workValue(size_t dom, size_t comp, size_t localPoint) const;
/// Specify a profile for one component of one domain.
/**
* Specify a profile for one component of one domain.
* @param dom domain number, beginning with 0 for the leftmost domain.
* @param comp component number
* @param pos A vector of relative positions, beginning with 0.0 at the
* left of the domain, and ending with 1.0 at the right of the domain.
* @param values A vector of values corresponding to the relative position
* locations.
*
* Note that the vector pos and values can have lengths different than the
* number of grid points, but their lengths must be equal. The values at
* the grid points will be linearly interpolated based on the (pos,
* values) specification.
*/
void setProfile(size_t dom, size_t comp, const vector_fp& pos,
const vector_fp& values);
@ -96,10 +118,16 @@ public:
/// Refine the grid in all domains.
int refine(int loglevel=0);
//! Add node for fixed temperature point of freely propagating flame
int setFixedTemperature(doublereal t);
void setAdiabaticFlame(void);
/// Set the criteria for grid refinement.
/**
* Set grid refinement criteria. If dom >= 0, then the settings
* apply only to the specified domain. If dom < 0, the settings
* are applied to each domain. @see Refiner::setCriteria.
*/
void setRefineCriteria(int dom = -1, doublereal ratio = 10.0,
doublereal slope = 0.8, doublereal curve = 0.8, doublereal prune = -0.1);
void setMaxGridPoints(int dom = -1, int npoints = 300);
@ -112,7 +140,9 @@ public:
*/
void setGridMin(int dom, double gridmin);
//! Initialize the solution with a previously-saved solution.
void restore(const std::string& fname, const std::string& id, int loglevel=2);
void getInitialSoln();
void setSolution(const doublereal* soln) {
@ -128,17 +158,20 @@ public:
void evalSSJacobian();
protected:
//! the solution vector
vector_fp m_x;
vector_fp m_x; // the solution vector
vector_fp m_xnew; // a work array used to hold the residual
// or the new solution
doublereal m_tstep; // timestep
vector_int m_steps; // array of number of steps to take before
// re-attempting the steady-state solution
//! a work array used to hold the residual or the new solution
vector_fp m_xnew;
//! timestep
doublereal m_tstep;
//! array of number of steps to take before re-attempting the steady-state
//! solution
vector_int m_steps;
private:
/// Calls method _finalize in each domain.
void finalize();
@ -146,11 +179,7 @@ private:
* @return 0 if successful, -1 on failure
*/
int newtonSolve(int loglevel);
};
}
#endif

159
include/cantera/oneD/StFlow.h
View file

@ -1,6 +1,5 @@
/**
* @file StFlow.h
*
*/
// Copyright 2001 California Institute of Technology
@ -34,36 +33,26 @@ const int c_Mixav_Transport = 0;
const int c_Multi_Transport = 1;
const int c_Soret = 2;
//-----------------------------------------------------------
// Class StFlow
//-----------------------------------------------------------
/**
* This class represents 1D flow domains that satisfy the
* one-dimensional similarity solution for chemically-reacting,
* axisymmetric, flows.
* This class represents 1D flow domains that satisfy the one-dimensional
* similarity solution for chemically-reacting, axisymmetric, flows.
*/
class StFlow : public Domain1D
{
public:
//--------------------------------
// construction and destruction
//--------------------------------
/// Constructor. Create a new flow domain.
/// @param ph Object representing the gas phase. This object
/// will be used to evaluate all thermodynamic, kinetic, and transport
/// properties.
/// @param nsp Number of species.
//! Create a new flow domain.
//! @param ph Object representing the gas phase. This object will be used
//! to evaluate all thermodynamic, kinetic, and transport properties.
//! @param nsp Number of species.
//! @param points Initial number of grid points
StFlow(IdealGasPhase* ph = 0, size_t nsp = 1, size_t points = 1);
/**
* @name Problem Specification
*/
//@{
//! @name Problem Specification
//! @{
virtual void setupGrid(size_t n, const doublereal* z);
@ -85,31 +74,31 @@ public:
m_thermo = &th;
}
/// Set the kinetics manager. The kinetics manager must
//! Set the kinetics manager. The kinetics manager must
void setKinetics(Kinetics& kin) {
m_kin = &kin;
}
/// set the transport manager
//! set the transport manager
void setTransport(Transport& trans, bool withSoret = false);
void enableSoret(bool withSoret);
bool withSoret() const {
return m_do_soret;
}
/// Set the pressure. Since the flow equations are for the limit of
/// small Mach number, the pressure is very nearly constant
/// throughout the flow.
//! Set the pressure. Since the flow equations are for the limit of
//! small Mach number, the pressure is very nearly constant
//! throughout the flow.
void setPressure(doublereal p) {
m_press = p;
}
/// The current pressure [Pa].
//! The current pressure [Pa].
doublereal pressure() const {
return m_press;
}
/// @todo remove? may be unused
//! @todo remove? may be unused
virtual void setState(size_t point, const doublereal* state,
doublereal* x) {
setTemperature(point, state[2]);
@ -118,7 +107,7 @@ public:
}
}
/// Write the initial solution estimate into array x.
//! Write the initial solution estimate into array x.
virtual void _getInitialSoln(doublereal* x) {
for (size_t j = 0; j < m_points; j++) {
T(x,j) = m_thermo->temperature();
@ -128,64 +117,58 @@ public:
virtual void _finalize(const doublereal* x);
/// Sometimes it is desired to carry out the simulation
/// using a specified temperature profile, rather than
/// computing it by solving the energy equation. This
/// method specifies this profile.
//! Sometimes it is desired to carry out the simulation using a specified
//! temperature profile, rather than computing it by solving the energy
//! equation. This method specifies this profile.
void setFixedTempProfile(vector_fp& zfixed, vector_fp& tfixed) {
m_zfix = zfixed;
m_tfix = tfixed;
}
/**
* Set the temperature fixed point at grid point j, and
* disable the energy equation so that the solution will be
* held to this value.
/*!
* Set the temperature fixed point at grid point j, and disable the energy
* equation so that the solution will be held to this value.
*/
void setTemperature(size_t j, doublereal t) {
m_fixedtemp[j] = t;
m_do_energy[j] = false;
}
/**
* Set the mass fraction fixed point for species k at grid
* point j, and disable the species equation so that the
* solution will be held to this value.
* note: in practice, the species are hardly ever held fixed.
/*!
* Set the mass fraction fixed point for species k at grid point j, and
* disable the species equation so that the solution will be held to this
* value. Note: in practice, the species are hardly ever held fixed.
*/
void setMassFraction(size_t j, size_t k, doublereal y) {
m_fixedy(k,j) = y;
m_do_species[k] = true; // false;
}
/// The fixed temperature value at point j.
//! The fixed temperature value at point j.
doublereal T_fixed(size_t j) const {
return m_fixedtemp[j];
}
/// The fixed mass fraction value of species k at point j.
//! The fixed mass fraction value of species k at point j.
doublereal Y_fixed(size_t k, size_t j) const {
return m_fixedy(k,j);
}
// @}
virtual std::string componentName(size_t n) const;
size_t componentIndex(const std::string& name) const;
//! Print the solution.
virtual void showSolution(const doublereal* x);
//! Save the current solution for this domain into an XML_Node
/*!
*
* @param o XML_Node to save the solution to.
* @param sol Current value of the solution vector.
* The object will pick out which part of the solution
* vector pertains to this object.
* @param sol Current value of the solution vector. The object will pick
* out which part of the solution vector pertains to this
* object.
*/
virtual XML_Node& save(XML_Node& o, const doublereal* const sol);
@ -256,13 +239,18 @@ public:
void integrateChem(doublereal* x,doublereal dt);
//! Change the grid size. Called after grid refinement.
void resize(size_t components, size_t points);
virtual void setFixedPoint(int j0, doublereal t0) {}
void setJac(MultiJac* jac);
//! Set the gas object state to be consistent with the solution at point j.
void setGas(const doublereal* x, size_t j);
//! Set the gas state to be consistent with the solution at the midpoint
//! between j and j + 1.
void setGasAtMidpoint(const doublereal* x, size_t j);
doublereal density(size_t j) const {
@ -276,6 +264,13 @@ public:
m_dovisc = dovisc;
}
/*!
* Evaluate the residual function for axisymmetric stagnation flow. If
* jpt is less than zero, the residual function is evaluated at all grid
* points. If jpt >= 0, then the residual function is only evaluated at
* grid points jpt-1, jpt, and jpt+1. This option is used to efficiently
* evaluate the Jacobian numerically.
*/
virtual void eval(size_t j, doublereal* x, doublereal* r,
integer* mask, doublereal rdt);
@ -289,7 +284,6 @@ public:
integer* diag, doublereal rdt) = 0;
protected:
doublereal component(const doublereal* x, size_t i, size_t j) const {
return x[index(i,j)];
}
@ -306,15 +300,15 @@ protected:
return m_wdot(k,j);
}
/// write the net production rates at point j into array m_wdot
//! Write the net production rates at point `j` into array `m_wdot`
void getWdot(doublereal* x, size_t j) {
setGas(x,j);
m_kin->getNetProductionRates(&m_wdot(0,j));
}
/**
* update the thermodynamic properties from point
* j0 to point j1 (inclusive), based on solution x.
* Update the thermodynamic properties from point j0 to point j1
* (inclusive), based on solution x.
*/
void updateThermo(const doublereal* x, size_t j0, size_t j1) {
for (size_t j = j0; j <= j1; j++) {
@ -325,7 +319,6 @@ protected:
}
}
//--------------------------------
// central-differenced derivatives
//--------------------------------
@ -338,10 +331,8 @@ protected:
}
//--------------------------------
// solution components
//--------------------------------
//! @name Solution components
//! @{
doublereal T(const doublereal* x, size_t j) const {
return x[index(c_offset_T, j)];
@ -391,11 +382,11 @@ protected:
doublereal flux(size_t k, size_t j) const {
return m_flux(k, j);
}
//! @}
// convective spatial derivatives. These use upwind
// differencing, assuming u(z) is negative
//! @name convective spatial derivatives.
//! These use upwind differencing, assuming u(z) is negative
//! @{
doublereal dVdz(const doublereal* x, size_t j) const {
size_t jloc = (u(x,j) > 0.0 ? j : j + 1);
return (V(x,jloc) - V(x,jloc-1))/m_dz[jloc-1];
@ -410,6 +401,7 @@ protected:
size_t jloc = (u(x,j) > 0.0 ? j : j + 1);
return (T(x,jloc) - T(x,jloc-1))/m_dz[jloc-1];
}
//! @}
doublereal shear(const doublereal* x, size_t j) const {
doublereal c1 = m_visc[j-1]*(V(x,j) - V(x,j-1));
@ -427,13 +419,11 @@ protected:
return m*m_nsp*m_nsp + m_nsp*j + k;
}
//! Update the diffusive mass fluxes.
void updateDiffFluxes(const doublereal* x, size_t j0, size_t j1);
//---------------------------------------------------------
//
// member data
//
//---------------------------------------------------------
// inlet
@ -500,14 +490,15 @@ protected:
vector_fp m_tfix;
bool m_dovisc;
//! Update the transport properties at grid points in the range from `j0`
//! to `j1`, based on solution `x`.
void updateTransport(doublereal* x, size_t j0, size_t j1);
private:
vector_fp m_ybar;
};
/**
* A class for axisymmetric stagnation flows.
*/
@ -554,26 +545,14 @@ public:
}
};
/*
class OneDFlow : public StFlow {
public:
OneDFlow(igthermo_t* ph = 0, int nsp = 1, int points = 1) :
StFlow(ph, nsp, points) {
}
virtual ~OneDFlow() {}
virtual void eval(int j, doublereal* x, doublereal* r,
integer* mask, doublereal rdt);
virtual std::string flowType() { return "OneDFlow"; }
doublereal mdot(doublereal* x, int j) {
return x[index(c_offset_L,j)];
}
private:
void updateTransport(doublereal* x,int j0, int j1);
};
*/
/**
* Import a previous solution to use as an initial estimate. The
* previous solution may have been computed using a different
* reaction mechanism. Species in the old and new mechanisms are
* matched by name, and any species in the new mechanism that were
* not in the old one are set to zero. The new solution is created
* with the same number of grid points as in the old solution.
*/
void importSolution(size_t points, doublereal* oldSoln, IdealGasPhase& oldmech,
size_t size_new, doublereal* newSoln, IdealGasPhase& newmech);

16
include/cantera/oneD/refine.h
View file

@ -8,12 +8,25 @@ namespace Cantera
class Domain1D;
//! Refine Domain1D grids so that profiles satisfy adaptation tolerances
class Refiner
{
public:
Refiner(Domain1D& domain);
virtual ~Refiner() {}
//! Set grid refinement criteria
/*!
* @param ratio Maximum ratio between grid spacing at adjacent intervals.
* E.g. `(x[j+1] - x[j]) / (x[j] - x[j-1]) < ratio`
* @param slope Maximum fractional change in the value of each solution
* component between adjacent grid points
* @param curve Maximum fractional change in the derivative of each
* solution component between adjacent grid points.
* @param prune Threshold for removing unnecessary grid points. `prune`
* should be smaller than both `slope` and `curve`. Set `prune <= 0`
* to disable pruning.
*/
void setCriteria(doublereal ratio = 10.0,
doublereal slope = 0.8,
doublereal curve = 0.8,
@ -26,6 +39,8 @@ public:
void setActive(int comp, bool state = true) {
m_active[comp] = state;
}
//! Set the maximum number of points allowed in the domain
void setMaxPoints(int npmax) {
m_npmax = npmax;
}
@ -54,6 +69,7 @@ public:
return (m_keep[j] != -1);
}
double value(const double* x, size_t i, size_t j);
double maxRatio() {
return m_ratio;
}

10
src/oneD/Domain1D.cpp
View file

@ -1,6 +1,5 @@
/**
* @file Domain1D.cpp
*
*/
#include "cantera/oneD/Domain1D.h"
@ -165,7 +164,6 @@ void Domain1D::restore(const XML_Node& dom, doublereal* soln, int loglevel)
}
}
// called to set up initial grid, and after grid refinement
void Domain1D::setupGrid(size_t n, const doublereal* z)
{
if (n > 1) {
@ -176,17 +174,12 @@ void Domain1D::setupGrid(size_t n, const doublereal* z)
}
}
void drawline()
{
writelog("\n-------------------------------------"
"------------------------------------------");
}
/**
* Print the solution.
*/
void Domain1D::showSolution(const doublereal* x)
{
size_t nn = m_nv/5;
@ -235,8 +228,6 @@ void Domain1D::showSolution(const doublereal* x)
writelog("\n");
}
// initial solution
void Domain1D::_getInitialSoln(doublereal* x)
{
for (size_t j = 0; j < m_points; j++) {
@ -253,5 +244,4 @@ doublereal Domain1D::initialValue(size_t n, size_t j)
return 0.0;
}
} // namespace

10
src/oneD/MultiJac.cpp
View file

@ -1,7 +1,5 @@
/**
* @file MultiJac.cpp
*
* Implementation file for class MultiJac
* @file MultiJac.cpp Implementation file for class MultiJac
*/
/*
@ -47,10 +45,6 @@ void MultiJac::incrementDiagonal(int j, doublereal d)
value(j,j) = m_ssdiag[j];
}
/**
* Evaluate the Jacobian at x0. The array of residual values at x0
* is supplied as an input.
*/
void MultiJac::eval(doublereal* x0, doublereal* resid0, doublereal rdt)
{
m_nevals++;
@ -99,5 +93,3 @@ void MultiJac::eval(doublereal* x0, doublereal* resid0, doublereal rdt)
}
} // namespace
// $Log: MultiJac.cpp,v

55
src/oneD/MultiNewton.cpp
View file

@ -1,7 +1,5 @@
/**
* @file MultiNewton.cpp
*
* Damped Newton solver for 1D multi-domain problems
* @file MultiNewton.cpp Damped Newton solver for 1D multi-domain problems
*/
/*
@ -9,7 +7,6 @@
*/
#include <vector>
using namespace std;
#include "cantera/oneD/MultiNewton.h"
@ -40,7 +37,6 @@ public:
}
};
/**
* Return a damping coefficient that keeps the solution after taking one
* Newton step between specified lower and upper bounds. This function only
@ -49,7 +45,6 @@ public:
doublereal bound_step(const doublereal* x, const doublereal* step,
Domain1D& r, int loglevel)
{
char buf[100];
size_t np = r.nPoints();
size_t nv = r.nComponents();
@ -102,7 +97,6 @@ doublereal bound_step(const doublereal* x, const doublereal* step,
return fbound;
}
/**
* This function computes the square of a weighted norm of a step
* vector for one domain.
@ -125,7 +119,6 @@ doublereal bound_step(const doublereal* x, const doublereal* step,
* solution component n in the domain. The second term,
* \f$\epsilon_{a,n}\f$, is the absolute error tolerance for component
* n.
*
*/
doublereal norm_square(const doublereal* x,
const doublereal* step, Domain1D& r)
@ -165,13 +158,10 @@ const string dashedline =
const doublereal DampFactor = sqrt(2.0);
const size_t NDAMP = 7;
//-----------------------------------------------------------
// MultiNewton methods
//-----------------------------------------------------------
MultiNewton::MultiNewton(int sz)
: m_maxAge(5)
{
@ -186,9 +176,6 @@ MultiNewton::~MultiNewton()
}
}
/**
* Prepare for a new solution vector length.
*/
void MultiNewton::resize(size_t sz)
{
m_n = sz;
@ -198,10 +185,6 @@ void MultiNewton::resize(size_t sz)
m_workarrays.clear();
}
/**
* Compute the weighted 2-norm of 'step'.
*/
doublereal MultiNewton::norm2(const doublereal* x,
const doublereal* step, OneDim& r) const
{
@ -216,11 +199,6 @@ doublereal MultiNewton::norm2(const doublereal* x,
return sqrt(sum);
}
/**
* Compute the undamped Newton step. The residual function is
* evaluated at x, but the Jacobian is not recomputed.
*/
void MultiNewton::step(doublereal* x, doublereal* step,
OneDim& r, MultiJac& jac, int loglevel)
{
@ -281,12 +259,6 @@ void MultiNewton::step(doublereal* x, doublereal* step,
#endif
}
/**
* Return the factor by which the undamped Newton step 'step0'
* must be multiplied in order to keep all solution components in
* all domains between their specified lower and upper bounds.
*/
doublereal MultiNewton::boundStep(const doublereal* x0,
const doublereal* step0, const OneDim& r, int loglevel)
{
@ -299,20 +271,10 @@ doublereal MultiNewton::boundStep(const doublereal* x0,
return fbound;
}
/**
* On entry, step0 must contain an undamped Newton step for the
* solution x0. This method attempts to find a damping coefficient
* such that the next undamped step would have a norm smaller than
* that of step0. If successful, the new solution after taking the
* damped step is returned in x1, and the undamped step at x1 is
* returned in step1.
*/
int MultiNewton::dampStep(const doublereal* x0, const doublereal* step0,
doublereal* x1, doublereal* step1, doublereal& s1,
OneDim& r, MultiJac& jac, int loglevel, bool writetitle)
{
// write header
if (loglevel > 0 && writetitle) {
writelog("\n\nDamped Newton iteration:\n");
@ -404,12 +366,6 @@ int MultiNewton::dampStep(const doublereal* x0, const doublereal* step0,
}
}
/**
* Find the solution to F(X) = 0 by damped Newton iteration. On
* entry, x0 contains an initial estimate of the solution. On
* successful return, x1 contains the converged solution.
*/
int MultiNewton::solve(doublereal* x0, doublereal* x1,
OneDim& r, MultiJac& jac, int loglevel)
{
@ -510,11 +466,6 @@ int MultiNewton::solve(doublereal* x0, doublereal* x1,
return m;
}
/**
* Get a pointer to an array of length m_n for temporary work
* space.
*/
doublereal* MultiNewton::getWorkArray()
{
doublereal* w = 0;
@ -528,10 +479,6 @@ doublereal* MultiNewton::getWorkArray()
return w;
}
/**
* Release a work array by pushing its pointer onto the stack of
* available arrays.
*/
void MultiNewton::releaseWorkArray(doublereal* work)
{
m_workarrays.push_back(work);

84
src/oneD/OneDim.cpp
View file

@ -13,9 +13,6 @@ using namespace std;
namespace Cantera
{
/**
* Default constructor. Create an empty object.
*/
OneDim::OneDim()
: m_tmin(1.0e-16), m_tmax(10.0), m_tfactor(0.5),
m_jac(0), m_newt(0),
@ -30,11 +27,6 @@ OneDim::OneDim()
//m_solve_time = 0.0;
}
/**
* Construct a OneDim container for the domains pointed at by the
* input vector of pointers.
*/
OneDim::OneDim(vector<Domain1D*> domains) :
m_tmin(1.0e-16), m_tmax(10.0), m_tfactor(0.5),
m_jac(0), m_newt(0),
@ -57,7 +49,6 @@ OneDim::OneDim(vector<Domain1D*> domains) :
resize();
}
size_t OneDim::domainIndex(const std::string& name)
{
for (size_t n = 0; n < m_nd; n++) {
@ -69,13 +60,8 @@ size_t OneDim::domainIndex(const std::string& name)
return npos;
}
/**
* Domains are added left-to-right.
*/
void OneDim::addDomain(Domain1D* d)
{
// if 'd' is not the first domain, link it to the last domain
// added (the rightmost one)
int n = static_cast<int>(m_dom.size());
@ -98,7 +84,6 @@ void OneDim::addDomain(Domain1D* d)
resize();
}
OneDim::~OneDim()
{
delete m_jac;
@ -114,7 +99,6 @@ MultiNewton& OneDim::newton()
return *m_newt;
}
//==============================================================================================================
void OneDim::writeStats(int printTime)
{
saveStats();
@ -134,20 +118,7 @@ void OneDim::writeStats(int printTime)
writelog(buf);
}
}
//==============================================================================================================
/**
* Save statistics on function and Jacobian evaluation, and reset
* the counters. Statistics are saved only if the number of
* Jacobian evaluations is greater than zero. The statistics saved
* are
*
* - number of grid points
* - number of Jacobian evaluations
* - CPU time spent evaluating Jacobians
* - number of non-Jacobian function evaluations
* - CPU time spent evaluating functions
*/
void OneDim::saveStats()
{
if (m_jac) {
@ -164,10 +135,6 @@ void OneDim::saveStats()
}
}
/**
* Call after one or more grids has been refined.
*/
void OneDim::resize()
{
m_bw = 0;
@ -232,7 +199,6 @@ void OneDim::resize()
}
}
int OneDim::solve(doublereal* x, doublereal* xnew, int loglevel)
{
if (!m_jac_ok) {
@ -254,15 +220,6 @@ void OneDim::evalSSJacobian(doublereal* x, doublereal* xnew)
m_rdt = rdt_save;
}
/**
* Return a pointer to the domain that contains component i of the
* global solution vector. The domains are scanned right-to-left,
* and the first one with starting location less or equal to i is
* returned.
*
* 8/26/02 changed '<' to '<=' DGG
*
*/
Domain1D* OneDim::pointDomain(size_t i)
{
Domain1D* d = right();
@ -275,11 +232,6 @@ Domain1D* OneDim::pointDomain(size_t i)
return 0;
}
/**
* Evaluate the multi-domain residual function, and return the
* result in array r.
*/
void OneDim::eval(size_t j, double* x, double* r, doublereal rdt, int count)
{
clock_t t0 = clock();
@ -312,11 +264,6 @@ void OneDim::eval(size_t j, double* x, double* r, doublereal rdt, int count)
}
}
/**
* The 'infinity' (maximum magnitude) norm of the steady-state
* residual. Used only for diagnostic output.
*/
doublereal OneDim::ssnorm(doublereal* x, doublereal* r)
{
eval(npos, x, r, 0.0, 0);
@ -327,10 +274,6 @@ doublereal OneDim::ssnorm(doublereal* x, doublereal* r)
return ss;
}
/**
* Prepare for time stepping with timestep dt.
*/
void OneDim::initTimeInteg(doublereal dt, doublereal* x)
{
doublereal rdt_old = m_rdt;
@ -351,12 +294,6 @@ void OneDim::initTimeInteg(doublereal dt, doublereal* x)
}
}
/**
* Prepare to solve the steady-state problem. Set the reciprocal
* of the time step to zero, and, if it was previously non-zero,
* signal that a new Jacobian will be needed.
*/
void OneDim::setSteadyMode()
{
m_rdt = 0.0;
@ -370,11 +307,6 @@ void OneDim::setSteadyMode()
}
}
/**
* Initialize all domains. On the first call, this methods calls
* the init method of each domain, proceeding from left to right.
* Subsequent calls do nothing.
*/
void OneDim::init()
{
if (!m_init) {
@ -387,10 +319,6 @@ void OneDim::init()
m_init = true;
}
/**
* Signal that the current Jacobian is no longer valid.
*/
void Domain1D::needJacUpdate()
{
if (m_container) {
@ -399,18 +327,9 @@ void Domain1D::needJacUpdate()
}
}
/**
* Take time steps using Backward Euler.
*
* nsteps -- number of steps
* dt -- initial step size
* loglevel -- controls amount of printed diagnostics
*/
doublereal OneDim::timeStep(int nsteps, doublereal dt, doublereal* x,
doublereal* r, int loglevel)
{
// set the Jacobian age parameter to the transient value
newton().setOptions(m_ts_jac_age);
@ -468,12 +387,10 @@ doublereal OneDim::timeStep(int nsteps, doublereal dt, doublereal* x,
return dt;
}
void OneDim::save(const std::string& fname, std::string id,
const std::string& desc, doublereal* sol,
int loglevel)
{
struct tm* newtime;
time_t aclock;
::time(&aclock); /* Get time in seconds */
@ -519,7 +436,6 @@ void OneDim::save(const std::string& fname, std::string id,
writelog("Solution saved to file "+fname+" as solution "+id+".\n", loglevel);
}
void Domain1D::setGrid(size_t n, const doublereal* z)
{
m_z.resize(n);

67
src/oneD/Sim1D.cpp
View file

@ -20,20 +20,18 @@ static void sim1D_drawline()
s += '\n';
writelog(s.c_str());
}
//====================================================================================================================
Sim1D::Sim1D() :
OneDim()
{
//writelog("Sim1D default constructor\n");
}
//====================================================================================================================
Sim1D::Sim1D(vector<Domain1D*>& domains) :
OneDim(domains)
{
// resize the internal solution vector and the wprk array,
// and perform domain-specific initialization of the
// solution vector.
// resize the internal solution vector and the wprk array, and perform
// domain-specific initialization of the solution vector.
m_x.resize(size(), 0.0);
m_xnew.resize(size(), 0.0);
@ -49,13 +47,10 @@ Sim1D::Sim1D(vector<Domain1D*>& domains) :
m_steps.push_back(2);
m_steps.push_back(5);
m_steps.push_back(10);
}
//====================================================================================================================
void Sim1D::setInitialGuess(const std::string& component, vector_fp& locs, vector_fp& vals)
{
for (size_t dom=0; dom<m_nd; dom++) {
Domain1D& d = domain(dom);
size_t ncomp = d.nComponents();
@ -67,15 +62,6 @@ void Sim1D::setInitialGuess(const std::string& component, vector_fp& locs, vecto
}
}
/**
* Set a single value in the solution vector.
* @param dom domain number, beginning with 0 for the leftmost domain.
* @param comp component number
* @param localPoint grid point within the domain, beginning with 0 for
* the leftmost grid point in the domain.
* @param value the value.
*/
void Sim1D::setValue(size_t dom, size_t comp, size_t localPoint, doublereal value)
{
size_t iloc = domain(dom).loc() + domain(dom).index(comp, localPoint);
@ -85,13 +71,6 @@ void Sim1D::setValue(size_t dom, size_t comp, size_t localPoint, doublereal val
m_x[iloc] = value;
}
/**
* @param dom domain number, beginning with 0 for the leftmost domain.
* @param comp component number
* @param localPoint grid point within the domain, beginning with 0 for
* the leftmost grid point in the domain.
*/
doublereal Sim1D::value(size_t dom, size_t comp, size_t localPoint) const
{
size_t iloc = domain(dom).loc() + domain(dom).index(comp, localPoint);
@ -110,21 +89,6 @@ doublereal Sim1D::workValue(size_t dom, size_t comp, size_t localPoint) const
return m_xnew[iloc];
}
/**
* @param dom domain number, beginning with 0 for the leftmost domain.
* @param comp component number
* @param pos A vector of relative positions, beginning with 0.0 at the
* left of the domain, and ending with 1.0 at the right of the domain.
* @param values A vector of values corresponding to the relative position
* locations.
*
* Note that the vector pos and values can have lengths
* different than the number of grid points, but their lengths
* must be equal. The values at the grid points will be
* linearly interpolated based on the (pos, values)
* specification.
*/
void Sim1D::setProfile(size_t dom, size_t comp,
const vector_fp& pos, const vector_fp& values)
{
@ -141,7 +105,6 @@ void Sim1D::setProfile(size_t dom, size_t comp,
}
}
void Sim1D::save(const std::string& fname, const std::string& id,
const std::string& desc, int loglevel)
{
@ -156,9 +119,6 @@ void Sim1D::saveResidual(const std::string& fname, const std::string& id,
OneDim::save(fname, id, desc, &res[0], loglevel);
}
/**
* Initialize the solution with a previously-saved solution.
*/
void Sim1D::restore(const std::string& fname, const std::string& id,
int loglevel)
{
@ -203,7 +163,6 @@ void Sim1D::restore(const std::string& fname, const std::string& id,
finalize();
}
void Sim1D::setFlatProfile(size_t dom, size_t comp, doublereal v)
{
size_t np = domain(dom).nPoints();
@ -213,7 +172,6 @@ void Sim1D::setFlatProfile(size_t dom, size_t comp, doublereal v)
}
}
void Sim1D::showSolution(ostream& s)
{
for (size_t n = 0; n < m_nd; n++) {
@ -248,7 +206,6 @@ void Sim1D::finalize()
}
}
void Sim1D::setTimeStep(doublereal stepsize, size_t n, integer* tsteps)
{
m_tstep = stepsize;
@ -258,7 +215,6 @@ void Sim1D::setTimeStep(doublereal stepsize, size_t n, integer* tsteps)
}
}
int Sim1D::newtonSolve(int loglevel)
{
int m = OneDim::solve(DATA_PTR(m_x), DATA_PTR(m_xnew), loglevel);
@ -273,7 +229,6 @@ int Sim1D::newtonSolve(int loglevel)
}
}
void Sim1D::solve(int loglevel, bool refine_grid)
{
int new_points = 1;
@ -388,10 +343,6 @@ void Sim1D::solve(int loglevel, bool refine_grid)
}
}
/**
* Refine the grid in all domains.
*/
int Sim1D::refine(int loglevel)
{
int ianalyze, np = 0;
@ -485,10 +436,6 @@ int Sim1D::refine(int loglevel)
return np;
}
/**
* Add node for fixed temperature point of freely propagating flame
*/
int Sim1D::setFixedTemperature(doublereal t)
{
int np = 0;
@ -598,11 +545,6 @@ void Sim1D::setAdiabaticFlame(void)
}
}
/**
* Set grid refinement criteria. If dom >= 0, then the settings
* apply only to the specified domain. If dom < 0, the settings
* are applied to each domain. @see Refiner::setCriteria.
*/
void Sim1D::setRefineCriteria(int dom, doublereal ratio,
doublereal slope, doublereal curve, doublereal prune)
{
@ -630,7 +572,6 @@ void Sim1D::setGridMin(int dom, double gridmin)
}
}
void Sim1D::setMaxGridPoints(int dom, int npoints)
{
if (dom >= 0) {

69
src/oneD/StFlow.cpp
View file

@ -18,22 +18,10 @@ using namespace std;
namespace Cantera
{
//------------------- importSolution ------------------------
/**
* Import a previous solution to use as an initial estimate. The
* previous solution may have been computed using a different
* reaction mechanism. Species in the old and new mechanisms are
* matched by name, and any species in the new mechanism that were
* not in the old one are set to zero. The new solution is created
* with the same number of grid points as in the old solution.
*/
void importSolution(size_t points,
doublereal* oldSoln, IdealGasPhase& oldmech,
size_t size_new, doublereal* newSoln, IdealGasPhase& newmech)
{
// Number of components in old and new solutions
size_t nv_old = oldmech.nSpecies() + 4;
size_t nv_new = newmech.nSpecies() + 4;
@ -81,7 +69,6 @@ void importSolution(size_t points,
}
}
static void st_drawline()
{
writelog("\n-------------------------------------"
@ -193,10 +180,6 @@ StFlow::StFlow(IdealGasPhase* ph, size_t nsp, size_t points) :
setID("stagnation flow");
}
/**
* Change the grid size. Called after grid refinement.
*/
void StFlow::resize(size_t ncomponents, size_t points)
{
Domain1D::resize(ncomponents, points);
@ -225,7 +208,6 @@ void StFlow::resize(size_t ncomponents, size_t points)
m_z.resize(m_points);
}
void StFlow::setupGrid(size_t n, const doublereal* z)
{
resize(m_nv, n);
@ -238,10 +220,6 @@ void StFlow::setupGrid(size_t n, const doublereal* z)
}
}
/**
* Install a transport manager.
*/
void StFlow::setTransport(Transport& trans, bool withSoret)
{
m_trans = &trans;
@ -275,11 +253,6 @@ void StFlow::enableSoret(bool withSoret)
}
}
/**
* Set the gas object state to be consistent with the solution at
* point j.
*/
void StFlow::setGas(const doublereal* x, size_t j)
{
m_thermo->setTemperature(T(x,j));
@ -288,11 +261,6 @@ void StFlow::setGas(const doublereal* x, size_t j)
m_thermo->setPressure(m_press);
}
/**
* Set the gas state to be consistent with the solution at the
* midpoint between j and j + 1.
*/
void StFlow::setGasAtMidpoint(const doublereal* x, size_t j)
{
m_thermo->setTemperature(0.5*(T(x,j)+T(x,j+1)));
@ -305,7 +273,6 @@ void StFlow::setGasAtMidpoint(const doublereal* x, size_t j)
m_thermo->setPressure(m_press);
}
void StFlow::_finalize(const doublereal* x)
{
size_t k, j;
@ -353,23 +320,9 @@ void StFlow::_finalize(const doublereal* x)
}
}
//------------------------------------------------------
/**
* Evaluate the residual function for axisymmetric stagnation
* flow. If jpt is less than zero, the residual function is
* evaluated at all grid points. If jpt >= 0, then the residual
* function is only evaluated at grid points jpt-1, jpt, and
* jpt+1. This option is used to efficiently evaluate the
* Jacobian numerically.
*
*/
void StFlow::eval(size_t jg, doublereal* xg,
doublereal* rg, integer* diagg, doublereal rdt)
{
// if evaluating a Jacobian, and the global point is outside
// the domain of influence for this domain, then skip
// evaluating the residual
@ -555,10 +508,6 @@ void StFlow::eval(size_t jg, doublereal* xg,
}
}
/**
* Update the transport properties at grid points in the range
* from j0 to j1, based on solution x.
*/
void StFlow::updateTransport(doublereal* x, size_t j0, size_t j1)
{
if (m_transport_option == c_Mixav_Transport) {
@ -602,13 +551,6 @@ void StFlow::updateTransport(doublereal* x, size_t j0, size_t j1)
}
}
/**
* Print the solution.
*/
void StFlow::showSolution(const doublereal* x)
{
size_t nn = m_nv/5;
@ -660,10 +602,6 @@ void StFlow::showSolution(const doublereal* x)
writelog("\n");
}
/**
* Update the diffusive mass fluxes.
*/
void StFlow::updateDiffFluxes(const doublereal* x, size_t j0, size_t j1)
{
size_t j, k, m;
@ -706,7 +644,6 @@ void StFlow::updateDiffFluxes(const doublereal* x, size_t j0, size_t j1)
}
}
string StFlow::componentName(size_t n) const
{
switch (n) {
@ -727,11 +664,8 @@ string StFlow::componentName(size_t n) const
}
}
size_t StFlow::componentIndex(const std::string& name) const
{
if (name=="u") {
return 0;
} else if (name=="V") {
@ -751,7 +685,6 @@ size_t StFlow::componentIndex(const std::string& name) const
return npos;
}
void StFlow::restore(const XML_Node& dom, doublereal* soln, int loglevel)
{
Domain1D::restore(dom, soln, loglevel);
@ -922,7 +855,6 @@ void StFlow::restore(const XML_Node& dom, doublereal* soln, int loglevel)
}
}
XML_Node& StFlow::save(XML_Node& o, const doublereal* const sol)
{
size_t k;
@ -1061,7 +993,6 @@ void FreeFlame::evalRightBoundary(doublereal* x, doublereal* rsd,
diag[index(4,j)] = 0;
}
void FreeFlame::evalContinuity(size_t j, doublereal* x, doublereal* rsd,
integer* diag, doublereal rdt)
{

40
src/oneD/boundaries1D.cpp
View file

@ -23,7 +23,6 @@ Bdry1D::Bdry1D() : Domain1D(1, 1, 0.0),
m_type = cConnectorType;
}
void Bdry1D::
_init(size_t n)
{
@ -73,16 +72,10 @@ _init(size_t n)
}
}
//----------------------------------------------------------
//
// Inlet1D methods
//
//----------------------------------------------------------
void Inlet1D::
setMoleFractions(const std::string& xin)
{
@ -121,7 +114,6 @@ componentName(size_t n) const
void Inlet1D::
init()
{
_init(2);
// set bounds (mdot, T)
@ -158,7 +150,6 @@ init()
}
}
void Inlet1D::
eval(size_t jg, doublereal* xg, doublereal* rg,
integer* diagg, doublereal rdt)
@ -234,8 +225,6 @@ eval(size_t jg, doublereal* xg, doublereal* rg,
rb[4+k] += x[0]*(m_yin[k]);
}
}
}
XML_Node& Inlet1D::
@ -275,8 +264,6 @@ restore(const XML_Node& dom, doublereal* soln, int loglevel)
resize(2,1);
}
//--------------------------------------------------
// Empty1D
//--------------------------------------------------
@ -339,8 +326,6 @@ restore(const XML_Node& dom, doublereal* soln, int loglevel)
resize(1,1);
}
//--------------------------------------------------
// Symm1D
//--------------------------------------------------
@ -413,7 +398,6 @@ eval(size_t jg, doublereal* xg, doublereal* rg,
}
}
XML_Node& Symm1D::
save(XML_Node& o, const doublereal* const soln)
{
@ -429,7 +413,6 @@ restore(const XML_Node& dom, doublereal* soln, int loglevel)
resize(1,1);
}
//--------------------------------------------------
// Outlet1D
//--------------------------------------------------
@ -466,7 +449,6 @@ init()
}
}
void Outlet1D::
eval(size_t jg, doublereal* xg, doublereal* rg,
integer* diagg, doublereal rdt)
@ -522,7 +504,6 @@ eval(size_t jg, doublereal* xg, doublereal* rg,
}
}
XML_Node& Outlet1D::
save(XML_Node& o, const doublereal* const soln)
{
@ -538,14 +519,10 @@ restore(const XML_Node& dom, doublereal* soln, int loglevel)
resize(1,1);
}
//--------------------------------------------------
// OutletRes1D
//--------------------------------------------------
void OutletRes1D::
setMoleFractions(const std::string& xres)
{
@ -609,7 +586,6 @@ init()
}
}
void OutletRes1D::
eval(size_t jg, doublereal* xg, doublereal* rg,
integer* diagg, doublereal rdt)
@ -670,7 +646,6 @@ eval(size_t jg, doublereal* xg, doublereal* rg,
}
}
XML_Node& OutletRes1D::
save(XML_Node& o, const doublereal* const soln)
{
@ -704,15 +679,10 @@ restore(const XML_Node& dom, doublereal* soln, int loglevel)
resize(1,1);
}
//-----------------------------------------------------------
//
// Surf1D
//
//-----------------------------------------------------------
string Surf1D::componentName(size_t n) const
{
switch (n) {
@ -739,7 +709,6 @@ init()
setTolerances(1, &rtol, 1, &atol);
}
void Surf1D::
eval(size_t jg, doublereal* xg, doublereal* rg,
integer* diagg, doublereal rdt)
@ -793,17 +762,10 @@ restore(const XML_Node& dom, doublereal* soln, int loglevel)
resize(1,1);
}
//-----------------------------------------------------------
//
// ReactingSurf1D
//
//-----------------------------------------------------------
string ReactingSurf1D::componentName(size_t n) const
{
if (n == 0) {
@ -842,7 +804,6 @@ init()
setTolerances(m_nv, DATA_PTR(rtol), m_nv, DATA_PTR(atol));
}
void ReactingSurf1D::
eval(size_t jg, doublereal* xg, doublereal* rg,
integer* diagg, doublereal rdt)
@ -967,4 +928,3 @@ restore(const XML_Node& dom, doublereal* soln, int loglevel)
resize(m_nsp+1,1);
}
}

4
src/oneD/refine.cpp
View file

@ -8,7 +8,6 @@ using namespace std;
namespace Cantera
{
static void r_drawline()
{
string s(78,'#');
@ -26,11 +25,9 @@ Refiner::Refiner(Domain1D& domain) :
m_thresh = std::sqrt(std::numeric_limits<double>::epsilon());
}
int Refiner::analyze(size_t n, const doublereal* z,
const doublereal* x)
{
if (n >= m_npmax) {
writelog("max number of grid points reached ("+int2str(m_npmax)+".\n");
return -2;
@ -214,7 +211,6 @@ void Refiner::show()
}
}
int Refiner::getNewGrid(int n, const doublereal* z,
int nn, doublereal* zn)
{