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