Cantera  2.0
vcs_internal.h
Go to the documentation of this file.
1 /**
2  * @file vcs_internal.h
3  * Internal declarations for the VCSnonideal package
4  */
5 /*
6  * Copyright (2005) 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 #ifndef _VCS_INTERNAL_H
12 #define _VCS_INTERNAL_H
13 
14 #include <cstring>
15 
16 #include "cantera/equil/vcs_defs.h"
17 
18 #include "cantera/base/global.h"
19 
20 #ifndef ALTLINPROG
21 #define ALTLINPROG
22 #endif
23 
24 namespace VCSnonideal
25 {
26 using Cantera::npos;
27 
28 //! Points to the data in a std::vector<> object
29 #define VCS_DATA_PTR(vvv) (&(vvv[0]))
30 
31 //! define this Cantera function to replace printf
32 /*!
33  * We can replace this with printf easily
34  */
35 #define plogf Cantera::writelogf
36 
37 //! define this Cantera function to replace cout << endl;
38 /*!
39  * We use this to place an endl in the log file, and
40  * ensure that the IO buffers are flushed.
41  */
42 #define plogendl() Cantera::writelogendl()
43 
44 //! Global hook for turning on and off time printing.
45 /*!
46  * Default is to allow printing. But, you can assign this to zero
47  * globally to turn off all time printing.
48  * This is helpful for test suite purposes where you are interested
49  * in differences in text files.
50  */
51 extern int vcs_timing_print_lvl;
52 
53 /*
54  * Forward references
55  */
56 class VCS_SPECIES_THERMO;
57 class VCS_PROB;
58 
59 //! Amount of extra printing that is done while in debug mode.
60 /*!
61  * 0 -> none
62  * 1 -> some
63  * 2 -> alot (default)
64  * 3 -> everything
65  */
66 
67 //! Class to keep track of time and iterations
68 /*!
69  * class keeps all of the counters together.
70  */
72 {
73 public:
74  //! Total number of iterations in the main loop
75  //! of vcs_TP() to solve for thermo equilibrium
76  int T_Its;
77 
78  //! Current number of iterations in the main loop
79  //! of vcs_TP() to solve for thermo equilibrium
80  int Its;
81 
82  //! Total number of optimizations of the
83  //! components basis set done
85 
86  //! number of optimizations of the components basis set done
88 
89  //! Current number of times the initial thermo
90  //! equilibrium estimator has been called
92 
93  //! Current number of calls to vcs_TP
95 
96  //! Current time spent in vcs_TP
97  double T_Time_vcs_TP;
98 
99  //! Current time spent in vcs_TP
100  double Time_vcs_TP;
101 
102  //! Total Time spent in basopt
104 
105  //! Current Time spent in basopt
106  double Time_basopt;
107 
108  //! Time spent in initial estimator
109  double T_Time_inest;
110 
111  //! Time spent in the vcs suite of programs
112  double T_Time_vcs;
113 };
114 
115 //! Returns the value of the gas constant in the units specified by parameter
116 /*!
117  * @param mu_units Specifies the units.
118  * - VCS_UNITS_KCALMOL: kcal gmol-1 K-1
119  * - VCS_UNITS_UNITLESS: 1.0 K-1
120  * - VCS_UNITS_KJMOL: kJ gmol-1 K-1
121  * - VCS_UNITS_KELVIN: 1.0 K-1
122  * - VCS_UNITS_MKS: joules kmol-1 K-1 = kg m2 s-2 kmol-1 K-1
123  */
124 double vcsUtil_gasConstant(int mu_units);
125 
126 //! Invert an n x n matrix and solve m rhs's
127 /*!
128  * Solve a square matrix with multiple right hand sides
129  *
130  * \f[
131  * C X + B = 0;
132  * \f]
133  *
134  * This routine uses Gauss elimination and is optimized for the solution
135  * of lots of rhs's. A crude form of row pivoting is used here.
136  * The matrix C is destroyed during the solve.
137  *
138  * @return The solution x[] is returned in the matrix <I>B</I>.
139  * Routine returns an integer representing success:
140  * - 1 : Matrix is singular
141  * - 0 : solution is OK
142  *
143  *
144  * @param c Matrix to be inverted. c is in fortran format, i.e., rows
145  * are the inner loop. Row numbers equal to idem.
146  * c[i+j*idem] = c_i_j = Matrix to be inverted:
147  * - i = row number
148  * - j = column number
149  *
150  * @param idem number of row dimensions in c
151  * @param n Number of rows and columns in c
152  * @param b Multiple RHS. Note, b is actually the negative of
153  * most formulations. Row numbers equal to idem.
154  * b[i+j*idem] = b_i_j = vectors of rhs's:
155  * - i = row number
156  * - j = column number
157  * (each column is a new rhs)
158  * @param m number of rhs's
159  */
160 int vcsUtil_mlequ(double* c, size_t idem, size_t n, double* b, size_t m);
161 
162 //! Invert an n x n matrix and solve m rhs's
163 /*!
164  * Solve a square matrix with multiple right hand sides
165  *
166  * \f[
167  * C X + B = 0;
168  * \f]
169  *
170  * This routine uses Gauss-Jordan elimination and is optimized for the solution
171  * of lots of rhs's. Full row and column pivoting is used here. It's been
172  * shown to be necessary in at least one case.
173  * The matrix C is destroyed during the solve.
174  *
175  * @return The solution x[] is returned in the matrix <I>B</I>.
176  * Routine returns an integer representing success:
177  * - 1 : Matrix is singular
178  * - 0 : solution is OK
179  *
180  * @param c Matrix to be inverted. c is in fortran format, i.e., rows
181  * are the inner loop. Row numbers equal to idem.
182  * c[i+j*idem] = c_i_j = Matrix to be inverted:
183  * - i = row number
184  * - j = column number
185  *
186  * @param idem number of row dimensions in c
187  * @param n Number of rows and columns in c
188  * @param b Multiple RHS. Note, b is actually the negative of
189  * most formulations. Row numbers equal to idem.
190  * b[i+j*idem] = b_i_j = vectors of rhs's:
191  * - i = row number
192  * - j = column number
193  * (each column is a new rhs)
194  * @param m number of rhs's
195  */
196 int vcsUtil_gaussj(double* c, size_t idem, size_t n, double* b, size_t m);
197 
198 
199 //! Definition of the function pointer for the root finder
200 /*!
201  * see vcsUtil_root1d for a definition of how to use this.
202  */
203 typedef double(*VCS_FUNC_PTR)(double xval, double Vtarget,
204  int varID, void* fptrPassthrough,
205  int* err);
206 
207 //! One dimensional root finder
208 /*!
209  *
210  * This root finder will find the root of a one dimensional
211  * equation
212  *
213  * \f[
214  * f(x) = 0
215  * \f]
216  * where x is a bounded quantity: \f$ x_{min} < x < x_max \f$
217  *
218  * The functional to be minimized must have the following call
219  * structure:
220  *
221  * @verbatim
222  typedef double (*VCS_FUNC_PTR)(double xval, double Vtarget,
223  int varID, void *fptrPassthrough,
224  int *err); @endverbatim
225  *
226  * xval is the current value of the x variable. Vtarget is the
227  * requested value of f(x), usually 0. varID is an integer
228  * that is passed through. fptrPassthrough is a void pointer
229  * that is passed through. err is a return error indicator.
230  * err = 0 is the norm. anything else is considered a fatal
231  * error.
232  * The return value of the function is the current value of
233  * f(xval).
234  *
235  * @param xmin Minimum permissible value of the x variable
236  * @param xmax Maximum permissible value of the x parameter
237  * @param itmax Maximum number of iterations
238  * @param func function pointer, pointing to the function to be
239  * minimized
240  * @param fptrPassthrough Pointer to void that gets passed through
241  * the rootfinder, unchanged, to the func.
242  * @param FuncTargVal Target value of the function. This is usually set
243  * to zero.
244  * @param varID Variable ID. This is usually set to zero.
245  * @param xbest Pointer to the initial value of x on input. On output
246  * This contains the root value.
247  * @param printLvl Print level of the routine.
248  *
249  *
250  * Following is a nontrial example for vcs_root1d() in which the position of a
251  * cylinder floating on the water is calculated.
252  *
253  * @verbatim
254  #include <cmath>
255  #include <cstdlib>
256 
257  #include "equil/vcs_internal.h"
258 
259  const double g_cgs = 980.;
260  const double mass_cyl = 0.066;
261  const double diam_cyl = 0.048;
262  const double rad_cyl = diam_cyl / 2.0;
263  const double len_cyl = 5.46;
264  const double vol_cyl = Pi * diam_cyl * diam_cyl / 4 * len_cyl;
265  const double rho_cyl = mass_cyl / vol_cyl;
266  const double rho_gas = 0.0;
267  const double rho_liq = 1.0;
268  const double sigma = 72.88;
269  // Contact angle in radians
270  const double alpha1 = 40.0 / 180. * Pi;
271 
272  double func_vert(double theta1, double h_2, double rho_c) {
273  double f_grav = - Pi * rad_cyl * rad_cyl * rho_c * g_cgs;
274  double tmp = rad_cyl * rad_cyl * g_cgs;
275  double tmp1 = theta1 + sin(theta1) * cos(theta1) - 2.0 * h_2 / rad_cyl * sin(theta1);
276  double f_buoy = tmp * (Pi * rho_gas + (rho_liq - rho_gas) * tmp1);
277  double f_sten = 2 * sigma * sin(theta1 + alpha1 - Pi);
278  double f_net = f_grav + f_buoy + f_sten;
279  return f_net;
280  }
281  double calc_h2_farfield(double theta1) {
282  double rhs = sigma * (1.0 + cos(alpha1 + theta1));
283  rhs *= 2.0;
284  rhs = rhs / (rho_liq - rho_gas) / g_cgs;
285  double sign = -1.0;
286  if (alpha1 + theta1 < Pi) sign = 1.0;
287  double res = sign * sqrt(rhs);
288  double h2 = res + rad_cyl * cos(theta1);
289  return h2;
290  }
291  double funcZero(double xval, double Vtarget, int varID, void *fptrPassthrough, int *err) {
292  double theta = xval;
293  double h2 = calc_h2_farfield(theta);
294  double fv = func_vert(theta, h2, rho_cyl);
295  return fv;
296  }
297  int main () {
298  double thetamax = Pi;
299  double thetamin = 0.0;
300  int maxit = 1000;
301  int iconv;
302  double thetaR = Pi/2.0;
303  int printLvl = 4;
304 
305  iconv = VCSnonideal::vcsUtil_root1d(thetamin, thetamax, maxit,
306  funcZero,
307  (void *) 0, 0.0, 0,
308  &thetaR, printLvl);
309  printf("theta = %g\n", thetaR);
310  double h2Final = calc_h2_farfield(thetaR);
311  printf("h2Final = %g\n", h2Final);
312  return 0;
313  } @endverbatim
314  *
315  */
316 int vcsUtil_root1d(double xmin, double xmax, size_t itmax, VCS_FUNC_PTR func,
317  void* fptrPassthrough,
318  double FuncTargVal, int varID, double* xbest,
319  int printLvl = 0);
320 
321 //! Returns the system wall clock time in seconds
322 /*!
323  * @return time in seconds.
324  */
325 double vcs_second();
326 
327 //! This define turns on using memset and memcpy. I have not run into
328 //! any systems where this is a problem. It's the fastest way to do
329 //! low lvl operations where applicable. There are alternative routines
330 //! available if this ever fails.
331 #define USE_MEMSET
332 #ifdef USE_MEMSET
333 
334 //! Zero a double vector
335 /*!
336  * @param vec_to vector of doubles
337  * @param length length of the vector to zero.
338  */
339 inline void vcs_dzero(double* const vec_to, const size_t length)
340 {
341  (void) memset((void*) vec_to, 0, length * sizeof(double));
342 }
343 
344 //! Zero an int vector
345 /*!
346  * @param vec_to vector of ints
347  * @param length length of the vector to zero.
348  */
349 inline void vcs_izero(int* const vec_to, const size_t length)
350 {
351  (void) memset((void*) vec_to, 0, length * sizeof(int));
352 }
353 
354 //! Copy a double vector
355 /*!
356  * @param vec_to Vector to copy into. This vector must be dimensioned
357  * at least as large as the vec_from vector.
358  * @param vec_from Vector to copy from
359  * @param length Number of doubles to copy.
360  */
361 inline void vcs_dcopy(double* const vec_to,
362  const double* const vec_from, const size_t length)
363 {
364  (void) memcpy((void*) vec_to, (const void*) vec_from,
365  (length) * sizeof(double));
366 }
367 
368 
369 //! Copy an int vector
370 /*!
371  * @param vec_to Vector to copy into. This vector must be dimensioned
372  * at least as large as the vec_from vector.
373  * @param vec_from Vector to copy from
374  * @param length Number of int to copy.
375  */
376 inline void vcs_icopy(int* const vec_to,
377  const int* const vec_from, const size_t length)
378 {
379  (void) memcpy((void*) vec_to, (const void*) vec_from,
380  (length) * sizeof(int));
381 }
382 
383 //! Zero a std double vector
384 /*!
385  * @param vec_to vector of doubles
386  * @param length length of the vector to zero.
387  */
388 inline void vcs_vdzero(std::vector<double> &vec_to, const size_t length)
389 {
390  (void) memset((void*)VCS_DATA_PTR(vec_to), 0, (length) * sizeof(double));
391 }
392 
393 //! Zero a std int vector
394 /*!
395  * @param vec_to vector of ints
396  * @param length length of the vector to zero.
397  */
398 inline void vcs_vizero(std::vector<int> &vec_to, const size_t length)
399 {
400  (void) memset((void*)VCS_DATA_PTR(vec_to), 0, (length) * sizeof(int));
401 }
402 
403 //! Copy one std double vector into another
404 /*!
405  * This is an inlined function that uses memcpy. memcpy is probably
406  * the fastest way to do this. This routine requires the vectors to be
407  * previously dimensioned appropriately. No error checking is done.
408  *
409  * @param vec_to Vector to copy into. This vector must be dimensioned
410  * at least as large as the vec_from vector.
411  * @param vec_from Vector to copy from
412  * @param length Number of doubles to copy.
413  */
414 inline void vcs_vdcopy(std::vector<double> & vec_to,
415  const std::vector<double> & vec_from, size_t length)
416 {
417  (void) memcpy((void*)&(vec_to[0]), (const void*) &(vec_from[0]),
418  (length) * sizeof(double));
419 }
420 
421 //! Copy one std integer vector into another
422 /*!
423  * This is an inlined function that uses memcpy. memcpy is probably
424  * the fastest way to do this. This routine requires the
425  *
426  * @param vec_to Vector to copy into. This vector must be dimensioned
427  * at least as large as the vec_from vector.
428  * @param vec_from Vector to copy from
429  * @param length Number of integers to copy.
430  */
431 inline void vcs_vicopy(std::vector<int> & vec_to,
432  const std::vector<int> & vec_from, const int length)
433 {
434  (void) memcpy((void*)&(vec_to[0]), (const void*) &(vec_from[0]),
435  (length) * sizeof(int));
436 }
437 #else
438 extern void vcs_dzero(double* const, const int);
439 extern void vcs_izero(int* const , const int);
440 extern void vcs_dcopy(double* const, const double* const, const int);
441 extern void vcs_icopy(int* const, const int* const, const int);
442 extern void vcs_vdzero(std::vector<double> &vvv, const int len = -1);
443 extern void vcs_vizero(std::vector<double> &vvv, const int len = -1);
444 void vcs_vdcopy(std::vector<double> &vec_to,
445  const std::vector<double> vec_from, const int len = -1);
446 void vcs_vicopy(std::vector<int> &vec_to,
447  const std::vector<int> vec_from, const int len = -1);
448 #endif
449 
450 //! determine the l2 norm of a vector of doubles
451 /*!
452  * @param vec vector of doubles
453  *
454  * @return Returns the l2 norm of the vector
455  */
456 double vcs_l2norm(const std::vector<double> vec);
457 
458 //! Finds the location of the maximum component in a double vector
459 /*!
460  * @param x pointer to a vector of doubles
461  * @param xSize pointer to a vector of doubles used as a multiplier
462  * to x[]
463  * @param j lowest index to search from
464  * @param n highest index to search from
465  * @return Return index of the greatest value on X(i) searched
466  * j <= i < n
467  */
468 size_t vcs_optMax(const double* x, const double* xSize, size_t j, size_t n);
469 
470 //! Returns the maximum integer in a list
471 /*!
472  * @param vector pointer to a vector of ints
473  * @param length length of the integer vector
474  *
475  * @return returns the max integer value in the list
476  */
477 int vcs_max_int(const int* vector, int length);
478 
479 //! Prints a line consisting of multiple occurrences of the same string
480 /*!
481  * This prints a string num times, and then terminate with a
482  * end of line character
483  *
484  * @param str C string that is null terminated
485  * @param num number of times the string is to be printed
486  */
487 void vcs_print_line(const char* str, int num);
488 
489 //! Returns a const char string representing the type of the
490 //! species given by the first argument
491 /*!
492  * @param speciesStatus Species status integer representing the type
493  * of the species.
494  * @param length Maximum length of the string to be returned.
495  * Shorter values will yield abbreviated strings.
496  * Defaults to a value of 100.
497  */
498 const char* vcs_speciesType_string(int speciesStatus, int length = 100);
499 
500 //! Print a string within a given space limit
501 /*!
502  * This routine limits the amount of the string that will be printed to a
503  * maximum of "space" characters. Printing is done to
504  * to Cantera's writelog() function.
505  *
506  * @param str String, which must be null terminated.
507  * @param space space limit for the printing.
508  * @param alignment Alignment of string within the space:
509  * - 0 centered
510  * - 1 right aligned
511  * - 2 left aligned
512  */
513 void vcs_print_stringTrunc(const char* str, size_t space, int alignment);
514 
515 //! Simple routine to check whether two doubles are equal up to
516 //! roundoff error
517 /*!
518  * Currently it's set to check for 10 digits of
519  * relative accuracy.
520  *
521  * @param d1 first double
522  * @param d2 second double
523  *
524  * @return returns true if the doubles are "equal" and false otherwise
525  */
526 bool vcs_doubleEqual(double d1, double d2);
527 
528 }
529 
530 #endif