6 #include "cantera/oneD/refine.h" 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),
18 m_nv = m_domain->nComponents();
19 m_active.resize(m_nv,
true);
20 m_thresh = std::sqrt(std::numeric_limits<double>::epsilon());
24 doublereal curve, doublereal prune)
28 "'ratio' must be greater than 2.0 ({} was specified).", ratio);
29 }
else if (slope < 0.0 || slope > 1.0) {
31 "'slope' must be between 0.0 and 1.0 ({} was specified).", slope);
32 }
else if (curve < 0.0 || curve > 1.0) {
34 "'curve' must be between 0.0 and 1.0 ({} was specified).", curve);
35 }
else if (prune > curve || prune > slope) {
37 "'prune' must be less than 'curve' and 'slope' ({} was specified).",
46 int Refiner::analyze(
size_t n,
const doublereal* z,
50 throw CanteraError(
"Refiner::analyze",
"max number of grid points reached ({}).", m_npmax);
68 throw CanteraError(
"Refiner::analyze",
"inconsistent");
75 for (
size_t j = 0; j < n-1; j++) {
76 dz[j] = z[j+1] - z[j];
79 for (
size_t i = 0; i < m_nv; i++) {
83 for (
size_t j = 0; j < n; j++) {
84 v[j] = value(x, i, j);
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]);
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());
99 doublereal aa = std::max(fabs(vmax), fabs(vmin));
100 doublereal ss = std::max(fabs(smax), fabs(smin));
106 if ((vmax - vmin) > m_min_range*aa) {
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;
119 }
else if (m_keep[j] == 0) {
129 if ((smax - smin) > m_min_range*ss) {
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]);
143 }
else if (m_keep[j+1] == 0) {
151 StFlow* fflame =
dynamic_cast<StFlow*
>(m_domain);
154 for (
size_t j = 1; j < n-1; j++) {
156 if (dz[j] > m_ratio*dz[j-1]) {
158 m_c[fmt::format(
"point {}", j)] = 1;
166 if (dz[j] < dz[j-1]/m_ratio) {
168 m_c[fmt::format(
"point {}", j-1)] = 1;
177 if (j > 1 && z[j+1]-z[j-1] > m_ratio * dz[j-2]) {
183 if (j < n-2 && z[j+1]-z[j-1] > m_ratio * dz[j+1]) {
188 if (fflame->domainType() == cFreeFlow && z[j] == fflame->m_zfixed) {
195 for (
size_t j = 2; j < n-1; j++) {
196 if (m_keep[j] == -1 && m_keep[j-1] == -1) {
201 return int(m_loc.size());
204 double Refiner::value(
const double* x,
size_t i,
size_t j)
206 return x[m_domain->index(i,j)];
211 if (!m_loc.empty()) {
213 writelog(
string(
"Refining grid in ") +
215 +
" New points inserted after grid points ");
216 for (
const auto& loc : m_loc) {
221 for (
const auto& c : m_c) {
226 }
else if (m_domain->
nPoints() > 1) {
227 writelog(
"no new points needed in "+m_domain->id()+
"\n");
231 int Refiner::getNewGrid(
int n,
const doublereal* z,
232 int nn, doublereal* zn)
234 int nnew =
static_cast<int>(m_loc.size());
236 throw CanteraError(
"Refine::getNewGrid",
"array size too small.");
245 for (
int j = 0; j < n - 1; j++) {
248 if (m_loc.count(j)) {
249 zn[jn] = 0.5*(z[j] + z[j+1]);
void writelog(const std::string &fmt, const Args &... args)
Write a formatted message to the screen.
doublereal m_gridmin
minimum grid spacing [m]
void setCriteria(doublereal ratio=10.0, doublereal slope=0.8, doublereal curve=0.8, doublereal prune=-0.1)
Set grid refinement criteria.
Base class for exceptions thrown by Cantera classes.
virtual std::string componentName(size_t n) const
Name of the nth component. May be overloaded.
size_t nComponents() const
Number of components at each grid point.
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
size_t nPoints() const
Number of grid points in this domain.
Namespace for the Cantera kernel.