Update of /cvsroot/cantera/cantera/Cantera/src
In directory sc8-pr-cvs4.sourceforge.net:/tmp/cvs-serv23382
Modified Files:
ChemEquil.cpp ChemEquil.h
Log Message:
Reworking of the solver algorithm within ChemEquil:
Much of the algorithm for relaxing the element potentials
to equilibrium employs the Brinkley algorithm (const T, P)
(p. 129 Smith and Missen). This algorithm does better at
relaxing the system towards steady state from a bad initial
guess.
The old algorithm now does the last little bit and it does
the calculation for fixed pairs other than const T, P, using
the Brinkley algorithm's initial guess.
Modified the damping algorithm in the old method. It now doesn't
create false convergence failures when it is close to the
end of the calculation.
The new solver has been tested against a matrix of initial conditions
for the gri30 mechanism and a 200 species jp8 mechanism.
Yet to do:
This routine should work but has not been checked out against:
phases with electric charge
nonideal phases.
This routine should probably not yet work with
molality based phases.
More work may have to be done with the Brinkley algorithm for
problems which suffer from stoichiometric constraints. A full
Gauss-Jordon elimination step should probably be added to
flag cases where there are effective stoichiometric constraints.
Index: ChemEquil.h
===================================================================
RCS file: /cvsroot/cantera/cantera/Cantera/src/ChemEquil.h,v
retrieving revision 1.18
retrieving revision 1.19
diff -C2 -d -r1.18 -r1.19
*** ChemEquil.h 25 Jul 2005 03:55:33 -0000 1.18
--- ChemEquil.h 20 Oct 2006 21:34:00 -0000 1.19
***************
*** 28,171 ****
namespace Cantera {
! int _equilflag(const char* xy);
/**
! * Chemical equilibrium options. Used internally by class ChemEquil.
*/
! class EquilOpt {
! public:
! EquilOpt() : relTolerance(1.e-10), maxIterations(1000), iterations(0),
! maxStepSize(10.0), propertyPair(TP), contin(false) {}
!
! doublereal relTolerance; ///< Relative tolerance
! int maxIterations; ///< Maximum number of iterations
! int iterations; ///< Iteration counter
! /**
! * Maximum step size. Largest change in any element potential or
! * in log(T) allowed in one Newton step. Default: 10.0
! */
! doublereal maxStepSize;
! /**
! * Property pair flag. Determines which two thermodynamic properties
! * are fixed.
! */
! int propertyPair;
! /**
! * Continuation flag. Set true if the calculation should be
! * initialized from the last calculation. Otherwise, the
! * calculation will be started from scratch and the initial
! * composition and element potentials estimated.
! */
! bool contin;
! };
! template<class M>
! class PropertyCalculator;
! /**
! * @defgroup equil Chemical Equilibrium
! *
! */
/**
! * Class ChemEquil implements a chemical equilibrium solver for
! * single-phase solutions. It is a "non-stoichiometric" solver in
! * the terminology of Smith and Missen, meaning that every
! * intermediate state is a valid chemical equilibrium state, but
! * does not necessarily satisfy the element constraints. In
! * contrast, the solver implemented in class MultiPhaseEquil uses
! * a "stoichiometric" algorithm, in which each intermediate state
! * satisfies the element constraints but is not a state of
! * chemical equilibrium. Non-stoichiometric methods are faster
! * when they converge, but stoichiometric ones tend to be more
! * robust and can be used also for problems with multiple
! * condensed phases. As expected, the ChemEquil solver is faster
! * than MultiPhaseEquil for many single-phase equilibrium
! * problems (particularly if there are only a few elements but
! * vvery many species), but can be less stable. Problem
! * situations include low temperatures where only a few species
! * have non-zero mole fractions, precisely stoichiometric
! * compositions (e.g. 2 H2 + O2). In general, if speed is
! * important, this solver should be tried first, and if it fails
! * then use MultiPhaseEquil.
! * @ingroup equil
! */
! class ChemEquil {
- public:
- ChemEquil();
- virtual ~ChemEquil();
! int equilibrate(thermo_t& s, const char* XY);
! int equilibrate(thermo_t& s, const char* XY, vector_fp& elMoles);
! const vector_fp& elementPotentials() const { return m_lambda; }
! /**
! * Options controlling how the calculation is carried out.
! * @see EquilOptions
! */
! EquilOpt options;
! protected:
! thermo_t* m_phase;
! thermo_t* m_thermo;
! /// number of atoms of element m in species k.
! doublereal nAtoms(int k, int m) const { return m_comp[k*m_mm + m]; }
! void initialize(thermo_t& s);
! void setToEquilState(thermo_t& s,
! const vector_fp& x, doublereal t);
! int setInitialMoles(thermo_t& s);
! int estimateElementPotentials(thermo_t& s,
! vector_fp& lambda);
!
! int dampStep(thermo_t& s, vector_fp& oldx,
! double oldf, vector_fp& grad, vector_fp& step, vector_fp& x,
! double& f, vector_fp& elmols, double xval, double yval );
! void equilResidual(thermo_t& s, const vector_fp& x,
! const vector_fp& elmtotal, vector_fp& resid,
! double xval, double yval);
! void equilJacobian(thermo_t& s, vector_fp& x,
! const vector_fp& elmols, DenseMatrix& jac,
! double xval, double yval);
! void update(const thermo_t& s);
! int m_mm;
! int m_kk;
! int m_skip;
! PropertyCalculator<thermo_t> *m_p1, *m_p2;
! vector_fp m_molefractions;
! vector_fp m_lambda;
! vector_fp m_elementmolefracs;
! vector_fp m_reswork;
! vector_fp m_jwork1;
! vector_fp m_jwork2;
! vector_fp m_comp;
! doublereal m_temp, m_dens;
! doublereal m_p0;
! int m_eloc;
! doublereal m_abscharge;
! doublereal m_startTemp, m_startDens;
! vector_fp m_startSoln;
! vector_fp m_grt;
! vector_fp m_mu_RT;
! vector_int m_component;
! };
}
--- 28,205 ----
namespace Cantera {
! int _equilflag(const char* xy);
!
! /**
! * Chemical equilibrium options. Used internally by class ChemEquil.
! */
! class EquilOpt {
! public:
! EquilOpt() : relTolerance(1.e-8), absElemTol(1.0E-70),maxIterations(1000),
! iterations(0),
! maxStepSize(10.0), propertyPair(TP), contin(false) {}
!
! doublereal relTolerance; ///< Relative tolerance
! doublereal absElemTol; ///< Abs Tol in element number
! int maxIterations; ///< Maximum number of iterations
! int iterations; ///< Iteration counter
/**
! * Maximum step size. Largest change in any element potential or
! * in log(T) allowed in one Newton step. Default: 10.0
*/
! doublereal maxStepSize;
! /**
! * Property pair flag. Determines which two thermodynamic properties
! * are fixed.
! */
! int propertyPair;
! /**
! * Continuation flag. Set true if the calculation should be
! * initialized from the last calculation. Otherwise, the
! * calculation will be started from scratch and the initial
! * composition and element potentials estimated.
! */
! bool contin;
! };
! template<class M>
! class PropertyCalculator;
! /**
! * @defgroup equil Chemical Equilibrium
! *
! */
! /**
! * Class ChemEquil implements a chemical equilibrium solver for
! * single-phase solutions. It is a "non-stoichiometric" solver in
! * the terminology of Smith and Missen, meaning that every
! * intermediate state is a valid chemical equilibrium state, but
! * does not necessarily satisfy the element constraints. In
! * contrast, the solver implemented in class MultiPhaseEquil uses
! * a "stoichiometric" algorithm, in which each intermediate state
! * satisfies the element constraints but is not a state of
! * chemical equilibrium. Non-stoichiometric methods are faster
! * when they converge, but stoichiometric ones tend to be more
! * robust and can be used also for problems with multiple
! * condensed phases. As expected, the ChemEquil solver is faster
! * than MultiPhaseEquil for many single-phase equilibrium
! * problems (particularly if there are only a few elements but
! * vvery many species), but can be less stable. Problem
! * situations include low temperatures where only a few species
! * have non-zero mole fractions, precisely stoichiometric
! * compositions (e.g. 2 H2 + O2). In general, if speed is
! * important, this solver should be tried first, and if it fails
! * then use MultiPhaseEquil.
! * @ingroup equil
! */
! class ChemEquil {
!
! public:
! ChemEquil();
! virtual ~ChemEquil();
!
! int equilibrate(thermo_t& s, const char* XY,
! bool useThermoPhaseElementPotentials = false);
! int equilibrate(thermo_t& s, const char* XY, vector_fp& elMoles,
! bool useThermoPhaseElementPotentials = false);
! const vector_fp& elementPotentials() const { return m_lambda; }
/**
! * Options controlling how the calculation is carried out.
! * @see EquilOptions
! */
! EquilOpt options;
! protected:
! thermo_t* m_phase;
! thermo_t* m_thermo;
+ /// number of atoms of element m in species k.
+ doublereal nAtoms(int k, int m) const { return m_comp[k*m_mm + m]; }
! void initialize(thermo_t& s);
! void setToEquilState(thermo_t& s,
! const vector_fp& x, doublereal t);
! int setInitialMoles(thermo_t& s);
! int estimateElementPotentials(thermo_t& s, vector_fp& lambda);
!
! int estimateEP_Brinkley(thermo_t&s, vector_fp& lambda, vector_fp& elMoles);
! int dampStep(thermo_t& s, vector_fp& oldx,
! double oldf, vector_fp& grad, vector_fp& step, vector_fp& x,
! double& f, vector_fp& elmols, double xval, double yval );
! void equilResidual(thermo_t& s, const vector_fp& x,
! const vector_fp& elmtotal, vector_fp& resid,
! double xval, double yval);
! void equilJacobian(thermo_t& s, vector_fp& x,
! const vector_fp& elmols, DenseMatrix& jac,
! double xval, double yval);
! void update(const thermo_t& s);
! int m_mm;
! int m_kk;
! int m_skip;
! PropertyCalculator<thermo_t> *m_p1, *m_p2;
! /**
! * Current value of the mole fractions in the single phase.
! * -> length = m_kk.
! */
! vector_fp m_molefractions;
! /**
! * Current value of the dimensional element potentials
! * -> length = m_mm
! */
! vector_fp m_lambda;
! /*
! * Current value of the sum of the element abundances given the
! * current element potentials.
! */
! doublereal m_elementTotalSum;
! /*
! * Current value of the element mole fractions. Note these aren't
! * the goal element mole fractions.
! */
! vector_fp m_elementmolefracs;
! vector_fp m_reswork;
! vector_fp m_jwork1;
! vector_fp m_jwork2;
! /*
! * Storage of the element compositions
! * natom(k,m) = m_comp[k*m_mm+ m];
! */
! vector_fp m_comp;
! doublereal m_temp, m_dens;
! doublereal m_p0;
! int m_eloc;
! doublereal m_abscharge;
! doublereal m_startTemp, m_startDens;
! vector_fp m_startSoln;
! vector_fp m_grt;
! vector_fp m_mu_RT;
! vector_int m_component;
! /*
! * element fractional cutoff, below which the element will be
! * zeroed.
! */
! double m_elemFracCutoff;
! bool m_doResPerturb;
! };
}
***************
*** 173,175 ****
#endif
-
--- 207,208 ----
Index: ChemEquil.cpp
===================================================================
RCS file: /cvsroot/cantera/cantera/Cantera/src/ChemEquil.cpp,v
retrieving revision 1.23
retrieving revision 1.24
diff -C2 -d -r1.23 -r1.24
*** ChemEquil.cpp 28 Apr 2006 17:22:23 -0000 1.23
--- ChemEquil.cpp 20 Oct 2006 21:34:00 -0000 1.24
***************
*** 6,9 ****
--- 6,12 ----
* ChemEquil.
*
+ *
+ * $Id$
+ *
* Copyright 2001 California Institute of Technology
*
***************
*** 28,833 ****
[...2366 lines suppressed...]
! }
! n_t *= exp(beta * resid[m_mm]);
!
!
! #ifdef DEBUG_HKM_EPEQUIL
! if (debug_prnt_lvl > 0) {
! printf("(it %d) OLD_SOLUTION NEW SOLUTION (undamped updated)\n", iter);
! for (m = 0; m < m_mm; m++) {
! string eee = eNames[m];
! printf(" %5s %10.5g %10.5g %10.5g\n", eee.c_str(), x_old[m], x[m], resid[m]);
! }
! printf(" n_t %10.5g %10.5g %10.5g \n", x_old[m_mm], n_t, exp(resid[m_mm]));
! }
! #endif
! }
! exit:
! return retn;
! }
!
! } // namespace
|