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

Bug fixes related to deleting multispecies phases.

Browse source
This commit is contained in:
Harry Moffat 2008-07-25 20:12:05 +00:00
parent 3df90ebff6
commit e91bdbe60d
7 changed files with 281 additions and 207 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

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

@ -433,8 +433,12 @@ namespace Cantera {
}
m_moles[ip] = phasemoles;
if (nsp > 1) {
if (phasemoles > 0.0) {
p->setState_TPX(m_temp, m_press, n + loc);
p->getMoleFractions(DATA_PTR(m_moleFractions) + loc);
} else {
p->getMoleFractions(DATA_PTR(m_moleFractions) + loc);
}
}
else {
m_moleFractions[loc] = 1.0;

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

@ -887,177 +887,7 @@ namespace VCSnonideal {
fclose(FP);
}
/*!
* This function decides whether a phase has charged species
* or not.
*/
static bool hasChargedSpecies(Cantera::ThermoPhase *tPhase) {
int nSpPhase = tPhase->nSpecies();
for (int k = 0; k < nSpPhase; k++) {
if (tPhase->charge(k) != 0.0) {
return true;
}
}
return false;
}
/**********************************************************************
*
* chargeNeutralityElement():
*
* This utility routine decides whether a Cantera ThermoPhase needs
* a constraint equation representing the charge neutrality of the
* phase. It does this by searching for charged species. If it
* finds one, and if the phase needs one, then it returns true.
*/
static bool chargeNeutralityElement(Cantera::ThermoPhase *tPhase) {
int hasCharge = hasChargedSpecies(tPhase);
if (tPhase->chargeNeutralityNecessary()) {
if (hasCharge) {
return true;
}
}
return false;
}
/*************************************************************************
*
*
*
*/
static int setUpElements(vcs_VolPhase *VolPhase, Cantera::ThermoPhase *tPhase) {
int e, k, eT;
string ename;
int eFound = -2;
/*
*
*/
int nebase = tPhase->nElements();
int ne = nebase;
int ns = tPhase->nSpecies();
/*
* Decide whether we need an extra element constraint for charge
* neutrality of the phase
*/
bool cne = chargeNeutralityElement(tPhase);
if (cne) {
VolPhase->ChargeNeutralityElement = ne;
ne++;
}
/*
* Assign and malloc structures
*/
VolPhase->elemResize(ne);
if (VolPhase->ChargeNeutralityElement >= 0) {
VolPhase->m_elType[VolPhase->ChargeNeutralityElement] =
VCS_ELEM_TYPE_CHARGENEUTRALITY;
}
if (hasChargedSpecies(tPhase)) {
if (cne) {
/*
* We need a charge neutrality constraint.
* We also have an Electron Element. These are
* duplicates of each other. To avoid trouble with
* possible range error conflicts, sometimes we eliminate
* the Electron condition. Flag that condition for elimination
* by toggling the ElActive variable. If we find we need it
* later, we will retoggle ElActive to true.
*/
for (eT = 0; eT < nebase; eT++) {
ename = tPhase->elementName(eT);
if (ename == "E") {
eFound = eT;
VolPhase->ElActive[eT] = 0;
VolPhase->m_elType[eT] = VCS_ELEM_TYPE_ELECTRONCHARGE;
}
}
} else {
for (eT = 0; eT < nebase; eT++) {
ename = tPhase->elementName(eT);
if (ename == "E") {
eFound = eT;
VolPhase->m_elType[eT] = VCS_ELEM_TYPE_ELECTRONCHARGE;
}
}
}
if (eFound == -2) {
eFound = ne;
VolPhase->m_elType[ne] = VCS_ELEM_TYPE_ELECTRONCHARGE;
VolPhase->ElActive[ne] = 0;
string ename = "E";
VolPhase->ElName[ne] = ename;
ne++;
VolPhase->elemResize(ne);
}
}
VolPhase->FormulaMatrix.resize(ne, ns, 0.0);
VolPhase->m_speciesUnknownType.resize(ns, VCS_SPECIES_TYPE_MOLNUM);
VolPhase->elemResize(ne);
//ElGlobalIndex.resize(ne, -1);
e = 0;
for (eT = 0; eT < nebase; eT++) {
ename = tPhase->elementName(eT);
VolPhase->ElName[e] = ename;
e++;
}
if (cne) {
string pname = tPhase->id();
if (pname == "") {
char sss[50];
sprintf(sss, "phase%d", VolPhase->VP_ID);
pname = sss;
}
ename = "cn_" + pname;
e = VolPhase->ChargeNeutralityElement;
VolPhase->ElName[e] = ename;
}
double * const * const fm = VolPhase->FormulaMatrix.baseDataAddr();
for (k = 0; k < ns; k++) {
e = 0;
for (eT = 0; eT < nebase; eT++) {
fm[e][k] = tPhase->nAtoms(k, eT);
e++;
}
if (eFound >= 0) {
fm[eFound][k] = - tPhase->charge(k);
}
}
if (cne) {
for (k = 0; k < ns; k++) {
fm[VolPhase->ChargeNeutralityElement][k] = tPhase->charge(k);
}
}
/*
* Here, we figure out what is the species types are
* The logic isn't set in stone, and is just for a particular type
* of problem that I'm solving first.
*/
if (ns == 1) {
if (tPhase->charge(0) != 0.0) {
VolPhase->m_speciesUnknownType[0] = VCS_SPECIES_TYPE_INTERFACIALVOLTAGE;
VolPhase->setPhiVarIndex(0);
}
}
return ne;
}
static void print_char(const char letter, const int num) {
for (int i = 0; i < num; i++) plogf("%c", letter);
}
@ -1207,7 +1037,7 @@ namespace VCSnonideal {
* element in the phase to enforce a charge neutrality
* constraint.
*/
setUpElements(VolPhase, tPhase);
VolPhase->transferElementsFM(tPhase);
/*
* Combine the element information in the vcs_VolPhase

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

@ -18,6 +18,7 @@
#include "ThermoPhase.h"
#include "mix_defs.h"
#include <string>
#include <cstdio>
#include <cstdlib>
@ -1322,7 +1323,8 @@ namespace VCSnonideal {
//! Returns the global index of the local element index for the phase
void vcs_VolPhase::setElemGlobalIndex(const int eLocal, const int eGlobal) {
DebugAssertThrowVCS(eLocal >= 0, "vcs_VolPhase::setElemGlobalIndex");
DebugAssertThrowVCS(eLocal < m_numElemConstraints, "vcs_VolPhase::setElemGlobalIndex");
DebugAssertThrowVCS(eLocal < m_numElemConstraints,
"vcs_VolPhase::setElemGlobalIndex");
m_elemGlobalIndex[eLocal] = eGlobal;
}
@ -1330,5 +1332,173 @@ namespace VCSnonideal {
return m_numElemConstraints;
}
std::string vcs_VolPhase::elementName(const int e) const {
return ElName[e];
}
/*!
* This function decides whether a phase has charged species
* or not.
*/
static bool hasChargedSpecies(const Cantera::ThermoPhase * const tPhase) {
int nSpPhase = tPhase->nSpecies();
for (int k = 0; k < nSpPhase; k++) {
if (tPhase->charge(k) != 0.0) {
return true;
}
}
return false;
}
/**********************************************************************
*
* chargeNeutralityElement():
*
* This utility routine decides whether a Cantera ThermoPhase needs
* a constraint equation representing the charge neutrality of the
* phase. It does this by searching for charged species. If it
* finds one, and if the phase needs one, then it returns true.
*/
static bool chargeNeutralityElement(const Cantera::ThermoPhase * const tPhase) {
int hasCharge = hasChargedSpecies(tPhase);
if (tPhase->chargeNeutralityNecessary()) {
if (hasCharge) {
return true;
}
}
return false;
}
int vcs_VolPhase::transferElementsFM(const Cantera::ThermoPhase * const tPhase) {
int e, k, eT;
std::string ename;
int eFound = -2;
/*
*
*/
int nebase = tPhase->nElements();
int ne = nebase;
int ns = tPhase->nSpecies();
/*
* Decide whether we need an extra element constraint for charge
* neutrality of the phase
*/
bool cne = chargeNeutralityElement(tPhase);
if (cne) {
ChargeNeutralityElement = ne;
ne++;
}
/*
* Assign and malloc structures
*/
elemResize(ne);
if (ChargeNeutralityElement >= 0) {
m_elType[ChargeNeutralityElement] =
VCS_ELEM_TYPE_CHARGENEUTRALITY;
}
if (hasChargedSpecies(tPhase)) {
if (cne) {
/*
* We need a charge neutrality constraint.
* We also have an Electron Element. These are
* duplicates of each other. To avoid trouble with
* possible range error conflicts, sometimes we eliminate
* the Electron condition. Flag that condition for elimination
* by toggling the ElActive variable. If we find we need it
* later, we will retoggle ElActive to true.
*/
for (eT = 0; eT < nebase; eT++) {
ename = tPhase->elementName(eT);
if (ename == "E") {
eFound = eT;
ElActive[eT] = 0;
m_elType[eT] = VCS_ELEM_TYPE_ELECTRONCHARGE;
}
}
} else {
for (eT = 0; eT < nebase; eT++) {
ename = tPhase->elementName(eT);
if (ename == "E") {
eFound = eT;
m_elType[eT] = VCS_ELEM_TYPE_ELECTRONCHARGE;
}
}
}
if (eFound == -2) {
eFound = ne;
m_elType[ne] = VCS_ELEM_TYPE_ELECTRONCHARGE;
ElActive[ne] = 0;
std::string ename = "E";
ElName[ne] = ename;
ne++;
elemResize(ne);
}
}
FormulaMatrix.resize(ne, ns, 0.0);
m_speciesUnknownType.resize(ns, VCS_SPECIES_TYPE_MOLNUM);
elemResize(ne);
//ElGlobalIndex.resize(ne, -1);
e = 0;
for (eT = 0; eT < nebase; eT++) {
ename = tPhase->elementName(eT);
ElName[e] = ename;
e++;
}
if (cne) {
std::string pname = tPhase->id();
if (pname == "") {
char sss[50];
sprintf(sss, "phase%d", VP_ID);
pname = sss;
}
ename = "cn_" + pname;
e = ChargeNeutralityElement;
ElName[e] = ename;
}
double * const * const fm = FormulaMatrix.baseDataAddr();
for (k = 0; k < ns; k++) {
e = 0;
for (eT = 0; eT < nebase; eT++) {
fm[e][k] = tPhase->nAtoms(k, eT);
e++;
}
if (eFound >= 0) {
fm[eFound][k] = - tPhase->charge(k);
}
}
if (cne) {
for (k = 0; k < ns; k++) {
fm[ChargeNeutralityElement][k] = tPhase->charge(k);
}
}
/*
* Here, we figure out what is the species types are
* The logic isn't set in stone, and is just for a particular type
* of problem that I'm solving first.
*/
if (ns == 1) {
if (tPhase->charge(0) != 0.0) {
m_speciesUnknownType[0] = VCS_SPECIES_TYPE_INTERFACIALVOLTAGE;
setPhiVarIndex(0);
}
}
return ne;
}
}

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

@ -480,6 +480,22 @@ namespace VCSnonideal {
int nElemConstraints() const;
std::string elementName(const int e) const;
//! Transfer all of the element information from the
//! ThermoPhase object to the vcs_VolPhase object.
/*!
* Also decide whether we need a new charge neutrality
* element in the phase to enforce a charge neutrality
* constraint.
*
* @param tPhase Pointer to the thermophase object
*/
int transferElementsFM(const Cantera::ThermoPhase * const tPhase);
private:
//! Evaluate the activity coefficients at the current conditions
@ -549,11 +565,12 @@ namespace VCSnonideal {
*/
void _updateMoleFractionDependencies();
/*************************************************************************
* MEMBER DATA *
************************************************************************/
public:
private:
//! Backtrack value of VCS_SOLVE *
/*!
* Note the default for this is 0. That's a valid value too, since
@ -561,6 +578,7 @@ namespace VCSnonideal {
*/
VCS_SOLVE *m_owningSolverObject;
public:
//! Original ID of the phase in the problem.
/*!
* If a non-ideal phase splits into two due to a
@ -610,12 +628,13 @@ namespace VCSnonideal {
*/
int ChargeNeutralityElement;
private:
//! vector of strings containing the element names
/*!
* Length = nElemConstraints
*/
std::vector<std::string> ElName;
public:
//! boolean indicating whether an element constraint is active
//! for the current problem
std::vector<int> ElActive;

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

@ -369,7 +369,7 @@ namespace VCSnonideal {
*/
for (eVP = 0; eVP < neVP; eVP++) {
foundPos = -1;
enVP = volPhase->ElName[eVP];
enVP = volPhase->elementName(eVP);
/*
* Search for matches with the existing elements.
* If found, then fill in the entry in the global

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

@ -1059,8 +1059,10 @@ private:
* loop.
*
* @param iph Phase to be deleted
*
* @return Returns whether the operation was successful or not
*/
void vcs_delete_multiphase(const int iph);
bool vcs_delete_multiphase(const int iph);
//! Change the concentration of a species by delta moles.
/*!

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

@ -2116,8 +2116,8 @@ namespace VCSnonideal {
*
*
* @param kspec The species index
* @delta_ptr pointer to the delta for the species. This may change during
* the calculation
* @param delta_ptr pointer to the delta for the species. This may change during
* the calculation
*
* @return
* 1: succeeded without change of dx
@ -2129,6 +2129,13 @@ namespace VCSnonideal {
int j;
double tmp;
double delta = *delta_ptr;
#ifdef DEBUG_MODE
if (irxn < 0) {
plogf(" --- delete_species() ERROR: called for a component %d", kspec);
plogendl();
std::exit(-1);
}
#endif
if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
/*
* Attempt the given dx. If it doesn't work, try to see if a smaller
@ -2406,10 +2413,11 @@ namespace VCSnonideal {
*
* @param iph Phase to be deleted
*/
void VCS_SOLVE::vcs_delete_multiphase(const int iph) {
int kspec, j, irxn;
bool VCS_SOLVE::vcs_delete_multiphase(const int iph) {
int kspec, irxn;
double dx;
vcs_VolPhase *Vphase = m_VolPhaseList[iph];
bool successful = true;
/*
* set the phase existence flag to dead
*/
@ -2426,6 +2434,8 @@ namespace VCSnonideal {
m_tPhaseMoles_new[iph] = 0.0;
m_deltaPhaseMoles[iph] = 0.0;
/*
* Loop over all of the active species in the phase.
*/
@ -2433,33 +2443,71 @@ namespace VCSnonideal {
if (m_phaseID[kspec] == iph) {
if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
irxn = kspec - m_numComponents;
/*
* calculate an extent of rxn, dx, that zeroes out the species.
*/
dx = - (m_molNumSpecies_old[kspec]);
/*
* Set the mole numbers of that species to zero.
*/
m_molNumSpecies_old[kspec] = 0.0;
m_molNumSpecies_new[kspec] = 0.0;
m_deltaMolNumSpecies[kspec] = 0.0;
/*
* Change the status flag of the species to that of an
* zeroed phase
*/
m_rxnStatus[irxn] = VCS_SPECIES_ZEROEDPHASE;
/*
* Changed the component mole numbers to account for the
* final extent of reaction. Make sure to keep component
* mole numbers above zero
*
*/
for (j = 0; j < m_numComponents; ++j) {
m_molNumSpecies_old[j] += m_stoichCoeffRxnMatrix[irxn][j] * dx;
if (m_speciesUnknownType[j] == VCS_SPECIES_TYPE_MOLNUM) {
if (m_molNumSpecies_old[j] < 0.0) {
m_molNumSpecies_old[j] = 0.0;
if (irxn >= 0) {
/*
* calculate an extent of rxn, dx, that zeroes out the species.
*/
dx = - (m_molNumSpecies_old[kspec]);
double dxTent = dx;
int retn = delta_species(kspec, &dxTent);
if (retn != 1) {
successful = false;
#ifdef DEBUG_MODE
if (m_debug_print_lvl >= 2) {
plogf(" --- delete_multiphase %d, %s ERROR problems deleting species %s\n",
iph, Vphase->PhaseName.c_str(), m_speciesName[kspec].c_str() );
plogf(" --- delta attempted: %g achieved: %g "
" Zeroing it manually\n", dx, dxTent);
}
#endif
m_molNumSpecies_old[kspec] = 0.0;
m_molNumSpecies_new[kspec] = 0.0;
m_deltaMolNumSpecies[kspec] = 0.0;
// recover the total phase moles.
vcs_tmoles();
} else {
/*
* Set the mole number of that species to zero.
*/
m_molNumSpecies_old[kspec] = 0.0;
m_molNumSpecies_new[kspec] = 0.0;
m_deltaMolNumSpecies[kspec] = 0.0;
}
/*
* Change the status flag of the species to that of an
* zeroed phase
*/
m_rxnStatus[irxn] = VCS_SPECIES_ZEROEDPHASE;
/*
* Changed the component mole numbers to account for the
* final extent of reaction. Make sure to keep component
* mole numbers above zero
*
*/
// for (j = 0; j < m_numComponents; ++j) {
// m_molNumSpecies_old[j] += m_stoichCoeffRxnMatrix[irxn][j] * dx;
//if (m_speciesUnknownType[j] == VCS_SPECIES_TYPE_MOLNUM) {
//if (m_molNumSpecies_old[j] < 0.0) {
// m_molNumSpecies_old[j] = 0.0;
//}
//}
//}
}
else {
#ifdef DEBUG_MODE
if (m_debug_print_lvl >= 2) {
plogf(" --- delete_multiphase One of the species is a component %d - %s with mole number %g\n",
kspec, m_speciesName[kspec].c_str(), m_molNumSpecies_old[kspec]);
}
#endif
if (m_molNumSpecies_old[kspec] > VCS_RELDELETE_SPECIES_CUTOFF * VCS_DELETE_PHASE_CUTOFF ) {
plogf(" --- delete_multiphase unknown situation error exit");
plogendl();
std::exit(-1);
} else {
m_molNumSpecies_old[kspec] = 0.0;
m_molNumSpecies_new[kspec] = 0.0;
}
}
}
@ -2505,6 +2553,7 @@ namespace VCSnonideal {
Vphase->setMolesFromVCSCheck(VCS_STATECALC_OLD,
VCS_DATA_PTR(m_molNumSpecies_old),
VCS_DATA_PTR(m_tPhaseMoles_old));
return successful;
}
/**********************************************************************************/
@ -4647,7 +4696,7 @@ namespace VCSnonideal {
if (! vcs_doubleEqual(tlogMoles[iph], tPhMoles_ptr[iph])) {
plogf("phase Moles may be off, iph = %d, %20.14g %20.14g \n",
iph, tlogMoles[iph], tPhMoles_ptr[iph]);
exit(0);
std::exit(0);
}
}
#endif