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

modifications to simplify solving boundary value problems

Browse source
This commit is contained in:
Dave Goodwin 2005-09-13 17:10:08 +00:00
parent 45c5ab1f98
commit 7718361c94
5 changed files with 245 additions and 65 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

190
Cantera/src/oneD/Domain1D.cpp Normal file
View file

@ -0,0 +1,190 @@
/**
* @file Domain1D.cpp
*
*/
#include "Domain1D.h"
namespace Cantera {
void Domain1D::
setTolerances(int nr, const doublereal* rtol,
int na, const doublereal* atol, int ts) {
if (nr < m_nv || na < m_nv)
throw CanteraError("Domain1D::setTolerances",
"wrong array size for solution error tolerances. "
"Size should be at least "+int2str(m_nv));
if (ts >= 0) {
copy(rtol, rtol + m_nv, m_rtol_ss.begin());
copy(atol, atol + m_nv, m_atol_ss.begin());
}
if (ts <= 0) {
copy(rtol, rtol + m_nv, m_rtol_ts.begin());
copy(atol, atol + m_nv, m_atol_ts.begin());
}
}
void Domain1D::
setTolerances(int n, doublereal rtol, doublereal atol, int ts) {
if (ts >= 0) {
m_rtol_ss[n] = rtol;
m_atol_ss[n] = atol;
}
if (ts <= 0) {
m_rtol_ts[n] = rtol;
m_atol_ts[n] = atol;
}
}
void Domain1D::
setTolerances(doublereal rtol, doublereal atol,int ts) {
for (int n = 0; n < m_nv; n++){
if(ts >= 0) {
m_rtol_ss[n] = rtol;
m_atol_ss[n] = atol;
}
if (ts <= 0) {
m_rtol_ts[n] = rtol;
m_atol_ts[n] = atol;
}
}
}
void Domain1D::
setTolerancesTS(doublereal rtol, doublereal atol) {
for (int n = 0; n < m_nv; n++){
m_rtol_ts[n] = rtol;
m_atol_ts[n] = atol;
}
}
void Domain1D::
setTolerancesSS(doublereal rtol, doublereal atol) {
for (int n = 0; n < m_nv; n++){
m_rtol_ss[n] = rtol;
m_atol_ss[n] = atol;
}
}
void Domain1D::
eval(int jg, doublereal* xg, doublereal* rg,
integer* mask, doublereal rdt) {
if (jg >=0 && (jg < firstPoint() - 1 || jg > lastPoint() + 1)) return;
// if evaluating a Jacobian, compute the steady-state residual
if (jg >= 0) rdt = 0.0;
// start of local part of global arrays
doublereal* x = xg + loc();
doublereal* rsd = rg + loc();
integer* diag = mask + loc();
int jmin, jmax, jpt, j, i;
jpt = jg - firstPoint();
if (jg < 0) { // evaluate all points
jmin = 0;
jmax = m_points - 1;
}
else { // evaluate points for Jacobian
jmin = max(jpt-1, 0);
jmax = min(jpt+1,m_points-1);
}
for (j = jmin; j <= jmax; j++) {
if (j == 0 || j == m_points - 1) {
for (i = 0; i < m_nv; i++) {
rsd[index(i,j)] = residual(x,i,j);
diag[index(i,j)] = 0;
}
}
else {
for (i = 0; i < m_nv; i++) {
rsd[index(i,j)] = residual(x,i,j)
- timeDerivativeFlag(i)*rdt*(value(x,i,j) - prevSoln(i,j));
diag[index(i,j)] = timeDerivativeFlag(i);
}
}
}
}
// called to set up initial grid, and after grid refinement
void Domain1D::setupGrid(int n, const doublereal* z) {
resize(m_nv, n);
int j;
for (j = 0; j < m_points; j++) m_z[j] = z[j];
}
void drawline() {
writelog("\n-------------------------------------"
"------------------------------------------");
}
/**
* Print the solution.
*/
void Domain1D::showSolution(const doublereal* x) {
int nn = m_nv/5;
int i, j, n;
//char* buf = new char[100];
char buf[100];
doublereal v;
for (i = 0; i < nn; i++) {
drawline();
sprintf(buf, "\n z ");
writelog(buf);
for (n = 0; n < 5; n++) {
sprintf(buf, " %10s ",componentName(i*5 + n).c_str());
writelog(buf);
}
drawline();
for (j = 0; j < m_points; j++) {
sprintf(buf, "\n %10.4g ",m_z[j]);
writelog(buf);
for (n = 0; n < 5; n++) {
v = value(x, i*5+n, j);
sprintf(buf, " %10.4g ",v);
writelog(buf);
}
}
writelog("\n");
}
int nrem = m_nv - 5*nn;
drawline();
sprintf(buf, "\n z ");
writelog(buf);
for (n = 0; n < nrem; n++) {
sprintf(buf, " %10s ", componentName(nn*5 + n).c_str());
writelog(buf);
}
drawline();
for (j = 0; j < m_points; j++) {
sprintf(buf, "\n %10.4g ",m_z[j]);
writelog(buf);
for (n = 0; n < nrem; n++) {
v = value(x, nn*5+n, j);
sprintf(buf, " %10.4g ", v);
writelog(buf);
}
}
writelog("\n");
}
// initial solution
void Domain1D::_getInitialSoln(doublereal* x) {
for (int j = 0; j < m_points; j++) {
for (int n = 0; n < m_nv; n++) {
x[index(n,j)] = initialValue(n,j);
}
}
}
doublereal Domain1D::initialValue(int n, int j) {
throw CanteraError("Domain1D::initialValue",
"base class method called!");
}
} // namespace

107
Cantera/src/oneD/Domain1D.h
View file

@ -102,12 +102,16 @@ namespace Cantera {
* actions required to resize the domain.
*/
virtual void resize(int nv, int np) {
// if the number of components is being changed, then a
// new grid refiner is required.
if (nv != m_nv || !m_refiner) {
m_nv = nv;
delete m_refiner;
m_refiner = new Refiner(*this);
}
m_nv = nv;
m_td.resize(m_nv, 1);
m_name.resize(m_nv,"");
m_max.resize(m_nv, 0.0);
m_min.resize(m_nv, 0.0);
m_rtol_ss.resize(m_nv, 1.0e-8);
@ -120,6 +124,7 @@ namespace Cantera {
locate();
}
/// Return a reference to the grid refiner.
Refiner& refiner() { return *m_refiner; }
/// Number of components at each grid point.
@ -130,9 +135,19 @@ namespace Cantera {
/// Name of the nth component. May be overloaded.
virtual string componentName(int n) const {
return "component " + int2str(n);
if (m_name[n] != "") return m_name[n];
else return "component " + int2str(n);
}
void setComponentName(int n, string name) {
m_name[n] = name;
}
void setComponentType(int n, int ctype) {
if (ctype == 0) setAlgebraic(n);
}
/// index of component with name \a name.
int componentIndex(string name) {
int nc = nComponents();
for (int n = 0; n < nc; n++) {
@ -160,60 +175,23 @@ namespace Cantera {
m_max[n] = upper;
}
/// set the error tolerances for all solution components.
void setTolerances(int nr, const doublereal* rtol,
int na, const doublereal* atol, int ts = 0) {
if (nr < m_nv || na < m_nv)
throw CanteraError("Domain1D::setTolerances",
"wrong array size for solution error tolerances. "
"Size should be at least "+int2str(m_nv));
if (ts >= 0) {
copy(rtol, rtol + m_nv, m_rtol_ss.begin());
copy(atol, atol + m_nv, m_atol_ss.begin());
}
if (ts <= 0) {
copy(rtol, rtol + m_nv, m_rtol_ts.begin());
copy(atol, atol + m_nv, m_atol_ts.begin());
}
}
int na, const doublereal* atol, int ts = 0);
void setTolerances(int n, doublereal rtol, doublereal atol, int ts = 0) {
if (ts >= 0) {
m_rtol_ss[n] = rtol;
m_atol_ss[n] = atol;
}
if (ts <= 0) {
m_rtol_ts[n] = rtol;
m_atol_ts[n] = atol;
}
}
/// set the error tolerances for solution component \a n.
void setTolerances(int n, doublereal rtol, doublereal atol, int ts = 0);
//added by Karl Meredith
void setTolerances(doublereal rtol, doublereal atol,int ts=0) {
for (int n=0;n<m_nv;n++){
if(ts>=0) {
m_rtol_ss[n] = rtol;
m_atol_ss[n] = atol;
}
if (ts <= 0) {
m_rtol_ts[n] = rtol;
m_atol_ts[n] = atol;
}
}
}
/// set scalar error tolerances. All solution components will
/// have the same relative and absolute error tolerances.
void setTolerances(doublereal rtol, doublereal atol,int ts=0);
//added by Karl Meredith
void setTolerancesTS(doublereal rtol, doublereal atol) {
for (int n=0;n<m_nv;n++){
m_rtol_ts[n] = rtol;
m_atol_ts[n] = atol;
}
}
void setTolerancesTS(doublereal rtol, doublereal atol);
//added by Karl Meredith
void setTolerancesSS(doublereal rtol, doublereal atol) {
for (int n=0;n<m_nv;n++){
m_rtol_ss[n] = rtol;
m_atol_ss[n] = atol;
}
}
void setTolerancesSS(doublereal rtol, doublereal atol);
/// Relative tolerance of the nth component.
doublereal rtol(int n) { return (m_rdt == 0.0 ? m_rtol_ss[n] : m_rtol_ts[n]); }
@ -271,11 +249,15 @@ namespace Cantera {
* evaluate the residual function at all points.
*/
virtual void eval(int j, doublereal* x, doublereal* r,
integer* mask, doublereal rdt=0.0) {
throw CanteraError("Domain1D::eval",
"residual function not defined.");
integer* mask, doublereal rdt=0.0);
virtual doublereal residual(doublereal* x, int n, int j) {
throw CanteraError("Domain1D::residual","residual function must be overloaded in derived class");
}
int timeDerivativeFlag(int n) { return m_td[n];}
void setAlgebraic(int n) { m_td[n] = 0; }
/**
* Does nothing.
*/
@ -284,7 +266,7 @@ namespace Cantera {
doublereal time() const { return m_time;}
void incrementTime(doublereal dt) { m_time += dt; }
size_t index(int n, int j) const { return m_nv*j + n; }
doublereal value(doublereal* x, int n, int j) const {
doublereal value(const doublereal* x, int n, int j) const {
return x[index(n,j)];
}
@ -396,7 +378,7 @@ namespace Cantera {
virtual void getTransientMask(integer* mask){}
virtual void showSolution(ostream& s, const doublereal* x) {}
virtual void showSolution(const doublereal* x) {}
virtual void showSolution(const doublereal* x);
virtual void restore(const XML_Node& dom, doublereal* soln) {}
@ -425,7 +407,8 @@ namespace Cantera {
const vector_fp& grid() const { return m_z; }
doublereal grid(int point) { return m_z[point]; }
virtual void setupGrid(int n, const doublereal* z) {}
virtual void setupGrid(int n, const doublereal* z);
void setGrid(int n, const doublereal* z);
/**
@ -436,10 +419,12 @@ namespace Cantera {
* prior to installing this domain into the container to be
* written to the global solution vector.
*/
virtual void _getInitialSoln(doublereal* x) {
throw CanteraError("Domain1D::_getInitialSoln",
"base class method _getInitialSoln called!");
}
virtual void _getInitialSoln(doublereal* x);
/**
* Initial value of solution component \a n at grid point \a j.
*/
virtual doublereal initialValue(int n, int j);
/**
* In some cases, a domain may need to set parameters that
@ -479,8 +464,8 @@ namespace Cantera {
Domain1D *m_left, *m_right;
string m_id, m_desc;
Refiner* m_refiner;
vector_int m_td;
vector<string> m_name;
private:

1
Cantera/src/oneD/Inlet1D.h
View file

@ -202,6 +202,7 @@ namespace Cantera {
virtual ~Empty1D(){}
virtual string componentName(int n) const;
virtual void showSolution(const doublereal* x) {}
virtual void init();

11
Cantera/src/oneD/Sim1D.cpp
View file

@ -186,15 +186,18 @@ namespace Cantera {
void Sim1D::showSolution(ostream& s) {
for (int n = 0; n < m_nd; n++) {
domain(n).showSolution(s, m_x.begin() + start(n));
if (domain(n).domainType() != cEmptyType)
domain(n).showSolution(s, m_x.begin() + start(n));
}
}
void Sim1D::showSolution() {
for (int n = 0; n < m_nd; n++) {
writelog("\n\n>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> "+domain(n).id()
+" <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n\n");
domain(n).showSolution(m_x.begin() + start(n));
if (domain(n).domainType() != cEmptyType) {
writelog("\n\n>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> "+domain(n).id()
+" <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n\n");
domain(n).showSolution(m_x.begin() + start(n));
}
}
}

1
Cantera/src/oneD/oneD_files.cpp
View file

@ -6,3 +6,4 @@
#include "boundaries1D.cpp"
#include "refine.cpp"
#include "Sim1D.cpp"
#include "Domain1D.cpp"