Cantera  2.0
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 
12 #include "cantera/equil/vcs_internal.h"
13 
14 #include <cstdio>
15 #include <cstdlib>
16 #include <cmath>
17 
18 namespace VCSnonideal
19 {
20 
21 #define TOL_CONV 1.0E-5
22 /*****************************************************************************/
23 /*****************************************************************************/
24 /*****************************************************************************/
25 #ifdef DEBUG_MODE
26 static void print_funcEval(FILE* fp, double xval, double fval, int its)
27 {
28  fprintf(fp,"\n");
29  fprintf(fp,"...............................................................\n");
30  fprintf(fp,".................. vcs_root1d Function Evaluation .............\n");
31  fprintf(fp,".................. iteration = %5d ........................\n", its);
32  fprintf(fp,".................. value = %12.5g ......................\n", xval);
33  fprintf(fp,".................. funct = %12.5g ......................\n", fval);
34  fprintf(fp,"...............................................................\n");
35  fprintf(fp,"\n");
36 }
37 #endif
38 /*****************************************************************************/
39 /*****************************************************************************/
40 /*****************************************************************************/
41 
42 // One Dimensional Root Finder
43 /*
44  *
45  * vcs_root1d:
46  *
47  *
48  *
49  * Following is a nontrial example for vcs_root1d() where the buoyancy of a
50  * cylinder floating on water is calculated.
51  *
52  * @verbatim
53  * #include <cmath>
54  * #include <cstdlib>
55  *
56  * #include "equil/vcs_internal.h"
57  *
58  * const double g_cgs = 980.;
59  * const double mass_cyl = 0.066;
60  * const double diam_cyl = 0.048;
61  * const double rad_cyl = diam_cyl / 2.0;
62  * const double len_cyl = 5.46;
63  * const double vol_cyl = Pi * diam_cyl * diam_cyl / 4 * len_cyl;
64  * const double rho_cyl = mass_cyl / vol_cyl;
65  * const double rho_gas = 0.0;
66  * const double rho_liq = 1.0;
67  * const double sigma = 72.88;
68  * // Contact angle in radians
69  * const double alpha1 = 40.0 / 180. * Pi;
70  *
71  * using namespace Cantera;
72  * using namespace VCSnonideal;
73  *
74  * double func_vert(double theta1, double h_2, double rho_c) {
75  * double f_grav = - Pi * rad_cyl * rad_cyl * rho_c * g_cgs;
76  * double tmp = rad_cyl * rad_cyl * g_cgs;
77  * double tmp1 = theta1 + sin(theta1) * cos(theta1) - 2.0 * h_2 / rad_cyl * sin(theta1);
78  * double f_buoy = tmp * (Pi * rho_gas + (rho_liq - rho_gas) * tmp1);
79  * double f_sten = 2 * sigma * sin(theta1 + alpha1 - Pi);
80  * double f_net = f_grav + f_buoy + f_sten;
81  * return f_net;
82  * }
83  * double calc_h2_farfield(double theta1) {
84  * double rhs = sigma * (1.0 + cos(alpha1 + theta1));
85  * rhs *= 2.0;
86  * rhs = rhs / (rho_liq - rho_gas) / g_cgs;
87  * double sign = -1.0;
88  * if (alpha1 + theta1 < Pi) sign = 1.0;
89  * double res = sign * sqrt(rhs);
90  * double h2 = res + rad_cyl * cos(theta1);
91  * return h2;
92  * }
93  * double funcZero(double xval, double Vtarget, int varID, void *fptrPassthrough, int *err) {
94  * double theta = xval;
95  * double h2 = calc_h2_farfield(theta);
96  * double fv = func_vert(theta, h2, rho_cyl);
97  * return fv;
98  * }
99  *
100  * int main () {
101  *
102  * double thetamax = Pi;
103  * double thetamin = 0.0;
104  * int maxit = 1000;
105  * int iconv;
106  * double thetaR = Pi/2.0;
107  * int printLvl = 4;
108  *
109  * iconv = VCSnonideal::vcsUtil_root1d(thetamin, thetamax, maxit, funcZero,
110  * (void *) 0, 0.0, 0, &thetaR, printLvl);
111  * printf("theta = %g\n", thetaR);
112  * double h2Final = calc_h2_farfield(thetaR);
113  * printf("h2Final = %g\n", h2Final);
114  * return 0;
115  * }
116  * @endverbatim
117  *
118  */
119 int vcsUtil_root1d(double xmin, double xmax, size_t itmax,
120  VCS_FUNC_PTR func, void* fptrPassthrough,
121  double FuncTargVal, int varID,
122  double* xbest, int printLvl)
123 {
124  static int callNum = 0;
125  const char* stre = "vcs_root1d ERROR: ";
126  const char* strw = "vcs_root1d WARNING: ";
127  bool converged = false;
128  int err = 0;
129 #ifdef DEBUG_MODE
130  char fileName[80];
131  FILE* fp = 0;
132 #endif
133  double x1, x2, xnew, f1, f2, fnew, slope;
134  size_t its = 0;
135  int posStraddle = 0;
136  int retn = VCS_SUCCESS;
137  bool foundPosF = false;
138  bool foundNegF = false;
139  bool foundStraddle = false;
140  double xPosF = 0.0;
141  double xNegF = 0.0;
142  double fnorm; /* A valid norm for the making the function value
143  * dimensionless */
144  double c[9], f[3], xn1, xn2, x0 = 0.0, f0 = 0.0, root, theta, xquad;
145 
146  callNum++;
147 #ifdef DEBUG_MODE
148  if (printLvl >= 3) {
149  sprintf(fileName, "rootfd_%d.log", callNum);
150  fp = fopen(fileName, "w");
151  fprintf(fp, " Iter TP_its xval Func_val | Reasoning\n");
152  fprintf(fp, "-----------------------------------------------------"
153  "-------------------------------\n");
154  }
155 #else
156  if (printLvl >= 3) {
157  plogf("WARNING: vcsUtil_root1d: printlvl >= 3, but debug mode not turned on\n");
158  }
159 #endif
160  if (xmax <= xmin) {
161  plogf("%sxmin and xmax are bad: %g %g\n", stre, xmin, xmax);
162  return VCS_PUB_BAD;
163  }
164  x1 = *xbest;
165  if (x1 < xmin || x1 > xmax) {
166  x1 = (xmin + xmax) / 2.0;
167  }
168  f1 = func(x1, FuncTargVal, varID, fptrPassthrough, &err);
169 #ifdef DEBUG_MODE
170  if (printLvl >= 3) {
171  print_funcEval(fp, x1, f1, its);
172  fprintf(fp, "%-5d %-5d %-15.5E %-15.5E\n", -2, 0, x1, f1);
173  }
174 #endif
175  if (f1 == 0.0) {
176  *xbest = x1;
177  return VCS_SUCCESS;
178  } else if (f1 > 0.0) {
179  foundPosF = true;
180  xPosF = x1;
181  } else {
182  foundNegF = true;
183  xNegF = x1;
184  }
185  x2 = x1 * 1.1;
186  if (x2 > xmax) {
187  x2 = x1 - (xmax - xmin) / 100.;
188  }
189  f2 = func(x2, FuncTargVal, varID, fptrPassthrough, &err);
190 #ifdef DEBUG_MODE
191  if (printLvl >= 3) {
192  print_funcEval(fp, x2, f2, its);
193  fprintf(fp, "%-5d %-5d %-15.5E %-15.5E", -1, 0, x2, f2);
194  }
195 #endif
196 
197  if (FuncTargVal != 0.0) {
198  fnorm = fabs(FuncTargVal) + 1.0E-13;
199  } else {
200  fnorm = 0.5*(fabs(f1) + fabs(f2)) + fabs(FuncTargVal);
201  }
202 
203  if (f2 == 0.0) {
204  return retn;
205  } else if (f2 > 0.0) {
206  if (!foundPosF) {
207  foundPosF = true;
208  xPosF = x2;
209  }
210  } else {
211  if (!foundNegF) {
212  foundNegF = true;
213  xNegF = x2;
214  }
215  }
216  foundStraddle = foundPosF && foundNegF;
217  if (foundStraddle) {
218  if (xPosF > xNegF) {
219  posStraddle = true;
220  } else {
221  posStraddle = false;
222  }
223  }
224 
225  do {
226  /*
227  * Find an estimate of the next point to try based on
228  * a linear approximation.
229  */
230  slope = (f2 - f1) / (x2 - x1);
231  if (slope == 0.0) {
232  plogf("%s functions evals produced the same result, %g, at %g and %g\n",
233  strw, f2, x1, x2);
234  xnew = 2*x2 - x1 + 1.0E-3;
235  } else {
236  xnew = x2 - f2 / slope;
237  }
238 #ifdef DEBUG_MODE
239  if (printLvl >= 3) {
240  fprintf(fp, " | xlin = %-9.4g", xnew);
241  }
242 #endif
243 
244  /*
245  * Do a quadratic fit -> Note this algorithm seems
246  * to work OK. The quadratic approximation doesn't kick in until
247  * the end of the run, when it becomes reliable.
248  */
249  if (its > 0) {
250  c[0] = 1.;
251  c[1] = 1.;
252  c[2] = 1.;
253  c[3] = x0;
254  c[4] = x1;
255  c[5] = x2;
256  c[6] = SQUARE(x0);
257  c[7] = SQUARE(x1);
258  c[8] = SQUARE(x2);
259  f[0] = - f0;
260  f[1] = - f1;
261  f[2] = - f2;
262  retn = vcsUtil_mlequ(c, 3, 3, f, 1);
263  if (retn == 1) {
264  goto QUAD_BAIL;
265  }
266  root = f[1]* f[1] - 4.0 * f[0] * f[2];
267  if (root >= 0.0) {
268  xn1 = (- f[1] + sqrt(root)) / (2.0 * f[2]);
269  xn2 = (- f[1] - sqrt(root)) / (2.0 * f[2]);
270  if (fabs(xn2 - x2) < fabs(xn1 - x2) && xn2 > 0.0) {
271  xquad = xn2;
272  } else {
273  xquad = xn1;
274  }
275  theta = fabs(xquad - xnew) / fabs(xnew - x2);
276  theta = std::min(1.0, theta);
277  xnew = theta * xnew + (1.0 - theta) * xquad;
278 #ifdef DEBUG_MODE
279  if (printLvl >= 3) {
280  if (theta != 1.0) {
281  fprintf(fp, " | xquad = %-9.4g", xnew);
282  }
283  }
284 #endif
285  } else {
286  /*
287  * Pick out situations where the convergence may be
288  * accelerated.
289  */
290  if ((DSIGN(xnew - x2) == DSIGN(x2 - x1)) &&
291  (DSIGN(x2 - x1) == DSIGN(x1 - x0))) {
292  xnew += xnew - x2;
293 #ifdef DEBUG_MODE
294  if (printLvl >= 3) {
295  fprintf(fp, " | xquada = %-9.4g", xnew);
296  }
297 #endif
298  }
299  }
300  }
301 QUAD_BAIL:
302  ;
303 
304 
305  /*
306  *
307  * Put heuristic bounds on the step jump
308  */
309  if ((xnew > x1 && xnew < x2) || (xnew < x1 && xnew > x2)) {
310  /*
311  *
312  * If we are doing a jump in between two points, make sure
313  * the new trial is between 10% and 90% of the distance
314  * between the old points.
315  */
316  slope = fabs(x2 - x1) / 10.;
317  if (fabs(xnew - x1) < slope) {
318  xnew = x1 + DSIGN(xnew-x1) * slope;
319 #ifdef DEBUG_MODE
320  if (printLvl >= 3) {
321  fprintf(fp, " | x10%% = %-9.4g", xnew);
322  }
323 #endif
324  }
325  if (fabs(xnew - x2) < slope) {
326  xnew = x2 + DSIGN(xnew-x2) * slope;
327 #ifdef DEBUG_MODE
328  if (printLvl >= 3) {
329  fprintf(fp, " | x10%% = %-9.4g", xnew);
330  }
331 #endif
332  }
333  } else {
334  /*
335  * If we are venturing into new ground, only allow the step jump
336  * to increase by 100% at each iteration
337  */
338  slope = 2.0 * fabs(x2 - x1);
339  if (fabs(slope) < fabs(xnew - x2)) {
340  xnew = x2 + DSIGN(xnew-x2) * slope;
341 #ifdef DEBUG_MODE
342  if (printLvl >= 3) {
343  fprintf(fp, " | xlimitsize = %-9.4g", xnew);
344  }
345 #endif
346  }
347  }
348 
349 
350  if (xnew > xmax) {
351  xnew = x2 + (xmax - x2) / 2.0;
352 #ifdef DEBUG_MODE
353  if (printLvl >= 3) {
354  fprintf(fp, " | xlimitmax = %-9.4g", xnew);
355  }
356 #endif
357  }
358  if (xnew < xmin) {
359  xnew = x2 + (x2 - xmin) / 2.0;
360 #ifdef DEBUG_MODE
361  if (printLvl >= 3) {
362  fprintf(fp, " | xlimitmin = %-9.4g", xnew);
363  }
364 #endif
365  }
366  if (foundStraddle) {
367 #ifdef DEBUG_MODE
368  slope = xnew;
369 #endif
370  if (posStraddle) {
371  if (f2 > 0.0) {
372  if (xnew > x2) {
373  xnew = (xNegF + x2)/2;
374  }
375  if (xnew < xNegF) {
376  xnew = (xNegF + x2)/2;
377  }
378  } else {
379  if (xnew < x2) {
380  xnew = (xPosF + x2)/2;
381  }
382  if (xnew > xPosF) {
383  xnew = (xPosF + x2)/2;
384  }
385  }
386  } else {
387  if (f2 > 0.0) {
388  if (xnew < x2) {
389  xnew = (xNegF + x2)/2;
390  }
391  if (xnew > xNegF) {
392  xnew = (xNegF + x2)/2;
393  }
394  } else {
395  if (xnew > x2) {
396  xnew = (xPosF + x2)/2;
397  }
398  if (xnew < xPosF) {
399  xnew = (xPosF + x2)/2;
400  }
401  }
402  }
403 #ifdef DEBUG_MODE
404  if (printLvl >= 3) {
405  if (slope != xnew) {
406  fprintf(fp, " | xstraddle = %-9.4g", xnew);
407  }
408  }
409 #endif
410  }
411 
412  fnew = func(xnew, FuncTargVal, varID, fptrPassthrough, &err);
413 #ifdef DEBUG_MODE
414  if (printLvl >= 3) {
415  fprintf(fp,"\n");
416  print_funcEval(fp, xnew, fnew, its);
417  fprintf(fp, "%-5d %-5d %-15.5E %-15.5E", its, 0, xnew, fnew);
418  }
419 #endif
420 
421  if (foundStraddle) {
422  if (posStraddle) {
423  if (fnew > 0.0) {
424  if (xnew < xPosF) {
425  xPosF = xnew;
426  }
427  } else {
428  if (xnew > xNegF) {
429  xNegF = xnew;
430  }
431  }
432  } else {
433  if (fnew > 0.0) {
434  if (xnew > xPosF) {
435  xPosF = xnew;
436  }
437  } else {
438  if (xnew < xNegF) {
439  xNegF = xnew;
440  }
441  }
442  }
443  }
444 
445  if (! foundStraddle) {
446  if (fnew > 0.0) {
447  if (!foundPosF) {
448  foundPosF = true;
449  xPosF = xnew;
450  foundStraddle = true;
451  posStraddle = (xPosF > xNegF);
452  }
453  } else {
454  if (!foundNegF) {
455  foundNegF = true;
456  xNegF = xnew;
457  foundStraddle = true;
458  posStraddle = (xPosF > xNegF);
459  }
460  }
461  }
462 
463  x0 = x1;
464  f0 = f1;
465  x1 = x2;
466  f1 = f2;
467  x2 = xnew;
468  f2 = fnew;
469  if (fabs(fnew / fnorm) < 1.0E-5) {
470  converged = true;
471  }
472  its++;
473  } while (! converged && its < itmax);
474  if (converged) {
475  if (printLvl >= 1) {
476  plogf("vcs_root1d success: convergence achieved\n");
477  }
478 #ifdef DEBUG_MODE
479  if (printLvl >= 3) {
480  fprintf(fp, " | vcs_root1d success in %d its, fnorm = %g\n", its, fnorm);
481  }
482 #endif
483  } else {
484  retn = VCS_FAILED_CONVERGENCE;
485  if (printLvl >= 1) {
486  plogf("vcs_root1d ERROR: maximum iterations exceeded without convergence\n");
487  }
488 #ifdef DEBUG_MODE
489  if (printLvl >= 3) {
490  fprintf(fp, "\nvcs_root1d failure in %d its\n", its);
491  }
492 #endif
493  }
494  *xbest = x2;
495 #ifdef DEBUG_MODE
496  if (printLvl >= 3) {
497  fclose(fp);
498  }
499 #endif
500  return retn;
501 }
502 /*****************************************************************************/
503 }
504 
Generated by   doxygen 1.8.2