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

Update on the NonlinearSolver:

Browse source 
Update to :calcSolnToResNormVector() to align it with documentation.
This commit is contained in:
Harry Moffat 2011-05-30 17:57:53 +00:00
parent 8710440a1f
commit 4d245501cf
2 changed files with 46 additions and 19 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

51
Cantera/src/numerics/NonlinearSolver.cpp
View file

@ -348,17 +348,12 @@ namespace Cantera {
* The program always assumes that atol is specific
* to the solution component
*
* param y vector of the current solution values
* @param y vector of the current solution values
*/
void NonlinearSolver::createSolnWeights(const doublereal * const y) {
for (int i = 0; i < neq_; i++) {
m_ewt[i] = rtol_ * fabs(y[i]) + atolk_[i];
}
#ifdef DEBUG_DOGLEG
// for (int i = 0; i < neq_; i++) {
// m_ewt[i] = 1.0E-4;
// }
#endif
}
//====================================================================================================================
// set bounds constraints for all variables in the problem
@ -665,8 +660,10 @@ namespace Cantera {
}
for (jcol = 0; jcol < neq_; jcol++) {
for (irow = 0; irow < neq_; irow++) {
m_rowScales[irow] += fabs(*jptr);
//m_rowScales[irow] += fabs(*jptr);
if (m_colScaling) {
// This is needed in order to mitgate the change in J_ij carried out just above this loop.
// Alternatively, we could move this loop up to the top
m_rowWtScales[irow] += fabs(*jptr) * m_ewt[jcol] / m_colScales[jcol];
} else {
m_rowWtScales[irow] += fabs(*jptr) * m_ewt[jcol];
@ -699,16 +696,40 @@ namespace Cantera {
}
//====================================================================================================================
// Calculate the scaling factor for translating residual norms into solution norms.
/*
* This routine calls computeResidWts() a couple of times in the calculation of m_ScaleSolnNormToResNorm.
* A more sophisticated routine may do more with signs to get a better value. Perhaps, a series of calculations
* with different signs attached may be in order. Then, m_ScaleSolnNormToResNorm would be calculated
* as the minimum of a series of calculations.
*/
void NonlinearSolver::calcSolnToResNormVector()
{
// double oldVal = m_ScaleSolnNormToResNorm;
// if (m_normSolnFRaw > 1.0E-13) {
// m_ScaleSolnNormToResNorm = 1.0E-2 * m_normResid0 / m_normSolnFRaw * oldVal;
//}
// m_normResid0 = m_normSolnFRaw;
//computeResidWts();
m_ScaleSolnNormToResNorm = 1.0;
//double tmp = residErrorNorm(DATA_PTR(m_resid));
if (! jacCopy_.m_factored) {
m_ScaleSolnNormToResNorm = 1.0;
computeResidWts();
for (int n = 0; n < neq_; n++) {
m_wksp[n] = 0.0;
}
doublereal *jptr = &(*(jacCopy_.begin()));
for (int jcol = 0; jcol < neq_; jcol++) {
for (int irow = 0; irow < neq_; irow++) {
m_wksp[irow] += (*jptr) * m_ewt[jcol];
jptr++;
}
}
double resNormOld = residErrorNorm(DATA_PTR(m_wksp));
if (resNormOld > 0.0) {
m_ScaleSolnNormToResNorm = m_ScaleSolnNormToResNorm * resNormOld;
}
if (m_ScaleSolnNormToResNorm < 1.0E-8) {
m_ScaleSolnNormToResNorm = 1.0E-8;
}
// Recalculate the residual weights now that we know the value of m_ScaleSolnNormToResNorm
computeResidWts();
} else {
throw CanteraError("NonlinearSolver::calcSolnToResNormVector()" , "Logic error");
}
}
//====================================================================================================================
// Compute the undamped Newton step based on the current jacobian and an input rhs

14
Cantera/src/numerics/NonlinearSolver.h
View file

@ -579,8 +579,13 @@ namespace Cantera {
*/
void setMaxNewtIts(const int maxNewtIts);
//! Calculate the scaling factor for translating residual norms into
//! solution norms.
//! Calculate the scaling factor for translating residual norms into solution norms.
/*!
* This routine calls computeResidWts() a couple of times in the calculation of m_ScaleSolnNormToResNorm.
* A more sophisticated routine may do more with signs to get a better value. Perhaps, a series of calculations
* with different signs attached may be in order. Then, m_ScaleSolnNormToResNorm would be calculated
* as the minimum of a series of calculations.
*/
void calcSolnToResNormVector();
//! Calculate the steepest descent direction and the Cauchy Point where the quadratic formulation
@ -758,8 +763,9 @@ namespace Cantera {
//! Weights for normalizing the values of the residuals
/*!
* These are computed if row scaling, m_rowScaling, is turned on. They are calculated currently as the
* sum of the absolute values jacobian multiplied by the solution weight function
* They are calculated as the sum of the absolute values of the jacobian
* multiplied by the solution weight function.
* This is carried out in scaleMatrix().
*/
std::vector<doublereal> m_rowWtScales;