1 // This file is part of Cantera. See License.txt in the top-level directory or
2 // at https://cantera.org/license.txt for license and copyright information.
3
4 #include "cantera/equil/vcs_MultiPhaseEquil.h"
5
6 #include "cantera/thermo/ThermoFactory.h"
7 #include "cantera/thermo/IdealGasPhase.h"
8 #include "cantera/thermo/HMWSoln.h"
9 #include "cantera/thermo/StoichSubstance.h"
10
11 using namespace Cantera;
12 using namespace std;
13
printUsage()14 void printUsage()
15 {
16 cout << "usage: nacl_equil [-h] [-help_cmdfile] [-d #] [HMW_NaCl.yaml] "
17 << endl;
18 cout << " -h help" << endl;
19 cout << " -d # : level of debug printing" << endl;
20 cout << " [HMW_NaCl.yaml] - Optionally change the name of the input file " << endl;
21 cout << endl;
22 cout << endl;
23 }
24
main(int argc,char ** argv)25 int main(int argc, char** argv)
26 {
27 #if defined(_MSC_VER) && _MSC_VER < 1900
28 _set_output_format(_TWO_DIGIT_EXPONENT);
29 #endif
30 suppress_deprecation_warnings();
31 int numSucc = 0;
32 int numFail = 0;
33 int printLvl = 1;
34 string inputFile = "HMW_NaCl.yaml";
35 VCS_SOLVE::disableTiming();
36
37 /*
38 * Process the command line arguments
39 */
40 if (argc > 1) {
41 string tok;
42 for (int j = 1; j < argc; j++) {
43 tok = string(argv[j]);
44 if (tok[0] == '-') {
45 int nopt = static_cast<int>(tok.size());
46 for (int n = 1; n < nopt; n++) {
47 if (!strcmp(tok.c_str() + 1, "help_cmdfile")) {
48 } else if (tok[n] == 'h') {
49 printUsage();
50 exit(1);
51 } else if (tok[n] == 'd') {
52 printLvl = 2;
53 int lvl = 2;
54 if (j < (argc - 1)) {
55 string tokla = string(argv[j+1]);
56 if (strlen(tokla.c_str()) > 0) {
57 lvl = atoi(tokla.c_str());
58 n = nopt - 1;
59 j += 1;
60 if (lvl >= 0) {
61 printLvl = lvl;
62 }
63 }
64 }
65 } else {
66 printUsage();
67 exit(1);
68 }
69 }
70 } else if (inputFile == "HMW_NaCl.yaml") {
71 inputFile = tok;
72 } else {
73 printUsage();
74 exit(1);
75 }
76 }
77 }
78
79
80
81 try {
82 int estimateEquil = -1;
83 double T = 298.15;
84 double pres = OneAtm;
85
86 // Initialize the individual phases
87
88 HMWSoln hmw(inputFile, "NaCl_electrolyte_complex_shomate");
89 hmw.setName("NaCl_electrolyte");
90 size_t kk = hmw.nSpecies();
91 vector_fp Xmol(kk, 0.0);
92 size_t iH2OL = hmw.speciesIndex("H2O(L)");
93 Xmol[iH2OL] = 1.0;
94 hmw.setState_TPX(T, pres, Xmol.data());
95
96 ThermoPhase* gas = newPhase("NaCl_gas.yaml");
97
98 kk = gas->nSpecies();
99 Xmol.resize(kk, 0.0);
100 for (size_t i = 0; i < kk; i++) {
101 Xmol[i] = 0.0;
102 }
103 size_t iN2 = gas->speciesIndex("N2");
104 Xmol[iN2] = 1.0;
105 gas->setState_TPX(T, pres, Xmol.data());
106
107
108 StoichSubstance ss("NaCl_Solid.yaml");
109 ss.setState_TP(T, pres);
110
111
112 // Construct the multiphase object
113 MultiPhase* mp = new MultiPhase();
114
115 mp->addPhase(&hmw, 2.0);
116 mp->addPhase(gas, 4.0);
117 mp->addPhase(&ss, 5.0);
118
119
120 try {
121 mp->equilibrate("TP", "vcs", 1e-9, 50000, 100, estimateEquil, printLvl);
122 cout << *mp;
123 numSucc++;
124 } catch (CanteraError& err) {
125 cout << *mp;
126 std::cerr << err.what() << std::endl;
127 cerr << "ERROR: MultiEquil equilibration step failed at "
128 << " T = " << T
129 << " Pres = " << pres
130 << endl;
131 cout << "ERROR: MultiEqiul equilibration step failed at "
132 << " T = " << T
133 << " Pres = " << pres
134 << endl;
135 exit(-1);
136 }
137
138 cout << "NUMBER OF SUCCESSES = " << numSucc << endl;
139 cout << "NUMBER OF FAILURES = " << numFail << endl;
140
141 return numFail;
142 } catch (CanteraError& err) {
143 std::cerr << err.what() << std::endl;
144 cerr << "ERROR: program terminating due to unforeseen circumstances." << endl;
145 return -1;
146 }
147 }
148