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

Eliminate the use of SquareMatrix in favor of DenseMatrix

Browse source
This commit is contained in:
Ray Speth 2016-04-10 21:25:55 -04:00
parent b49c2e4c2d
commit 961ca06cd5
4 changed files with 10 additions and 60 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

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

@ -13,7 +13,7 @@
#define SOLVESP_H
#include "cantera/kinetics/InterfaceKinetics.h"
#include "cantera/numerics/SquareMatrix.h"
#include "cantera/numerics/DenseMatrix.h"
//! @defgroup solvesp_methods Surface Problem Solver Methods
//! @{
@ -301,7 +301,7 @@ private:
* @param do_time Calculate a time dependent residual
* @param deltaT Delta time for time dependent problem.
*/
void resjac_eval(SquareMatrix& jac, doublereal* resid,
void resjac_eval(DenseMatrix& jac, doublereal* resid,
doublereal* CSolnSP,
const doublereal* CSolnSPOld, const bool do_time,
const doublereal deltaT);
@ -504,7 +504,7 @@ private:
//! Jacobian. m_neq by m_neq computed Jacobian matrix for the local
//! Newton's method.
SquareMatrix m_Jac;
DenseMatrix m_Jac;
public:
int m_ioflag;

3
include/cantera/transport/MultiTransport.h
View file

@ -10,7 +10,6 @@
// Cantera includes
#include "GasTransport.h"
#include "cantera/numerics/SquareMatrix.h"
namespace Cantera
{
@ -151,7 +150,7 @@ protected:
// L matrix quantities
DenseMatrix m_Lmatrix;
SquareMatrix m_aa;
DenseMatrix m_aa;
vector_fp m_a;
vector_fp m_b;

37
src/kinetics/solveSP.cpp
View file

@ -223,39 +223,8 @@ int solveSP::solveSurfProb(int ifunc, doublereal time_scale, doublereal TKelvin,
// Find the weighted norm of the residual
double resid_norm = calcWeightedNorm(m_wtResid.data(), m_resid.data(), m_neq);
// Solve Linear system. The solution is in resid[]
int info = m_Jac.factor();
if (info==0) {
m_Jac.solve(&m_resid[0]);
} else {
// Force convergence if residual is small to avoid "nan" results
// from the linear solve.
if (m_ioflag) {
writelogf("solveSurfSS: Zero pivot, assuming converged: %g (%d)\n",
resid_norm, info);
}
for (size_t jcol = 0; jcol < m_neq; jcol++) {
m_resid[jcol] = 0.0;
}
// print out some helpful info
if (m_ioflag > 1) {
writelog("-----\n");
writelogf("solveSurfProb: iter %d t_real %g delta_t %g\n\n",
iter,t_real, 1.0/inv_t);
writelog("solveSurfProb: init guess, current concentration,"
"and prod rate:\n");
for (size_t jcol = 0; jcol < m_neq; jcol++) {
writelog("\t%d %g %g %g\n", jcol,
m_CSolnSPInit[jcol], m_CSolnSP[jcol],
m_netProductionRatesSave[m_kinSpecIndex[jcol]]);
}
writelog("-----\n");
}
if (do_time) {
t_real += time_scale;
}
}
// Solve Linear system. The solution is in m_resid
solve(m_Jac, m_resid.data());
// Calculate the Damping factor needed to keep all unknowns between 0
// and 1, and not allow too large a change (factor of 2) in any unknown.
@ -476,7 +445,7 @@ void solveSP::fun_eval(doublereal* resid, const doublereal* CSoln,
}
}
void solveSP::resjac_eval(SquareMatrix& jac,
void solveSP::resjac_eval(DenseMatrix& jac,
doublereal resid[], doublereal CSoln[],
const doublereal CSolnOld[], const bool do_time,
const doublereal deltaT)

24
src/transport/MultiTransport.cpp
View file

@ -226,17 +226,8 @@ void MultiTransport::getSpeciesFluxes(size_t ndim, const doublereal* const grad_
fluxes[jmax + n*ldf] = 0.0;
}
// use LAPACK to solve the equations
int info = m_aa.factor();
if (info) {
throw CanteraError("MultiTransport::getSpeciesFluxes",
"Error factorizing matrix.");
}
info = m_aa.solve(fluxes, ndim, ldf);
if (info) {
throw CanteraError("MultiTransport::getSpeciesFluxes",
"Error solving linear system.");
}
// solve the equations
solve(m_aa, fluxes, ndim, ldf);
size_t offset;
doublereal pp = pressure_ig();
@ -332,16 +323,7 @@ void MultiTransport::getMassFluxes(const doublereal* state1, const doublereal* s
fluxes[jmax] = 0.0;
// Solve the equations
int info = m_aa.factor();
if (info) {
throw CanteraError("MultiTransport::getMassFluxes",
"Error in factorization. Info = {}", info);
}
info = m_aa.solve(fluxes);
if (info) {
throw CanteraError("MultiTransport::getMassFluxes",
"Error in linear solve. Info = {}", info);
}
solve(m_aa, fluxes);
doublereal pp = pressure_ig();
// multiply diffusion velocities by rho * Y_k to create