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

Moved Tortuosity.h to include file structure so that 1Delectrode code can work.

Browse source 
It works and produces almost the same results as before.

Added DAE_solvers.cpp back in. There was a missing factory function.

Fixed bandsolver indexing calculations in NonlinearSolver and Bandmatrix due
   to int to size_t conversion. numerics test programs now work.
This commit is contained in:
Harry Moffat 2012-04-07 00:06:00 +00:00
parent 9c06d1e07c
commit 232b861284
7 changed files with 124 additions and 46 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

6
include/cantera/base/ct_defs.h
View file

@ -164,7 +164,13 @@ const doublereal MaxExp = 690.775527898;
//! Fairly random number to be used to initialize variables against
//! to see if they are subsequently defined.
const doublereal Undef = -999.1234;
//! Small number to compare differences of mole fractions against.
/*!
* This number is used for the interconversion of mole fraction and mass fraction quantities
* when the molecuar weight of a species is zero. It's also used for the matrix inversion
* of transport properties when mole fractions must be positive.
*/
const doublereal Tiny = 1.e-20;
//! inline function to return the max value of two doubles.

2
include/cantera/thermo/Phase.h
View file

@ -9,7 +9,7 @@
#include "cantera/base/vec_functions.h"
#include "cantera/base/ctml.h"
#include "Elements.h"
#include "cantera/thermo/Elements.h"
namespace Cantera
{

0
src/transport/Tortuosity.h → include/cantera/transport/Tortuosity.h
View file

47
src/numerics/BandMatrix.cpp
View file

@ -160,7 +160,9 @@ doublereal BandMatrix::value(size_t i, size_t j) const
//====================================================================================================================
size_t BandMatrix::index(size_t i, size_t j) const
{
size_t rw = m_kl + m_ku + i - j;
int jj = j;
int ii = i;
size_t rw = (int) m_kl + (int) m_ku + (int) ii - jj;
return (2*m_kl + m_ku + 1)*j + rw;
}
//====================================================================================================================
@ -218,12 +220,14 @@ vector_int& BandMatrix::ipiv()
*/
void BandMatrix::mult(const doublereal* b, doublereal* prod) const
{
size_t nr = nRows();
int kl = m_kl;
int ku = m_ku;
int nr = nRows();
doublereal sum = 0.0;
for (size_t m = 0; m < nr; m++) {
for (int m = 0; m < nr; m++) {
sum = 0.0;
for (size_t j = m - m_kl; j <= m + m_ku; j++) {
if (j < m_n) {
for (int j = m - kl; j <= m + ku; j++) {
if (j >= 0 && j < (int) m_n) {
sum += _value(m,j) * b[j];
}
}
@ -236,13 +240,16 @@ void BandMatrix::mult(const doublereal* b, doublereal* prod) const
*/
void BandMatrix::leftMult(const doublereal* const b, doublereal* const prod) const
{
size_t nc = nColumns();
int kl = m_kl;
int ku = m_ku;
int nc = nColumns();
doublereal sum = 0.0;
for (size_t n = 0; n < nc; n++) {
for (int n = 0; n < nc; n++) {
sum = 0.0;
for (size_t i = n - m_ku; i <= n + m_kl; i++) {
if (i < m_n) {
sum += _value(i,n) * b[i];
for (int i = n - ku; i <= n + kl; i++) {
if (i >= 0 && i < (int) m_n) {
size_t ii = i;
sum += _value(ii,n) * b[ii];
}
}
prod[n] = sum;
@ -421,12 +428,14 @@ int BandMatrix::factorAlgorithm() const
// Returns the one norm of the matrix
doublereal BandMatrix::oneNorm() const
{
int ku = m_ku;
int kl = m_kl;
doublereal value = 0.0;
for (size_t j = 0; j < m_n; j++) {
for (int j = 0; j < (int) m_n; j++) {
doublereal sum = 0.0;
doublereal* colP = m_colPtrs[j];
for (size_t i = j - m_ku; i <= j + m_kl; i++) {
sum += fabs(colP[m_kl + m_ku + i - j]);
for (int i = j - ku; i <= j + kl; i++) {
sum += fabs(colP[kl + ku + i - j]);
}
if (sum > value) {
value = sum;
@ -440,10 +449,10 @@ size_t BandMatrix::checkRows(doublereal& valueSmall) const
valueSmall = 1.0E300;
size_t iSmall = npos;
double vv;
for (size_t i = 0; i < m_n; i++) {
for (int i = 0; i < (int) m_n; i++) {
double valueS = 0.0;
for (size_t j = i - m_kl; j <= i + m_ku; j++) {
if (j < m_n) {
for (int j = i - (int) m_kl; j <= i + (int) m_ku; j++) {
if (j >= 0 && j < (int) m_n) {
vv = fabs(value(i,j));
if (vv > valueS) {
valueS = vv;
@ -466,10 +475,10 @@ size_t BandMatrix::checkColumns(doublereal& valueSmall) const
valueSmall = 1.0E300;
size_t jSmall = npos;
double vv;
for (size_t j = 0; j < m_n; j++) {
for (int j = 0; j < (int) m_n; j++) {
double valueS = 0.0;
for (size_t i = j - m_ku; i <= j + m_kl; i++) {
if (i < m_n) {
for (int i = j - (int) m_ku; i <= j + (int) m_kl; i++) {
if (i >= 0 && i < (int) m_n) {
vv = fabs(value(i,j));
if (vv > valueS) {
valueS = vv;

45
src/numerics/DAE_solvers.cpp Normal file
View file

@ -0,0 +1,45 @@
/**
* @file DAE_solvers.cpp
* Factory routine for picking the DAE solver package
*/
/*
* $Revision: 725 $
* $Date: 2011-05-16 18:45:08 -0600 (Mon, 16 May 2011) $
*/
/*
* Copyright 2004 Sandia Corporation. Under the terms of Contract
* DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
* retains certain rights in this software.
* See file License.txt for licensing information.
*/
#include "cantera/base/ct_defs.h"
#include "cantera/numerics/DAE_Solver.h"
#include "cantera/numerics/IDA_Solver.h"
// DAE_DEVEL is turned off at the current time
#define DAE_DEVEL
#ifdef DAE_DEVEL
namespace Cantera {
DAE_Solver* newDAE_Solver(std::string itype, ResidJacEval& f) {
if (itype == "IDA") {
#ifdef HAS_SUNDIALS
return new IDA_Solver(f);
#else
throw CanteraError("newDAE_Solver","IDA solver requires sundials"
" package, but Cantera was not built with sundials.");
#endif
}
else {
throw CanteraError("newDAE_Solver",
"unknown DAE solver: "+itype);
}
}
}
#
#endif

45
src/numerics/NonlinearSolver.cpp
View file

@ -754,7 +754,6 @@ void NonlinearSolver::scaleMatrix(GeneralMatrix& jac, doublereal* const y_comm,
doublereal time_curr, int num_newt_its)
{
size_t irow, jcol;
size_t ku, kl;
size_t ivec[2];
jac.nRowsAndStruct(ivec);
double* colP_j;
@ -783,12 +782,12 @@ void NonlinearSolver::scaleMatrix(GeneralMatrix& jac, doublereal* const y_comm,
}
}
} else if (jac.matrixType_ == 1) {
kl = ivec[0];
ku = ivec[1];
for (jcol = 0; jcol < neq_; jcol++) {
int kl = ivec[0];
int ku = ivec[1];
for (int jcol = 0; jcol < (int) neq_; jcol++) {
colP_j = (doublereal*) jac.ptrColumn(jcol);
for (irow = jcol - ku; irow <= jcol + kl; irow++) {
if (irow < neq_) {
for (int irow = jcol - ku; irow <= jcol + kl; irow++) {
if (irow >= 0 && irow < (int) neq_) {
colP_j[kl + ku + irow - jcol] *= m_colScales[jcol];
}
}
@ -828,12 +827,12 @@ void NonlinearSolver::scaleMatrix(GeneralMatrix& jac, doublereal* const y_comm,
}
}
} else if (jac.matrixType_ == 1) {
kl = ivec[0];
ku = ivec[1];
for (jcol = 0; jcol < neq_; jcol++) {
int kl = ivec[0];
int ku = ivec[1];
for (int jcol = 0; jcol < (int) neq_; jcol++) {
colP_j = (doublereal*) jac.ptrColumn(jcol);
for (irow = jcol - ku; irow <= jcol + kl; irow++) {
if (irow < neq_) {
for (int irow = jcol - ku; irow <= jcol + kl; irow++) {
if (irow >= 0 && irow < (int) neq_) {
double vv = fabs(colP_j[kl + ku + irow - jcol]);
if (m_rowScaling) {
m_rowScales[irow] += vv;
@ -871,12 +870,12 @@ void NonlinearSolver::scaleMatrix(GeneralMatrix& jac, doublereal* const y_comm,
}
}
} else if (jac.matrixType_ == 1) {
kl = ivec[0];
ku = ivec[1];
for (jcol = 0; jcol < neq_; jcol++) {
int kl = ivec[0];
int ku = ivec[1];
for (int jcol = 0; jcol < (int) neq_; jcol++) {
colP_j = (doublereal*) jac.ptrColumn(jcol);
for (irow = jcol - ku; irow <= jcol + kl; irow++) {
if (irow < neq_) {
for (int irow = jcol - ku; irow <= jcol + kl; irow++) {
if (irow >= 0 && irow < (int) neq_) {
colP_j[kl + ku + irow - jcol] *= m_rowScales[irow];
}
}
@ -3881,14 +3880,24 @@ int NonlinearSolver::beuler_jac(GeneralMatrix& J, doublereal* const f,
doublereal diff;
int ileft = (int) j - (int) ku;
int iright= j + kl;
for (int i = ileft; i <= iright; i++) {
if (i >= 0 && i < (int) neq_) {
size_t ii = i;
size_t index = (int) kl + (int) ku + i - (int) j;
diff = subtractRD(m_wksp[ii], f[ii]);
col_j[index] = diff / dy;
}
}
/*
for (size_t i = j - ku; i <= j + kl; i++) {
if (i < neq_) {
diff = subtractRD(m_wksp[i], f[i]);
col_j[kl + ku + i - j] = diff / dy;
}
}
*/
y[j] = ysave;
if (solnType_ != NSOLN_TYPE_STEADY_STATE) {
ydot[j] = ydotsave;

25
src/thermo/Phase.cpp
View file

@ -318,6 +318,9 @@ void Phase::setMoleFractions(const doublereal* const x)
{
// Use m_y as a temporary work vector for the non-negative mole fractions
doublereal norm = 0.0;
/*
* sum is calculated below as the unnormalized molecular weight
*/
doublereal sum = 0;
for (size_t k = 0; k < m_kk; k++) {
double xk = std::max(x[k], 0.0); // Ignore negative mole fractions
@ -325,13 +328,19 @@ void Phase::setMoleFractions(const doublereal* const x)
norm += xk;
sum += m_molwts[k] * xk;
}
transform(m_y.begin(), m_y.end(), m_ym.begin(),
timesConstant<double>(1.0/sum));
// Now set m_y to the mass fractions
transform(m_ym.begin(), m_ym.begin() + m_kk, m_molwts.begin(),
m_y.begin(), multiplies<double>());
/*
* Set m_ym_ to the normalized mole fractions divided by the normalized mean molecular weight:
* m_ym_k = X_k / (sum_k X_k M_k)
*/
transform(m_y.begin(), m_y.end(), m_ym.begin(), timesConstant<double>(1.0/sum));
/*
* Now set m_y to the normalized mass fractions
* m_y = X_k M_k / (sum_k X_k M_k)
*/
transform(m_ym.begin(), m_ym.begin() + m_kk, m_molwts.begin(), m_y.begin(), multiplies<double>());
/*
* Calculate the normalized molecular weight
*/
m_mmw = sum/norm;
// Call a routine to determine whether state has changed.
@ -961,7 +970,7 @@ void Phase::init(const vector_fp& mw)
}
// Some surface phases may define species representing empty sites
// hat have zero molecular weight. Give them a very small molecular
// that have zero molecular weight. Give them a very small molecular
// weight to avoid dividing by zero.
if (m_molwts[k] < Tiny) {
m_molwts[k] = Tiny;