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

Modifications to sovleProb to get it to cross the origin

Browse source
This commit is contained in:
Harry Moffat 2010-07-27 20:07:24 +00:00
parent 71167df80f
commit 6490892f95
2 changed files with 45 additions and 27 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

45
Cantera/src/numerics/solveProb.cpp
View file

@ -57,7 +57,7 @@ namespace Cantera {
solveProb::solveProb(ResidEval* resid) :
m_residFunc(resid),
m_neq(0),
m_atol(1.0E-15),
m_atol(0),
m_rtol(1.0E-4),
m_maxstep(1000),
m_ioflag(0)
@ -67,7 +67,7 @@ namespace Cantera {
// Dimension solution vector
int dim1 = MAX(1, m_neq);
m_atol.resize(dim1, 0.0);
m_atol.resize(dim1, 1.0E-9);
m_netProductionRatesSave.resize(dim1, 0.0);
m_numEqn1.resize(dim1, 0.0);
m_numEqn2.resize(dim1, 0.0);
@ -102,7 +102,7 @@ namespace Cantera {
* proportional to their production rates.
*/
int solveProb::solve(int ifunc, doublereal time_scale,
doublereal reltol, doublereal abstol)
doublereal reltol)
{
doublereal EXTRA_ACCURACY = 0.001;
if (ifunc == SOLVEPROB_JACOBIAN) {
@ -164,7 +164,7 @@ namespace Cantera {
if (m_ioflag) {
print_header(m_ioflag, ifunc, time_scale, reltol, abstol,
print_header(m_ioflag, ifunc, time_scale, reltol,
DATA_PTR(m_netProductionRatesSave));
}
@ -323,9 +323,7 @@ namespace Cantera {
for (irow = 0; irow < m_neq; irow++) {
m_CSolnSP[irow] -= damp * m_resid[irow];
}
for (irow = 0; irow < m_neq; irow++) {
m_CSolnSP[irow] = MAX(0.0, m_CSolnSP[irow]);
}
if (do_time) t_real += damp/inv_t;
@ -513,6 +511,10 @@ namespace Cantera {
* - Only going to allow x[i] to converge to the top and bottom bounds by a
* single order of magnitude at one time
*/
bool canCrossOrigin = false;
if (topBounds > 0.0 && botBounds < 0.0) {
canCrossOrigin = true;
}
xtop = topBounds - 0.1 * fabs(topBounds - x[i]);
@ -525,13 +527,21 @@ namespace Cantera {
else if (xnew < xbot) {
damp = APPROACH * (x[i] - xbot) / dxneg[i];
*label = i;
} else if (xnew > 3.0*MAX(x[i], 1.0E-10)) {
damp = - 2.0 * MAX(x[i], 1.0E-10) / dxneg[i];
*label = i;
}
double denom = fabs(x[i]) + m_atol[i];
// else if (fabs(xnew) > 2.0*MAX(fabs(x[i]), 1.0E-10)) {
// damp = 0.5 * MAX(fabs(x[i]), 1.0E-9)/ fabs(xnew);
// *label = i;
// }
double denom = fabs(x[i]) + 1.0E5 * m_atol[i];
if ((fabs(delta_x) / denom) > 0.3) {
double newdamp = 0.3 * denom / delta_x;
double newdamp = 0.3 * denom / fabs(delta_x);
if (canCrossOrigin) {
if (xnew * x[i] < 0.0) {
if (fabs(x[i]) < 1.0E8 * m_atol[i]) {
newdamp = 2.0 * fabs(x[i]) / fabs(delta_x);
}
}
}
damp = MIN(damp, newdamp);
}
@ -698,7 +708,7 @@ namespace Cantera {
* Optional printing at the start of the solveProb problem
*/
void solveProb::print_header(int ioflag, int ifunc, doublereal time_scale,
doublereal reltol, doublereal abstol,
doublereal reltol,
doublereal netProdRate[]) {
int damping = 1;
if (ioflag) {
@ -732,7 +742,7 @@ namespace Cantera {
else
printf(" Damping is OFF \n");
printf(" Reltol = %9.3e, Abstol = %9.3e\n", reltol, abstol);
printf(" Reltol = %9.3e, Abstol = %9.3e\n", reltol, m_atol[0]);
}
/*
@ -974,4 +984,11 @@ namespace Cantera {
}
#endif
//================================================================================================
void solveProb::setAtol(const doublereal atol[]) {
for (int k = 0; k < m_neq; k++, k++) {
m_atol[k] = atol[k];
}
}
}

27
Cantera/src/numerics/solveProb.h
View file

@ -188,15 +188,13 @@ namespace Cantera {
* where applicable
*
* @param reltol Relative tolerance to use
* @param abstol absolute tolerance.
*
* @return Returns 1 if the surface problem is successfully solved.
* Returns -1 if the surface problem wasn't solved successfully.
* Note the actual converged solution is returned as part of the
* internal state of the InterfaceKinetics objects.
*/
int solve(int ifunc, doublereal time_scale,
doublereal reltol, doublereal abstol);
int solve(int ifunc, doublereal time_scale, doublereal reltol);
//! Report the current state of the solution
/*!
@ -204,13 +202,25 @@ namespace Cantera {
*/
virtual void reportState(doublereal * const CSoln) const;
//! Set the bottom and top bounds on the solution vector
/*!
* The default is for the bottom is 0.0, while the default for the top is 1.0
*
* @param botBounds Vector of bottom bounds
* @param topBounds vector of top bounds
*/
virtual void setBounds(const doublereal botBounds[], const doublereal topBounds[]);
void setAtol(const doublereal atol[]);
private:
//! Printing routine that gets called at the start of every
//! invocation
virtual void print_header(int ioflag, int ifunc, doublereal time_scale,
doublereal reltol, doublereal abstol,
doublereal reltol,
doublereal netProdRate[]);
#ifdef DEBUG_SOLVEPROB
@ -354,15 +364,6 @@ namespace Cantera {
*/
virtual doublereal calc_damping(doublereal x[], doublereal dxneg[], int dim, int *label);
//! Set the bottom and top bounds on the solution vector
/*!
* The default is for the bottom is 0.0, while the default for the top is 1.0
*
* @param botBounds Vector of bottom bounds
* @param topBounds vector of top bounds
*/
virtual void setBounds(const doublereal botBounds[], const doublereal topBounds[]);
//! residual function pointer to be solved.
ResidEval *m_residFunc;