Cantera  2.5.1
refine.cpp
Go to the documentation of this file.
1 //! @file refine.cpp
2 
3 // This file is part of Cantera. See License.txt in the top-level directory or
4 // at https://cantera.org/license.txt for license and copyright information.
5 
6 #include "cantera/oneD/refine.h"
7 #include "cantera/oneD/StFlow.h"
8 
9 using namespace std;
10 
11 namespace Cantera
12 {
13 Refiner::Refiner(Domain1D& domain) :
14  m_ratio(10.0), m_slope(0.8), m_curve(0.8), m_prune(-0.001),
15  m_min_range(0.01), m_domain(&domain), m_npmax(1000),
16  m_gridmin(1e-10)
17 {
18  m_nv = m_domain->nComponents();
19  m_active.resize(m_nv, true);
20  m_thresh = std::sqrt(std::numeric_limits<double>::epsilon());
21 }
22 
23 void Refiner::setCriteria(doublereal ratio, doublereal slope,
24  doublereal curve, doublereal prune)
25 {
26  if (ratio < 2.0) {
27  throw CanteraError("Refiner::setCriteria",
28  "'ratio' must be greater than 2.0 ({} was specified).", ratio);
29  } else if (slope < 0.0 || slope > 1.0) {
30  throw CanteraError("Refiner::setCriteria",
31  "'slope' must be between 0.0 and 1.0 ({} was specified).", slope);
32  } else if (curve < 0.0 || curve > 1.0) {
33  throw CanteraError("Refiner::setCriteria",
34  "'curve' must be between 0.0 and 1.0 ({} was specified).", curve);
35  } else if (prune > curve || prune > slope) {
36  throw CanteraError("Refiner::setCriteria",
37  "'prune' must be less than 'curve' and 'slope' ({} was specified).",
38  prune);
39  }
40  m_ratio = ratio;
41  m_slope = slope;
42  m_curve = curve;
43  m_prune = prune;
44 }
45 
46 int Refiner::analyze(size_t n, const doublereal* z,
47  const doublereal* x)
48 {
49  if (n >= m_npmax) {
50  throw CanteraError("Refiner::analyze", "max number of grid points reached ({}).", m_npmax);
51  }
52 
53  if (m_domain->nPoints() <= 1) {
54  return 0;
55  }
56 
57  m_loc.clear();
58  m_c.clear();
59  m_keep.clear();
60 
61  m_keep[0] = 1;
62  m_keep[n-1] = 1;
63 
64  m_nv = m_domain->nComponents();
65 
66  // check consistency
67  if (n != m_domain->nPoints()) {
68  throw CanteraError("Refiner::analyze", "inconsistent");
69  }
70 
71  // find locations where cell size ratio is too large.
72  vector_fp v(n), s(n-1);
73 
74  vector_fp dz(n-1);
75  for (size_t j = 0; j < n-1; j++) {
76  dz[j] = z[j+1] - z[j];
77  }
78 
79  for (size_t i = 0; i < m_nv; i++) {
80  if (m_active[i]) {
81  string name = m_domain->componentName(i);
82  // get component i at all points
83  for (size_t j = 0; j < n; j++) {
84  v[j] = value(x, i, j);
85  }
86 
87  // slope of component i
88  for (size_t j = 0; j < n-1; j++) {
89  s[j] = (value(x, i, j+1) - value(x, i, j))/(z[j+1] - z[j]);
90  }
91 
92  // find the range of values and slopes
93  doublereal vmin = *min_element(v.begin(), v.end());
94  doublereal vmax = *max_element(v.begin(), v.end());
95  doublereal smin = *min_element(s.begin(), s.end());
96  doublereal smax = *max_element(s.begin(), s.end());
97 
98  // max absolute values of v and s
99  doublereal aa = std::max(fabs(vmax), fabs(vmin));
100  doublereal ss = std::max(fabs(smax), fabs(smin));
101 
102  // refine based on component i only if the range of v is
103  // greater than a fraction 'min_range' of max |v|. This
104  // eliminates components that consist of small fluctuations
105  // on a constant background.
106  if ((vmax - vmin) > m_min_range*aa) {
107  // maximum allowable difference in value between adjacent
108  // points.
109  doublereal dmax = m_slope*(vmax - vmin) + m_thresh;
110  for (size_t j = 0; j < n-1; j++) {
111  doublereal r = fabs(v[j+1] - v[j])/dmax;
112  if (r > 1.0 && dz[j] >= 2 * m_gridmin) {
113  m_loc[j] = 1;
114  m_c[name] = 1;
115  }
116  if (r >= m_prune) {
117  m_keep[j] = 1;
118  m_keep[j+1] = 1;
119  } else if (m_keep[j] == 0) {
120  m_keep[j] = -1;
121  }
122  }
123  }
124 
125  // refine based on the slope of component i only if the
126  // range of s is greater than a fraction 'min_range' of max
127  // |s|. This eliminates components that consist of small
128  // fluctuations on a constant slope background.
129  if ((smax - smin) > m_min_range*ss) {
130  // maximum allowable difference in slope between
131  // adjacent points.
132  doublereal dmax = m_curve*(smax - smin);
133  for (size_t j = 0; j < n-2; j++) {
134  doublereal r = fabs(s[j+1] - s[j]) / (dmax + m_thresh/dz[j]);
135  if (r > 1.0 && dz[j] >= 2 * m_gridmin &&
136  dz[j+1] >= 2 * m_gridmin) {
137  m_c[name] = 1;
138  m_loc[j] = 1;
139  m_loc[j+1] = 1;
140  }
141  if (r >= m_prune) {
142  m_keep[j+1] = 1;
143  } else if (m_keep[j+1] == 0) {
144  m_keep[j+1] = -1;
145  }
146  }
147  }
148  }
149  }
150 
151  StFlow* fflame = dynamic_cast<StFlow*>(m_domain);
152 
153  // Refine based on properties of the grid itself
154  for (size_t j = 1; j < n-1; j++) {
155  // Add a new point if the ratio with left interval is too large
156  if (dz[j] > m_ratio*dz[j-1]) {
157  m_loc[j] = 1;
158  m_c[fmt::format("point {}", j)] = 1;
159  m_keep[j-1] = 1;
160  m_keep[j] = 1;
161  m_keep[j+1] = 1;
162  m_keep[j+2] = 1;
163  }
164 
165  // Add a point if the ratio with right interval is too large
166  if (dz[j] < dz[j-1]/m_ratio) {
167  m_loc[j-1] = 1;
168  m_c[fmt::format("point {}", j-1)] = 1;
169  m_keep[j-2] = 1;
170  m_keep[j-1] = 1;
171  m_keep[j] = 1;
172  m_keep[j+1] = 1;
173  }
174 
175  // Keep the point if removing would make the ratio with the left
176  // interval too large.
177  if (j > 1 && z[j+1]-z[j-1] > m_ratio * dz[j-2]) {
178  m_keep[j] = 1;
179  }
180 
181  // Keep the point if removing would make the ratio with the right
182  // interval too large.
183  if (j < n-2 && z[j+1]-z[j-1] > m_ratio * dz[j+1]) {
184  m_keep[j] = 1;
185  }
186 
187  // Keep the point where the temperature is fixed
188  if (fflame && fflame->domainType() == cFreeFlow && z[j] == fflame->m_zfixed) {
189  m_keep[j] = 1;
190  }
191  }
192 
193  // Don't allow pruning to remove multiple adjacent grid points
194  // in a single pass.
195  for (size_t j = 2; j < n-1; j++) {
196  if (m_keep[j] == -1 && m_keep[j-1] == -1) {
197  m_keep[j] = 1;
198  }
199  }
200 
201  return int(m_loc.size());
202 }
203 
204 double Refiner::value(const double* x, size_t i, size_t j)
205 {
206  return x[m_domain->index(i,j)];
207 }
208 
209 void Refiner::show()
210 {
211  if (!m_loc.empty()) {
212  writeline('#', 78);
213  writelog(string("Refining grid in ") +
214  m_domain->id()+".\n"
215  +" New points inserted after grid points ");
216  for (const auto& loc : m_loc) {
217  writelog("{} ", loc.first);
218  }
219  writelog("\n");
220  writelog(" to resolve ");
221  for (const auto& c : m_c) {
222  writelog(string(c.first)+" ");
223  }
224  writelog("\n");
225  writeline('#', 78);
226  } else if (m_domain->nPoints() > 1) {
227  writelog("no new points needed in "+m_domain->id()+"\n");
228  }
229 }
230 
231 int Refiner::getNewGrid(int n, const doublereal* z,
232  int nn, doublereal* zn)
233 {
234  int nnew = static_cast<int>(m_loc.size());
235  if (nnew + n > nn) {
236  throw CanteraError("Refine::getNewGrid", "array size too small.");
237  }
238 
239  if (m_loc.empty()) {
240  copy(z, z + n, zn);
241  return 0;
242  }
243 
244  int jn = 0;
245  for (int j = 0; j < n - 1; j++) {
246  zn[jn] = z[j];
247  jn++;
248  if (m_loc.count(j)) {
249  zn[jn] = 0.5*(z[j] + z[j+1]);
250  jn++;
251  }
252  }
253  zn[jn] = z[n-1];
254  return 0;
255 }
256 }
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:61
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
Definition: ct_defs.h:180
void writelog(const std::string &fmt, const Args &... args)
Write a formatted message to the screen.
Definition: global.h:158
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:264