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

Changed variable names

Browse source
This commit is contained in:
Harry Moffat 2008-04-23 16:08:16 +00:00
parent 4769e38ac0
commit 0f4c0a3f7c
11 changed files with 703 additions and 682 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

74
Cantera/src/equil/vcs_Gibbs.cpp
View file

@ -19,12 +19,12 @@
namespace VCSnonideal {
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
double VCS_SOLVE::vcs_Total_Gibbs(double *molesSp, double *chemPot,
double *tPhMoles)
double VCS_SOLVE::vcs_Total_Gibbs(double *molesSp, double *chemPot,
double *tPhMoles)
/*************************************************************************
*
@ -36,14 +36,14 @@ double VCS_SOLVE::vcs_Total_Gibbs(double *molesSp, double *chemPot,
* Note, for this algorithm this function should be MONOTONICALLY
* DECREASING.
*************************************************************************/
{
{
double g = 0.0;
for (int iph = 0; iph < NPhase; iph++) {
vcs_VolPhase *Vphase = VPhaseList[iph];
if ((TPhInertMoles[iph] > 0.0) && (tPhMoles[iph] > 0.0)) {
g += TPhInertMoles[iph] *
log(TPhInertMoles[iph] / tPhMoles[iph]);
log(TPhInertMoles[iph] / tPhMoles[iph]);
if (Vphase->GasPhase) {
g += TPhInertMoles[iph] * log(Pres);
}
@ -55,48 +55,40 @@ double VCS_SOLVE::vcs_Total_Gibbs(double *molesSp, double *chemPot,
}
return g;
}
}
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
double VCS_SOLVE::vcs_GibbsPhase(int iphase, double *w, double *fe)
/*************************************************************************
*
* vcs_Total_Gibbs:
*
* Calculate the total dimensionless Gibbs free energy
* -> Inert species are handled as if they had a standard free
* energy of zero.
* Note, for this algorithm this function should be MONOTONICALLY
* DECREASING.
*************************************************************************/
{
// Calculate the total dimensionless Gibbs free energy of a single phase
/*
* -> Inert species are handled as if they had a standard free
* energy of zero and if they obeyed ideal solution/gas theory
*
* @param iphase ID of the phase
* @param w Species mole number vector
* @param fe vector of partial molar free energies of the species.
*/
double VCS_SOLVE::vcs_GibbsPhase(int iphase, const double * const w,
const double * const fe) {
double g = 0.0;
vcs_VolPhase *Vphase = VPhaseList[iphase];
if ((TPhInertMoles[iphase] > 0.0) && (TPhMoles[iphase] > 0.0)) {
g += TPhInertMoles[iphase] *
log(TPhInertMoles[iphase] / TPhMoles[iphase]);
double phaseMols = 0.0;
for (int kspec = 0; kspec < m_numSpeciesRdc; ++kspec) {
if (PhaseID[kspec] == iphase) {
g += w[kspec] * fe[kspec];
phaseMols += w[kspec];
}
}
if (TPhInertMoles[iphase] > 0.0) {
phaseMols += TPhInertMoles[iphase];
g += TPhInertMoles[iphase] * log(TPhInertMoles[iphase] / phaseMols);
vcs_VolPhase *Vphase = VPhaseList[iphase];
if (Vphase->GasPhase == iphase) {
g += TPhInertMoles[iphase] * log(Pres);
}
}
for (int kspec = 0; kspec < m_numSpeciesRdc; ++kspec) {
if (PhaseID[kspec] == iphase) {
g += w[kspec] * fe[kspec];
}
}
return g;
}
}
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
}

1043
Cantera/src/equil/vcs_elem.cpp
View file

File diff suppressed because it is too large Load diff

1
Cantera/src/equil/vcs_elem_rearrange.cpp
View file

@ -7,7 +7,6 @@
* $Date$
* $Revision$
*/
/*
* Copywrite (2005) Sandia Corporation. Under the terms of
* Contract DE-AC04-94AL85000 with Sandia Corporation, the

67
Cantera/src/equil/vcs_inest.cpp
View file

@ -1,6 +1,6 @@
/**
* @file vcs_inest.cpp
* Methods for obtaining a good initial guess
* Implementation methods for obtaining a good initial guess
*/
/* $Author$
* $Date$
@ -13,36 +13,35 @@
* U.S. Government retains certain rights in this software.
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "vcs_solve.h"
#include "vcs_internal.h"
#include "vcs_VolPhase.h"
#include "clockWC.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
namespace VCSnonideal {
static char pprefix[20] = " --- vcs_inest: ";
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
void VCS_SOLVE::inest(double *aw, double *sa, double *sm,
double *ss, double test)
/**************************************************************************
*
* inest:
*
* Estimates equilibrium compositions.
* Algorithm covered in a section of Smith and Missen's Book.
*
* Linear programming module is based on using dbolm.
***************************************************************************/
{
// Estimate equilibrium compositions
/*
* Estimates equilibrium compositions.
* Algorithm covered in a section of Smith and Missen's Book.
*
* Linear programming module is based on using dbolm.
*
* @param aw aw[i[ Mole fraction work space (ne in length)
* @param sa sa[j] = Gramm-Schmidt orthog work space (ne in length)
* @param sm sm[i+j*ne] = QR matrix work space (ne*ne in length)
* @param ss ss[j] = Gramm-Schmidt orthog work space (ne in length)
* @param test This is a small negative number.
*/
void VCS_SOLVE::inest(double * const aw, double * const sa, double * const sm,
double * const ss, double test) {
int conv, k, lt, ikl, kspec, iph, irxn;
double s;
double s1 = 0.0;
@ -169,24 +168,24 @@ namespace VCSnonideal {
/* **** CHEMICAL POTENTIALS OF BASIS ****************** */
/* ***************************************************************** */
/*
* Calculate TMoles and TPhMoles[]
* Calculate TMoles and m_tPhaseMoles_old[]
*/
vcs_tmoles();
/*
* TPhMoles1[] will consist of just the component moles
* m_tPhaseMoles_new[] will consist of just the component moles
*/
for (iph = 0; iph < NPhase; iph++) {
TPhMoles1[iph] = TPhInertMoles[iph] + 1.0E-20;
m_tPhaseMoles_new[iph] = TPhInertMoles[iph] + 1.0E-20;
}
for (kspec = 0; kspec < m_numComponents; ++kspec) {
if (SpeciesUnknownType[kspec] == VCS_SPECIES_TYPE_MOLNUM) {
TPhMoles1[PhaseID[kspec]] += molNum[kspec];
m_tPhaseMoles_new[PhaseID[kspec]] += molNum[kspec];
}
}
TMolesMultiphase = 0.0;
for (iph = 0; iph < NPhase; iph++) {
if (! VPhaseList[iph]->SingleSpecies) {
TMolesMultiphase += TPhMoles1[iph];
TMolesMultiphase += m_tPhaseMoles_new[iph];
}
}
vcs_dcopy(VCS_DATA_PTR(m_molNumSpecies_new), molNum, nspecies);
@ -202,7 +201,7 @@ namespace VCSnonideal {
if (SpeciesUnknownType[kspec] == VCS_SPECIES_TYPE_MOLNUM) {
if (! SSPhase[kspec]) {
iph = PhaseID[kspec];
m_feSpecies_curr[kspec] += log(m_molNumSpecies_new[kspec] / TPhMoles[iph]);
m_feSpecies_curr[kspec] += log(m_molNumSpecies_new[kspec] / m_tPhaseMoles_old[iph]);
}
} else {
m_molNumSpecies_new[kspec] = 0.0;
@ -227,8 +226,8 @@ namespace VCSnonideal {
/* ********************************************************** */
vcs_dzero(VCS_DATA_PTR(DelTPhMoles), NPhase);
for (iph = 0; iph < NPhase; iph++) {
xtphMax[iph] = log(TPhMoles1[iph] * 1.0E32);
xtphMin[iph] = log(TPhMoles1[iph] * 1.0E-32);
xtphMax[iph] = log(m_tPhaseMoles_new[iph] * 1.0E32);
xtphMin[iph] = log(m_tPhaseMoles_new[iph] * 1.0E-32);
}
for (irxn = 0; irxn < nrxn; ++irxn) {
kspec = ir[irxn];
@ -251,7 +250,7 @@ namespace VCSnonideal {
* phase.
* It cut diamond4.vin iterations down from 62 to 14.
*/
m_deltaMolNumSpecies[kspec] = 0.5 * (TPhMoles1[iph] + TMolesMultiphase)
m_deltaMolNumSpecies[kspec] = 0.5 * (m_tPhaseMoles_new[iph] + TMolesMultiphase)
* exp(-m_deltaGRxn_new[irxn]);
for (k = 0; k < m_numComponents; ++k) {
@ -317,7 +316,7 @@ namespace VCSnonideal {
}
/*
* We have a new w[] estimate, go get the
* TMoles and TPhMoles[] values
* TMoles and m_tPhaseMoles_old[] values
*/
vcs_tmoles();
if (lt > 0) goto finished;
@ -505,7 +504,7 @@ namespace VCSnonideal {
if (vcs_debug_print_lvl >= 2) {
plogf("%sTotal Dimensionless Gibbs Free Energy = %15.7E", pprefix,
vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_feSpecies_curr),
VCS_DATA_PTR(TPhMoles)));
VCS_DATA_PTR(m_tPhaseMoles_old)));
plogendl();
}
#endif
@ -517,7 +516,7 @@ namespace VCSnonideal {
m_VCount->T_Time_inest += tsecond;
(m_VCount->T_Calls_Inest)++;
return retn;
}/**** vcs_inest() ***********************************************************/
}
}

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

@ -290,7 +290,7 @@ int VCS_SOLVE::vcs_prep(void) {
vcs_dzero(&(DnPhase[0][0]), m_numSpeciesTot*NPhase);
vcs_izero(&(PhaseParticipation[0][0]), m_numSpeciesTot*NPhase);
vcs_dzero(VCS_DATA_PTR(DelTPhMoles), NPhase);
vcs_dzero(VCS_DATA_PTR(TPhMoles1), NPhase);
vcs_dzero(VCS_DATA_PTR(m_tPhaseMoles_new), NPhase);
/*
* Calculate the total number of moles in all phases.
*/

14
Cantera/src/equil/vcs_report.cpp
View file

@ -153,7 +153,7 @@ int VCS_SOLVE::vcs_report(int iconv)
(VPhaseList[i])->PhaseName.c_str());
}
plogf("%14.7E %14.7E %12.4E\n", TPhInertMoles[i],
TPhInertMoles[i] / TPhMoles[i], 0.0);
TPhInertMoles[i] / m_tPhaseMoles_old[i], 0.0);
}
}
if (m_numSpeciesRdc != nspecies) {
@ -242,10 +242,10 @@ int VCS_SOLVE::vcs_report(int iconv)
plogf(" %3d ", iphase);
vcs_VolPhase *VPhase = VPhaseList[iphase];
plogf("%-12.12s |",VPhase->PhaseName.c_str());
plogf("%10.3e |", TPhMoles[iphase]);
totalMoles += TPhMoles[iphase];
if (TPhMoles[iphase] != VPhase->TotalMoles()) {
if (! vcs_doubleEqual(TPhMoles[iphase], VPhase->TotalMoles())) {
plogf("%10.3e |", m_tPhaseMoles_old[iphase]);
totalMoles += m_tPhaseMoles_old[iphase];
if (m_tPhaseMoles_old[iphase] != VPhase->TotalMoles()) {
if (! vcs_doubleEqual(m_tPhaseMoles_old[iphase], VPhase->TotalMoles())) {
plogf("We have a problem\n");
exit(-1);
}
@ -281,7 +281,7 @@ int VCS_SOLVE::vcs_report(int iconv)
*/
g = vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_feSpecies_curr),
VCS_DATA_PTR(TPhMoles));
VCS_DATA_PTR(m_tPhaseMoles_old));
plogf("\n\tTotal Dimensionless Gibbs Free Energy = G/RT = %15.7E\n", g);
if (inertYes)
plogf("\t\t(Inert species have standard free energy of zero)\n");
@ -315,7 +315,7 @@ int VCS_SOLVE::vcs_report(int iconv)
plogf(" %14.7E ", m_molNumSpecies_old[l]);
plogf("%14.7E ", m_SSfeSpecies[l]);
plogf("%14.7E ", log(ActCoeff[l]));
double tpmoles = TPhMoles[pid];
double tpmoles = m_tPhaseMoles_old[pid];
double phi = phasePhi[pid];
double eContrib = phi * Charge[l] * Faraday_dim;
double lx = 0.0;

10
Cantera/src/equil/vcs_rxnadj.cpp
View file

@ -142,8 +142,8 @@ int VCS_SOLVE::vcs_rxn_adj_cg(void)
}
for (j = 0; j < NPhase; j++) {
if (! (VPhaseList[j])->SingleSpecies) {
if (TPhMoles[j] > 0.0)
s -= SQUARE(dnPhase_irxn[j]) / TPhMoles[j];
if (m_tPhaseMoles_old[j] > 0.0)
s -= SQUARE(dnPhase_irxn[j]) / m_tPhaseMoles_old[j];
}
}
if (s != 0.0) {
@ -195,13 +195,13 @@ int VCS_SOLVE::vcs_rxn_adj_cg(void)
*/
if (dss != 0.0) {
m_molNumSpecies_old[kspec] += dss;
TPhMoles[PhaseID[kspec]] += dss;
m_tPhaseMoles_old[PhaseID[kspec]] += dss;
for (j = 0; j < m_numComponents; ++j) {
m_molNumSpecies_old[j] += dss * m_stoichCoeffRxnMatrix[irxn][j];
TPhMoles[PhaseID[j]] += dss * m_stoichCoeffRxnMatrix[irxn][j];
m_tPhaseMoles_old[PhaseID[j]] += dss * m_stoichCoeffRxnMatrix[irxn][j];
}
m_molNumSpecies_old[k] = 0.0;
TPhMoles[PhaseID[k]] = 0.0;
m_tPhaseMoles_old[PhaseID[k]] = 0.0;
#ifdef DEBUG_MODE
plogf(" --- vcs_st2 Special section to delete ");
plogf("%-12.12s", SpName[k].c_str());

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

@ -131,8 +131,8 @@ namespace VCSnonideal {
m_elemAbundances.resize(nelements, 0.0);
m_elemAbundancesGoal.resize(nelements, 0.0);
TPhMoles.resize(nphase0, 0.0);
TPhMoles1.resize(nphase0, 0.0);
m_tPhaseMoles_old.resize(nphase0, 0.0);
m_tPhaseMoles_new.resize(nphase0, 0.0);
DelTPhMoles.resize(nphase0, 0.0);
TmpPhase.resize(nphase0, 0.0);
TmpPhase2.resize(nphase0, 0.0);

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

@ -270,7 +270,13 @@ public:
void vcs_redim_TP(void);
void vcs_printChemPotUnits(int unitsFormat);
void vcs_elab(void);
//! Computes the current elemental abundances vector
/*!
* Computes the elemental abundances vector, m_elemAbundances[], and stores it
* back into the global structure
*/
void vcs_elab();
int vcs_elabcheck(int ibound);
void vcs_elabPhase(int iphase, double * const elemAbundPhase);
int vcs_elcorr(double aa[], double x[]);
@ -294,7 +300,19 @@ public:
#endif
double vcs_Total_Gibbs(double *w, double *fe, double *tPhMoles);
double vcs_GibbsPhase(int iphase, double *w, double *fe);
//! Calculate the total dimensionless Gibbs free energy of a single phase
/*!
* -> Inert species are handled as if they had a standard free
* energy of zero and if they obeyed ideal solution/gas theory
*
* @param iphase ID of the phase
* @param w Species mole number vector for all species
* @param fe vector of partial molar free energies of all of the
* species
*/
double vcs_GibbsPhase(int iphase, const double * const w,
const double * const fe);
double vcs_Gxs_phase_calc(vcs_VolPhase *Vphase, double *mf_PO);
double vcs_Gxs_calc(int iphase);
@ -358,8 +376,25 @@ private:
void prneav(void);
void checkDelta1(double * const ds, double * const delTPhMoles, int kspec);
#endif
void inest(double *aw, double *sa, double *sm,
double *ss, double test);
//! Estimate equilibrium compositions
/*!
* Estimates equilibrium compositions.
* Algorithm covered in a section of Smith and Missen's Book.
*
* Linear programming module is based on using dbolm.
*
* @param aw aw[i[ Mole fraction work space (ne in length)
* @param sa sa[j] = Gramm-Schmidt orthog work space (ne in length)
* @param sm sm[i+j*ne] = QR matrix work space (ne*ne in length)
* @param ss ss[j] = Gramm-Schmidt orthog work space (ne in length)
* @param test This is a small negative number.
*/
void inest(double * const aw, double * const sa, double * const sm,
double * const ss, double test);
void vcs_SSPhase(void);
double deltaG_Recalc_Rxn(int irxn, const double *const molNum,
double * const ac, double * const mu_i);
@ -612,7 +647,7 @@ public:
*
* Length = number of phases
*/
std::vector<double> TPhMoles;
std::vector<double> m_tPhaseMoles_old;
//! total gmols of species in each phase in the tentative soln vector
/*!
@ -621,7 +656,7 @@ public:
*
* Length = number of phases
*/
std::vector<double> TPhMoles1;
std::vector<double> m_tPhaseMoles_new;
//! Temporary vector of length NPhase
std::vector<double> TmpPhase;

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

@ -814,7 +814,7 @@ namespace VCSnonideal {
*/
if (m_molNumSpecies_new[kspec] < 0.005 * TMoles) {
iph = PhaseID[kspec];
if (m_molNumSpecies_new[kspec] < (TPhMoles[iph] * 0.01)) {
if (m_molNumSpecies_new[kspec] < (m_tPhaseMoles_old[iph] * 0.01)) {
#ifdef DEBUG_MODE
if (vcs_debug_print_lvl >= 2) {
plogf(" --- Major species changed to minor: ");
@ -887,9 +887,9 @@ namespace VCSnonideal {
*/
dnPhase_irxn = DnPhase[irxn];
for (int iphase = 0; iphase < NPhase; iphase++) {
TPhMoles[iphase] += dnPhase_irxn[iphase] * dx;
m_tPhaseMoles_old[iphase] += dnPhase_irxn[iphase] * dx;
}
TPhMoles[iph] = 0.0;
m_tPhaseMoles_old[iph] = 0.0;
vcs_updateVP(0);
/*
* Recalcuate the chemical potentials, FE(), and the
@ -1086,7 +1086,7 @@ namespace VCSnonideal {
* Calculate the tentative total mole numbers for each phase
*/
for (iph = 0; iph < NPhase; iph++) {
TPhMoles1[iph] = TPhMoles[iph] + DelTPhMoles[iph];
m_tPhaseMoles_new[iph] = m_tPhaseMoles_old[iph] + DelTPhMoles[iph];
}
/*
* Calculate the new chemical potentials using the tentative
@ -1119,10 +1119,10 @@ namespace VCSnonideal {
if (printDetails) {
plogf(" --- Total Old Dimensionless Gibbs Free Energy = %20.13E\n",
vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_feSpecies_old),
VCS_DATA_PTR(TPhMoles)));
VCS_DATA_PTR(m_tPhaseMoles_old)));
plogf(" --- Total tentative Dimensionless Gibbs Free Energy = %20.13E",
vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_new), VCS_DATA_PTR(m_feSpecies_curr),
VCS_DATA_PTR(TPhMoles1)));
VCS_DATA_PTR(m_tPhaseMoles_new)));
plogendl();
}
@ -1153,15 +1153,15 @@ namespace VCSnonideal {
plogf("Norms of Delta G():%14.6E%14.6E\n",
l2normdg(VCS_DATA_PTR(m_deltaGRxn_old)),
l2normdg(VCS_DATA_PTR(m_deltaGRxn_new)));
plogf(" Total moles of gas = %15.7E\n", TPhMoles[0]);
plogf(" Total moles of gas = %15.7E\n", m_tPhaseMoles_old[0]);
if ((NPhase > 1) && (! (VPhaseList[1])->SingleSpecies)) {
plogf(" Total moles of liquid = %15.7E\n", TPhMoles[1]);
plogf(" Total moles of liquid = %15.7E\n", m_tPhaseMoles_old[1]);
} else {
plogf(" Total moles of liquid = %15.7E\n", 0.0);
}
plogf(" Total New Dimensionless Gibbs Free Energy = %20.13E\n",
vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_new), VCS_DATA_PTR(m_feSpecies_curr),
VCS_DATA_PTR(TPhMoles1)));
VCS_DATA_PTR(m_tPhaseMoles_new)));
plogf(" -----------------------------------------------------");
plogendl();
}
@ -1211,15 +1211,15 @@ namespace VCSnonideal {
plogf(" --- "); vcs_print_line("-", 50);
for (iph = 0; iph < NPhase; iph++) {
Vphase = VPhaseList[iph];
plogf(" --- %18s = %15.7E\n", Vphase->PhaseName.c_str(), TPhMoles1[iph]);
plogf(" --- %18s = %15.7E\n", Vphase->PhaseName.c_str(), m_tPhaseMoles_new[iph]);
}
plogf(" "); vcs_print_line("-", 103);
plogf(" --- Total Old Dimensionless Gibbs Free Energy = %20.13E\n",
vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_feSpecies_old),
VCS_DATA_PTR(TPhMoles)));
VCS_DATA_PTR(m_tPhaseMoles_old)));
plogf(" --- Total New Dimensionless Gibbs Free Energy = %20.13E",
vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_new), VCS_DATA_PTR(m_feSpecies_curr),
VCS_DATA_PTR(TPhMoles1)));
VCS_DATA_PTR(m_tPhaseMoles_new)));
plogendl();
if (m_VCount->Its > 550) {
plogf(" --- Troublesome solve");
@ -1243,7 +1243,7 @@ namespace VCSnonideal {
* we have already done this inside the FORCED
* loop.
*/
vcs_dcopy(VCS_DATA_PTR(TPhMoles), VCS_DATA_PTR(TPhMoles1), NPhase);
vcs_dcopy(VCS_DATA_PTR(m_tPhaseMoles_old), VCS_DATA_PTR(m_tPhaseMoles_new), NPhase);
vcs_dcopy(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_molNumSpecies_new), m_numSpeciesRdc);
vcs_dcopy(VCS_DATA_PTR(m_deltaGRxn_old), VCS_DATA_PTR(m_deltaGRxn_new), m_numRxnRdc);
vcs_dcopy(VCS_DATA_PTR(m_feSpecies_old), VCS_DATA_PTR(m_feSpecies_curr), m_numSpeciesRdc);
@ -1275,8 +1275,8 @@ namespace VCSnonideal {
for (iph = 0; iph < NPhase; iph++) {
Vphase = VPhaseList[iph];
if (!(Vphase->SingleSpecies)) {
if (TPhMoles[iph] != 0.0 &&
TPhMoles[iph]/TMoles <= VCS_DELETE_PHASE_CUTOFF) {
if (m_tPhaseMoles_old[iph] != 0.0 &&
m_tPhaseMoles_old[iph]/TMoles <= VCS_DELETE_PHASE_CUTOFF) {
soldel = 1;
for (kspec = 0; kspec < m_numSpeciesRdc; kspec++) {
if (PhaseID[kspec] == iph && m_molNumSpecies_old[kspec] > 0.0) {
@ -1922,8 +1922,8 @@ namespace VCSnonideal {
m_molNumSpecies_new[kspec] = 1.0;
} else {
iph = PhaseID[kspec];
if (TPhMoles[iph] != 0.0) {
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] / TPhMoles[iph];
if (m_tPhaseMoles_old[iph] != 0.0) {
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] / m_tPhaseMoles_old[iph];
} else {
/*
* For MultiSpecies phases that are zeroed out,
@ -2166,12 +2166,12 @@ namespace VCSnonideal {
*delta_ptr = dx;
m_molNumSpecies_old[kspec] += dx;
int iph = PhaseID[kspec];
TPhMoles[iph] += dx;
m_tPhaseMoles_old[iph] += dx;
for (j = 0; j < m_numComponents; ++j) {
iph = PhaseID[j];
tmp = sc_irxn[j] * dx;
m_molNumSpecies_old[j] += tmp;
TPhMoles[iph] += tmp;
m_tPhaseMoles_old[iph] += tmp;
if (m_molNumSpecies_old[j] < 0.0) {
m_molNumSpecies_old[j] = 0.0;
}
@ -2193,7 +2193,7 @@ namespace VCSnonideal {
* Zero out the concentration of a species. Make sure to conserve
* elements and keep track of the total moles in all phases.
* w[]
* TPhMoles[]
* m_tPhaseMoles_old[]
*
* return:
* 1: succeeded
@ -2247,7 +2247,7 @@ namespace VCSnonideal {
int irxn = kspec - m_numComponents; /* This is the noncomponent rxn index */
/*
* Zero the concentration of the species.
* -> This zeroes w[kspec] and modifies TPhMoles[]
* -> This zeroes w[kspec] and modifies m_tPhaseMoles_old[]
*/
int retn = zero_species(kspec);
if (! retn) {
@ -2275,7 +2275,7 @@ namespace VCSnonideal {
/*
* Adjust the total moles in a phase downwards.
*/
Vphase->setMolesFromVCSCheck(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(TPhMoles));
Vphase->setMolesFromVCSCheck(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_tPhaseMoles_old));
/*
* Adjust the current number of active species and reactions counters
@ -2340,7 +2340,7 @@ namespace VCSnonideal {
#endif
/*
* Set the species back to minor species status
* this adjusts m_molNumSpecies_old[] and TPhMoles[]
* this adjusts m_molNumSpecies_old[] and m_tPhaseMoles_old[]
* HKM -> make this a relative mole number!
*/
dx = VCS_DELETE_SPECIES_CUTOFF * 10.;
@ -2353,7 +2353,7 @@ namespace VCSnonideal {
}
int iph = PhaseID[kspec];
vcs_VolPhase *Vphase = VPhaseList[iph];
Vphase->setMolesFromVCSCheck(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(TPhMoles));
Vphase->setMolesFromVCSCheck(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_tPhaseMoles_old));
/*
* We may have popped a multispecies phase back
* into existence. If we did, we have to check
@ -2421,8 +2421,8 @@ namespace VCSnonideal {
/*
* Zero out the total moles counters for the phase
*/
TPhMoles[iph] = 0.0;
TPhMoles1[iph] = 0.0;
m_tPhaseMoles_old[iph] = 0.0;
m_tPhaseMoles_new[iph] = 0.0;
DelTPhMoles[iph] = 0.0;
/*
@ -2499,7 +2499,7 @@ namespace VCSnonideal {
/*
* Upload the state to the VP object
*/
Vphase->setMolesFromVCSCheck(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(TPhMoles), iph);
Vphase->setMolesFromVCSCheck(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_tPhaseMoles_old), iph);
} /* delete_multiphase() *****************************************************/
@ -2536,8 +2536,8 @@ namespace VCSnonideal {
vcs_deltag(0, true);
for (iph = 0; iph < NPhase; iph++) {
if (TPhMoles[iph] > 0.0)
xtcutoff[iph] = log (TPhMoles[iph] / VCS_DELETE_SPECIES_CUTOFF);
if (m_tPhaseMoles_old[iph] > 0.0)
xtcutoff[iph] = log (m_tPhaseMoles_old[iph] / VCS_DELETE_SPECIES_CUTOFF);
else
xtcutoff[iph] = 0.0;
}
@ -2546,7 +2546,7 @@ namespace VCSnonideal {
* We are checking the equation:
*
* sum_u = sum_j_comp [ sigma_i_j * u_j ]
* = u_i_O + log((AC_i * W_i)/TPhMoles)
* = u_i_O + log((AC_i * W_i)/m_tPhaseMoles_old)
*
* by first evaluating:
*
@ -2571,14 +2571,14 @@ namespace VCSnonideal {
for (irxn = m_numRxnRdc; irxn < m_numRxnTot; ++irxn) {
kspec = ir[irxn];
iph = PhaseID[kspec];
if (TPhMoles[iph] == 0.0) {
if (m_tPhaseMoles_old[iph] == 0.0) {
if (m_deltaGRxn_new[irxn] < 0.0) {
vcs_reinsert_deleted(kspec);
npb++;
} else {
m_molNumSpecies_old[kspec] = 0.0;
}
} else if (TPhMoles[iph] > 0.0) {
} else if (m_tPhaseMoles_old[iph] > 0.0) {
if (m_deltaGRxn_new[irxn] < xtcutoff[iph]) {
vcs_reinsert_deleted(kspec);
npb++;
@ -2624,9 +2624,9 @@ namespace VCSnonideal {
for (int irxn = m_numRxnRdc; irxn < m_numRxnTot; ++irxn) {
kspec = ir[irxn];
iph = PhaseID[kspec];
if (TPhMoles[iph] > 0.0) {
if (m_tPhaseMoles_old[iph] > 0.0) {
double maxDG = MIN(m_deltaGRxn_new[irxn], 300);
double dx = TPhMoles[iph] * exp(- maxDG);
double dx = m_tPhaseMoles_old[iph] * exp(- maxDG);
retn = delta_species(kspec, &dx);
}
}
@ -2750,7 +2750,7 @@ namespace VCSnonideal {
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + al * m_deltaMolNumSpecies[kspec];
}
for (iph = 0; iph < NPhase; iph++) {
TPhMoles1[iph] = TPhMoles[iph] + al * DelTPhMoles[iph];
m_tPhaseMoles_new[iph] = m_tPhaseMoles_old[iph] + al * DelTPhMoles[iph];
}
vcs_updateVP(1);
@ -2865,7 +2865,7 @@ namespace VCSnonideal {
* is nontrivial in size.
*/
iph = PhaseID[kspec];
double tphmoles = TPhMoles[iph];
double tphmoles = m_tPhaseMoles_old[iph];
double trphmoles = tphmoles / TMoles;
if (trphmoles > VCS_DELETE_PHASE_CUTOFF) {
m_deltaMolNumSpecies[kspec] = TMoles * VCS_SMALL_MULTIPHASE_SPECIES;
@ -2946,8 +2946,8 @@ namespace VCSnonideal {
for (j = 0; j < NPhase; j++) {
Vphase = VPhaseList[j];
if (! Vphase->SingleSpecies) {
if (TPhMoles[j] > 0.0)
s -= SQUARE(dnPhase_irxn[j]) / TPhMoles[j];
if (m_tPhaseMoles_old[j] > 0.0)
s -= SQUARE(dnPhase_irxn[j]) / m_tPhaseMoles_old[j];
}
}
if (s != 0.0) {
@ -3053,16 +3053,16 @@ namespace VCSnonideal {
*/
if (dss != 0.0) {
m_molNumSpecies_old[kspec] += dss;
TPhMoles[PhaseID[kspec]] += dss;
m_tPhaseMoles_old[PhaseID[kspec]] += dss;
for (j = 0; j < m_numComponents; ++j) {
m_molNumSpecies_old[j] += dss * m_stoichCoeffRxnMatrix[irxn][j];
TPhMoles[PhaseID[j]] += dss * m_stoichCoeffRxnMatrix[irxn][j];
m_tPhaseMoles_old[PhaseID[j]] += dss * m_stoichCoeffRxnMatrix[irxn][j];
}
m_molNumSpecies_old[k] = 0.0;
iph = PhaseID[k];
Vphase = VPhaseList[iph];
Vphase->Existence = 0;
TPhMoles[iph] = 0.0;
m_tPhaseMoles_old[iph] = 0.0;
#ifdef DEBUG_MODE
if (vcs_debug_print_lvl >= 2) {
plogf(" --- vcs_RxnStepSizes Special section to delete %s\n",
@ -3901,7 +3901,7 @@ namespace VCSnonideal {
if (SSPhase[kspec]) {
return VCS_SPECIES_ZEROEDSS;
} else {
if (TPhMoles[iph] == 0.0) return VCS_SPECIES_ZEROEDPHASE;
if (m_tPhaseMoles_old[iph] == 0.0) return VCS_SPECIES_ZEROEDPHASE;
else return VCS_SPECIES_ZEROEDMS;
}
}
@ -3973,7 +3973,7 @@ namespace VCSnonideal {
* Check to see whether the current species is a major component
* of its phase. If it is, it is a major component
*/
if (m_molNumSpecies_old[kspec] > (TPhMoles[iph] * 0.1)) return VCS_SPECIES_MAJOR;
if (m_molNumSpecies_old[kspec] > (m_tPhaseMoles_old[iph] * 0.1)) return VCS_SPECIES_MAJOR;
/*
* Main check in the loop:
* Check to see if there is a component with a mole number that is
@ -4243,9 +4243,9 @@ namespace VCSnonideal {
}
#endif
if (kk <= 0) {
tPhMoles_ptr = VCS_DATA_PTR(TPhMoles);
tPhMoles_ptr = VCS_DATA_PTR(m_tPhaseMoles_old);
} else {
tPhMoles_ptr = VCS_DATA_PTR(TPhMoles1);
tPhMoles_ptr = VCS_DATA_PTR(m_tPhaseMoles_new);
}
tlogMoles = VCS_DATA_PTR(TmpPhase);
/*
@ -4535,20 +4535,20 @@ namespace VCSnonideal {
double sum;
vcs_VolPhase *Vphase;
for (i = 0; i < NPhase; i++) {
TPhMoles[i] = TPhInertMoles[i];
m_tPhaseMoles_old[i] = TPhInertMoles[i];
}
for (i = 0; i < m_numSpeciesTot; i++) {
if (SpeciesUnknownType[i] == VCS_SPECIES_TYPE_MOLNUM) {
TPhMoles[PhaseID[i]] += m_molNumSpecies_old[i];
m_tPhaseMoles_old[PhaseID[i]] += m_molNumSpecies_old[i];
}
}
sum = 0.0;
for (i = 0; i < NPhase; i++) {
sum += TPhMoles[i];
sum += m_tPhaseMoles_old[i];
Vphase = VPhaseList[i];
// Took out because we aren't updating mole fractions in Vphase
// Vphase->TMoles = TPhMoles[i];
if (TPhMoles[i] == 0.0) {
// Vphase->TMoles = m_tPhaseMoles_old[i];
if (m_tPhaseMoles_old[i] == 0.0) {
Vphase->Existence = 0;
} else {
if (TPhInertMoles[i] > 0.0) {
@ -4581,10 +4581,10 @@ namespace VCSnonideal {
Vphase = VPhaseList[i];
if (place == 0) {
Vphase->setMolesFromVCSCheck(VCS_DATA_PTR(m_molNumSpecies_old),
VCS_DATA_PTR(TPhMoles), i);
VCS_DATA_PTR(m_tPhaseMoles_old), i);
} else if (place == 1) {
Vphase->setMolesFromVCSCheck(VCS_DATA_PTR(m_molNumSpecies_new),
VCS_DATA_PTR(TPhMoles1), i);
VCS_DATA_PTR(m_tPhaseMoles_new), i);
} else {
plogf("we shouldn't be here\n");
exit(-1);

19
Cantera/src/equil/vcs_species_thermo.cpp
View file

@ -1,16 +1,17 @@
/**
* @file vcs_species_thermo.cpp
* Implementation for the VCS_SPECIES_THERMO object.
*/
/*
* $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 <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "vcs_solve.h"
#include "vcs_species_thermo.h"
@ -20,14 +21,14 @@
#include "vcs_Exception.h"
#include "vcs_internal.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
using namespace std;
namespace VCSnonideal {
/*****************************************************************************
*
* constructor():
*/
VCS_SPECIES_THERMO::VCS_SPECIES_THERMO(int indexPhase,
int indexSpeciesPhase) :
@ -436,7 +437,7 @@ double VCS_SOLVE::vcs_Gxs_calc(int iphase)
{
int kspec;
double Gxs = 0.0, ac;
double totmol = TPhMoles[iphase];
double totmol = m_tPhaseMoles_old[iphase];
vcs_VolPhase *Vphase = VPhaseList[iphase];
VCS_SPECIES_THERMO *ts_ptr;