research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Incremental update on solver.

Browse source 
Changing some names
This commit is contained in:
Harry Moffat 2010-07-30 18:16:31 +00:00
parent 14a349a14a
commit 43798deb62
8 changed files with 79 additions and 33 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

4
Cantera/src/equil/vcs_MultiPhaseEquil.cpp
View file

@ -1280,7 +1280,7 @@ namespace VCSnonideal {
vcs_VolPhase *VolPhase = vprob->VPhaseList[iphase];
std::string sEOS = string16_EOSType(VolPhase->m_eqnState);
plogf("%16s %5d %5d %8d %16s %8d %16e ", VolPhase->PhaseName.c_str(),
VolPhase->VP_ID, VolPhase->m_singleSpecies,
VolPhase->VP_ID_, VolPhase->m_singleSpecies,
VolPhase->m_gasPhase, sEOS.c_str(),
VolPhase->nSpecies(), VolPhase->totalMolesInert() );
plogf("%16e\n", VolPhase->totalMoles());
@ -1402,7 +1402,7 @@ namespace VCSnonideal {
vcs_VolPhase *VolPhase = vprob->VPhaseList[iphase];
std::string sEOS = string16_EOSType(VolPhase->m_eqnState);
plogf("%16s %5d %5d %8d %16s %8d %16e ", VolPhase->PhaseName.c_str(),
VolPhase->VP_ID, VolPhase->m_singleSpecies,
VolPhase->VP_ID_, VolPhase->m_singleSpecies,
VolPhase->m_gasPhase, sEOS.c_str(),
VolPhase->nSpecies(), VolPhase->totalMolesInert() );
plogf("%16e\n", VolPhase->totalMoles() );

18
Cantera/src/equil/vcs_VolPhase.cpp
View file

@ -32,7 +32,7 @@ namespace VCSnonideal {
*/
vcs_VolPhase::vcs_VolPhase(VCS_SOLVE * owningSolverObject) :
m_owningSolverObject(0),
VP_ID(-1),
VP_ID_(-1),
Domain_ID(-1),
m_singleSpecies(true),
m_gasPhase(false),
@ -94,7 +94,7 @@ namespace VCSnonideal {
*/
vcs_VolPhase::vcs_VolPhase(const vcs_VolPhase& b) :
m_owningSolverObject(b.m_owningSolverObject),
VP_ID(b.VP_ID),
VP_ID_(b.VP_ID_),
Domain_ID(b.Domain_ID),
m_singleSpecies(b.m_singleSpecies),
m_gasPhase(b.m_gasPhase),
@ -147,7 +147,7 @@ namespace VCSnonideal {
// operator but is true for a copy constructor
// m_owningSolverObject = b.m_owningSolverObject;
VP_ID = b.VP_ID;
VP_ID_ = b.VP_ID_;
Domain_ID = b.Domain_ID;
m_singleSpecies = b.m_singleSpecies;
m_gasPhase = b.m_gasPhase;
@ -261,17 +261,17 @@ namespace VCSnonideal {
m_phi = 0.0;
m_phiVarIndex = -1;
if (phaseNum == VP_ID) {
if (phaseNum == VP_ID_) {
if (strcmp(PhaseName.c_str(), phaseName)) {
plogf("Strings are different: %s %s :unknown situation\n",
PhaseName.c_str(), phaseName);
exit(EXIT_FAILURE);
}
} else {
VP_ID = phaseNum;
VP_ID_ = phaseNum;
if (!phaseName) {
char itmp[40];
sprintf(itmp, "Phase_%d", VP_ID);
sprintf(itmp, "Phase_%d", VP_ID_);
PhaseName = itmp;
} else {
PhaseName = phaseName;
@ -711,7 +711,7 @@ namespace VCSnonideal {
/*
* Check for consistency with TPhMoles[]
*/
double Tcheck = TPhMoles[VP_ID];
double Tcheck = TPhMoles[VP_ID_];
if (Tcheck != v_totalMoles) {
if (vcs_doubleEqual(Tcheck, v_totalMoles)) {
Tcheck = v_totalMoles;
@ -1104,7 +1104,7 @@ namespace VCSnonideal {
if (m_numSpecies != 0) {
plogf("Warning Nsp != NVolSpeces: %d %d \n", nsp, m_numSpecies);
}
resize(VP_ID, nsp, nelem, PhaseName.c_str());
resize(VP_ID_, nsp, nelem, PhaseName.c_str());
}
TP_ptr->getMoleFractions(VCS_DATA_PTR(Xmol));
fractionCreationDelta_ = Xmol;
@ -1567,7 +1567,7 @@ namespace VCSnonideal {
std::string pname = tPhase->id();
if (pname == "") {
char sss[50];
sprintf(sss, "phase%d", VP_ID);
sprintf(sss, "phase%d", VP_ID_);
pname = sss;
}
ename = "cn_" + pname;

2
Cantera/src/equil/vcs_VolPhase.h
View file

@ -673,7 +673,7 @@ namespace VCSnonideal {
* miscibility gap, these numbers will stay the
* same after the split.
*/
int VP_ID;
int VP_ID_;
//! ID of the surface or volume domain in which the
//! this phase exists

69
Cantera/src/equil/vcs_phaseStability.cpp
View file

@ -1,8 +1,17 @@
/* ======================================================================= */
/* $RCSfile: vcs_phaseStability.cpp,v $ */
/* $Date$ */
/* $Revision$ */
/* ======================================================================= */
/**
* @file vcs_phaseStability.cpp
* Implementation class for functions associated with determining the stability of a phase
* (see Class \link Cantera::VCS_SOLVE VCS_SOLVE\endlink and \ref equilfunctions ).
*/
/*
* $Id$
*/
/*
* Copywrite (2005) Sandia Corporation. Under the terms of
* Contract DE-AC04-94AL85000 with Sandia Corporation, the
* U.S. Government retains certain rights in this software.
*/
#include "vcs_solve.h"
@ -138,7 +147,7 @@ namespace VCSnonideal {
* @return returns the phase id of the phase that pops back into
* existence. Returns -1 if there are no phases
*/
int VCS_SOLVE::vcs_popPhaseID() {
int VCS_SOLVE::vcs_popPhaseID(std::vector<int> & phasePopPhaseIDs) {
int iphasePop = -1;
int iph;
int irxn, kspec;
@ -151,8 +160,8 @@ namespace VCSnonideal {
char anote[128];
if (m_debug_print_lvl >= 2) {
plogf(" --- vcs_popPhaseID() called\n");
plogf(" --- Phase Status F_e MoleNum\n");
plogf(" --------------------------------------------------------------\n");
plogf(" --- Phase Status F_e MoleNum\n");
plogf(" --------------------------------------------------------------------------\n");
}
#endif
for (iph = 0; iph < m_numPhases; iph++) {
@ -165,7 +174,7 @@ namespace VCSnonideal {
#ifdef DEBUG_MODE
if (m_debug_print_lvl >= 2) {
plogf(" --- %18s %5d NA %11.3e\n",
plogf(" --- %18s %5d NA %11.3e\n",
Vphase->PhaseName.c_str(),
existence,
m_tPhaseMoles_old[iph]);
@ -194,10 +203,19 @@ namespace VCSnonideal {
#endif
}
}
#ifdef DEBUG_MODE
if (Fephase < 0.0) {
strcpy(anote," (not stable)");
if (m_tPhaseMoles_old[iph] > 0.0) {
printf("shouldn't be here\n");
exit(-1);
}
}
#endif
#ifdef DEBUG_MODE
if (m_debug_print_lvl >= 2) {
plogf(" --- %18s %5d NA %11.3g %s\n",
plogf(" --- %18s %5d %10.3g %10.3g %s\n",
Vphase->PhaseName.c_str(),
existence, Fephase,
m_tPhaseMoles_old[iph], anote);
@ -242,11 +260,19 @@ namespace VCSnonideal {
}
}
}
phasePopPhaseIDs.resize(0);
if (iphasePop >= 0) {
phasePopPhaseIDs.push_back(iphasePop);
}
/*
* Insert logic here to figure out if phase pops are linked together. Only do one linked
* pop at a time.
*/
#ifdef DEBUG_MODE
if (m_debug_print_lvl >= 2) {
plogf(" --------------------------------------------------------------\n");
plogf(" ---------------------------------------------------------------------\n");
}
#endif
return iphasePop;
@ -280,6 +306,7 @@ namespace VCSnonideal {
int kspec = Vphase->spGlobalIndexVCS(0);
// Identify the formation reaction for that species
int irxn = kspec - m_numComponents;
std::vector<int> creationGlobalRxnNumbers;
doublereal s;
int j, k;
@ -378,8 +405,10 @@ namespace VCSnonideal {
} else {
vector<doublereal> fracDelta(Vphase->nSpecies());
vector<doublereal> X_est(Vphase->nSpecies());
// fracDelta = Vphase->creationMoleNumbers(creationGlobalRxnNumbers);
fracDelta = Vphase->fractionCreationDeltas();
double sumFrac = 0.0;
for (k = 0; k < Vphase->nSpecies(); k++) {
sumFrac += fracDelta[k];
@ -490,6 +519,8 @@ namespace VCSnonideal {
vector<doublereal> fracDelta_new(Vphase->nSpecies(), 0.0);
vector<doublereal> fracDelta_old(Vphase->nSpecies(), 0.0);
vector<doublereal> fracDelta_raw(Vphase->nSpecies(), 0.0);
vector<int> creationGlobalRxnNumbers(Vphase->nSpecies(), -1);
vector<doublereal> m_feSpecies_Deficient(m_numComponents, 0.0);
doublereal damp = 1.0;
@ -505,8 +536,8 @@ namespace VCSnonideal {
// Get the storred estimate for the composition of the phase if
// it gets created
// fracDelta_new = Vphase->creationMoleNumbers(creationGlobalRxnNumbers);
fracDelta_new = Vphase->fractionCreationDeltas();
bool oneIsComponent = false;
std::vector<int> componentList;
@ -539,6 +570,12 @@ namespace VCSnonideal {
plogf(" --- vcs_phaseStabilityTest() called for phase %d\n", iph);
}
#endif
for (k = 0; k < Vphase->nSpecies(); k++) {
if (fracDelta_new[k] < 1.0E-13) {
fracDelta_new[k] = 1.0E-13;
}
}
bool converged = false;
for (int its = 0; its < 200 && (!converged); its++) {
@ -547,6 +584,8 @@ namespace VCSnonideal {
fracDelta_old = fracDelta_new;
dirProdOld = dirProd;
// Given a set of fracDelta's, we calculate the fracDelta's
// for the component species, if any
for (i = 0; i < (int) componentList.size(); i++) {
@ -753,7 +792,9 @@ namespace VCSnonideal {
if (converged) {
Vphase->setMoleFractionsState(0.0, VCS_DATA_PTR(X_est),
VCS_STATECALC_PHASESTABILITY);
Vphase->setFractionCreationDeltas( VCS_DATA_PTR(fracDelta_new));
// Vphase->setCreationMoleNumbers(VCS_DATA_PTR(fracDelta_new), creationGlobalRxnNumbers);
Vphase->setFractionCreationDeltas(VCS_DATA_PTR(fracDelta_new));
}

2
Cantera/src/equil/vcs_prob.cpp
View file

@ -284,7 +284,7 @@ namespace VCSnonideal {
Vphase = VPhaseList[iphase];
std::string EOS_cstr = string16_EOSType(Vphase->m_eqnState);
plogf("%16s %5d %5d %8d ", Vphase->PhaseName.c_str(),
Vphase->VP_ID, Vphase->m_singleSpecies, Vphase->m_gasPhase);
Vphase->VP_ID_, Vphase->m_singleSpecies, Vphase->m_gasPhase);
plogf("%16s %8d %16e ", EOS_cstr.c_str(),
Vphase->nSpecies(), Vphase->totalMolesInert());
if (iest >= 0) plogf("%16e\n", Vphase->totalMoles());

4
Cantera/src/equil/vcs_solve.cpp
View file

@ -869,9 +869,9 @@ namespace VCSnonideal {
vcs_VolPhase *vPhase = m_VolPhaseList[iph];
vcs_VolPhase *pub_phase_ptr = pub->VPhaseList[iph];
if (vPhase->VP_ID != pub_phase_ptr->VP_ID) {
if (vPhase->VP_ID_ != pub_phase_ptr->VP_ID_) {
plogf("%sPhase numbers have changed:%d %d\n", yo.c_str(),
vPhase->VP_ID, pub_phase_ptr->VP_ID);
vPhase->VP_ID_, pub_phase_ptr->VP_ID_);
retn = VCS_PUB_BAD;
}

9
Cantera/src/equil/vcs_solve.h
View file

@ -1,6 +1,7 @@
/**
* @file vcs_solve.h
* Header file for the internal object that holds the problem
* Header file for the internal object that holds the vcs equilibrium problem
* (see Class \link Cantera::VCS_SOLVE VCS_SOLVE\endlink and \ref equilfunctions ).
*/
/*
* $Id$
@ -501,7 +502,6 @@ public:
*/
void vcs_updateVP(const int stateCalc);
//! Utility function that evaluates whether a phase can be popped
//! into existence
/*!
@ -514,10 +514,13 @@ public:
//! Decision as to whether a phase pops back into existence
/*!
* @param phasePopPhaseIDs Vector containing the phase ids of the phases
* that will be popped this step.
*
* @return returns the phase id of the phase that pops back into
* existence. Returns -1 if there are no phases
*/
int vcs_popPhaseID();
int vcs_popPhaseID(std::vector<int> &phasePopPhaseIDs);
//! Calculates the deltas of the reactions due to phases popping
//! into existence

4
Cantera/src/equil/vcs_solve_TP.cpp
View file

@ -122,6 +122,7 @@ namespace VCSnonideal {
double *dnPhase_irxn;
double atomComp;
int iphasePop;
std::vector<int> phasePopPhaseIDs(0);
#ifdef DEBUG_MODE
char ANOTE[128];
/*
@ -399,7 +400,8 @@ namespace VCSnonideal {
* First step is a major branch in the algorithm.
* We first determine if a phase pops into existence.
*/
iphasePop = vcs_popPhaseID();
phasePopPhaseIDs.clear();
iphasePop = vcs_popPhaseID(phasePopPhaseIDs);
/*
*
*/