Second iteration for solveProp. Added a routine to get the

final answer out.
This commit is contained in:
Harry Moffat 2010-05-22 21:50:26 +00:00
parent 45a427fe54
commit d8ae74198d
2 changed files with 30 additions and 46 deletions

View file

@ -48,9 +48,6 @@ namespace Cantera {
# define MIN(x,y) (( (x) < (y) ) ? (x) : (y)) /* min function */
#endif
#ifndef DAMPING
# define DAMPING true
#endif
/***************************************************************************
* solveSP Class Definitinos
@ -152,27 +149,17 @@ namespace Cantera {
* upload the initial conditions
*/
m_residFunc->getInitialConditions(t_real, DATA_PTR(m_CSolnSP), DATA_PTR(m_numEqn1));
/*
* Store the initial guess in the soln vector,
* CSoln, and in an separate vector CSolnInit.
* CSolnSP, and in an separate vector CSolnSPInit.
*/
for (int n = 0; n < m_neq; n++) {
// m_CSolnSP[loc] = m_numEqn1[k];
}
std::copy(m_CSolnSP.begin(), m_CSolnSP.end(), m_CSolnSPInit.begin());
// Calculate the largest species in each phase
// evalSurfLarge(DATA_PTR(m_CSolnSP));
/*
* Get the net production rate of all species in the kinetics manager.
*/
// m_kin->getNetProductionRates(DATA_PTR(m_netProductionRatesSave));
if (m_ioflag) {
print_header(m_ioflag, ifunc, time_scale, DAMPING, reltol, abstol,
print_header(m_ioflag, ifunc, time_scale, reltol, abstol,
DATA_PTR(m_netProductionRatesSave));
}
@ -313,9 +300,9 @@ namespace Cantera {
* in any unknown.
*/
#ifdef DAMPING
damp = calc_damping( DATA_PTR(m_CSolnSP), DATA_PTR(m_resid), m_neq, &label_d);
#endif
damp = calc_damping(DATA_PTR(m_CSolnSP), DATA_PTR(m_resid), m_neq, &label_d);
/*
* Calculate the weighted norm of the update vector
@ -334,7 +321,6 @@ namespace Cantera {
for (irow = 0; irow < m_neq; irow++) {
m_CSolnSP[irow] = MAX(0.0, m_CSolnSP[irow]);
}
updateState(DATA_PTR(m_CSolnSP));
if (do_time) t_real += damp/inv_t;
@ -411,17 +397,14 @@ namespace Cantera {
if (update_norm > 1.0) {
return -1;
}
return 1;
return 0;
}
#undef DAMPING
//================================================================================================
/*
* Update the surface states of the surface phases.
*/
void solveProb::updateState(const doublereal *CSolnSP) {
void solveProb::reportState(doublereal * const CSolnSP) const {
std::copy(m_CSolnSP.begin(), m_CSolnSP.end(), CSolnSP);
}
//================================================================================================
/*
@ -434,8 +417,8 @@ namespace Cantera {
* This routine uses the m_numEqn1 and m_netProductionRatesSave vectors
* as temporary internal storage.
*/
void solveProb::fun_eval(doublereal * resid, const doublereal *CSoln,
const doublereal *CSolnOld, const bool do_time,
void solveProb::fun_eval(doublereal * const resid, const doublereal * const CSoln,
const doublereal * const CSolnOld, const bool do_time,
const doublereal deltaT)
{
if (do_time) {
@ -690,8 +673,9 @@ namespace Cantera {
* Optional printing at the start of the solveProb problem
*/
void solveProb::print_header(int ioflag, int ifunc, doublereal time_scale,
int damping, doublereal reltol, doublereal abstol,
doublereal netProdRate[]) {
doublereal reltol, doublereal abstol,
doublereal netProdRate[]) {
int damping = 1;
if (ioflag) {
printf("\n================================ SOLVEPROB CALL SETUP "
"========================================\n");
@ -864,7 +848,7 @@ namespace Cantera {
#endif
} /* printIteration */
//================================================================================================
//================================================================================================
void solveProb::printFinal(int ioflag, doublereal damp, int label_d, int label_t,
doublereal inv_t, doublereal t_real, int iter,
doublereal update_norm, doublereal resid_norm,

View file

@ -86,12 +86,12 @@ namespace Cantera {
* 4 values:
*
*
* 1: SFLUX_INITIALIZE = This assumes that the initial guess supplied to the
* 1: SOLVEPROB_INITIALIZE = This assumes that the initial guess supplied to the
* routine is far from the correct one. Substantial
* work plus transient time-stepping is to be expected
* to find a solution.
*
* 2: SFLUX_RESIDUAL = Need to solve the nonlinear problem in order to
* 2: SOLVEPROB_RESIDUAL = Need to solve the nonlinear problem in order to
* calculate quantities for a residual calculation
* (Can expect a moderate change in the solution
* vector -> try to solve the system by direct methods
@ -101,14 +101,14 @@ namespace Cantera {
* algorithm to determine when to shut off
* time-stepping.
*
* 3: SFLUX_JACOBIAN = Calculation of the surface problem is due to the
* 3: SOLVEPROB_JACOBIAN = Calculation of the surface problem is due to the
* need for a numerical jacobian for the gas-problem.
* The solution is expected to be very close to the
* initial guess, and extra accuracy is needed because
* solution variables have been delta'd from
* nominal values to create jacobian entries.
*
* 4: SFLUX_TRANSIENT = The transient calculation is performed here for an
* 4: SOLVEPROB_TRANSIENT = The transient calculation is performed here for an
* amount of time specified by "time_scale". It is
* not garraunted to be time-accurate - just stable
* and fairly fast. The solution after del_t time is
@ -193,12 +193,19 @@ namespace Cantera {
int solve(int ifunc, doublereal time_scale,
doublereal reltol, doublereal abstol);
//! Report the current state of the solution
/*!
* @param Report the solution vector for the nonlinear problem
*/
virtual void reportState(doublereal * const CSoln) const;
private:
//! Printing routine that gets called at the start of every
//! invocation
virtual void print_header(int ioflag, int ifunc, doublereal time_scale,
int damping, doublereal reltol, doublereal abstol,
doublereal reltol, doublereal abstol,
doublereal netProdRate[]);
#ifdef DEBUG_SOLVEPROB
@ -289,13 +296,6 @@ namespace Cantera {
doublereal inv_t, doublereal t_real, int iter,
bool do_time);
#endif
/**
* Update the surface states of the surface phases.
*/
virtual void updateState(const doublereal *cSurfSpec);
//! Main Function evalulation
/*!
@ -308,8 +308,8 @@ namespace Cantera {
* @param do_time Calculate a time dependent residual
* @param deltaT Delta time for time dependent problem.
*/
virtual void fun_eval(doublereal* resid, const doublereal *CSolnSP,
const doublereal *CSolnOldSP, const bool do_time, const doublereal deltaT);
virtual void fun_eval(doublereal* const resid, const doublereal * const CSolnSP,
const doublereal * const CSolnOldSP, const bool do_time, const doublereal deltaT);
//! Main routine that calculates the current residual and Jacobian
/*!