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

Added a way to specify weighting norms for the residuals.

Browse source
This commit is contained in:
Harry Moffat 2012-02-08 15:59:11 +00:00
parent fdd37f5f2d
commit cb4048249c
2 changed files with 145 additions and 27 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

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

@ -146,6 +146,9 @@ namespace Cantera {
atolBase_(1.0E-10),
m_ydot_nm1(0),
atolk_(0),
userResidAtol_(0),
userResidRtol_(1.0E-3),
checkUserResidualTols_(0),
m_print_flag(0),
m_ScaleSolnNormToResNorm(0.001),
jacCopyPtr_(0),
@ -267,6 +270,9 @@ namespace Cantera {
atolBase_(1.0E-10),
m_ydot_nm1(0),
atolk_(0),
userResidAtol_(0),
userResidRtol_(1.0E-3),
checkUserResidualTols_(0),
m_print_flag(0),
m_ScaleSolnNormToResNorm(0.001),
jacCopyPtr_(0),
@ -368,6 +374,9 @@ namespace Cantera {
rtol_ = right.rtol_;
atolBase_ = right.atolBase_;
atolk_ = right.atolk_;
userResidAtol_ = right.userResidAtol_;
userResidRtol_ = right.userResidRtol_;
checkUserResidualTols_ = right.checkUserResidualTols_;
m_print_flag = right.m_print_flag;
m_ScaleSolnNormToResNorm = right.m_ScaleSolnNormToResNorm;
@ -736,7 +745,7 @@ namespace Cantera {
int irow, jcol;
int ku, kl;
int ivec[2];
int n = jac.nRowsAndStruct(ivec);
jac.nRowsAndStruct(ivec);
double *colP_j;
/*
@ -882,15 +891,25 @@ namespace Cantera {
{
if (! jacCopyPtr_->factored()) {
doublereal sum = 0.0;
for (int irow = 0; irow < neq_; irow++) {
m_residWts[irow] = m_rowWtScales[irow] / neq_;
sum += m_residWts[irow];
}
sum /= neq_;
for (int irow = 0; irow < neq_; irow++) {
m_residWts[irow] = (m_residWts[irow] + atolBase_ * atolBase_ * sum);
if (checkUserResidualTols_ != 1) {
doublereal sum = 0.0;
for (int irow = 0; irow < neq_; irow++) {
m_residWts[irow] = m_rowWtScales[irow] / neq_;
sum += m_residWts[irow];
}
sum /= neq_;
for (int irow = 0; irow < neq_; irow++) {
m_residWts[irow] = (m_residWts[irow] + atolBase_ * atolBase_ * sum);
}
if (checkUserResidualTols_ == 2) {
for (int irow = 0; irow < neq_; irow++) {
m_residWts[irow] = MIN(m_residWts[irow], userResidAtol_[irow] + userResidRtol_ * m_rowWtScales[irow] / neq_);
}
}
} else {
for (int irow = 0; irow < neq_; irow++) {
m_residWts[irow] = userResidAtol_[irow] + userResidRtol_ * m_rowWtScales[irow] / neq_;
}
}
@ -919,6 +938,7 @@ namespace Cantera {
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 {
@ -2232,7 +2252,7 @@ namespace Cantera {
* @param alpha Relative length along the dog length that you are on.
* @param deltaX Vector to be filled up
*/
void NonlinearSolver::fillDogLegStep(int leg, doublereal alpha, std::vector<doublereal> & deltaX) const {
void NonlinearSolver::fillDogLegStep(int leg, doublereal alpha, std::vector<doublereal> & deltaX) const {
if (leg == 0) {
for (int i = 0; i < neq_; i++) {
deltaX[i] = alpha * deltaX_CP_[i];
@ -2343,7 +2363,8 @@ namespace Cantera {
* Maximum decrease in variable in any one newton iteration:
* factor of 5
*/
doublereal NonlinearSolver::boundStep(const doublereal * const y, const doublereal * const step0) {
doublereal NonlinearSolver::boundStep(const doublereal * const y, const doublereal * const step0)
{
int i, i_lower = -1;
doublereal fbound = 1.0, f_bounds = 1.0;
doublereal ff, y_new;
@ -2611,7 +2632,8 @@ namespace Cantera {
}
if (stepNorm_2 < 1.0) {
if (m_print_flag >= 4 ) {
printf("\t dampStep(): current trial step accepted and soln converged retnTrial = %d, its = %d, damp = %g\n", 0, m+1, ff);
printf("\t dampStep(): current trial step accepted and soln converged retnTrial ="
"%d, its = %d, damp = %g\n", 0, m+1, ff);
}
return 0;
}
@ -2715,7 +2737,8 @@ namespace Cantera {
* OK, we have the step0. Now, ask the question whether it satisfies the acceptance criteria
* as a good step. The overall outcome is returned in the variable info.
*/
info = decideStep(time_curr, dogLegID_, dogLegAlpha_, y_n_curr, ydot_n_curr, step_1, y_n_1, ydot_n_1, trustDeltaOld);
info = decideStep(time_curr, dogLegID_, dogLegAlpha_, y_n_curr, ydot_n_curr, step_1,
y_n_1, ydot_n_1, trustDeltaOld);
m_normResid_Bound = m_normResid_1;
/*
@ -2825,9 +2848,11 @@ namespace Cantera {
* -2 Current value of the solution vector caused a residual error in its evaluation.
* Step is a failure, and the step size must be reduced in order to proceed further.
*/
int NonlinearSolver::decideStep(const doublereal time_curr, int leg, doublereal alpha, const doublereal * const y_n_curr,
int NonlinearSolver::decideStep(const doublereal time_curr, int leg, doublereal alpha,
const doublereal * const y_n_curr,
const doublereal * const ydot_n_curr, const std::vector<doublereal> & step_1,
const doublereal * const y_n_1, const doublereal * const ydot_n_1, doublereal trustDeltaOld)
const doublereal * const y_n_1, const doublereal * const ydot_n_1,
doublereal trustDeltaOld)
{
int retn = 2;
bool goodStep = false;
@ -2889,11 +2914,13 @@ namespace Cantera {
m_normResid_1 = m_normResidTrial;
retn = 0;
if (m_print_flag >= 4) {
printf("\t\t decideStep: Norm Residual(leg=%1d, alpha=%10.2E) = %11.4E passes\n", dogLegID_, dogLegAlpha_, m_normResidTrial);
printf("\t\t decideStep: Norm Residual(leg=%1d, alpha=%10.2E) = %11.4E passes\n",
dogLegID_, dogLegAlpha_, m_normResidTrial);
}
} else {
if (m_print_flag >= 4) {
printf("\t\t decideStep: Norm Residual(leg=%1d, alpha=%10.2E) = %11.4E failes\n", dogLegID_, dogLegAlpha_, m_normResidTrial);
printf("\t\t decideStep: Norm Residual(leg=%1d, alpha=%10.2E) = %11.4E failes\n",
dogLegID_, dogLegAlpha_, m_normResidTrial);
}
trustDelta_ *= 0.33;
CurrentTrustFactor_ *= 0.33;
@ -3994,14 +4021,26 @@ namespace Cantera {
NonlinearSolver::computeResidWts()
{
ResidWtsReevaluated_ = true;
doublereal sum = 0.0;
for (int i = 0; i < neq_; i++) {
m_residWts[i] = m_rowWtScales[i] / neq_;
sum += m_residWts[i];
}
sum /= neq_;
for (int i = 0; i < neq_; i++) {
m_residWts[i] = m_ScaleSolnNormToResNorm * (m_residWts[i] + atolBase_ * atolBase_ * sum);
if (checkUserResidualTols_ == 1) {
for (int i = 0; i < neq_; i++) {
m_residWts[i] = userResidAtol_[i] + userResidRtol_ * m_rowWtScales[i] / neq_;
}
} else {
doublereal sum = 0.0;
for (int i = 0; i < neq_; i++) {
m_residWts[i] = m_rowWtScales[i] / neq_;
sum += m_residWts[i];
}
sum /= neq_;
for (int i = 0; i < neq_; i++) {
m_residWts[i] = m_ScaleSolnNormToResNorm * (m_residWts[i] + atolBase_ * atolBase_ * sum);
}
if (checkUserResidualTols_ == 2) {
for (int i = 0; i < neq_; i++) {
double uR = userResidAtol_[i] + userResidRtol_ * m_rowWtScales[i] / neq_;
m_residWts[i] = MIN(m_residWts[i], uR);
}
}
}
}
//=====================================================================================================================
@ -4093,7 +4132,7 @@ namespace Cantera {
//=====================================================================================================================
// Set the relative tolerances for the solution variables
/*
* Set the relative tolerances used in the calculation
* Set the relative tolerances used in the calculation for the solution variables.
*
* @param rtol single double
*/
@ -4102,6 +4141,40 @@ namespace Cantera {
rtol_ = rtol;
}
//=====================================================================================================================
// Set the relative and absolute tolerances for the Residual norm comparisons, if used
/*
*
* residWeightNorm[i] = residAtol[i] + residRtol * m_rowWtScales[i] / neq
*
* @param residNormHandling Parameter that sets the default handling of the residual norms
* 0 The residual weighting vector is calculated to make sure that the solution
* norms are roughly 1 when the residual norm is roughly 1.
* This is the default if this routine is not called.
* 1 Use the user residual norm specified by the parameters in this routine
* 2 Use the minimum value of the residual weights calculcated by method 1 and 2.
* This is the default if this routine is called and this parameter isn't specified.
*/
void NonlinearSolver::setResidualTols(double residRtol, double * residATol, int residNormHandling)
{
if (residNormHandling < 0 || residNormHandling > 2) {
throw CanteraError("NonlinearSolver::setResidualTols()",
"Unknown int for residNormHandling");
}
checkUserResidualTols_ = residNormHandling;
userResidRtol_ = residRtol;
if (residATol) {
userResidAtol_.resize(neq_);
for (int i = 0; i < neq_; i++) {
userResidAtol_[i] = residATol[i];
}
} else {
if (residNormHandling ==1 || residNormHandling == 2) {
throw CanteraError("NonlinearSolver::setResidualTols()",
"Must set residATol vector");
}
}
}
//=====================================================================================================================
void NonlinearSolver::setPrintLvl(int printLvl)
{
m_print_flag = printLvl;

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

@ -667,6 +667,35 @@ namespace Cantera {
*/
void setRtol(const doublereal rtol);
//! Set the relative and absolute tolerances for the Residual norm comparisons, if used
/*!
* Residual norms are used to calculate convergence within the nonlinear solver, since
* these are the norms that are associated with convergence proofs, especially for ill-conditioned systems.
* Usually the residual weights for each row are calculated by the program such that they
* correlate with the convergence requirements on the solution variables input by the user using
* the routines setAtol() and setRtol().
* The residual weights are essentially calculated from the value
*
* residWeightNorm[i] = m_ScaleSolnNormToResNorm * sum_j ( fabs(A_i,j) ewt(j))
*
* The factor, m_ScaleSolnNormToResNorm, is computed periodically to ensure that the solution norms
* and the residual norms are converging at the same time and thus accounts for some-illconditioning issues
* but not all.
*
* The user specified tolerance for the residual is given by the following quantity
*
* residWeightNorm[i] = residAtol[i] + residRtol * m_rowWtScales[i] / neq
*
* @param residNormHandling Parameter that sets the default handling of the residual norms
* 0 The residual weighting vector is calculated to make sure that the solution
* norms are roughly 1 when the residual norm is roughly 1.
* This is the default if this routine is not called.
* 1 Use the user residual norm specified by the parameters in this routine
* 2 Use the minimum value of the residual weights calculcated by method 1 and 2.
* This is the default if this routine is called and this parameter isn't specified.
*/
void setResidualTols(double residRtol, double * residATol, int residNormHandling = 2);
//! Set the value of the maximum # of newton iterations
/*!
* @param maxNewtIts Maximum number of newton iterations
@ -878,6 +907,8 @@ namespace Cantera {
*/
void setSolverScheme(int doDogLeg, int doAffineSolve);
/*
* -----------------------------------------------------------------------------------------------------------------
* MEMBER DATA
@ -1082,6 +1113,20 @@ namespace Cantera {
*/
std::vector<doublereal> atolk_;
//! absolute tolerance in the unscaled solution unknowns
std::vector<doublereal> userResidAtol_;
//! absolute tolerance in the unscaled solution unknowns
doublereal userResidRtol_;
//! Check the residual tolerances explictly against user input
/*!
* 0 Don't calculate residual weights from residual tolerance inputs
* 1 Calculate residual weights from residual tolerance inputs only
* 2 Calculate residual weights from a minimum of the solution error weights process and the direct residual tolerance inputs
*/
int checkUserResidualTols_;
//! Determines the level of printing for each time step.
/*!
* 0 -> absolutely nothing is printed for a single time step.