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

More cleanup

Browse source
This commit is contained in:
Harry Moffat 2011-09-09 21:32:50 +00:00
parent 4d2e6b65bb
commit 93d5d43a19
2 changed files with 139 additions and 69 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

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

@ -70,7 +70,7 @@ namespace Cantera {
printf("\n"); printf("\n");
} }
bool NonlinearSolver::m_TurnOffTiming(false); bool NonlinearSolver::s_TurnOffTiming(false);
#ifdef DEBUG_NUMJAC #ifdef DEBUG_NUMJAC
bool NonlinearSolver::s_print_NumJac(true); bool NonlinearSolver::s_print_NumJac(true);
@ -98,7 +98,8 @@ namespace Cantera {
m_ewt(0), m_ewt(0),
m_manualDeltaStepSet(0), m_manualDeltaStepSet(0),
m_deltaStepMinimum(0), m_deltaStepMinimum(0),
m_y_n(0), m_y_n_curr(0),
m_ydot_n_curr(0),
m_y_nm1(0), m_y_nm1(0),
ydot_new(0), ydot_new(0),
m_colScales(0), m_colScales(0),
@ -106,6 +107,7 @@ namespace Cantera {
m_rowWtScales(0), m_rowWtScales(0),
m_resid(0), m_resid(0),
m_wksp(0), m_wksp(0),
m_wksp_2(0),
m_residWts(0), m_residWts(0),
m_normResid0(0.0), m_normResid0(0.0),
m_normResidFRaw(0.0), m_normResidFRaw(0.0),
@ -162,7 +164,8 @@ namespace Cantera {
m_ewt.resize(neq_, rtol_); m_ewt.resize(neq_, rtol_);
m_deltaStepMinimum.resize(neq_, 0.001); m_deltaStepMinimum.resize(neq_, 0.001);
m_deltaStepMaximum.resize(neq_, 1.0E10); m_deltaStepMaximum.resize(neq_, 1.0E10);
m_y_n.resize(neq_, 0.0); m_y_n_curr.resize(neq_, 0.0);
m_ydot_n_curr.resize(neq_, 0.0);
m_y_nm1.resize(neq_, 0.0); m_y_nm1.resize(neq_, 0.0);
ydot_new.resize(neq_, 0.0); ydot_new.resize(neq_, 0.0);
m_colScales.resize(neq_, 1.0); m_colScales.resize(neq_, 1.0);
@ -170,6 +173,7 @@ namespace Cantera {
m_rowWtScales.resize(neq_, 1.0); m_rowWtScales.resize(neq_, 1.0);
m_resid.resize(neq_, 0.0); m_resid.resize(neq_, 0.0);
m_wksp.resize(neq_, 0.0); m_wksp.resize(neq_, 0.0);
m_wksp_2.resize(neq_, 0.0);
m_residWts.resize(neq_, 0.0); m_residWts.resize(neq_, 0.0);
atolk_.resize(neq_, atolBase_); atolk_.resize(neq_, atolBase_);
deltaX_Newton_.resize(neq_, 0.0); deltaX_Newton_.resize(neq_, 0.0);
@ -198,7 +202,8 @@ namespace Cantera {
m_ewt(0), m_ewt(0),
m_manualDeltaStepSet(0), m_manualDeltaStepSet(0),
m_deltaStepMinimum(0), m_deltaStepMinimum(0),
m_y_n(0), m_y_n_curr(0),
m_ydot_n_curr(0),
m_y_nm1(0), m_y_nm1(0),
ydot_new(0), ydot_new(0),
m_colScales(0), m_colScales(0),
@ -206,6 +211,7 @@ namespace Cantera {
m_rowWtScales(0), m_rowWtScales(0),
m_resid(0), m_resid(0),
m_wksp(0), m_wksp(0),
m_wksp_2(0),
m_residWts(0), m_residWts(0),
m_normResid0(0.0), m_normResid0(0.0),
m_normResidFRaw(0.0), m_normResidFRaw(0.0),
@ -277,7 +283,8 @@ namespace Cantera {
m_ewt = right.m_ewt; m_ewt = right.m_ewt;
m_manualDeltaStepSet = right.m_manualDeltaStepSet; m_manualDeltaStepSet = right.m_manualDeltaStepSet;
m_deltaStepMinimum = right.m_deltaStepMinimum; m_deltaStepMinimum = right.m_deltaStepMinimum;
m_y_n = right.m_y_n; m_y_n_curr = right.m_y_n_curr;
m_ydot_n_curr = right.m_ydot_n_curr;
m_y_nm1 = right.m_y_nm1; m_y_nm1 = right.m_y_nm1;
ydot_new = right.ydot_new; ydot_new = right.ydot_new;
m_colScales = right.m_colScales; m_colScales = right.m_colScales;
@ -285,6 +292,7 @@ namespace Cantera {
m_rowWtScales = right.m_rowWtScales; m_rowWtScales = right.m_rowWtScales;
m_resid = right.m_resid; m_resid = right.m_resid;
m_wksp = right.m_wksp; m_wksp = right.m_wksp;
m_wksp_2 = right.m_wksp_2;
m_residWts = right.m_residWts; m_residWts = right.m_residWts;
m_normResid0 = right.m_normResid0; m_normResid0 = right.m_normResid0;
m_normResidFRaw = right.m_normResidFRaw; m_normResidFRaw = right.m_normResidFRaw;
@ -472,7 +480,7 @@ namespace Cantera {
error = delta_y[i] / m_ewt[i]; error = delta_y[i] / m_ewt[i];
normContrib = sqrt(error * error); normContrib = sqrt(error * error);
printf("\t\t %4d %12.4e | %12.4e %12.4e %12.4e %12.4e\n", i, normContrib/sqrt((double)neq_), printf("\t\t %4d %12.4e | %12.4e %12.4e %12.4e %12.4e\n", i, normContrib/sqrt((double)neq_),
delta_y[i], m_y_n[i], m_y_n[i] + dampFactor * delta_y[i], m_ewt[i]); delta_y[i], m_y_n_curr[i], m_y_n_curr[i] + dampFactor * delta_y[i], m_ewt[i]);
} }
} }
@ -492,7 +500,7 @@ namespace Cantera {
* out to standard output. * out to standard output.
*/ */
doublereal NonlinearSolver::residErrorNorm(const doublereal * const resid, const char * title, const int printLargest, doublereal NonlinearSolver::residErrorNorm(const doublereal * const resid, const char * title, const int printLargest,
const doublereal * const y) const doublereal * const y) const
{ {
int i; int i;
doublereal sum_norm = 0.0, error; doublereal sum_norm = 0.0, error;
@ -580,7 +588,7 @@ namespace Cantera {
m_colScales[i] = 1.0; m_colScales[i] = 1.0;
} }
} }
m_func->calcSolnScales(time_n, DATA_PTR(m_y_n), DATA_PTR(m_y_nm1), DATA_PTR(m_colScales)); m_func->calcSolnScales(time_n, DATA_PTR(m_y_n_curr), DATA_PTR(m_y_nm1), DATA_PTR(m_colScales));
} }
//==================================================================================================================== //====================================================================================================================
// Compute the current residual // Compute the current residual
@ -595,7 +603,7 @@ namespace Cantera {
* -0 or neg value Means an unsuccessful operation * -0 or neg value Means an unsuccessful operation
*/ */
int NonlinearSolver::doResidualCalc(const doublereal time_curr, const int typeCalc, const doublereal * const y_curr, int NonlinearSolver::doResidualCalc(const doublereal time_curr, const int typeCalc, const doublereal * const y_curr,
const doublereal * const ydot_curr, const ResidEval_Type_Enum evalType) const doublereal * const ydot_curr, const ResidEval_Type_Enum evalType) const
{ {
int retn = m_func->evalResidNJ(time_curr, delta_t_n, y_curr, ydot_curr, DATA_PTR(m_resid), evalType); int retn = m_func->evalResidNJ(time_curr, delta_t_n, y_curr, ydot_curr, DATA_PTR(m_resid), evalType);
m_nfe++; m_nfe++;
@ -755,7 +763,7 @@ namespace Cantera {
/* /*
* Compute the undamped Newton step. The residual function is * Compute the undamped Newton step. The residual function is
* evaluated at the current time, t_n, at the current values of the * evaluated at the current time, t_n, at the current values of the
* solution vector, m_y_n, and the solution time derivative, m_ydot_n. * solution vector, m_y_n_curr, and the solution time derivative, m_ydot_n.
* The Jacobian is not recomputed. * The Jacobian is not recomputed.
* *
* A factored jacobian is reused, if available. If a factored jacobian * A factored jacobian is reused, if available. If a factored jacobian
@ -1276,7 +1284,7 @@ namespace Cantera {
double cauchyDistanceNorm = solnErrorNorm(DATA_PTR(deltaX_CP_)); double cauchyDistanceNorm = solnErrorNorm(DATA_PTR(deltaX_CP_));
for (int i = 0; i < neq_; i++) { for (int i = 0; i < neq_; i++) {
mdp::checkFinite(deltaX_CP_[i]); mdp::checkFinite(deltaX_CP_[i]);
y1[i] = m_y_n[i] + ff * deltaX_CP_[i]; y1[i] = m_y_n_curr[i] + ff * deltaX_CP_[i];
} }
/* /*
* Calculate the residual that would result if y1[] were the new solution vector * Calculate the residual that would result if y1[] were the new solution vector
@ -1295,7 +1303,7 @@ namespace Cantera {
double sNewt = solnErrorNorm(DATA_PTR(newtDir)); double sNewt = solnErrorNorm(DATA_PTR(newtDir));
for (int i = 0; i < neq_; i++) { for (int i = 0; i < neq_; i++) {
y1[i] = m_y_n[i] + ff * newtDir[i]; y1[i] = m_y_n_curr[i] + ff * newtDir[i];
} }
/* /*
* Calculate the residual that would result if y1[] were the new solution vector * Calculate the residual that would result if y1[] were the new solution vector
@ -1477,15 +1485,15 @@ namespace Cantera {
//==================================================================================================================== //====================================================================================================================
// Here we print out the residual at various points along the double dogleg, comparing against the quadratic model // Here we print out the residual at various points along the double dogleg, comparing against the quadratic model
// in a table format // in a table format
/*! /*
* @param time_curr INPUT current time * @param time_curr INPUT current time
* @param ydot0 INPUT Current value of the derivative of the solution vector for non-time dependent * @param ydot0 INPUT Current value of the derivative of the solution vector for non-time dependent
* determinations * determinations
* @param ydot1 INPUT Time derivate of solution at the conditions which are evalulated * @param ydot1 INPUT Time derivate of solution at the conditions which are evalulated
*/ */
void NonlinearSolver::residualComparisonLeg(const double time_curr, const double * const ydot0, void NonlinearSolver::residualComparisonLeg(const double time_curr, const double * const ydot0) const {
double * const ydot1) {
double *y1 = DATA_PTR(m_wksp); double *y1 = DATA_PTR(m_wksp);
double *ydot1 = DATA_PTR(m_wksp_2);
double sLen; double sLen;
if (s_print_DogLeg || (doDogLeg_ && m_print_flag > 6)) { if (s_print_DogLeg || (doDogLeg_ && m_print_flag > 6)) {
printf(" residualComparisonLeg() \n"); printf(" residualComparisonLeg() \n");
@ -1503,7 +1511,7 @@ namespace Cantera {
for (int iteration = 0; iteration < (int) alphaT.size(); iteration++) { for (int iteration = 0; iteration < (int) alphaT.size(); iteration++) {
double alpha = alphaT[iteration]; double alpha = alphaT[iteration];
for (int i = 0; i < neq_; i++) { for (int i = 0; i < neq_; i++) {
y1[i] = m_y_n[i] + alpha * deltaX_CP_[i]; y1[i] = m_y_n_curr[i] + alpha * deltaX_CP_[i];
} }
if (solnType_ != NSOLN_TYPE_STEADY_STATE) { if (solnType_ != NSOLN_TYPE_STEADY_STATE) {
calc_ydot(m_order, y1, ydot1); calc_ydot(m_order, y1, ydot1);
@ -1532,7 +1540,7 @@ namespace Cantera {
for (int iteration = 0; iteration < (int) alphaT.size(); iteration++) { for (int iteration = 0; iteration < (int) alphaT.size(); iteration++) {
double alpha = alphaT[iteration]; double alpha = alphaT[iteration];
for (int i = 0; i < neq_; i++) { for (int i = 0; i < neq_; i++) {
y1[i] = m_y_n[i] + (1.0 - alpha) * deltaX_CP_[i]; y1[i] = m_y_n_curr[i] + (1.0 - alpha) * deltaX_CP_[i];
y1[i] += alpha * Nuu_ * deltaX_Newton_[i]; y1[i] += alpha * Nuu_ * deltaX_Newton_[i];
} }
if (solnType_ != NSOLN_TYPE_STEADY_STATE) { if (solnType_ != NSOLN_TYPE_STEADY_STATE) {
@ -1549,7 +1557,7 @@ namespace Cantera {
} }
for (int i = 0; i < neq_; i++) { for (int i = 0; i < neq_; i++) {
y1[i] -= m_y_n[i]; y1[i] -= m_y_n_curr[i];
} }
sLen = solnErrorNorm(DATA_PTR(y1)); sLen = solnErrorNorm(DATA_PTR(y1));
@ -1565,7 +1573,7 @@ namespace Cantera {
for (int iteration = 0; iteration < (int) alphaT.size(); iteration++) { for (int iteration = 0; iteration < (int) alphaT.size(); iteration++) {
double alpha = alphaT[iteration]; double alpha = alphaT[iteration];
for (int i = 0; i < neq_; i++) { for (int i = 0; i < neq_; i++) {
y1[i] = m_y_n[i] + ( Nuu_ + alpha * (1.0 - Nuu_))* deltaX_Newton_[i]; y1[i] = m_y_n_curr[i] + ( Nuu_ + alpha * (1.0 - Nuu_))* deltaX_Newton_[i];
} }
if (solnType_ != NSOLN_TYPE_STEADY_STATE) { if (solnType_ != NSOLN_TYPE_STEADY_STATE) {
calc_ydot(m_order, y1, ydot1); calc_ydot(m_order, y1, ydot1);
@ -1605,7 +1613,7 @@ namespace Cantera {
{ {
for (int i = 0; i < neq_; i++) { for (int i = 0; i < neq_; i++) {
m_deltaStepMinimum[i] = 1000. * atolk_[i]; m_deltaStepMinimum[i] = 1000. * atolk_[i];
m_deltaStepMinimum[i] = MAX(m_deltaStepMinimum[i], 0.1 * fabs(m_y_n[i])); m_deltaStepMinimum[i] = MAX(m_deltaStepMinimum[i], 0.1 * fabs(m_y_n_curr[i]));
} }
} }
//==================================================================================================================== //====================================================================================================================
@ -1729,9 +1737,9 @@ namespace Cantera {
return f_delta_bounds; return f_delta_bounds;
} }
//==================================================================================================================== //====================================================================================================================
//! Calculate the trust region vectors // Calculate the trust region vectors
/*! /*
* The trust region is made up of the trust region vector calculation and the trustDelta_ value * The trust region is made up of the trust region vector calculation and the trustDelta_ value
* We periodically recalculate the trustVector_ values so that they renormalize to the * We periodically recalculate the trustVector_ values so that they renormalize to the
* correct length. * correct length.
@ -1753,7 +1761,7 @@ namespace Cantera {
// we use the old value of the trust region as an indicator // we use the old value of the trust region as an indicator
for (int i = 0; i < neq_; i++) { for (int i = 0; i < neq_; i++) {
oldVal = deltaX_trust_[i]; oldVal = deltaX_trust_[i];
fabsy = fabs(m_y_n[i]); fabsy = fabs(m_y_n_curr[i]);
// First off make sure that each trust region vector is 1/2 the size of each variable or smaller // First off make sure that each trust region vector is 1/2 the size of each variable or smaller
// unless overridden by the deltaStepMininum value. // unless overridden by the deltaStepMininum value.
double newValue = trustDeltaEach * m_ewt[i] / wtSum; double newValue = trustDeltaEach * m_ewt[i] / wtSum;
@ -1859,7 +1867,15 @@ namespace Cantera {
return sum; return sum;
} }
//==================================================================================================================== //====================================================================================================================
int NonlinearSolver::calcTrustIntersection(double trustDelta, double &lambda, double &alpha) const // Given a trust distance, this routine calculates the intersection of the this distance with the
// double dogleg curve
/*
* @param trustDelta (INPUT) Value of the trust distance
* @param lambda (OUTPUT) Returns the internal coordinate of the double dogleg
* @param alpha (OUTPUT) Returns the relative distance along the appropriate leg
* @return leg (OUTPUT) Returns the leg ID (0, 1, or 2)
*/
int NonlinearSolver::calcTrustIntersection(double trustDelta, double &lambda, double &alpha) const
{ {
double dist; double dist;
if (normTrust_Newton_ < trustDelta) { if (normTrust_Newton_ < trustDelta) {
@ -2557,20 +2573,20 @@ namespace Cantera {
// std::vector<doublereal> y_curr(neq_, 0.0); // std::vector<doublereal> y_curr(neq_, 0.0);
std::vector<doublereal> ydot_curr(neq_, 0.0); // std::vector<doublereal> ydot_curr(neq_, 0.0);
std::vector<doublereal> stp(neq_, 0.0); std::vector<doublereal> stp(neq_, 0.0);
std::vector<doublereal> stp1(neq_, 0.0); std::vector<doublereal> stp1(neq_, 0.0);
std::vector<doublereal> y_new(neq_, 0.0); std::vector<doublereal> y_new(neq_, 0.0);
mdp::mdp_copy_dbl_1(DATA_PTR(m_y_n), DATA_PTR(y_comm), neq_); mdp::mdp_copy_dbl_1(DATA_PTR(m_y_n_curr), DATA_PTR(y_comm), neq_);
if (SolnType != NSOLN_TYPE_STEADY_STATE || ydot_comm) { if (SolnType != NSOLN_TYPE_STEADY_STATE || ydot_comm) {
mdp::mdp_copy_dbl_1(DATA_PTR(ydot_curr), ydot_comm, neq_); mdp::mdp_copy_dbl_1(DATA_PTR(m_ydot_n_curr), ydot_comm, neq_);
mdp::mdp_copy_dbl_1(DATA_PTR(ydot_new), ydot_comm, neq_); mdp::mdp_copy_dbl_1(DATA_PTR(ydot_new), ydot_comm, neq_);
} }
// Redo the solution weights every time we enter the function // Redo the solution weights every time we enter the function
createSolnWeights(DATA_PTR(m_y_n)); createSolnWeights(DATA_PTR(m_y_n_curr));
m_normDeltaSoln_Newton = 1.0E1; m_normDeltaSoln_Newton = 1.0E1;
bool frst = true; bool frst = true;
num_newt_its = 0; num_newt_its = 0;
@ -2607,7 +2623,7 @@ namespace Cantera {
* If we are far enough away from the solution, redo the solution weights and the trust vectors. * If we are far enough away from the solution, redo the solution weights and the trust vectors.
*/ */
if (m_normDeltaSoln_Newton > 1.0E2) { if (m_normDeltaSoln_Newton > 1.0E2) {
createSolnWeights(DATA_PTR(m_y_n)); createSolnWeights(DATA_PTR(m_y_n_curr));
#ifdef DEBUG_DOGLEG #ifdef DEBUG_DOGLEG
calcTrustVector(); calcTrustVector();
#else #else
@ -2618,7 +2634,7 @@ namespace Cantera {
} else { } else {
// Do this stuff every 5 iterations // Do this stuff every 5 iterations
if ( (num_newt_its % 5) == 1) { if ( (num_newt_its % 5) == 1) {
createSolnWeights(DATA_PTR(m_y_n)); createSolnWeights(DATA_PTR(m_y_n_curr));
#ifdef DEBUG_DOGLEG #ifdef DEBUG_DOGLEG
calcTrustVector(); calcTrustVector();
#else #else
@ -2629,7 +2645,7 @@ namespace Cantera {
} }
} }
//mdp::mdp_copy_dbl_1(DATA_PTR(m_y_n), DATA_PTR(y_curr), neq_); //mdp::mdp_copy_dbl_1(DATA_PTR(m_y_n_curr), DATA_PTR(y_curr), neq_);
/* /*
* Set default values of Delta bounds constraints * Set default values of Delta bounds constraints
*/ */
@ -2651,7 +2667,8 @@ namespace Cantera {
if (m_print_flag > 3) { if (m_print_flag > 3) {
printf("\tsolve_nonlinear_problem(): Getting a new Jacobian and solving system\n"); printf("\tsolve_nonlinear_problem(): Getting a new Jacobian and solving system\n");
} }
info = beuler_jac(jac, DATA_PTR(m_resid), time_curr, CJ, DATA_PTR(m_y_n), DATA_PTR(ydot_curr), num_newt_its); info = beuler_jac(jac, DATA_PTR(m_resid), time_curr, CJ, DATA_PTR(m_y_n_curr),
DATA_PTR(m_ydot_n_curr), num_newt_its);
if (info == 0) { if (info == 0) {
m = -4; m = -4;
goto done; goto done;
@ -2672,7 +2689,7 @@ namespace Cantera {
/* /*
* Calculate the base residual * Calculate the base residual
*/ */
info = doResidualCalc(time_curr, NSOLN_TYPE_STEADY_STATE, DATA_PTR(m_y_n), DATA_PTR(ydot_curr)); info = doResidualCalc(time_curr, NSOLN_TYPE_STEADY_STATE, DATA_PTR(m_y_n_curr), DATA_PTR(m_ydot_n_curr));
if (info != 1) { if (info != 1) {
if (m_print_flag > 0) { if (m_print_flag > 0) {
printf("\t\t\tsolve_nonlinear_problem(): Residual Calc ERROR %d. Bailing\n", info); printf("\t\t\tsolve_nonlinear_problem(): Residual Calc ERROR %d. Bailing\n", info);
@ -2685,18 +2702,18 @@ namespace Cantera {
* Scale the matrix and the rhs, if they aren't already scaled * Scale the matrix and the rhs, if they aren't already scaled
* Figure out and store the residual scaling factors. * Figure out and store the residual scaling factors.
*/ */
scaleMatrix(jac, DATA_PTR(m_y_n), DATA_PTR(ydot_curr), time_curr); scaleMatrix(jac, DATA_PTR(m_y_n_curr), DATA_PTR(m_ydot_n_curr), time_curr);
/* /*
* Optional print out the initial residual * Optional print out the initial residual
*/ */
if (m_print_flag >= 6) { if (m_print_flag >= 6) {
m_normResid0 = residErrorNorm(DATA_PTR(m_resid), "Initial norm of the residual", 10, DATA_PTR(m_y_n)); m_normResid0 = residErrorNorm(DATA_PTR(m_resid), "Initial norm of the residual", 10, DATA_PTR(m_y_n_curr));
} else if (m_print_flag == 4 || m_print_flag == 5) { } else if (m_print_flag == 4 || m_print_flag == 5) {
m_normResid0 = residErrorNorm(DATA_PTR(m_resid), "Initial norm of the residual", 0, DATA_PTR(m_y_n)); m_normResid0 = residErrorNorm(DATA_PTR(m_resid), "Initial norm of the residual", 0, DATA_PTR(m_y_n_curr));
} else { } else {
m_normResid0 = residErrorNorm(DATA_PTR(m_resid), "Initial norm of the residual", 0, DATA_PTR(m_y_n)); m_normResid0 = residErrorNorm(DATA_PTR(m_resid), "Initial norm of the residual", 0, DATA_PTR(m_y_n_curr));
} }
#ifdef DEBUG_DOGLEG #ifdef DEBUG_DOGLEG
@ -2715,9 +2732,9 @@ namespace Cantera {
// compute the undamped Newton step // compute the undamped Newton step
if (doAffineSolve_) { if (doAffineSolve_) {
info = doAffineNewtonSolve(DATA_PTR(m_y_n), DATA_PTR(ydot_curr), DATA_PTR(deltaX_Newton_), jac); info = doAffineNewtonSolve(DATA_PTR(m_y_n_curr), DATA_PTR(m_ydot_n_curr), DATA_PTR(deltaX_Newton_), jac);
} else { } else {
info = doNewtonSolve(time_curr, DATA_PTR(m_y_n), DATA_PTR(ydot_curr), DATA_PTR(deltaX_Newton_), jac); info = doNewtonSolve(time_curr, DATA_PTR(m_y_n_curr), DATA_PTR(m_ydot_n_curr), DATA_PTR(deltaX_Newton_), jac);
} }
if (info) { if (info) {
@ -2750,28 +2767,28 @@ namespace Cantera {
/* /*
* Filter out bad directions * Filter out bad directions
*/ */
filterNewStep(time_curr, DATA_PTR(m_y_n), DATA_PTR(stp)); filterNewStep(time_curr, DATA_PTR(m_y_n_curr), DATA_PTR(stp));
#ifdef DEBUG_DOGLEG #ifdef DEBUG_DOGLEG
descentComparison(time_curr, DATA_PTR(ydot_curr), DATA_PTR(ydot_new), DATA_PTR(stp)); descentComparison(time_curr, DATA_PTR(m_ydot_n_curr), DATA_PTR(ydot_new), DATA_PTR(stp));
#endif #endif
if (doDogLeg_) { if (doDogLeg_) {
setupDoubleDogleg(); setupDoubleDogleg();
#ifdef DEBUG_DOGLEG #ifdef DEBUG_DOGLEG
residualComparisonLeg(time_curr, DATA_PTR(ydot_curr), DATA_PTR(ydot_new)); residualComparisonLeg(time_curr, DATA_PTR(m_ydot_n_curr));
#endif #endif
m = dampDogLeg(time_curr, DATA_PTR(m_y_n), DATA_PTR(ydot_curr), m = dampDogLeg(time_curr, DATA_PTR(m_y_n_curr), DATA_PTR(m_ydot_n_curr),
stp, DATA_PTR(y_new), DATA_PTR(ydot_new), stp, DATA_PTR(y_new), DATA_PTR(ydot_new),
DATA_PTR(stp1), s1, jac, frst, i_backtracks); DATA_PTR(stp1), s1, jac, frst, i_backtracks);
} }
#ifdef DEBUG_DOGLEG #ifdef DEBUG_DOGLEG
else { else {
residualComparisonLeg(time_curr, DATA_PTR(ydot_curr), DATA_PTR(ydot_new)); residualComparisonLeg(time_curr, DATA_PTR(m_ydot_n_curr));
} }
#endif #endif
@ -2786,7 +2803,7 @@ namespace Cantera {
* s1 * s1
*/ */
if (!doDogLeg_) { if (!doDogLeg_) {
m = dampStep(time_curr, DATA_PTR(m_y_n), DATA_PTR(ydot_curr), m = dampStep(time_curr, DATA_PTR(m_y_n_curr), DATA_PTR(m_ydot_n_curr),
DATA_PTR(stp), DATA_PTR(y_new), DATA_PTR(ydot_new), DATA_PTR(stp), DATA_PTR(y_new), DATA_PTR(ydot_new),
DATA_PTR(stp1), s1, jac, frst, i_backtracks); DATA_PTR(stp1), s1, jac, frst, i_backtracks);
frst = false; frst = false;
@ -2856,10 +2873,10 @@ namespace Cantera {
// Exchange new for curr solutions // Exchange new for curr solutions
if (m >= 0) { if (m >= 0) {
mdp::mdp_copy_dbl_1(DATA_PTR(m_y_n), CONSTD_DATA_PTR(y_new), neq_); mdp::mdp_copy_dbl_1(DATA_PTR(m_y_n_curr), CONSTD_DATA_PTR(y_new), neq_);
if (solnType_ != NSOLN_TYPE_STEADY_STATE) { if (solnType_ != NSOLN_TYPE_STEADY_STATE) {
calc_ydot(m_order, DATA_PTR(m_y_n), DATA_PTR(ydot_curr)); calc_ydot(m_order, DATA_PTR(m_y_n_curr), DATA_PTR(m_ydot_n_curr));
} }
} }
@ -2909,9 +2926,9 @@ namespace Cantera {
} }
mdp::mdp_copy_dbl_1(y_comm, CONSTD_DATA_PTR(m_y_n), neq_); mdp::mdp_copy_dbl_1(y_comm, CONSTD_DATA_PTR(m_y_n_curr), neq_);
if (solnType_ != NSOLN_TYPE_STEADY_STATE) { if (solnType_ != NSOLN_TYPE_STEADY_STATE) {
mdp::mdp_copy_dbl_1(ydot_comm, CONSTD_DATA_PTR(ydot_curr), neq_); mdp::mdp_copy_dbl_1(ydot_comm, CONSTD_DATA_PTR(m_ydot_n_curr), neq_);
} }
num_linear_solves += m_numTotalLinearSolves; num_linear_solves += m_numTotalLinearSolves;
@ -2919,7 +2936,7 @@ namespace Cantera {
doublereal time_elapsed = wc.secondsWC(); doublereal time_elapsed = wc.secondsWC();
if (m_print_flag > 1) { if (m_print_flag > 1) {
if (m > 0) { if (m > 0) {
if (NonlinearSolver::m_TurnOffTiming) { if (NonlinearSolver::s_TurnOffTiming) {
printf("\t\tNonlinear problem solved successfully in %d its\n", printf("\t\tNonlinear problem solved successfully in %d its\n",
num_newt_its); num_newt_its);
} else { } else {
@ -3208,14 +3225,17 @@ namespace Cantera {
return retn; return retn;
} }
//==================================================================================================================== //====================================================================================================================
// Internal function to calculate the time derivative at the new step // Internal function to calculate the time derivative of the solution at the new step
/* /*
* Previously, the user must have supplied information about the previous time step for this routine to
* work as intended.
*
* @param order of the BDF method * @param order of the BDF method
* @param y_curr current value of the solution * @param y_curr current value of the solution
* @param ydot_curr Calculated value of the solution derivative that is consistent with y_curr * @param ydot_curr Calculated value of the solution derivative that is consistent with y_curr
*/ */
void NonlinearSolver:: void NonlinearSolver::
calc_ydot(const int order, const doublereal * const y_curr, doublereal * const ydot_curr) calc_ydot(const int order, const doublereal * const y_curr, doublereal * const ydot_curr) const
{ {
if (!ydot_curr) { if (!ydot_curr) {
return; return;
@ -3235,8 +3255,10 @@ namespace Cantera {
for (i = 0; i < neq_; i++) { for (i = 0; i < neq_; i++) {
ydot_curr[i] = c1 * (y_curr[i] - m_y_nm1[i]) - m_ydot_nm1[i]; ydot_curr[i] = c1 * (y_curr[i] - m_y_nm1[i]) - m_ydot_nm1[i];
} }
throw CanteraError("", "not implemented");
return; return;
default:
throw CanteraError("calc_ydot()", "Case not covered");
} }
} }
//==================================================================================================================== //====================================================================================================================
@ -3391,8 +3413,7 @@ namespace Cantera {
rtol_ = rtol; rtol_ = rtol;
} }
//===================================================================================================================== //=====================================================================================================================
void NonlinearSolver::setPrintLvl(int printLvl)
void NonlinearSolver::setPrintLvl( int printLvl)
{ {
m_print_flag = printLvl; m_print_flag = printLvl;
} }

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

@ -61,6 +61,29 @@ namespace Cantera {
* depending upon whether the direction is straight. * depending upon whether the direction is straight.
* *
* *
* TIME STEP TYPE
*
* The code solves a nonlinear problem. Frequently the nonlinear problem is created from time-dependent
* residual. Whenever you change the solution vector, you are also changing the derivative of the
* solution vector. Therefore, the code has the option of altering ydot, a vector of time derivatives
* of the solution in tandem with the solution vector and then feeding a residual and Jacobian routine
* with the time derivatives as well as the solution. The code has support for a backwards euler method
* and a second order Adams-Bashforth or Trapezoidal Rule.
*
* In order to use these methods, the solver must be initialized with delta_t and m_y_nm1[i] to specify
* the conditions at the previous time step. For second order methods, the time derivative at t_nm1 must
* also be supplied, m_ydot_nm1[i]. Then the solution type NSOLN_TYPE_TIME_DEPENDENT may be used to
* solve the problem.
*
* For steady state problem whose residual doesn't have a solution time derivative in it, you should
* use the NSOLN_TYPE_STEADY_STATE problem type.
*
* We have a NSOLN_TYPE_PSEUDO_TIME_DEPENDENT defined. However, this is not implemented yet. This would
* be a pseudo time dependent calculation, where an optional time derivative could be added in order to
* help equilibrate a nonlinear steady state system. The time transient is not important in and of
* itself. Many physical systems have a time dependence to them that provides a natural way to relax
* the nonlinear system.
*
* *
* @code * @code
* *
@ -159,11 +182,12 @@ namespace Cantera {
* @return Returns the L2 norm of the delta * @return Returns the L2 norm of the delta
*/ */
doublereal residErrorNorm(const doublereal * const resid, const char * title = 0, const int printLargest = 0, doublereal residErrorNorm(const doublereal * const resid, const char * title = 0, const int printLargest = 0,
const doublereal * const y = 0); const doublereal * const y = 0) const;
//! Compute the current residual //! Compute the current residual
/*! /*!
* The current value of the residual is storred in the internal work array m_resid. * The current value of the residual is storred in the internal work array m_resid, which is defined
* as mutable
* *
* @param time_curr Value of the time * @param time_curr Value of the time
* @param typeCalc Type of the calculation * @param typeCalc Type of the calculation
@ -178,7 +202,7 @@ namespace Cantera {
*/ */
int doResidualCalc(const doublereal time_curr, const int typeCalc, const doublereal * const y_curr, int doResidualCalc(const doublereal time_curr, const int typeCalc, const doublereal * const y_curr,
const doublereal * const ydot_curr, const doublereal * const ydot_curr,
const ResidEval_Type_Enum evalType = Base_ResidEval); const ResidEval_Type_Enum evalType = Base_ResidEval) const;
//! Compute the undamped Newton step //! Compute the undamped Newton step
/*! /*!
@ -351,13 +375,16 @@ namespace Cantera {
//! Return an editable vector of the high bounds constraints //! Return an editable vector of the high bounds constraints
std::vector<double> & highBoundsConstraintVector(); std::vector<double> & highBoundsConstraintVector();
//! Internal function to calculate the time derivative at the new step //! Internal function to calculate the time derivative of the solution at the new step
/*! /*!
* Previously, the user must have supplied information about the previous time step for this routine to
* work as intended.
*
* @param order of the BDF method * @param order of the BDF method
* @param y_curr current value of the solution * @param y_curr current value of the solution
* @param ydot_curr Calculated value of the solution derivative that is consistent with y_curr * @param ydot_curr Calculated value of the solution derivative that is consistent with y_curr
*/ */
void calc_ydot(const int order, const doublereal * const y_curr, doublereal * const ydot_curr); void calc_ydot(const int order, const doublereal * const y_curr, doublereal * const ydot_curr) const;
//! Function called to evaluate the jacobian matrix and the current //! Function called to evaluate the jacobian matrix and the current
//! residual vector at the current time step //! residual vector at the current time step
@ -630,6 +657,14 @@ namespace Cantera {
*/ */
int lambdaToLeg(const double lambda, double &alpha) const; int lambdaToLeg(const double lambda, double &alpha) const;
//! Given a trust distance, this routine calculates the intersection of the this distance with the
//! double dogleg curve
/*!
* @param trustDelta (INPUT) Value of the trust distance
* @param lambda (OUTPUT) Returns the internal coordinate of the double dogleg
* @param alpha (OUTPUT) Returns the relative distance along the appropriate leg
* @return leg (OUTPUT) Returns the leg ID (0, 1, or 2)
*/
int calcTrustIntersection(double trustVal, double &lambda, double &alpha) const; int calcTrustIntersection(double trustVal, double &lambda, double &alpha) const;
//! Initialize the size of the trust vector. //! Initialize the size of the trust vector.
@ -686,7 +721,14 @@ namespace Cantera {
*/ */
double expectedResidLeg(int leg, doublereal alpha) const; double expectedResidLeg(int leg, doublereal alpha) const;
void residualComparisonLeg(const double time_curr, const double * const ydot0, double * const ydot1); //! Here we print out the residual at various points along the double dogleg, comparing against the quadratic model
//! in a table format
/*
* @param time_curr INPUT current time
* @param ydot0 INPUT Current value of the derivative of the solution vector for non-time dependent
* determinations
*/
void residualComparisonLeg(const double time_curr, const double * const ydot0) const;
//! Set the print level from the rootfinder //! Set the print level from the rootfinder
/*! /*!
@ -750,7 +792,11 @@ namespace Cantera {
std::vector<doublereal> m_deltaStepMaximum; std::vector<doublereal> m_deltaStepMaximum;
//! Vector containing the current solution vector within the nonlinear solver //! Vector containing the current solution vector within the nonlinear solver
std::vector<doublereal> m_y_n; std::vector<doublereal> m_y_n_curr;
//! Vector containing the time derivative of the current solution vector within the nonlinear solver
//! (where applicable)
std::vector<doublereal> m_ydot_n_curr;
//! Vector containing the solution at the previous time step //! Vector containing the solution at the previous time step
std::vector<doublereal> m_y_nm1; std::vector<doublereal> m_y_nm1;
@ -777,11 +823,14 @@ namespace Cantera {
std::vector<doublereal> m_rowWtScales; std::vector<doublereal> m_rowWtScales;
//! Value of the residual for the nonlinear problem //! Value of the residual for the nonlinear problem
std::vector<doublereal> m_resid; mutable std::vector<doublereal> m_resid;
//! Workspace of length neq_ //! Workspace of length neq_
mutable std::vector<doublereal> m_wksp; mutable std::vector<doublereal> m_wksp;
//! Workspace of length neq_
mutable std::vector<doublereal> m_wksp_2;
/***************************************************************************************** /*****************************************************************************************
* INTERNAL WEIGHTS FOR TAKING SOLUTION NORMS * INTERNAL WEIGHTS FOR TAKING SOLUTION NORMS
******************************************************************************************/ ******************************************************************************************/
@ -809,7 +858,7 @@ namespace Cantera {
doublereal m_normResidPoints[15]; doublereal m_normResidPoints[15];
//! Boolean indicating whether we should scale the residual //! Boolean indicating whether we should scale the residual
bool m_resid_scaled; mutable bool m_resid_scaled;
/***************************************************************************************** /*****************************************************************************************
* INTERNAL BOUNDARY INFO FOR SOLUTIONS * INTERNAL BOUNDARY INFO FOR SOLUTIONS
@ -831,7 +880,7 @@ namespace Cantera {
doublereal delta_t_n; doublereal delta_t_n;
//! Counter for the total number of function evaluations //! Counter for the total number of function evaluations
int m_nfe; mutable int m_nfe;
/*********************************************************************************************** /***********************************************************************************************
* MATRIX INFORMATION * MATRIX INFORMATION
@ -999,7 +1048,7 @@ namespace Cantera {
/*! /*!
* Necessary to do for test suites * Necessary to do for test suites
*/ */
static bool m_TurnOffTiming; static bool s_TurnOffTiming;
//! Turn on or off printing of the Jacobian //! Turn on or off printing of the Jacobian
static bool s_print_NumJac; static bool s_print_NumJac;