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

ResidJacEval:

Browse source 
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.
This commit is contained in:
Harry Moffat 2010-12-07 18:05:25 +00:00
parent 44e03fcf79
commit d5f44a5e6b
11 changed files with 153 additions and 17 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

4
Cantera/cxx/include/electrolyteThermo.h
View file

@ -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

29
Cantera/cxx/include/thermo.h
View file

@ -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

24
Cantera/src/base/mdp_allo.cpp
View file

@ -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));
}
}
/****************************************************************************/
/****************************************************************************/
/****************************************************************************/

7
Cantera/src/base/mdp_allo.h
View file

@ -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.
/*!

31
Cantera/src/numerics/DenseMatrix.h
View file

@ -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
};
//==================================================================================================================

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

@ -281,6 +281,14 @@ namespace Cantera {
}
}
//====================================================================================================================
std::vector<double> & NonlinearSolver::lowBoundsConstraintVector() {
return m_y_low_bounds;
}
//====================================================================================================================
std::vector<double> & 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

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

@ -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<double> & lowBoundsConstraintVector();
//! return an editable vector of the high bounds constraints
std::vector<double> & 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);

7
Cantera/src/numerics/ResidJacEval.cpp
View file

@ -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.

12
Cantera/src/numerics/ResidJacEval.h
View file

@ -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,

4
docs/Cantera.cfg.in
View file

@ -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.

3
docs/install_examples/linux.64_sierra_gcc444_python264_numpy
View file

@ -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"