research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions
cantera/include/cantera/numerics/NonlinearSolver.h
Ray Speth 9109888882 Fixed regressions in NonlinearSolver introduced in r1978 
Fixed signed/unsigned type consistency, spelling mistakes, and removed
unnecessary macros.
2013-02-01 23:39:58 +00:00

1324 lines
55 KiB
C++
Raw Blame History

/**
* @file NonlinearSolver.h
* Class that calculates the solution to a nonlinear, dense, set
* of equations (see \ref numerics
* and class \link Cantera::NonlinearSolver NonlinearSolver\endlink).
*/
/*
* $Date$
* $Revision$
*/
/*
* Copywrite 2004 Sandia Corporation. Under the terms of Contract
* DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
* retains certain rights in this software.
* See file License.txt for licensing information.
*/
#ifndef CT_NONLINEARSOLVER_H
#define CT_NONLINEARSOLVER_H
#include "cantera/numerics/ResidJacEval.h"
#include "cantera/numerics/SquareMatrix.h"
namespace Cantera {
//@{
/// @name Constant which determines the type of the nonlinear solve
/*!
* I think steady state is the only option I'm gunning for
*/
//! The nonlinear problem is part of a pseudo time dependent calculation (NOT TESTED)
#define NSOLN_TYPE_PSEUDO_TIME_DEPENDENT 2
//! The nonlinear problem is part of a time dependent calculation
#define NSOLN_TYPE_TIME_DEPENDENT 1
//! The nonlinear problem is part of a steady state calculation
#define NSOLN_TYPE_STEADY_STATE 0
//@}
//@{
/// @name Constant which determines the Return int from the nonlinear solver
/*!
* This int is returned from the nonlinear solver
*/
//! The nonlinear solve is successful.
#define NSOLN_RETN_SUCCESS 1
//! Problem isn't solved yet
#define NSOLN_RETN_CONTINUE 0
//! The nonlinear problem started to take too small an update step. This indicates that either the
//! Jacobian is bad, or a constraint is being bumped up against.
#define NSOLN_RETN_FAIL_STEPTOOSMALL -1
//! The nonlinear problem didn't solve the problem
#define NSOLN_RETN_FAIL_DAMPSTEP -2
//! The nonlinear problem's jacobian is singular
#define NSOLN_RETN_MATRIXINVERSIONERROR -3
//! The nonlinear problem's jacobian formation produced an error
#define NSOLN_RETN_JACOBIANFORMATIONERROR -4
//! The nonlinear problem's base residual produced an error
#define NSOLN_RETN_RESIDUALFORMATIONERROR -5
//! The nonlinear problem's max number of iterations has been exceeded
#define NSOLN_RETN_MAXIMUMITERATIONSEXCEEDED -7
//@}
//@}
//@{
/// @name Constant which determines the type of the Jacobian
//! The jacobian will be calculated from a numerical method
#define NSOLN_JAC_NUM 1
//! The jacobian is calculated from an analytical function
#define NSOLN_JAC_ANAL 2
//@}
//! Class that calculates the solution to a nonlinear system
/*!
* This is a small nonlinear solver that can solve highly nonlinear problems that
* must use a dense matrix to relax the system.
*
* Newton's method is used.
*
* Damping is used extensively when relaxing the system.
*
*
* The basic idea is that we predict a direction that is parameterized by an overall coordinate
* value, beta, from zero to one, This may or may not be the same as the value, damp,
* 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.
*
* MATRIX SCALING
*
*
*
*
* @code
*
*
* NonlinearSolver *nls = new NonlinearSolver(&r1);
*
* int solnType = NSOLN_TYPE_STEADY_STATE ;
*
* nls->setDeltaBoundsMagnitudes(deltaBounds);
*
* nls->solve_nonlinear_problem(solnType, y_comm, ydot_comm, CJ, time_curr, jac,
* num_newt_its, num_linear_solves, numBacktracks,
* loglevelInput);
*
* @endcode
*
*
* @ingroup numerics
*/
class NonlinearSolver {
public:
//! Default constructor
/*!
* @param func Residual and jacobian evaluator function object
*/
NonlinearSolver(ResidJacEval *func);
//!Copy Constructor for the %ThermoPhase object.
/*!
* @param right Item to be copied
*/
NonlinearSolver(const NonlinearSolver &right);
//! Destructor
~NonlinearSolver();
//! Assignment operator
/*!
* This is NOT a virtual function.
*
* @param right Reference to %NonlinearSolver object to be
* copied into the
* current one.
*/
NonlinearSolver& operator=(const NonlinearSolver &right);
//! Create solution weights for convergence criteria
/*!
* We create soln weights from the following formula
*
* wt[i] = rtol * abs(y[i]) + atol[i]
*
* The program always assumes that atol is specific
* to the solution component
*
* @param y vector of the current solution values
*/
void createSolnWeights(const doublereal * const y);
//! L2 norm of the delta of the solution vector
/*!
* calculate the norm of the solution vector. This will
* involve the column scaling of the matrix
*
* The third argument has a default of false. However,
* if true, then a table of the largest values is printed
* out to standard output.
*
* @param delta_y Vector to take the norm of
* @param title Optional title to be printed out
* @param printLargest int indicating how many specific lines should be printed out
* @param dampFactor Current value of the damping factor. Defaults to 1.
* only used for printout out a table.
*
* @return Returns the L2 norm of the delta
*/
doublereal solnErrorNorm(const doublereal * const delta_y, const char * title = 0, int printLargest = 0,
const doublereal dampFactor = 1.0) const;
//! L2 norm of the residual of the equation system
/*!
* Calculate the norm of the residual vector. This may
* involve using the row sum scaling from the matrix problem.
*
* The second argument has a default of false. However,
* if true, then a table of the largest values is printed
* out to standard output.
*
* @param resid Vector of the residuals
* @param title Optional title to be printed out
* @param printLargest Number of specific entries to be printed
* @param y Current value of y - only used for printouts
*
*
* @return Returns the L2 norm of the delta
*/
doublereal residErrorNorm(const doublereal * const resid, const char * title = 0, const int printLargest = 0,
const doublereal * const y = 0) const;
//! Compute the current residual
/*!
* The current value of the residual is stored in the internal work array m_resid, which is defined
* as mutable
*
* @param time_curr Value of the time
* @param typeCalc Type of the calculation
* @param y_curr Current value of the solution vector
* @param ydot_curr Current value of the time derivative of the solution vector
* @param evalType Base evalulation type
* Defaults to Base_ResidEval
*
* @return Returns a flag to indicate that operation is successful.
* 1 Means a successful operation
* -0 or neg value Means an unsuccessful operation
*/
int doResidualCalc(const doublereal time_curr, const int typeCalc, const doublereal * const y_curr,
const doublereal * const ydot_curr,
const ResidEval_Type_Enum evalType = Base_ResidEval) const;
//! Compute the undamped Newton step
/*!
*
* Compute the undamped Newton step. The residual function is
* 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.
* The Jacobian is not recomputed.
*
* A factored jacobian is reused, if available. If a factored jacobian
* is not available, then the jacobian is factored. Before factoring,
* the jacobian is row and column-scaled. Column scaling is not
* recomputed. The row scales are recomputed here, after column
* scaling has been implemented.
*
* @param time_curr Current value of the time
* @param y_curr Current value of the solution
* @param ydot_curr Current value of the solution derivative.
* @param delta_y return value of the raw change in y
* @param jac Jacobian
*
* @return Returns the result code from lapack. A zero means success. Anything
* else indicates a failure.
*/
int doNewtonSolve(const doublereal time_curr, const doublereal * const y_curr,
const doublereal * const ydot_curr, doublereal * const delta_y,
GeneralMatrix& jac);
//! Compute the newton step, either by direct newton's or by solving a close problem that is represented
//! by a Hessian (
/*!
* This is algorith A.6.5.1 in Dennis / Schnabel
*
* Compute the QR decomposition
*
* Compute the undamped Newton step. The residual function is
* 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.
* The Jacobian is not recomputed.
*
* A factored jacobian is reused, if available. If a factored jacobian
* is not available, then the jacobian is factored. Before factoring,
* the jacobian is row and column-scaled. Column scaling is not
* recomputed. The row scales are recomputed here, after column
* scaling has been implemented.
*
* @param y_curr Current value of the solution
* @param ydot_curr Current value of the solution derivative.
* @param delta_y return value of the raw change in y
* @param jac Jacobian
*
* Internal input
* ---------------
* internal m_resid Stored residual is used as input
*
*
* @return Returns the result code from lapack. A zero means success. Anything
* else indicates a failure.
*/
int doAffineNewtonSolve(const doublereal * const y_curr, const doublereal * const ydot_curr,
doublereal * const delta_y, GeneralMatrix& jac);
//! Calculate the length of the current trust region in terms of the solution error norm
/*!
* We carry out a norm of deltaX_trust_ first. Then, we multiply that value
* by trustDelta_
*/
doublereal trustRegionLength() const;
//! Set default deulta bounds amounts
/*!
* Delta bounds are set to 0.01 for all unknowns arbitrarily and capriciously
* Then, for each call to the nonlinear solver
* Then, they are increased to 1000 x atol
* then, they are increased to 0.1 fab(y[i])
*/
void setDefaultDeltaBoundsMagnitudes();
//! Adjust the step minimums
void adjustUpStepMinimums();
//! Set the delta Bounds magnitudes by hand
/*!
* @param deltaBoundsMagnitudes set the deltaBoundsMagnitude vector
*/
void setDeltaBoundsMagnitudes(const doublereal * const deltaBoundsMagnitudes);
protected:
//! Readjust the trust region vectors
/*!
* 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
* correct length. We change the trustDelta_ values regularly
*
* The trust region calculate is based on
*
* || delta_x dot 1/trustDeltaX_ || <= trustDelta_
*
*/
void readjustTrustVector();
//! Fill a dogleg solution step vector
/*!
* Previously, we have filled up deltaX_Newton_[], deltaX_CP_[], and Nuu_, so that
* this routine is straightforward.
*
* @param leg Leg of the dog leg you are on (0, 1, or 2)
* @param alpha Relative length along the dog length that you are on.
* @param deltaX Vector to be filled up
*/
void fillDogLegStep(int leg, doublereal alpha, std::vector<doublereal> & deltaX) const;
//! Calculate the trust distance of a step in the solution variables
/*!
* The trust distance is defined as the length of the step according to the norm wrt to the trust region.
* We calculate the trust distance by the following method.
*
* trustDist = || delta_x dot 1/trustDeltaX_ ||
*
* @param deltaX Current value of deltaX
*/
doublereal calcTrustDistance(std::vector<doublereal> const & deltaX) const;
public:
//! Bound the step
/*!
*
* Return the factor by which the undamped Newton step 'step0'
* must be multiplied in order to keep all solution components in
* all domains between their specified lower and upper bounds.
* Other bounds may be applied here as well.
*
* Currently the bounds are hard coded into this routine:
*
* Minimum value for all variables: - 0.01 * m_ewt[i]
* Maximum value = none.
*
* Thus, this means that all solution components are expected
* to be numerical greater than zero in the limit of time step
* truncation errors going to zero.
*
* Delta bounds: The idea behind these is that the Jacobian
* couldn't possibly be representative if the
* variable is changed by a lot. (true for
* nonlinear systems, false for linear systems)
* Maximum increase in variable in any one newton iteration:
* factor of 2
* Maximum decrease in variable in any one newton iteration:
* factor of 5
*
* @param y Current solution value of the old step
* @param step0 Proposed step change in the solution
*
* @return Returns the damping factor determined by the bounds calculation
*/
doublereal boundStep(const doublereal * const y, const doublereal * const step0);
//! Set bounds constraints for all variables in the problem
/*!
*
* @param y_low_bounds Vector of lower bounds
* @param y_high_bounds Vector of high bounds
*/
void setBoundsConstraints(const doublereal * const y_low_bounds,
const doublereal * const y_high_bounds);
//! Return an editable vector of the low bounds constraints
std::vector<doublereal> & lowBoundsConstraintVector();
//! Return an editable vector of the high bounds constraints
std::vector<doublereal> & highBoundsConstraintVector();
//! 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 y_curr current value of the solution
* @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) const;
//! Function called to evaluate the jacobian matrix and the current
//! residual vector at the current time step
/*!
*
*
* @param J Jacobian matrix to be filled in
* @param f Right hand side. This routine returns the current
* value of the rhs (output), so that it does
* not have to be computed again.
* @param time_curr Current time
* @param CJ inverse of the value of deltaT
* @param y value of the solution vector
* @param ydot value of the time derivative of the solution vector
* @param num_newt_its Number of newton iterations
*
* @return Returns a flag to indicate that operation is successful.
* 1 Means a successful operation
* 0 Means an unsuccessful operation
*/
int beuler_jac(GeneralMatrix &J, doublereal * const f,
doublereal time_curr, doublereal CJ, doublereal * const y,
doublereal * const ydot, int num_newt_its);
//! Apply a filtering process to the step
/*!
* @param timeCurrent Current value of the time
* @param ybase current value of the solution
* @param step0 Proposed step change in the solution
*
* @return Returns the norm of the value of the amount filtered
*/
doublereal filterNewStep(const doublereal timeCurrent, const doublereal * const ybase, doublereal * const step0);
//! Apply a filter to the solution
/*!
* @param timeCurrent Current value of the time
* @param y_current current value of the solution
* @param ydot_current Current value of the solution derivative.
*
* @return Returns the norm of the value of the amount filtered
*/
doublereal filterNewSolution(const doublereal timeCurrent, doublereal * const y_current,
doublereal * const ydot_current);
//! Return the factor by which the undamped Newton step 'step0'
//! must be multiplied in order to keep the update within the bounds of an accurate jacobian.
/*!
*
* The idea behind these is that the Jacobian couldn't possibly be representative, if the
* variable is changed by a lot. (true for nonlinear systems, false for linear systems)
* Maximum increase in variable in any one newton iteration:
* factor of 1.5
* Maximum decrease in variable in any one newton iteration:
* factor of 2
*
* @param y Initial value of the solution vector
* @param step0 initial proposed step size
*
* @return returns the damping factor
*/
doublereal deltaBoundStep(const doublereal * const y, const doublereal * const step0);
//! Find a damping coefficient through a look-ahead mechanism
/*!
* On entry, step_1 must contain an undamped Newton step for the
* solution y_n_curr. This method attempts to find a damping coefficient
* such that all components stay in bounds, and the next
* undamped step would have a norm smaller than
* that of step_1. If successful, the new solution after taking the
* damped step is returned in y_n_1, and the undamped step at y_n_1 is
* returned in step_2.
*
* @param time_curr Current physical time
* @param y_n_curr Base value of the solution before any steps
* are taken
* @param ydot_n_curr Base value of the time derivative of the
* solution
* @param step_1 Initial step suggested.
* @param y_n_1 Value of y1, the suggested solution after damping
* @param ydot_n_1 Value of the time derivative of the solution at y_n_1
* @param step_2 Value of the step change from y_n_1 to y_n_2
* @param stepNorm_2 norm of the step change in going from y_n_1 to y_n_2
* @param jac Jacobian
* @param writetitle Write a title line
* @param num_backtracks Number of backtracks taken
*
* @return returns an integer indicating what happened.
*/
int dampStep(const doublereal time_curr, const doublereal * const y_n_curr,
const doublereal * const ydot_n_curr, doublereal * const step_1,
doublereal * const y_n_1, doublereal * const ydot_n_1, doublereal * step_2,
doublereal & stepNorm_2, GeneralMatrix& jac, bool writetitle,
int& num_backtracks);
//! Find the solution to F(X) = 0 by damped Newton iteration.
/*!
* On
* entry, x0 contains an initial estimate of the solution. On
* successful return, x1 contains the converged solution.
*
* SolnType = TRANSIENT -> we will assume we are relaxing a transient
* equation system for now. Will make it more general later,
* if an application comes up.
*
* @param SolnType Solution type
* @param y_comm Initial value of the solution. On return this is the converged
* value of the solution
* @param ydot_comm Initial value of the solution derivative. On return this is the
* converged value of the solution derivative.
* @param CJ Inverse of the value of deltaT
* @param time_curr Current value of the time
* @param jac Matrix that will be used to store the jacobian
* @param num_newt_its Number of newton iterations taken
* @param num_linear_solves Number of linear solves taken
* @param num_backtracks Number of backtracking steps taken
* @param loglevelInput Input log level determines the amount of printing.
*
*
* @return A positive value indicates a successful convergence
* -1 Failed convergence
*/
int solve_nonlinear_problem(int SolnType, doublereal * const y_comm, doublereal * const ydot_comm, doublereal CJ,
doublereal time_curr, GeneralMatrix & jac, int &num_newt_its,
int &num_linear_solves, int &num_backtracks, int loglevelInput);
//! Set the values for the previous time step
/*!
* We set the values for the previous time step here. These are used in the nonlinear
* solve because they affect the calculation of ydot.
*
* @param y_nm1 Value of the solution vector at the previous time step
* @param ydot_nm1 Value of the solution vector derivative at the previous time step
*/
virtual void
setPreviousTimeStep(const std::vector<doublereal>& y_nm1, const std::vector<doublereal>& ydot_nm1);
private:
//! Set the column scales
void calcColumnScales();
public:
//! Set the column scaling that are used for the inversion of the matrix
/*!
* There are three ways to do this.
*
* The first method is to set the bool useColScaling to true, leaving the scaling factors unset.
* Then, the column scales will be set to the solution error weighting factors. This has the
* effect of ensuring that all delta variables will have the same order of magnitude at convergence
* end.
*
* The second way is the explicitly set the column factors in the second argument of this function call.
*
* The final way to input the scales is to override the ResidJacEval member function call,
*
* calcSolnScales(double time_n, const double *m_y_n_curr, const double *m_y_nm1, double *m_colScales)
*
* Overriding this function call will trump all other ways to specify the column scaling factors.
*
* @param useColScaling Turn this on if you want to use column scaling in the calculations
* @param scaleFactors A vector of doubles that specifies the column factors.
*/
void setColumnScaling(bool useColScaling, const double * const scaleFactors = 0);
//! Set the rowscaling that are used for the inversion of the matrix
/*!
* Row scaling is set here. Right now the row scaling is set internally in the code.
*
* @param useRowScaling Turn row scaling on or off.
*/
void setRowScaling(bool useRowScaling);
//! Scale the matrix
/*!
* @param jac Jacobian
* @param y_comm Current value of the solution vector
* @param ydot_comm Current value of the time derivative of the solution vector
* @param time_curr current value of the time
* @param num_newt_its Current value of the number of newt its
*/
void scaleMatrix(GeneralMatrix& jac, doublereal * const y_comm, doublereal * const ydot_comm,
doublereal time_curr, int num_newt_its);
//! Print solution norm contribution
/*!
* Prints out the most important entries to the update to the solution vector for the current step
*
* @param step_1 Raw update vector for the current nonlinear step
* @param stepNorm_1 Norm of the vector step_1
* @param step_2 Raw update vector for the next solution value based on the old matrix
* @param stepNorm_2 Norm of the vector step_2
* @param title title of the printout
* @param y_n_curr Old value of the solution
* @param y_n_1 New value of the solution after damping corrections
* @param damp Value of the damping factor
* @param num_entries Number of entries to print out
*/
void
print_solnDelta_norm_contrib(const doublereal * const step_1, const char * const stepNorm_1,
const doublereal * const step_2, const char * const stepNorm_2,
const char * const title, const doublereal * const y_n_curr,
const doublereal * const y_n_1, doublereal damp, size_t num_entries);
//! Compute the Residual Weights
/*!
* The residual weights are defined here to be equal to the inverse of the row scaling factors used to
* row scale the matrix, after column scaling is used. They are multiplied by 10-3 because the column
* weights are also multiplied by that same quantity.
*
* The basic idea is that a change in the solution vector on the order of the convergence tolerance
* multiplied by [RJC] which is of order one after row scaling should give you the relative weight
* of the row. Values of the residual for that row can then be normalized by the value of this weight.
* When the tolerance in delta x is achieved, the tolerance in the residual is also achieved.
*/
void computeResidWts();
//! Return the residual weights
/*!
* @param residWts Vector of length neq_
*/
void getResidWts(doublereal * const residWts) const;
//! Check to see if the nonlinear problem has converged
/*!
*
* @param dampCode Code from the damping routine
* @param s1 Value of the norm of the step change
*
* @return integer is returned. If positive, then the problem has converged
* 1 Successful step was taken: Next step's norm is less than 1.0.
* The final residual norm is less than 1.0.
* 2 Successful step: Next step's norm is less than 0.8.
* This step's norm is less than 1.0.
* The residual norm can be anything.
* 3 Success: The final residual is less than 1.0
* The predicted deltaSoln is below 1.0.
* 0 Not converged yet
*/
int convergenceCheck(int dampCode, doublereal s1);
//! Set the absolute tolerances for the solution variables
/*!
* Set the absolute tolerances used in the calculation
*
* @param atol Vector of length neq_ that contains the tolerances to be used for the solution variables
*/
void setAtol(const doublereal * const atol);
//! Set the relative tolerances for the solution variables
/*!
* Set the relative tolerances used in the calculation
*
* @param rtol single double
*/
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.
*
* With this routine the user can override or add to the residual weighting norm evaluation by specifying
* their own vector of residual absolute and relative tolerances.
*
* 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 calculated 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
*/
void setMaxNewtIts(const int maxNewtIts);
//! 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
//! of the nonlinear problem expects a minimum along the descent direction.
/*!
* @param jac Jacobian matrix: must be unfactored.
*
* @return Returns the norm of the solution update
*/
doublereal doCauchyPointSolve(GeneralMatrix& jac);
//! This is a utility routine that can be used to print out the rates of the initial residual decline
/*!
* The residual**2 decline for various directions is printed out. The rate of decline of the
* square of the residuals multiplied by the number of equations along each direction is printed out
* This quantity can be directly related to the theory, and may be calculated from derivatives at the
* original point.
*
* ( (r)**2 * neq - (r0)**2 * neq ) / distance
*
* What's printed out:
*
* The theoretical linearized residual decline
* The actual residual decline in the steepest descent direction determined by numerical differencing
* The actual residual decline in the newton direction determined by numerical differencing
*
* This routine doesn't need to be called for the solution of the nonlinear problem.
*
* @param time_curr Current time
* @param ydot0 INPUT Current value of the derivative of the solution vector
* @param ydot1 INPUT Time derivatives of solution at the conditions which are evaluated for success
* @param numTrials OUTPUT Counter for the number of residual evaluations
*/
void descentComparison(doublereal time_curr ,doublereal * ydot0, doublereal * ydot1, int &numTrials);
//! Setup the parameters for the double dog leg
/*!
* The calls to the doCauchySolve() and doNewtonSolve() routines are done at the main level. This routine comes
* after those calls. We calculate the point Nuu_ here, the distances of the dog-legs,
* and the norms of the CP and Newton points in terms of the trust vectors.
*/
void setupDoubleDogleg();
//! Change the global lambda coordinate into the (leg,alpha) coordinate for the double dogleg
/*!
* @param lambda Global value of the distance along the double dogleg
* @param alpha relative value along the particular leg
*
* @return Returns the leg number ( 0, 1, or 2).
*/
int lambdaToLeg(const doublereal lambda, doublereal &alpha) const;
//! Given a trust distance, this routine calculates the intersection of the this distance with the
//! double dogleg curve
/*!
* @param trustVal (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(doublereal trustVal, doublereal &lambda, doublereal &alpha) const;
//! Initialize the size of the trust vector.
/*!
* The algorithm we use is to set it equal to the length of the Distance to the Cauchy point.
*/
void initializeTrustRegion();
//! Set Trust region initialization strategy
/*!
* The default is use method 2 with a factor of 1.
* Then, on subsequent invocations of solve_nonlinear_problem() the strategy flips to method 0.
*
* @param method Method to set the strategy
* 0 No strategy - Use the previous strategy
* 1 Factor of the solution error weights
* 2 Factor of the first Cauchy Point distance
* 3 Factor of the first Newton step distance
*
* @param factor Factor to use in combination with the method
*
*/
void setTrustRegionInitializationMethod(int method, doublereal factor);
//! Damp using the dog leg approach
/*!
*
* @param time_curr INPUT Current value of the time
* @param y_n_curr INPUT Current value of the solution vector
* @param ydot_n_curr INPUT Current value of the derivative of the solution vector
* @param step_1 INPUT First trial step for the first iteration
* @param y_n_1 INPUT First trial value of the solution vector
* @param ydot_n_1 INPUT First trial value of the derivative of the solution vector
* @param stepNorm_1 OUTPUT Norm of the vector step_1
* @param stepNorm_2 OUTPUT Estimated norm of the vector step_2
* @param jac INPUT jacobian
* @param num_backtracks OUTPUT number of backtracks taken in the current damping step
*
* @return 1 Successful step was taken. The predicted residual norm is less than one
* 2 Successful step: Next step's norm is less than 0.8
* 3 Success: The final residual is less than 1.0
* A predicted deltaSoln1 is not produced however. s1 is estimated.
* 4 Success: The final residual is less than the residual
* from the previous step.
* A predicted deltaSoln1 is not produced however. s1 is estimated.
* 0 Uncertain Success: s1 is about the same as s0
* -2 Unsuccessful step.
*/
int dampDogLeg(const doublereal time_curr, const doublereal* y_n_curr,
const doublereal *ydot_n_curr, std::vector<doublereal> & step_1,
doublereal* const y_n_1, doublereal* const ydot_n_1,
doublereal& stepNorm_1, doublereal& stepNorm_2, GeneralMatrix& jac, int& num_backtracks);
//! Decide whether the current step is acceptable and adjust the trust region size
/*!
* This is an extension of algorithm 6.4.5 of Dennis and Schnabel.
*
* Here we decide whether to accept the current step
* At the end of the calculation a new estimate of the trust region is calculated
*
* @param time_curr INPUT Current value of the time
* @param leg INPUT Leg of the dogleg that we are on
* @param alpha INPUT Distance down that leg that we are on
* @param y_n_curr INPUT Current value of the solution vector
* @param ydot_n_curr INPUT Current value of the derivative of the solution vector
* @param step_1 INPUT Trial step
* @param y_n_1 OUTPUT Solution values at the conditions which are evaluated for success
* @param ydot_n_1 OUTPUT Time derivatives of solution at the conditions which are evaluated for success
* @param trustDeltaOld INPUT Value of the trust length at the old conditions
*
*
* @return This function returns a code which indicates whether the step will be accepted or not.
* 3 Step passed with flying colors. Try redoing the calculation with a bigger trust region.
* 2 Step didn't pass deltaF requirement. Decrease the size of the next trust region for a retry and return
* 0 The step passed.
* -1 The step size is now too small (||d || < 0.1). A really small step isn't decreasing the function.
* This is an error condition.
* -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 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);
//! Calculated the expected residual along the double dogleg curve.
/*!
* @param leg 0, 1, or 2 representing the curves of the dogleg
* @param alpha Relative distance along the particular curve.
*
* @return Returns the expected value of the residual at that point according to the quadratic model.
* The residual at the newton point will always be zero.
*/
doublereal expectedResidLeg(int leg, doublereal alpha) const;
//! 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
* @param legBest OUTPUT leg of the dogleg that gives the lowest residual
* @param alphaBest OUTPUT distance along dogleg for best result.
*/
void residualComparisonLeg(const doublereal time_curr, const doublereal * const ydot0, int & legBest,
doublereal & alphaBest) const;
//! Set the print level from the nonlinear solver
/*!
*
* 0 -> absolutely nothing is printed for a single time step.
* 1 -> One line summary per solve_nonlinear call
* 2 -> short description, points of interest: Table of nonlinear solve - one line per iteration
* 3 -> Table is included -> More printing per nonlinear iteration (default) that occurs during the table
* 4 -> Summaries of the nonlinear solve iteration as they are occurring -> table no longer printed
* 5 -> Algorithm information on the nonlinear iterates are printed out
* 6 -> Additional info on the nonlinear iterates are printed out
* 7 -> Additional info on the linear solve is printed out.
* 8 -> Info on a per iterate of the linear solve is printed out.
*
* @param printLvl integer value
*/
void setPrintLvl(int printLvl);
//! Parameter to turn on solution solver schemes
/*!
* @param doDogLeg Parameter to turn on the double dog leg scheme
* Default is to always use a damping scheme in the Newton Direction.
* When this is nonzero, a model trust region approach is used using a double dog leg
* with the steepest descent direction used for small step sizes.
*
* @param doAffineSolve Parameter to turn on or off the solution of the system using a Hessian
* if the matrix has a bad condition number.
*/
void setSolverScheme(int doDogLeg, int doAffineSolve);
/*
* -----------------------------------------------------------------------------------------------------------------
* MEMBER DATA
* ------------------------------------------------------------------------------------------------
*/
private:
//! Pointer to the residual and jacobian evaluator for the
//! function
/*!
* See ResidJacEval.h for an evaluator.
*/
ResidJacEval *m_func;
//! Solution type
int solnType_;
//! Local copy of the number of equations
size_t neq_;
//! Soln error weights
std::vector<doublereal> m_ewt;
//! Boolean indicating whether a manual delta bounds has been input.
int m_manualDeltaStepSet;
//! Soln Delta bounds magnitudes
std::vector<doublereal> m_deltaStepMinimum;
//! Value of the delta step magnitudes
std::vector<doublereal> m_deltaStepMaximum;
//! Vector containing the current solution vector within the nonlinear solver
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
std::vector<doublereal> m_y_nm1;
//! Vector containing the solution derivative at the previous time step
std::vector<doublereal> m_ydot_nm1;
//! Vector containing the solution at the new point that is to be considered
std::vector<doublereal> m_y_n_trial;
//! Value of the solution time derivative at the new point that is to be considered
std::vector<doublereal> m_ydot_trial;
//! Value of the step to be taken in the solution
std::vector<doublereal> m_step_1;
//! Vector of column scaling factors
std::vector<doublereal> m_colScales;
//! 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 of the rows of the jacobian.
*/
std::vector<doublereal> m_rowScales;
//! Weights for normalizing the values of the residuals
/*!
* 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;
//! Value of the residual for the nonlinear problem
mutable std::vector<doublereal> m_resid;
//! Workspace of length neq_
mutable std::vector<doublereal> m_wksp;
//! Workspace of length neq_
mutable std::vector<doublereal> m_wksp_2;
/*****************************************************************************************
* INTERNAL WEIGHTS FOR TAKING SOLUTION NORMS
******************************************************************************************/
//! Vector of residual weights
/*!
* These are used to establish useful and informative weighted norms of the residual vector.
*/
std::vector<doublereal> m_residWts;
//! Norm of the residual at the start of each nonlinear iteration
doublereal m_normResid_0;
//! Norm of the residual after it has been bounded
doublereal m_normResid_Bound;
//! Norm of the residual at the end of the first leg of the current iteration
doublereal m_normResid_1;
//! Norm of the residual at the end of the first leg of the current iteration
doublereal m_normResid_full;
//! Norm of the solution update created by the iteration in its raw, undamped form, using the solution norm
doublereal m_normDeltaSoln_Newton;
//! Norm of the distance to the cauchy point using the solution norm
doublereal m_normDeltaSoln_CP;
//! Norm of the residual for a trial calculation which may or may not be used
doublereal m_normResidTrial;
//! Vector of the norm
doublereal m_normResidPoints[15];
//! Boolean indicating whether we should scale the residual
mutable bool m_resid_scaled;
/*****************************************************************************************
* INTERNAL BOUNDARY INFO FOR SOLUTIONS
*****************************************************************************************/
//! Bounds vector for each species
std::vector<doublereal> m_y_high_bounds;
//! Lower bounds vector for each species
std::vector<doublereal> m_y_low_bounds;
//! Damping factor imposed by hard bounds and by delta bounds
doublereal m_dampBound;
//! Additional damping factor due to bounds on the residual and solution norms
doublereal m_dampRes;
//! Delta t for the current step
doublereal delta_t_n;
//! Counter for the total number of function evaluations
mutable int m_nfe;
/***********************************************************************************************
* MATRIX INFORMATION
**************************************************************************************/
//! The type of column scaling used in the matrix inversion of the problem
/*!
* If 1 then colScaling = m_ewt[]
* If 2 then colScaling = user set
* if 0 then colScaling = 1.0
*/
int m_colScaling;
//! int indicating whether row scaling is turned on (1) or not (0)
int m_rowScaling;
//! Total number of linear solves taken by the solver object
int m_numTotalLinearSolves;
//! Number of local linear solves done during the current iteration
int m_numLocalLinearSolves;
//! Total number of newton iterations
int m_numTotalNewtIts;
public:
//! Minimum number of newton iterations to use
int m_min_newt_its;
private:
//! Maximum number of newton iterations
int maxNewtIts_;
//! Jacobian formation method
/*!
* 1 = numerical (default)
* 2 = analytical
*/
int m_jacFormMethod;
//! Number of Jacobian evaluations
int m_nJacEval;
//! Current system time
/*!
* Note, we assume even for steady state problems that the residual
* is a function of a system time.
*/
doublereal time_n;
//! Boolean indicating matrix conditioning
int m_matrixConditioning;
//! Order of the time step method = 1
int m_order;
//! value of the relative tolerance to use in solving the equation set
doublereal rtol_;
//! Base value of the absolute tolerance
doublereal atolBase_;
//! absolute tolerance in the solution unknown
/*!
* This is used to evaluating the weighting factor
*/
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.
* 1 -> One line summary per solve_nonlinear call
* 2 -> short description, points of interest: Table of nonlinear solve - one line per iteration
* 3 -> Table is included -> More printing per nonlinear iteration (default) that occurs during the table
* 4 -> Summaries of the nonlinear solve iteration as they are occurring -> table no longer printed
* Base_ShowSolution Residual called for residual printing at the end of convergence.
* 5 -> Algorithm information on the nonlinear iterates are printed out
* 6 -> Additional info on the nonlinear iterates are printed out
* Base_ShowSolution Residual called for residual printing at the end of each step.
* 7 -> Additional info on the linear solve is printed out.
* 8 -> Info on a per iterate of the linear solve is printed out.
*/
int m_print_flag;
//! Scale factor for turning residual norms into solution norms
doublereal m_ScaleSolnNormToResNorm;
//! Copy of the jacobian that doesn't get overwritten when the inverse is determined
/*!
* The jacobian stored here is the raw matrix, before any row or column scaling is carried out
*/
Cantera::GeneralMatrix * jacCopyPtr_;
//! Hessian
Cantera::GeneralMatrix * HessianPtr_;
/*********************************************************************************************
* VARIABLES ASSOCIATED WITH STEPS AND ASSOCIATED DOUBLE DOGLEG PARAMETERS
*********************************************************************************************/
//! Steepest descent direction. This is also the distance to the Cauchy Point
std::vector<doublereal> deltaX_CP_;
//! Newton Step - This is the newton step determined from the straight Jacobian
/*
* Newton step for the current step only
*/
std::vector<doublereal> deltaX_Newton_;
//! Expected value of the residual norm at the Cauchy point if the quadratic model
//! were valid
doublereal residNorm2Cauchy_;
//! Current leg
int dogLegID_;
//! Current Alpha param along the leg
doublereal dogLegAlpha_;
//! Residual dot Jd norm
/*!
* This is equal to R_hat dot J_hat d_y_descent
*/
doublereal RJd_norm_;
//! Value of lambdaStar_ which is used to calculate the Cauchy point
doublereal lambdaStar_;
//! Jacobian times the steepest descent direction in the normalized coordinates.
/*!
* This is equal to [ Jhat d^y_{descent} ] in the notes, Eqn. 18.
*/
std::vector<doublereal> Jd_;
//! Vector of trust region values.
std::vector<doublereal> deltaX_trust_;
//! Current norm of the vector deltaX_trust_ in terms of the solution norm
mutable doublereal norm_deltaX_trust_;
//! Current value of trust radius. This is used with deltaX_trust_ to
//! calculate the max step size.
doublereal trustDelta_;
//! Method for handling the trust region initialization
/*!
* Then, on subsequent invocations of solve_nonlinear_problem() the strategy flips to method 0.
*
* method Method to set the strategy
* 0 No strategy - Use the previous strategy
* 1 Factor of the solution error weights
* 2 Factor of the first Cauchy Point distance
* 3 Factor of the first Newton step distance
*/
int trustRegionInitializationMethod_;
//! Factor used to set the initial trust region
doublereal trustRegionInitializationFactor_;
//! Relative distance down the Newton step that the second dogleg starts
doublereal Nuu_;
//! Distance of the zeroeth leg of the dogleg in terms of the solution norm
doublereal dist_R0_;
//! Distance of the first leg of the dogleg in terms of the solution norm
doublereal dist_R1_;
//! Distance of the second leg of the dogleg in terms of the solution norm
doublereal dist_R2_;
//! Distance of the sum of all legs of the doglegs in terms of the solution norm
doublereal dist_Total_;
//! Dot product of the Jd_ variable defined above with itself.
doublereal JdJd_norm_;
//! Norm of the Newton Step wrt trust region
doublereal normTrust_Newton_;
//! Norm of the Cauchy Step direction wrt trust region
doublereal normTrust_CP_;
//! General toggle for turning on dog leg damping.
int doDogLeg_;
//! General toggle for turning on Affine solve with Hessian
int doAffineSolve_;
//! Condition number of the matrix
doublereal m_conditionNumber;
//! Factor indicating how much trust region has been changed this iteration - output variable
doublereal CurrentTrustFactor_;
//! Factor indicating how much trust region has been changed next iteration - output variable
doublereal NextTrustFactor_;
//! Boolean indicating that the residual weights have been reevaluated this iteration - output variable
bool ResidWtsReevaluated_;
//! Expected DResid_dS for the steepest descent path - output variable
doublereal ResidDecreaseSDExp_;
//! Actual DResid_dS for the steepest descent path - output variable
doublereal ResidDecreaseSD_;
//! Expected DResid_dS for the Newton path - output variable
doublereal ResidDecreaseNewtExp_;
//! Actual DResid_dS for the newton path - output variable
doublereal ResidDecreaseNewt_;
/*******************************************************************************************
* STATIC VARIABLES
*****************************************************************************************/
public:
//! Turn off printing of time
/*!
* Necessary to do for test suites
*/
static bool s_TurnOffTiming;
//! Turn on or off printing of the Jacobian
static bool s_print_NumJac;
//! Turn on extra printing of dogleg information
static bool s_print_DogLeg;
//! Turn on solving both the Newton and Hessian systems and comparing the results
/*!
* This is off by default
*/
static bool s_doBothSolvesAndCompare;
//! This toggle turns off the use of the Hessian when it is warranted by the condition number.
/*!
* This is a debugging option.
*/
static bool s_alwaysAssumeNewtonGood;
};
}
#endif
Reference in a new issue View git blame Copy permalink