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