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

[1D] General cleanup of class Refiner

Browse source
This commit is contained in:
Ray Speth 2013-02-07 23:41:49 +00:00
parent dc1ffee1a2
commit e48bd48c4e
2 changed files with 37 additions and 80 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

19
include/cantera/oneD/refine.h
View file

@ -10,9 +10,7 @@ class Domain1D;
class Refiner class Refiner
{ {
public: public:
Refiner(Domain1D& domain); Refiner(Domain1D& domain);
virtual ~Refiner() {} virtual ~Refiner() {}
@ -45,7 +43,6 @@ public:
int analyze(size_t n, const doublereal* z, const doublereal* x); int analyze(size_t n, const doublereal* z, const doublereal* x);
int getNewGrid(int n, const doublereal* z, int nn, doublereal* znew); int getNewGrid(int n, const doublereal* z, int nn, doublereal* znew);
//int getNewSoln(int n, const doublereal* x, doublereal* xnew);
int nNewPoints() { int nNewPoints() {
return static_cast<int>(m_loc.size()); return static_cast<int>(m_loc.size());
} }
@ -54,7 +51,7 @@ public:
return m_loc.find(j) != m_loc.end(); return m_loc.find(j) != m_loc.end();
} }
bool keepPoint(size_t j) { bool keepPoint(size_t j) {
return (m_keep[j] != -1); // m_keep.find(j) != m_keep.end(); return (m_keep[j] != -1);
} }
double value(const double* x, size_t i, size_t j); double value(const double* x, size_t i, size_t j);
double maxRatio() { double maxRatio() {
@ -71,18 +68,16 @@ public:
} }
protected: protected:
std::map<size_t, int> m_loc; std::map<size_t, int> m_loc;
std::map<size_t, int> m_keep; std::map<size_t, int> m_keep;
std::map<std::string, int> m_c; std::map<std::string, int> m_c;
std::vector<bool> m_active; std::vector<bool> m_active;
doublereal m_ratio, m_slope, m_curve, m_prune; doublereal m_ratio, m_slope, m_curve, m_prune;
doublereal m_min_range; doublereal m_min_range;
Domain1D* m_domain; Domain1D* m_domain;
size_t m_nv, m_npmax; size_t m_nv, m_npmax;
doublereal m_thresh; doublereal m_thresh;
doublereal m_gridmin; //!< minimum grid spacing [m] doublereal m_gridmin; //!< minimum grid spacing [m]
}; };
} }

98
src/oneD/refine.cpp
View file

@ -37,7 +37,6 @@ int Refiner::analyze(size_t n, const doublereal* z,
} }
if (m_domain->nPoints() <= 1) { if (m_domain->nPoints() <= 1) {
//writelog("can't refine a domain with 1 point: "+m_domain->id()+"\n");
return 0; return 0;
} }
@ -55,114 +54,87 @@ int Refiner::analyze(size_t n, const doublereal* z,
throw CanteraError("analyze","inconsistent"); throw CanteraError("analyze","inconsistent");
} }
// find locations where cell size ratio is too large.
/**
* find locations where cell size ratio is too large.
*/
size_t j;
string name;
doublereal vmin, vmax, smin, smax, aa, ss;
doublereal dmax, r;
vector_fp v(n), s(n-1); vector_fp v(n), s(n-1);
vector_fp dz(n-1); vector_fp dz(n-1);
for (j = 0; j < n-1; j++) { for (size_t j = 0; j < n-1; j++) {
dz[j] = z[j+1] - z[j]; dz[j] = z[j+1] - z[j];
} }
for (size_t i = 0; i < m_nv; i++) { for (size_t i = 0; i < m_nv; i++) {
if (m_active[i]) { if (m_active[i]) {
name = m_domain->componentName(i); string name = m_domain->componentName(i);
//writelog("refine: examining "+name+"\n");
// get component i at all points // get component i at all points
for (j = 0; j < n; j++) { for (size_t j = 0; j < n; j++) {
v[j] = value(x, i, j); v[j] = value(x, i, j);
} }
// slope of component i // slope of component i
for (j = 0; j < n-1; j++) for (size_t j = 0; j < n-1; j++) {
s[j] = (value(x, i, j+1) - value(x, i, j))/ s[j] = (value(x, i, j+1) - value(x, i, j))/(z[j+1] - z[j]);
(z[j+1] - z[j]); }
// find the range of values and slopes // find the range of values and slopes
doublereal vmin = *min_element(v.begin(), v.end());
vmin = *min_element(v.begin(), v.end()); doublereal vmax = *max_element(v.begin(), v.end());
vmax = *max_element(v.begin(), v.end()); doublereal smin = *min_element(s.begin(), s.end());
smin = *min_element(s.begin(), s.end()); doublereal smax = *max_element(s.begin(), s.end());
smax = *max_element(s.begin(), s.end());
// max absolute values of v and s // max absolute values of v and s
aa = std::max(fabs(vmax), fabs(vmin)); doublereal aa = std::max(fabs(vmax), fabs(vmin));
ss = std::max(fabs(smax), fabs(smin)); doublereal ss = std::max(fabs(smax), fabs(smin));
// refine based on component i only if the range of v is // refine based on component i only if the range of v is
// greater than a fraction 'min_range' of max |v|. This // greater than a fraction 'min_range' of max |v|. This
// eliminates components that consist of small fluctuations // eliminates components that consist of small fluctuations
// on a constant background. // on a constant background.
if ((vmax - vmin) > m_min_range*aa) { if ((vmax - vmin) > m_min_range*aa) {
// maximum allowable difference in value between adjacent
// maximum allowable difference in value between // points.
// adjacent points. doublereal dmax = m_slope*(vmax - vmin) + m_thresh;
for (size_t j = 0; j < n-1; j++) {
dmax = m_slope*(vmax - vmin) + m_thresh; doublereal r = fabs(v[j+1] - v[j])/dmax;
for (j = 0; j < n-1; j++) {
r = fabs(v[j+1] - v[j])/dmax;
if (r > 1.0 && dz[j] >= 2 * m_gridmin) { if (r > 1.0 && dz[j] >= 2 * m_gridmin) {
m_loc[j] = 1; m_loc[j] = 1;
m_c[name] = 1; m_c[name] = 1;
//if (int(m_loc.size()) + n > m_npmax) goto done;
} }
if (r >= m_prune) { if (r >= m_prune) {
m_keep[j] = 1; m_keep[j] = 1;
m_keep[j+1] = 1; m_keep[j+1] = 1;
} else { } else if (m_keep[j] == 0) {
//writelog(string("r = ")+fp2str(r)+"\n"); m_keep[j] = -1;
if (m_keep[j] == 0) {
//if (m_keep[j-1] > -1 && m_keep[j+1] > -1)
m_keep[j] = -1;
}
//if (m_keep[j+1] == 0) m_keep[j+1] = -1;
} }
} }
} }
// refine based on the slope of component i only if the // refine based on the slope of component i only if the
// range of s is greater than a fraction 'min_range' of max // range of s is greater than a fraction 'min_range' of max
// |s|. This eliminates components that consist of small // |s|. This eliminates components that consist of small
// fluctuations on a constant slope background. // fluctuations on a constant slope background.
if ((smax - smin) > m_min_range*ss) { if ((smax - smin) > m_min_range*ss) {
// maximum allowable difference in slope between // maximum allowable difference in slope between
// adjacent points. // adjacent points.
dmax = m_curve*(smax - smin); // + 0.5*m_curve*(smax + smin); doublereal dmax = m_curve*(smax - smin);
for (j = 0; j < n-2; j++) { for (size_t j = 0; j < n-2; j++) {
r = fabs(s[j+1] - s[j]) / (dmax + m_thresh/dz[j]); doublereal r = fabs(s[j+1] - s[j]) / (dmax + m_thresh/dz[j]);
if (r > 1.0 && dz[j] >= 2 * m_gridmin && if (r > 1.0 && dz[j] >= 2 * m_gridmin &&
dz[j+1] >= 2 * m_gridmin) { dz[j+1] >= 2 * m_gridmin) {
m_c[name] = 1; m_c[name] = 1;
m_loc[j] = 1; m_loc[j] = 1;
m_loc[j+1] = 1; m_loc[j+1] = 1;
//if (int(m_loc.size()) + n > m_npmax) goto done;
} }
if (r >= m_prune) { if (r >= m_prune) {
m_keep[j+1] = 1; m_keep[j+1] = 1;
} else { } else if (m_keep[j+1] == 0) {
//writelog(string("r slope = ")+fp2str(r)+"\n"); m_keep[j+1] = -1;
if (m_keep[j+1] == 0) {
//if (m_keep[j] > -1 && m_keep[j+2] > -1)
m_keep[j+1] = -1;
}
} }
} }
} }
} }
} }
for (j = 1; j < n-1; j++) { for (size_t j = 1; j < n-1; j++) {
if (dz[j] > m_ratio*dz[j-1]) { if (dz[j] > m_ratio*dz[j-1]) {
m_loc[j] = 1; m_loc[j] = 1;
m_c["point "+int2str(j)] = 1; m_c["point "+int2str(j)] = 1;
@ -177,11 +149,8 @@ int Refiner::analyze(size_t n, const doublereal* z,
if (j < n-2 && z[j+1]-z[j] > m_ratio * dz[j+1]) { if (j < n-2 && z[j+1]-z[j] > m_ratio * dz[j+1]) {
m_keep[j] = 1; m_keep[j] = 1;
} }
//if (m_loc.size() + n > m_npmax) goto done;
} }
//done:
//m_did_analysis = true;
return int(m_loc.size()); return int(m_loc.size());
} }
@ -192,8 +161,7 @@ double Refiner::value(const double* x, size_t i, size_t j)
void Refiner::show() void Refiner::show()
{ {
int nnew = static_cast<int>(m_loc.size()); if (!m_loc.empty()) {
if (nnew > 0) {
r_drawline(); r_drawline();
writelog(string("Refining grid in ") + writelog(string("Refining grid in ") +
m_domain->id()+".\n" m_domain->id()+".\n"
@ -212,9 +180,6 @@ void Refiner::show()
r_drawline(); r_drawline();
} else if (m_domain->nPoints() > 1) { } else if (m_domain->nPoints() > 1) {
writelog("no new points needed in "+m_domain->id()+"\n"); writelog("no new points needed in "+m_domain->id()+"\n");
//writelog("curve = "+fp2str(m_curve)+"\n");
//writelog("slope = "+fp2str(m_slope)+"\n");
//writelog("prune = "+fp2str(m_prune)+"\n");
} }
} }
@ -222,21 +187,18 @@ void Refiner::show()
int Refiner::getNewGrid(int n, const doublereal* z, int Refiner::getNewGrid(int n, const doublereal* z,
int nn, doublereal* zn) int nn, doublereal* zn)
{ {
int j;
int nnew = static_cast<int>(m_loc.size()); int nnew = static_cast<int>(m_loc.size());
if (nnew + n > nn) { if (nnew + n > nn) {
throw CanteraError("Refine::getNewGrid", throw CanteraError("Refine::getNewGrid", "array size too small.");
"array size too small.");
return -1;
} }
int jn = 0;
if (m_loc.empty()) { if (m_loc.empty()) {
copy(z, z + n, zn); copy(z, z + n, zn);
return 0; return 0;
} }
for (j = 0; j < n - 1; j++) { int jn = 0;
for (int j = 0; j < n - 1; j++) {
zn[jn] = z[j]; zn[jn] = z[j];
jn++; jn++;
if (m_loc.count(j)) { if (m_loc.count(j)) {