From 74f33b882938880c7e6858df955ed5e5a5fe88d2 Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Tue, 27 Jul 2010 01:10:58 +0000 Subject: [PATCH] Added a topBounds and botBounds capability --- Cantera/src/numerics/solveProb.cpp | 62 ++++++++++++++++++++---------- Cantera/src/numerics/solveProb.h | 53 +++++++++++++++++++++++-- 2 files changed, 91 insertions(+), 24 deletions(-) diff --git a/Cantera/src/numerics/solveProb.cpp b/Cantera/src/numerics/solveProb.cpp index 5bf2b766b..2da41b13c 100644 --- a/Cantera/src/numerics/solveProb.cpp +++ b/Cantera/src/numerics/solveProb.cpp @@ -79,6 +79,8 @@ namespace Cantera { m_wtSpecies.resize(dim1, 0.0); m_resid.resize(dim1, 0.0); m_ipiv.resize(dim1, 0); + m_topBounds.resize(dim1, 1.0); + m_botBounds.resize(dim1, 0.0); m_Jac.resize(dim1, dim1, 0.0); m_JacCol.resize(dim1, 0); @@ -473,47 +475,55 @@ namespace Cantera { } //================================================================================================ #define APPROACH 0.50 - /* This function calculates a damping factor for the Newton iteration update - * vector, dxneg, to insure that all site and bulk fractions, x, remain - * bounded between zero and one. + // This function calculates a damping factor for the Newton iteration update + // vector, dxneg, to insure that all solution components stay within perscribed 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 0-1, then Newton's method is allowed to operate normally. + * 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[], int dim, int *label) { - int i; - doublereal damp = 1.0, xnew, xtop, xbot; + doublereal damp = 1.0, xnew, xtop, xbot; static doublereal damp_old = 1.0; - *label = -1; - for (i = 0; i < dim; i++) { - + for (int i = 0; i < dim; i++) { + doublereal topBounds = m_topBounds[i]; + doublereal botBounds = m_botBounds[i]; /* * Calculate the new suggested new value of x[i] */ - // x_raw = x[i] - dxneg[i]; double delta_x = - dxneg[i]; xnew = x[i] - damp * dxneg[i]; /* * Calculate the allowed maximum and minimum values of x[i] - * - Only going to allow x[i] to converge to zero by a - * single order of magnitude at a time + * - Only going to allow x[i] to converge to the top and bottom bounds by a + * single order of magnitude at one time */ - xtop = 1.0 - 0.1*fabs(1.0-x[i]); - xbot = fabs(x[i]*0.1) - 1.0e-16; + xtop = topBounds - 0.1 * fabs(topBounds - x[i]); + + xbot = botBounds + 0.1 * fabs(x[i] - botBounds); + if (xnew > xtop) { - damp = - APPROACH * (1.0 - x[i]) / dxneg[i]; + damp = - APPROACH * (xtop - x[i]) / dxneg[i]; *label = i; } else if (xnew < xbot) { - damp = APPROACH * x[i] / dxneg[i]; + damp = APPROACH * (x[i] - xbot) / dxneg[i]; *label = i; } else if (xnew > 3.0*MAX(x[i], 1.0E-10)) { damp = - 2.0 * MAX(x[i], 1.0E-10) / dxneg[i]; @@ -527,7 +537,6 @@ namespace Cantera { } - // if (damp < 1.0e-2) damp = 1.0e-2; /* * Only allow the damping parameter to increase by a factor of three each * iteration. Heuristic to avoid oscillations in the value of damp @@ -541,7 +550,6 @@ namespace Cantera { * Save old value of the damping parameter for use * in subsequent calls. */ - damp_old = damp; return damp; @@ -628,7 +636,6 @@ namespace Cantera { } } - /* * Increase time step exponentially as same species repeatedly * controls time step @@ -658,7 +665,22 @@ namespace Cantera { 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 (int k = 0; k < m_neq; k++) { + m_botBounds[k] = botBounds[k]; + m_topBounds[k] = topBounds[k]; + } + } + //==================================================================================================================== /* * printResJac(): prints out the residual and Jacobian. * diff --git a/Cantera/src/numerics/solveProb.h b/Cantera/src/numerics/solveProb.h index 658fff0e3..1e7cf479c 100644 --- a/Cantera/src/numerics/solveProb.h +++ b/Cantera/src/numerics/solveProb.h @@ -66,12 +66,17 @@ namespace Cantera { //! Method to solve a pseudo steady state of a nonlinear problem /*! - * The following class handles solving nonlinear problem.s + * The following class handles the solution of a nonlinear problem. + * + * Res_ss(C) = - Res(C) = 0 * + * Optionally a pseudo transient algorithm may be used to relax the residual if + * it is available. * - * Note there are a couple of different types of species indecices - * floating around in the formulation of this object. - * + * Res_td(C) = dC/dt - Res(C) = 0; + * + * 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 @@ -329,8 +334,36 @@ namespace Cantera { const doublereal *CSolnSPOld, const bool do_time, const doublereal deltaT); + //! This function calculates a damping factor for the Newton iteration update + //! vector, dxneg, to insure that all solution components stay within perscribed 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. + */ virtual doublereal calc_damping(doublereal x[], doublereal dxneg[], int dim, int *label); + //! 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 + */ + virtual void setBounds(const doublereal botBounds[], const doublereal topBounds[]); + + //! residual function pointer to be solved. ResidEval *m_residFunc; //! Total number of equations to solve in the implicit problem. @@ -427,6 +460,18 @@ namespace Cantera { */ Array2D m_Jac; + //! Top bounds for the solution vector + /*! + * This defaults to 1.0 + */ + vector_fp m_topBounds; + + //! Bottom bounds for the solution vector + /*! + * This defaults to 0.0 + */ + vector_fp m_botBounds; + public: int m_ioflag;