Significant Update to the Equilibrium Solver
Added a capability to solve for the mathematical condition for the birth of a multispecies phase. Added the logic into the algorithm.
This commit is contained in:
Harry Moffat
parent
f5e6b47e7b
commit
ef278fc885
8 changed files with 1159 additions and 461 deletions
|
|
@ -1107,6 +1107,7 @@ namespace VCSnonideal {
|
|||
resize(VP_ID, nsp, nelem, PhaseName.c_str());
|
||||
}
|
||||
TP_ptr->getMoleFractions(VCS_DATA_PTR(Xmol));
|
||||
fractionCreationDelta_ = Xmol;
|
||||
_updateMoleFractionDependencies();
|
||||
|
||||
/*
|
||||
|
|
|
|||
|
|
@ -104,7 +104,7 @@ namespace VCSnonideal {
|
|||
//! Cutoff relative moles below which a phase is deleted
|
||||
//! from the equilibrium problem.
|
||||
#ifndef VCS_DELETE_PHASE_CUTOFF
|
||||
#define VCS_DELETE_PHASE_CUTOFF 1.0e-11
|
||||
#define VCS_DELETE_PHASE_CUTOFF 1.0e-12
|
||||
#endif
|
||||
|
||||
//@}
|
||||
|
|
|
|||
|
|
@ -22,90 +22,683 @@
|
|||
using namespace std;
|
||||
|
||||
namespace VCSnonideal {
|
||||
|
||||
|
||||
// Utility function that evaluates whether a phase can be popped
|
||||
// into existence
|
||||
/*
|
||||
* @param iphasePop id of the phase, which is currently zeroed,
|
||||
*
|
||||
* @return Returns true if the phase can come into existence
|
||||
* and false otherwise.
|
||||
*/
|
||||
bool VCS_SOLVE::vcs_popPhasePossible(const int iphasePop) const {
|
||||
|
||||
vcs_VolPhase *Vphase = m_VolPhaseList[iphasePop];
|
||||
|
||||
#ifdef DEBUG_MODE
|
||||
int existence = Vphase->exists();
|
||||
if (existence > 0) {
|
||||
printf("ERROR vcs_popPhasePossible called for a phase that exists!");
|
||||
exit(-1);
|
||||
}
|
||||
#endif
|
||||
|
||||
/*
|
||||
* section to do damping of the m_deltaMolNumSpecies[]
|
||||
*/
|
||||
for (int k = 0; k < Vphase->nSpecies(); k++) {
|
||||
int kspec = Vphase->spGlobalIndexVCS(k);
|
||||
int irxn = kspec - m_numComponents;
|
||||
if (irxn >= 0) {
|
||||
for (int j = 0; j < m_numComponents; ++j) {
|
||||
if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) {
|
||||
double stoicC = m_stoichCoeffRxnMatrix[irxn][j];
|
||||
if (stoicC != 0.0) {
|
||||
double negChangeComp = - stoicC * 1.0;
|
||||
if (negChangeComp > 0.0) {
|
||||
// TODO: We may have to come up with a tolerance here
|
||||
if (m_molNumSpecies_old[j] <= 1.0E-300) {
|
||||
#ifdef DEBUG_MODE
|
||||
if (m_debug_print_lvl >= 3) {
|
||||
plogf(" --- vcs_popPhasePosssible() Phase %d (%s) can't be popped\n", iphasePop,
|
||||
Vphase->PhaseName.c_str());
|
||||
plogf(" --- Component %d (%s)will go negative\n", j, m_speciesName[j].c_str());
|
||||
}
|
||||
#endif
|
||||
return false;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
// Decision as to whether a phase pops back into existence
|
||||
/*
|
||||
* @return returns the phase id of the phase that pops back into
|
||||
* existence. Returns -1 if there are no phases
|
||||
*/
|
||||
int VCS_SOLVE::vcs_popPhaseID() {
|
||||
int iphasePop = -1;
|
||||
int iph;
|
||||
int irxn, kspec;
|
||||
doublereal FephaseMax = -1.0E30;
|
||||
doublereal Fephase = -1.0E30;
|
||||
vcs_VolPhase *Vphase = 0;
|
||||
|
||||
|
||||
int VCS_SOLVE::vcs_phaseStabilityTest(const int iph) {
|
||||
#ifdef DEBUG_MODE
|
||||
char anote[128];
|
||||
if (m_debug_print_lvl >= 2) {
|
||||
plogf(" --- vcs_popPhaseID() called\n");
|
||||
plogf(" --- Phase Status F_e MoleNum\n");
|
||||
plogf(" --------------------------------------------------------------\n");
|
||||
}
|
||||
#endif
|
||||
for (iph = 0; iph < m_numPhases; iph++) {
|
||||
Vphase = m_VolPhaseList[iph];
|
||||
int existence = Vphase->exists();
|
||||
strcpy(anote, "");
|
||||
if (existence > 0) {
|
||||
|
||||
#ifdef DEBUG_MODE
|
||||
if (m_debug_print_lvl >= 2) {
|
||||
plogf(" --- %18s %5d NA %11.3e\n",
|
||||
Vphase->PhaseName.c_str(),
|
||||
existence,
|
||||
m_tPhaseMoles_old[iph]);
|
||||
}
|
||||
#endif
|
||||
} else {
|
||||
if (Vphase->m_singleSpecies) {
|
||||
/***********************************************************************
|
||||
*
|
||||
* Single Phase Stability Resolution
|
||||
*
|
||||
***********************************************************************/
|
||||
kspec = Vphase->spGlobalIndexVCS(0);
|
||||
irxn = kspec - m_numComponents;
|
||||
doublereal deltaGRxn = m_deltaGRxn_old[irxn];
|
||||
Fephase = exp(-deltaGRxn) - 1.0;
|
||||
if (Fephase > 0.0) {
|
||||
#ifdef DEBUG_MODE
|
||||
strcpy(anote," (ready to be birthed)");
|
||||
#endif
|
||||
if (Fephase > FephaseMax) {
|
||||
iphasePop = iph;
|
||||
FephaseMax = Fephase;
|
||||
#ifdef DEBUG_MODE
|
||||
strcpy(anote," (chosen to be birthed)");
|
||||
#endif
|
||||
}
|
||||
}
|
||||
|
||||
#ifdef DEBUG_MODE
|
||||
if (m_debug_print_lvl >= 2) {
|
||||
plogf(" --- %18s %5d NA %11.3g %s\n",
|
||||
Vphase->PhaseName.c_str(),
|
||||
existence, Fephase,
|
||||
m_tPhaseMoles_old[iph], anote);
|
||||
}
|
||||
#endif
|
||||
|
||||
} else {
|
||||
/***********************************************************************
|
||||
*
|
||||
* MultiSpecies Phase Stability Resolution
|
||||
*
|
||||
***********************************************************************/
|
||||
if (vcs_popPhasePossible(iph)) {
|
||||
Fephase = vcs_phaseStabilityTest(iph);
|
||||
if (Fephase > 0.0) {
|
||||
if (Fephase > FephaseMax) {
|
||||
iphasePop = iph;
|
||||
FephaseMax = Fephase;
|
||||
}
|
||||
} else {
|
||||
if (Fephase > FephaseMax) {
|
||||
FephaseMax = Fephase;
|
||||
}
|
||||
}
|
||||
#ifdef DEBUG_MODE
|
||||
if (m_debug_print_lvl >= 2) {
|
||||
plogf(" --- %18s %5d %11.3g %11.3g\n",
|
||||
Vphase->PhaseName.c_str(),
|
||||
existence, Fephase,
|
||||
m_tPhaseMoles_old[iph]);
|
||||
}
|
||||
#endif
|
||||
} else {
|
||||
#ifdef DEBUG_MODE
|
||||
if (m_debug_print_lvl >= 2) {
|
||||
plogf(" --- %18s %5d blocked %11.3g\n",
|
||||
Vphase->PhaseName.c_str(),
|
||||
existence, m_tPhaseMoles_old[iph]);
|
||||
}
|
||||
#endif
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
#ifdef DEBUG_MODE
|
||||
if (m_debug_print_lvl >= 2) {
|
||||
plogf(" --------------------------------------------------------------\n");
|
||||
}
|
||||
#endif
|
||||
return iphasePop;
|
||||
}
|
||||
|
||||
// Calculates the deltas of the reactions due to phases popping
|
||||
// into existence
|
||||
/*
|
||||
* @param iphasePop Phase id of the phase that will come into existence
|
||||
*
|
||||
* Output
|
||||
* -------
|
||||
* m_deltaMolNumSpecies(irxn) : reaction adjustments, where irxn refers
|
||||
* to the irxn'th species
|
||||
* formation reaction. This adjustment
|
||||
* is for species
|
||||
* irxn + M, where M is the number
|
||||
* of components.
|
||||
*
|
||||
* @return Returns an int representing the status of the step
|
||||
* - 0 : normal return
|
||||
* - 1 : A single species phase species has been zeroed out
|
||||
* in this routine. The species is a noncomponent
|
||||
* - 2 : Same as one but, the zeroed species is a component.
|
||||
* - 3 : Nothing was done because the phase couldn't be birthed
|
||||
* because a needed component is zero.
|
||||
*/
|
||||
int VCS_SOLVE::vcs_popPhaseRxnStepSizes(const int iphasePop) {
|
||||
vcs_VolPhase *Vphase = m_VolPhaseList[iphasePop];
|
||||
// Identify the first species in the phase
|
||||
int kspec = Vphase->spGlobalIndexVCS(0);
|
||||
// Identify the formation reaction for that species
|
||||
int irxn = kspec - m_numComponents;
|
||||
|
||||
doublereal s;
|
||||
int j, k;
|
||||
// Calculate the initial moles of the phase being born.
|
||||
// Here we set it to 10x of the value which would cause the phase to be
|
||||
// zeroed out within the algorithm. We may later adjust the value.
|
||||
doublereal tPhaseMoles = 10. * m_totalMolNum * VCS_DELETE_PHASE_CUTOFF;
|
||||
|
||||
|
||||
#ifdef DEBUG_MODE
|
||||
int existence = Vphase->exists();
|
||||
if (existence > 0) {
|
||||
printf("ERROR vcs_popPhaseRxnStepSizes called for a phase that exists!");
|
||||
exit(-1);
|
||||
}
|
||||
char anote[256];
|
||||
if (m_debug_print_lvl >= 2) {
|
||||
plogf(" --- vcs_popPhaseRxnStepSizes() called to pop phase %s %d into existence\n",
|
||||
Vphase->PhaseName.c_str(), iphasePop);
|
||||
}
|
||||
#endif
|
||||
// Section for a single-species phase
|
||||
//
|
||||
if (Vphase->m_singleSpecies) {
|
||||
s = 0.0;
|
||||
double *dnPhase_irxn = m_deltaMolNumPhase[irxn];
|
||||
for (j = 0; j < m_numComponents; ++j) {
|
||||
if (!m_SSPhase[j]) {
|
||||
if (m_molNumSpecies_old[j] > 0.0) {
|
||||
s += SQUARE(m_stoichCoeffRxnMatrix[irxn][j]) / m_molNumSpecies_old[j];
|
||||
}
|
||||
}
|
||||
}
|
||||
for (j = 0; j < m_numPhases; j++) {
|
||||
Vphase = m_VolPhaseList[j];
|
||||
if (! Vphase->m_singleSpecies) {
|
||||
if (m_tPhaseMoles_old[j] > 0.0)
|
||||
s -= SQUARE(dnPhase_irxn[j]) / m_tPhaseMoles_old[j];
|
||||
}
|
||||
}
|
||||
if (s != 0.0) {
|
||||
double s_old = s;
|
||||
s = vcs_Hessian_diag_adj(irxn, s_old);
|
||||
#ifdef DEBUG_MODE
|
||||
if (s_old != s) {
|
||||
sprintf(anote, "Normal calc: diag adjusted from %g "
|
||||
"to %g due to act coeff", s_old, s);
|
||||
}
|
||||
#endif
|
||||
m_deltaMolNumSpecies[kspec] = -m_deltaGRxn_new[irxn] / s;
|
||||
} else {
|
||||
// Ok, s is equal to zero. We can not apply a sophisticated theory
|
||||
// to birth the phase. Just pick a small delta and go with it.
|
||||
m_deltaMolNumSpecies[kspec] = tPhaseMoles;
|
||||
}
|
||||
|
||||
/*
|
||||
* section to do damping of the m_deltaMolNumSpecies[]
|
||||
*/
|
||||
for (j = 0; j < m_numComponents; ++j) {
|
||||
double stoicC = m_stoichCoeffRxnMatrix[irxn][j];
|
||||
if (stoicC != 0.0) {
|
||||
if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) {
|
||||
double negChangeComp = - stoicC * m_deltaMolNumSpecies[kspec];
|
||||
if (negChangeComp > m_molNumSpecies_old[j]) {
|
||||
if (m_molNumSpecies_old[j] > 0.0) {
|
||||
#ifdef DEBUG_MODE
|
||||
sprintf(anote, "Delta damped from %g "
|
||||
"to %g due to component %d (%10s) going neg", m_deltaMolNumSpecies[kspec],
|
||||
-m_molNumSpecies_old[j]/stoicC, j, m_speciesName[j].c_str());
|
||||
#endif
|
||||
m_deltaMolNumSpecies[kspec] = - 0.5 * m_molNumSpecies_old[j] / stoicC;
|
||||
} else {
|
||||
#ifdef DEBUG_MODE
|
||||
sprintf(anote, "Delta damped from %g "
|
||||
"to %g due to component %d (%10s) zero", m_deltaMolNumSpecies[kspec],
|
||||
-m_molNumSpecies_old[j]/stoicC, j, m_speciesName[j].c_str());
|
||||
#endif
|
||||
m_deltaMolNumSpecies[kspec] = 0.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
// Implement a damping term that limits m_deltaMolNumSpecies to the size of the mole number
|
||||
if (-m_deltaMolNumSpecies[kspec] > m_molNumSpecies_old[kspec]) {
|
||||
#ifdef DEBUG_MODE
|
||||
sprintf(anote, "Delta damped from %g "
|
||||
"to %g due to %s going negative", m_deltaMolNumSpecies[kspec],
|
||||
-m_molNumSpecies_old[kspec], m_speciesName[kspec].c_str());
|
||||
#endif
|
||||
m_deltaMolNumSpecies[kspec] = -m_molNumSpecies_old[kspec];
|
||||
}
|
||||
|
||||
|
||||
} else {
|
||||
vector<doublereal> fracDelta(Vphase->nSpecies());
|
||||
vector<doublereal> X_est(Vphase->nSpecies());
|
||||
fracDelta = Vphase->fractionCreationDeltas();
|
||||
|
||||
double sumFrac = 0.0;
|
||||
for (k = 0; k < Vphase->nSpecies(); k++) {
|
||||
sumFrac += fracDelta[k];
|
||||
}
|
||||
for (k = 0; k < Vphase->nSpecies(); k++) {
|
||||
X_est[k] = fracDelta[k] / sumFrac;
|
||||
}
|
||||
|
||||
doublereal deltaMolNumPhase = tPhaseMoles;
|
||||
doublereal damp = 1.0;
|
||||
m_deltaGRxn_tmp = m_molNumSpecies_old;
|
||||
double * molNumSpecies_tmp = DATA_PTR(m_deltaGRxn_tmp);
|
||||
|
||||
|
||||
for (k = 0; k < Vphase->nSpecies(); k++) {
|
||||
kspec = Vphase->spGlobalIndexVCS(k);
|
||||
double delmol = deltaMolNumPhase * X_est[k];
|
||||
irxn = kspec - m_numComponents;
|
||||
if (kspec >= m_numComponents) {
|
||||
for (j = 0; j < m_numComponents; ++j) {
|
||||
double stoicC = m_stoichCoeffRxnMatrix[irxn][j];
|
||||
if (stoicC != 0.0) {
|
||||
if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) {
|
||||
molNumSpecies_tmp[j] += stoicC * delmol;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
doublereal ratioComp = 0.0;
|
||||
for (j = 0; j < m_numComponents; ++j) {
|
||||
double deltaJ = m_molNumSpecies_old[j] - molNumSpecies_tmp[j];
|
||||
if (molNumSpecies_tmp[j] < 0.0) {
|
||||
ratioComp = 1.0;
|
||||
if (deltaJ > 0.0) {
|
||||
double delta0 = m_molNumSpecies_old[j];
|
||||
double dampj = delta0 / deltaJ * 0.9;
|
||||
if (dampj < damp) {
|
||||
damp = dampj;
|
||||
}
|
||||
}
|
||||
} else {
|
||||
if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) {
|
||||
int jph = m_phaseID[j];
|
||||
if ((jph != iphasePop) && (!m_SSPhase[j])) {
|
||||
double fdeltaJ = fabs(deltaJ);
|
||||
if ( m_molNumSpecies_old[j] > 0.0) {
|
||||
ratioComp = MAX(ratioComp, fdeltaJ/ m_molNumSpecies_old[j]);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// We may have greatly underestimated the deltaMoles for the phase pop
|
||||
// Here we create a damp > 1 to account for this possibility.
|
||||
// We adjust upwards to make sure that a component in an existing multispecies
|
||||
// phase is modified by a factor of 1/1000.
|
||||
if (ratioComp > 1.0E-30) {
|
||||
if (ratioComp < 0.001) {
|
||||
damp = 0.001 / ratioComp;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
if (damp <= 1.0E-6) {
|
||||
return 3;
|
||||
}
|
||||
|
||||
for (k = 0; k < Vphase->nSpecies(); k++) {
|
||||
kspec = Vphase->spGlobalIndexVCS(k);
|
||||
if (kspec < m_numComponents) {
|
||||
m_speciesStatus[kspec] = VCS_SPECIES_COMPONENT;
|
||||
} else {
|
||||
m_deltaMolNumSpecies[kspec] = deltaMolNumPhase * X_est[k] * damp;
|
||||
if (X_est[k] > 1.0E-3) {
|
||||
m_speciesStatus[kspec] = VCS_SPECIES_MAJOR;
|
||||
} else {
|
||||
m_speciesStatus[kspec] = VCS_SPECIES_MINOR;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
|
||||
//
|
||||
|
||||
double VCS_SOLVE::vcs_phaseStabilityTest(const int iph) {
|
||||
|
||||
/*
|
||||
* We will use the _new state calc here
|
||||
*/
|
||||
int kspec, irxn, k;
|
||||
int kspec, irxn, k, i, kc, kc_spec;
|
||||
vcs_VolPhase *Vphase = m_VolPhaseList[iph];
|
||||
double deltaGRxn;
|
||||
doublereal deltaGRxn;
|
||||
|
||||
// We will do a full newton calculation later, but for now, ...
|
||||
bool doSuccessiveSubstitution = true;
|
||||
int res = 0;
|
||||
double funcPhaseStability;
|
||||
vector<doublereal> X_est(Vphase->nSpecies(), 0.0);
|
||||
vector<doublereal> delFrac(Vphase->nSpecies(), 0.0);
|
||||
vector<doublereal> E_phi(Vphase->nSpecies(), 0.0);
|
||||
vector<doublereal> fracDelta_new(Vphase->nSpecies(), 0.0);
|
||||
vector<doublereal> fracDelta_old(Vphase->nSpecies(), 0.0);
|
||||
vector<doublereal> fracDelta_raw(Vphase->nSpecies(), 0.0);
|
||||
|
||||
vector<double> X_est(Vphase->nSpecies(), 0.0);
|
||||
vector<double> X_est_old(Vphase->nSpecies(), 0.0);
|
||||
vector<double> delX(Vphase->nSpecies(), 0.0);
|
||||
vector<double> E_phi(Vphase->nSpecies(), 0.0);
|
||||
double damp = 1.0;
|
||||
double normUpdate = 1.0;
|
||||
double normUpdateOld = 1.0;
|
||||
vector<doublereal> m_feSpecies_Deficient(m_numComponents, 0.0);
|
||||
doublereal damp = 1.0;
|
||||
doublereal dampOld = 1.0;
|
||||
doublereal normUpdate = 1.0;
|
||||
doublereal normUpdateOld = 1.0;
|
||||
doublereal sum = 0.0;
|
||||
doublereal dirProd = 0.0;
|
||||
doublereal dirProdOld = 0.0;
|
||||
|
||||
// get the activity coefficients
|
||||
Vphase->sendToVCS_ActCoeff(VCS_STATECALC_OLD, VCS_DATA_PTR(m_actCoeffSpecies_new));
|
||||
|
||||
|
||||
|
||||
// Get the storred estimate for the composition of the phase if
|
||||
// it gets created
|
||||
fracDelta_new = Vphase->fractionCreationDeltas();
|
||||
|
||||
#ifdef DEBUG_MODE
|
||||
if (m_temperature < 380.) {
|
||||
// fracDelta_new[0] = 0.8;
|
||||
// fracDelta_new[1] = 0.1;
|
||||
// fracDelta_new[2] = 1.0E-8;
|
||||
// fracDelta_new[3] = 0.1;
|
||||
//fracDelta_new[4] = 1.0E-8;
|
||||
}
|
||||
if (m_temperature < 390.) {
|
||||
printf("we are here\n");
|
||||
}
|
||||
#endif
|
||||
|
||||
bool oneIsComponent = false;
|
||||
std::vector<int> componentList;
|
||||
|
||||
for (k = 0; k < Vphase->nSpecies(); k++) {
|
||||
kspec = Vphase->spGlobalIndexVCS(k);
|
||||
if (kspec < m_numComponents) {
|
||||
oneIsComponent = true;
|
||||
componentList.push_back(k);
|
||||
}
|
||||
}
|
||||
|
||||
for (k = 0; k < m_numComponents; k++) {
|
||||
m_feSpecies_Deficient[k] = m_feSpecies_old[k];
|
||||
}
|
||||
normUpdate = 0.1 * vcs_l2norm(fracDelta_new);
|
||||
damp = 1.0E-2;
|
||||
|
||||
if (doSuccessiveSubstitution) {
|
||||
|
||||
for (int its = 0; its < 20; its++) {
|
||||
|
||||
|
||||
bool converged = false;
|
||||
for (int its = 0; its < 200 && (!converged); its++) {
|
||||
|
||||
dampOld = damp;
|
||||
normUpdateOld = normUpdate;
|
||||
for (k = 0; k < Vphase->nSpecies(); k++) {
|
||||
X_est_old[k] = X_est[k];
|
||||
fracDelta_old = fracDelta_new;
|
||||
dirProdOld = dirProd;
|
||||
|
||||
// Given a set of fracDelta's, we calculate the fracDelta's
|
||||
// for the component species, if any
|
||||
for (i = 0; i < (int) componentList.size(); i++) {
|
||||
kc = componentList[i];
|
||||
kc_spec = Vphase->spGlobalIndexVCS(kc);
|
||||
fracDelta_old[kc] = 0.0;
|
||||
for (k = 0; k < Vphase->nSpecies(); k++) {
|
||||
kspec = Vphase->spGlobalIndexVCS(k);
|
||||
irxn = kspec - m_numComponents;
|
||||
if (irxn >= 0) {
|
||||
fracDelta_old[kc] += m_stoichCoeffRxnMatrix[irxn][kc_spec] * fracDelta_old[k];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
double poly = -1.0;
|
||||
// Now, calculate the predicted mole fractions, X_est[k]
|
||||
double sumFrac = 0.0;
|
||||
for (k = 0; k < Vphase->nSpecies(); k++) {
|
||||
kspec = Vphase->spGlobalIndexVCS(k);
|
||||
irxn = kspec - m_numComponents;
|
||||
deltaGRxn = m_deltaGRxn_old[irxn];
|
||||
// We may need to look at deltaGRxn for components!
|
||||
if (irxn >= 0) {
|
||||
if (deltaGRxn > 50.0) deltaGRxn = 50.0;
|
||||
if (deltaGRxn < -50.0) deltaGRxn = -50.0;
|
||||
E_phi[k] = exp(-deltaGRxn)/m_actCoeffSpecies_new[kspec];
|
||||
poly += E_phi[k];
|
||||
sumFrac += fracDelta_old[k];
|
||||
}
|
||||
double sum_Xcomp = 0.0;
|
||||
for (k = 0; k < Vphase->nSpecies(); k++) {
|
||||
X_est[k] = fracDelta_old[k] / sumFrac;
|
||||
kc_spec = Vphase->spGlobalIndexVCS(k);
|
||||
if (kc_spec < m_numComponents) {
|
||||
sum_Xcomp += X_est[k];
|
||||
}
|
||||
}
|
||||
double sum = poly + 1.0;
|
||||
|
||||
for (k = 0; k < Vphase->nSpecies(); k++) {
|
||||
delX[k] = E_phi[k]/sum - X_est_old[k];
|
||||
}
|
||||
normUpdate = vcs_l2norm(delX);
|
||||
|
||||
// Figure out the damping coefficient
|
||||
double ratio = normUpdate / normUpdateOld;
|
||||
if (ratio < 0.4) {
|
||||
damp = 1.0;
|
||||
} else if (ratio > 1.0) {
|
||||
damp = 0.03;
|
||||
} else {
|
||||
damp = 0.1;
|
||||
}
|
||||
|
||||
|
||||
|
||||
/*
|
||||
* Feed the newly formed estimate of the mole fractions back into the
|
||||
* ThermoPhase object
|
||||
*/
|
||||
Vphase->setMoleFractionsState(0.0, VCS_DATA_PTR(X_est), VCS_STATECALC_PHASESTABILITY);
|
||||
|
||||
/*
|
||||
* get the activity coefficients
|
||||
*/
|
||||
Vphase->sendToVCS_ActCoeff(VCS_STATECALC_OLD, VCS_DATA_PTR(m_actCoeffSpecies_new));
|
||||
|
||||
/*
|
||||
* first Calculate altered chemical potentials for component species
|
||||
* belonging to this phase.
|
||||
*/
|
||||
for (i = 0; i < (int) componentList.size(); i++) {
|
||||
kc = componentList[i];
|
||||
kc_spec = Vphase->spGlobalIndexVCS(kc);
|
||||
if ( X_est[kc] > VCS_DELETE_MINORSPECIES_CUTOFF) {
|
||||
m_feSpecies_Deficient[kc_spec] = m_feSpecies_old[kc_spec]
|
||||
+ log(m_actCoeffSpecies_new[kc_spec] * X_est[kc]);
|
||||
} else {
|
||||
m_feSpecies_Deficient[kc_spec] = m_feSpecies_old[kc_spec]
|
||||
+ log(m_actCoeffSpecies_new[kc_spec] * VCS_DELETE_MINORSPECIES_CUTOFF);
|
||||
}
|
||||
}
|
||||
|
||||
for (i = 0; i < (int) componentList.size(); i++) {
|
||||
kc = componentList[i];
|
||||
kc_spec = Vphase->spGlobalIndexVCS(kc);
|
||||
|
||||
for (k = 0; k < Vphase->nSpecies(); k++) {
|
||||
kspec = Vphase->spGlobalIndexVCS(k);
|
||||
irxn = kspec - m_numComponents;
|
||||
if (irxn >= 0) {
|
||||
if (i == 0) {
|
||||
m_deltaGRxn_Deficient[irxn] = m_deltaGRxn_old[irxn];
|
||||
}
|
||||
double *dtmp_ptr = m_stoichCoeffRxnMatrix[irxn];
|
||||
if (dtmp_ptr[kc_spec] != 0.0) {
|
||||
m_deltaGRxn_Deficient[irxn] +=
|
||||
dtmp_ptr[kc_spec] * (m_feSpecies_Deficient[kc_spec]- m_feSpecies_old[kc_spec]);
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
/*
|
||||
* Calculate the E_phi's
|
||||
*/
|
||||
sum = 0.0;
|
||||
funcPhaseStability = sum_Xcomp - 1.0;
|
||||
for (k = 0; k < Vphase->nSpecies(); k++) {
|
||||
X_est[k] = X_est_old[k] + damp * delX[k];
|
||||
kspec = Vphase->spGlobalIndexVCS(k);
|
||||
irxn = kspec - m_numComponents;
|
||||
if (irxn >= 0) {
|
||||
deltaGRxn = m_deltaGRxn_Deficient[irxn];
|
||||
if (deltaGRxn > 50.0) deltaGRxn = 50.0;
|
||||
if (deltaGRxn < -50.0) deltaGRxn = -50.0;
|
||||
E_phi[k] = std::exp(-deltaGRxn) / m_actCoeffSpecies_new[kspec];
|
||||
sum += E_phi[k];
|
||||
funcPhaseStability += E_phi[k];
|
||||
} else {
|
||||
E_phi[k] = 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
/*
|
||||
* Calculate the raw estimate of the new fracs
|
||||
*/
|
||||
for (k = 0; k < Vphase->nSpecies(); k++) {
|
||||
kspec = Vphase->spGlobalIndexVCS(k);
|
||||
irxn = kspec - m_numComponents;
|
||||
double b = E_phi[k] / sum * (1.0 - sum_Xcomp);
|
||||
if (irxn >= 0) {
|
||||
fracDelta_raw[k] = (sumFrac - fracDelta_old[k]) * b / (1.0 - b);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Given a set of fracDelta's, we calculate the fracDelta's
|
||||
// for the component species, if any
|
||||
for (i = 0; i < (int) componentList.size(); i++) {
|
||||
kc = componentList[i];
|
||||
kc_spec = Vphase->spGlobalIndexVCS(kc);
|
||||
fracDelta_raw[kc] = 0.0;
|
||||
for (k = 0; k < Vphase->nSpecies(); k++) {
|
||||
kspec = Vphase->spGlobalIndexVCS(k);
|
||||
irxn = kspec - m_numComponents;
|
||||
if (irxn >= 0) {
|
||||
fracDelta_raw[kc] += m_stoichCoeffRxnMatrix[irxn][kc_spec] * fracDelta_raw[k];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
/*
|
||||
* Now possibly dampen the estimate.
|
||||
*/
|
||||
doublereal sumADel = 0.0;
|
||||
for (k = 0; k < Vphase->nSpecies(); k++) {
|
||||
delFrac[k] = fracDelta_raw[k] - fracDelta_old[k];
|
||||
sumADel += fabs(delFrac[k]);
|
||||
}
|
||||
normUpdate = vcs_l2norm(delFrac);
|
||||
|
||||
dirProd = 0.0;
|
||||
for (k = 0; k < Vphase->nSpecies(); k++) {
|
||||
dirProd += fracDelta_old[k] * delFrac[k];
|
||||
}
|
||||
bool crossedSign = false;
|
||||
if (dirProd * dirProdOld < 0.0) {
|
||||
crossedSign = true;
|
||||
}
|
||||
|
||||
|
||||
damp = 0.5;
|
||||
if (dampOld < 0.25) {
|
||||
damp = 2.0 * dampOld;
|
||||
}
|
||||
if (crossedSign) {
|
||||
if (normUpdate *1.5 > normUpdateOld) {
|
||||
damp = 0.5 * dampOld;
|
||||
} else if (normUpdate *2.0 > normUpdateOld) {
|
||||
damp = 0.8 * dampOld;
|
||||
}
|
||||
} else {
|
||||
if (normUpdate > normUpdateOld * 2.0) {
|
||||
damp = 0.6 * dampOld;
|
||||
} else if (normUpdate > normUpdateOld * 1.2) {
|
||||
damp = 0.9 * dampOld;
|
||||
}
|
||||
}
|
||||
|
||||
for (k = 0; k < Vphase->nSpecies(); k++) {
|
||||
kspec = Vphase->spGlobalIndexVCS(k);
|
||||
m_molNumSpecies_new[kspec] = X_est[k];
|
||||
if (fabs(damp * delFrac[k]) > 0.3*fabs(fracDelta_old[k])) {
|
||||
damp = MAX(0.3*fabs(fracDelta_old[k]) / fabs( delFrac[k]), 1.0E-8/fabs( delFrac[k]));
|
||||
}
|
||||
if (delFrac[k] < 0.0) {
|
||||
if (2.0 * damp * (-delFrac[k]) > fracDelta_old[k]) {
|
||||
damp = fracDelta_old[k] / (2.0 * (-delFrac[k]));
|
||||
}
|
||||
}
|
||||
if (delFrac[k] > 0.0) {
|
||||
if (2.0 * damp * delFrac[k] > fracDelta_old[k]) {
|
||||
damp = fracDelta_old[k] / (2.0 * delFrac[k]);
|
||||
}
|
||||
}
|
||||
}
|
||||
if (damp < 0.000001) {
|
||||
damp = 0.000001;
|
||||
}
|
||||
|
||||
for (k = 0; k < Vphase->nSpecies(); k++) {
|
||||
fracDelta_new[k] = fracDelta_old[k] + damp * (delFrac[k]);
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
if (normUpdate < 1.0E-5) {
|
||||
converged = true;
|
||||
}
|
||||
Vphase->setMolesFromVCS(VCS_STATECALC_NEW);
|
||||
|
||||
|
||||
}
|
||||
|
||||
if (converged) {
|
||||
Vphase->setMoleFractionsState(0.0, VCS_DATA_PTR(X_est), VCS_STATECALC_PHASESTABILITY);
|
||||
Vphase->setFractionCreationDeltas( VCS_DATA_PTR(fracDelta_new));
|
||||
}
|
||||
|
||||
|
||||
} else {
|
||||
printf("not done yet\n");
|
||||
exit(-1);
|
||||
}
|
||||
|
||||
return res;
|
||||
return funcPhaseStability;
|
||||
}
|
||||
|
||||
}
|
||||
|
|
|
|||
|
|
@ -27,10 +27,12 @@ namespace VCSnonideal {
|
|||
*
|
||||
* Output
|
||||
* -------
|
||||
* m_deltaMolNumSpecies(irxn) : reaction adjustments, where irxn refers
|
||||
* m_deltaMolNumSpecies[kspec] : reaction adjustments, where irxn refers
|
||||
* to the irxn'th species
|
||||
* formation reaction. This adjustment is for species
|
||||
* irxn + M, where M is the number of components.
|
||||
* formation reaction. This adjustment
|
||||
* is for species
|
||||
* irxn + M, where M is the number
|
||||
* of components.
|
||||
*
|
||||
* Special branching occurs sometimes. This causes the component basis
|
||||
* to be reevaluated
|
||||
|
|
@ -127,7 +129,7 @@ namespace VCSnonideal {
|
|||
}
|
||||
} else {
|
||||
/********************************************************************/
|
||||
/************************* REGULAR PROCESSING ************/
|
||||
/************************* REGULAR PROCESSING ***********************/
|
||||
/********************************************************************/
|
||||
/*
|
||||
* First take care of cases where we want to bail out
|
||||
|
|
|
|||
|
|
@ -498,7 +498,37 @@ public:
|
|||
*
|
||||
*/
|
||||
void vcs_updateVP(const int stateCalc);
|
||||
|
||||
|
||||
|
||||
//! Utility function that evaluates whether a phase can be popped
|
||||
//! into existence
|
||||
/*!
|
||||
* @param iphasePop id of the phase, which is currently zeroed,
|
||||
*
|
||||
* @return Returns true if the phase can come into existence
|
||||
* and false otherwise.
|
||||
*/
|
||||
bool vcs_popPhasePossible(const int iphasePop) const;
|
||||
|
||||
//! Decision as to whether a phase pops back into existence
|
||||
/*!
|
||||
* @return returns the phase id of the phase that pops back into
|
||||
* existence. Returns -1 if there are no phases
|
||||
*/
|
||||
int vcs_popPhaseID();
|
||||
|
||||
//! Calculates the deltas of the reactions due to phases popping
|
||||
//! into existence
|
||||
/*!
|
||||
* @param iphasePop Phase id of the phase that will come into existence
|
||||
*
|
||||
* @return Returns an int representing the status of the step
|
||||
* - 0 : normal return
|
||||
* - 1 : A single species phase species has been zeroed out
|
||||
* in this routine. The species is a noncomponent
|
||||
* - 2 : Same as one but, the zeroed species is a component.
|
||||
*/
|
||||
int vcs_popPhaseRxnStepSizes(const int iphasePop);
|
||||
|
||||
//! Calculates formation reaction step sizes.
|
||||
/*!
|
||||
|
|
@ -633,8 +663,13 @@ public:
|
|||
double vcs_birthGuess(const int kspec);
|
||||
|
||||
|
||||
|
||||
int vcs_phaseStabilityTest(const int iph);
|
||||
//! Main program to test whether a deleted phase should be brought
|
||||
//! back into existence
|
||||
/*!
|
||||
*
|
||||
* @param iph Phase id of the deleted phase
|
||||
*/
|
||||
double vcs_phaseStabilityTest(const int iph);
|
||||
|
||||
//! Solve an equilibrium problem at a particular fixed temperature
|
||||
//! and pressure
|
||||
|
|
@ -1441,7 +1476,7 @@ public:
|
|||
* stoichiometric coefficient of one is assumed for the
|
||||
* species K in this mechanism.
|
||||
*
|
||||
* NOTE: K = IRXN + NC
|
||||
* NOTE: kspec = Irxn + m_numComponents
|
||||
*
|
||||
* sc[irxn][j] :
|
||||
* j refers to the component number, and irxn
|
||||
|
|
|
|||
|
|
@ -447,45 +447,74 @@ namespace VCSnonideal {
|
|||
vcs_dzero(VCS_DATA_PTR(m_deltaMolNumSpecies), m_numSpeciesTot);
|
||||
|
||||
/*
|
||||
* Figure out whether we will calculate new reaction step sizes
|
||||
* for the major species.
|
||||
* -> We won't if all species are minors (im), OR
|
||||
* all major species have already converged
|
||||
* First step is a major branch in the algorithm.
|
||||
* We first determine if a phase pops into existence.
|
||||
*/
|
||||
if (!(MajorSpeciesHaveConverged) && ! allMinorZeroedSpecies) {
|
||||
soldel = vcs_RxnStepSizes();
|
||||
/* - If SOLDEL is true then we encountered a reaction between */
|
||||
/* - single-species-phase species, only, and have adjusted */
|
||||
/* - the mole number vector, W(), directly. In this case, */
|
||||
/* - we should immediately go back and recompute a new */
|
||||
/* - component basis, if the species that was zeroed was */
|
||||
/* - a component. SOLDEL is true when this is so. */
|
||||
if (soldel > 0) {
|
||||
/* - We have changed the base mole number amongst single- */
|
||||
/* - species-phase species. However, we don't need to */
|
||||
/* - recaculate their chemical potentials because they */
|
||||
/* - are constant, anyway! */
|
||||
if (soldel == 2) {
|
||||
goto L_COMPONENT_CALC;
|
||||
}
|
||||
/* - We have not changed the actual DG values for */
|
||||
/* - any species, even the one we deleted. Thus, */
|
||||
/* - we don't need to start over. */
|
||||
}
|
||||
} else {
|
||||
int iphasePop = vcs_popPhaseID();
|
||||
/*
|
||||
*
|
||||
*/
|
||||
soldel = -1;
|
||||
if (iphasePop >= 0) {
|
||||
soldel = vcs_popPhaseRxnStepSizes(iphasePop);
|
||||
if (soldel == 3) {
|
||||
iphasePop = -1;
|
||||
#ifdef DEBUG_MODE
|
||||
if (m_debug_print_lvl >= 2) {
|
||||
if (allMinorZeroedSpecies) {
|
||||
plogf(" --- vcs_RxnStepSizes not called because all"
|
||||
"species are minors\n");
|
||||
} else {
|
||||
plogf(" --- vcs_RxnStepSizes not called because "
|
||||
"all majors have converged\n");
|
||||
if (m_debug_print_lvl >= 2) {
|
||||
plogf(" --- vcs_popPhaseRxnStepSizes() was called but stoich prevented phase %d popping\n");
|
||||
}
|
||||
}
|
||||
#endif
|
||||
}
|
||||
}
|
||||
|
||||
if (iphasePop < 0) {
|
||||
/*
|
||||
* Figure out whether we will calculate new reaction step sizes
|
||||
* for the major species.
|
||||
* -> We won't if all species are minors (im), OR
|
||||
* all major species have already converged
|
||||
*/
|
||||
if (!(MajorSpeciesHaveConverged) && ! allMinorZeroedSpecies) {
|
||||
soldel = vcs_RxnStepSizes();
|
||||
/* - If SOLDEL is true then we encountered a reaction between */
|
||||
/* - single-species-phase species, only, and have adjusted */
|
||||
/* - the mole number vector, W(), directly. In this case, */
|
||||
/* - we should immediately go back and recompute a new */
|
||||
/* - component basis, if the species that was zeroed was */
|
||||
/* - a component. SOLDEL is true when this is so. */
|
||||
if (soldel > 0) {
|
||||
/* - We have changed the base mole number amongst single- */
|
||||
/* - species-phase species. However, we don't need to */
|
||||
/* - recaculate their chemical potentials because they */
|
||||
/* - are constant, anyway! */
|
||||
if (soldel == 2) {
|
||||
goto L_COMPONENT_CALC;
|
||||
}
|
||||
/* - We have not changed the actual DG values for */
|
||||
/* - any species, even the one we deleted. Thus, */
|
||||
/* - we don't need to start over. */
|
||||
}
|
||||
} else {
|
||||
#ifdef DEBUG_MODE
|
||||
if (m_debug_print_lvl >= 2) {
|
||||
if (allMinorZeroedSpecies) {
|
||||
plogf(" --- vcs_RxnStepSizes not called because all"
|
||||
"species are minors\n");
|
||||
} else {
|
||||
plogf(" --- vcs_RxnStepSizes not called because "
|
||||
"all majors have converged\n");
|
||||
}
|
||||
}
|
||||
#endif
|
||||
}
|
||||
}
|
||||
#ifdef DEBUG_MODE
|
||||
else {
|
||||
if (m_debug_print_lvl >= 2) {
|
||||
plogf(" --- vcs_RxnStepSizes not called because alternative"
|
||||
"phase creation delta was used instead\n");
|
||||
}
|
||||
}
|
||||
#endif
|
||||
lec = FALSE;
|
||||
doPhaseDeleteIph = -1;
|
||||
doPhaseDeleteKspec = -1;
|
||||
|
|
@ -544,389 +573,403 @@ namespace VCSnonideal {
|
|||
#ifdef DEBUG_MODE
|
||||
ANOTE[0] = '\0';
|
||||
#endif
|
||||
|
||||
if (m_speciesStatus[kspec] == VCS_SPECIES_INTERFACIALVOLTAGE) {
|
||||
/********************************************************************/
|
||||
/************************ VOLTAGE SPECIES ***************************/
|
||||
/********************************************************************/
|
||||
#ifdef DEBUG_MODE
|
||||
dx = vcs_minor_alt_calc(kspec, irxn, &soldel, ANOTE);
|
||||
#else
|
||||
dx = vcs_minor_alt_calc(kspec, irxn, &soldel);
|
||||
#endif
|
||||
m_deltaMolNumSpecies[kspec] = dx;
|
||||
}
|
||||
else if (m_speciesStatus[kspec] < VCS_SPECIES_MINOR) {
|
||||
/********************************************************************/
|
||||
/********************** ZEROED OUT SPECIES **************************/
|
||||
/********************************************************************/
|
||||
bool resurrect = true;
|
||||
#ifdef DEBUG_MODE
|
||||
if (m_debug_print_lvl >= 3) {
|
||||
plogf(" --- %s currently zeroed (SpStatus=%-2d):",
|
||||
m_speciesName[kspec].c_str(), m_speciesStatus[kspec]);
|
||||
plogf("%3d DG = %11.4E WT = %11.4E W = %11.4E DS = %11.4E\n",
|
||||
irxn, m_deltaGRxn_new[irxn], m_molNumSpecies_new[kspec],
|
||||
m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec]);
|
||||
}
|
||||
#endif
|
||||
// HKM Alternative is to not allow ds[] = 0.0 phases
|
||||
// to pop back into existence. For esthetics, I'm allowing this.
|
||||
// so that dg < 0.0 phases with zero mole numbers become components.
|
||||
// This is also better, because that component will be the first
|
||||
// one to pop into existence if there is a minute quantity of the element.
|
||||
// This could change in the future.
|
||||
//if (dg[irxn] >= 0.0 || ds[kspec] <= 0.0) {
|
||||
if (m_deltaGRxn_new[irxn] >= 0.0 ) {
|
||||
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec];
|
||||
m_deltaMolNumSpecies[kspec] = 0.0;
|
||||
resurrect = false;
|
||||
#ifdef DEBUG_MODE
|
||||
sprintf(ANOTE, "Species stays zeroed: DG = %11.4E",
|
||||
m_deltaGRxn_new[irxn]);
|
||||
if (m_deltaGRxn_new[irxn] < 0.0) {
|
||||
sprintf(ANOTE, "Species stays zeroed even though dg neg:DG = %11.4E, ds zeroed ",
|
||||
m_deltaGRxn_new[irxn]);
|
||||
}
|
||||
#endif
|
||||
if (iphasePop >= 0) {
|
||||
if (iph == iphasePop) {
|
||||
dx = m_deltaMolNumSpecies[kspec];
|
||||
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + m_deltaMolNumSpecies[kspec];
|
||||
#ifdef DEBUG_MODE
|
||||
sprintf(ANOTE, "Phase pop");
|
||||
#endif
|
||||
} else {
|
||||
for (int j = 0; j < m_numElemConstraints; ++j) {
|
||||
int elType = m_elType[j];
|
||||
if (elType == VCS_ELEM_TYPE_ABSPOS) {
|
||||
double atomComp = m_formulaMatrix[j][kspec];
|
||||
if (atomComp > 0.0) {
|
||||
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",
|
||||
m_elementName[j].c_str());
|
||||
#endif
|
||||
resurrect = false;
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
/*
|
||||
* Resurrect the species
|
||||
*/
|
||||
if (resurrect) {
|
||||
bool phaseResurrected = false;
|
||||
if (Vphase->exists() == VCS_PHASE_EXIST_NO) {
|
||||
//Vphase->setExistence(1);
|
||||
phaseResurrected = true;
|
||||
}
|
||||
--m_numRxnMinorZeroed;
|
||||
|
||||
if (phaseResurrected) {
|
||||
#ifdef DEBUG_MODE
|
||||
if (m_debug_print_lvl >= 2) {
|
||||
plogf(" --- Zeroed species changed to major: ");
|
||||
plogf("%-12s\n", m_speciesName[kspec].c_str());
|
||||
}
|
||||
#endif
|
||||
m_speciesStatus[kspec] = VCS_SPECIES_MAJOR;
|
||||
MajorSpeciesHaveConverged = false;
|
||||
allMinorZeroedSpecies = false;
|
||||
} else {
|
||||
#ifdef DEBUG_MODE
|
||||
if (m_debug_print_lvl >= 2) {
|
||||
plogf(" --- Zeroed species changed to minor: ");
|
||||
plogf("%-12s\n", m_speciesName[kspec].c_str());
|
||||
}
|
||||
#endif
|
||||
m_speciesStatus[kspec] = VCS_SPECIES_MINOR;
|
||||
}
|
||||
if (m_deltaMolNumSpecies[kspec] > 0.0) {
|
||||
dx = m_deltaMolNumSpecies[kspec] * 0.01;
|
||||
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx;
|
||||
} else {
|
||||
m_molNumSpecies_new[kspec] = m_totalMolNum * VCS_DELETE_PHASE_CUTOFF * 10.;
|
||||
dx = m_molNumSpecies_new[kspec] - m_molNumSpecies_old[kspec];
|
||||
}
|
||||
m_deltaMolNumSpecies[kspec] = dx;
|
||||
#ifdef DEBUG_MODE
|
||||
sprintf(ANOTE, "Born:IC=-1 to IC=1:DG=%11.4E", m_deltaGRxn_new[irxn]);
|
||||
#endif
|
||||
} else {
|
||||
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec];
|
||||
m_deltaMolNumSpecies[kspec] = 0.0;
|
||||
dx = 0.0;
|
||||
}
|
||||
} else if (m_speciesStatus[kspec] == VCS_SPECIES_MINOR) {
|
||||
/********************************************************************/
|
||||
/***************************** MINOR SPECIES ************************/
|
||||
/********************************************************************/
|
||||
/*
|
||||
* Unless ITI isn't equal to zero we zero out changes
|
||||
* to minor species.
|
||||
*/
|
||||
if (iti != 0) {
|
||||
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec];
|
||||
m_deltaMolNumSpecies[kspec] = 0.0;
|
||||
dx = 0.0;
|
||||
#ifdef DEBUG_MODE
|
||||
sprintf(ANOTE,"minor species not considered");
|
||||
if (m_debug_print_lvl >= 2) {
|
||||
plogf(" --- "); plogf("%-12s", m_speciesName[kspec].c_str());
|
||||
plogf("%3d%11.4E%11.4E%11.4E | %s",
|
||||
m_speciesStatus[kspec], m_molNumSpecies_old[kspec], m_molNumSpecies_new[kspec],
|
||||
m_deltaMolNumSpecies[kspec], ANOTE);
|
||||
plogendl();
|
||||
}
|
||||
#endif
|
||||
continue;
|
||||
}
|
||||
/*
|
||||
* Minor species alternative calculation
|
||||
* ---------------------------------------
|
||||
* This is based upon the following approximation:
|
||||
* The mole fraction changes due to these reactions don't affect
|
||||
* the mole numbers of the component species. Therefore the
|
||||
* following approximation is valid for an ideal solution
|
||||
* 0 = DG(I) + log(WT(I)/W(I))
|
||||
* (DG contains the contribution from FF(I) + log(W(I)/TL) )
|
||||
* Thus,
|
||||
* WT(I) = W(I) EXP(-DG(I))
|
||||
* If soldel is true on return, then we branch to the section
|
||||
* that deletes a species from the current set of active species.
|
||||
*/
|
||||
#ifdef DEBUG_MODE
|
||||
dx = vcs_minor_alt_calc(kspec, irxn, &soldel, ANOTE);
|
||||
#else
|
||||
dx = vcs_minor_alt_calc(kspec, irxn, &soldel);
|
||||
#endif
|
||||
m_deltaMolNumSpecies[kspec] = dx;
|
||||
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx;
|
||||
|
||||
if (soldel) {
|
||||
/*******************************************************************/
|
||||
/***** DELETE MINOR SPECIES LESS THAN VCS_DELETE_SPECIES_CUTOFF */
|
||||
/***** MOLE NUMBER */
|
||||
/*******************************************************************/
|
||||
#ifdef DEBUG_MODE
|
||||
if (m_debug_print_lvl >= 2) {
|
||||
plogf(" --- Delete minor species in multispec phase: %-12s",
|
||||
m_speciesName[kspec].c_str());
|
||||
plogendl();
|
||||
}
|
||||
#endif
|
||||
m_deltaMolNumSpecies[kspec] = 0.0;
|
||||
/*
|
||||
* Delete species, kspec. The alternate return is for the case
|
||||
* where all species become deleted. Then, we need to
|
||||
* branch to the code where we reevaluate the deletion
|
||||
* of all species.
|
||||
*/
|
||||
lnospec = vcs_delete_species(kspec);
|
||||
if (lnospec) goto L_RECHECK_DELETED;
|
||||
/*
|
||||
* Go back to consider the next species in the list.
|
||||
* Note, however, that the next species in the list is now
|
||||
* in slot l. In deleting the previous species L, We have
|
||||
* exchanged slot MR with slot l, and then have
|
||||
* decremented MR.
|
||||
* Therefore, we will decrement the species counter, here.
|
||||
*/
|
||||
--irxn;
|
||||
#ifdef DEBUG_MODE
|
||||
goto L_MAIN_LOOP_END_NO_PRINT;
|
||||
#else
|
||||
goto L_MAIN_LOOP_END;
|
||||
#endif
|
||||
}
|
||||
} else {
|
||||
/********************************************************************/
|
||||
/*********************** MAJOR SPECIES ******************************/
|
||||
/********************************************************************/
|
||||
#ifdef DEBUG_MODE
|
||||
sprintf(ANOTE, "Normal Major Calc");
|
||||
|
||||
|
||||
if (m_speciesStatus[kspec] == VCS_SPECIES_INTERFACIALVOLTAGE) {
|
||||
/********************************************************************/
|
||||
/************************ VOLTAGE SPECIES ***************************/
|
||||
/********************************************************************/
|
||||
#ifdef DEBUG_MODE
|
||||
dx = vcs_minor_alt_calc(kspec, irxn, &soldel, ANOTE);
|
||||
#else
|
||||
dx = vcs_minor_alt_calc(kspec, irxn, &soldel);
|
||||
#endif
|
||||
/*
|
||||
* Check for superconvergence of the formation reaction. Do
|
||||
* nothing if it is superconverged. Skip to the end of the
|
||||
* irxn loop if it is superconverged.
|
||||
*/
|
||||
if (fabs(m_deltaGRxn_new[irxn]) <= m_tolmaj2) {
|
||||
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec];
|
||||
m_deltaMolNumSpecies[kspec] = 0.0;
|
||||
dx = 0.0;
|
||||
m_deltaMolNumSpecies[kspec] = dx;
|
||||
}
|
||||
else if (m_speciesStatus[kspec] < VCS_SPECIES_MINOR) {
|
||||
/********************************************************************/
|
||||
/********************** ZEROED OUT SPECIES **************************/
|
||||
/********************************************************************/
|
||||
bool resurrect = true;
|
||||
#ifdef DEBUG_MODE
|
||||
sprintf(ANOTE, "major species is converged");
|
||||
if (m_debug_print_lvl >= 2) {
|
||||
plogf(" --- "); plogf("%-12s", m_speciesName[kspec].c_str());
|
||||
plogf("%3d%11.4E%11.4E%11.4E | %s",
|
||||
m_speciesStatus[kspec], m_molNumSpecies_old[kspec], m_molNumSpecies_new[kspec],
|
||||
m_deltaMolNumSpecies[kspec], ANOTE);
|
||||
plogendl();
|
||||
if (m_debug_print_lvl >= 3) {
|
||||
plogf(" --- %s currently zeroed (SpStatus=%-2d):",
|
||||
m_speciesName[kspec].c_str(), m_speciesStatus[kspec]);
|
||||
plogf("%3d DG = %11.4E WT = %11.4E W = %11.4E DS = %11.4E\n",
|
||||
irxn, m_deltaGRxn_new[irxn], m_molNumSpecies_new[kspec],
|
||||
m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec]);
|
||||
}
|
||||
#endif
|
||||
continue;
|
||||
}
|
||||
/*
|
||||
* Set the initial step size, dx, equal to the value produced
|
||||
* by the routine, vcs_RxnStepSize().
|
||||
*
|
||||
* Note the multiplition logic is to make sure that
|
||||
* dg[] didn't change sign due to w[] changing in the
|
||||
* middle of the iteration. (it can if a single species
|
||||
* phase goes out of existence).
|
||||
*/
|
||||
if ((m_deltaGRxn_new[irxn] * m_deltaMolNumSpecies[kspec]) <= 0.0) {
|
||||
dx = m_deltaMolNumSpecies[kspec];
|
||||
} else {
|
||||
dx = 0.0;
|
||||
m_deltaMolNumSpecies[kspec] = 0.0;
|
||||
// HKM Alternative is to not allow ds[] = 0.0 phases
|
||||
// to pop back into existence. For esthetics, I'm allowing this.
|
||||
// so that dg < 0.0 phases with zero mole numbers become components.
|
||||
// This is also better, because that component will be the first
|
||||
// one to pop into existence if there is a minute quantity of the element.
|
||||
// This could change in the future.
|
||||
//if (dg[irxn] >= 0.0 || ds[kspec] <= 0.0) {
|
||||
if (m_deltaGRxn_new[irxn] >= 0.0 ) {
|
||||
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec];
|
||||
m_deltaMolNumSpecies[kspec] = 0.0;
|
||||
resurrect = false;
|
||||
#ifdef DEBUG_MODE
|
||||
sprintf(ANOTE, "dx set to 0, DG flipped sign due to "
|
||||
"changed initial point");
|
||||
sprintf(ANOTE, "Species stays zeroed: DG = %11.4E",
|
||||
m_deltaGRxn_new[irxn]);
|
||||
if (m_deltaGRxn_new[irxn] < 0.0) {
|
||||
sprintf(ANOTE, "Species stays zeroed even though dg neg:DG = %11.4E, ds zeroed ",
|
||||
m_deltaGRxn_new[irxn]);
|
||||
}
|
||||
#endif
|
||||
}
|
||||
/*
|
||||
* Form a tentative value of the new species moles
|
||||
*/
|
||||
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx;
|
||||
|
||||
/*
|
||||
* Check for non-positive mole fraction of major species.
|
||||
* If we find one, we branch to a section below. Then,
|
||||
* depending upon the outcome, we branch to sections below,
|
||||
* or we restart the entire iteration.
|
||||
*/
|
||||
if (m_molNumSpecies_new[kspec] <= 0.0) {
|
||||
#ifdef DEBUG_MODE
|
||||
sprintf(ANOTE, "initial nonpos kmoles= %11.3E",
|
||||
m_molNumSpecies_new[kspec]);
|
||||
#endif
|
||||
/* ************************************************* */
|
||||
/* *** NON-POSITIVE MOLES OF MAJOR SPECIES ********* */
|
||||
/* ************************************************* */
|
||||
/*
|
||||
* We are here when a tentative value of a mole fraction
|
||||
* created by a tentative value of M_DELTAMOLNUMSPECIES(*) is negative.
|
||||
* We branch from here depending upon whether this
|
||||
* species is in a single species phase or in
|
||||
* a multispecies phase.
|
||||
*/
|
||||
if (! (m_SSPhase[kspec])) {
|
||||
/*
|
||||
* Section for multispecies phases:
|
||||
* - Cut reaction adjustment for positive kmoles of
|
||||
* major species in multispecies phases.
|
||||
* Decrease its concentration by a factor of 10.
|
||||
*/
|
||||
dx = -0.9 * m_molNumSpecies_old[kspec];
|
||||
m_deltaMolNumSpecies[kspec] = dx;
|
||||
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx;
|
||||
} else {
|
||||
/*
|
||||
* Section for single species phases:
|
||||
* Calculate a dx that will wipe out the
|
||||
* moles in the phase.
|
||||
*/
|
||||
dx = -m_molNumSpecies_old[kspec];
|
||||
/*
|
||||
* Calculate an update that doesn't create a negative mole
|
||||
* number for a component species. Actually, restrict this
|
||||
* a little more so that the component values can only be
|
||||
* reduced by two 99%,
|
||||
*/
|
||||
for (j = 0; j < m_numComponents; ++j) {
|
||||
if (sc_irxn[j] != 0.0) {
|
||||
wx[j] = m_molNumSpecies_old[j] + sc_irxn[j] * dx;
|
||||
if (wx[j] <= m_molNumSpecies_old[j] * 0.01 - 1.0E-150) {
|
||||
dx = MAX(dx, m_molNumSpecies_old[j] * -0.99 / sc_irxn[j]);
|
||||
for (int j = 0; j < m_numElemConstraints; ++j) {
|
||||
int elType = m_elType[j];
|
||||
if (elType == VCS_ELEM_TYPE_ABSPOS) {
|
||||
double atomComp = m_formulaMatrix[j][kspec];
|
||||
if (atomComp > 0.0) {
|
||||
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",
|
||||
m_elementName[j].c_str());
|
||||
#endif
|
||||
resurrect = false;
|
||||
break;
|
||||
}
|
||||
}
|
||||
} else {
|
||||
wx[j] = m_molNumSpecies_old[j];
|
||||
}
|
||||
}
|
||||
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx;
|
||||
if (m_molNumSpecies_new[kspec] > 0.0) {
|
||||
m_deltaMolNumSpecies[kspec] = dx;
|
||||
#ifdef DEBUG_MODE
|
||||
sprintf(ANOTE,
|
||||
"zeroing SS phase created a neg component species "
|
||||
"-> reducing step size instead");
|
||||
#endif
|
||||
} else {
|
||||
/*
|
||||
* We are going to zero the single species phase.
|
||||
* Set the existence flag
|
||||
*/
|
||||
iph = m_phaseID[kspec];
|
||||
Vphase = m_VolPhaseList[iph];
|
||||
//Vphase->setExistence(0);
|
||||
#ifdef DEBUG_MODE
|
||||
sprintf(ANOTE, "zeroing out SS phase: ");
|
||||
#endif
|
||||
/*
|
||||
* Change the base mole numbers for the iteration.
|
||||
* We need to do this here, because we have decided
|
||||
* to eliminate the phase in this special section
|
||||
* outside the main loop.
|
||||
*/
|
||||
m_molNumSpecies_new[kspec] = 0.0;
|
||||
doPhaseDeleteIph = iph;
|
||||
doPhaseDeleteKspec = kspec;
|
||||
}
|
||||
/*
|
||||
* Resurrect the species
|
||||
*/
|
||||
if (resurrect) {
|
||||
bool phaseResurrected = false;
|
||||
if (Vphase->exists() == VCS_PHASE_EXIST_NO) {
|
||||
//Vphase->setExistence(1);
|
||||
phaseResurrected = true;
|
||||
}
|
||||
--m_numRxnMinorZeroed;
|
||||
|
||||
if (phaseResurrected) {
|
||||
#ifdef DEBUG_MODE
|
||||
if (m_debug_print_lvl >= 2) {
|
||||
if (m_speciesStatus[kspec] >= 0) {
|
||||
plogf(" --- SS species changed to zeroedss: ");
|
||||
plogf("%-12s", m_speciesName[kspec].c_str());
|
||||
plogendl();
|
||||
}
|
||||
plogf(" --- Zeroed species changed to major: ");
|
||||
plogf("%-12s\n", m_speciesName[kspec].c_str());
|
||||
}
|
||||
#endif
|
||||
m_speciesStatus[kspec] = VCS_SPECIES_ZEROEDSS;
|
||||
++m_numRxnMinorZeroed;
|
||||
allMinorZeroedSpecies = (m_numRxnMinorZeroed == m_numRxnRdc);
|
||||
|
||||
for (int kk = 0; kk < m_numSpeciesTot; kk++) {
|
||||
m_deltaMolNumSpecies[kk] = 0.0;
|
||||
m_molNumSpecies_new[kk] = m_molNumSpecies_old[kk];
|
||||
m_speciesStatus[kspec] = VCS_SPECIES_MAJOR;
|
||||
MajorSpeciesHaveConverged = false;
|
||||
allMinorZeroedSpecies = false;
|
||||
} else {
|
||||
#ifdef DEBUG_MODE
|
||||
if (m_debug_print_lvl >= 2) {
|
||||
plogf(" --- Zeroed species changed to minor: ");
|
||||
plogf("%-12s\n", m_speciesName[kspec].c_str());
|
||||
}
|
||||
m_deltaMolNumSpecies[kspec] = dx;
|
||||
m_molNumSpecies_new[kspec] = 0.0;
|
||||
|
||||
for (k = 0; k < m_numComponents; ++k) {
|
||||
m_deltaMolNumSpecies[k] = 0.0;
|
||||
}
|
||||
for (iph = 0; iph < m_numPhases; iph++) {
|
||||
m_deltaPhaseMoles[iph] = 0.0;
|
||||
}
|
||||
|
||||
#endif
|
||||
m_speciesStatus[kspec] = VCS_SPECIES_MINOR;
|
||||
}
|
||||
if (m_deltaMolNumSpecies[kspec] > 0.0) {
|
||||
dx = m_deltaMolNumSpecies[kspec] * 0.01;
|
||||
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx;
|
||||
} else {
|
||||
m_molNumSpecies_new[kspec] = m_totalMolNum * VCS_DELETE_PHASE_CUTOFF * 10.;
|
||||
dx = m_molNumSpecies_new[kspec] - m_molNumSpecies_old[kspec];
|
||||
}
|
||||
m_deltaMolNumSpecies[kspec] = dx;
|
||||
#ifdef DEBUG_MODE
|
||||
sprintf(ANOTE, "Born:IC=-1 to IC=1:DG=%11.4E", m_deltaGRxn_new[irxn]);
|
||||
#endif
|
||||
} else {
|
||||
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec];
|
||||
m_deltaMolNumSpecies[kspec] = 0.0;
|
||||
dx = 0.0;
|
||||
}
|
||||
} else if (m_speciesStatus[kspec] == VCS_SPECIES_MINOR) {
|
||||
/********************************************************************/
|
||||
/***************************** MINOR SPECIES ************************/
|
||||
/********************************************************************/
|
||||
/*
|
||||
* Unless ITI isn't equal to zero we zero out changes
|
||||
* to minor species.
|
||||
*/
|
||||
if (iti != 0) {
|
||||
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec];
|
||||
m_deltaMolNumSpecies[kspec] = 0.0;
|
||||
dx = 0.0;
|
||||
#ifdef DEBUG_MODE
|
||||
sprintf(ANOTE,"minor species not considered");
|
||||
if (m_debug_print_lvl >= 2) {
|
||||
plogf(" --- "); plogf("%-12s", m_speciesName[kspec].c_str());
|
||||
plogf("%3d%11.4E%11.4E%11.4E | %s",
|
||||
m_speciesStatus[kspec], m_molNumSpecies_old[kspec], m_molNumSpecies_new[kspec],
|
||||
m_deltaMolNumSpecies[kspec], ANOTE);
|
||||
plogendl();
|
||||
}
|
||||
#endif
|
||||
continue;
|
||||
}
|
||||
/*
|
||||
* Minor species alternative calculation
|
||||
* ---------------------------------------
|
||||
* This is based upon the following approximation:
|
||||
* The mole fraction changes due to these reactions don't affect
|
||||
* the mole numbers of the component species. Therefore the
|
||||
* following approximation is valid for an ideal solution
|
||||
* 0 = DG(I) + log(WT(I)/W(I))
|
||||
* (DG contains the contribution from FF(I) + log(W(I)/TL) )
|
||||
* Thus,
|
||||
* WT(I) = W(I) EXP(-DG(I))
|
||||
* If soldel is true on return, then we branch to the section
|
||||
* that deletes a species from the current set of active species.
|
||||
*/
|
||||
#ifdef DEBUG_MODE
|
||||
dx = vcs_minor_alt_calc(kspec, irxn, &soldel, ANOTE);
|
||||
#else
|
||||
dx = vcs_minor_alt_calc(kspec, irxn, &soldel);
|
||||
#endif
|
||||
m_deltaMolNumSpecies[kspec] = dx;
|
||||
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx;
|
||||
|
||||
if (soldel) {
|
||||
/*******************************************************************/
|
||||
/***** DELETE MINOR SPECIES LESS THAN VCS_DELETE_SPECIES_CUTOFF */
|
||||
/***** MOLE NUMBER */
|
||||
/*******************************************************************/
|
||||
#ifdef DEBUG_MODE
|
||||
if (m_debug_print_lvl >= 2) {
|
||||
plogf(" --- Delete minor species in multispec phase: %-12s",
|
||||
m_speciesName[kspec].c_str());
|
||||
plogendl();
|
||||
}
|
||||
#endif
|
||||
m_deltaMolNumSpecies[kspec] = 0.0;
|
||||
/*
|
||||
* Delete species, kspec. The alternate return is for the case
|
||||
* where all species become deleted. Then, we need to
|
||||
* branch to the code where we reevaluate the deletion
|
||||
* of all species.
|
||||
*/
|
||||
lnospec = vcs_delete_species(kspec);
|
||||
if (lnospec) goto L_RECHECK_DELETED;
|
||||
/*
|
||||
* Go back to consider the next species in the list.
|
||||
* Note, however, that the next species in the list is now
|
||||
* in slot l. In deleting the previous species L, We have
|
||||
* exchanged slot MR with slot l, and then have
|
||||
* decremented MR.
|
||||
* Therefore, we will decrement the species counter, here.
|
||||
*/
|
||||
--irxn;
|
||||
#ifdef DEBUG_MODE
|
||||
goto L_MAIN_LOOP_END_NO_PRINT;
|
||||
#else
|
||||
goto L_MAIN_LOOP_END;
|
||||
#endif
|
||||
}
|
||||
} else {
|
||||
/********************************************************************/
|
||||
/*********************** MAJOR SPECIES ******************************/
|
||||
/********************************************************************/
|
||||
#ifdef DEBUG_MODE
|
||||
sprintf(ANOTE, "Normal Major Calc");
|
||||
#endif
|
||||
/*
|
||||
* Check for superconvergence of the formation reaction. Do
|
||||
* nothing if it is superconverged. Skip to the end of the
|
||||
* irxn loop if it is superconverged.
|
||||
*/
|
||||
if (fabs(m_deltaGRxn_new[irxn]) <= m_tolmaj2) {
|
||||
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec];
|
||||
m_deltaMolNumSpecies[kspec] = 0.0;
|
||||
dx = 0.0;
|
||||
#ifdef DEBUG_MODE
|
||||
sprintf(ANOTE, "major species is converged");
|
||||
if (m_debug_print_lvl >= 2) {
|
||||
plogf(" --- "); plogf("%-12s", m_speciesName[kspec].c_str());
|
||||
plogf("%3d%11.4E%11.4E%11.4E | %s",
|
||||
m_speciesStatus[kspec], m_molNumSpecies_old[kspec], m_molNumSpecies_new[kspec],
|
||||
m_deltaMolNumSpecies[kspec], ANOTE);
|
||||
plogendl();
|
||||
}
|
||||
#endif
|
||||
continue;
|
||||
}
|
||||
/*
|
||||
* Set the initial step size, dx, equal to the value produced
|
||||
* by the routine, vcs_RxnStepSize().
|
||||
*
|
||||
* Note the multiplition logic is to make sure that
|
||||
* dg[] didn't change sign due to w[] changing in the
|
||||
* middle of the iteration. (it can if a single species
|
||||
* phase goes out of existence).
|
||||
*/
|
||||
if ((m_deltaGRxn_new[irxn] * m_deltaMolNumSpecies[kspec]) <= 0.0) {
|
||||
dx = m_deltaMolNumSpecies[kspec];
|
||||
} else {
|
||||
dx = 0.0;
|
||||
m_deltaMolNumSpecies[kspec] = 0.0;
|
||||
#ifdef DEBUG_MODE
|
||||
sprintf(ANOTE, "dx set to 0, DG flipped sign due to "
|
||||
"changed initial point");
|
||||
#endif
|
||||
}
|
||||
/*
|
||||
* Form a tentative value of the new species moles
|
||||
*/
|
||||
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx;
|
||||
|
||||
/*
|
||||
* Check for non-positive mole fraction of major species.
|
||||
* If we find one, we branch to a section below. Then,
|
||||
* depending upon the outcome, we branch to sections below,
|
||||
* or we restart the entire iteration.
|
||||
*/
|
||||
if (m_molNumSpecies_new[kspec] <= 0.0) {
|
||||
#ifdef DEBUG_MODE
|
||||
sprintf(ANOTE, "initial nonpos kmoles= %11.3E",
|
||||
m_molNumSpecies_new[kspec]);
|
||||
#endif
|
||||
/* ************************************************* */
|
||||
/* *** NON-POSITIVE MOLES OF MAJOR SPECIES ********* */
|
||||
/* ************************************************* */
|
||||
/*
|
||||
* We are here when a tentative value of a mole fraction
|
||||
* created by a tentative value of M_DELTAMOLNUMSPECIES(*) is negative.
|
||||
* We branch from here depending upon whether this
|
||||
* species is in a single species phase or in
|
||||
* a multispecies phase.
|
||||
*/
|
||||
if (! (m_SSPhase[kspec])) {
|
||||
/*
|
||||
* Section for multispecies phases:
|
||||
* - Cut reaction adjustment for positive kmoles of
|
||||
* major species in multispecies phases.
|
||||
* Decrease its concentration by a factor of 10.
|
||||
*/
|
||||
dx = -0.9 * m_molNumSpecies_old[kspec];
|
||||
m_deltaMolNumSpecies[kspec] = dx;
|
||||
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx;
|
||||
} else {
|
||||
/*
|
||||
* Section for single species phases:
|
||||
* Calculate a dx that will wipe out the
|
||||
* moles in the phase.
|
||||
*/
|
||||
dx = -m_molNumSpecies_old[kspec];
|
||||
/*
|
||||
* Calculate an update that doesn't create a negative mole
|
||||
* number for a component species. Actually, restrict this
|
||||
* a little more so that the component values can only be
|
||||
* reduced by two 99%,
|
||||
*/
|
||||
for (j = 0; j < m_numComponents; ++j) {
|
||||
if (sc_irxn[j] != 0.0) {
|
||||
wx[j] = m_molNumSpecies_old[j] + sc_irxn[j] * dx;
|
||||
if (wx[j] <= m_molNumSpecies_old[j] * 0.01 - 1.0E-150) {
|
||||
dx = MAX(dx, m_molNumSpecies_old[j] * -0.99 / sc_irxn[j]);
|
||||
}
|
||||
} else {
|
||||
wx[j] = m_molNumSpecies_old[j];
|
||||
}
|
||||
}
|
||||
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx;
|
||||
if (m_molNumSpecies_new[kspec] > 0.0) {
|
||||
m_deltaMolNumSpecies[kspec] = dx;
|
||||
#ifdef DEBUG_MODE
|
||||
sprintf(ANOTE,
|
||||
"zeroing SS phase created a neg component species "
|
||||
"-> reducing step size instead");
|
||||
#endif
|
||||
} else {
|
||||
/*
|
||||
* We are going to zero the single species phase.
|
||||
* Set the existence flag
|
||||
*/
|
||||
iph = m_phaseID[kspec];
|
||||
Vphase = m_VolPhaseList[iph];
|
||||
//Vphase->setExistence(0);
|
||||
#ifdef DEBUG_MODE
|
||||
sprintf(ANOTE, "zeroing out SS phase: ");
|
||||
#endif
|
||||
/*
|
||||
* Change the base mole numbers for the iteration.
|
||||
* We need to do this here, because we have decided
|
||||
* to eliminate the phase in this special section
|
||||
* outside the main loop.
|
||||
*/
|
||||
m_molNumSpecies_new[kspec] = 0.0;
|
||||
doPhaseDeleteIph = iph;
|
||||
doPhaseDeleteKspec = kspec;
|
||||
|
||||
#ifdef DEBUG_MODE
|
||||
if (m_debug_print_lvl >= 2) {
|
||||
if (m_speciesStatus[kspec] >= 0) {
|
||||
plogf(" --- SS species changed to zeroedss: ");
|
||||
plogf("%-12s", m_speciesName[kspec].c_str());
|
||||
plogendl();
|
||||
}
|
||||
}
|
||||
#endif
|
||||
m_speciesStatus[kspec] = VCS_SPECIES_ZEROEDSS;
|
||||
++m_numRxnMinorZeroed;
|
||||
allMinorZeroedSpecies = (m_numRxnMinorZeroed == m_numRxnRdc);
|
||||
|
||||
for (int kk = 0; kk < m_numSpeciesTot; kk++) {
|
||||
m_deltaMolNumSpecies[kk] = 0.0;
|
||||
m_molNumSpecies_new[kk] = m_molNumSpecies_old[kk];
|
||||
}
|
||||
m_deltaMolNumSpecies[kspec] = dx;
|
||||
m_molNumSpecies_new[kspec] = 0.0;
|
||||
|
||||
for (k = 0; k < m_numComponents; ++k) {
|
||||
m_deltaMolNumSpecies[k] = 0.0;
|
||||
}
|
||||
for (iph = 0; iph < m_numPhases; iph++) {
|
||||
m_deltaPhaseMoles[iph] = 0.0;
|
||||
}
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
#ifdef VCS_LINE_SEARCH
|
||||
/*********************************************************************/
|
||||
/*** LINE SEARCH ALGORITHM FOR MAJOR SPECIES IN NON-IDEAL PHASES *****/
|
||||
/*********************************************************************/
|
||||
/*
|
||||
* Skip the line search if we are birthing a species
|
||||
*/
|
||||
if ((dx != 0.0) &&
|
||||
(m_molNumSpecies_old[kspec] > 0.0) &&
|
||||
(doPhaseDeleteIph == -1) &&
|
||||
(m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE)) {
|
||||
double dx_old = dx;
|
||||
/*********************************************************************/
|
||||
/*** LINE SEARCH ALGORITHM FOR MAJOR SPECIES IN NON-IDEAL PHASES *****/
|
||||
/*********************************************************************/
|
||||
/*
|
||||
* Skip the line search if we are birthing a species
|
||||
*/
|
||||
if ((dx != 0.0) &&
|
||||
(m_molNumSpecies_old[kspec] > 0.0) &&
|
||||
(doPhaseDeleteIph == -1) &&
|
||||
(m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE)) {
|
||||
double dx_old = dx;
|
||||
|
||||
#ifdef DEBUG_MODE
|
||||
dx = vcs_line_search(irxn, dx_old, ANOTE);
|
||||
dx = vcs_line_search(irxn, dx_old, ANOTE);
|
||||
#else
|
||||
dx = vcs_line_search(irxn, dx_old);
|
||||
dx = vcs_line_search(irxn, dx_old);
|
||||
#endif
|
||||
vcs_setFlagsVolPhases(false, VCS_STATECALC_NEW);
|
||||
}
|
||||
m_deltaMolNumSpecies[kspec] = dx;
|
||||
vcs_setFlagsVolPhases(false, VCS_STATECALC_NEW);
|
||||
}
|
||||
m_deltaMolNumSpecies[kspec] = dx;
|
||||
#endif
|
||||
|
||||
} /* End of Loop on ic[irxn] -> the type of species */
|
||||
}/* End of Loop on ic[irxn] -> the type of species */
|
||||
}
|
||||
/***********************************************************************/
|
||||
/****** CALCULATE KMOLE NUMBER CHANGE FOR THE COMPONENT BASIS **********/
|
||||
/***********************************************************************/
|
||||
|
|
@ -4748,7 +4791,8 @@ namespace VCSnonideal {
|
|||
* This should be implemented.
|
||||
*/
|
||||
int k;
|
||||
if (alterZeroedPhases) {
|
||||
//alterZeroedPhases = false;
|
||||
if (alterZeroedPhases && false) {
|
||||
for (iph = 0; iph < m_numPhases; iph++) {
|
||||
lneed = FALSE;
|
||||
vcs_VolPhase *Vphase = m_VolPhaseList[iph];
|
||||
|
|
|
|||
|
|
@ -75,8 +75,10 @@ namespace Cantera {
|
|||
MC_apCut_(0.0),
|
||||
MC_bpCut_(0.0),
|
||||
MC_cpCut_(0.0),
|
||||
CROP_ln_gamma_o_min(-25.0),
|
||||
CROP_ln_gamma_k_max(23.0),
|
||||
CROP_ln_gamma_o_min(-10.0),
|
||||
CROP_ln_gamma_o_max(3.0),
|
||||
CROP_ln_gamma_k_min(-5.0),
|
||||
CROP_ln_gamma_k_max(15.0),
|
||||
m_debugCalc(0)
|
||||
{
|
||||
for (int i = 0; i < 17; i++) {
|
||||
|
|
@ -130,8 +132,10 @@ namespace Cantera {
|
|||
MC_apCut_(0.0),
|
||||
MC_bpCut_(0.0),
|
||||
MC_cpCut_(0.0),
|
||||
CROP_ln_gamma_o_min(-25.0),
|
||||
CROP_ln_gamma_k_max(23.0),
|
||||
CROP_ln_gamma_o_min(-10.0),
|
||||
CROP_ln_gamma_o_max(3.0),
|
||||
CROP_ln_gamma_k_min(-5.0),
|
||||
CROP_ln_gamma_k_max(15.0),
|
||||
m_debugCalc(0)
|
||||
{
|
||||
for (int i = 0; i < 17; i++) {
|
||||
|
|
@ -179,8 +183,10 @@ namespace Cantera {
|
|||
MC_apCut_(0.0),
|
||||
MC_bpCut_(0.0),
|
||||
MC_cpCut_(0.0),
|
||||
CROP_ln_gamma_o_min(-25.0),
|
||||
CROP_ln_gamma_k_max(23.0),
|
||||
CROP_ln_gamma_o_min(-10.0),
|
||||
CROP_ln_gamma_o_max(3.0),
|
||||
CROP_ln_gamma_k_min(-5.0),
|
||||
CROP_ln_gamma_k_max(15.0),
|
||||
m_debugCalc(0)
|
||||
{
|
||||
for (int i = 0; i < 17; i++) {
|
||||
|
|
@ -234,8 +240,10 @@ namespace Cantera {
|
|||
MC_apCut_(0.0),
|
||||
MC_bpCut_(0.0),
|
||||
MC_cpCut_(0.0),
|
||||
CROP_ln_gamma_o_min(-25.0),
|
||||
CROP_ln_gamma_k_max(23.0),
|
||||
CROP_ln_gamma_o_min(-10.0),
|
||||
CROP_ln_gamma_o_max(3.0),
|
||||
CROP_ln_gamma_k_min(-5.0),
|
||||
CROP_ln_gamma_k_max(15.0),
|
||||
m_debugCalc(0)
|
||||
{
|
||||
/*
|
||||
|
|
@ -389,6 +397,8 @@ namespace Cantera {
|
|||
MC_bpCut_ = b.MC_bpCut_;
|
||||
MC_cpCut_ = b.MC_cpCut_;
|
||||
CROP_ln_gamma_o_min = b.CROP_ln_gamma_o_min;
|
||||
CROP_ln_gamma_o_max = b.CROP_ln_gamma_o_max;
|
||||
CROP_ln_gamma_k_min = b.CROP_ln_gamma_k_min;
|
||||
CROP_ln_gamma_k_max = b.CROP_ln_gamma_k_max;
|
||||
m_CounterIJ = b.m_CounterIJ;
|
||||
m_molalitiesCropped = b.m_molalitiesCropped;
|
||||
|
|
@ -463,8 +473,10 @@ namespace Cantera {
|
|||
MC_apCut_(0.0),
|
||||
MC_bpCut_(0.0),
|
||||
MC_cpCut_(0.0),
|
||||
CROP_ln_gamma_o_min(-25.0),
|
||||
CROP_ln_gamma_k_max(23.0),
|
||||
CROP_ln_gamma_o_min(-10.0),
|
||||
CROP_ln_gamma_o_max(3.0),
|
||||
CROP_ln_gamma_k_min(-5.0),
|
||||
CROP_ln_gamma_k_max(15.0),
|
||||
m_debugCalc(0)
|
||||
{
|
||||
if (testProb != 1) {
|
||||
|
|
@ -1812,14 +1824,22 @@ namespace Cantera {
|
|||
|
||||
for (int k = 1; k < m_kk; k++) {
|
||||
m_lnActCoeffMolal_Unscaled[k] += IMS_lnActCoeffMolal_[k];
|
||||
if (m_lnActCoeffMolal_Unscaled[k] > (CROP_ln_gamma_k_max + lnxs)) {
|
||||
m_lnActCoeffMolal_Unscaled[k] = CROP_ln_gamma_k_max + lnxs;
|
||||
if (m_lnActCoeffMolal_Unscaled[k] > (CROP_ln_gamma_k_max)) {
|
||||
m_lnActCoeffMolal_Unscaled[k] = CROP_ln_gamma_k_max;
|
||||
}
|
||||
if (m_lnActCoeffMolal_Unscaled[k] < (CROP_ln_gamma_k_min - 2.5 *lnxs)) {
|
||||
// -1.0 and -1.5 caused multiple solutions
|
||||
m_lnActCoeffMolal_Unscaled[k] = CROP_ln_gamma_k_min - 2.5 * lnxs;
|
||||
}
|
||||
}
|
||||
|
||||
m_lnActCoeffMolal_Unscaled[0] += (IMS_lnActCoeffMolal_[0] - lnActCoeffMolal0);
|
||||
if (m_lnActCoeffMolal_Unscaled[0] < CROP_ln_gamma_o_min - lnxs) {
|
||||
m_lnActCoeffMolal_Unscaled[0] = CROP_ln_gamma_o_min - lnxs;
|
||||
if (m_lnActCoeffMolal_Unscaled[0] < CROP_ln_gamma_o_min) {
|
||||
m_lnActCoeffMolal_Unscaled[0] = CROP_ln_gamma_o_min;
|
||||
}
|
||||
if (m_lnActCoeffMolal_Unscaled[0] > CROP_ln_gamma_o_max) {
|
||||
// -0.5 caused multiple solutions
|
||||
m_lnActCoeffMolal_Unscaled[0] = CROP_ln_gamma_o_max;
|
||||
}
|
||||
|
||||
/*
|
||||
|
|
|
|||
|
|
@ -3194,9 +3194,12 @@ namespace Cantera {
|
|||
doublereal MC_cpCut_;
|
||||
|
||||
doublereal CROP_ln_gamma_o_min;
|
||||
doublereal CROP_ln_gamma_o_max;
|
||||
|
||||
doublereal CROP_ln_gamma_k_min;
|
||||
doublereal CROP_ln_gamma_k_max;
|
||||
|
||||
|
||||
//! Local error routine
|
||||
/*!
|
||||
* @param msg print out a message and error exit
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue