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

[Kinetics] Abstract use of LAPACK in solveSP and solveProb

Browse source
This commit is contained in:
Ray Speth 2014-03-24 04:02:18 +00:00
parent 8933203a20
commit 44c53a6f27
4 changed files with 9 additions and 30 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

7
include/cantera/kinetics/solveSP.h
View file

@ -14,7 +14,7 @@
#include "ImplicitSurfChem.h"
#include "cantera/kinetics/InterfaceKinetics.h"
#include "cantera/base/Array.h"
#include "cantera/numerics/SquareMatrix.h"
//! @defgroup solvesp_methods Surface Problem Solver Methods
//! @{
@ -533,9 +533,6 @@ private:
//! Vector of mole fractions. length m_maxTotSpecies
vector_fp m_XMolKinSpecies;
//! pivots. length MAX(1, m_neq)
vector_int m_ipiv;
//! Vector of pointers to the top of the columns of the Jacobian
/*!
* The "dim" by "dim" computed Jacobian matrix for the
@ -545,7 +542,7 @@ private:
//! Jacobian. m_neq by m_neq computed Jacobian matrix for the local
//! Newton's method.
Array2D m_Jac;
SquareMatrix m_Jac;
public:
int m_ioflag;

10
include/cantera/numerics/solveProb.h
View file

@ -17,7 +17,7 @@
* @defgroup solverGroup Solvers for Equation Systems
*/
#include "cantera/base/Array.h"
#include "cantera/numerics/SquareMatrix.h"
#include "ResidEval.h"
//! Solution Methods
@ -415,12 +415,6 @@ private:
*/
vector_fp m_resid;
//! pivots
/*!
* length MAX(1, m_neq)
*/
vector_int m_ipiv;
//! Vector of pointers to the top of the columns of the jacobians
/*!
* The "dim" by "dim" computed Jacobian matrix for the
@ -432,7 +426,7 @@ private:
/*!
* m_neq by m_neq computed Jacobian matrix for the local Newton's method.
*/
Array2D m_Jac;
SquareMatrix m_Jac;
//! Top bounds for the solution vector
/*!

12
src/kinetics/solveSP.cpp
View file

@ -10,7 +10,6 @@
#include "cantera/kinetics/solveSP.h"
#include "cantera/base/clockWC.h"
#include "cantera/numerics/ctlapack.h"
#include <cstdio>
@ -125,8 +124,6 @@ solveSP::solveSP(ImplicitSurfChem* surfChemPtr, int bulkFunc) :
m_wtResid.resize(dim1, 0.0);
m_wtSpecies.resize(dim1, 0.0);
m_resid.resize(dim1, 0.0);
m_ipiv.resize(dim1, 0);
m_Jac.resize(dim1, dim1, 0.0);
m_JacCol.resize(dim1, 0);
for (size_t k = 0; k < dim1; k++) {
@ -152,7 +149,6 @@ int solveSP::solveSurfProb(int ifunc, doublereal time_scale, doublereal TKelvin,
doublereal label_factor = 1.0;
int iter=0; // iteration number on numlinear solver
int iter_max=1000; // maximum number of nonlinear iterations
int nrhs=1;
doublereal deltaT = 1.0E-10; // Delta time step
doublereal damp=1.0, tmp;
// Weighted L2 norm of the residual. Currently, this is only
@ -291,13 +287,11 @@ int solveSP::solveSurfProb(int ifunc, doublereal time_scale, doublereal TKelvin,
DATA_PTR(m_resid), m_neq);
/*
* Solve Linear system (with LAPACK). The solution is in resid[]
* Solve Linear system. The solution is in resid[]
*/
ct_dgetrf(m_neq, m_neq, m_JacCol[0], m_neq, DATA_PTR(m_ipiv), info);
info = m_Jac.factor();
if (info==0) {
ct_dgetrs(ctlapack::NoTranspose, m_neq, nrhs, m_JacCol[0],
m_neq, DATA_PTR(m_ipiv), DATA_PTR(m_resid), m_neq,
info);
m_Jac.solve(&m_resid[0]);
}
/*
* Force convergence if residual is small to avoid

10
src/numerics/solveProb.cpp
View file

@ -11,7 +11,6 @@
#include "cantera/numerics/solveProb.h"
#include "cantera/base/clockWC.h"
#include "cantera/numerics/ctlapack.h"
#include "cantera/base/stringUtils.h"
using namespace std;
@ -48,7 +47,6 @@ solveProb::solveProb(ResidEval* resid) :
m_wtResid.resize(dim1, 0.0);
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);
@ -84,7 +82,6 @@ int solveProb::solve(int ifunc, doublereal time_scale,
doublereal label_factor = 1.0;
int iter=0; // iteration number on numlinear solver
int iter_max=1000; // maximum number of nonlinear iterations
int nrhs=1;
doublereal deltaT = 1.0E-10; // Delta time step
doublereal damp=1.0, tmp;
// Weighted L2 norm of the residual. Currently, this is only
@ -238,12 +235,9 @@ int solveProb::solve(int ifunc, doublereal time_scale,
/*
* Solve Linear system (with LAPACK). The solution is in resid[]
*/
ct_dgetrf(m_neq, m_neq, m_JacCol[0], m_neq, DATA_PTR(m_ipiv), info);
info = m_Jac.factor();
if (info==0) {
ct_dgetrs(ctlapack::NoTranspose, m_neq, nrhs, m_JacCol[0],
m_neq, DATA_PTR(m_ipiv), DATA_PTR(m_resid), m_neq,
info);
m_Jac.solve(&m_resid[0]);
}
/*
* Force convergence if residual is small to avoid