From d5f44a5e6b959dc8a20effb07cb252e115b3a7fd Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Tue, 7 Dec 2010 18:05:25 +0000 Subject: [PATCH] ResidJacEval: Added in support for returning a failure condition from calcDeltaSolnVariables() NonLinearSolver: Added in support for returning a failure condition from bueler_jac(). Added in a call to return an editable bounds constraint vector. DenseMatrix: Added in comments about error handling from LAPACK routines. mdp_allo: Added in a new routine, mdp_zero_int_1 thermo.h electrolyteThermo.h Started to add in ifdef blocks. These are not handled correctly currently. --- Cantera/cxx/include/electrolyteThermo.h | 4 ++- Cantera/cxx/include/thermo.h | 29 +++++++++++++++++ Cantera/src/base/mdp_allo.cpp | 24 ++++++++++++++ Cantera/src/base/mdp_allo.h | 7 +++++ Cantera/src/numerics/DenseMatrix.h | 31 ++++++++++++++++++- Cantera/src/numerics/NonlinearSolver.cpp | 29 ++++++++++++----- Cantera/src/numerics/NonlinearSolver.h | 20 ++++++++++-- Cantera/src/numerics/ResidJacEval.cpp | 7 ++++- Cantera/src/numerics/ResidJacEval.h | 12 +++++-- docs/Cantera.cfg.in | 4 ++- .../linux.64_sierra_gcc444_python264_numpy | 3 +- 11 files changed, 153 insertions(+), 17 deletions(-) diff --git a/Cantera/cxx/include/electrolyteThermo.h b/Cantera/cxx/include/electrolyteThermo.h index cdd1db33c..aef557f24 100644 --- a/Cantera/cxx/include/electrolyteThermo.h +++ b/Cantera/cxx/include/electrolyteThermo.h @@ -9,8 +9,8 @@ #ifndef CT_ELECTROLYTETHERMO_INCL #define CT_ELECTROLYTETHERMO_INCL -#include "thermo.h" +#ifdef WITH_ELECTROLYTES #include "kernel/electrolytes.h" #include "kernel/MolalityVPSSTP.h" #include "kernel/VPStandardStateTP.h" @@ -26,3 +26,5 @@ #include "kernel/VPSSMgr_Water_HKFT.h" #include "kernel/VPSSMgr_Water_ConstVol.h" #endif + +#endif diff --git a/Cantera/cxx/include/thermo.h b/Cantera/cxx/include/thermo.h index a27b23336..5eff5eaea 100755 --- a/Cantera/cxx/include/thermo.h +++ b/Cantera/cxx/include/thermo.h @@ -14,4 +14,33 @@ #include "kernel/SurfPhase.h" #include "kernel/EdgePhase.h" + +#ifdef WITH_IDEAL_SOLUTIONS + +#include "kernel/GibbsExcessVPSSTP.h" +#include "kernel/MargulesVPSSTP.h" + +#endif + +#ifdef WITH_ELECTROLYTES + +#include "electrolyteThermo.h" + +#endif + + +#ifdef WITH_LATTICE_SOLID + +#include "kernel/LatticePhase.h" +#include "kernel/LatticeSolidPhase.h" + +#endif + +#ifdef WITH_PURE_FLUIDS + + +#endif + + + #endif diff --git a/Cantera/src/base/mdp_allo.cpp b/Cantera/src/base/mdp_allo.cpp index d650ed919..c1308049a 100644 --- a/Cantera/src/base/mdp_allo.cpp +++ b/Cantera/src/base/mdp_allo.cpp @@ -1697,6 +1697,30 @@ namespace mdp { (void) memset((void *)v, 0, len * sizeof(double)); } } + + /****************************************************************************/ + /****************************************************************************/ + /****************************************************************************/ + + void mdp_zero_int_1(int * const v, const int len) + + /************************************************************************** + * + * mdp_zero_int_1: + * + * Zeroes out an int vector + * + * Input + * ------------- + * v = Vector of values to be set to zero + * len = Length of the vector + **************************************************************************/ + { + if (len > 0) { + (void) memset((void *)v, 0, len * sizeof(int)); + } + } + /****************************************************************************/ /****************************************************************************/ /****************************************************************************/ diff --git a/Cantera/src/base/mdp_allo.h b/Cantera/src/base/mdp_allo.h index 756b6dce4..1b196bf3f 100644 --- a/Cantera/src/base/mdp_allo.h +++ b/Cantera/src/base/mdp_allo.h @@ -674,6 +674,13 @@ namespace mdp { */ extern void mdp_zero_dbl_1(double * const v , const int len); + //! Zeroes an int vector + /*! + * @param v = Vector of values to be assigned + * @param len = Length of the vector + */ + extern void mdp_zero_int_1(int * const v , const int len); + //! Assigns a single value to a double matrix. Contiguous data for the //! matrix is assumed. /*! diff --git a/Cantera/src/numerics/DenseMatrix.h b/Cantera/src/numerics/DenseMatrix.h index c9d38af47..04aac8f48 100644 --- a/Cantera/src/numerics/DenseMatrix.h +++ b/Cantera/src/numerics/DenseMatrix.h @@ -33,6 +33,10 @@ namespace Cantera { //! Exception thrown when an LAPACK error is encountered associated with inverting or solving a matrix + /*! + * A named error condition is used so that the calling code may differentiate this type of error + * from other error conditions. + */ class CELapackError : public CanteraError { public: @@ -48,12 +52,25 @@ namespace Cantera { }; - //! A class for full (non-sparse) matrices with Fortran-compatible + //! A class for full (non-sparse) matrices with Fortran-compatible //! data storage, which adds matrix operations to class Array2D. /*! * The dense matrix class adds matrix operations onto the Array2D class. * These matrix operations are carried out by the appropriate BLAS and LAPACK routines * + * Error handling from BLAS and LAPACK are handled via the following formulation. + * Depending on a variable, a singular matrix or other terminal error condition from + * LAPACK is handled by either throwing an exception of type, CELapackError, or by + * returning the error code condition to the calling routine. + * + * The int variable, m_useReturnErrorCode, determines which method is used. + * The default value of zero means that an exception is thrown. A value of 1 + * means that a return code is used. + * + * Reporting of these LAPACK error conditions is handled by the class variable + * m_printLevel. The default is for no reporting. If m_printLevel is nonzero, + * the error condition is reported to Cantera's log file. + * * @ingroup numerics */ class DenseMatrix : public Array2D { @@ -134,6 +151,7 @@ namespace Cantera { vector_int m_ipiv; public: + //! Error Handling Flag /*! * The default is to set this to 0. In this case, if a factorization is requested and can't be achieved, @@ -150,6 +168,17 @@ namespace Cantera { * Level of printing that is carried out. Only error conditions are printed out, if this value is nonzero. */ int m_printLevel; + + + // Listing of friend functions which are defined below + + friend int solve(DenseMatrix& A, double* b); + friend int solve(DenseMatrix& A, DenseMatrix& b); + friend int invert(DenseMatrix& A, int nn); +#ifdef INCL_LEAST_SQUARES + friend int leastSquares(DenseMatrix& A, double* b); +#endif + }; //================================================================================================================== diff --git a/Cantera/src/numerics/NonlinearSolver.cpp b/Cantera/src/numerics/NonlinearSolver.cpp index c0e65774c..b8df9939d 100644 --- a/Cantera/src/numerics/NonlinearSolver.cpp +++ b/Cantera/src/numerics/NonlinearSolver.cpp @@ -281,6 +281,14 @@ namespace Cantera { } } //==================================================================================================================== + 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 @@ -1224,8 +1232,11 @@ namespace Cantera { if (m_print_flag > 3) { printf("\tsolve_nonlinear_problem(): Getting a new Jacobian and solving system\n"); } - beuler_jac(jac, DATA_PTR(m_resid), time_curr, CJ, DATA_PTR(m_y_n), DATA_PTR(ydot_curr), - num_newt_its); + info = beuler_jac(jac, DATA_PTR(m_resid), time_curr, CJ, DATA_PTR(m_y_n), DATA_PTR(ydot_curr), num_newt_its); + if (info == 0) { + m = -1; + goto done; + } m_residCurrent = true; } else { if (m_print_flag > 1) { @@ -1526,9 +1537,12 @@ namespace Cantera { * @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 */ - void NonlinearSolver::beuler_jac(SquareMatrix &J, doublereal * const f, + int NonlinearSolver::beuler_jac(SquareMatrix &J, doublereal * const f, doublereal time_curr, doublereal CJ, doublereal * const y, doublereal * const ydot, @@ -1537,7 +1551,7 @@ namespace Cantera { int i, j; double* col_j; doublereal ysave, ydotsave, dy; - + int retn = 1; /* * Clear the factor flag */ @@ -1568,7 +1582,8 @@ namespace Cantera { * derivative. */ doublereal *dyVector = mdp::mdp_alloc_dbl_1(neq_, MDP_DBL_NOINIT); - m_func->calcDeltaSolnVariables(time_curr, y, ydot, dyVector, DATA_PTR(m_ewt)); + retn = m_func->calcDeltaSolnVariables(time_curr, y, ydot, dyVector, DATA_PTR(m_ewt)); + if (s_print_NumJac) { @@ -1668,7 +1683,7 @@ namespace Cantera { printf("--------------\n"); } } - + return retn; } //==================================================================================================================== // Internal function to calculate the time derivative at the new step diff --git a/Cantera/src/numerics/NonlinearSolver.h b/Cantera/src/numerics/NonlinearSolver.h index 21062cc8e..08699226e 100644 --- a/Cantera/src/numerics/NonlinearSolver.h +++ b/Cantera/src/numerics/NonlinearSolver.h @@ -229,6 +229,11 @@ namespace Cantera { void setBoundsConstraints(const doublereal * const y_low_bounds, const doublereal * const y_high_bounds); + //! return an editable vector of the low bounds constraints + std::vector & lowBoundsConstraintVector(); + + //! return an editable vector of the high bounds constraints + std::vector & highBoundsConstraintVector(); //! Internal function to calculate the time derivative at the new step /*! @@ -238,13 +243,22 @@ namespace Cantera { */ void calc_ydot(const int order, const doublereal * const y_curr, doublereal * const ydot_curr); - //! Function called to evaluate the jacobian matrix and the curent - //! residual vector. + //! Function called to evaluate the jacobian matrix and the current + //! residual vector at the current time step /*! + * * + * @param J = Jacobian matrix to be filled in + * @param f = Right hand side. This routine returns the current + * value of the rhs (output), so that it does + * not have to be computed again. + * + * @return Returns a flag to indicate that operation is successful. + * 1 Means a successful operation + * 0 Means an unsuccessful operation * */ - void beuler_jac(SquareMatrix &J, doublereal * const f, + int beuler_jac(SquareMatrix &J, doublereal * const f, doublereal time_curr, doublereal CJ, doublereal * const y, doublereal * const ydot, int num_newt_its); diff --git a/Cantera/src/numerics/ResidJacEval.cpp b/Cantera/src/numerics/ResidJacEval.cpp index e27aa8db9..b1ae4edf1 100644 --- a/Cantera/src/numerics/ResidJacEval.cpp +++ b/Cantera/src/numerics/ResidJacEval.cpp @@ -172,8 +172,12 @@ namespace Cantera { * @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 */ - void ResidJacEval:: + int ResidJacEval:: calcDeltaSolnVariables(const doublereal t, const doublereal * const ySoln, const doublereal * const ySolnDot, doublereal * const deltaYSoln, const doublereal *const solnWeights) @@ -187,6 +191,7 @@ namespace Cantera { deltaYSoln[i] = fmaxx(1.0E-2 * solnWeights[i], 1.0E-6 * fabs(ySoln[i])); } } + return 1; } //==================================================================================================================== // Returns a vector of column scale factors that can be used to column scale Jacobians. diff --git a/Cantera/src/numerics/ResidJacEval.h b/Cantera/src/numerics/ResidJacEval.h index 7ddb4e0a9..e1b1fc069 100644 --- a/Cantera/src/numerics/ResidJacEval.h +++ b/Cantera/src/numerics/ResidJacEval.h @@ -138,8 +138,10 @@ namespace Cantera { * solution vector eliminating illegal directions. * * @param t Time (input) - * @param y Solution vector (input, output) + * @param ybase Solution vector (input, output) * @param step Proposed step in the solution that will be cropped + * + * @return Return the norm of the amount of filtering */ virtual doublereal filterNewStep(const doublereal t, const doublereal * const ybase, doublereal * const step); @@ -151,6 +153,8 @@ namespace Cantera { * * @param t Time (input) * @param y Solution vector (input, output) + * + * @return Return the norm of the amount of filtering */ virtual doublereal filterSolnPrediction(const doublereal t, doublereal * const y); @@ -205,8 +209,12 @@ namespace Cantera { * @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 or neg value Means an unsuccessful operation */ - virtual void + virtual int calcDeltaSolnVariables(const doublereal t, const doublereal * const y, const doublereal * const ydot, diff --git a/docs/Cantera.cfg.in b/docs/Cantera.cfg.in index e73623793..a732ce226 100755 --- a/docs/Cantera.cfg.in +++ b/docs/Cantera.cfg.in @@ -792,7 +792,9 @@ FILE_PATTERNS = Kinetics.h Kinetics.cpp \ WaterTransport.h \ WaterTransport.cpp \ DenseMatrix.h \ - DenseMatrix.cpp + DenseMatrix.cpp \ + ResidJacEval.h \ + ResidJacEval.cpp # The RECURSIVE tag can be used to turn specify whether or not subdirectories # should be searched for input files as well. Possible values are YES and NO. diff --git a/docs/install_examples/linux.64_sierra_gcc444_python264_numpy b/docs/install_examples/linux.64_sierra_gcc444_python264_numpy index dc66fcce8..30e9e8181 100755 --- a/docs/install_examples/linux.64_sierra_gcc444_python264_numpy +++ b/docs/install_examples/linux.64_sierra_gcc444_python264_numpy @@ -86,7 +86,8 @@ export F77 CFLAGS="-g -Wall" export CFLAGS -CXXFLAGS="-g -Wall -Woverloaded-virtual -DDEBUG_HKM -DDEBUG_HKM_EPEQUIL" +#CXXFLAGS="-g -Wall -Woverloaded-virtual -DDEBUG_HKM -DDEBUG_HKM_EPEQUIL -DDEBUG_NUMJAC" +CXXFLAGS="-g -Wall -Woverloaded-virtual -DDEBUG_HKM -DDEBUG_HKM_EPEQUIL " export CXXFLAGS FFLAGS="-g -DDEBUG_HKM -fno-second-underscore"