Cantera  2.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
vcs_root1d.cpp
Go to the documentation of this file.
1 /**
2  * @file vcs_root1d.cpp
3  * Code for a one dimensional root finder program.
4  */
5 /*
6  * Copyright (2006) Sandia Corporation. Under the terms of
7  * Contract DE-AC04-94AL85000 with Sandia Corporation, the
8  * U.S. Government retains certain rights in this software.
9  */
10 
11 #include "cantera/equil/vcs_internal.h"
12 #include "cantera/equil/vcs_defs.h"
13 #include "cantera/numerics/ctlapack.h"
14 
15 #include <cstdio>
16 
17 namespace Cantera
18 {
19 
20 #define TOL_CONV 1.0E-5
21 
22 static void print_funcEval(FILE* fp, double xval, double fval, int its)
23 {
24  fprintf(fp,"\n");
25  fprintf(fp,"...............................................................\n");
26  fprintf(fp,".................. vcs_root1d Function Evaluation .............\n");
27  fprintf(fp,".................. iteration = %5d ........................\n", its);
28  fprintf(fp,".................. value = %12.5g ......................\n", xval);
29  fprintf(fp,".................. funct = %12.5g ......................\n", fval);
30  fprintf(fp,"...............................................................\n");
31  fprintf(fp,"\n");
32 }
33 
34 int vcsUtil_root1d(double xmin, double xmax, size_t itmax,
35  VCS_FUNC_PTR func, void* fptrPassthrough,
36  double FuncTargVal, int varID,
37  double* xbest, int printLvl)
38 {
39  warn_deprecated("vcsUtil_root1d", "To be removed after Cantera 2.2.");
40  static int callNum = 0;
41  const char* stre = "vcs_root1d ERROR: ";
42  const char* strw = "vcs_root1d WARNING: ";
43  bool converged = false;
44  int err = 0;
45 
46 #ifdef DEBUG_MODE
47  char fileName[80];
48 #else
49  char* fileName;
50 #endif
51  FILE* fp = 0;
52 
53  double x1, x2, xnew, f1, f2, fnew, slope;
54  size_t its = 0;
55  int posStraddle = 0;
56  int retn = VCS_SUCCESS;
57  bool foundPosF = false;
58  bool foundNegF = false;
59  bool foundStraddle = false;
60  double xPosF = 0.0;
61  double xNegF = 0.0;
62  double fnorm; /* A valid norm for the making the function value
63  * dimensionless */
64  double c[9], f[3], xn1, xn2, x0 = 0.0, f0 = 0.0, root, theta, xquad;
65 
66  callNum++;
67  if (DEBUG_MODE_ENABLED && printLvl >= 3) {
68  sprintf(fileName, "rootfd_%d.log", callNum);
69  fp = fopen(fileName, "w");
70  fprintf(fp, " Iter TP_its xval Func_val | Reasoning\n");
71  fprintf(fp, "-----------------------------------------------------"
72  "-------------------------------\n");
73  } else if (printLvl >= 3) {
74  plogf("WARNING: vcsUtil_root1d: printlvl >= 3, but debug mode not turned on\n");
75  }
76  if (xmax <= xmin) {
77  plogf("%sxmin and xmax are bad: %g %g\n", stre, xmin, xmax);
78  return VCS_PUB_BAD;
79  }
80  x1 = *xbest;
81  if (x1 < xmin || x1 > xmax) {
82  x1 = (xmin + xmax) / 2.0;
83  }
84  f1 = func(x1, FuncTargVal, varID, fptrPassthrough, &err);
85  if (DEBUG_MODE_ENABLED && printLvl >= 3) {
86  print_funcEval(fp, x1, f1, static_cast<int>(its));
87  fprintf(fp, "%-5d %-5d %-15.5E %-15.5E\n", -2, 0, x1, f1);
88  }
89  if (f1 == 0.0) {
90  *xbest = x1;
91  return VCS_SUCCESS;
92  } else if (f1 > 0.0) {
93  foundPosF = true;
94  xPosF = x1;
95  } else {
96  foundNegF = true;
97  xNegF = x1;
98  }
99  x2 = x1 * 1.1;
100  if (x2 > xmax) {
101  x2 = x1 - (xmax - xmin) / 100.;
102  }
103  f2 = func(x2, FuncTargVal, varID, fptrPassthrough, &err);
104  if (DEBUG_MODE_ENABLED && printLvl >= 3) {
105  print_funcEval(fp, x2, f2, static_cast<int>(its));
106  fprintf(fp, "%-5d %-5d %-15.5E %-15.5E", -1, 0, x2, f2);
107  }
108 
109  if (FuncTargVal != 0.0) {
110  fnorm = fabs(FuncTargVal) + 1.0E-13;
111  } else {
112  fnorm = 0.5*(fabs(f1) + fabs(f2)) + fabs(FuncTargVal);
113  }
114 
115  if (f2 == 0.0) {
116  return retn;
117  } else if (f2 > 0.0) {
118  if (!foundPosF) {
119  foundPosF = true;
120  xPosF = x2;
121  }
122  } else {
123  if (!foundNegF) {
124  foundNegF = true;
125  xNegF = x2;
126  }
127  }
128  foundStraddle = foundPosF && foundNegF;
129  if (foundStraddle) {
130  if (xPosF > xNegF) {
131  posStraddle = true;
132  } else {
133  posStraddle = false;
134  }
135  }
136  int ipiv[3];
137  int info;
138 
139  do {
140  /*
141  * Find an estimate of the next point to try based on
142  * a linear approximation.
143  */
144  slope = (f2 - f1) / (x2 - x1);
145  if (slope == 0.0) {
146  plogf("%s functions evals produced the same result, %g, at %g and %g\n",
147  strw, f2, x1, x2);
148  xnew = 2*x2 - x1 + 1.0E-3;
149  } else {
150  xnew = x2 - f2 / slope;
151  }
152  if (DEBUG_MODE_ENABLED && printLvl >= 3) {
153  fprintf(fp, " | xlin = %-9.4g", xnew);
154  }
155 
156  /*
157  * Do a quadratic fit -> Note this algorithm seems
158  * to work OK. The quadratic approximation doesn't kick in until
159  * the end of the run, when it becomes reliable.
160  */
161  if (its > 0) {
162  c[0] = 1.;
163  c[1] = 1.;
164  c[2] = 1.;
165  c[3] = x0;
166  c[4] = x1;
167  c[5] = x2;
168  c[6] = pow(x0, 2);
169  c[7] = pow(x1, 2);
170  c[8] = pow(x2, 2);
171  f[0] = f0;
172  f[1] = f1;
173  f[2] = f2;
174 
175  ct_dgetrf(3, 3, c, 3, ipiv, info);
176  if (info) {
177  goto QUAD_BAIL;
178  }
179  ct_dgetrs(ctlapack::NoTranspose, 3, 1, c, 3, ipiv, f, 3, info);
180  root = f[1]* f[1] - 4.0 * f[0] * f[2];
181  if (root >= 0.0) {
182  xn1 = (- f[1] + sqrt(root)) / (2.0 * f[2]);
183  xn2 = (- f[1] - sqrt(root)) / (2.0 * f[2]);
184  if (fabs(xn2 - x2) < fabs(xn1 - x2) && xn2 > 0.0) {
185  xquad = xn2;
186  } else {
187  xquad = xn1;
188  }
189  theta = fabs(xquad - xnew) / fabs(xnew - x2);
190  theta = std::min(1.0, theta);
191  xnew = theta * xnew + (1.0 - theta) * xquad;
192  if (DEBUG_MODE_ENABLED && printLvl >= 3) {
193  if (theta != 1.0) {
194  fprintf(fp, " | xquad = %-9.4g", xnew);
195  }
196  }
197  } else {
198  /*
199  * Pick out situations where the convergence may be
200  * accelerated.
201  */
202  if ((sign(xnew - x2) == sign(x2 - x1)) &&
203  (sign(x2 - x1) == sign(x1 - x0))) {
204  xnew += xnew - x2;
205  if (DEBUG_MODE_ENABLED && printLvl >= 3) {
206  fprintf(fp, " | xquada = %-9.4g", xnew);
207  }
208  }
209  }
210  }
211 QUAD_BAIL:
212  ;
213 
214 
215  /*
216  *
217  * Put heuristic bounds on the step jump
218  */
219  if ((xnew > x1 && xnew < x2) || (xnew < x1 && xnew > x2)) {
220  /*
221  *
222  * If we are doing a jump in between two points, make sure
223  * the new trial is between 10% and 90% of the distance
224  * between the old points.
225  */
226  slope = fabs(x2 - x1) / 10.;
227  if (fabs(xnew - x1) < slope) {
228  xnew = x1 + sign(xnew-x1) * slope;
229  if (DEBUG_MODE_ENABLED && printLvl >= 3) {
230  fprintf(fp, " | x10%% = %-9.4g", xnew);
231  }
232  }
233  if (fabs(xnew - x2) < slope) {
234  xnew = x2 + sign(xnew-x2) * slope;
235  if (DEBUG_MODE_ENABLED && printLvl >= 3) {
236  fprintf(fp, " | x10%% = %-9.4g", xnew);
237  }
238  }
239  } else {
240  /*
241  * If we are venturing into new ground, only allow the step jump
242  * to increase by 100% at each iteration
243  */
244  slope = 2.0 * fabs(x2 - x1);
245  if (fabs(slope) < fabs(xnew - x2)) {
246  xnew = x2 + sign(xnew-x2) * slope;
247  if (DEBUG_MODE_ENABLED && printLvl >= 3) {
248  fprintf(fp, " | xlimitsize = %-9.4g", xnew);
249  }
250  }
251  }
252 
253 
254  if (xnew > xmax) {
255  xnew = x2 + (xmax - x2) / 2.0;
256  if (DEBUG_MODE_ENABLED && printLvl >= 3) {
257  fprintf(fp, " | xlimitmax = %-9.4g", xnew);
258  }
259  }
260  if (xnew < xmin) {
261  xnew = x2 + (x2 - xmin) / 2.0;
262  if (DEBUG_MODE_ENABLED && printLvl >= 3) {
263  fprintf(fp, " | xlimitmin = %-9.4g", xnew);
264  }
265  }
266  if (foundStraddle) {
267 #ifdef DEBUG_MODE
268  slope = xnew;
269 #endif
270  if (posStraddle) {
271  if (f2 > 0.0) {
272  if (xnew > x2) {
273  xnew = (xNegF + x2)/2;
274  }
275  if (xnew < xNegF) {
276  xnew = (xNegF + x2)/2;
277  }
278  } else {
279  if (xnew < x2) {
280  xnew = (xPosF + x2)/2;
281  }
282  if (xnew > xPosF) {
283  xnew = (xPosF + x2)/2;
284  }
285  }
286  } else {
287  if (f2 > 0.0) {
288  if (xnew < x2) {
289  xnew = (xNegF + x2)/2;
290  }
291  if (xnew > xNegF) {
292  xnew = (xNegF + x2)/2;
293  }
294  } else {
295  if (xnew > x2) {
296  xnew = (xPosF + x2)/2;
297  }
298  if (xnew < xPosF) {
299  xnew = (xPosF + x2)/2;
300  }
301  }
302  }
303  if (DEBUG_MODE_ENABLED && printLvl >= 3) {
304  if (slope != xnew) {
305  fprintf(fp, " | xstraddle = %-9.4g", xnew);
306  }
307  }
308  }
309 
310  fnew = func(xnew, FuncTargVal, varID, fptrPassthrough, &err);
311  if (DEBUG_MODE_ENABLED && printLvl >= 3) {
312  fprintf(fp,"\n");
313  print_funcEval(fp, xnew, fnew, static_cast<int>(its));
314  fprintf(fp, "%-5d %-5d %-15.5E %-15.5E", (int) its, 0, xnew, fnew);
315  }
316 
317  if (foundStraddle) {
318  if (posStraddle) {
319  if (fnew > 0.0) {
320  if (xnew < xPosF) {
321  xPosF = xnew;
322  }
323  } else {
324  if (xnew > xNegF) {
325  xNegF = xnew;
326  }
327  }
328  } else {
329  if (fnew > 0.0) {
330  if (xnew > xPosF) {
331  xPosF = xnew;
332  }
333  } else {
334  if (xnew < xNegF) {
335  xNegF = xnew;
336  }
337  }
338  }
339  }
340 
341  if (! foundStraddle) {
342  if (fnew > 0.0) {
343  if (!foundPosF) {
344  foundPosF = true;
345  xPosF = xnew;
346  foundStraddle = true;
347  posStraddle = (xPosF > xNegF);
348  }
349  } else {
350  if (!foundNegF) {
351  foundNegF = true;
352  xNegF = xnew;
353  foundStraddle = true;
354  posStraddle = (xPosF > xNegF);
355  }
356  }
357  }
358 
359  x0 = x1;
360  f0 = f1;
361  x1 = x2;
362  f1 = f2;
363  x2 = xnew;
364  f2 = fnew;
365  if (fabs(fnew / fnorm) < 1.0E-5) {
366  converged = true;
367  }
368  its++;
369  } while (! converged && its < itmax);
370  if (converged) {
371  if (printLvl >= 1) {
372  plogf("vcs_root1d success: convergence achieved\n");
373  }
374  if (DEBUG_MODE_ENABLED && printLvl >= 3) {
375  fprintf(fp, " | vcs_root1d success in %d its, fnorm = %g\n", (int) its, fnorm);
376  }
377  } else {
378  retn = VCS_FAILED_CONVERGENCE;
379  if (printLvl >= 1) {
380  plogf("vcs_root1d ERROR: maximum iterations exceeded without convergence\n");
381  }
382  if (DEBUG_MODE_ENABLED && printLvl >= 3) {
383  fprintf(fp, "\nvcs_root1d failure in %lu its\n", its);
384  }
385  }
386  *xbest = x2;
387  if (DEBUG_MODE_ENABLED && printLvl >= 3) {
388  fclose(fp);
389  }
390  return retn;
391 }
392 
393 }
Cantera::warn_deprecated
void warn_deprecated(const std::string &method, const std::string &extra)
Print a warning indicating that method is deprecated.
Definition: global.cpp:78
Cantera::sign
int sign(T x)
Sign of a number. Returns -1 if x < 0, 1 if x > 0 and 0 if x == 0.
Definition: global.h:276
Cantera::vcsUtil_root1d
int vcsUtil_root1d(double xmin, double xmax, size_t itmax, VCS_FUNC_PTR func, void *fptrPassthrough, double FuncTargVal, int varID, double *xbest, int printLvl)
One dimensional root finder.
Definition: vcs_root1d.cpp:34
vcs_defs.h
Defines and definitions within the vcs package.
vcs_internal.h
Internal declarations for the VCSnonideal package.
VCS_SUCCESS
#define VCS_SUCCESS
Definition: vcs_defs.h:21
Cantera::VCS_FUNC_PTR
double(* VCS_FUNC_PTR)(double xval, double Vtarget, int varID, void *fptrPassthrough, int *err)
Definition of the function pointer for the root finder.
Definition: vcs_internal.h:107
plogf
#define plogf
define this Cantera function to replace printf
Definition: vcs_internal.h:24
Generated by   doxygen 1.8.6