00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016 #ifdef WIN32
00017 #pragma warning(disable:4786)
00018 #endif
00019
00020 #include "ThermoFactory.h"
00021
00022 #include "speciesThermoTypes.h"
00023 #include "SpeciesThermoFactory.h"
00024 #include "IdealGasPhase.h"
00025 #include "VPSSMgr.h"
00026 #include "VPSSMgrFactory.h"
00027
00028 #ifdef WITH_IDEAL_SOLUTIONS
00029 #include "IdealSolidSolnPhase.h"
00030 #include "MargulesVPSSTP.h"
00031 #include "IonsFromNeutralVPSSTP.h"
00032 #endif
00033
00034 #ifdef WITH_PURE_FLUIDS
00035 #include "PureFluidPhase.h"
00036 #endif
00037
00038 #include "ConstDensityThermo.h"
00039 #include "SurfPhase.h"
00040 #include "EdgePhase.h"
00041
00042 #ifdef WITH_METAL
00043 #include "MetalPhase.h"
00044 #endif
00045
00046 #ifdef WITH_SEMICONDUCTOR
00047 #include "SemiconductorPhase.h"
00048 #endif
00049
00050 #undef USE_SSTP
00051 #ifdef WITH_STOICH_SUBSTANCE
00052 #ifdef USE_SSTP
00053 #include "StoichSubstanceSSTP.h"
00054
00055 #else
00056 #include "StoichSubstance.h"
00057 #endif
00058 #endif
00059
00060 #ifdef WITH_STOICH_SUBSTANCE
00061 #include "MineralEQ3.h"
00062 #include "MetalSHEelectrons.h"
00063 #endif
00064
00065
00066
00067 #ifdef WITH_LATTICE_SOLID
00068 #include "LatticeSolidPhase.h"
00069 #include "LatticePhase.h"
00070 #endif
00071
00072 #ifdef WITH_ELECTROLYTES
00073 #include "HMWSoln.h"
00074 #include "DebyeHuckel.h"
00075 #include "IdealMolalSoln.h"
00076 #endif
00077
00078 #include "IdealSolnGasVPSS.h"
00079
00080 #include <cstdlib>
00081
00082 using namespace std;
00083
00084 namespace Cantera {
00085
00086 ThermoFactory* ThermoFactory::s_factory = 0;
00087 #if defined(THREAD_SAFE_CANTERA)
00088 boost::mutex ThermoFactory::thermo_mutex;
00089 #endif
00090
00091 static int ntypes = 18;
00092 static string _types[] = {"IdealGas", "Incompressible",
00093 "Surface", "Edge", "Metal", "StoichSubstance",
00094 "PureFluid", "LatticeSolid", "Lattice",
00095 "HMW", "IdealSolidSolution", "DebyeHuckel",
00096 "IdealMolalSolution", "IdealGasVPSS",
00097 "MineralEQ3", "MetalSHEelectrons", "Margules",
00098 "IonsFromNeutralMolecule"
00099 };
00100
00101 static int _itypes[] = {cIdealGas, cIncompressible,
00102 cSurf, cEdge, cMetal, cStoichSubstance,
00103 cPureFluid, cLatticeSolid, cLattice,
00104 cHMW, cIdealSolidSolnPhase, cDebyeHuckel,
00105 cIdealMolalSoln, cVPSS_IdealGas,
00106 cMineralEQ3, cMetalSHEelectrons,
00107 cMargulesVPSSTP, cIonsFromNeutral
00108 };
00109
00110
00111
00112
00113 ThermoPhase* ThermoFactory::newThermoPhase(std::string model) {
00114
00115 int ieos=-1;
00116
00117 for (int n = 0; n < ntypes; n++) {
00118 if (model == _types[n]) ieos = _itypes[n];
00119 }
00120
00121 ThermoPhase* th=0;
00122 switch (ieos) {
00123
00124 case cIdealGas:
00125 th = new IdealGasPhase;
00126 break;
00127
00128 case cIncompressible:
00129 th = new ConstDensityThermo;
00130 break;
00131
00132 case cSurf:
00133 th = new SurfPhase;
00134 break;
00135
00136 case cEdge:
00137 th = new EdgePhase;
00138 break;
00139
00140 #ifdef WITH_IDEAL_SOLUTIONS
00141 case cIdealSolidSolnPhase:
00142 th = new IdealSolidSolnPhase();
00143 break;
00144
00145 case cMargulesVPSSTP:
00146 th = new MargulesVPSSTP();
00147 break;
00148
00149 case cIonsFromNeutral:
00150 th = new IonsFromNeutralVPSSTP();
00151 break;
00152 #endif
00153
00154 #ifdef WITH_METAL
00155 case cMetal:
00156 th = new MetalPhase;
00157 break;
00158 #endif
00159
00160 #ifdef WITH_STOICH_SUBSTANCE
00161 case cStoichSubstance:
00162 #ifdef USE_SSTP
00163 th = new StoichSubstanceSSTP;
00164 #else
00165 th = new StoichSubstance;
00166 #endif
00167 break;
00168 #endif
00169
00170 #ifdef WITH_STOICH_SUBSTANCE
00171 case cMineralEQ3:
00172 th = new MineralEQ3();
00173 break;
00174 #endif
00175
00176 #ifdef WITH_STOICH_SUBSTANCE
00177 case cMetalSHEelectrons:
00178 th = new MetalSHEelectrons();
00179 break;
00180 #endif
00181
00182 #ifdef WITH_LATTICE_SOLID
00183 case cLatticeSolid:
00184 th = new LatticeSolidPhase;
00185 break;
00186
00187 case cLattice:
00188 th = new LatticePhase;
00189 break;
00190 #endif
00191
00192 #ifdef WITH_PURE_FLUIDS
00193 case cPureFluid:
00194 th = new PureFluidPhase;
00195 break;
00196 #endif
00197 #ifdef WITH_ELECTROLYTES
00198 case cHMW:
00199 th = new HMWSoln;
00200 break;
00201
00202 case cDebyeHuckel:
00203 th = new DebyeHuckel;
00204 break;
00205
00206 case cIdealMolalSoln:
00207 th = new IdealMolalSoln;
00208 break;
00209 #endif
00210
00211 case cVPSS_IdealGas:
00212 th = new IdealSolnGasVPSS;
00213 break;
00214
00215 default:
00216 throw UnknownThermoPhaseModel("ThermoFactory::newThermoPhase",
00217 model);
00218 }
00219 return th;
00220 }
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230 std::string eosTypeString(int ieos, int length)
00231 {
00232 std::string ss = "UnknownPhaseType";
00233
00234 for (int n = 0; n < ntypes; n++) {
00235 if (_itypes[n] == ieos) {
00236 ss = _types[n];
00237
00238 }
00239 }
00240 return ss;
00241 }
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255 ThermoPhase* newPhase(XML_Node& xmlphase) {
00256 const XML_Node& th = xmlphase.child("thermo");
00257 string model = th["model"];
00258 ThermoPhase* t = newThermoPhase(model);
00259 if (model == "singing cows") {
00260 throw CanteraError(" newPhase", "Cows don't sing");
00261 }
00262 #ifdef WITH_ELECTROLYTES
00263 else if (model == "HMW") {
00264 HMWSoln* p = dynamic_cast<HMWSoln*>(t);
00265 p->constructPhaseXML(xmlphase,"");
00266 }
00267 #endif
00268 #ifdef WITH_IDEAL_SOLUTIONS
00269 else if (model == "IonsFromNeutralMolecule") {
00270 IonsFromNeutralVPSSTP* p = dynamic_cast<IonsFromNeutralVPSSTP*>(t);
00271 p->constructPhaseXML(xmlphase,"");
00272 }
00273 #endif
00274 else {
00275 importPhase(xmlphase, t);
00276 }
00277 return t;
00278 }
00279
00280 ThermoPhase* newPhase(std::string infile, std::string id) {
00281 XML_Node* root = get_XML_File(infile);
00282 if (id == "-") id = "";
00283 XML_Node* xphase = get_XML_NameID("phase", std::string("#")+id, root);
00284 if (!xphase) {
00285 throw CanteraError("newPhase",
00286 "Couldn't find phase named \"" + id + "\" in file, " + infile);
00287 }
00288 if (xphase)
00289 return newPhase(*xphase);
00290 else
00291 return (ThermoPhase *) 0;
00292 }
00293
00294
00295 static void formSpeciesXMLNodeList(std::vector<XML_Node *> &spDataNodeList,
00296 std::vector<std::string> &spNamesList,
00297 std::vector<int> &spRuleList,
00298 const std::vector<XML_Node *> spArray_names,
00299 const std::vector<XML_Node *> spArray_dbases,
00300 const vector_int sprule) {
00301
00302
00303 std::map<std::string, bool> declared;
00304
00305 int nspa = spArray_dbases.size();
00306 int nSpecies = 0;
00307 bool skip;
00308
00309 for (int jsp = 0; jsp < nspa; jsp++) {
00310 const XML_Node& speciesArray = *spArray_names[jsp];
00311
00312
00313 const XML_Node *db = spArray_dbases[jsp];
00314
00315
00316 std::vector<std::string> spnames;
00317 getStringArray(speciesArray, spnames);
00318 int nsp = static_cast<int>(spnames.size());
00319
00320
00321
00322
00323 if (nsp == 1 && spnames[0] == "all") {
00324 std::vector<XML_Node *> allsp;
00325 db->getChildren("species", allsp);
00326 nsp = static_cast<int>(allsp.size());
00327 spnames.resize(nsp);
00328 for (int nn = 0; nn < nsp; nn++) {
00329 string stemp = (*allsp[nn])["name"];
00330 bool skip = false;
00331 if (declared[stemp]) {
00332 if (sprule[jsp] >= 10) {
00333 skip = true;
00334 } else {
00335 throw CanteraError("ThermoFactory::formSpeciesXMLNodeList()",
00336 "duplicate species: \"" + stemp + "\"");
00337 }
00338 }
00339 if (!skip) {
00340 declared[stemp] = true;
00341 nSpecies++;
00342 spNamesList.resize(nSpecies);
00343 spDataNodeList.resize(nSpecies, 0);
00344 spRuleList.resize(nSpecies, 0);
00345 spNamesList[nSpecies-1] = stemp;
00346 spDataNodeList[nSpecies-1] = allsp[nn];
00347 spRuleList[nSpecies-1] = sprule[jsp];
00348 }
00349 }
00350 }
00351 else if (nsp == 1 && spnames[0] == "unique") {
00352 std::vector<XML_Node *> allsp;
00353 db->getChildren("species", allsp);
00354 nsp = static_cast<int>(allsp.size());
00355 spnames.resize(nsp);
00356 for (int nn = 0; nn < nsp; nn++) {
00357 string stemp = (*allsp[nn])["name"];
00358 bool skip = false;
00359 if (declared[stemp]) {
00360 skip = true;
00361 }
00362 if (!skip) {
00363 declared[stemp] = true;
00364 nSpecies++;
00365 spNamesList.resize(nSpecies);
00366 spDataNodeList.resize(nSpecies, 0);
00367 spRuleList.resize(nSpecies, 0);
00368 spNamesList[nSpecies-1] = stemp;
00369 spDataNodeList[nSpecies-1] = allsp[nn];
00370 spRuleList[nSpecies-1] = sprule[jsp];
00371 }
00372 }
00373 } else {
00374 for (int k = 0; k < nsp; k++) {
00375 string stemp = spnames[k];
00376 skip = false;
00377 if (declared[stemp]) {
00378 if (sprule[jsp] >= 10) {
00379 skip = true;
00380 } else {
00381 throw CanteraError("ThermoFactory::formSpeciesXMLNodeList()",
00382 "duplicate species: \"" + stemp + "\"");
00383 }
00384 }
00385 if (!skip) {
00386 declared[stemp] = true;
00387
00388 XML_Node* s = db->findByAttr("name", stemp);
00389 if (!s) {
00390 throw CanteraError("importPhase","no data for species, \""
00391 + stemp + "\"");
00392 }
00393 nSpecies++;
00394 spNamesList.resize(nSpecies);
00395 spDataNodeList.resize(nSpecies, 0);
00396 spRuleList.resize(nSpecies, 0);
00397 spNamesList[nSpecies-1] = stemp;
00398 spDataNodeList[nSpecies-1] = s;
00399 spRuleList[nSpecies-1] = sprule[jsp];
00400 }
00401 }
00402 }
00403 }
00404 }
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431 bool importPhase(XML_Node& phase, ThermoPhase* th,
00432 SpeciesThermoFactory* spfactory) {
00433
00434
00435
00436 if (phase.name() != "phase") {
00437 throw CanteraError("importPhase",
00438 "Current const XML_Node named, " + phase.name() +
00439 ", is not a phase element.");
00440 }
00441
00442
00443
00444
00445
00446
00447
00448
00449 XML_Node &phaseNode_XML = th->xml();
00450 phaseNode_XML.clear();
00451 phase.copy(&phaseNode_XML);
00452
00453
00454
00455 th->setID(phase.id());
00456 th->setName(phase.id());
00457
00458
00459 if (phase.hasAttrib("dim")) {
00460 int idim = intValue(phase["dim"]);
00461 if (idim < 1 || idim > 3)
00462 throw CanteraError("importPhase",
00463 "phase, " + th->id() +
00464 ", has unphysical number of dimensions: " + phase["dim"]);
00465 th->setNDim(idim);
00466 }
00467 else {
00468 th->setNDim(3);
00469 }
00470
00471
00472
00473
00474 if (phase.hasChild("thermo")) {
00475 const XML_Node& eos = phase.child("thermo");
00476 th->setParametersFromXML(eos);
00477 } else {
00478 throw CanteraError("importPhase",
00479 " phase, " + th->id() +
00480 ", XML_Node does not have a \"thermo\" XML_Node");
00481 }
00482
00483 VPStandardStateTP *vpss_ptr = 0;
00484 int ssConvention = th->standardStateConvention();
00485 if (ssConvention == cSS_CONVENTION_VPSS) {
00486 vpss_ptr = dynamic_cast <VPStandardStateTP *>(th);
00487 if (vpss_ptr == 0) {
00488 throw CanteraError("importPhase",
00489 "phase, " + th->id() + ", was VPSS, but dynamic cast failed");
00490 }
00491 }
00492
00493
00494
00495 if (!spfactory) {
00496 spfactory = SpeciesThermoFactory::factory();
00497 }
00498
00499
00500
00501
00502 th->addElementsFromXML(phase);
00503
00504
00505
00506
00507
00508
00509
00510
00511 XML_Node* db = 0;
00512 vector<XML_Node*> sparrays;
00513 phase.getChildren("speciesArray", sparrays);
00514 int jsp, nspa = static_cast<int>(sparrays.size());
00515 if (nspa == 0) {
00516 throw CanteraError("importPhase",
00517 "phase, " + th->id() + ", has zero \"speciesArray\" XML nodes.\n"
00518 + " There must be at least one speciesArray nodes "
00519 "with one or more species");
00520 }
00521 vector<XML_Node*> dbases;
00522 vector_int sprule(nspa,0);
00523
00524
00525 for (jsp = 0; jsp < nspa; jsp++) {
00526
00527 const XML_Node& speciesArray = *sparrays[jsp];
00528
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541
00542
00543 if (speciesArray.hasChild("skip")) {
00544 const XML_Node& sk = speciesArray.child("skip");
00545 string eskip = sk["element"];
00546 if (eskip == "undeclared") {
00547 sprule[jsp] = 1;
00548 }
00549 string dskip = sk["species"];
00550 if (dskip == "duplicate") {
00551 sprule[jsp] += 10;
00552 }
00553 }
00554
00555 string fname, idstr;
00556
00557
00558
00559
00560
00561
00562 db = get_XML_Node(speciesArray["datasrc"], &phase.root());
00563 if (db == 0) {
00564 throw CanteraError("importPhase",
00565 " Can not find XML node for species database: "
00566 + speciesArray["datasrc"]);
00567 }
00568
00569
00570 dbases.push_back(db);
00571 }
00572
00573
00574
00575
00576
00577 std::vector<XML_Node *> spDataNodeList;
00578 std::vector<std::string> spNamesList;
00579 std::vector<int> spRuleList;
00580 formSpeciesXMLNodeList(spDataNodeList, spNamesList, spRuleList,
00581 sparrays, dbases, sprule);
00582
00583
00584
00585 delete &th->speciesThermo();
00586
00587
00588 SpeciesThermo* spth = 0;
00589 VPSSMgr* vp_spth = 0;
00590 if (ssConvention == cSS_CONVENTION_TEMPERATURE) {
00591
00592
00593
00594
00595
00596 spth = newSpeciesThermoMgr(spDataNodeList);
00597
00598
00599 th->setSpeciesThermo(spth);
00600 } else {
00601 vp_spth = newVPSSMgr(vpss_ptr, &phase, spDataNodeList);
00602 vpss_ptr->setVPSSMgr(vp_spth);
00603 spth = vp_spth->SpeciesThermoMgr();
00604 th->setSpeciesThermo(spth);
00605 }
00606
00607
00608 int k = 0;
00609
00610 int nsp = spDataNodeList.size();
00611 for (int i = 0; i < nsp; i++) {
00612 XML_Node *s = spDataNodeList[i];
00613 AssertTrace(s != 0);
00614 bool ok = installSpecies(k, *s, *th, spth, spRuleList[i],
00615 &phase, vp_spth, spfactory);
00616 if (ok) {
00617 th->saveSpeciesData(k, s);
00618 ++k;
00619 }
00620 }
00621
00622
00623 th->freezeSpecies();
00624
00625
00626 th->initThermo();
00627
00628
00629
00630 string id = "";
00631 th->initThermoXML(phase, id);
00632
00633 return true;
00634 }
00635
00636
00637
00638
00639
00640
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655
00656
00657
00658
00659
00660
00661
00662
00663
00664
00665
00666
00667
00668
00669
00670
00671
00672
00673
00674
00675
00676 bool installSpecies(int k, const XML_Node& s, thermo_t& th,
00677 SpeciesThermo *spthermo_ptr, int rule,
00678 XML_Node *phaseNode_ptr,
00679 VPSSMgr *vpss_ptr,
00680 SpeciesThermoFactory* factory) {
00681
00682 std::string xname = s.name();
00683 if (xname != "species") {
00684 throw CanteraError("installSpecies",
00685 "Unexpected XML name of species XML_Node: " + xname);
00686 }
00687
00688 const XML_Node& a = s.child("atomArray");
00689 map<string,string> comp;
00690 getMap(a, comp);
00691
00692
00693
00694
00695
00696 map<string,string>::const_iterator _b = comp.begin();
00697 for (; _b != comp.end(); ++_b) {
00698 if (th.elementIndex(_b->first) < 0) {
00699 if (rule == 0) {
00700 throw CanteraError("installSpecies",
00701 "Species " + s["name"] +
00702 " contains undeclared element " + _b->first);
00703 }
00704 else
00705 return false;
00706 }
00707 }
00708
00709
00710
00711
00712
00713 int m, nel = th.nElements();
00714 vector_fp ecomp(nel, 0.0);
00715 for (m = 0; m < nel; m++) {
00716 const char *es = comp[th.elementName(m)].c_str();
00717 if (strlen(es) > 0) {
00718 ecomp[m] = atofCheck(es);
00719 }
00720 }
00721
00722
00723
00724
00725
00726 doublereal chrg = 0.0;
00727 if (s.hasChild("charge")) chrg = getFloat(s, "charge");
00728
00729
00730
00731 doublereal sz = 1.0;
00732 if (s.hasChild("size")) sz = getFloat(s, "size");
00733
00734
00735 th.addUniqueSpecies(s["name"], &ecomp[0], chrg, sz);
00736
00737 if (vpss_ptr) {
00738 VPStandardStateTP *vp_ptr = dynamic_cast<VPStandardStateTP *>(&th);
00739 factory->installVPThermoForSpecies(k, s, vp_ptr, vpss_ptr, spthermo_ptr,
00740 phaseNode_ptr);
00741 } else {
00742
00743
00744 factory->installThermoForSpecies(k, s, &th, *spthermo_ptr, phaseNode_ptr);
00745 }
00746
00747
00748 return true;
00749 }
00750
00751
00752
00753
00754
00755
00756
00757
00758
00759
00760
00761
00762
00763
00764 const XML_Node *speciesXML_Node(std::string kname,
00765 const XML_Node *phaseSpeciesData) {
00766 if (!phaseSpeciesData) return ((const XML_Node *) 0);
00767 string jname = phaseSpeciesData->name();
00768 if (jname != "speciesData") {
00769 throw CanteraError("speciesXML_Node()",
00770 "Unexpected phaseSpeciesData name: " + jname);
00771 }
00772 vector<XML_Node*> xspecies;
00773 phaseSpeciesData->getChildren("species", xspecies);
00774 int jj = xspecies.size();
00775 for (int j = 0; j < jj; j++) {
00776 const XML_Node& sp = *xspecies[j];
00777 jname = sp["name"];
00778 if (jname == kname) {
00779 return &sp;
00780 }
00781 }
00782 return ((const XML_Node *) 0);
00783 }
00784
00785 }