diff --git a/Cantera/src/equil/vcs_inest.cpp b/Cantera/src/equil/vcs_inest.cpp index 2e144e6b9..87141ac96 100644 --- a/Cantera/src/equil/vcs_inest.cpp +++ b/Cantera/src/equil/vcs_inest.cpp @@ -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 diff --git a/Cantera/src/equil/vcs_prep.cpp b/Cantera/src/equil/vcs_prep.cpp index ec951e170..e5eaae50f 100644 --- a/Cantera/src/equil/vcs_prep.cpp +++ b/Cantera/src/equil/vcs_prep.cpp @@ -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; } } diff --git a/Cantera/src/equil/vcs_report.cpp b/Cantera/src/equil/vcs_report.cpp index 916c4067f..ba592dace 100644 --- a/Cantera/src/equil/vcs_report.cpp +++ b/Cantera/src/equil/vcs_report.cpp @@ -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 sortindex(nspecies,0); - std::vector xy(nspecies,0.0); + std::vector sortindex(nspecies,0); + std::vector 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 gaPhase(m_numElemConstraints, 0.0); - std::vector 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 gaPhase(m_numElemConstraints, 0.0); + std::vector 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 "); + } } -} -/*****************************************************************************/ -/*****************************************************************************/ -/*****************************************************************************/ + /*****************************************************************************/ + /*****************************************************************************/ + /*****************************************************************************/ } diff --git a/Cantera/src/equil/vcs_rxnadj.cpp b/Cantera/src/equil/vcs_rxnadj.cpp index e60bcbfe0..1f1c51ead 100644 --- a/Cantera/src/equil/vcs_rxnadj.cpp +++ b/Cantera/src/equil/vcs_rxnadj.cpp @@ -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]; diff --git a/Cantera/src/equil/vcs_solve.cpp b/Cantera/src/equil/vcs_solve.cpp index ebc2ff553..da617939b 100644 --- a/Cantera/src/equil/vcs_solve.cpp +++ b/Cantera/src/equil/vcs_solve.cpp @@ -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); diff --git a/Cantera/src/equil/vcs_solve.h b/Cantera/src/equil/vcs_solve.h index b9fd68e4e..b98013494 100644 --- a/Cantera/src/equil/vcs_solve.h +++ b/Cantera/src/equil/vcs_solve.h @@ -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 ir; + std::vector m_indexRxnToSpecies; //! Major -Minor status vector for the formation reaction /*! diff --git a/Cantera/src/equil/vcs_solve_TP.cpp b/Cantera/src/equil/vcs_solve_TP.cpp index 5e20da319..da4b79a4b 100644 --- a/Cantera/src/equil/vcs_solve_TP.cpp +++ b/Cantera/src/equil/vcs_solve_TP.cpp @@ -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;