diff --git a/Cantera/src/equil/vcs_elem.cpp b/Cantera/src/equil/vcs_elem.cpp index f92c646e4..275dbbfbe 100644 --- a/Cantera/src/equil/vcs_elem.cpp +++ b/Cantera/src/equil/vcs_elem.cpp @@ -21,15 +21,15 @@ void VCS_SOLVE::vcs_elab(void) * * vcs_elab: * - * Computes the elemental abundances vector, ga[], and stores it + * Computes the elemental abundances vector, m_elemAbundances[], and stores it * back into the global structure *************************************************************************/ { for (int j = 0; j < m_numElemConstraints; ++j) { - ga[j] = 0.0; + m_elemAbundances[j] = 0.0; for (int i = 0; i < m_numSpeciesTot; ++i) { if (SpeciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - ga[j] += FormulaMatrix[j][i] * m_molNumSpecies_old[i]; + m_elemAbundances[j] += FormulaMatrix[j][i] * m_molNumSpecies_old[i]; } } } @@ -79,14 +79,14 @@ int VCS_SOLVE::vcs_elabcheck(int ibound) { * Require 12 digits of accuracy on non-zero constraints. */ for (i = 0; i < top; ++i) { - if (fabs(ga[i] - gai[i]) > (fabs(gai[i]) * 1.0e-12)) { + if (fabs(m_elemAbundances[i] - m_elemAbundancesGoal[i]) > (fabs(m_elemAbundancesGoal[i]) * 1.0e-12)) { /* * This logic is for charge neutrality condition */ if (m_elType[i] == VCS_ELEM_TYPE_CHARGENEUTRALITY) { - AssertThrowVCS(gai[i] == 0.0, "vcs_elabcheck"); + AssertThrowVCS(m_elemAbundancesGoal[i] == 0.0, "vcs_elabcheck"); } - if (gai[i] == 0.0 || (m_elType[i] == VCS_ELEM_TYPE_ELECTRONCHARGE)) { + if (m_elemAbundancesGoal[i] == 0.0 || (m_elType[i] == VCS_ELEM_TYPE_ELECTRONCHARGE)) { scale = VCS_DELETE_MINORSPECIES_CUTOFF; /* * Find out if the constraint is a multisign constraint. @@ -107,11 +107,11 @@ int VCS_SOLVE::vcs_elabcheck(int ibound) { } } if (multisign) { - if (fabs(ga[i] - gai[i]) > 1e-11 * scale) { + if (fabs(m_elemAbundances[i] - m_elemAbundancesGoal[i]) > 1e-11 * scale) { return FALSE; } } else { - if (fabs(ga[i] - gai[i]) > VCS_DELETE_MINORSPECIES_CUTOFF) { + if (fabs(m_elemAbundances[i] - m_elemAbundancesGoal[i]) > VCS_DELETE_MINORSPECIES_CUTOFF) { return FALSE; } } @@ -193,7 +193,7 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) * Internal data to be worked on:: * * ga Current element abundances - * gai Required elemental abundances + * m_elemAbundancesGoal Required elemental abundances * m_molNumSpecies_old Current mole number of species. * FormulaMatrix[][] Formular matrix of the species * ne Number of elements @@ -213,7 +213,7 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) #ifdef DEBUG_MODE double l2before = 0.0, l2after = 0.0; std::vector ga_save(m_numElemConstraints, 0.0); - vcs_dcopy(VCS_DATA_PTR(ga_save), VCS_DATA_PTR(ga), m_numElemConstraints); + vcs_dcopy(VCS_DATA_PTR(ga_save), VCS_DATA_PTR(m_elemAbundances), m_numElemConstraints); if (vcs_debug_print_lvl >= 2) { plogf(" --- vcsc_elcorr: Element abundances correction routine"); if (m_numElemConstraints != m_numComponents) { @@ -223,7 +223,7 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) } for (i = 0; i < m_numElemConstraints; ++i) { - x[i] = ga[i] - gai[i]; + x[i] = m_elemAbundances[i] - m_elemAbundancesGoal[i]; } l2before = 0.0; for (i = 0; i < m_numElemConstraints; ++i) { @@ -261,7 +261,7 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) if (SpeciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { double eval = FormulaMatrix[i][kspec]; if (eval > 0.0) { - m_molNumSpecies_old[kspec] = gai[i] / eval; + m_molNumSpecies_old[kspec] = m_elemAbundancesGoal[i] / eval; changed = true; } } @@ -279,7 +279,7 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) } } if (numCompNonZero == 1) { - double diff = gai[i]; + double diff = m_elemAbundancesGoal[i]; for (kspec = m_numComponents; kspec < m_numSpeciesTot; kspec++) { if (SpeciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { double eval = FormulaMatrix[i][kspec]; @@ -313,7 +313,7 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) if (SpeciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { double atomComp = FormulaMatrix[i][kspec]; if (atomComp > 0.0) { - double maxPermissible = gai[i] / atomComp; + double maxPermissible = m_elemAbundancesGoal[i] / atomComp; if (m_molNumSpecies_old[kspec] > maxPermissible) { #ifdef DEBUG_MODE @@ -356,7 +356,7 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) * nc .ne. ne. */ for (i = 0; i < m_numComponents; ++i) { - x[i] = ga[i] - gai[i]; + x[i] = m_elemAbundances[i] - m_elemAbundancesGoal[i]; if (fabs(x[i]) > 1.0E-13) retn = 1; for (j = 0; j < m_numComponents; ++j) { aa[j + i*m_numElemConstraints] = FormulaMatrix[j][i]; @@ -435,7 +435,7 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) saveDir = 0.0; goodSpec = TRUE; for (i = 0; i < m_numComponents; ++i) { - dir = FormulaMatrix[i][kspec] * (gai[i] - ga[i]); + dir = FormulaMatrix[i][kspec] * (m_elemAbundancesGoal[i] - m_elemAbundances[i]); if (fabs(dir) > 1.0E-10) { if (dir > 0.0) { if (saveDir < 0.0) { @@ -461,7 +461,7 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) xx = 0.0; for (i = 0; i < m_numComponents; ++i) { if (FormulaMatrix[i][kspec] != 0.0) { - xx += (gai[i] - ga[i]) / FormulaMatrix[i][kspec]; + xx += (m_elemAbundancesGoal[i] - m_elemAbundances[i]) / FormulaMatrix[i][kspec]; its++; } } @@ -489,11 +489,11 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) for (i = 0; i < m_numElemConstraints; ++i) { if (m_elType[i] == VCS_ELEM_TYPE_CHARGENEUTRALITY || - (m_elType[i] == VCS_ELEM_TYPE_ABSPOS && gai[i] == 0.0)) { + (m_elType[i] == VCS_ELEM_TYPE_ABSPOS && m_elemAbundancesGoal[i] == 0.0)) { for (kspec = 0; kspec < m_numSpeciesRdc; kspec++) { - if (ga[i] > 0.0) { + if (m_elemAbundances[i] > 0.0) { if (FormulaMatrix[i][kspec] < 0.0) { - m_molNumSpecies_old[kspec] -= ga[i] / FormulaMatrix[i][kspec] ; + m_molNumSpecies_old[kspec] -= m_elemAbundances[i] / FormulaMatrix[i][kspec] ; if (m_molNumSpecies_old[kspec] < 0.0) { m_molNumSpecies_old[kspec] = 0.0; } @@ -501,9 +501,9 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) break; } } - if (ga[i] < 0.0) { + if (m_elemAbundances[i] < 0.0) { if (FormulaMatrix[i][kspec] > 0.0) { - m_molNumSpecies_old[kspec] -= ga[i] / FormulaMatrix[i][kspec]; + m_molNumSpecies_old[kspec] -= m_elemAbundances[i] / FormulaMatrix[i][kspec]; if (m_molNumSpecies_old[kspec] < 0.0) { m_molNumSpecies_old[kspec] = 0.0; } @@ -525,7 +525,7 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) * electron charge exactly. */ for (i = 0; i < m_numElemConstraints; ++i) { - double dev = gai[i] - ga[i]; + double dev = m_elemAbundancesGoal[i] - m_elemAbundances[i]; if (m_elType[i] == VCS_ELEM_TYPE_ELECTRONCHARGE && (fabs(dev) > 1.0E-300)) { bool useZeroed = true; for (kspec = 0; kspec < m_numSpeciesRdc; kspec++) { @@ -581,7 +581,7 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) #ifdef DEBUG_MODE l2after = 0.0; for (i = 0; i < m_numElemConstraints; ++i) { - l2after += SQUARE(ga[i] - gai[i]); + l2after += SQUARE(m_elemAbundances[i] - m_elemAbundancesGoal[i]); } l2after = sqrt(l2after/m_numElemConstraints); if (vcs_debug_print_lvl >= 2) { @@ -589,7 +589,7 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) " Final\n"); for (i = 0; i < m_numElemConstraints; ++i) { plogf(" --- "); plogf("%-2.2s", ElName[i].c_str()); - plogf(" %20.12E %20.12E %20.12E\n", gai[i], ga_save[i], ga[i]); + plogf(" %20.12E %20.12E %20.12E\n", m_elemAbundancesGoal[i], ga_save[i], m_elemAbundances[i]); } plogf(" --- Diff_Norm: %20.12E %20.12E\n", l2before, l2after); diff --git a/Cantera/src/equil/vcs_elem_rearrange.cpp b/Cantera/src/equil/vcs_elem_rearrange.cpp index bb9fcbcdd..a11c25ffa 100644 --- a/Cantera/src/equil/vcs_elem_rearrange.cpp +++ b/Cantera/src/equil/vcs_elem_rearrange.cpp @@ -1,3 +1,8 @@ +/** + * @file vcs_elem_rearrange.cpp + * Contains implementations for rearranging the element columns, and + * it contains the algorithm for choosing the rearrangement. + */ /* $Author$ * $Date$ * $Revision$ @@ -9,18 +14,19 @@ * U.S. Government retains certain rights in this software. */ -#include -#include -#include - #include "vcs_solve.h" #include "vcs_internal.h" #include "vcs_VolPhase.h" +#include +#include +#include + namespace VCSnonideal { + // Rearrange the constraint equations represented by the Formula + // Matrix so that the operational ones are in the front /* - * vcs_elem_rearrange: * * This subroutine handles the rearrangement of the constraint * equations represented by the Formula Matrix. Rearrangement is only @@ -55,125 +61,127 @@ namespace VCSnonideal { * sm[i+j*ne] = QR matrix work space (ne*ne in length) * */ -int VCS_SOLVE::vcs_elem_rearrange(double *aw, double *sa, double *sm, - double *ss) { - int j, k, l, i, jl, ml, jr, lindep, ielem; - int ncomponents = m_numComponents; - double test = -1.0E10; + int VCS_SOLVE::vcs_elem_rearrange(double * const aw, double * const sa, + double * const sm, double * const ss) { + int j, k, l, i, jl, ml, jr, lindep, ielem; + int ncomponents = m_numComponents; + double test = -1.0E10; #ifdef DEBUG_MODE - if (vcs_debug_print_lvl >= 2) { - plogf(" "); for(i=0; i<77; i++) plogf("-"); plogf("\n"); - plogf(" --- Subroutine elem_rearrange() called to "); - plogf("check stoich. coefficent matrix\n"); - plogf(" --- and to rearrange the element ordering once\n"); - } + if (vcs_debug_print_lvl >= 2) { + plogf(" "); for(i=0; i<77; i++) plogf("-"); plogf("\n"); + plogf(" --- Subroutine elem_rearrange() called to "); + plogf("check stoich. coefficent matrix\n"); + plogf(" --- and to rearrange the element ordering once"); + plogendl(); + } #endif - /* - * Use a temporary work array for the element numbers - * Also make sure the value of test is unique. - */ - lindep = FALSE; - do { + /* + * Use a temporary work array for the element numbers + * Also make sure the value of test is unique. + */ + lindep = FALSE; + do { lindep = FALSE; for (i = 0; i < m_numElemConstraints; ++i) { test -= 1.0; - aw[i] = gai[i]; + aw[i] = m_elemAbundancesGoal[i]; if (test == aw[i]) lindep = TRUE; } - } while (lindep); + } while (lindep); - /* - * Top of a loop of some sort based on the index JR. JR is the - * current number independent elements found. - */ - jr = -1; - do { + /* + * Top of a loop of some sort based on the index JR. JR is the + * current number independent elements found. + */ + jr = -1; + do { ++jr; /* - * Top of another loop point based on finding a linearly - * independent species - */ + * Top of another loop point based on finding a linearly + * independent species + */ do { - /* + /* * Search the remaining part of the mole fraction vector, AW, * for the largest remaining species. Return its identity in K. */ - k = m_numElemConstraints; - for (ielem = jr; ielem < m_numElemConstraints; ielem++) { - if (ElActive[ielem]) { - if (aw[ielem] != test) { - k = ielem; - break; - } - } - } - if (k == m_numElemConstraints) { - plogf("Shouldn't be here\n"); - exit(-1); - } + k = m_numElemConstraints; + for (ielem = jr; ielem < m_numElemConstraints; ielem++) { + if (ElActive[ielem]) { + if (aw[ielem] != test) { + k = ielem; + break; + } + } + } + if (k == m_numElemConstraints) { + plogf("vcs_elem_rearrange::Shouldn't be here. Algorithm misfired."); + plogendl(); + exit(-1); + } - /* + /* * Assign a large negative number to the element that we have * just found, in order to take it out of further consideration. */ - aw[k] = test; + aw[k] = test; - /* *********************************************************** */ - /* **** CHECK LINEAR INDEPENDENCE OF CURRENT FORMULA MATRIX */ - /* **** LINE WITH PREVIOUS LINES OF THE FORMULA MATRIX ****** */ - /* *********************************************************** */ - /* + /* *********************************************************** */ + /* **** CHECK LINEAR INDEPENDENCE OF CURRENT FORMULA MATRIX */ + /* **** LINE WITH PREVIOUS LINES OF THE FORMULA MATRIX ****** */ + /* *********************************************************** */ + /* * Modified Gram-Schmidt Method, p. 202 Dalquist * QR factorization of a matrix without row pivoting. */ - jl = jr; - /* + jl = jr; + /* * Fill in the row for the current element, k, under consideration * The row will contain the Formula matrix value for that element * from the current component. */ - for (j = 0; j < ncomponents; ++j) { - sm[j + jr*ncomponents] = FormulaMatrix[k][j]; - } - if (jl > 0) { - /* - * Compute the coefficients of JA column of the - * the upper triangular R matrix, SS(J) = R_J_JR - * (this is slightly different than Dalquist) - * R_JA_JA = 1 - */ - for (j = 0; j < jl; ++j) { - ss[j] = 0.0; - for (i = 0; i < ncomponents; ++i) { - ss[j] += sm[i + jr*ncomponents] * sm[i + j*ncomponents]; - } - ss[j] /= sa[j]; + for (j = 0; j < ncomponents; ++j) { + sm[j + jr*ncomponents] = FormulaMatrix[k][j]; + } + if (jl > 0) { + /* + * Compute the coefficients of JA column of the + * the upper triangular R matrix, SS(J) = R_J_JR + * (this is slightly different than Dalquist) + * R_JA_JA = 1 + */ + for (j = 0; j < jl; ++j) { + ss[j] = 0.0; + for (i = 0; i < ncomponents; ++i) { + ss[j] += sm[i + jr*ncomponents] * sm[i + j*ncomponents]; } - /* - * Now make the new column, (*,JR), orthogonal to the - * previous columns - */ - for (j = 0; j < jl; ++j) { - for (l = 0; l < ncomponents; ++l) { - sm[l + jr*ncomponents] -= ss[j] * sm[l + j*ncomponents]; - } + ss[j] /= sa[j]; + } + /* + * Now make the new column, (*,JR), orthogonal to the + * previous columns + */ + for (j = 0; j < jl; ++j) { + for (l = 0; l < ncomponents; ++l) { + sm[l + jr*ncomponents] -= ss[j] * sm[l + j*ncomponents]; } - } + } + } - /* + /* * Find the new length of the new column in Q. * It will be used in the denominator in future row calcs. */ - sa[jr] = 0.0; - for (ml = 0; ml < ncomponents; ++ml) { - sa[jr] += SQUARE(sm[ml + jr*ncomponents]); - } - /* **************************************************** */ - /* **** IF NORM OF NEW ROW .LT. 1E-6 REJECT ********** */ - /* **************************************************** */ - if (sa[jr] < 1.0e-6) lindep = TRUE; - else lindep = FALSE; + sa[jr] = 0.0; + for (ml = 0; ml < ncomponents; ++ml) { + sa[jr] += SQUARE(sm[ml + jr*ncomponents]); + } + /* **************************************************** */ + /* **** IF NORM OF NEW ROW .LT. 1E-6 REJECT ********** */ + /* **************************************************** */ + if (sa[jr] < 1.0e-6) lindep = TRUE; + else lindep = FALSE; } while(lindep); /* ****************************************** */ /* **** REARRANGE THE DATA ****************** */ @@ -182,69 +190,72 @@ int VCS_SOLVE::vcs_elem_rearrange(double *aw, double *sa, double *sm, #ifdef DEBUG_MODE if (vcs_debug_print_lvl >= 2) { plogf(" --- "); plogf("%-2.2s", (ElName[k]).c_str()); - plogf("(%9.2g) replaces ", gai[k]); + plogf("(%9.2g) replaces ", m_elemAbundancesGoal[k]); plogf("%-2.2s", (ElName[jr]).c_str()); - plogf("(%9.2g) as element %3d\n", gai[jr], jr); + plogf("(%9.2g) as element %3d", m_elemAbundancesGoal[jr], jr); + plogendl(); } #endif - vcs_switch_elem_pos(jr, k); - vcsUtil_dsw(aw, jr, k); + vcs_switch_elem_pos(jr, k); + vcsUtil_dsw(aw, jr, k); } /* - * If we haven't found enough components, go back - * and find some more. (nc -1 is used below, because - * jr is counted from 0, via the C convention. - */ - } while (jr < (ncomponents-1)); - return VCS_SUCCESS; -} /* vcs_elem_rearrange() ****************************************************/ + * If we haven't found enough components, go back + * and find some more. (nc -1 is used below, because + * jr is counted from 0, via the C convention. + */ + } while (jr < (ncomponents-1)); + return VCS_SUCCESS; + } -// Swaps the indecises for all of the global data for two elements, ipos -// and jpos. -/* - * This function knows all of the element information with VCS_SOLVE, and - * can therefore switch element positions - * - * @param ipos first global element index - * @param jpos second global element index - */ -void VCS_SOLVE::vcs_switch_elem_pos(int ipos, int jpos) { - if (ipos == jpos) return; - int j; - double dtmp; - vcs_VolPhase *volPhase; + // Swaps the indecises for all of the global data for two elements, ipos + // and jpos. + /* + * This function knows all of the element information with VCS_SOLVE, and + * can therefore switch element positions + * + * @param ipos first global element index + * @param jpos second global element index + */ + void VCS_SOLVE::vcs_switch_elem_pos(int ipos, int jpos) { + if (ipos == jpos) return; + int j; + double dtmp; + vcs_VolPhase *volPhase; #ifdef DEBUG_MODE - if (ipos < 0 || ipos > (m_numElemConstraints - 1) || - jpos < 0 || jpos > (m_numElemConstraints - 1) ) { + if (ipos < 0 || ipos > (m_numElemConstraints - 1) || + jpos < 0 || jpos > (m_numElemConstraints - 1) ) { plogf("vcs_switch_elem_pos: ifunc = 0: inappropriate args: %d %d\n", - ipos, jpos); - } + ipos, jpos); + plogendl(); + exit(-1); + } #endif - /* - * Change the element Global Index list in each phase object - * to reflect the switch in the element positions. - */ - for (int iph = 0; iph < NPhase; iph++) { - volPhase = VPhaseList[iph]; - for (int e = 0; e < volPhase->nElemConstraints; e++) { - if (volPhase->ElGlobalIndex[e] == ipos) { + /* + * Change the element Global Index list in each phase object + * to reflect the switch in the element positions. + */ + for (int iph = 0; iph < NPhase; iph++) { + volPhase = VPhaseList[iph]; + for (int e = 0; e < volPhase->nElemConstraints; e++) { + if (volPhase->ElGlobalIndex[e] == ipos) { volPhase->ElGlobalIndex[e] = jpos; - } - if (volPhase->ElGlobalIndex[e] == jpos) { + } + if (volPhase->ElGlobalIndex[e] == jpos) { volPhase->ElGlobalIndex[e] =ipos; - } - } - } - vcsUtil_dsw(VCS_DATA_PTR(gai), ipos, jpos); - vcsUtil_dsw(VCS_DATA_PTR(ga), ipos, jpos); - vcsUtil_isw(VCS_DATA_PTR(IndEl), ipos, jpos); - vcsUtil_isw(VCS_DATA_PTR(m_elType), ipos, jpos); - vcsUtil_isw(VCS_DATA_PTR(ElActive), ipos, jpos); - for (j = 0; j < m_numSpeciesTot; ++j) { + } + } + } + vcsUtil_dsw(VCS_DATA_PTR(m_elemAbundancesGoal), ipos, jpos); + vcsUtil_dsw(VCS_DATA_PTR(m_elemAbundances), ipos, jpos); + vcsUtil_isw(VCS_DATA_PTR(IndEl), ipos, jpos); + vcsUtil_isw(VCS_DATA_PTR(m_elType), ipos, jpos); + vcsUtil_isw(VCS_DATA_PTR(ElActive), ipos, jpos); + for (j = 0; j < m_numSpeciesTot; ++j) { SWAP(FormulaMatrix[ipos][j], FormulaMatrix[jpos][j], dtmp); - } - vcsUtil_stsw(ElName, ipos, jpos); -} /* vcs_switch_elem_pos() ***************************************************/ + } + vcsUtil_stsw(ElName, ipos, jpos); + } } diff --git a/Cantera/src/equil/vcs_inest.cpp b/Cantera/src/equil/vcs_inest.cpp index 93aa4e6ac..f5fa65f15 100644 --- a/Cantera/src/equil/vcs_inest.cpp +++ b/Cantera/src/equil/vcs_inest.cpp @@ -120,7 +120,7 @@ namespace VCSnonideal { tmp += FormulaMatrix[j][kspec] * molNum[kspec]; } plogf("%s ", pprefix); plogf(" %-9.9s", (ElName[j]).c_str()); - plogf(" %12.3g %12.3g\n", gai[j], tmp); + plogf(" %12.3g %12.3g\n", m_elemAbundancesGoal[j], tmp); jj++; } } diff --git a/Cantera/src/equil/vcs_nondim.cpp b/Cantera/src/equil/vcs_nondim.cpp index 03dfc1dc8..feff08e42 100644 --- a/Cantera/src/equil/vcs_nondim.cpp +++ b/Cantera/src/equil/vcs_nondim.cpp @@ -117,7 +117,7 @@ void VCS_SOLVE::vcs_nondim_TP(void) { } } for (i = 0; i < m_numElemConstraints; ++i) { - gai[i] *= 1.0E3; + m_elemAbundancesGoal[i] *= 1.0E3; } } } @@ -157,7 +157,7 @@ void VCS_SOLVE::vcs_redim_TP(void) } } for (i = 0; i < m_numElemConstraints; ++i) { - gai[i] /= 1.0E3; + m_elemAbundancesGoal[i] /= 1.0E3; } } diff --git a/Cantera/src/equil/vcs_report.cpp b/Cantera/src/equil/vcs_report.cpp index 880717de2..69ff0afaa 100644 --- a/Cantera/src/equil/vcs_report.cpp +++ b/Cantera/src/equil/vcs_report.cpp @@ -234,7 +234,7 @@ int VCS_SOLVE::vcs_report(int iconv) plogf(" | |\n"); plogf(" PhaseName | MolTarget |"); for (j = 0; j < m_numElemConstraints; j++) { - plogf(" %10.3g", gai[j]); + plogf(" %10.3g", m_elemAbundancesGoal[j]); } plogf(" | Gibbs Total |\n"); print_line("-", m_numElemConstraints*10 + 58); @@ -290,7 +290,7 @@ int VCS_SOLVE::vcs_report(int iconv) plogf(" Actual Target Type ElActive\n"); for (i = 0; i < m_numElemConstraints; ++i) { print_space(26); plogf("%-2.2s", (ElName[i]).c_str()); - plogf("%20.12E %20.12E", ga[i], gai[i]); + plogf("%20.12E %20.12E", m_elemAbundances[i], m_elemAbundancesGoal[i]); plogf(" %3d %3d\n", m_elType[i], ElActive[i]); } plogf("\n"); diff --git a/Cantera/src/equil/vcs_solve.cpp b/Cantera/src/equil/vcs_solve.cpp index 1499bde1e..6e6ed42c3 100644 --- a/Cantera/src/equil/vcs_solve.cpp +++ b/Cantera/src/equil/vcs_solve.cpp @@ -128,9 +128,8 @@ namespace VCSnonideal { m_deltaMolNumSpecies.resize(nspecies0, 0.0); m_feSpecies_old.resize(nspecies0, 0.0); - ga.resize(nelements, 0.0); - gai.resize(nelements, 0.0); - + m_elemAbundances.resize(nelements, 0.0); + m_elemAbundancesGoal.resize(nelements, 0.0); TPhMoles.resize(nphase0, 0.0); TPhMoles1.resize(nphase0, 0.0); @@ -550,22 +549,22 @@ namespace VCSnonideal { } /* - * Formulate the Goal Element Abundance Vector, gai[] + * Formulate the Goal Element Abundance Vector */ if (pub->gai.size() != 0) { - for (i = 0; i < nelements; i++) gai[i] = pub->gai[i]; + for (i = 0; i < nelements; i++) m_elemAbundancesGoal[i] = pub->gai[i]; } else { if (iest == 0) { for (j = 0; j < nelements; j++) { - gai[j] = 0.0; + m_elemAbundancesGoal[j] = 0.0; for (kspec = 0; kspec < nspecies; kspec++) { if (SpeciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - gai[j] += FormulaMatrix[j][kspec] * m_molNumSpecies_old[kspec]; + m_elemAbundancesGoal[j] += FormulaMatrix[j][kspec] * m_molNumSpecies_old[kspec]; } } } } else { - plogf("%sElement Abundances, gai[], not specified\n", ser); + plogf("%sElement Abundances, m_elemAbundancesGoal[], not specified\n", ser); return VCS_PUB_BAD; } } @@ -797,7 +796,7 @@ namespace VCSnonideal { */ for (i = 0; i < m_numElemConstraints; i++) { j = IndEl[i]; - gai[i] = pub->gai[j]; + m_elemAbundancesGoal[i] = pub->gai[j]; } /* diff --git a/Cantera/src/equil/vcs_solve.h b/Cantera/src/equil/vcs_solve.h index 8f1f4a56c..f6c65671e 100644 --- a/Cantera/src/equil/vcs_solve.h +++ b/Cantera/src/equil/vcs_solve.h @@ -197,7 +197,55 @@ public: bool vcs_wellPosed(VCS_PROB *vprob); - int vcs_elem_rearrange(double *aw, double *sa, double *sm, double *ss); + //! Rearrange the constraint equations represented by the Formula + //! Matrix so that the operational ones are in the front + /*! + * + * This subroutine handles the rearrangement of the constraint + * equations represented by the Formula Matrix. Rearrangement is only + * necessary when the number of components is less than the number of + * elements. For this case, some constraints can never be satisfied + * exactly, because the range space represented by the Formula + * Matrix of the components can't span the extra space. These + * constraints, which are out of the range space of the component + * Formula matrix entries, are migrated to the back of the Formula + * matrix. + * + * A prototypical example is an extra element column in + * FormulaMatrix[], + * which is identically zero. For example, let's say that argon is + * has an element column in FormulaMatrix[], but no species in the + * mechanism + * actually contains argon. Then, nc < ne. Also, without perturbation + * of FormulaMatrix[] vcs_basopt[] would produce a zero pivot + * because the matrix + * would be singular (unless the argon element column was already the + * last column of FormulaMatrix[]. + * This routine borrows heavily from vcs_basopt's algorithm. It + * finds nc constraints which span the range space of the Component + * Formula matrix, and assigns them as the first nc components in the + * formular matrix. This guarrantees that vcs_basopt[] has a + * nonsingular matrix to invert. + * + * Other Variables + * @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) + * + */ + int vcs_elem_rearrange(double *const aw, double * const sa, + double * const sm, double * const ss); + + //! Swaps the indecises for all of the global data for two elements, ipos + //! and jpos. + /*! + * This function knows all of the element information with VCS_SOLVE, and + * can therefore switch element positions + * + * @param ipos first global element index + * @param jpos second global element index + */ void vcs_switch_elem_pos(int ipos, int jpos); int vcs_rxn_adj_cg(void); @@ -212,6 +260,7 @@ public: int vcs_report(int); + int vcs_rearrange(void); @@ -516,6 +565,10 @@ public: //! Last deltag[irxn] from the previous step std::vector m_deltaGRxn_old; + //! Temporary vector of Rxn DeltaG's + /*! + * This is used from time to time, for printing purposes + */ std::vector m_deltaGRxn_tmp; //! Reaction Adjustments for each species during the current step @@ -525,18 +578,29 @@ public: */ std::vector m_deltaMolNumSpecies; + //! Element abundances vector + /*! + * Vector of moles of each element actually in the solution + * vector. Except for certain parts of the algorithm, + * this is a constant. + * Note other constraint conditions are added to this vector. + * This is input from the input file and + * is considered a constant from thereon. + * units = gmoles + */ + std::vector m_elemAbundances; + + //! Element abundances vector Goals + /*! + * Vector of moles of each element that are the goals of the + * simulation. This is a constant in the problem. + * Note other constraint conditions are added to this vector. + * This is input from the input file and + * is considered a constant from thereon. + * units = gmoles + */ + std::vector m_elemAbundancesGoal; - std::vector ga; /* ga[j] = Element abundances for jth element from - * estimate - * -> this is calculated from the current mole - * fraction vector and BM, the formula - * vector. - * units = gmoles */ - std::vector gai; /* gai[j] = Element abundances for jth element - * -> corrected - * -> this is input from the input file and - * is considered a constant from thereon. - * units = gmoles */ double TMoles; /* TMoles = Total number of moles in all phases * This number includes the inerts. * -> Don't use this except for scaling diff --git a/Cantera/src/equil/vcs_solve_TP.cpp b/Cantera/src/equil/vcs_solve_TP.cpp index d3e442bf4..b0fa33433 100644 --- a/Cantera/src/equil/vcs_solve_TP.cpp +++ b/Cantera/src/equil/vcs_solve_TP.cpp @@ -250,7 +250,7 @@ namespace VCSnonideal { plogf(" FROM ESTIMATE Type\n\n"); for (i = 0; i < m_numElemConstraints; ++i) { print_space(26); plogf("%-2.2s", (ElName[i]).c_str()); - plogf("%20.12E%20.12E %3d\n", gai[i], ga[i], m_elType[i]); + plogf("%20.12E%20.12E %3d\n", m_elemAbundancesGoal[i], m_elemAbundances[i], m_elType[i]); } if (iest < 0) { plogf("\n MODIFIED LINEAR PROGRAMMING ESTIMATE OF EQUILIBRIUM\n"); @@ -596,7 +596,7 @@ namespace VCSnonideal { if (elType == VCS_ELEM_TYPE_ABSPOS) { double atomComp = FormulaMatrix[j][kspec]; if (atomComp > 0.0) { - double maxPermissible = gai[j] / atomComp; + double maxPermissible = m_elemAbundancesGoal[j] / atomComp; if (maxPermissible < VCS_DELETE_MINORSPECIES_CUTOFF) { #ifdef DEBUG_MODE sprintf(ANOTE, "Species stays zeroed even though dG neg, because of %s elemAbund", @@ -3484,7 +3484,7 @@ namespace VCSnonideal { double nu = FormulaMatrix[j][kspec]; if (nu != 0.0) { nonZeroesKspec++; - maxConcPossKspec = MIN(gai[j] / nu, maxConcPossKspec); + maxConcPossKspec = MIN(m_elemAbundancesGoal[j] / nu, maxConcPossKspec); } } } @@ -3657,7 +3657,7 @@ namespace VCSnonideal { * Then, the first row in sm[], below will be indentically * zero. bleh. * What needs to be done is to perform a rearrangement - * of the ELEMENTS -> i.e. rearrange, FormulaMatrix, sp, and gai, such + * of the ELEMENTS -> i.e. rearrange, FormulaMatrix, sp, and m_elemAbundancesGoal, such * that the first nc elements form in combination with the * nc components create an invertible sm[]. not a small * project, but very doable. @@ -3944,7 +3944,7 @@ namespace VCSnonideal { if (elType == VCS_ELEM_TYPE_ABSPOS) { double atomComp = FormulaMatrix[j][kspec]; if (atomComp > 0.0) { - double maxPermissible = gai[j] / atomComp; + double maxPermissible = m_elemAbundancesGoal[j] / atomComp; if (maxPermissible < VCS_DELETE_MINORSPECIES_CUTOFF) { #ifdef DEBUG_MODE if (vcs_debug_print_lvl >= 2) { @@ -4474,9 +4474,9 @@ namespace VCSnonideal { for (j = 0; j < m_numElemConstraints; ++j) { plogf(" "); plogf("%-2.2s", (ElName[j]).c_str()); plogf(" = %15.6E %15.6E %15.6E %3d\n", - eav[j], gai[j], eav[j] - gai[j], m_elType[j]); - if (gai[j] != 0.) { - if (fabs(eav[j] - gai[j]) > gai[j] * 5.0e-9) + eav[j], m_elemAbundancesGoal[j], eav[j] - m_elemAbundancesGoal[j], m_elType[j]); + if (m_elemAbundancesGoal[j] != 0.) { + if (fabs(eav[j] - m_elemAbundancesGoal[j]) > m_elemAbundancesGoal[j] * 5.0e-9) kerr = TRUE; } else { if (fabs(eav[j]) > 1.0e-10) kerr = TRUE;