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

Changed names of variables.

Browse source
This commit is contained in:
Harry Moffat 2008-05-16 21:18:29 +00:00
parent ca3bc7c5be
commit c76d8b7504
7 changed files with 342 additions and 336 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

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

@ -230,7 +230,7 @@ namespace VCSnonideal {
xtphMin[iph] = log(m_tPhaseMoles_new[iph] * 1.0E-32);
}
for (irxn = 0; irxn < nrxn; ++irxn) {
kspec = ir[irxn];
kspec = m_indexRxnToSpecies[irxn];
/*
* For single species phases, we will not estimate the
* mole numbers. If the phase exists, it stays. If it

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

@ -137,7 +137,7 @@ namespace VCSnonideal {
m_numRxnRdc = m_numRxnTot;
m_numSpeciesRdc = m_numSpeciesTot;
for (i = 0; i < m_numRxnRdc; ++i) {
ir[i] = m_numElemConstraints + i;
m_indexRxnToSpecies[i] = m_numElemConstraints + i;
}
for (kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
@ -226,7 +226,7 @@ namespace VCSnonideal {
if (m_numElemConstraints != m_numComponents) {
m_numRxnTot = m_numRxnRdc = m_numSpeciesTot - m_numComponents;
for (i = 0; i < m_numRxnRdc; ++i) {
ir[i] = m_numComponents + i;
m_indexRxnToSpecies[i] = m_numComponents + i;
}
}

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

@ -18,117 +18,117 @@
namespace VCSnonideal {
/*****************************************************************************/
static void print_space(int num)
{
for (int j = 0; j < num; j++) {
plogf(" ");
}
}
/*****************************************************************************/
static void print_space(int num)
{
for (int j = 0; j < num; j++) {
plogf(" ");
}
}
static void print_line(std::string schar, int num) {
static void print_line(std::string schar, int num) {
for (int j = 0; j < num; j++) plogf("%s", schar.c_str());
plogf("\n");
}
}
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
int VCS_SOLVE::vcs_report(int iconv)
int VCS_SOLVE::vcs_report(int iconv)
/**************************************************************************
*
* vcs_report:
*
* Print out a report on the state of the equilibrium problem to
* standard output.
* This prints out the current contents of the VCS_SOLVE class, V.
* The "old" solution vector is printed out.
***************************************************************************/
{
int i, j, l, k, inertYes = FALSE, kspec;
int nspecies = m_numSpeciesTot;
double g;
/**************************************************************************
*
* vcs_report:
*
* Print out a report on the state of the equilibrium problem to
* standard output.
* This prints out the current contents of the VCS_SOLVE class, V.
* The "old" solution vector is printed out.
***************************************************************************/
{
int i, j, l, k, inertYes = FALSE, kspec;
int nspecies = m_numSpeciesTot;
double g;
char originalUnitsState = UnitsState;
char originalUnitsState = UnitsState;
std::vector<int> sortindex(nspecies,0);
std::vector<double> xy(nspecies,0.0);
std::vector<int> sortindex(nspecies,0);
std::vector<double> xy(nspecies,0.0);
/* ************************************************************** */
/* **** SORT DEPENDENT SPECIES IN DECREASING ORDER OF MOLES ***** */
/* ************************************************************** */
/* ************************************************************** */
/* **** SORT DEPENDENT SPECIES IN DECREASING ORDER OF MOLES ***** */
/* ************************************************************** */
for (i = 0; i < nspecies; ++i) {
for (i = 0; i < nspecies; ++i) {
sortindex[i] = i;
xy[i] = m_molNumSpecies_old[i];
}
/*
* Sort the XY vector, the mole fraction vector,
* and the sort index vector, sortindex, according to
* the magnitude of the mole fraction vector.
*/
for (l = m_numComponents; l < m_numSpeciesRdc; ++l) {
}
/*
* Sort the XY vector, the mole fraction vector,
* and the sort index vector, sortindex, according to
* the magnitude of the mole fraction vector.
*/
for (l = m_numComponents; l < m_numSpeciesRdc; ++l) {
k = vcs_optMax(VCS_DATA_PTR(xy), 0, l, m_numSpeciesRdc);
if (k != l) {
vcsUtil_dsw(VCS_DATA_PTR(xy), k, l);
vcsUtil_isw(VCS_DATA_PTR(sortindex), k, l);
vcsUtil_dsw(VCS_DATA_PTR(xy), k, l);
vcsUtil_isw(VCS_DATA_PTR(sortindex), k, l);
}
}
}
/*
* Decide whether we have to nondimensionalize the equations.
* -> For the printouts from this routine, we will use nondimensional
* representations. This may be expanded in the future.
*/
if (UnitsState == VCS_DIMENSIONAL_G) {
/*
* Decide whether we have to nondimensionalize the equations.
* -> For the printouts from this routine, we will use nondimensional
* representations. This may be expanded in the future.
*/
if (UnitsState == VCS_DIMENSIONAL_G) {
vcs_nondim_TP();
}
}
/* ******************************************************** */
/* *** PRINT OUT RESULTS ********************************** */
/* ******************************************************** */
/* ******************************************************** */
/* *** PRINT OUT RESULTS ********************************** */
/* ******************************************************** */
plogf("\n\n\n\n");
print_line("-", 80);
print_line("-", 80);
plogf("\t\t VCS_TP REPORT\n");
print_line("-", 80);
print_line("-", 80);
if (iconv < 0) {
plogf("\n\n\n\n");
print_line("-", 80);
print_line("-", 80);
plogf("\t\t VCS_TP REPORT\n");
print_line("-", 80);
print_line("-", 80);
if (iconv < 0) {
plogf(" ERROR: CONVERGENCE CRITERION NOT SATISFIED.\n");
} else if (iconv == 1) {
plogf(" RANGE SPACE ERROR: Equilibrium Found but not all Element Abundances are Satisfied\n");
}
/*
* Calculate some quantities that may need updating
*/
vcs_tmoles();
Vol = vcs_VolTotal(m_temperature, m_pressurePA,
VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(VolPM));
} else if (iconv == 1) {
plogf(" RANGE SPACE ERROR: Equilibrium Found but not all Element Abundances are Satisfied\n");
}
/*
* Calculate some quantities that may need updating
*/
vcs_tmoles();
Vol = vcs_VolTotal(m_temperature, m_pressurePA,
VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(VolPM));
plogf("\t\tTemperature = %15.2g Kelvin\n", m_temperature);
plogf("\t\tPressure = %15.5g Pa \n", m_pressurePA);
plogf("\t\tVolume = %15.5g m**3\n", Vol);
plogf("\t\tTemperature = %15.2g Kelvin\n", m_temperature);
plogf("\t\tPressure = %15.5g Pa \n", m_pressurePA);
plogf("\t\tVolume = %15.5g m**3\n", Vol);
/*
* -------- TABLE OF SPECIES IN DECREASING MOLE NUMBERS --------------
*/
plogf("\n\n");
print_line("-", 80);
plogf(" Species Equilibrium kmoles ");
plogf("Mole Fraction ChemPot/RT SpecUnkType\n");
print_line("-", 80);
for (i = 0; i < m_numComponents; ++i) {
/*
* -------- TABLE OF SPECIES IN DECREASING MOLE NUMBERS --------------
*/
plogf("\n\n");
print_line("-", 80);
plogf(" Species Equilibrium kmoles ");
plogf("Mole Fraction ChemPot/RT SpecUnkType\n");
print_line("-", 80);
for (i = 0; i < m_numComponents; ++i) {
plogf(" %-12.12s", m_speciesName[i].c_str());
print_space(13);
plogf("%14.7E %14.7E %12.4E", m_molNumSpecies_old[i], m_molNumSpecies_new[i], m_feSpecies_curr[i]);
plogf(" %3d", m_speciesUnknownType[i]);
plogf("\n");
}
for (i = m_numComponents; i < m_numSpeciesRdc; ++i) {
}
for (i = m_numComponents; i < m_numSpeciesRdc; ++i) {
l = sortindex[i];
plogf(" %-12.12s", m_speciesName[l].c_str());
print_space(13);
@ -144,173 +144,175 @@ int VCS_SOLVE::vcs_report(int iconv)
exit(-1);
}
plogf("\n");
}
for (i = 0; i < m_numPhases; i++) {
}
for (i = 0; i < m_numPhases; i++) {
if (TPhInertMoles[i] > 0.0) {
inertYes = TRUE;
if (i == 0) {
plogf(" Inert Gas Species ");
} else {
plogf(" Inert Species in phase %16s ",
(VPhaseList[i])->PhaseName.c_str());
}
plogf("%14.7E %14.7E %12.4E\n", TPhInertMoles[i],
TPhInertMoles[i] / m_tPhaseMoles_old[i], 0.0);
inertYes = TRUE;
if (i == 0) {
plogf(" Inert Gas Species ");
} else {
plogf(" Inert Species in phase %16s ",
(VPhaseList[i])->PhaseName.c_str());
}
plogf("%14.7E %14.7E %12.4E\n", TPhInertMoles[i],
TPhInertMoles[i] / m_tPhaseMoles_old[i], 0.0);
}
}
if (m_numSpeciesRdc != nspecies) {
}
if (m_numSpeciesRdc != nspecies) {
plogf("\n SPECIES WITH LESS THAN 1.0E-32 KMOLES:\n\n");
for (kspec = m_numSpeciesRdc; kspec < nspecies; ++kspec) {
plogf(" %-12.12s", m_speciesName[kspec].c_str());
plogf(" %14.7E %14.7E %12.4E",
m_molNumSpecies_old[kspec], m_molNumSpecies_new[kspec], m_deltaGRxn_new[kspec]);
if (m_speciesUnknownType[i] == VCS_SPECIES_TYPE_MOLNUM) {
plogf(" KMol_Num");
} else if (m_speciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
plogf(" Voltage");
} else {
plogf(" Unknown");
}
plogf("\n");
}
}
print_line("-", 80);
plogf("\n");
/*
* ---------- TABLE OF SPECIES FORMATION REACTIONS ------------------
*/
plogf("\n");
print_line("-", m_numComponents*10 + 45);
plogf(" |ComponentID|");
for (j = 0; j < m_numComponents; j++) {
plogf(" %3d", j);
}
plogf(" | |\n");
plogf(" | Components|");
for (j = 0; j < m_numComponents; j++) {
plogf(" %10.10s", m_speciesName[j].c_str());
}
plogf(" | |\n");
plogf(" NonComponent | Moles |");
for (j = 0; j < m_numComponents; j++) {
plogf(" %10.3g", m_molNumSpecies_old[j]);
}
plogf(" | DG/RT Rxn |\n");
print_line("-", m_numComponents*10 + 45);
for (i = 0; i < m_numRxnTot; i++) {
int irxn = ir[i];
plogf(" %3d ", irxn);
plogf("%-10.10s", m_speciesName[irxn].c_str());
plogf("|%10.3g |", m_molNumSpecies_old[irxn]);
for (j = 0; j < m_numComponents; j++) {
plogf(" %6.2f", m_stoichCoeffRxnMatrix[i][j]);
plogf(" %-12.12s", m_speciesName[kspec].c_str());
// Note m_deltaGRxn_new[] stores in kspec slot not irxn slot, after solve
plogf(" %14.7E %14.7E %12.4E",
m_molNumSpecies_old[kspec], m_molNumSpecies_new[kspec], m_deltaGRxn_new[kspec]);
if (m_speciesUnknownType[i] == VCS_SPECIES_TYPE_MOLNUM) {
plogf(" KMol_Num");
} else if (m_speciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
plogf(" Voltage");
} else {
plogf(" Unknown");
}
plogf(" |%10.3g |", m_deltaGRxn_new[irxn]);
plogf("\n");
}
print_line("-", m_numComponents*10 + 45);
plogf("\n");
}
}
print_line("-", 80);
plogf("\n");
/*
* ------------------ TABLE OF PHASE INFORMATION ---------------------
*/
std::vector<double> gaPhase(m_numElemConstraints, 0.0);
std::vector<double> gaTPhase(m_numElemConstraints, 0.0);
double totalMoles = 0.0;
double gibbsPhase = 0.0;
double gibbsTotal = 0.0;
plogf("\n\n");
plogf("\n");
print_line("-", m_numElemConstraints*10 + 58);
plogf(" | ElementID |");
for (j = 0; j < m_numElemConstraints; j++) {
plogf(" %3d", j);
}
plogf(" | |\n");
plogf(" | Element |");
for (j = 0; j < m_numElemConstraints; j++) {
plogf(" %10.10s", (m_elementName[j]).c_str());
}
plogf(" | |\n");
plogf(" PhaseName |KMolTarget |");
for (j = 0; j < m_numElemConstraints; j++) {
plogf(" %10.3g", m_elemAbundancesGoal[j]);
}
plogf(" | Gibbs Total |\n");
print_line("-", m_numElemConstraints*10 + 58);
for (int iphase = 0; iphase < m_numPhases; iphase++) {
plogf(" %3d ", iphase);
vcs_VolPhase *VPhase = VPhaseList[iphase];
plogf("%-12.12s |",VPhase->PhaseName.c_str());
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);
}
}
vcs_elabPhase(iphase, VCS_DATA_PTR(gaPhase));
for (j = 0; j < m_numElemConstraints; j++) {
plogf(" %10.3g", gaPhase[j]);
gaTPhase[j] += gaPhase[j];
}
gibbsPhase = vcs_GibbsPhase(iphase, VCS_DATA_PTR(m_molNumSpecies_old),
VCS_DATA_PTR(m_feSpecies_curr));
gibbsTotal += gibbsPhase;
plogf(" | %18.11E |\n", gibbsPhase);
}
print_line("-", m_numElemConstraints*10 + 58);
plogf(" TOTAL |%10.3e |", totalMoles);
for (j = 0; j < m_numElemConstraints; j++) {
plogf(" %10.3g", gaTPhase[j]);
}
plogf(" | %18.11E |\n", gibbsTotal);
/*
* ---------- TABLE OF SPECIES FORMATION REACTIONS ------------------
*/
plogf("\n");
print_line("-", m_numComponents*10 + 45);
plogf(" |ComponentID|");
for (j = 0; j < m_numComponents; j++) {
plogf(" %3d", j);
}
plogf(" | |\n");
plogf(" | Components|");
for (j = 0; j < m_numComponents; j++) {
plogf(" %10.10s", m_speciesName[j].c_str());
}
plogf(" | |\n");
plogf(" NonComponent | Moles |");
for (j = 0; j < m_numComponents; j++) {
plogf(" %10.3g", m_molNumSpecies_old[j]);
}
plogf(" | DG/RT Rxn |\n");
print_line("-", m_numComponents*10 + 45);
for (i = 0; i < m_numRxnTot; i++) {
int kspec = m_indexRxnToSpecies[i];
plogf(" %3d ", kspec);
plogf("%-10.10s", m_speciesName[kspec].c_str());
plogf("|%10.3g |", m_molNumSpecies_old[kspec]);
for (j = 0; j < m_numComponents; j++) {
plogf(" %6.2f", m_stoichCoeffRxnMatrix[i][j]);
}
// Note m_deltaGRxn_new[] stores in kspec slot not irxn slot, after solve
plogf(" |%10.3g |", m_deltaGRxn_new[kspec]);
plogf("\n");
}
print_line("-", m_numComponents*10 + 45);
plogf("\n");
print_line("-", m_numElemConstraints*10 + 58);
plogf("\n");
/*
* ------------------ TABLE OF PHASE INFORMATION ---------------------
*/
std::vector<double> gaPhase(m_numElemConstraints, 0.0);
std::vector<double> gaTPhase(m_numElemConstraints, 0.0);
double totalMoles = 0.0;
double gibbsPhase = 0.0;
double gibbsTotal = 0.0;
plogf("\n\n");
plogf("\n");
print_line("-", m_numElemConstraints*10 + 58);
plogf(" | ElementID |");
for (j = 0; j < m_numElemConstraints; j++) {
plogf(" %3d", j);
}
plogf(" | |\n");
plogf(" | Element |");
for (j = 0; j < m_numElemConstraints; j++) {
plogf(" %10.10s", (m_elementName[j]).c_str());
}
plogf(" | |\n");
plogf(" PhaseName |KMolTarget |");
for (j = 0; j < m_numElemConstraints; j++) {
plogf(" %10.3g", m_elemAbundancesGoal[j]);
}
plogf(" | Gibbs Total |\n");
print_line("-", m_numElemConstraints*10 + 58);
for (int iphase = 0; iphase < m_numPhases; iphase++) {
plogf(" %3d ", iphase);
vcs_VolPhase *VPhase = VPhaseList[iphase];
plogf("%-12.12s |",VPhase->PhaseName.c_str());
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);
}
}
vcs_elabPhase(iphase, VCS_DATA_PTR(gaPhase));
for (j = 0; j < m_numElemConstraints; j++) {
plogf(" %10.3g", gaPhase[j]);
gaTPhase[j] += gaPhase[j];
}
gibbsPhase = vcs_GibbsPhase(iphase, VCS_DATA_PTR(m_molNumSpecies_old),
VCS_DATA_PTR(m_feSpecies_curr));
gibbsTotal += gibbsPhase;
plogf(" | %18.11E |\n", gibbsPhase);
}
print_line("-", m_numElemConstraints*10 + 58);
plogf(" TOTAL |%10.3e |", totalMoles);
for (j = 0; j < m_numElemConstraints; j++) {
plogf(" %10.3g", gaTPhase[j]);
}
plogf(" | %18.11E |\n", gibbsTotal);
/*
* ----------- GLOBAL SATISFACTION INFORMATION -----------------------
*/
print_line("-", m_numElemConstraints*10 + 58);
plogf("\n");
/*
* Calculate the total dimensionless Gibbs Free Energy
* -> Inert species are handled as if they had a standard free
* energy of zero
*/
/*
* ----------- GLOBAL SATISFACTION INFORMATION -----------------------
*/
/*
* Calculate the total dimensionless Gibbs Free Energy
* -> Inert species are handled as if they had a standard free
* energy of zero
*/
g = vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_feSpecies_curr),
VCS_DATA_PTR(m_tPhaseMoles_old));
plogf("\n\tTotal Dimensionless Gibbs Free Energy = G/RT = %15.7E\n", g);
if (inertYes)
g = vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_feSpecies_curr),
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");
plogf("\nElemental Abundances (kmol): ");
plogf(" Actual Target Type ElActive\n");
for (i = 0; i < m_numElemConstraints; ++i) {
plogf("\nElemental Abundances (kmol): ");
plogf(" Actual Target Type ElActive\n");
for (i = 0; i < m_numElemConstraints; ++i) {
print_space(26); plogf("%-2.2s", (m_elementName[i]).c_str());
plogf("%20.12E %20.12E", m_elemAbundances[i], m_elemAbundancesGoal[i]);
plogf(" %3d %3d\n", m_elType[i], m_elementActive[i]);
}
plogf("\n");
}
plogf("\n");
/*
* ------------------ TABLE OF SPECIES CHEM POTS ---------------------
*/
plogf("\n"); print_line("-", 93);
plogf("Chemical Potentials of the Species: (dimensionless)\n");
/*
* ------------------ TABLE OF SPECIES CHEM POTS ---------------------
*/
plogf("\n"); print_line("-", 93);
plogf("Chemical Potentials of the Species: (dimensionless)\n");
double rt = vcs_nondimMult_TP(m_VCS_UnitsFormat, m_temperature);
plogf("\t\t(RT = %g ", rt);
vcs_printChemPotUnits(m_VCS_UnitsFormat);
plogf(")\n");
plogf(" Name TKMoles StandStateChemPot "
double rt = vcs_nondimMult_TP(m_VCS_UnitsFormat, m_temperature);
plogf("\t\t(RT = %g ", rt);
vcs_printChemPotUnits(m_VCS_UnitsFormat);
plogf(")\n");
plogf(" Name TKMoles StandStateChemPot "
" ln(AC) ln(X_i) | F z_i phi | ChemPot | (-lnMnaught)\n");
print_line("-", 115);
for (i = 0; i < nspecies; ++i) {
print_line("-", 115);
for (i = 0; i < nspecies; ++i) {
l = sortindex[i];
int pid = m_phaseID[l];
plogf(" %-12.12s", m_speciesName[l].c_str());
@ -344,80 +346,80 @@ int VCS_SOLVE::vcs_report(int iconv)
plogf(" (%14.7E)", - SpecLnMnaught[l]);
}
plogf("\n");
}
print_line("-", 115);
}
print_line("-", 115);
/*
* ------------- TABLE OF SOLUTION COUNTERS --------------------------
*/
plogf("\n");
plogf("\nCounters: Iterations Time (seconds)\n");
if (m_timing_print_lvl > 0) {
plogf(" vcs_basopt: %5d %11.5E\n",
m_VCount->Basis_Opts, m_VCount->Time_basopt);
plogf(" vcs_TP: %5d %11.5E\n",
m_VCount->Its, m_VCount->Time_vcs_TP);
} else {
plogf(" vcs_basopt: %5d %11s\n",
m_VCount->Basis_Opts," NA ");
plogf(" vcs_TP: %5d %11s\n",
m_VCount->Its," NA " );
}
print_line("-", 80);
print_line("-", 80);
/*
* ------------- TABLE OF SOLUTION COUNTERS --------------------------
*/
plogf("\n");
plogf("\nCounters: Iterations Time (seconds)\n");
if (m_timing_print_lvl > 0) {
plogf(" vcs_basopt: %5d %11.5E\n",
m_VCount->Basis_Opts, m_VCount->Time_basopt);
plogf(" vcs_TP: %5d %11.5E\n",
m_VCount->Its, m_VCount->Time_vcs_TP);
} else {
plogf(" vcs_basopt: %5d %11s\n",
m_VCount->Basis_Opts," NA ");
plogf(" vcs_TP: %5d %11s\n",
m_VCount->Its," NA " );
}
print_line("-", 80);
print_line("-", 80);
/*
* Set the Units state of the system back to where it was when we
* entered the program.
*/
if (originalUnitsState != UnitsState) {
/*
* Set the Units state of the system back to where it was when we
* entered the program.
*/
if (originalUnitsState != UnitsState) {
if (originalUnitsState == VCS_DIMENSIONAL_G ) vcs_redim_TP();
else vcs_nondim_TP();
}
/*
* Return a successful completion flag
*/
return VCS_SUCCESS;
} /* vcs_report() ************************************************************/
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
}
/*
* Return a successful completion flag
*/
return VCS_SUCCESS;
} /* vcs_report() ************************************************************/
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
void VCS_SOLVE::vcs_TCounters_report(int timing_print_lvl)
void VCS_SOLVE::vcs_TCounters_report(int timing_print_lvl)
/**************************************************************************
*
* vcs_TCounters_report:
*
* Print out the total Its and time counters to standard output
***************************************************************************/
{
plogf("\nTCounters: Num_Calls Total_Its Total_Time (seconds)\n");
if (timing_print_lvl > 0) {
plogf(" vcs_basopt: %5d %5d %11.5E\n",
m_VCount->T_Basis_Opts, m_VCount->T_Basis_Opts,
m_VCount->T_Time_basopt);
plogf(" vcs_TP: %5d %5d %11.5E\n",
m_VCount->T_Calls_vcs_TP, m_VCount->T_Its,
m_VCount->T_Time_vcs_TP);
plogf(" vcs_inest: %5d %11.5E\n",
m_VCount->T_Calls_Inest, m_VCount->T_Time_inest);
plogf(" vcs_TotalTime: %11.5E\n",
m_VCount->T_Time_vcs);
} else {
plogf(" vcs_basopt: %5d %5d %11s\n",
m_VCount->T_Basis_Opts, m_VCount->T_Basis_Opts," NA ");
plogf(" vcs_TP: %5d %5d %11s\n",
m_VCount->T_Calls_vcs_TP, m_VCount->T_Its," NA ");
plogf(" vcs_inest: %5d %11s\n",
m_VCount->T_Calls_Inest, " NA ");
plogf(" vcs_TotalTime: %11s\n",
" NA ");
/**************************************************************************
*
* vcs_TCounters_report:
*
* Print out the total Its and time counters to standard output
***************************************************************************/
{
plogf("\nTCounters: Num_Calls Total_Its Total_Time (seconds)\n");
if (timing_print_lvl > 0) {
plogf(" vcs_basopt: %5d %5d %11.5E\n",
m_VCount->T_Basis_Opts, m_VCount->T_Basis_Opts,
m_VCount->T_Time_basopt);
plogf(" vcs_TP: %5d %5d %11.5E\n",
m_VCount->T_Calls_vcs_TP, m_VCount->T_Its,
m_VCount->T_Time_vcs_TP);
plogf(" vcs_inest: %5d %11.5E\n",
m_VCount->T_Calls_Inest, m_VCount->T_Time_inest);
plogf(" vcs_TotalTime: %11.5E\n",
m_VCount->T_Time_vcs);
} else {
plogf(" vcs_basopt: %5d %5d %11s\n",
m_VCount->T_Basis_Opts, m_VCount->T_Basis_Opts," NA ");
plogf(" vcs_TP: %5d %5d %11s\n",
m_VCount->T_Calls_vcs_TP, m_VCount->T_Its," NA ");
plogf(" vcs_inest: %5d %11s\n",
m_VCount->T_Calls_Inest, " NA ");
plogf(" vcs_TotalTime: %11s\n",
" NA ");
}
}
}
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
}

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

@ -70,7 +70,7 @@ namespace VCSnonideal {
sprintf(ANOTE,"Normal Calc");
#endif
kspec = ir[irxn];
kspec = m_indexRxnToSpecies[irxn];
dnPhase_irxn = m_deltaMolNumPhase[irxn];
if (m_molNumSpecies_old[kspec] == 0.0 && (! m_SSPhase[kspec])) {
@ -279,7 +279,7 @@ namespace VCSnonideal {
int kspec, k, l, kph;
double s;
double *sc_irxn;
kspec = ir[irxn];
kspec = m_indexRxnToSpecies[irxn];
kph = m_phaseID[kspec];
sc_irxn = m_stoichCoeffRxnMatrix[irxn];
/*
@ -413,7 +413,7 @@ namespace VCSnonideal {
{
int its = 0;
int k;
int kspec = ir[irxn];
int kspec = m_indexRxnToSpecies[irxn];
const int MAXITS = 10;
double dx = dx_orig;
double *sc_irxn = m_stoichCoeffRxnMatrix[irxn];

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

@ -158,7 +158,7 @@ namespace VCSnonideal {
* mapping. We can't fill it in until we know the number of c
* components in the problem
*/
ir.resize(nspecies0, 0);
m_indexRxnToSpecies.resize(nspecies0, 0);
/* Initialize all species to be major species */
m_rxnStatus.resize(nspecies0, 1);

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

@ -1004,7 +1004,7 @@ public:
//! Delta G(I) for the noncomponent species in the mechanism.
/*!
* Computed by the subroutine DELTAG. DG is the free
* Computed by the subroutine DELTAG. DG is the free
* energy change for the reaction which
* forms species K from the
* component species. This vector has length
@ -1157,16 +1157,18 @@ public:
//! Mapping between the species index for noncomponent species and the
//! full species index.
/*!
* ir[irxn] = Mapping between the species index for
* noncomponent species and the full species
* ir[irxn] = Mapping between the reaction index for
* noncomponent formation reaction of a species
* and the full species
* index.
* - Initially set to a value of K = NC + I
* This vector has length equal to number
* of noncomponent species in the mechanism.
* It starts with the first current
* noncomponent species in the mechanism.
* kspec = ir[irxn]
*/
std::vector<int> ir;
std::vector<int> m_indexRxnToSpecies;
//! Major -Minor status vector for the formation reaction
/*!

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

@ -315,7 +315,7 @@ namespace VCSnonideal {
/*************************************************************************/
m_numRxnMinorZeroed = 0;
for (irxn = 0; irxn < m_numRxnRdc; ++irxn) {
kspec = ir[irxn];
kspec = m_indexRxnToSpecies[irxn];
m_rxnStatus[irxn] = vcs_species_type(kspec);
if (m_rxnStatus[irxn] == VCS_SPECIES_MINOR) {
m_rxnStatus[irxn] = VCS_SPECIES_MAJOR;
@ -500,7 +500,7 @@ namespace VCSnonideal {
#endif
for (irxn = 0; irxn < m_numRxnRdc; irxn++) {
kspec = ir[irxn];
kspec = m_indexRxnToSpecies[irxn];
sc_irxn = m_stoichCoeffRxnMatrix[irxn];
iph = m_phaseID[kspec];
Vphase = VPhaseList[iph];
@ -1366,7 +1366,7 @@ namespace VCSnonideal {
dofast = false;
if (dofast) {
for (i = 0; i < m_numRxnRdc; ++i) {
l = ir[i];
l = m_indexRxnToSpecies[i];
for (j = m_numComponents - 1; j >= 0; j--) {
bool doSwap = false;
if (m_SSPhase[j]) {
@ -1426,7 +1426,7 @@ namespace VCSnonideal {
}
} else {
for (i = 0; i < m_numRxnRdc; ++i) {
l = ir[i];
l = m_indexRxnToSpecies[i];
for (j = 0; j < m_numComponents; ++j) {
bool doSwap = false;
if (m_SSPhase[j]) {
@ -1506,7 +1506,7 @@ namespace VCSnonideal {
#endif
m_numRxnMinorZeroed = 0;
for (irxn = 0; irxn < m_numRxnRdc; irxn++) {
kspec = ir[irxn];
kspec = m_indexRxnToSpecies[irxn];
int speciesType = vcs_species_type(kspec);
if (speciesType < VCS_SPECIES_MINOR) {
@ -1590,7 +1590,7 @@ namespace VCSnonideal {
} else {
#ifdef DEBUG_MODE
if (vcs_debug_print_lvl >= 2) {
plogf("%s failed\n", m_speciesName[ir[irxn]].c_str());
plogf("%s failed\n", m_speciesName[m_indexRxnToSpecies[irxn]].c_str());
}
#endif
/*
@ -1657,7 +1657,7 @@ namespace VCSnonideal {
}
#ifdef DEBUG_MODE
if (vcs_debug_print_lvl >= 2) {
plogf("%s failed\n", m_speciesName[ir[irxn]].c_str());
plogf("%s failed\n", m_speciesName[m_indexRxnToSpecies[irxn]].c_str());
}
#endif
/*
@ -2543,7 +2543,7 @@ namespace VCSnonideal {
*/
npb = 0;
for (irxn = m_numRxnRdc; irxn < m_numRxnTot; ++irxn) {
kspec = ir[irxn];
kspec = m_indexRxnToSpecies[irxn];
iph = m_phaseID[kspec];
if (m_tPhaseMoles_old[iph] == 0.0) {
if (m_deltaGRxn_new[irxn] < 0.0) {
@ -2596,7 +2596,7 @@ namespace VCSnonideal {
for (int irxn = m_numRxnRdc; irxn < m_numRxnTot; ++irxn) {
kspec = ir[irxn];
kspec = m_indexRxnToSpecies[irxn];
iph = m_phaseID[kspec];
if (m_tPhaseMoles_old[iph] > 0.0) {
double maxDG = MIN(m_deltaGRxn_new[irxn], 300);
@ -2818,7 +2818,7 @@ namespace VCSnonideal {
sprintf(ANOTE,"Normal Calc");
#endif
kspec = ir[irxn];
kspec = m_indexRxnToSpecies[irxn];
if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
@ -3130,7 +3130,7 @@ namespace VCSnonideal {
for (irxn = 0; irxn < m_numRxnRdc; ++irxn) {
if (m_rxnStatus[irxn] != VCS_SPECIES_MINOR) {
icase = 0;
m_deltaGRxn_new[irxn] = m_feSpecies_curr[ir[irxn]];
m_deltaGRxn_new[irxn] = m_feSpecies_curr[m_indexRxnToSpecies[irxn]];
dtmp_ptr = m_stoichCoeffRxnMatrix[irxn];
for (kspec = 0; kspec < m_numComponents; ++kspec) {
m_deltaGRxn_new[irxn] += dtmp_ptr[kspec] * m_feSpecies_curr[kspec];
@ -3149,11 +3149,12 @@ namespace VCSnonideal {
/* ************************************************* */
for (irxn = 0; irxn < irxnl; ++irxn) {
icase = 0;
m_deltaGRxn_new[irxn] = m_feSpecies_curr[ir[irxn]];
m_deltaGRxn_new[irxn] = m_feSpecies_curr[m_indexRxnToSpecies[irxn]];
dtmp_ptr = m_stoichCoeffRxnMatrix[irxn];
for (kspec = 0; kspec < m_numComponents; ++kspec) {
m_deltaGRxn_new[irxn] += dtmp_ptr[kspec] * m_feSpecies_curr[kspec];
if (m_molNumSpecies_old[kspec] < VCS_DELETE_MINORSPECIES_CUTOFF && dtmp_ptr[kspec] < 0.0) {
if (m_molNumSpecies_old[kspec] < VCS_DELETE_MINORSPECIES_CUTOFF &&
dtmp_ptr[kspec] < 0.0) {
icase = 1;
}
}
@ -3168,11 +3169,12 @@ namespace VCSnonideal {
for (irxn = 0; irxn < m_numRxnRdc; ++irxn) {
if (m_rxnStatus[irxn] <= VCS_SPECIES_MINOR) {
icase = 0;
m_deltaGRxn_new[irxn] = m_feSpecies_curr[ir[irxn]];
m_deltaGRxn_new[irxn] = m_feSpecies_curr[m_indexRxnToSpecies[irxn]];
dtmp_ptr = m_stoichCoeffRxnMatrix[irxn];
for (kspec = 0; kspec < m_numComponents; ++kspec) {
m_deltaGRxn_new[irxn] += dtmp_ptr[kspec] * m_feSpecies_curr[kspec];
if (m_molNumSpecies_old[kspec] < VCS_DELETE_MINORSPECIES_CUTOFF && dtmp_ptr[kspec] < 0.0) {
if (m_molNumSpecies_old[kspec] < VCS_DELETE_MINORSPECIES_CUTOFF &&
dtmp_ptr[kspec] < 0.0) {
icase = 1;
}
}
@ -3513,7 +3515,7 @@ namespace VCSnonideal {
m_numRxnRdc = m_numRxnTot - numPreDeleted;
m_numSpeciesRdc = m_numSpeciesTot - numPreDeleted;
for (i = 0; i < m_numSpeciesTot; ++i) {
ir[i] = ncTrial + i;
m_indexRxnToSpecies[i] = ncTrial + i;
}
#ifdef DEBUG_MODE
if (vcs_debug_print_lvl >= 2) {
@ -3649,7 +3651,7 @@ namespace VCSnonideal {
}
}
for (i = 0; i < m_numRxnTot; ++i) {
k = ir[i];
k = m_indexRxnToSpecies[i];
for (j = 0; j < ncTrial; ++j) {
m_stoichCoeffRxnMatrix[i][j] = m_formulaMatrix[j][k];
}
@ -3698,7 +3700,7 @@ namespace VCSnonideal {
}
}
for (i = 0; i < m_numRxnTot; ++i) {
k = ir[i];
k = m_indexRxnToSpecies[i];
for (j = 0; j < ncTrial; ++j) {
if (j == jlose) {
aw[j] = m_formulaMatrix[juse][k];
@ -3749,9 +3751,9 @@ namespace VCSnonideal {
//plogf("| m_scSize");
plogf("\n");
for (i = 0; i < m_numRxnTot; i++) {
plogf(" --- %3d ", ir[i]);
plogf("%-10.10s", m_speciesName[ir[i]].c_str());
plogf("|%10.3g|", m_molNumSpecies_old[ir[i]]);
plogf(" --- %3d ", m_indexRxnToSpecies[i]);
plogf("%-10.10s", m_speciesName[m_indexRxnToSpecies[i]].c_str());
plogf("|%10.3g|", m_molNumSpecies_old[m_indexRxnToSpecies[i]]);
for (j = 0; j < ncTrial; j++) {
plogf(" %6.2f", m_stoichCoeffRxnMatrix[i][j]);
}
@ -3782,7 +3784,7 @@ namespace VCSnonideal {
for (irxn = 0; irxn < m_numRxnTot; ++irxn) {
scrxn_ptr = m_stoichCoeffRxnMatrix[irxn];
dptr = m_deltaMolNumPhase[irxn];
kspec = ir[irxn];
kspec = m_indexRxnToSpecies[irxn];
int iph = m_phaseID[kspec];
int *pp_ptr = m_phaseParticipation[irxn];
dptr[iph] = 1.0;
@ -4375,7 +4377,7 @@ namespace VCSnonideal {
if (ll < 0) {
for (irxn = 0; irxn < m_numRxnRdc; ++irxn) {
if (m_rxnStatus[irxn] != VCS_SPECIES_MINOR) {
kspec = ir[irxn];
kspec = m_indexRxnToSpecies[irxn];
iphase = m_phaseID[kspec];
if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
#ifdef DEBUG_MODE
@ -4419,7 +4421,7 @@ namespace VCSnonideal {
} else if (ll > 0) {
for (irxn = 0; irxn < m_numRxnRdc; ++irxn) {
if (m_rxnStatus[irxn] == VCS_SPECIES_MINOR) {
kspec = ir[irxn];
kspec = m_indexRxnToSpecies[irxn];
iphase = m_phaseID[kspec];
if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
#ifdef DEBUG_MODE
@ -4791,7 +4793,7 @@ namespace VCSnonideal {
bool zeroedPhase = TRUE;
for (irxn = 0; irxn < irxnl; ++irxn) {
kspec = ir[irxn];
kspec = m_indexRxnToSpecies[irxn];
if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
iph = m_phaseID[kspec];
if (iph == iphase ) {
@ -4848,7 +4850,7 @@ namespace VCSnonideal {
if (zeroedPhase) {
double phaseDG = 1.0;
for (irxn = 0; irxn < irxnl; ++irxn) {
kspec = ir[irxn];
kspec = m_indexRxnToSpecies[irxn];
iph = m_phaseID[kspec];
if (iph == iphase) {
if (m_deltaGRxn_new[irxn] > 50.0) m_deltaGRxn_new[irxn] = 50.0;
@ -4860,7 +4862,7 @@ namespace VCSnonideal {
* Overwrite the individual dg's with the phase DG.
*/
for (irxn = 0; irxn < irxnl; ++irxn) {
kspec = ir[irxn];
kspec = m_indexRxnToSpecies[irxn];
iph = m_phaseID[kspec];
if (iph == iphase) {
m_deltaGRxn_new[irxn] = 1.0 - phaseDG;