Changed variable names
This commit is contained in:
Harry Moffat
parent
6b660b30c9
commit
4769e38ac0
8 changed files with 275 additions and 201 deletions
|
|
@ -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<double> 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);
|
||||
|
|
|
|||
|
|
@ -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 <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <math.h>
|
||||
|
||||
#include "vcs_solve.h"
|
||||
#include "vcs_internal.h"
|
||||
#include "vcs_VolPhase.h"
|
||||
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <math.h>
|
||||
|
||||
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);
|
||||
}
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -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++;
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
}
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -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");
|
||||
|
|
|
|||
|
|
@ -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];
|
||||
}
|
||||
|
||||
/*
|
||||
|
|
|
|||
|
|
@ -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<double> m_deltaGRxn_old;
|
||||
|
||||
//! Temporary vector of Rxn DeltaG's
|
||||
/*!
|
||||
* This is used from time to time, for printing purposes
|
||||
*/
|
||||
std::vector<double> m_deltaGRxn_tmp;
|
||||
|
||||
//! Reaction Adjustments for each species during the current step
|
||||
|
|
@ -525,18 +578,29 @@ public:
|
|||
*/
|
||||
std::vector<double> 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<double> 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<double> m_elemAbundancesGoal;
|
||||
|
||||
std::vector<double> 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<double> 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
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue