diff --git a/include/cantera/numerics/DAE_Solver.h b/include/cantera/numerics/DAE_Solver.h index 465298418..99b89110d 100644 --- a/include/cantera/numerics/DAE_Solver.h +++ b/include/cantera/numerics/DAE_Solver.h @@ -169,7 +169,7 @@ public: * derivatives of all differential variables. * @param y Calculated value of the solution vector after the procedure ends * @param yp Calculated value of the solution derivative after the procedure - * @param The first value of t at which a soluton will be + * @param tout The first value of t at which a soluton will be * requested (from IDASolve). (This is needed here to * determine the direction of integration and rough scale * in the independent variable t. diff --git a/include/cantera/numerics/NonlinearSolver.h b/include/cantera/numerics/NonlinearSolver.h index a28445813..ffe03656b 100644 --- a/include/cantera/numerics/NonlinearSolver.h +++ b/include/cantera/numerics/NonlinearSolver.h @@ -82,12 +82,10 @@ namespace Cantera * * 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 @@ -113,30 +111,19 @@ namespace Cantera * * 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 /*! @@ -144,23 +131,13 @@ public: */ NonlinearSolver(ResidJacEval* func); - //!Copy Constructor for the %ThermoPhase object. - /*! - * @param right Item to be copied - */ + //!Copy Constructor 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 @@ -176,14 +153,13 @@ public: */ void createSolnWeights(const doublereal* const y); - //! L2 norm of the delta of the solution vector + //! 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. + * 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 @@ -201,16 +177,14 @@ public: * 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. + * 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 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, @@ -218,19 +192,18 @@ public: //! Compute the current residual /*! - * The current value of the residual is stored in the internal work array m_resid, which is defined - * as mutable + * 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 + * @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 + * 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, @@ -238,11 +211,9 @@ public: //! 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. + * 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, @@ -256,15 +227,15 @@ public: * @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. + * @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 ( + //! 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 * @@ -290,7 +261,6 @@ public: * --------------- * 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. */ @@ -323,7 +293,6 @@ public: 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 @@ -333,7 +302,6 @@ protected: * The trust region calculate is based on * * || delta_x dot 1/trustDeltaX_ || <= trustDelta_ - * */ void readjustTrustVector(); @@ -359,12 +327,9 @@ protected: */ doublereal calcTrustDistance(std::vector 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. @@ -397,7 +362,6 @@ public: //! 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 */ @@ -424,8 +388,6 @@ public: //! 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 @@ -468,13 +430,10 @@ public: //! 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 + * 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 @@ -507,7 +466,17 @@ public: * @param writetitle Write a title line * @param num_backtracks Number of backtracks taken * - * @return returns an integer indicating what happened. + * @return 1 Successful step was taken: Next step was less than previous step. + * s1 is calculated + * 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 + * NSOLN_RETN_FAIL_DAMPSTEP + * Unsuccessful step. We can not find a damping factor that is suitable. */ int dampStep(const doublereal time_curr, const doublereal* const y_n_curr, const doublereal* const ydot_n_curr, doublereal* const step_1, @@ -517,8 +486,7 @@ public: //! Find the solution to F(X) = 0 by damped Newton iteration. /*! - * On - * entry, x0 contains an initial estimate of the solution. On + * 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 @@ -538,7 +506,6 @@ public: * @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 */ @@ -558,7 +525,7 @@ public: setPreviousTimeStep(const std::vector& y_nm1, const std::vector& ydot_nm1); private: - //! Set the column scales + //! Set the column scaling vector at the current time void calcColumnScales(); public: @@ -585,7 +552,6 @@ public: */ 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. @@ -644,11 +610,8 @@ public: */ 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 * @@ -702,6 +665,9 @@ public: * * residWeightNorm[i] = residAtol[i] + residRtol * m_rowWtScales[i] / neq * + * @param residRtol scalar residual relative tolerance + * @param residAtol vector of residual absolute tolerances + * * @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. @@ -710,7 +676,7 @@ public: * 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); + void setResidualTols(double residRtol, double* residAtol, int residNormHandling = 2); //! Set the value of the maximum # of newton iterations /*! @@ -760,7 +726,6 @@ public: */ 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 @@ -806,14 +771,11 @@ public: * 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 @@ -857,7 +819,6 @@ public: * @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 @@ -896,7 +857,6 @@ public: //! 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 @@ -923,8 +883,6 @@ public: */ void setSolverScheme(int doDogLeg, int doAffineSolve); - - /* * ----------------------------------------------------------------------------------------------------------------- * MEMBER DATA diff --git a/include/cantera/numerics/ResidEval.h b/include/cantera/numerics/ResidEval.h index fc1b4a78c..06650fd0e 100644 --- a/include/cantera/numerics/ResidEval.h +++ b/include/cantera/numerics/ResidEval.h @@ -78,10 +78,8 @@ public: return (m_alg[k] == 1); } - /** - * Evaluate the residual function. Called by the - * integrator. + * Evaluate the residual function. Called by the integrator. * @param t time. (input) * @param y solution vector. (input) * @param ydot rate of change of solution vector. (input) @@ -109,11 +107,17 @@ public: return eval(t, y, DATA_PTR(ydot), r); } - /** - * Fill the solution and derivative vectors with the initial - * conditions at initial time t0. - * @return 1 Everything is fine - * 0 or neg Something went wrong + //! Fill in the initial conditions + /*! + * Values for both the solution and the value of ydot may be provided. + * + * @param[in] t0 Time + * @param[out] y Solution vector + * @param[out] ydot Rate of change of solution vector. + * + * @return Returns a flag to indicate that operation is successful. + * 1 Means a successful operation + * -0 or neg value Means an unsuccessful operation */ virtual int getInitialConditions(const doublereal t0, doublereal* const y, doublereal* const ydot) { diff --git a/include/cantera/numerics/ResidJacEval.h b/include/cantera/numerics/ResidJacEval.h index e6ea18ad6..19d10ef54 100644 --- a/include/cantera/numerics/ResidJacEval.h +++ b/include/cantera/numerics/ResidJacEval.h @@ -47,41 +47,27 @@ enum ResidEval_Type_Enum { * A class for full (non-sparse dense matrices with Fortran-compatible data storage. * The class adds support for identifying what types of calls are made to the residual * evaluator by adding the ResidEval_Type_Enum class. - * - * */ class ResidJacEval : public ResidEval { - public: - //!Default constructor /*! * @param atol Initial value of the global tolerance (defaults to 1.0E-13) */ ResidJacEval(doublereal atol = 1.0e-13); - //!Copy Constructor for the %ResidJacEval object - /*! - * @param right Item to be copied - */ + //!Copy Constructor ResidJacEval(const ResidJacEval& right); //! Assignment operator - /*! - * This is NOT a virtual function. - * - * @param right Reference to %ResidJacEval object to be copied into the - * current one. - */ ResidJacEval& operator=(const ResidJacEval& right); - //! Duplication routine for objects which inherit from - //! residJacEval + //! Duplication routine for objects derived from residJacEval /*! - * This virtual routine can be used to duplicate %ResidJacEval objects - * inherited from %ResidJacEval even if the application only has - * a pointer to %ResidJacEval to work with. + * This virtual routine can be used to duplicate objects which inherit + * from ResidJacEval even if the application only has a pointer to + * ResidJacEval to work with. * * These routines are basically wrappers around the derived copy * constructor. @@ -116,33 +102,10 @@ public: const int id_x = -1, const doublereal delta_x = 0.0); - - /** - * Evaluate the residual function. Called by the - * integrator. - * @param t time. (input) - * @param y solution vector. (input) - * @param ydot rate of change of solution vector. (input) - * @param r residual vector (output) - */ virtual int eval(const doublereal t, const doublereal* const y, const doublereal* const ydot, doublereal* const r); - - - //! Fill in the initial conditions - /*! - * Values for both the solution and the value of ydot may be provided. - * - * @param t0 Time (input) - * @param y Solution vector (output) - * @param ydot Rate of change of solution vector. (output) - * - * @return Returns a flag to indicate that operation is successful. - * 1 Means a successful operation - * -0 or neg value Means an unsuccessful operation - */ virtual int getInitialConditions(const doublereal t0, doublereal* const y, doublereal* const ydot); //! Filter the solution predictions @@ -251,7 +214,6 @@ public: //! This function may be used to create output at various points in the execution of an application. /*! - * * @param ifunc identity of the call * 0 Initial call * 1 Called at the end of every successful time step @@ -324,7 +286,6 @@ public: const doublereal* const y, const doublereal* const ydot, GeneralMatrix& J, doublereal* const resid); - //! Calculate an analytical jacobian and the residual at the current time and values. /*! * Only called if the jacFormation method is set to analytical @@ -349,7 +310,6 @@ public: doublereal* const resid); protected: - //! constant value of atol doublereal m_atol; @@ -359,4 +319,3 @@ protected: } #endif - diff --git a/include/cantera/numerics/RootFind.h b/include/cantera/numerics/RootFind.h index 15958c098..3de0743a5 100644 --- a/include/cantera/numerics/RootFind.h +++ b/include/cantera/numerics/RootFind.h @@ -43,18 +43,11 @@ namespace Cantera //! This means that the rootfinder believes the solution is lower than xmin #define ROOTFIND_SOLNLOWERTHANXMIN -3 //! This means that the rootfinder believes the solution is higher than xmax -/*! - * - */ #define ROOTFIND_SOLNHIGHERTHANXMAX -4 //@} - - - //! Root finder for 1D problems /*! - * * The root finder solves a single nonlinear equation described below. * * \f[ @@ -85,13 +78,10 @@ namespace Cantera * * Supplying Hints with General Function Behavior Flags * - * - * * Stopping Criteria * * Specification of the Stopping Criteria * - * * Additional constraints * * Bounds Criteria For the Routine @@ -141,34 +131,21 @@ namespace Cantera */ class RootFind { - public: - //! Constructor for the object /*! - * * @param resid Pointer to the residual function to be used to calculate f(x) */ RootFind(ResidEval* resid); //! Copy constructor - /*! - * @param r object to be copied - */ RootFind(const RootFind& r); - //! Destructor. Deletes the integrator. ~RootFind(); //! Assignment operator - /*! - * @param right object to be copied - * - * @return returns a reference to the current object - */ RootFind& operator=(const RootFind& right); - private: //! Calculate a deltaX from an input value of x /*! @@ -216,7 +193,6 @@ private: bool theSame(doublereal x2, doublereal x1, doublereal factor = 1.0) const; public: - //! Using a line search method, find the root of a 1D function /*! * This routine solves the following equation. @@ -333,7 +309,6 @@ public: void printTable(); public: - //! Pointer to the residual function evaluator ResidEval* m_residFunc; @@ -356,7 +331,6 @@ public: doublereal m_maxstep; protected: - //! Print level /*! * 0 No printing of any kind @@ -368,7 +342,6 @@ protected: int printLvl; public: - //! Boolean to turn on the possibility of writing a log file. bool writeLogAllowed_; @@ -413,7 +386,6 @@ protected: //! Internal variable tracking f(x) of smallest x tried. doublereal fx_minTried_; - //! Structure containing the iteration history struct rfTable { //@{ @@ -428,7 +400,6 @@ protected: double deltaFConverged; double delX; - std::string reasoning; void clear() { @@ -453,16 +424,11 @@ protected: delX(-1.0E300), reasoning("") { }; - - //@} }; //! Vector of iteration histories std::vector rfHistory_; - - - }; } #endif diff --git a/include/cantera/numerics/solveProb.h b/include/cantera/numerics/solveProb.h index 9d1e81086..78f034ee6 100644 --- a/include/cantera/numerics/solveProb.h +++ b/include/cantera/numerics/solveProb.h @@ -55,13 +55,8 @@ const int SOLVEPROB_RESIDUAL = 2; const int SOLVEPROB_JACOBIAN = 3; const int SOLVEPROB_TRANSIENT = 4; - - - namespace Cantera { - - //! Method to solve a pseudo steady state of a nonlinear problem /*! * The following class handles the solution of a nonlinear problem. @@ -76,7 +71,6 @@ namespace Cantera * Res_ss(C) is the steady state residual to be solved. Res_td(C) is the * time dependent residual which leads to the steady state residual. * - * * Solution Method * * This routine is typically used within a residual calculation in a large code. @@ -84,11 +78,10 @@ namespace Cantera * work every time. Therefore, requirements demand that it be robust but also * efficient. * - * The solution methodology is largely determined by the ifunc<\TT> parameter, + * The solution methodology is largely determined by the ifunc parameter, * that is input to the solution object. This parameter may have the following * 4 values: * - * * 1: SOLVEPROB_INITIALIZE = This assumes that the initial guess supplied to the * routine is far from the correct one. Substantial * work plus transient time-stepping is to be expected @@ -127,7 +120,6 @@ namespace Cantera * Damping is based on a "delta damping" technique. The solution unknowns * are not allowed to vary too much between iterations. * - * * EXTRA_ACCURACY:A constant that is the ratio of the required update norm in * this Newton iteration compared to that in the nonlinear solver. * A value of 0.1 is used so surface species are safely overconverged. @@ -146,21 +138,11 @@ namespace Cantera */ class solveProb { - public: //! Constructor for the object - /*! - * @param surfChemPtr Pointer to the ImplicitSurfChem object that - * defines the surface problem to be solved. - * - * @param bulkFunc Integer representing how the bulk phases - * should be handled. Currently, only the - * default value of BULK_ETCH is supported. - */ solveProb(ResidEval* resid); - //! Destructor. Deletes the integrator. virtual ~solveProb(); private: @@ -172,7 +154,6 @@ private: solveProb& operator=(const solveProb& right); public: - //! Main routine that actually calculates the pseudo steady state //! of the surface problem /*! @@ -197,7 +178,7 @@ public: //! Report the current state of the solution /*! - * @param Report the solution vector for the nonlinear problem + * @param[out] CSoln solution vector for the nonlinear problem */ virtual void reportState(doublereal* const CSoln) const; @@ -210,20 +191,17 @@ public: */ virtual void setBounds(const doublereal botBounds[], const doublereal topBounds[]); - void setAtol(const doublereal atol[]); void setAtolConst(const doublereal atolconst); private: - - //! Printing routine that gets called at the start of every - //! invocation + //! Printing routine that gets called at the start of every invocation virtual void print_header(int ioflag, int ifunc, doublereal time_scale, doublereal reltol, doublereal netProdRate[]); #ifdef DEBUG_SOLVEPROB - + //! Prints out the residual and Jacobian virtual void printResJac(int ioflag, int neq, const Array2D& Jac, doublereal resid[], doublereal wtResid[], doublereal norm); #endif @@ -236,11 +214,7 @@ private: doublereal resid[], doublereal wtSpecies[], size_t dim, bool do_time); - //! Print a summary of the solution - /*! - * - */ virtual void printFinal(int ioflag, doublereal damp, size_t label_d, size_t label_t, doublereal inv_t, doublereal t_real, int iter, doublereal update_norm, doublereal resid_norm, @@ -264,7 +238,7 @@ private: * * @param netProdRateSolnSP Output variable. Net production rate * of all of the species in the solution vector. - * @param XMolSolnSP output variable. + * @param Csoln output variable. * Mole fraction of all of the species in the solution vector * @param label Output variable. Pointer to the value of the * species index (kindexSP) that is controlling @@ -286,7 +260,9 @@ private: //! Calculate the solution and residual weights /*! - * @param wtSpecies Weights to use for the soln unknowns. These + * Calculate the weighting factors for norms wrt both the species + * concentration unknowns and the residual unknowns. + * @param wtSpecies Weights to use for the soln unknowns. These * are in concentration units * @param wtResid Weights to sue for the residual unknowns. * @@ -313,6 +289,7 @@ private: //! Main Function evaluation /*! + * This calculates the net production rates of all species * * @param resid output Vector of residuals, length = m_neq * @param CSolnSP Vector of species concentrations, unknowns in the @@ -323,7 +300,7 @@ private: * @param deltaT Delta time for time dependent problem. */ virtual void fun_eval(doublereal* const resid, const doublereal* const CSolnSP, - const doublereal* const CSolnOldSP, const bool do_time, const doublereal deltaT); + const doublereal* const CSolnSPOld, const bool do_time, const doublereal deltaT); //! Main routine that calculates the current residual and Jacobian /*! @@ -445,8 +422,7 @@ private: */ vector_int m_ipiv; - //! Vector of pointers to the top of the columns of the - //! jacobians + //! Vector of pointers to the top of the columns of the jacobians /*! * The "dim" by "dim" computed Jacobian matrix for the * local Newton's method. @@ -455,8 +431,7 @@ private: //! Jacobian /*! - * m_neq by m_neq computed Jacobian matrix for the - * local Newton's method. + * m_neq by m_neq computed Jacobian matrix for the local Newton's method. */ Array2D m_Jac; @@ -472,7 +447,6 @@ private: */ vector_fp m_botBounds; - public: int m_ioflag; }; diff --git a/src/numerics/NonlinearSolver.cpp b/src/numerics/NonlinearSolver.cpp index 33885c31d..ffaea2075 100644 --- a/src/numerics/NonlinearSolver.cpp +++ b/src/numerics/NonlinearSolver.cpp @@ -1,8 +1,5 @@ /** - * - * @file NonlinearSolver.cpp - * - * Damped Newton solver for 0D and 1D problems + * @file NonlinearSolver.cpp Damped Newton solver for 0D and 1D problems */ /* @@ -34,16 +31,16 @@ using namespace std; namespace Cantera { - -//==================================================================================================================== //----------------------------------------------------------- // Constants //----------------------------------------------------------- -//! Dampfactor is the factor by which the damping factor is reduced by when a reduction in step length is warranted +//! Dampfactor is the factor by which the damping factor is reduced by when a +//! reduction in step length is warranted const doublereal DampFactor = 4.0; -//! Number of damping steps that are carried out before the solution is deemed a failure +//! Number of damping steps that are carried out before the solution is deemed +//! a failure const int NDAMP = 7; -//==================================================================================================================== + //! Print a line of a single repeated character string /*! * @param str Character string @@ -80,11 +77,6 @@ bool NonlinearSolver::s_doBothSolvesAndCompare(false); */ bool NonlinearSolver::s_alwaysAssumeNewtonGood(false); -//==================================================================================================================== -// Default constructor -/* - * @param func Residual and jacobian evaluator function object - */ NonlinearSolver::NonlinearSolver(ResidJacEval* func) : m_func(func), solnType_(NSOLN_TYPE_STEADY_STATE), @@ -206,9 +198,8 @@ NonlinearSolver::NonlinearSolver(ResidJacEval* func) : deltaX_CP_.resize(neq_, 0.0); Jd_.resize(neq_, 0.0); deltaX_trust_.resize(neq_, 1.0); - } -//==================================================================================================================== + NonlinearSolver::NonlinearSolver(const NonlinearSolver& right) : m_func(right.m_func), solnType_(NSOLN_TYPE_STEADY_STATE), @@ -298,13 +289,12 @@ NonlinearSolver::NonlinearSolver(const NonlinearSolver& right) : *this =operator=(right); } -//==================================================================================================================== NonlinearSolver::~NonlinearSolver() { delete jacCopyPtr_; delete HessianPtr_; } -//==================================================================================================================== + NonlinearSolver& NonlinearSolver::operator=(const NonlinearSolver& right) { if (this == &right) { @@ -405,18 +395,7 @@ NonlinearSolver& NonlinearSolver::operator=(const NonlinearSolver& right) return *this; } -//==================================================================================================================== -// 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 NonlinearSolver::createSolnWeights(const doublereal* const y) { for (size_t i = 0; i < neq_; i++) { @@ -428,13 +407,7 @@ void NonlinearSolver::createSolnWeights(const doublereal* const y) #endif } } -//==================================================================================================================== -// 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 NonlinearSolver::setBoundsConstraints(const doublereal* const y_low_bounds, const doublereal* const y_high_bounds) { @@ -443,38 +416,23 @@ void NonlinearSolver::setBoundsConstraints(const doublereal* const y_low_bounds, m_y_high_bounds[i] = y_high_bounds[i]; } } -//==================================================================================================================== + void NonlinearSolver::setSolverScheme(int doDogLeg, int doAffineSolve) { doDogLeg_ = doDogLeg; doAffineSolve_ = doAffineSolve; } -//==================================================================================================================== + std::vector & NonlinearSolver::lowBoundsConstraintVector() { return m_y_low_bounds; } -//==================================================================================================================== + std::vector & NonlinearSolver::highBoundsConstraintVector() { return m_y_high_bounds; } -//==================================================================================================================== -// 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. - */ + doublereal NonlinearSolver::solnErrorNorm(const doublereal* const delta_y, const char* title, int printLargest, const doublereal dampFactor) const { @@ -551,14 +509,7 @@ doublereal NonlinearSolver::solnErrorNorm(const doublereal* const delta_y, const } return sum_norm; } -//==================================================================================================================== -/* - * L2 Norm of the residual - * - * The second argument has a default of false. However, - * if true, then a table of the largest values is printed - * out to standard output. - */ + doublereal NonlinearSolver::residErrorNorm(const doublereal* const resid, const char* title, const int printLargest, const doublereal* const y) const { @@ -639,27 +590,7 @@ doublereal NonlinearSolver::residErrorNorm(const doublereal* const resid, const } return sum_norm; } -//==================================================================================================================== -// 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 NonlinearSolver::setColumnScaling(bool useColScaling, const double* const scaleFactors) { if (useColScaling) { @@ -678,23 +609,12 @@ void NonlinearSolver::setColumnScaling(bool useColScaling, const double* const s m_colScaling = 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 NonlinearSolver::setRowScaling(bool useRowScaling) { m_rowScaling = useRowScaling; } -//==================================================================================================================== -/* - * calcColumnScales(): - * - * Set the column scaling vector at the current time - */ + void NonlinearSolver::calcColumnScales() { if (m_colScaling == 1) { @@ -710,18 +630,7 @@ void NonlinearSolver::calcColumnScales() m_func->calcSolnScales(time_n, DATA_PTR(m_y_n_curr), DATA_PTR(m_y_nm1), DATA_PTR(m_colScales)); } } -//==================================================================================================================== -// Compute the current residual -/* - * @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 - * - * @return Returns a flag to indicate that operation is successful. - * 1 Means a successful operation - * -0 or neg value Means an unsuccessful operation - */ + 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 { @@ -730,14 +639,7 @@ int NonlinearSolver::doResidualCalc(const doublereal time_curr, const int typeCa m_resid_scaled = false; return retn; } -//==================================================================================================================== -// 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 - */ + void NonlinearSolver::scaleMatrix(GeneralMatrix& jac, doublereal* const y_comm, doublereal* const ydot_comm, doublereal time_curr, int num_newt_its) { @@ -883,14 +785,7 @@ void NonlinearSolver::scaleMatrix(GeneralMatrix& jac, doublereal* const y_comm, } } -//==================================================================================================================== -// Calculate the scaling factor for translating residual norms into solution norms. -/* - * This routine calls computeResidWts() a couple of times in the calculation of m_ScaleSolnNormToResNorm. - * A more sophisticated routine may do more with signs to get a better value. Perhaps, a series of calculations - * with different signs attached may be in order. Then, m_ScaleSolnNormToResNorm would be calculated - * as the minimum of a series of calculations. - */ + void NonlinearSolver::calcSolnToResNormVector() { if (! jacCopyPtr_->factored()) { @@ -949,20 +844,7 @@ void NonlinearSolver::calcSolnToResNormVector() throw CanteraError("NonlinearSolver::calcSolnToResNormVector()" , "Logic error"); } } -//==================================================================================================================== -// Compute the undamped Newton step based on the current jacobian and an input rhs -/* - * 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_curr, 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. - */ + int NonlinearSolver::doNewtonSolve(const doublereal time_curr, const doublereal* const y_curr, const doublereal* const ydot_curr, doublereal* const delta_y, GeneralMatrix& jac) @@ -1045,22 +927,16 @@ int NonlinearSolver::doNewtonSolve(const doublereal time_curr, const doublereal* m_numLocalLinearSolves++; return info; } -//==================================================================================================================== -// 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 - * - * Notes on banded Hessian solve: - * The matrix for jT j has a larger band width. Both the top and bottom band widths - * are doubled, going from KU to KU+KL and KL to KU+KL in size. This is not an impossible increase in cost, but - * has to be considered. - */ + int NonlinearSolver::doAffineNewtonSolve(const doublereal* const y_curr, const doublereal* const ydot_curr, doublereal* const delta_y, GeneralMatrix& jac) { + /* + * Notes on banded Hessian solve: The matrix for jT j has a larger band + * width. Both the top and bottom band widths are doubled, going from KU + * to KU+KL and KL to KU+KL in size. This is not an impossible increase + * in cost, but has to be considered. + */ bool newtonGood = true; // We can default to QR here ( or not ) jac.useFactorAlgorithm(1); @@ -1397,11 +1273,7 @@ int NonlinearSolver::doAffineNewtonSolve(const doublereal* const y_curr, const return info; } -//==================================================================================================================== -// Do a steepest descent calculation -/* - * This call must be made on the unfactored jacobian! - */ + doublereal NonlinearSolver::doCauchyPointSolve(GeneralMatrix& jac) { doublereal rowFac = 1.0; @@ -1531,7 +1403,7 @@ doublereal NonlinearSolver::doCauchyPointSolve(GeneralMatrix& jac) } return normSoln; } -//=================================================================================================================== + void NonlinearSolver::descentComparison(doublereal time_curr, doublereal* ydot0, doublereal* ydot1, int& numTrials) { doublereal ff = 1.0E-5; @@ -1651,21 +1523,10 @@ void NonlinearSolver::descentComparison(doublereal time_curr, doublereal* ydot0 funcDecreaseNewt2 = 0.5 * (residNewt2 - normResid02) / (ff * sNewt); printf("\t\t %10.3E %13.5E %22.16E\n", ff, funcDecreaseNewt2, residNewt); } - } - - } - } -//==================================================================================================================== -// 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 NonlinearSolver::setupDoubleDogleg() { /* @@ -1720,16 +1581,8 @@ void NonlinearSolver::setupDoubleDogleg() */ normTrust_Newton_ = calcTrustDistance(deltaX_Newton_); normTrust_CP_ = calcTrustDistance(deltaX_CP_); - } -//==================================================================================================================== -// 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 NonlinearSolver::lambdaToLeg(const doublereal lambda, doublereal& alpha) const { @@ -1743,18 +1596,9 @@ int NonlinearSolver::lambdaToLeg(const doublereal lambda, doublereal& alpha) con alpha = (lambda * dist_Total_ - dist_R0_ - dist_R1_) / dist_R2_; return 2; } -//==================================================================================================================== -// 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 NonlinearSolver::expectedResidLeg(int leg, doublereal alpha) const { - doublereal resD2, res2, resNorm; doublereal normResid02 = m_normResid_0 * m_normResid_0 * neq_; @@ -1805,16 +1649,7 @@ doublereal NonlinearSolver::expectedResidLeg(int leg, doublereal alpha) const return resNorm; } -//==================================================================================================================== -// 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 NonlinearSolver::residualComparisonLeg(const doublereal time_curr, const doublereal* const ydot0, int& legBest, doublereal& alphaBest) const { @@ -1970,18 +1805,13 @@ void NonlinearSolver::residualComparisonLeg(const doublereal time_curr, const do } } -//==================================================================================================================== -// 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 NonlinearSolver::trustRegionLength() const { norm_deltaX_trust_ = solnErrorNorm(DATA_PTR(deltaX_trust_)); return trustDelta_ * norm_deltaX_trust_; } -//==================================================================================================================== + void NonlinearSolver::setDefaultDeltaBoundsMagnitudes() { for (size_t i = 0; i < neq_; i++) { @@ -1989,7 +1819,7 @@ void NonlinearSolver::setDefaultDeltaBoundsMagnitudes() m_deltaStepMinimum[i] = std::max(m_deltaStepMinimum[i], 0.1 * fabs(m_y_n_curr[i])); } } -//==================================================================================================================== + void NonlinearSolver::adjustUpStepMinimums() { for (size_t i = 0; i < neq_; i++) { @@ -2000,37 +1830,18 @@ void NonlinearSolver::adjustUpStepMinimums() } } -//==================================================================================================================== + void NonlinearSolver::setDeltaBoundsMagnitudes(const doublereal* const deltaStepMinimum) { - for (size_t i = 0; i < neq_; i++) { m_deltaStepMinimum[i] = deltaStepMinimum[i]; } m_manualDeltaStepSet = 1; } -//==================================================================================================================== -/* - * - * 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_n_curr Initial value of the solution vector - * @param step_1 initial proposed step size - * - * @return returns the damping factor - */ + double NonlinearSolver::deltaBoundStep(const doublereal* const y_n_curr, const doublereal* const step_1) { - size_t i_fbounds = 0; int ifbd = 0; int i_fbd = 0; @@ -2125,13 +1936,7 @@ NonlinearSolver::deltaBoundStep(const doublereal* const y_n_curr, const doublere return f_delta_bounds; } -//==================================================================================================================== -// 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. - */ + void NonlinearSolver::readjustTrustVector() { doublereal trustDeltaOld = trustDelta_; @@ -2194,11 +1999,7 @@ void NonlinearSolver::readjustTrustVector() trustNormGoal, deltaXSizeOld, trustDeltaOld, norm_deltaX_trust_, trustDelta_); } } -//==================================================================================================================== -//! 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 NonlinearSolver::initializeTrustRegion() { if (trustRegionInitializationMethod_ == 0) { @@ -2242,16 +2043,6 @@ void NonlinearSolver::initializeTrustRegion() } } -//==================================================================================================================== -// 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 NonlinearSolver::fillDogLegStep(int leg, doublereal alpha, std::vector & deltaX) const { if (leg == 0) { @@ -2268,16 +2059,7 @@ void NonlinearSolver::fillDogLegStep(int leg, doublereal alpha, std::vector const& deltaX) const { doublereal sum = 0.0; @@ -2288,15 +2070,7 @@ doublereal NonlinearSolver::calcTrustDistance(std::vector const& de } return sqrt(sum / neq_) / trustDelta_; } -//==================================================================================================================== -// 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(doublereal trustDelta, doublereal& lambda, doublereal& alpha) const { doublereal dist; @@ -2335,34 +2109,7 @@ int NonlinearSolver::calcTrustIntersection(doublereal trustDelta, doublereal& la lambda = dist / dist_Total_; return 1; } -//==================================================================================================================== -/* - * - * boundStep(): - * - * 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 - */ + doublereal NonlinearSolver::boundStep(const doublereal* const y, const doublereal* const step0) { size_t i_lower = npos; @@ -2412,30 +2159,7 @@ doublereal NonlinearSolver::boundStep(const doublereal* const y, const doublerea doublereal f_delta_bounds = deltaBoundStep(y, step0); return std::min(f_bounds, f_delta_bounds); } -//=================================================================================================================== -// Find a damping coefficient through a look-ahead mechanism -/* - * - * On entry, step0 must contain an undamped Newton step to the - * current solution y0. This method attempts to find a damping coefficient - * such that the next undamped step would have a norm smaller than - * that of step0. If successful, the new solution after taking the - * damped step is returned in y1, and the undamped step at y1 is - * returned in step1. - * - * - * @return 1 Successful step was taken: Next step was less than previous step. - * s1 is calculated - * 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 - * NSOLN_RETN_FAIL_DAMPSTEP - * Unsuccessful step. We can not find a damping factor that is suitable. - */ + int NonlinearSolver::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* const step_2, @@ -2643,31 +2367,7 @@ int NonlinearSolver::dampStep(const doublereal time_curr, const doublereal* cons } return NSOLN_RETN_FAIL_DAMPSTEP; } -//==================================================================================================================== -// 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 s1 OUTPUT Norm of the vector step_1 - * @param jac INPUT jacobian - * @param numTrials OUTPUT number of trials taken in the current damping step - * - * - * @return 1 Success: Good step was taken. The predicted residual norm is less than one - * 2 Success: Good 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 Unknown Uncertain Success: s1 is about the same as s0 - * NSOLN_RETN_FAIL_DAMPSTEP - * Unsuccessful step. Can not find a damping coefficient that is suitable - */ + int NonlinearSolver::dampDogLeg(const doublereal time_curr, const doublereal* y_n_curr, const doublereal* ydot_n_curr, std::vector & step_1, doublereal* const y_n_1, doublereal* const ydot_n_1, @@ -2816,34 +2516,7 @@ int NonlinearSolver::dampDogLeg(const doublereal time_curr, const doublereal* y_ } return NSOLN_RETN_FAIL_DAMPSTEP; } -//==================================================================================================================== -// 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 y0 INPUT Current value of the solution vector - * @param ydot0 INPUT Current value of the derivative of the solution vector - * @param step0 INPUT Trial step - * @param y1 OUTPUT Solution values at the conditions which are evaluated for success - * @param ydot1 OUTPUT Time derivatives of solution at the conditions which are evalulated 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 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 & step_1, @@ -3010,21 +2683,7 @@ int NonlinearSolver::decideStep(const doublereal time_curr, int leg, doublereal } return retn; } -//==================================================================================================================== -/* - * solve_nonlinear_problem(): - * - * 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. - * - * @return A positive value indicates a successful convergence - * -1 Failed convergence - */ + int NonlinearSolver::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, @@ -3585,36 +3244,14 @@ done: return retnCode; } -//==================================================================================================================== -//! 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 - */ + void NonlinearSolver:: setPreviousTimeStep(const std::vector& y_nm1, const std::vector& ydot_nm1) { m_y_nm1 = y_nm1; m_ydot_nm1 = ydot_nm1; } -//==================================================================================================================== -// 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 NonlinearSolver:: print_solnDelta_norm_contrib(const doublereal* const step_1, const char* const stepNorm_1, @@ -3702,21 +3339,7 @@ static inline doublereal subtractRD(doublereal a, doublereal b) } return diff; } -//==================================================================================================================== -/* - * - * Function called by BEuler to evaluate the Jacobian matrix and the - * current residual at the current time step. - * @param N = The size of the equation system - * @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. - * - * @return Returns a flag to indicate that operation is successful. - * 1 Means a successful operation - * 0 Means an unsuccessful operation - */ + int NonlinearSolver::beuler_jac(GeneralMatrix& J, doublereal* const f, doublereal time_curr, doublereal CJ, doublereal* const y, doublereal* const ydot, @@ -3989,16 +3612,7 @@ int NonlinearSolver::beuler_jac(GeneralMatrix& J, doublereal* const f, return retn; } -//==================================================================================================================== -// 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 NonlinearSolver:: calc_ydot(const int order, const doublereal* const y_curr, doublereal* const ydot_curr) const { @@ -4025,49 +3639,20 @@ calc_ydot(const int order, const doublereal* const y_curr, doublereal* const ydo throw CanteraError("calc_ydot()", "Case not covered"); } } -//==================================================================================================================== -// Apply a filtering process to the new step -/* - * @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 NonlinearSolver::filterNewStep(const doublereal timeCurrent, const doublereal* const ybase, doublereal* const step0) { return m_func->filterNewStep(timeCurrent, ybase, step0); } -//==================================================================================================================== -// Apply a filtering process to the new 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 NonlinearSolver::filterNewSolution(const doublereal timeCurrent, doublereal* const y_current, doublereal* const ydot_current) { return m_func->filterSolnPrediction(timeCurrent, y_current); } -//==================================================================================================================== -// 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 rtol and an atol factor - * is added as well so that if the residual is less than 1, then the calculation is deemed to be converged. - * - * 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 should also be achieved - * and should be checked. - */ -void -NonlinearSolver::computeResidWts() + +void NonlinearSolver::computeResidWts() { ResidWtsReevaluated_ = true; if (checkUserResidualTols_ == 1) { @@ -4098,34 +3683,15 @@ NonlinearSolver::computeResidWts() } } } -//===================================================================================================================== -// return the residual weights -/* - * @param residWts Vector of length neq_ - */ -void -NonlinearSolver::getResidWts(doublereal* const residWts) const + +void NonlinearSolver::getResidWts(doublereal* const residWts) const { for (size_t i = 0; i < neq_; i++) { residWts[i] = (m_residWts)[i]; } } -//===================================================================================================================== -// Check to see if the nonlinear problem has converged -/* - * - * @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.0E-2 - * The predicted deltaSoln is below 1.0E-2. - * 0 Not converged yet - */ -int -NonlinearSolver::convergenceCheck(int dampCode, doublereal s1) + +int NonlinearSolver::convergenceCheck(int dampCode, doublereal s1) { int retn = 0; if (m_dampBound < 0.9999) { @@ -4175,13 +3741,7 @@ NonlinearSolver::convergenceCheck(int dampCode, doublereal s1) } return retn; } -//===================================================================================================================== -// 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 NonlinearSolver::setAtol(const doublereal* const atol) { for (size_t i = 0; i < neq_; i++) { @@ -4192,13 +3752,7 @@ void NonlinearSolver::setAtol(const doublereal* const atol) atolk_[i]= atol[i]; } } -//===================================================================================================================== -// Set the relative tolerances for the solution variables -/* - * Set the relative tolerances used in the calculation for the solution variables. - * - * @param rtol single double - */ + void NonlinearSolver::setRtol(const doublereal rtol) { if (rtol <= 0.0) { @@ -4207,20 +3761,7 @@ void NonlinearSolver::setRtol(const doublereal rtol) } 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 calculated by method 0 and 1. - * 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) { @@ -4241,11 +3782,10 @@ void NonlinearSolver::setResidualTols(double residRtol, double* residATol, int } } } -//===================================================================================================================== + void NonlinearSolver::setPrintLvl(int printLvl) { m_print_flag = printLvl; } -//===================================================================================================================== -} +} diff --git a/src/numerics/ResidJacEval.cpp b/src/numerics/ResidJacEval.cpp index 8d8944002..9c4cc3de2 100644 --- a/src/numerics/ResidJacEval.cpp +++ b/src/numerics/ResidJacEval.cpp @@ -20,23 +20,18 @@ using namespace std; namespace Cantera { -//==================================================================================================================== - ResidJacEval::ResidJacEval(doublereal atol) : ResidEval(), m_atol(atol) { } -//==================================================================================================================== -// Copy Constructor for the %ResidJacEval object -/* - */ + ResidJacEval::ResidJacEval(const ResidJacEval& right) : ResidEval() { *this = operator=(right); } -//==================================================================================================================== + ResidJacEval& ResidJacEval::operator=(const ResidJacEval& right) { if (this == &right) { @@ -50,30 +45,17 @@ ResidJacEval& ResidJacEval::operator=(const ResidJacEval& right) return *this; } -//==================================================================================================================== -// Duplication routine for objects which inherit from %ResidJacEval -/* - * This virtual routine can be used to duplicate %ResidJacEval objects - * inherited from %ResidJacEval even if the application only has - * a pointer to %ResidJacEval to work with. - * - * These routines are basically wrappers around the derived copy - * constructor. - */ + ResidJacEval* ResidJacEval::duplMyselfAsResidJacEval() const { return new ResidJacEval(*this); } -//==================================================================================================================== + int ResidJacEval::nEquations() const { return neq_; } -//==================================================================================================================== -// Set a global value of the absolute tolerance -/* - * @param atol Value of atol - */ + void ResidJacEval::setAtol(doublereal atol) { m_atol = atol; @@ -82,15 +64,7 @@ void ResidJacEval::setAtol(doublereal atol) "atol must be greater than zero"); } } -//==================================================================================================================== -//! Fill in the initial conditions -/*! - * Values for both the solution and the value of ydot may be provided. - * - * @param t0 Time (input) - * @param y Solution vector (output) - * @param ydot Rate of change of solution vector. (output) - */ + int ResidJacEval:: getInitialConditions(doublereal t0, doublereal* const y, doublereal* const ydot) { @@ -104,79 +78,27 @@ getInitialConditions(doublereal t0, doublereal* const y, doublereal* const ydot) } return 1; } -//==================================================================================================================== -// This function may be used to create output at various points in the execution of an application. -/* - * - * @param ifunc identity of the call - * 0 Initial call - * 1 Called at the end of every successful time step - * -1 Called at the end of every unsuccessful time step - * 2 Called at the end of every call to integrateRJE() - * - * @param t Time (input) - * @param delta_t The current value of the time step (input) - * @param y Solution vector (input, do not modify) - * @param ydot Rate of change of solution vector. (input) - */ + void ResidJacEval:: user_out2(const int ifunc, const doublereal t, const doublereal deltaT, const doublereal* y, const doublereal* ydot) { } -//==================================================================================================================== -// This function may be used to create output at various points in the execution of an application. -/* - * This routine calls user_out2(). - * - * @param ifunc identity of the call - * @param t Time (input) - * @param y Solution vector (input, do not modify) - * @param ydot Rate of change of solution vector. (input) - */ void ResidJacEval:: user_out(const int ifunc, const doublereal t, const doublereal* y, const doublereal* ydot) { user_out2(ifunc, t, 0.0, y, ydot); } -//==================================================================================================================== -//! Evaluate the time tracking equations, if any -/*! - * Evaluate time integrated quantities that are calculated at the - * end of every successful time step. This call is made once at the end of every successful - * time step that advances the time. It's also made once at the start of the time stepping. - * - * @param t Time (input) - * @param delta_t The current value of the time step (input) - * @param y Solution vector (input, do not modify) - * @param ydot Rate of change of solution vector. (input, do not modify) - */ + int ResidJacEval:: evalTimeTrackingEqns(const doublereal t, const doublereal delta_t, const doublereal* y, const doublereal* ydot) { return 1; } -//==================================================================================================================== -// Return a vector of delta y's for calculation of the numerical Jacobian -/* - * There is a default algorithm provided. - * - * delta_y[i] = atol[i] + 1.0E-6 ysoln[i] - * delta_y[i] = atol[i] + MAX(1.0E-6 ysoln[i] * 0.01 * solnWeights[i]) - * - * @param t Time (input) - * @param y Solution vector (input, do not modify) - * @param ydot Rate of change of solution vector. (input, do not modify) - * @param delta_y Value of the delta to be used in calculating the numerical jacobian - * @param solnWeights Value of the solution weights that are used in determining convergence (default = 0) - * - * @return Returns a flag to indicate that operation is successful. - * 1 Means a successful operation - * 0 Means an unsuccessful operation - */ + int ResidJacEval:: calcDeltaSolnVariables(const doublereal t, const doublereal* const ySoln, const doublereal* const ySolnDot, doublereal* const deltaYSoln, @@ -193,16 +115,7 @@ calcDeltaSolnVariables(const doublereal t, const doublereal* const ySoln, } return 1; } -//==================================================================================================================== -// Returns a vector of column scale factors that can be used to column scale Jacobians. -/* - * Default to yScales[] = 1.0 - * - * @param t Time (input) - * @param y Solution vector (input, do not modify) - * @param y_old Old Solution vector (input, do not modify) - * @param yScales Value of the column scales - */ + void ResidJacEval:: calcSolnScales(const doublereal t, const doublereal* const ysoln, const doublereal* const ysolnOld, doublereal* const ysolnScales) @@ -215,47 +128,17 @@ calcSolnScales(const doublereal t, const doublereal* const ysoln, const doublere } } } -//==================================================================================================================== -// Filter the solution predictions -/* - * Codes might provide a predicted step change. This routine filters the predicted - * solution vector eliminating illegal directions. - * - * @param t Time (input) - * @param y Solution vector (input, output) - * @param step Proposed step in the solution that will be cropped - */ + doublereal ResidJacEval::filterNewStep(doublereal t, const doublereal* const ybase, doublereal* const step) { return 0.0; } -//==================================================================================================================== -// Filter the solution predictions -/* - * Codes might provide a predicted solution vector. This routine filters the predicted - * solution vector. - * - * @param t Time (input) - * @param y Solution vector (input, output) - */ + doublereal ResidJacEval::filterSolnPrediction(doublereal t, doublereal* const y) { return 0.0; } -//==================================================================================================================== -// Evaluate any stopping criteria other than a final time limit -/* - * If we are to stop the time integration for any reason other than reaching a final time limit, tout, - * provide a test here. This call is made at the end of every successful time step iteration - * - * @return If true, the the time stepping is stopped. If false, then time stepping is stopped if t >= tout - * Defaults to false. - * - * @param t Time (input) - * @param delta_t The current value of the time step (input) - * @param y Solution vector (input, do not modify) - * @param ydot Rate of change of solution vector. (input, do not modify) - */ + bool ResidJacEval:: evalStoppingCritera(const doublereal t, const doublereal delta_t, @@ -264,42 +147,13 @@ evalStoppingCritera(const doublereal t, { return false; } -//==================================================================================================================== -// Multiply the matrix by another matrix that leads to better conditioning -/* - * Provide a left sided matrix that will multiply the current jacobian, after scaling - * and lead to a better conditioned system. - * This routine is called just before the matrix is factored. - * - * Original Problem: - * J delta_x = - Resid - * - * New problem: - * M (J delta_x) = - M Resid - * - * @param matrix Pointer to the current jacobian (if zero, it's already been factored) - * @param nrows offsets for the matrix - * @param rhs residual vector. This also needs to be lhs multiplied by M - */ + int ResidJacEval:: matrixConditioning(doublereal* const matrix, const int nrows, doublereal* const rhs) { return 1; } -//==================================================================================================================== -// Evaluate the residual function -/* - * @param t Time (input) - * @param delta_t The current value of the time step (input) - * @param y Solution vector (input, do not modify) - * @param ydot Rate of change of solution vector. (input, do not modify) - * @param resid Value of the residual that is computed (output) - * @param evalType Type of the residual being computed (defaults to Base_ResidEval) - * @param id_x Index of the variable that is being numerically differenced to find - * the jacobian (defaults to -1, which indicates that no variable is being - * differenced or that the residual doesn't take this issue into account) - * @param delta_x Value of the delta used in the numerical differencing - */ + int ResidJacEval:: evalResidNJ(const doublereal t, const doublereal deltaT, const doublereal* y, const doublereal* ydot, doublereal* const resid, const ResidEval_Type_Enum evalType, @@ -308,25 +162,14 @@ evalResidNJ(const doublereal t, const doublereal deltaT, const doublereal* y, throw CanteraError("ResidJacEval::evalResidNJ()", "Not implemented\n"); return 1; } -//==================================================================================================================== + int ResidJacEval::eval(const doublereal t, const doublereal* const y, const doublereal* const ydot, doublereal* const r) { double deltaT = -1.0; return evalResidNJ(t, deltaT, y, ydot, r); } -//==================================================================================================================== -// Calculate an analytical jacobian and the residual at the current time and values. -/* - * Only called if the jacFormation method is set to analytical - * - * @param t Time (input) - * @param delta_t The current value of the time step (input) - * @param y Solution vector (input, do not modify) - * @param ydot Rate of change of solution vector. (input, do not modify) - * @param J Reference to the SquareMatrix object to be calculated (output) - * @param resid Value of the residual that is computed (output) - */ + int ResidJacEval:: evalJacobian(const doublereal t, const doublereal delta_t, doublereal cj, const doublereal* const y, @@ -337,19 +180,7 @@ evalJacobian(const doublereal t, const doublereal delta_t, doublereal cj, doublereal* const* jac_colPts = J.colPts(); return evalJacobianDP(t, delta_t, cj, y, ydot, jac_colPts, resid); } -//==================================================================================================================== -// Calculate an analytical jacobian and the residual at the current time and values. -/* - * Only called if the jacFormation method is set to analytical - * - * @param t Time (input) - * @param delta_t The current value of the time step (input) - * @param c_j The current value of the coefficient of the time derivative - * @param y Solution vector (input, do not modify) - * @param ydot Rate of change of solution vector. (input, do not modify) - * @param jac_colPts Reference to the SquareMatrix object to be calculated (output) - * @param resid Value of the residual that is computed (output) - */ + int ResidJacEval:: evalJacobianDP(const doublereal t, const doublereal delta_t, const doublereal c_j, @@ -361,6 +192,5 @@ evalJacobianDP(const doublereal t, const doublereal delta_t, throw CanteraError("ResidJacEval::evalJacobianDP()", "Not implemented\n"); return 1; } -//==================================================================================================================== } diff --git a/src/numerics/RootFind.cpp b/src/numerics/RootFind.cpp index 981b87dc6..bef5b8180 100644 --- a/src/numerics/RootFind.cpp +++ b/src/numerics/RootFind.cpp @@ -40,9 +40,6 @@ namespace Cantera #define DSIGN(x) (( (x) == (0.0) ) ? (0.0) : ( ((x) > 0.0) ? 1.0 : -1.0 )) #endif -/*****************************************************************************/ -/*****************************************************************************/ -/*****************************************************************************/ #ifdef DEBUG_MODE //! Print out a form for the current function evaluation /*! @@ -64,8 +61,6 @@ static void print_funcEval(FILE* fp, doublereal xval, doublereal fval, int its) } #endif -//================================================================================================ -// Main constructor RootFind::RootFind(ResidEval* resid) : m_residFunc(resid), m_funcTargetValue(0.0), @@ -88,9 +83,8 @@ RootFind::RootFind(ResidEval* resid) : x_minTried_(1.0E300), fx_minTried_(0.0) { - } -//================================================================================================ + RootFind::RootFind(const RootFind& r) : m_residFunc(r.m_residFunc), m_funcTargetValue(0.0), @@ -115,12 +109,11 @@ RootFind::RootFind(const RootFind& r) : { *this = r; } -//================================================================================================ -// Empty destructor + RootFind::~RootFind() { } -//==================================================================================================================== + RootFind& RootFind::operator=(const RootFind& right) { if (this == &right) { @@ -149,14 +142,7 @@ RootFind& RootFind::operator=(const RootFind& right) return *this; } -//================================================================================================ -// Calculate a deltaX from an input value of x -/* - * This routine ensure that the deltaX will be greater or equal to DeltaXNorm_ - * or 1.0E-14 x - * - * @param x1 input value of x - */ + doublereal RootFind::delXNonzero(doublereal x1) const { doublereal deltaX = 1.0E-14 * fabs(x1); @@ -166,14 +152,7 @@ doublereal RootFind::delXNonzero(doublereal x1) const } return deltaX; } -//================================================================================================ -// Calculate a deltaX from an input value of x -/* - * This routine ensure that the deltaX will be greater or equal to DeltaXNorm_ - * or 1.0E-14 x or deltaXConverged_. - * - * @param x1 input value of x - */ + doublereal RootFind::delXMeaningful(doublereal x1) const { doublereal del = delXNonzero(x1); @@ -182,16 +161,7 @@ doublereal RootFind::delXMeaningful(doublereal x1) const } return del; } -//================================================================================================ -// Calculate a controlled, nonzero delta between two numbers -/* - * The delta is designed to be greater than or equal to delXNonzero(x) defined above - * with the same sign as the original delta. Therefore if you subtract it from either - * of the two original numbers, you get a different number. - * - * @param x1 first number - * @param x2 second number - */ + double RootFind::deltaXControlled(doublereal x2, doublereal x1) const { doublereal sgnn = 1.0; @@ -206,18 +176,7 @@ double RootFind::deltaXControlled(doublereal x2, doublereal x1) const } return deltaX; } -//==================================================================================================================== -// Function to decide whether two real numbers are the same or not -/* - * A comparison is made between the two numbers to decide whether they - * are close to one another. This is defined as being within factor * delXMeaningful() of each other. - * - * @param x2 First number - * @param x1 second number - * @param factor Multiplicative factor for delta X. defaults to 1 - * - * @return Returns a boolean indicating whether the two numbers are the same or not. - */ + bool RootFind::theSame(doublereal x2, doublereal x1, doublereal factor) const { doublereal x = fabs(x2) + fabs(x1); @@ -229,26 +188,13 @@ bool RootFind::theSame(doublereal x2, doublereal x1, doublereal factor) const } return false; } -//==================================================================================================================== -/* - * The following calculation is a line search method to find the root of a function - * - * - * xbest Returns the x that satisfies the function - * On input, xbest should contain the best estimate - * - * return: - * 0 Found function - */ + int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& funcTargetValue, doublereal* xbest) { - /* * We store the function target and then actually calculate a modified functional * * func = eval(x1) - m_funcTargetValue = 0 - * - * */ m_funcTargetValue = funcTargetValue; @@ -1176,7 +1122,7 @@ done: return retn; } -//==================================================================================================================== + doublereal RootFind::func(doublereal x) { doublereal r; @@ -1198,15 +1144,7 @@ doublereal RootFind::func(doublereal x) } return ff; } -//==================================================================================================================== -// Set the tolerance parameters for the rootfinder -/* - * These tolerance parameters are used on the function value to determine convergence - * - * - * @param rtol Relative tolerance. The default is 10^-5 - * @param atol absolute tolerance. The default is 10^-11 - */ + void RootFind::setTol(doublereal rtolf, doublereal atolf, doublereal rtolx, doublereal atolx) { m_atolf = atolf; @@ -1222,38 +1160,12 @@ void RootFind::setTol(doublereal rtolf, doublereal atolf, doublereal rtolx, doub m_atolx = atolx; } } -//==================================================================================================================== -// Set the print level from the rootfinder -/* - * - * 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 RootFind::setPrintLvl(int printlvl) { printLvl = printlvl; } -//==================================================================================================================== -// Set the function behavior flag -/* - * If this is true, the function is generally an increasing function of x. - * In particular, if the algorithm is seeking a higher value of f, it will look - * in the positive x direction. - * - * This type of function is needed because this algorithm must deal with regions of f(x) where - * f is not changing with x. - * - * @param value boolean value - */ + void RootFind::setFuncIsGenerallyIncreasing(bool value) { if (value) { @@ -1261,18 +1173,7 @@ void RootFind::setFuncIsGenerallyIncreasing(bool value) } FuncIsGenerallyIncreasing_ = value; } -//==================================================================================================================== -// Set the function behavior flag -/* - * If this is true, the function is generally a decreasing function of x. - * In particular, if the algorithm is seeking a higher value of f, it will look - * in the negative x direction. - * - * This type of function is needed because this algorithm must deal with regions of f(x) where - * f is not changing with x. - * - * @param value boolean value - */ + void RootFind::setFuncIsGenerallyDecreasing(bool value) { if (value) { @@ -1280,33 +1181,19 @@ void RootFind::setFuncIsGenerallyDecreasing(bool value) } FuncIsGenerallyDecreasing_ = value; } -//==================================================================================================================== -// Set the nominal value of deltaX -/* - * This sets the value of deltaXNorm_ - * - * @param deltaXNorm - */ + void RootFind::setDeltaX(doublereal deltaXNorm) { DeltaXnorm_ = deltaXNorm; specifiedDeltaXnorm_ = 1; } -//==================================================================================================================== -// Set the maximum value of deltaX -/* - * This sets the value of deltaXMax_ - * - * @param deltaX - */ + void RootFind::setDeltaXMax(doublereal deltaX) { DeltaXMax_ = deltaX; specifiedDeltaXMax_ = 1; } -//==================================================================================================================== -//==================================================================================================================== void RootFind::printTable() { printf("\t----------------------------------------------------------------------------------------------------------------------------------------\n"); @@ -1322,6 +1209,5 @@ void RootFind::printTable() } printf("\t----------------------------------------------------------------------------------------------------------------------------------------\n"); } -//==================================================================================================================== } diff --git a/src/numerics/funcs.cpp b/src/numerics/funcs.cpp index adbd3acbb..72d7534b8 100644 --- a/src/numerics/funcs.cpp +++ b/src/numerics/funcs.cpp @@ -7,6 +7,8 @@ * See file License.txt for licensing information */ +#include "cantera/numerics/funcs.h" + #include #include @@ -16,7 +18,6 @@ using namespace std; #include "cantera/base/ct_defs.h" #include "cantera/base/ctexceptions.h" #include "cantera/base/stringUtils.h" -#include "cantera/numerics/funcs.h" #include "cantera/numerics/polyfit.h" #ifndef FTN_TRAILING_UNDERSCORE @@ -38,21 +39,6 @@ extern "C" { namespace Cantera { -// Linearly interpolate a function defined on a discrete grid. -/* - * Vector xpts contains a monotonic sequence of grid points, and - * vector fpts contains function values defined at these points. - * The value returned is the linear interpolate at point x. - * If x is outside the range of xpts, the value of fpts at the - * nearest end is returned. - * - * @param x value of the x coordinate - * @param xpts value of the grid points - * @param fpts value of the interpolant at the grid points - * - * @return Returned value is the value of of the interpolated - * function at x. - */ doublereal linearInterp(doublereal x, const vector_fp& xpts, const vector_fp& fpts) { diff --git a/src/numerics/solveProb.cpp b/src/numerics/solveProb.cpp index 3814f17dd..202383753 100644 --- a/src/numerics/solveProb.cpp +++ b/src/numerics/solveProb.cpp @@ -30,8 +30,6 @@ namespace Cantera static doublereal calcWeightedNorm(const doublereal [], const doublereal dx[], size_t); -//================================================================================================ -// Main constructor solveProb::solveProb(ResidEval* resid) : m_residFunc(resid), m_neq(0), @@ -67,22 +65,20 @@ solveProb::solveProb(ResidEval* resid) : } } -//================================================================================================ -// Empty destructor + solveProb::~solveProb() { } -//================================================================================================ -/* - * The following calculation is a Newton's method to - * get the surface fractions of the surface and bulk species by - * requiring that the - * surface species production rate = 0 and that the bulk fractions are - * proportional to their production rates. - */ + int solveProb::solve(int ifunc, doublereal time_scale, doublereal reltol) { + /* + * The following calculation is a Newton's method to get the surface fractions + * of the surface and bulk species by requiring that the surface species + * production rate = 0 and that the bulk fractions are proportional to their + * production rates. + */ doublereal EXTRA_ACCURACY = 0.001; if (ifunc == SOLVEPROB_JACOBIAN) { EXTRA_ACCURACY *= 0.001; @@ -392,43 +388,27 @@ int solveProb::solve(int ifunc, doublereal time_scale, } return 0; } -//================================================================================================ -/* - * Update the surface states of the surface phases. - */ + void solveProb::reportState(doublereal* const CSolnSP) const { std::copy(m_CSolnSP.begin(), m_CSolnSP.end(), CSolnSP); } -//================================================================================================ -/* - * This calculates the net production rates of all species - * - * This calculates the function eval. - * (should switch to special_species formulation for sum condition) - * - * @internal - * This routine uses the m_numEqn1 and m_netProductionRatesSave vectors - * as temporary internal storage. - */ + void solveProb::fun_eval(doublereal* const resid, const doublereal* const CSoln, const doublereal* const CSolnOld, const bool do_time, const doublereal deltaT) { + /* + * This routine uses the m_numEqn1 and m_netProductionRatesSave vectors + * as temporary internal storage. + */ if (do_time) { m_residFunc->evalSimpleTD(0.0, CSoln, CSolnOld, deltaT, resid); } else { m_residFunc->evalSS(0.0, CSoln, resid); } } -//================================================================================================ -/* - * Calculate the Jacobian and residual - * - * @internal - * This routine uses the m_numEqn2 vector - * as temporary internal storage. - */ + void solveProb::resjac_eval(std::vector &JacCol, doublereal resid[], doublereal CSoln[], const doublereal CSolnOld[], const bool do_time, @@ -452,6 +432,7 @@ void solveProb::resjac_eval(std::vector &JacCol, } dc = std::max(1.0E-11 * sd, fabs(cSave) * 1.0E-6); CSoln[kCol] += dc; + // Use the m_numEqn2 vector as temporary internal storage. fun_eval(DATA_PTR(m_numEqn2), CSoln, CSolnOld, do_time, deltaT); col_j = JacCol[kCol]; for (size_t i = 0; i < m_neq; i++) { @@ -461,26 +442,9 @@ void solveProb::resjac_eval(std::vector &JacCol, } } -//================================================================================================ + #define APPROACH 0.50 -// This function calculates a damping factor for the Newton iteration update -// vector, dxneg, to insure that all solution components stay within prescribed bounds -/* - * The default for this class is that all solution components are bounded between zero and one. - * this is because the original unknowns were mole fractions and surface site fractions. - * - * dxneg[] = negative of the update vector. - * - * The constant "APPROACH" sets the fraction of the distance to the boundary - * that the step can take. If the full step would not force any fraction - * outside of the bounds, then Newton's method is mostly allowed to operate normally. - * There is also some solution damping employed. - * - * @param x Vector of the current solution components - * @param dxneg Vector of the negative of the full solution update vector. - * @param dim Size of the solution vector - * @param label return int, stating which solution component caused the most damping. - */ + doublereal solveProb::calc_damping(doublereal x[], doublereal dxneg[], size_t dim, size_t* label) { doublereal damp = 1.0, xnew, xtop, xbot; @@ -554,7 +518,7 @@ doublereal solveProb::calc_damping(doublereal x[], doublereal dxneg[], size_t di } #undef APPROACH -//================================================================================================ + /* * This function calculates the norm of an update, dx[], * based on the weighted values of x. @@ -572,12 +536,7 @@ static doublereal calcWeightedNorm(const doublereal wtX[], const doublereal dx[] } return sqrt(norm/dim); } -//================================================================================================ -/* - * Calculate the weighting factors for norms wrt both the species - * concentration unknowns and the residual unknowns. - * - */ + void solveProb::calcWeights(doublereal wtSpecies[], doublereal wtResid[], const doublereal CSoln[]) { @@ -601,18 +560,7 @@ void solveProb::calcWeights(doublereal wtSpecies[], doublereal wtResid[], } } } -//================================================================================================ -/* - * This routine calculates a pretty conservative 1/del_t based - * on MAX_i(sdot_i/(X_i*SDen0)). This probably guarantees - * diagonal dominance. - * - * Small surface fractions are allowed to intervene in the del_t - * determination, no matter how small. This may be changed. - * Now minimum changed to 1.0e-12, - * - * Maximum time step set to time_scale. - */ + doublereal solveProb:: calc_t(doublereal netProdRateSolnSP[], doublereal Csoln[], size_t* label, size_t* label_old, doublereal* label_factor, int ioflag) @@ -665,14 +613,7 @@ calc_t(doublereal netProdRateSolnSP[], doublereal Csoln[], return inv_timeScale; } -//==================================================================================================================== -// Set the bottom and top bounds on the solution vector -/* - * The default is for the bottom is 0.0, while the default for the top is 1.0 - * - * @param botBounds Vector of bottom bounds - * @param topBounds vector of top bounds - */ + void solveProb::setBounds(const doublereal botBounds[], const doublereal topBounds[]) { for (size_t k = 0; k < m_neq; k++) { @@ -680,11 +621,7 @@ void solveProb::setBounds(const doublereal botBounds[], const doublereal topBoun m_topBounds[k] = topBounds[k]; } } -//==================================================================================================================== -/* - * printResJac(): prints out the residual and Jacobian. - * - */ + #ifdef DEBUG_SOLVEPROB void solveProb::printResJac(int ioflag, int neq, const Array2D& Jac, doublereal resid[], doublereal wtRes[], @@ -693,10 +630,7 @@ void solveProb::printResJac(int ioflag, int neq, const Array2D& Jac, } #endif -//================================================================================================ -/* - * Optional printing at the start of the solveProb problem - */ + void solveProb::print_header(int ioflag, int ifunc, doublereal time_scale, doublereal reltol, doublereal netProdRate[]) @@ -881,9 +815,8 @@ void solveProb::printIteration(int ioflag, doublereal damp, size_t label_d, "------------------------------\n"); } #endif -} /* printIteration */ +} -//================================================================================================ void solveProb::printFinal(int ioflag, doublereal damp, size_t label_d, size_t label_t, doublereal inv_t, doublereal t_real, int iter, doublereal update_norm, doublereal resid_norm, @@ -971,7 +904,7 @@ void solveProb::printFinal(int ioflag, doublereal damp, size_t label_d, size_t l } #endif } -//================================================================================================ + #ifdef DEBUG_SOLVEPROB void solveProb:: printIterationHeader(int ioflag, doublereal damp,doublereal inv_t, doublereal t_real, @@ -995,21 +928,19 @@ printIterationHeader(int ioflag, doublereal damp,doublereal inv_t, doublereal t_ } } #endif -//================================================================================================ + void solveProb::setAtol(const doublereal atol[]) { for (size_t k = 0; k < m_neq; k++, k++) { m_atol[k] = atol[k]; } } -//================================================================================================ + void solveProb::setAtolConst(const doublereal atolconst) { for (size_t k = 0; k < m_neq; k++, k++) { m_atol[k] = atolconst; } } -//================================================================================================ - }