Cantera  2.0
newEquilibrate.cpp
1 
2 #include "cantera/equil/equil.h"
4 #include "cantera/equil/MultiPhaseEquil.h"
6 #include "cantera/base/global.h"
8 
9 namespace Cantera
10 {
11 
12 int new_equilibrate(thermo_t& s, const char* XY, int solver,
13  doublereal rtol, int maxsteps, int maxiter, int loglevel)
14 {
15  // Decide which solvers to consider
16  int ixy = _equilflag(XY);
17  vector_int solvers;
18  if (solver == -1) {
19  solvers.push_back(2);
20  if (ixy == TP || ixy == HP || ixy == SP || ixy == TV) {
21  solvers.push_back(1);
22  }
23  solvers.push_back(0);
24  } else if (solver == 1) {
25  if (ixy == TP || ixy == HP || ixy == SP || ixy == TV) {
26  solvers.push_back(1);
27  } else {
28  throw CanteraError("equilibrate",
29  std::string("MultiPhase equilibrium solver does not support property pair ") + XY);
30  }
31  } else if (solver == 0 || solver == 2) {
32  solvers.push_back(0);
33  } else {
34  throw CanteraError("equilibrate",
35  "unknown solver:" + int2str(solver));
36  }
37 
38  int nAttempts = 0;
39  bool success = false;
40  while (!solvers.empty()) {
41  nAttempts++;
42  solver = solvers.back();
43  solvers.pop_back();
44 
45  if (solver == 0) {
46  ChemEquil e;
47  e.options.maxIterations = maxsteps;
48  e.options.relTolerance = rtol;
49  nAttempts++;
50  int retnSub = e.equilibrate(s, XY, true, loglevel-1);
51  if (retnSub < 0) {
52  if (loglevel > 0) {
53  addLogEntry("ChemEquil solver failed.");
54  }
55  } else {
56  if (loglevel > 0) {
57  addLogEntry("ChemEquil solver succeeded.");
58  }
59  return nAttempts;
60  }
61  s.setElementPotentials(e.elementPotentials());
62  } else if (solver == 1) {
63  try {
64  MultiPhase m;
65  m.addPhase(&s, 1.0);
66  m.init();
67  nAttempts++;
68  std::cout << "Calling equilibrate(Multiphase)" << std::endl;
69  equilibrate(m, XY, rtol, maxsteps, maxiter, loglevel-1);
70  if (loglevel > 0) {
71  addLogEntry("MultiPhaseEquil solver succeeded.");
72  }
73  return nAttempts;
74  } catch (CanteraError& err) {
75  err.save();
76  if (loglevel > 0) {
77  addLogEntry("MultiPhaseEquil solver failed.");
78  }
79  if (nAttempts < 2) {
80  if (loglevel > 0) {
81  addLogEntry("Trying single phase ChemEquil solver.");
82  }
83  solver = -1;
84  } else {
85  if (loglevel > 0) {
86  endLogGroup("equilibrate");
87  }
88  throw err;
89  }
90  }
91 
92  }
93 
94  }
95 
96  throw CanteraError("equilibrate", "All available solvers failed");
97 
98 }
99 
100 }