Cantera  2.0
vcs_rxnadj.cpp
Go to the documentation of this file.
1 /**
2  * @file vcs_rxnadj.cpp
3  * Routines for carrying out various adjustments to the reaction steps
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_solve.h"
12 #include "cantera/equil/vcs_internal.h"
13 #include "cantera/equil/vcs_VolPhase.h"
14 
15 #include <cstdio>
16 #include <cstdlib>
17 #include <cmath>
18 
19 namespace VCSnonideal
20 {
21 
22 // Calculates formation reaction step sizes.
23 /*
24  * This is equation 6.4-16, p. 143 in Smith and Missen.
25  *
26  * Output
27  * -------
28  * m_deltaMolNumSpecies[kspec] : reaction adjustments, where irxn refers
29  * to the irxn'th species
30  * formation reaction. This adjustment
31  * is for species
32  * irxn + M, where M is the number
33  * of components.
34  *
35  * Special branching occurs sometimes. This causes the component basis
36  * to be reevaluated
37  *
38  * @return Returns an int representing the status of the step
39  * - 0 : normal return
40  * - 1 : A single species phase species has been zeroed out
41  * in this routine. The species is a noncomponent
42  * - 2 : Same as one but, the zeroed species is a component.
43  */
44 size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial)
45 {
46  size_t kspec, iph;
47  size_t iphDel = npos;
48  double s, xx, dss;
49  size_t k = 0;
50  vcs_VolPhase* Vphase = 0;
51  double* dnPhase_irxn;
52 #ifdef DEBUG_MODE
53  char ANOTE[128];
54  if (m_debug_print_lvl >= 2) {
55  plogf(" ");
56  for (int j = 0; j < 82; j++) {
57  plogf("-");
58  }
59  plogf("\n");
60  plogf(" --- Subroutine vcs_RxnStepSizes called - Details:\n");
61  plogf(" ");
62  for (int j = 0; j < 82; j++) {
63  plogf("-");
64  }
65  plogf("\n");
66  plogf(" --- Species KMoles Rxn_Adjustment DeltaG"
67  " | Comment\n");
68  }
69 #endif
70  /*
71  * We update the matrix dlnActCoeffdmolNumber[][] at the
72  * top of the loop, when necessary
73  */
74  if (m_useActCoeffJac) {
75  vcs_CalcLnActCoeffJac(VCS_DATA_PTR(m_molNumSpecies_old));
76  }
77  /************************************************************************
78  ******** LOOP OVER THE FORMATION REACTIONS *****************************
79  ************************************************************************/
80 
81  for (size_t irxn = 0; irxn < m_numRxnRdc; ++irxn) {
82 #ifdef DEBUG_MODE
83  sprintf(ANOTE,"Normal Calc");
84 #endif
85 
86  kspec = m_indexRxnToSpecies[irxn];
87 
88  if (m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDPHASE) {
89  m_deltaMolNumSpecies[kspec] = 0.0;
90 #ifdef DEBUG_MODE
91  sprintf(ANOTE, "ZeroedPhase: Phase is artificially zeroed");
92 #endif
93  } else if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
94 
95  dnPhase_irxn = m_deltaMolNumPhase[irxn];
96 
97  if (m_molNumSpecies_old[kspec] == 0.0 && (! m_SSPhase[kspec])) {
98  /********************************************************************/
99  /******* MULTISPECIES PHASE WITH total moles equal to zero *********/
100  /*******************************************************************/
101  /*
102  * If dg[irxn] is negative, then the multispecies phase should
103  * come alive again. Add a small positive step size to
104  * make it come alive.
105  */
106  if (m_deltaGRxn_new[irxn] < -1.0e-4) {
107  /*
108  * First decide if this species is part of a multiphase that
109  * is nontrivial in size.
110  */
111  iph = m_phaseID[kspec];
112  double tphmoles = m_tPhaseMoles_old[iph];
113  double trphmoles = tphmoles / m_totalMolNum;
114  if (trphmoles > VCS_DELETE_PHASE_CUTOFF) {
115  m_deltaMolNumSpecies[kspec] = m_totalMolNum * VCS_SMALL_MULTIPHASE_SPECIES;
116  if (m_speciesStatus[kspec] == VCS_SPECIES_STOICHZERO) {
117  m_deltaMolNumSpecies[kspec] = 0.0;
118 #ifdef DEBUG_MODE
119  sprintf(ANOTE,
120  "MultSpec (%s): Species not born due to STOICH/PHASEPOP even though DG = %11.3E",
121  vcs_speciesType_string(m_speciesStatus[kspec], 15),
122  m_deltaGRxn_new[irxn]);
123 #endif
124  } else {
125  m_deltaMolNumSpecies[kspec] = m_totalMolNum * VCS_SMALL_MULTIPHASE_SPECIES;
126 #ifdef DEBUG_MODE
127  sprintf(ANOTE,
128  "MultSpec (%s): small species born again DG = %11.3E",
129  vcs_speciesType_string(m_speciesStatus[kspec], 15),
130  m_deltaGRxn_new[irxn]);
131 #endif
132  }
133  } else {
134 #ifdef DEBUG_MODE
135  sprintf(ANOTE, "MultSpec (%s): phase come alive DG = %11.3E",
136  vcs_speciesType_string(m_speciesStatus[kspec], 15),
137  m_deltaGRxn_new[irxn]);
138 #endif
139  Vphase = m_VolPhaseList[iph];
140  size_t numSpPhase = Vphase->nSpecies();
141  m_deltaMolNumSpecies[kspec] =
142  m_totalMolNum * 10.0 * VCS_DELETE_PHASE_CUTOFF / numSpPhase;
143  }
144 
145  } else {
146 #ifdef DEBUG_MODE
147  sprintf(ANOTE, "MultSpec (%s): still dead DG = %11.3E",
148  vcs_speciesType_string(m_speciesStatus[kspec], 15),
149  m_deltaGRxn_new[irxn]);
150 #endif
151  m_deltaMolNumSpecies[kspec] = 0.0;
152  }
153  } else {
154  /********************************************************************/
155  /************************* REGULAR PROCESSING ***********************/
156  /********************************************************************/
157  /*
158  * First take care of cases where we want to bail out
159  *
160  *
161  * Don't bother if superconvergence has already been achieved
162  * in this mode.
163  */
164  if (fabs(m_deltaGRxn_new[irxn]) <= m_tolmaj2) {
165 #ifdef DEBUG_MODE
166  sprintf(ANOTE,"Skipped: superconverged DG = %11.3E", m_deltaGRxn_new[irxn]);
167  if (m_debug_print_lvl >= 2) {
168  plogf(" --- %-12.12s", m_speciesName[kspec].c_str());
169  plogf(" %12.4E %12.4E %12.4E | %s\n",
170  m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec],
171  m_deltaGRxn_new[irxn], ANOTE);
172  }
173 #endif
174  continue;
175  }
176  /*
177  * Don't calculate for minor or nonexistent species if
178  * their values are to be decreasing anyway.
179  */
180  if ((m_speciesStatus[kspec] != VCS_SPECIES_MAJOR) && (m_deltaGRxn_new[irxn] >= 0.0)) {
181 #ifdef DEBUG_MODE
182  sprintf(ANOTE,"Skipped: IC = %3d and DG >0: %11.3E",
183  m_speciesStatus[kspec], m_deltaGRxn_new[irxn]);
184  if (m_debug_print_lvl >= 2) {
185  plogf(" --- %-12.12s", m_speciesName[kspec].c_str());
186  plogf(" %12.4E %12.4E %12.4E | %s\n",
187  m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec],
188  m_deltaGRxn_new[irxn], ANOTE);
189  }
190 #endif
191  continue;
192  }
193  /*
194  * Start of the regular processing
195  */
196  if (m_SSPhase[kspec]) {
197  s = 0.0;
198  } else {
199  s = 1.0 / m_molNumSpecies_old[kspec] ;
200  }
201  for (size_t j = 0; j < m_numComponents; ++j) {
202  if (!m_SSPhase[j]) {
203  if (m_molNumSpecies_old[j] > 0.0) {
204  s += SQUARE(m_stoichCoeffRxnMatrix[irxn][j]) / m_molNumSpecies_old[j];
205  }
206  }
207  }
208  for (size_t j = 0; j < m_numPhases; j++) {
209  Vphase = m_VolPhaseList[j];
210  if (! Vphase->m_singleSpecies) {
211  if (m_tPhaseMoles_old[j] > 0.0) {
212  s -= SQUARE(dnPhase_irxn[j]) / m_tPhaseMoles_old[j];
213  }
214  }
215  }
216  if (s != 0.0) {
217  /*
218  * Take into account of the
219  * derivatives of the activity coefficients with respect to the
220  * mole numbers, even in our diagonal approximation.
221  */
222  if (m_useActCoeffJac) {
223  double s_old = s;
224  s = vcs_Hessian_diag_adj(irxn, s_old);
225 #ifdef DEBUG_MODE
226  if (s_old != s) {
227  sprintf(ANOTE, "Normal calc: diag adjusted from %g "
228  "to %g due to act coeff", s_old, s);
229  }
230 #endif
231  }
232 
233  m_deltaMolNumSpecies[kspec] = -m_deltaGRxn_new[irxn] / s;
234  // New section to do damping of the m_deltaMolNumSpecies[]
235  for (size_t j = 0; j < m_numComponents; ++j) {
236  double stoicC = m_stoichCoeffRxnMatrix[irxn][j];
237  if (stoicC != 0.0) {
238  double negChangeComp = - stoicC * m_deltaMolNumSpecies[kspec];
239  if (negChangeComp > m_molNumSpecies_old[j]) {
240  if (m_molNumSpecies_old[j] > 0.0) {
241 #ifdef DEBUG_MODE
242  sprintf(ANOTE, "Delta damped from %g "
243  "to %g due to component %d (%10s) going neg", m_deltaMolNumSpecies[kspec],
244  -m_molNumSpecies_old[j]/stoicC, j, m_speciesName[j].c_str());
245 #endif
246  m_deltaMolNumSpecies[kspec] = - m_molNumSpecies_old[j] / stoicC;
247  } else {
248 #ifdef DEBUG_MODE
249  sprintf(ANOTE, "Delta damped from %g "
250  "to %g due to component %d (%10s) zero", m_deltaMolNumSpecies[kspec],
251  -m_molNumSpecies_old[j]/stoicC, j, m_speciesName[j].c_str());
252 #endif
253  m_deltaMolNumSpecies[kspec] = 0.0;
254  }
255  }
256  }
257  }
258  // Implement a damping term that limits m_deltaMolNumSpecies to the size of the mole number
259  if (-m_deltaMolNumSpecies[kspec] > m_molNumSpecies_old[kspec]) {
260 #ifdef DEBUG_MODE
261  sprintf(ANOTE, "Delta damped from %g "
262  "to %g due to %s going negative", m_deltaMolNumSpecies[kspec],
263  -m_molNumSpecies_old[kspec], m_speciesName[kspec].c_str());
264 #endif
265  m_deltaMolNumSpecies[kspec] = -m_molNumSpecies_old[kspec];
266  }
267 
268  } else {
269  /* ************************************************************ */
270  /* **** REACTION IS ENTIRELY AMONGST SINGLE SPECIES PHASES **** */
271  /* **** DELETE ONE OF THE PHASES AND RECOMPUTE BASIS ********* */
272  /* ************************************************************ */
273  /*
274  * Either the species L will disappear or one of the
275  * component single species phases will disappear. The sign
276  * of DG(I) will indicate which way the reaction will go.
277  * Then, we need to follow the reaction to see which species
278  * will zero out first.
279  * -> The species to be zeroed out will be "k".
280  */
281  if (m_deltaGRxn_new[irxn] > 0.0) {
282  dss = m_molNumSpecies_old[kspec];
283  k = kspec;
284  for (size_t j = 0; j < m_numComponents; ++j) {
285  if (m_stoichCoeffRxnMatrix[irxn][j] > 0.0) {
286  xx = m_molNumSpecies_old[j] / m_stoichCoeffRxnMatrix[irxn][j];
287  if (xx < dss) {
288  dss = xx;
289  k = j;
290  }
291  }
292  }
293  dss = -dss;
294  } else {
295  dss = 1.0e10;
296  for (size_t j = 0; j < m_numComponents; ++j) {
297  if (m_stoichCoeffRxnMatrix[irxn][j] < 0.0) {
298  xx = -m_molNumSpecies_old[j] / m_stoichCoeffRxnMatrix[irxn][j];
299  if (xx < dss) {
300  dss = xx;
301  k = j;
302  }
303  }
304  }
305  }
306  /*
307  * Here we adjust the mole fractions
308  * according to DSS and the stoichiometric array
309  * to take into account that we are eliminating
310  * the kth species. DSS contains the amount
311  * of moles of the kth species that needs to be
312  * added back into the component species.
313  */
314  if (dss != 0.0) {
315 
316  if ((k == kspec) && (m_SSPhase[kspec] != 1)) {
317  /*
318  * Found out that we can be in this spot, when components of multispecies phases
319  * are zeroed, leaving noncomponent species of the same phase having all of the
320  * mole numbers of that phases. it seems that we can suggest a zero of the species
321  * and the code will recover.
322  */
323 #ifdef DEBUG_MODE
324  sprintf(ANOTE, "Delta damped from %g to %g due to delete %s",
325  m_deltaMolNumSpecies[kspec],
326  -m_molNumSpecies_old[kspec], m_speciesName[kspec].c_str());
327 #endif
328  m_deltaMolNumSpecies[kspec] = -m_molNumSpecies_old[kspec];
329 #ifdef DEBUG_MODE
330  if (m_debug_print_lvl >= 2) {
331  plogf(" --- %-12.12s", m_speciesName[kspec].c_str());
332  plogf(" %12.4E %12.4E %12.4E | %s\n",
333  m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec],
334  m_deltaGRxn_new[irxn], ANOTE);
335  }
336 #endif
337  continue;
338  }
339  /*
340  * Delete the single species phase
341  */
342  for (size_t j = 0; j < m_numSpeciesTot; j++) {
343  m_deltaMolNumSpecies[j] = 0.0;
344  }
345  m_deltaMolNumSpecies[kspec] = dss;
346  for (size_t j = 0; j < m_numComponents; ++j) {
347  m_deltaMolNumSpecies[j] = dss * m_stoichCoeffRxnMatrix[irxn][j];
348  }
349 
350  iphDel = m_phaseID[k];
351  kSpecial = k;
352 
353 #ifdef DEBUG_MODE
354  if (k != kspec) {
355  sprintf(ANOTE, "Delete component SS phase %d named %s - SS phases only",
356  iphDel, m_speciesName[k].c_str());
357  } else {
358  sprintf(ANOTE, "Delete this SS phase %d - SS components only", iphDel);
359  }
360  if (m_debug_print_lvl >= 2) {
361  plogf(" --- %-12.12s", m_speciesName[kspec].c_str());
362  plogf(" %12.4E %12.4E %12.4E | %s\n",
363  m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec],
364  m_deltaGRxn_new[irxn], ANOTE);
365  plogf(" --- vcs_RxnStepSizes Special section to set up to delete %s",
366  m_speciesName[k].c_str());
367  plogendl();
368  }
369 #endif
370  if (k != kspec) {
371  forceComponentCalc = 1;
372 #ifdef DEBUG_MODE
373  if (m_debug_print_lvl >= 2) {
374  plogf(" --- Force a component recalculation \n");
375  plogendl();
376  }
377 #endif
378  }
379 #ifdef DEBUG_MODE
380  if (m_debug_print_lvl >= 2) {
381  plogf(" ");
382  vcs_print_line("-", 82);
383  }
384 #endif
385  return iphDel;
386  }
387  }
388  } /* End of regular processing */
389 #ifdef DEBUG_MODE
390  if (m_debug_print_lvl >= 2) {
391  plogf(" --- %-12.12s", m_speciesName[kspec].c_str());
392  plogf(" %12.4E %12.4E %12.4E | %s\n",
393  m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec],
394  m_deltaGRxn_new[irxn], ANOTE);
395  }
396 #endif
397  } /* End of loop over m_speciesUnknownType */
398  } /* End of loop over non-component stoichiometric formation reactions */
399 #ifdef DEBUG_MODE
400  if (m_debug_print_lvl >= 2) {
401  plogf(" ");
402  vcs_print_line("-", 82);
403  }
404 #endif
405  return iphDel;
406 }
407 
408 //====================================================================================================================
409 // Calculates reaction adjustments using a full Hessian approximation
410 /*
411  * Calculates reaction adjustments. This does what equation 6.4-16, p. 143
412  * in Smith and Missen is suppose to do. However, a full matrix is
413  * formed and then solved via a conjugate gradient algorithm. No
414  * preconditioning is done.
415  *
416  * If special branching is warranted, then the program bails out.
417  *
418  * Output
419  * -------
420  * DS(I) : reaction adjustment, where I refers to the Ith species
421  * Special branching occurs sometimes. This causes the component basis
422  * to be reevaluated
423  * return = 0 : normal return
424  * 1 : A single species phase species has been zeroed out
425  * in this routine. The species is a noncomponent
426  * 2 : Same as one but, the zeroed species is a component.
427  *
428  * Special attention is taken to flag cases where the direction of the
429  * update is contrary to the steepest descent rule. This is an important
430  * attribute of the regular vcs algorithm. We don't want to violate this.
431  *
432  * NOTE: currently this routine is not used.
433  */
434 int VCS_SOLVE::vcs_rxn_adj_cg()
435 {
436  size_t irxn, j;
437  size_t k = 0;
438  size_t kspec;
439  int soldel = 0;
440  double s, xx, dss;
441  double* dnPhase_irxn;
442 #ifdef DEBUG_MODE
443  char ANOTE[128];
444  plogf(" ");
445  for (j = 0; j < 77; j++) {
446  plogf("-");
447  }
448  plogf("\n --- Subroutine rxn_adj_cg() called\n");
449  plogf(" --- Species Moles Rxn_Adjustment | Comment\n");
450 #endif
451 
452  /*
453  * Precalculation loop -> we calculate quantities based on
454  * loops over the number of species.
455  * We also evaluate whether the matrix is appropriate for
456  * this algorithm. If not, we bail out.
457  */
458  for (irxn = 0; irxn < m_numRxnRdc; ++irxn) {
459 #ifdef DEBUG_MODE
460  sprintf(ANOTE,"Normal Calc");
461 #endif
462 
463  kspec = m_indexRxnToSpecies[irxn];
464  dnPhase_irxn = m_deltaMolNumPhase[irxn];
465 
466  if (m_molNumSpecies_old[kspec] == 0.0 && (! m_SSPhase[kspec])) {
467  /* *******************************************************************/
468  /* **** MULTISPECIES PHASE WITH total moles equal to zero ************/
469  /* *******************************************************************/
470  /*
471  * HKM -> the statment below presupposes units in m_deltaGRxn_new[]. It probably
472  * should be replaced with something more relativistic
473  */
474  if (m_deltaGRxn_new[irxn] < -1.0e-4) {
475 #ifdef DEBUG_MODE
476  (void) sprintf(ANOTE, "MultSpec: come alive DG = %11.3E", m_deltaGRxn_new[irxn]);
477 #endif
478  m_deltaMolNumSpecies[kspec] = 1.0e-10;
479  m_speciesStatus[kspec] = VCS_SPECIES_MAJOR;
480  --(m_numRxnMinorZeroed);
481  } else {
482 #ifdef DEBUG_MODE
483  (void) sprintf(ANOTE, "MultSpec: still dead DG = %11.3E", m_deltaGRxn_new[irxn]);
484 #endif
485  m_deltaMolNumSpecies[kspec] = 0.0;
486  }
487  } else {
488  /* ********************************************** */
489  /* **** REGULAR PROCESSING ********** */
490  /* ********************************************** */
491  /*
492  * First take care of cases where we want to bail out
493  *
494  *
495  * Don't bother if superconvergence has already been achieved
496  * in this mode.
497  */
498  if (fabs(m_deltaGRxn_new[irxn]) <= m_tolmaj2) {
499 #ifdef DEBUG_MODE
500  sprintf(ANOTE,"Skipped: converged DG = %11.3E\n", m_deltaGRxn_new[irxn]);
501  plogf(" --- ");
502  plogf("%-12.12s", m_speciesName[kspec].c_str());
503  plogf(" %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec],
504  m_deltaMolNumSpecies[kspec], ANOTE);
505 #endif
506  continue;
507  }
508  /*
509  * Don't calculate for minor or nonexistent species if
510  * their values are to be decreasing anyway.
511  */
512  if (m_speciesStatus[kspec] <= VCS_SPECIES_MINOR && m_deltaGRxn_new[irxn] >= 0.0) {
513 #ifdef DEBUG_MODE
514  sprintf(ANOTE,"Skipped: IC = %3d and DG >0: %11.3E\n",
515  m_speciesStatus[kspec], m_deltaGRxn_new[irxn]);
516  plogf(" --- ");
517  plogf("%-12.12s", m_speciesName[kspec].c_str());
518  plogf(" %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec],
519  m_deltaMolNumSpecies[kspec], ANOTE);
520 #endif
521  continue;
522  }
523  /*
524  * Start of the regular processing
525  */
526  if (m_SSPhase[kspec]) {
527  s = 0.0;
528  } else {
529  s = 1.0 / m_molNumSpecies_old[kspec];
530  }
531  for (j = 0; j < m_numComponents; ++j) {
532  if (! m_SSPhase[j]) {
533  s += SQUARE(m_stoichCoeffRxnMatrix[irxn][j]) / m_molNumSpecies_old[j];
534  }
535  }
536  for (j = 0; j < m_numPhases; j++) {
537  if (!(m_VolPhaseList[j])->m_singleSpecies) {
538  if (m_tPhaseMoles_old[j] > 0.0) {
539  s -= SQUARE(dnPhase_irxn[j]) / m_tPhaseMoles_old[j];
540  }
541  }
542  }
543  if (s != 0.0) {
544  m_deltaMolNumSpecies[kspec] = -m_deltaGRxn_new[irxn] / s;
545  } else {
546  /* ************************************************************ */
547  /* **** REACTION IS ENTIRELY AMONGST SINGLE SPECIES PHASES **** */
548  /* **** DELETE ONE SOLID AND RECOMPUTE BASIS ********* */
549  /* ************************************************************ */
550  /*
551  * Either the species L will disappear or one of the
552  * component single species phases will disappear. The sign
553  * of DG(I) will indicate which way the reaction will go.
554  * Then, we need to follow the reaction to see which species
555  * will zero out first.
556  */
557  if (m_deltaGRxn_new[irxn] > 0.0) {
558  dss = m_molNumSpecies_old[kspec];
559  k = kspec;
560  for (j = 0; j < m_numComponents; ++j) {
561  if (m_stoichCoeffRxnMatrix[irxn][j] > 0.0) {
562  xx = m_molNumSpecies_old[j] / m_stoichCoeffRxnMatrix[irxn][j];
563  if (xx < dss) {
564  dss = xx;
565  k = j;
566  }
567  }
568  }
569  dss = -dss;
570  } else {
571  dss = 1.0e10;
572  for (j = 0; j < m_numComponents; ++j) {
573  if (m_stoichCoeffRxnMatrix[irxn][j] < 0.0) {
574  xx = -m_molNumSpecies_old[j] / m_stoichCoeffRxnMatrix[irxn][j];
575  if (xx < dss) {
576  dss = xx;
577  k = j;
578  }
579  }
580  }
581  }
582  /*
583  * Here we adjust the mole fractions
584  * according to DSS and the stoichiometric array
585  * to take into account that we are eliminating
586  * the kth species. DSS contains the amount
587  * of moles of the kth species that needs to be
588  * added back into the component species.
589  */
590  if (dss != 0.0) {
591  m_molNumSpecies_old[kspec] += dss;
592  m_tPhaseMoles_old[m_phaseID[kspec]] += dss;
593  for (j = 0; j < m_numComponents; ++j) {
594  m_molNumSpecies_old[j] += dss * m_stoichCoeffRxnMatrix[irxn][j];
595  m_tPhaseMoles_old[m_phaseID[j]] += dss * m_stoichCoeffRxnMatrix[irxn][j];
596  }
597  m_molNumSpecies_old[k] = 0.0;
598  m_tPhaseMoles_old[m_phaseID[k]] = 0.0;
599 #ifdef DEBUG_MODE
600  plogf(" --- vcs_st2 Special section to delete ");
601  plogf("%-12.12s", m_speciesName[k].c_str());
602  plogf("\n --- Immediate return - Restart iteration\n");
603 #endif
604  /*
605  * We need to immediately recompute the
606  * component basis, because we just zeroed
607  * it out.
608  */
609  if (k != kspec) {
610  soldel = 2;
611  } else {
612  soldel = 1;
613  }
614  return soldel;
615  }
616  }
617  } /* End of regular processing */
618 #ifdef DEBUG_MODE
619  plogf(" --- ");
620  plogf("%-12.12s", m_speciesName[kspec].c_str());
621  plogf(" %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec],
622  m_deltaMolNumSpecies[kspec], ANOTE);
623 #endif
624  } /* End of loop over non-component stoichiometric formation reactions */
625 
626  /*
627  *
628  * When we form the Hessian we must be careful to ensure that it
629  * is a symmetric positive definate matrix, still. This means zeroing
630  * out columns when we zero out rows as well.
631  * -> I suggest writing a small program to make sure of this
632  * property.
633  */
634 
635 #ifdef DEBUG_MODE
636  plogf(" ");
637  for (j = 0; j < 77; j++) {
638  plogf("-");
639  }
640  plogf("\n");
641 #endif
642  return soldel;
643 }
644 
645 //====================================================================================================================
646 // Calculates the diagonal contribution to the Hessian due to
647 // the dependence of the activity coefficients on the mole numbers.
648 /*
649  * (See framemaker notes, Eqn. 20 - VCS Equations document)
650  *
651  * We allow the diagonal to be increased positively to any degree.
652  * We allow the diagonal to be decreased to 1/3 of the ideal solution
653  * value, but no more -> it must remain positive.
654  *
655  * NOTE: currently this routine is not used
656  */
657 double VCS_SOLVE::vcs_Hessian_diag_adj(size_t irxn, double hessianDiag_Ideal)
658 {
659  double diag = hessianDiag_Ideal;
660  double hessActCoef = vcs_Hessian_actCoeff_diag(irxn);
661  if (hessianDiag_Ideal <= 0.0) {
662  plogf("vcs_Hessian_diag_adj::We shouldn't be here\n");
663  exit(EXIT_FAILURE);
664  }
665  if (hessActCoef >= 0.0) {
666  diag += hessActCoef;
667  } else if (fabs(hessActCoef) < 0.6666 * hessianDiag_Ideal) {
668  diag += hessActCoef;
669  } else {
670  diag -= 0.6666 * hessianDiag_Ideal;
671  }
672  return diag;
673 }
674 
675 //====================================================================================================================
676 //! Calculates the diagonal contribution to the Hessian due to
677 //! the dependence of the activity coefficients on the mole numbers.
678 /*!
679  * (See framemaker notes, Eqn. 20 - VCS Equations document)
680  *
681  * NOTE: currently this routine is not used
682  */
683 double VCS_SOLVE::vcs_Hessian_actCoeff_diag(size_t irxn)
684 {
685  size_t kspec, k, l, kph;
686  double s;
687  double* sc_irxn;
688  kspec = m_indexRxnToSpecies[irxn];
689  kph = m_phaseID[kspec];
690  sc_irxn = m_stoichCoeffRxnMatrix[irxn];
691  /*
692  * First the diagonal term of the Jacobian
693  */
694  s = m_dLnActCoeffdMolNum[kspec][kspec];
695  /*
696  * Next, the other terms. Note this only a loop over the components
697  * So, it's not too expensive to calculate.
698  */
699  for (l = 0; l < m_numComponents; l++) {
700  if (!m_SSPhase[l]) {
701  for (k = 0; k < m_numComponents; ++k) {
702  if (m_phaseID[k] == m_phaseID[l]) {
703  s += sc_irxn[k] * sc_irxn[l] * m_dLnActCoeffdMolNum[k][l];
704  }
705  }
706  if (kph == m_phaseID[l]) {
707  s += sc_irxn[l] * (m_dLnActCoeffdMolNum[kspec][l] + m_dLnActCoeffdMolNum[l][kspec]);
708  }
709  }
710  }
711  return s;
712 }
713 //====================================================================================================================
714 // Recalculate all of the activity coefficients in all of the phases
715 // based on input mole numbers
716 /*
717  *
718  * @param moleSpeciesVCS kmol of species to be used in the update.
719  *
720  * NOTE: This routine needs to be regulated.
721  */
722 void VCS_SOLVE::vcs_CalcLnActCoeffJac(const double* const moleSpeciesVCS)
723 {
724  /*
725  * Loop over all of the phases in the problem
726  */
727  for (size_t iphase = 0; iphase < m_numPhases; iphase++) {
728  vcs_VolPhase* Vphase = m_VolPhaseList[iphase];
729  /*
730  * We don't need to call single species phases;
731  */
732  if (!Vphase->m_singleSpecies && !Vphase->isIdealSoln()) {
733  /*
734  * update the mole numbers
735  */
736  Vphase->setMolesFromVCS(VCS_STATECALC_OLD, moleSpeciesVCS);
737  /*
738  * Download the resulting calculation into the full vector
739  * -> This scatter calculation is carried out in the
740  * vcs_VolPhase object.
741  */
742  Vphase->sendToVCS_LnActCoeffJac(m_dLnActCoeffdMolNum.baseDataAddr());
743  }
744  }
745 }
746 /*****************************************************************************/
747 
748 //! This function recalculates the deltaG for reaction, irxn
749 /*!
750  * This function recalculates the deltaG for reaction irxn,
751  * given the mole numbers in molNum. It uses the temporary
752  * space mu_i, to hold the recalculated chemical potentials.
753  * It only recalculates the chemical potentials for species in phases
754  * which participate in the irxn reaction.
755  *
756  * This function is used by the vcs_line_search algorithm() and
757  * should not be used widely due to the unknown state it leaves the
758  * system.
759  *
760  * Input
761  * ------------
762  * @param irxn Reaction number
763  * @param molNum Current mole numbers of species to be used as
764  * input to the calculation (units = kmol)
765  * (length = totalNuMSpecies)
766  *
767  * Output
768  * ------------
769  * @param ac output Activity coefficients (length = totalNumSpecies)
770  * Note this is only partially formed. Only species in
771  * phases that participate in the reaction will be updated
772  * @param mu_i dimensionless chemical potentials (length - totalNumSpecies
773  * Note this is only partially formed. Only species in
774  * phases that participate in the reaction will be updated
775  *
776  * @return Returns the dimensionless deltaG of the reaction
777  *
778  * Note, this is a dangerous routine that leaves the underlying objects in
779  * an unknown state.
780  */
781 double VCS_SOLVE::deltaG_Recalc_Rxn(const int stateCalc,
782  const size_t irxn, const double* const molNum,
783  double* const ac, double* const mu_i)
784 {
785  size_t kspec = irxn + m_numComponents;
786  int* pp_ptr = m_phaseParticipation[irxn];
787  for (size_t iphase = 0; iphase < m_numPhases; iphase++) {
788  if (pp_ptr[iphase]) {
789  vcs_chemPotPhase(stateCalc, iphase, molNum, ac, mu_i);
790  }
791  }
792  double deltaG = mu_i[kspec];
793  double* sc_irxn = m_stoichCoeffRxnMatrix[irxn];
794  for (size_t k = 0; k < m_numComponents; k++) {
795  deltaG += sc_irxn[k] * mu_i[k];
796  }
797  return deltaG;
798 }
799 /*****************************************************************************/
800 
801 #ifdef DEBUG_MODE
802 // A line search algorithm is carried out on one reaction
803 /*
804  * In this routine we carry out a rough line search algorithm
805  * to make sure that the m_deltaGRxn_new doesn't switch signs prematurely.
806  *
807  * @param irxn Reaction number
808  * @param dx_orig Original step length
809  *
810  * @param ANOTE Output character string stating the conclusions of the
811  * line search
812  *
813  * @return Returns the optimized step length found by the search
814  */
815 double VCS_SOLVE::vcs_line_search(const size_t irxn, const double dx_orig,
816  char* const ANOTE)
817 #else
818 double VCS_SOLVE::vcs_line_search(const size_t irxn, const double dx_orig)
819 #endif
820 {
821  int its = 0;
822  size_t k;
823  size_t kspec = m_indexRxnToSpecies[irxn];
824  const int MAXITS = 10;
825  double dx = dx_orig;
826  double* sc_irxn = m_stoichCoeffRxnMatrix[irxn];
827  double* molNumBase = VCS_DATA_PTR(m_molNumSpecies_old);
828  double* acBase = VCS_DATA_PTR(m_actCoeffSpecies_old);
829  double* ac = VCS_DATA_PTR(m_actCoeffSpecies_new);
830  double molSum = 0.0;
831  double slope;
832  /*
833  * Calculate the deltaG value at the dx = 0.0 point
834  */
835  vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD);
836  double deltaGOrig = deltaG_Recalc_Rxn(VCS_STATECALC_OLD,
837  irxn, molNumBase, acBase,
838  VCS_DATA_PTR(m_feSpecies_old));
839  double forig = fabs(deltaGOrig) + 1.0E-15;
840  if (deltaGOrig > 0.0) {
841  if (dx_orig > 0.0) {
842  dx = 0.0;
843 #ifdef DEBUG_MODE
844  if (m_debug_print_lvl >= 2) {
845  //plogf(" --- %s :Warning possible error dx>0 dg > 0\n", SpName[kspec]);
846  }
847  sprintf(ANOTE,"Rxn reduced to zero step size in line search: dx>0 dg > 0");
848 #endif
849  return dx;
850  }
851  } else if (deltaGOrig < 0.0) {
852  if (dx_orig < 0.0) {
853  dx = 0.0;
854 #ifdef DEBUG_MODE
855  if (m_debug_print_lvl >= 2) {
856  //plogf(" --- %s :Warning possible error dx<0 dg < 0\n", SpName[kspec]);
857  }
858  sprintf(ANOTE,"Rxn reduced to zero step size in line search: dx<0 dg < 0");
859 #endif
860  return dx;
861  }
862  } else if (deltaGOrig == 0.0) {
863  return 0.0;
864  }
865  if (dx_orig == 0.0) {
866  return 0.0;
867  }
868 
869  vcs_dcopy(VCS_DATA_PTR(m_molNumSpecies_new), molNumBase, m_numSpeciesRdc);
870  molSum = molNumBase[kspec];
871  m_molNumSpecies_new[kspec] = molNumBase[kspec] + dx_orig;
872  for (k = 0; k < m_numComponents; k++) {
873  m_molNumSpecies_new[k] = molNumBase[k] + sc_irxn[k] * dx_orig;
874  molSum += molNumBase[k];
875  }
876  vcs_setFlagsVolPhases(false, VCS_STATECALC_NEW);
877 
878  double deltaG1 = deltaG_Recalc_Rxn(VCS_STATECALC_NEW,
879  irxn, VCS_DATA_PTR(m_molNumSpecies_new),
880  ac, VCS_DATA_PTR(m_feSpecies_new));
881 
882  /*
883  * If deltaG hasn't switched signs when going the full distance
884  * then we are heading in the appropriate direction, and
885  * we should accept the current full step size
886  */
887  if (deltaG1 * deltaGOrig > 0.0) {
888  dx = dx_orig;
889  goto finalize;
890  }
891  /*
892  * If we have decreased somewhat, the deltaG return after finding
893  * a better estimate for the line search.
894  */
895  if (fabs(deltaG1) < 0.8*forig) {
896  if (deltaG1 * deltaGOrig < 0.0) {
897  slope = (deltaG1 - deltaGOrig) / dx_orig;
898  dx = -deltaGOrig / slope;
899  } else {
900  dx = dx_orig;
901  }
902  goto finalize;
903  }
904 
905  dx = dx_orig;
906 
907  for (its = 0; its < MAXITS; its++) {
908  /*
909  * Calculate the approximation to the total Gibbs free energy at
910  * the dx *= 0.5 point
911  */
912  dx *= 0.5;
913  m_molNumSpecies_new[kspec] = molNumBase[kspec] + dx;
914  for (k = 0; k < m_numComponents; k++) {
915  m_molNumSpecies_new[k] = molNumBase[k] + sc_irxn[k] * dx;
916  }
917  vcs_setFlagsVolPhases(false, VCS_STATECALC_NEW);
918  double deltaG = deltaG_Recalc_Rxn(VCS_STATECALC_NEW,
919  irxn, VCS_DATA_PTR(m_molNumSpecies_new),
920  ac, VCS_DATA_PTR(m_feSpecies_new));
921  /*
922  * If deltaG hasn't switched signs when going the full distance
923  * then we are heading in the appropriate direction, and
924  * we should accept the current step
925  */
926  if (deltaG * deltaGOrig > 0.0) {
927  goto finalize;
928  }
929  /*
930  * If we have decreased somewhat, the deltaG return after finding
931  * a better estimate for the line search.
932  */
933  if (fabs(deltaG) / forig < (1.0 - 0.1 * dx / dx_orig)) {
934  if (deltaG * deltaGOrig < 0.0) {
935  slope = (deltaG - deltaGOrig) / dx;
936  dx = -deltaGOrig / slope;
937  }
938  goto finalize;
939  }
940  }
941 
942 finalize:
943  vcs_setFlagsVolPhases(false, VCS_STATECALC_NEW);
944  if (its >= MAXITS) {
945 #ifdef DEBUG_MODE
946  sprintf(ANOTE,"Rxn reduced to zero step size from %g to %g (MAXITS)",
947  dx_orig, dx);
948  return dx;
949 #endif
950  }
951 #ifdef DEBUG_MODE
952  if (dx != dx_orig) {
953  sprintf(ANOTE,"Line Search reduced step size from %g to %g",
954  dx_orig, dx);
955  }
956 #endif
957 
958  return dx;
959 }
960 /*****************************************************************************/
961 }
962 
Generated by   doxygen 1.8.2