diff --git a/ext/SConscript b/ext/SConscript index 2a4832c91..95e638bb8 100644 --- a/ext/SConscript +++ b/ext/SConscript @@ -26,8 +26,7 @@ def prep_f2c(env): if env['VERBOSE'] : print "INFO 2: prep_f2c: adding MSDOS to CPPDEFINES" - if env['VERBOSE'] : - localenv.Append(CPPDEFINES=['added']) + if env['VERBOSE'] : print "INFO 2: CPPDEFINES ", localenv['CPPDEFINES'] #print "localenv", localenv.Dump() #exit(0) diff --git a/include/cantera/equil/vcs_VolPhase.h b/include/cantera/equil/vcs_VolPhase.h index 3b7be7ad4..9b75e4e87 100644 --- a/include/cantera/equil/vcs_VolPhase.h +++ b/include/cantera/equil/vcs_VolPhase.h @@ -916,13 +916,13 @@ private: mutable std::vector ActCoeff; //! Vector of the derivatives of the ln activity coefficient wrt to the - //! current mole number + //! current mole number multiplied by the current phase moles /*! - * dLnActCoeffdMolNumber[j][k]; + * np_dLnActCoeffdMolNumber[j][k]; * j = id of the species mole number * k = id of the species activity coefficient */ - mutable DoubleStarStar dLnActCoeffdMolNumber; + mutable DoubleStarStar np_dLnActCoeffdMolNumber; //! Status /*! diff --git a/include/cantera/equil/vcs_internal.h b/include/cantera/equil/vcs_internal.h index da6ccbd08..4d0c83f4e 100644 --- a/include/cantera/equil/vcs_internal.h +++ b/include/cantera/equil/vcs_internal.h @@ -521,6 +521,23 @@ void vcs_print_stringTrunc(const char* str, size_t space, int alignment); */ bool vcs_doubleEqual(double d1, double d2); + +//! Sorts a vector of ints in place from lowest to the highest values +/*! + * The vector is returned sorted from lowest to highest. + * + * @param x Reference to a vector of ints. + */ +void vcs_heapsort(std::vector &x); + +//! Sorts a vector of ints and eliminates duplicates from the resulting list +/*! + * @param xOrderedUnique Ordered vector of unique ints that were part of the original list + * @param x Reference to a constant vector of ints. + */ +void vcs_orderedUnique(std::vector & xOrderedUnique, const std::vector & x); + + } #endif diff --git a/include/cantera/equil/vcs_solve.h b/include/cantera/equil/vcs_solve.h index 9bbeb1def..cfd931b1c 100644 --- a/include/cantera/equil/vcs_solve.h +++ b/include/cantera/equil/vcs_solve.h @@ -1958,16 +1958,15 @@ public: */ std::vector m_actCoeffSpecies_old; - //! Change in activity coefficient with mole number + //! Change in the log of the activity coefficient with respect to the mole number + //! multiplied by the phase mole number /*! * length = [nspecies][nspecies] * - * (This is a temporary array that - * gets regenerated every time it's - * needed. It is not swapped wrt species - * (unused atm) + * This is a temporary array that gets regenerated every time it's + * needed. It is not swapped wrt species. */ - DoubleStarStar m_dLnActCoeffdMolNum; + DoubleStarStar m_np_dLnActCoeffdMolNum; //! Molecular weight of each species /*! diff --git a/include/cantera/kinetics/InterfaceKinetics.h b/include/cantera/kinetics/InterfaceKinetics.h index 6335ad1d7..fe51c45c4 100644 --- a/include/cantera/kinetics/InterfaceKinetics.h +++ b/include/cantera/kinetics/InterfaceKinetics.h @@ -556,7 +556,7 @@ public: * @param iphase Index of the phase. This is the order within the internal thermo vector object * @param exists Boolean indicating whether the phase exists or not */ - void setPhaseExistence(const size_t iphase, const bool exists); + void setPhaseExistence(const size_t iphase, const int exists); //! Set the stability of a phase in the reaction object /*! diff --git a/include/cantera/thermo/Phase.h b/include/cantera/thermo/Phase.h index d64b7fda5..aacfbf01d 100644 --- a/include/cantera/thermo/Phase.h +++ b/include/cantera/thermo/Phase.h @@ -543,9 +543,8 @@ public: //! Note the density of a phase is an independent variable. //! @param[in] density density (kg/m^3). virtual void setDensity(const doublereal density_) { - if (density_ <= 0) { - throw CanteraError("Phase::setDensity", - "density must be positive"); + if (density_ <= 0.0) { + throw CanteraError("Phase::setDensity()", "density must be positive"); } m_dens = density_; } diff --git a/src/equil/vcs_VolPhase.cpp b/src/equil/vcs_VolPhase.cpp index 480835d67..ac98b79e0 100644 --- a/src/equil/vcs_VolPhase.cpp +++ b/src/equil/vcs_VolPhase.cpp @@ -213,7 +213,7 @@ vcs_VolPhase& vcs_VolPhase::operator=(const vcs_VolPhase& b) StarMolarVol = b.StarMolarVol; PartialMolarVol = b.PartialMolarVol; ActCoeff = b.ActCoeff; - dLnActCoeffdMolNumber = b.dLnActCoeffdMolNumber; + np_dLnActCoeffdMolNumber = b.np_dLnActCoeffdMolNumber; m_vcsStateStatus = b.m_vcsStateStatus; m_phi = b.m_phi; m_UpToDate = false; @@ -311,7 +311,7 @@ void vcs_VolPhase::resize(const size_t phaseNum, const size_t nspecies, StarMolarVol.resize(nspecies, -1.0); PartialMolarVol.resize(nspecies, -1.0); ActCoeff.resize(nspecies, 1.0); - dLnActCoeffdMolNumber.resize(nspecies, nspecies, 0.0); + np_dLnActCoeffdMolNumber.resize(nspecies, nspecies, 0.0); m_speciesUnknownType.resize(nspecies, VCS_SPECIES_TYPE_MOLNUM); @@ -914,6 +914,11 @@ double vcs_VolPhase::_updateVolPM() const void vcs_VolPhase::_updateLnActCoeffJac() { + double phaseTotalMoles = v_totalMoles; + if (phaseTotalMoles < 1.0E-14) { + phaseTotalMoles = 1.0; + } + /* * Evaluate the current base activity coefficients if necessary */ @@ -923,15 +928,15 @@ void vcs_VolPhase::_updateLnActCoeffJac() if (!TP_ptr) { return; } - TP_ptr->getdlnActCoeffdlnN(m_numSpecies, &dLnActCoeffdMolNumber[0][0]); + TP_ptr->getdlnActCoeffdlnN(m_numSpecies, &np_dLnActCoeffdMolNumber[0][0]); for (size_t j = 0; j < m_numSpecies; j++) { - double moles_j_base = v_totalMoles * Xmol_[j]; - double* const lnActCoeffCol = dLnActCoeffdMolNumber[j]; + double moles_j_base = phaseTotalMoles * Xmol_[j]; + double* const np_lnActCoeffCol = np_dLnActCoeffdMolNumber[j]; if (moles_j_base < 1.0E-200) { - moles_j_base = 1.0E-7 * moles_j_base + 1.0E-20 * v_totalMoles + 1.0E-150; + moles_j_base = 1.0E-7 * moles_j_base + 1.0E-13 * phaseTotalMoles + 1.0E-150; } for (size_t k = 0; k < m_numSpecies; k++) { - lnActCoeffCol[k] /= moles_j_base; + np_lnActCoeffCol[k] = np_lnActCoeffCol[k] * phaseTotalMoles / moles_j_base; } } @@ -939,7 +944,7 @@ void vcs_VolPhase::_updateLnActCoeffJac() // Make copies of ActCoeff and Xmol_ for use in taking differences std::vector ActCoeff_Base(ActCoeff); std::vector Xmol_Base(Xmol_); - double TMoles_base = v_totalMoles; + double TMoles_base = phaseTotalMoles; /* * Loop over the columns species to be deltad @@ -950,17 +955,17 @@ void vcs_VolPhase::_updateLnActCoeffJac() * -> Note Xmol_[] and Tmoles are always positive or zero * quantities. */ - double moles_j_base = v_totalMoles * Xmol_Base[j]; - deltaMoles_j = 1.0E-7 * moles_j_base + 1.0E-13 * v_totalMoles + 1.0E-150; + double moles_j_base = phaseTotalMoles * Xmol_Base[j]; + deltaMoles_j = 1.0E-7 * moles_j_base + 1.0E-13 * phaseTotalMoles + 1.0E-150; /* * Now, update the total moles in the phase and all of the * mole fractions based on this. */ - v_totalMoles = TMoles_base + deltaMoles_j; + phaseTotalMoles = TMoles_base + deltaMoles_j; for (size_t k = 0; k < m_numSpecies; k++) { - Xmol_[k] = Xmol_Base[k] * TMoles_base / v_totalMoles; + Xmol_[k] = Xmol_Base[k] * TMoles_base / phaseTotalMoles; } - Xmol_[j] = (moles_j_base + deltaMoles_j) / v_totalMoles; + Xmol_[j] = (moles_j_base + deltaMoles_j) / phaseTotalMoles; /* * Go get new values for the activity coefficients. @@ -971,12 +976,12 @@ void vcs_VolPhase::_updateLnActCoeffJac() /* * Calculate the column of the matrix */ - double* const lnActCoeffCol = dLnActCoeffdMolNumber[j]; + double* const np_lnActCoeffCol = np_dLnActCoeffdMolNumber[j]; for (size_t k = 0; k < m_numSpecies; k++) { double tmp; tmp = (ActCoeff[k] - ActCoeff_Base[k]) / ((ActCoeff[k] + ActCoeff_Base[k]) * 0.5 * deltaMoles_j); - if (fabs(tmp - lnActCoeffCol[k]) > 1.0E-4 * fabs(tmp) + fabs(lnActCoeffCol[k])) { + if (fabs(tmp - np_lnActCoeffCol[k]) > 1.0E-4 * fabs(tmp) + fabs(np_lnActCoeffCol[k])) { // printf(" we have an error\n"); } @@ -1014,7 +1019,7 @@ void vcs_VolPhase::_updateLnActCoeffJac() * k = id of the species activity coefficient */ void -vcs_VolPhase::sendToVCS_LnActCoeffJac(double* const* const LnACJac_VCS) +vcs_VolPhase::sendToVCS_LnActCoeffJac(double* const* const np_LnACJac_VCS) { /* * update the Ln Act Coeff jacobian entries with respect to the @@ -1028,11 +1033,11 @@ vcs_VolPhase::sendToVCS_LnActCoeffJac(double* const* const LnACJac_VCS) */ for (size_t j = 0; j < m_numSpecies; j++) { size_t jglob = IndSpecies[j]; - double* const lnACJacVCS_col = LnACJac_VCS[jglob]; - const double* const lnACJac_col = dLnActCoeffdMolNumber[j]; + double* const np_lnACJacVCS_col = np_LnACJac_VCS[jglob]; + const double* const np_lnACJac_col = np_dLnActCoeffdMolNumber[j]; for (size_t k = 0; k < m_numSpecies; k++) { size_t kglob = IndSpecies[k]; - lnACJacVCS_col[kglob] = lnACJac_col[k]; + np_lnACJacVCS_col[kglob] = np_lnACJac_col[k]; } } } diff --git a/src/equil/vcs_root1d.cpp b/src/equil/vcs_root1d.cpp index 7817ccd02..b0c5b7201 100644 --- a/src/equil/vcs_root1d.cpp +++ b/src/equil/vcs_root1d.cpp @@ -487,7 +487,7 @@ QUAD_BAIL: } #ifdef DEBUG_MODE if (printLvl >= 3) { - fprintf(fp, "\nvcs_root1d failure in %d its\n", its); + fprintf(fp, "\nvcs_root1d failure in %lu its\n", its); } #endif } diff --git a/src/equil/vcs_rxnadj.cpp b/src/equil/vcs_rxnadj.cpp index d95805cd4..f40ee2ca4 100644 --- a/src/equil/vcs_rxnadj.cpp +++ b/src/equil/vcs_rxnadj.cpp @@ -16,8 +16,7 @@ #include #include -namespace VCSnonideal -{ +namespace VCSnonideal { // Calculates formation reaction step sizes. /* @@ -64,7 +63,7 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) } plogf("\n"); plogf(" --- Species KMoles Rxn_Adjustment DeltaG" - " | Comment\n"); + " | Comment\n"); } #endif /* @@ -74,13 +73,13 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) if (m_useActCoeffJac) { vcs_CalcLnActCoeffJac(VCS_DATA_PTR(m_molNumSpecies_old)); } - /************************************************************************ - ******** LOOP OVER THE FORMATION REACTIONS ***************************** - ************************************************************************/ + /************************************************************************ + ******** LOOP OVER THE FORMATION REACTIONS ***************************** + ************************************************************************/ for (size_t irxn = 0; irxn < m_numRxnRdc; ++irxn) { #ifdef DEBUG_MODE - sprintf(ANOTE,"Normal Calc"); + sprintf(ANOTE, "Normal Calc"); #endif kspec = m_indexRxnToSpecies[irxn]; @@ -94,7 +93,7 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) dnPhase_irxn = m_deltaMolNumPhase[irxn]; - if (m_molNumSpecies_old[kspec] == 0.0 && (! m_SSPhase[kspec])) { + if (m_molNumSpecies_old[kspec] == 0.0 && (!m_SSPhase[kspec])) { /********************************************************************/ /******* MULTISPECIES PHASE WITH total moles equal to zero *********/ /*******************************************************************/ @@ -111,44 +110,45 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) iph = m_phaseID[kspec]; double tphmoles = m_tPhaseMoles_old[iph]; double trphmoles = tphmoles / m_totalMolNum; - if (trphmoles > VCS_DELETE_PHASE_CUTOFF) { + Vphase = m_VolPhaseList[iph]; + if (Vphase->exists() && (trphmoles > VCS_DELETE_PHASE_CUTOFF)) { m_deltaMolNumSpecies[kspec] = m_totalMolNum * VCS_SMALL_MULTIPHASE_SPECIES; if (m_speciesStatus[kspec] == VCS_SPECIES_STOICHZERO) { m_deltaMolNumSpecies[kspec] = 0.0; #ifdef DEBUG_MODE - sprintf(ANOTE, - "MultSpec (%s): Species not born due to STOICH/PHASEPOP even though DG = %11.3E", - vcs_speciesType_string(m_speciesStatus[kspec], 15), - m_deltaGRxn_new[irxn]); + sprintf(ANOTE, "MultSpec (%s): Species not born due to STOICH/PHASEPOP even though DG = %11.3E", + vcs_speciesType_string(m_speciesStatus[kspec], 15), m_deltaGRxn_new[irxn]); #endif } else { - m_deltaMolNumSpecies[kspec] = m_totalMolNum * VCS_SMALL_MULTIPHASE_SPECIES; + m_deltaMolNumSpecies[kspec] = m_totalMolNum * VCS_SMALL_MULTIPHASE_SPECIES * 10.0; #ifdef DEBUG_MODE - sprintf(ANOTE, - "MultSpec (%s): small species born again DG = %11.3E", - vcs_speciesType_string(m_speciesStatus[kspec], 15), - m_deltaGRxn_new[irxn]); + sprintf(ANOTE, "MultSpec (%s): small species born again DG = %11.3E", + vcs_speciesType_string(m_speciesStatus[kspec], 15), m_deltaGRxn_new[irxn]); #endif } } else { #ifdef DEBUG_MODE - sprintf(ANOTE, "MultSpec (%s): phase come alive DG = %11.3E", - vcs_speciesType_string(m_speciesStatus[kspec], 15), - m_deltaGRxn_new[irxn]); + sprintf(ANOTE, "MultSpec (%s):still dead, no phase pop, even though DG = %11.3E", + vcs_speciesType_string(m_speciesStatus[kspec], 15), m_deltaGRxn_new[irxn]); #endif - Vphase = m_VolPhaseList[iph]; - size_t numSpPhase = Vphase->nSpecies(); - m_deltaMolNumSpecies[kspec] = - m_totalMolNum * 10.0 * VCS_DELETE_PHASE_CUTOFF / numSpPhase; + m_deltaMolNumSpecies[kspec] = 0.0; + if (Vphase->exists() > 0 && trphmoles > 0.0) { + m_deltaMolNumSpecies[kspec] = m_totalMolNum * VCS_SMALL_MULTIPHASE_SPECIES * 10.; +#ifdef DEBUG_MODE + sprintf(ANOTE, + "MultSpec (%s): birthed species because it was zero in a small existing phase with DG = %11.3E", + vcs_speciesType_string(m_speciesStatus[kspec], 15), m_deltaGRxn_new[irxn]); +#endif + } } } else { #ifdef DEBUG_MODE - sprintf(ANOTE, "MultSpec (%s): still dead DG = %11.3E", - vcs_speciesType_string(m_speciesStatus[kspec], 15), + sprintf(ANOTE, "MultSpec (%s): still dead DG = %11.3E", vcs_speciesType_string(m_speciesStatus[kspec], 15), m_deltaGRxn_new[irxn]); #endif m_deltaMolNumSpecies[kspec] = 0.0; + } } else { /********************************************************************/ @@ -163,12 +163,12 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) */ if (fabs(m_deltaGRxn_new[irxn]) <= m_tolmaj2) { #ifdef DEBUG_MODE - sprintf(ANOTE,"Skipped: superconverged DG = %11.3E", m_deltaGRxn_new[irxn]); + sprintf(ANOTE, "Skipped: superconverged DG = %11.3E", m_deltaGRxn_new[irxn]); if (m_debug_print_lvl >= 2) { plogf(" --- %-12.12s", m_speciesName[kspec].c_str()); plogf(" %12.4E %12.4E %12.4E | %s\n", - m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec], - m_deltaGRxn_new[irxn], ANOTE); + m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec], + m_deltaGRxn_new[irxn], ANOTE); } #endif continue; @@ -179,13 +179,12 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) */ if ((m_speciesStatus[kspec] != VCS_SPECIES_MAJOR) && (m_deltaGRxn_new[irxn] >= 0.0)) { #ifdef DEBUG_MODE - sprintf(ANOTE,"Skipped: IC = %3d and DG >0: %11.3E", - m_speciesStatus[kspec], m_deltaGRxn_new[irxn]); + sprintf(ANOTE, "Skipped: IC = %3d and DG >0: %11.3E", m_speciesStatus[kspec], m_deltaGRxn_new[irxn]); if (m_debug_print_lvl >= 2) { plogf(" --- %-12.12s", m_speciesName[kspec].c_str()); plogf(" %12.4E %12.4E %12.4E | %s\n", - m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec], - m_deltaGRxn_new[irxn], ANOTE); + m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec], + m_deltaGRxn_new[irxn], ANOTE); } #endif continue; @@ -196,7 +195,7 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) if (m_SSPhase[kspec]) { s = 0.0; } else { - s = 1.0 / m_molNumSpecies_old[kspec] ; + s = 1.0 / m_molNumSpecies_old[kspec]; } for (size_t j = 0; j < m_numComponents; ++j) { if (!m_SSPhase[j]) { @@ -207,7 +206,7 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) } for (size_t j = 0; j < m_numPhases; j++) { Vphase = m_VolPhaseList[j]; - if (! Vphase->m_singleSpecies) { + if (!Vphase->m_singleSpecies) { if (m_tPhaseMoles_old[j] > 0.0) { s -= SQUARE(dnPhase_irxn[j]) / m_tPhaseMoles_old[j]; } @@ -225,7 +224,7 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) #ifdef DEBUG_MODE if (s_old != s) { sprintf(ANOTE, "Normal calc: diag adjusted from %g " - "to %g due to act coeff", s_old, s); + "to %g due to act coeff", s_old, s); } #endif } @@ -235,20 +234,20 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) for (size_t j = 0; j < m_numComponents; ++j) { double stoicC = m_stoichCoeffRxnMatrix[irxn][j]; if (stoicC != 0.0) { - double negChangeComp = - stoicC * m_deltaMolNumSpecies[kspec]; + 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()); + "to %g due to component %lu (%10s) going neg", m_deltaMolNumSpecies[kspec], + -m_molNumSpecies_old[j] / stoicC, j, m_speciesName[j].c_str()); #endif - m_deltaMolNumSpecies[kspec] = - m_molNumSpecies_old[j] / stoicC; + m_deltaMolNumSpecies[kspec] = -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()); + "to %g due to component %lu (%10s) zero", m_deltaMolNumSpecies[kspec], + -m_molNumSpecies_old[j] / stoicC, j, m_speciesName[j].c_str()); #endif m_deltaMolNumSpecies[kspec] = 0.0; } @@ -259,8 +258,8 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) 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()); + "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]; } @@ -321,17 +320,16 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) * and the code will recover. */ #ifdef DEBUG_MODE - sprintf(ANOTE, "Delta damped from %g to %g due to delete %s", - m_deltaMolNumSpecies[kspec], - -m_molNumSpecies_old[kspec], m_speciesName[kspec].c_str()); + sprintf(ANOTE, "Delta damped from %g to %g due to delete %s", m_deltaMolNumSpecies[kspec], + -m_molNumSpecies_old[kspec], m_speciesName[kspec].c_str()); #endif m_deltaMolNumSpecies[kspec] = -m_molNumSpecies_old[kspec]; #ifdef DEBUG_MODE if (m_debug_print_lvl >= 2) { plogf(" --- %-12.12s", m_speciesName[kspec].c_str()); plogf(" %12.4E %12.4E %12.4E | %s\n", - m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec], - m_deltaGRxn_new[irxn], ANOTE); + m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec], + m_deltaGRxn_new[irxn], ANOTE); } #endif continue; @@ -352,18 +350,18 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) #ifdef DEBUG_MODE if (k != kspec) { - sprintf(ANOTE, "Delete component SS phase %d named %s - SS phases only", - iphDel, m_speciesName[k].c_str()); + sprintf(ANOTE, "Delete component SS phase %lu named %s - SS phases only", iphDel, + m_speciesName[k].c_str()); } else { - sprintf(ANOTE, "Delete this SS phase %d - SS components only", iphDel); + sprintf(ANOTE, "Delete this SS phase %lu - SS components only", iphDel); } if (m_debug_print_lvl >= 2) { plogf(" --- %-12.12s", m_speciesName[kspec].c_str()); plogf(" %12.4E %12.4E %12.4E | %s\n", - m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec], - m_deltaGRxn_new[irxn], ANOTE); + m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec], + m_deltaGRxn_new[irxn], ANOTE); plogf(" --- vcs_RxnStepSizes Special section to set up to delete %s", - m_speciesName[k].c_str()); + m_speciesName[k].c_str()); plogendl(); } #endif @@ -390,8 +388,8 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) if (m_debug_print_lvl >= 2) { plogf(" --- %-12.12s", m_speciesName[kspec].c_str()); plogf(" %12.4E %12.4E %12.4E | %s\n", - m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec], - m_deltaGRxn_new[irxn], ANOTE); + m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec], + m_deltaGRxn_new[irxn], ANOTE); } #endif } /* End of loop over m_speciesUnknownType */ @@ -457,13 +455,13 @@ int VCS_SOLVE::vcs_rxn_adj_cg() */ for (irxn = 0; irxn < m_numRxnRdc; ++irxn) { #ifdef DEBUG_MODE - sprintf(ANOTE,"Normal Calc"); + sprintf(ANOTE, "Normal Calc"); #endif kspec = m_indexRxnToSpecies[irxn]; dnPhase_irxn = m_deltaMolNumPhase[irxn]; - if (m_molNumSpecies_old[kspec] == 0.0 && (! m_SSPhase[kspec])) { + if (m_molNumSpecies_old[kspec] == 0.0 && (!m_SSPhase[kspec])) { /* *******************************************************************/ /* **** MULTISPECIES PHASE WITH total moles equal to zero ************/ /* *******************************************************************/ @@ -497,11 +495,11 @@ int VCS_SOLVE::vcs_rxn_adj_cg() */ if (fabs(m_deltaGRxn_new[irxn]) <= m_tolmaj2) { #ifdef DEBUG_MODE - sprintf(ANOTE,"Skipped: converged DG = %11.3E\n", m_deltaGRxn_new[irxn]); + sprintf(ANOTE, "Skipped: converged DG = %11.3E\n", m_deltaGRxn_new[irxn]); plogf(" --- "); plogf("%-12.12s", m_speciesName[kspec].c_str()); - plogf(" %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec], - m_deltaMolNumSpecies[kspec], ANOTE); + plogf(" %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec], + m_deltaMolNumSpecies[kspec], ANOTE); #endif continue; } @@ -511,12 +509,11 @@ int VCS_SOLVE::vcs_rxn_adj_cg() */ if (m_speciesStatus[kspec] <= VCS_SPECIES_MINOR && m_deltaGRxn_new[irxn] >= 0.0) { #ifdef DEBUG_MODE - sprintf(ANOTE,"Skipped: IC = %3d and DG >0: %11.3E\n", - m_speciesStatus[kspec], m_deltaGRxn_new[irxn]); + sprintf(ANOTE, "Skipped: IC = %3d and DG >0: %11.3E\n", m_speciesStatus[kspec], m_deltaGRxn_new[irxn]); plogf(" --- "); plogf("%-12.12s", m_speciesName[kspec].c_str()); plogf(" %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec], - m_deltaMolNumSpecies[kspec], ANOTE); + m_deltaMolNumSpecies[kspec], ANOTE); #endif continue; } @@ -529,7 +526,7 @@ int VCS_SOLVE::vcs_rxn_adj_cg() s = 1.0 / m_molNumSpecies_old[kspec]; } for (j = 0; j < m_numComponents; ++j) { - if (! m_SSPhase[j]) { + if (!m_SSPhase[j]) { s += SQUARE(m_stoichCoeffRxnMatrix[irxn][j]) / m_molNumSpecies_old[j]; } } @@ -589,10 +586,10 @@ int VCS_SOLVE::vcs_rxn_adj_cg() */ if (dss != 0.0) { m_molNumSpecies_old[kspec] += dss; - m_tPhaseMoles_old[m_phaseID[kspec]] += dss; + m_tPhaseMoles_old[m_phaseID[kspec]] += dss; for (j = 0; j < m_numComponents; ++j) { m_molNumSpecies_old[j] += dss * m_stoichCoeffRxnMatrix[irxn][j]; - m_tPhaseMoles_old[m_phaseID[j]] += dss * m_stoichCoeffRxnMatrix[irxn][j]; + m_tPhaseMoles_old[m_phaseID[j]] += dss * m_stoichCoeffRxnMatrix[irxn][j]; } m_molNumSpecies_old[k] = 0.0; m_tPhaseMoles_old[m_phaseID[k]] = 0.0; @@ -619,7 +616,7 @@ int VCS_SOLVE::vcs_rxn_adj_cg() plogf(" --- "); plogf("%-12.12s", m_speciesName[kspec].c_str()); plogf(" %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec], - m_deltaMolNumSpecies[kspec], ANOTE); + m_deltaMolNumSpecies[kspec], ANOTE); #endif } /* End of loop over non-component stoichiometric formation reactions */ @@ -664,10 +661,10 @@ double VCS_SOLVE::vcs_Hessian_diag_adj(size_t irxn, double hessianDiag_Ideal) } if (hessActCoef >= 0.0) { diag += hessActCoef; - } else if (fabs(hessActCoef) < 0.6666 * hessianDiag_Ideal) { + } else if (fabs(hessActCoef) < 0.6666 * hessianDiag_Ideal) { diag += hessActCoef; } else { - diag -= 0.6666 * hessianDiag_Ideal; + diag -= 0.6666 * hessianDiag_Ideal; } return diag; } @@ -687,11 +684,15 @@ double VCS_SOLVE::vcs_Hessian_actCoeff_diag(size_t irxn) double* sc_irxn; kspec = m_indexRxnToSpecies[irxn]; kph = m_phaseID[kspec]; + double np_kspec = m_tPhaseMoles_old[kph]; + if (np_kspec < 1.0E-13) { + np_kspec = 1.0E-13; + } sc_irxn = m_stoichCoeffRxnMatrix[irxn]; /* * First the diagonal term of the Jacobian */ - s = m_dLnActCoeffdMolNum[kspec][kspec]; + s = m_np_dLnActCoeffdMolNum[kspec][kspec] / np_kspec; /* * Next, the other terms. Note this only a loop over the components * So, it's not too expensive to calculate. @@ -700,12 +701,16 @@ double VCS_SOLVE::vcs_Hessian_actCoeff_diag(size_t irxn) if (!m_SSPhase[l]) { for (k = 0; k < m_numComponents; ++k) { if (m_phaseID[k] == m_phaseID[l]) { - s += sc_irxn[k] * sc_irxn[l] * m_dLnActCoeffdMolNum[k][l]; + double np = m_tPhaseMoles_old[m_phaseID[k]]; + if (np > 0.0) { + s += sc_irxn[k] * sc_irxn[l] * m_np_dLnActCoeffdMolNum[k][l] / np; + } } } if (kph == m_phaseID[l]) { - s += sc_irxn[l] * (m_dLnActCoeffdMolNum[kspec][l] + m_dLnActCoeffdMolNum[l][kspec]); + s += sc_irxn[l] * (m_np_dLnActCoeffdMolNum[kspec][l] + m_np_dLnActCoeffdMolNum[l][kspec]) / np_kspec; } + } } return s; @@ -739,7 +744,7 @@ void VCS_SOLVE::vcs_CalcLnActCoeffJac(const double* const moleSpeciesVCS) * -> This scatter calculation is carried out in the * vcs_VolPhase object. */ - Vphase->sendToVCS_LnActCoeffJac(m_dLnActCoeffdMolNum.baseDataAddr()); + Vphase->sendToVCS_LnActCoeffJac(m_np_dLnActCoeffdMolNum.baseDataAddr()); } } } @@ -778,9 +783,8 @@ void VCS_SOLVE::vcs_CalcLnActCoeffJac(const double* const moleSpeciesVCS) * Note, this is a dangerous routine that leaves the underlying objects in * an unknown state. */ -double VCS_SOLVE::deltaG_Recalc_Rxn(const int stateCalc, - const size_t irxn, const double* const molNum, - double* const ac, double* const mu_i) +double VCS_SOLVE::deltaG_Recalc_Rxn(const int stateCalc, const size_t irxn, const double* const molNum, double* const ac, + double* const mu_i) { size_t kspec = irxn + m_numComponents; int* pp_ptr = m_phaseParticipation[irxn]; @@ -812,8 +816,7 @@ double VCS_SOLVE::deltaG_Recalc_Rxn(const int stateCalc, * * @return Returns the optimized step length found by the search */ -double VCS_SOLVE::vcs_line_search(const size_t irxn, const double dx_orig, - char* const ANOTE) +double VCS_SOLVE::vcs_line_search(const size_t irxn, const double dx_orig, char* const ANOTE) #else double VCS_SOLVE::vcs_line_search(const size_t irxn, const double dx_orig) #endif @@ -833,9 +836,7 @@ double VCS_SOLVE::vcs_line_search(const size_t irxn, const double dx_orig) * Calculate the deltaG value at the dx = 0.0 point */ vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD); - double deltaGOrig = deltaG_Recalc_Rxn(VCS_STATECALC_OLD, - irxn, molNumBase, acBase, - VCS_DATA_PTR(m_feSpecies_old)); + double deltaGOrig = deltaG_Recalc_Rxn(VCS_STATECALC_OLD, irxn, molNumBase, acBase, VCS_DATA_PTR(m_feSpecies_old)); double forig = fabs(deltaGOrig) + 1.0E-15; if (deltaGOrig > 0.0) { if (dx_orig > 0.0) { @@ -844,7 +845,7 @@ double VCS_SOLVE::vcs_line_search(const size_t irxn, const double dx_orig) if (m_debug_print_lvl >= 2) { //plogf(" --- %s :Warning possible error dx>0 dg > 0\n", SpName[kspec]); } - sprintf(ANOTE,"Rxn reduced to zero step size in line search: dx>0 dg > 0"); + sprintf(ANOTE, "Rxn reduced to zero step size in line search: dx>0 dg > 0"); #endif return dx; } @@ -855,7 +856,7 @@ double VCS_SOLVE::vcs_line_search(const size_t irxn, const double dx_orig) if (m_debug_print_lvl >= 2) { //plogf(" --- %s :Warning possible error dx<0 dg < 0\n", SpName[kspec]); } - sprintf(ANOTE,"Rxn reduced to zero step size in line search: dx<0 dg < 0"); + sprintf(ANOTE, "Rxn reduced to zero step size in line search: dx<0 dg < 0"); #endif return dx; } @@ -867,17 +868,16 @@ double VCS_SOLVE::vcs_line_search(const size_t irxn, const double dx_orig) } vcs_dcopy(VCS_DATA_PTR(m_molNumSpecies_new), molNumBase, m_numSpeciesRdc); - molSum = molNumBase[kspec]; + molSum = molNumBase[kspec]; m_molNumSpecies_new[kspec] = molNumBase[kspec] + dx_orig; for (k = 0; k < m_numComponents; k++) { m_molNumSpecies_new[k] = molNumBase[k] + sc_irxn[k] * dx_orig; - molSum += molNumBase[k]; + molSum += molNumBase[k]; } vcs_setFlagsVolPhases(false, VCS_STATECALC_NEW); - double deltaG1 = deltaG_Recalc_Rxn(VCS_STATECALC_NEW, - irxn, VCS_DATA_PTR(m_molNumSpecies_new), - ac, VCS_DATA_PTR(m_feSpecies_new)); + double deltaG1 = deltaG_Recalc_Rxn(VCS_STATECALC_NEW, irxn, VCS_DATA_PTR(m_molNumSpecies_new), + ac, VCS_DATA_PTR(m_feSpecies_new)); /* * If deltaG hasn't switched signs when going the full distance @@ -892,7 +892,7 @@ double VCS_SOLVE::vcs_line_search(const size_t irxn, const double dx_orig) * If we have decreased somewhat, the deltaG return after finding * a better estimate for the line search. */ - if (fabs(deltaG1) < 0.8*forig) { + if (fabs(deltaG1) < 0.8 * forig) { if (deltaG1 * deltaGOrig < 0.0) { slope = (deltaG1 - deltaGOrig) / dx_orig; dx = -deltaGOrig / slope; @@ -915,9 +915,8 @@ double VCS_SOLVE::vcs_line_search(const size_t irxn, const double dx_orig) m_molNumSpecies_new[k] = molNumBase[k] + sc_irxn[k] * dx; } vcs_setFlagsVolPhases(false, VCS_STATECALC_NEW); - double deltaG = deltaG_Recalc_Rxn(VCS_STATECALC_NEW, - irxn, VCS_DATA_PTR(m_molNumSpecies_new), - ac, VCS_DATA_PTR(m_feSpecies_new)); + double deltaG = deltaG_Recalc_Rxn(VCS_STATECALC_NEW, irxn, VCS_DATA_PTR(m_molNumSpecies_new), + ac, VCS_DATA_PTR(m_feSpecies_new)); /* * If deltaG hasn't switched signs when going the full distance * then we are heading in the appropriate direction, and @@ -939,19 +938,16 @@ double VCS_SOLVE::vcs_line_search(const size_t irxn, const double dx_orig) } } -finalize: - vcs_setFlagsVolPhases(false, VCS_STATECALC_NEW); + finalize: vcs_setFlagsVolPhases(false, VCS_STATECALC_NEW); if (its >= MAXITS) { #ifdef DEBUG_MODE - sprintf(ANOTE,"Rxn reduced to zero step size from %g to %g (MAXITS)", - dx_orig, dx); + sprintf(ANOTE, "Rxn reduced to zero step size from %g to %g (MAXITS)", dx_orig, dx); return dx; #endif } #ifdef DEBUG_MODE if (dx != dx_orig) { - sprintf(ANOTE,"Line Search reduced step size from %g to %g", - dx_orig, dx); + sprintf(ANOTE, "Line Search reduced step size from %g to %g", dx_orig, dx); } #endif diff --git a/src/equil/vcs_solve.cpp b/src/equil/vcs_solve.cpp index b0674f111..3db3f2505 100644 --- a/src/equil/vcs_solve.cpp +++ b/src/equil/vcs_solve.cpp @@ -218,7 +218,7 @@ void VCS_SOLVE::vcs_initSizes(const size_t nspecies0, const size_t nelements, */ m_useActCoeffJac = true; if (m_useActCoeffJac) { - m_dLnActCoeffdMolNum.resize(nspecies0, nspecies0, 0.0); + m_np_dLnActCoeffdMolNum.resize(nspecies0, nspecies0, 0.0); } m_PMVolumeSpecies.resize(nspecies0, 0.0); diff --git a/src/equil/vcs_solve_TP.cpp b/src/equil/vcs_solve_TP.cpp index 0cd9c256d..372399ad1 100644 --- a/src/equil/vcs_solve_TP.cpp +++ b/src/equil/vcs_solve_TP.cpp @@ -2075,7 +2075,7 @@ double VCS_SOLVE::vcs_minor_alt_calc(size_t kspec, size_t irxn, bool* do_delete /* * get the diagonal of the activity coefficient jacobian */ - s = m_dLnActCoeffdMolNum[kspec][kspec]; + s = m_np_dLnActCoeffdMolNum[kspec][kspec] / (m_tPhaseMoles_old[iph]); // s *= (m_tPhaseMoles_old[iph]); /* * We fit it to a power law approximation of the activity coefficient @@ -5463,7 +5463,7 @@ void VCS_SOLVE::vcs_switch_pos(const bool ifunc, const size_t k1, const size_t k std::swap(m_formulaMatrix[j][k1], m_formulaMatrix[j][k2]); } if (m_useActCoeffJac) { - vcs_switch2D(m_dLnActCoeffdMolNum.baseDataAddr(), k1, k2); + vcs_switch2D(m_np_dLnActCoeffdMolNum.baseDataAddr(), k1, k2); } std::swap(m_speciesStatus[k1], m_speciesStatus[k2]); /* diff --git a/src/equil/vcs_util.cpp b/src/equil/vcs_util.cpp index 311aae7fe..e13276dd5 100644 --- a/src/equil/vcs_util.cpp +++ b/src/equil/vcs_util.cpp @@ -17,8 +17,7 @@ using namespace std; -namespace VCSnonideal -{ +namespace VCSnonideal { /***************************************************************************/ /***************************************************************************/ @@ -33,10 +32,10 @@ void vcs_dzero(double* vector, int length) * Zeroes a double vector *************************************************************************/ { - int i; - for (i = 0; i < length; i++) { - vector[i] = 0.0; - } + int i; + for (i = 0; i < length; i++) { + vector[i] = 0.0; + } } /* vcs_dzero() ***********************************************************/ #endif /***************************************************************************/ @@ -52,10 +51,10 @@ void vcs_izero(int* vector, int length) * Zeroes an int vector *************************************************************************/ { - int i; - for (i = 0; i < length; i++) { - vector[i] = 0; - } + int i; + for (i = 0; i < length; i++) { + vector[i] = 0; + } } /* vcs_izero() ***********************************************************/ #endif /***************************************************************************/ @@ -72,10 +71,10 @@ void vcs_dcopy(double* const vec_to, const double* const vec_from, int length) * Copies a double vector ***************************************************************************/ { - int i; - for (i = 0; i < length; i++) { - vec_to[i] = vec_from[i]; - } + int i; + for (i = 0; i < length; i++) { + vec_to[i] = vec_from[i]; + } } /* vcs_dzero() *************************************************************/ #endif /*****************************************************************************/ @@ -92,10 +91,10 @@ void vcs_icopy(int* vec_to, int* vec_from, int length) * copies an int vector ***************************************************************************/ { - int i; - for (i = 0; i < length; i++) { - vec_to[i] = vec_from[i]; - } + int i; + for (i = 0; i < length; i++) { + vec_to[i] = vec_from[i]; + } } /* vcs_dzero() *************************************************************/ #endif @@ -111,26 +110,25 @@ void vcs_icopy(int* vec_to, int* vec_from, int length) */ void vcs_vdzero(std::vector &vvv, int len) { - if (len < 0) { - std::fill(vvv.begin(), vvv.end(), 0.0); - } else { - std::fill_n(vvv.begin(), len, 0.0); - } + if (len < 0) { + std::fill(vvv.begin(), vvv.end(), 0.0); + } else { + std::fill_n(vvv.begin(), len, 0.0); + } } #endif -double vcs_l2norm(const std::vector vec) -{ - size_t len = vec.size(); - if (len == 0) { - return 0.0; - } - double sum = 0.0; - std::vector::const_iterator pos; - for (pos = vec.begin(); pos != vec.end(); ++pos) { - sum += (*pos) * (*pos); - } - return std::sqrt(sum/len); +double vcs_l2norm(const std::vector vec) { + size_t len = vec.size(); + if (len == 0) { + return 0.0; + } + double sum = 0.0; + std::vector::const_iterator pos; + for (pos = vec.begin(); pos != vec.end(); ++pos) { + sum += (*pos) * (*pos); + } + return std::sqrt(sum / len); } /*****************************************************************************/ @@ -145,16 +143,14 @@ double vcs_l2norm(const std::vector vec) */ void vcs_vizero(std::vector &vvv, int len) { - if (len < 0) { - std::fill(vvv.begin(), vvv.end(), 0.0); - } else { - std::fill_n(vvv.begin(), len, 0.0); - } + if (len < 0) { + std::fill(vvv.begin(), vvv.end(), 0.0); + } else { + std::fill_n(vvv.begin(), len, 0.0); + } } #endif - - #ifndef USE_MEMSET /* * vcs_vdcopy @@ -166,13 +162,12 @@ void vcs_vizero(std::vector &vvv, int len) * @param length Length of the copy */ void vcs_vdcopy(std::vector &vec_to, - const std::vector & vec_from, int length) + const std::vector & vec_from, int length) { - std::copy(vec_from.begin(), vec_from.begin() + length, vec_to.begin()); + std::copy(vec_from.begin(), vec_from.begin() + length, vec_to.begin()); } #endif - #ifndef USE_MEMSET /* * vcs_vicopy @@ -184,9 +179,9 @@ void vcs_vdcopy(std::vector &vec_to, * @param length Length of the copy */ void vcs_vicopy(std::vector &vec_to, - const std::vector & vec_from, int length) + const std::vector & vec_from, int length) { - std::copy(vec_from.begin(), vec_from.begin() + length, vec_to.begin()); + std::copy(vec_from.begin(), vec_from.begin() + length, vec_to.begin()); } #endif @@ -202,30 +197,29 @@ void vcs_vicopy(std::vector &vec_to, * RETURN * return index of the greatest value on X(*) searched */ -size_t vcs_optMax(const double* x, const double* xSize, size_t j, size_t n) -{ - size_t i; - size_t largest = j; - double big = x[j]; - if (xSize) { - assert(xSize[j] > 0.0); - big *= xSize[j]; - for (i = j + 1; i < n; ++i) { - assert(xSize[i] > 0.0); - if ((x[i]*xSize[i]) > big) { - largest = i; - big = x[i]*xSize[i]; - } - } - } else { - for (i = j + 1; i < n; ++i) { - if (x[i] > big) { - largest = i; - big = x[i]; - } - } - } - return largest; +size_t vcs_optMax(const double* x, const double* xSize, size_t j, size_t n) { + size_t i; + size_t largest = j; + double big = x[j]; + if (xSize) { + assert(xSize[j] > 0.0); + big *= xSize[j]; + for (i = j + 1; i < n; ++i) { + assert(xSize[i] > 0.0); + if ((x[i] * xSize[i]) > big) { + largest = i; + big = x[i] * xSize[i]; + } + } + } else { + for (i = j + 1; i < n; ++i) { + if (x[i] > big) { + largest = i; + big = x[i]; + } + } + } + return largest; } int vcs_max_int(const int* vector, int length) @@ -237,45 +231,44 @@ int vcs_max_int(const int* vector, int length) * returns the maximum integer in a list. ***************************************************************************/ { - int i, retn; - if (vector == NULL || length <= 0) { - return 0; - } - retn = vector[0]; - for (i = 1; i < length; i++) { - retn = std::max(retn, vector[i]); - } - return retn; + int i, retn; + if (vector == NULL || length <= 0) { + return 0; + } + retn = vector[0]; + for (i = 1; i < length; i++) { + retn = std::max(retn, vector[i]); + } + return retn; } - //==================================================================================================================== #ifdef DEBUG_HKM static void mlequ_matrixDump(double* c, int idem, int n) { - int i, j; - printf("vcsUtil_mlequ() MATRIX DUMP --------------------------------------------------\n"); - printf(" "); - for (j = 0; j < n; ++j) { - printf(" % 3d ", j); - } - printf("\n"); - for (j = 0; j < n; ++j) { - printf("-----------"); - } - printf("\n"); - for (i = 0; i < n; ++i) { - printf(" %3d | ", i); - for (j = 0; j < n; ++j) { - printf("% 10.3e ", c[i + j * idem]); - } - printf("\n"); - } - for (j = 0; j < n; ++j) { - printf("-----------"); - } - printf("\n"); - printf("vcsUtil_mlequ() END MATRIX DUMP --------------------------------------------------\n"); + int i, j; + printf("vcsUtil_mlequ() MATRIX DUMP --------------------------------------------------\n"); + printf(" "); + for (j = 0; j < n; ++j) { + printf(" % 3d ", j); + } + printf("\n"); + for (j = 0; j < n; ++j) { + printf("-----------"); + } + printf("\n"); + for (i = 0; i < n; ++i) { + printf(" %3d | ", i); + for (j = 0; j < n; ++j) { + printf("% 10.3e ", c[i + j * idem]); + } + printf("\n"); + } + for (j = 0; j < n; ++j) { + printf("-----------"); + } + printf("\n"); + printf("vcsUtil_mlequ() END MATRIX DUMP --------------------------------------------------\n"); } #endif @@ -290,17 +283,17 @@ static void mlequ_matrixDump(double* c, int idem, int n) * @param irowa first row to swap * @param irowb second row to swap */ -static void vcsUtil_swapRows(double* c, size_t idem, size_t n, double* b, size_t m, size_t irowa, size_t irowb) -{ - if (irowa == irowb) { - return; - } - for (size_t j = 0; j < n; j++) { - std::swap(c[irowa + j * idem], c[irowb + j * idem]); - } - for (size_t j = 0; j < m; j++) { - std::swap(b[irowa + j * idem], b[irowb + j * idem]); - } +static void vcsUtil_swapRows(double* c, size_t idem, size_t n, double* b, + size_t m, size_t irowa, size_t irowb) { + if (irowa == irowb) { + return; + } + for (size_t j = 0; j < n; j++) { + std::swap(c[irowa + j * idem], c[irowb + j * idem]); + } + for (size_t j = 0; j < m; j++) { + std::swap(b[irowa + j * idem], b[irowb + j * idem]); + } } //==================================================================================================================== //! Swap rows in the c matrix and the b rhs matrix to lower the condition number of the matrix @@ -311,87 +304,89 @@ static void vcsUtil_swapRows(double* c, size_t idem, size_t n, double* b, size_t * @param b RHS of the Ax=b problem to solve * @param m Number of rhs to solve */ -static void vcsUtil_mlequ_preprocess(double* c, size_t idem, size_t n, double* b, size_t m) -{ - size_t j = 0; - std::vector irowUsed(n, 0); +static void vcsUtil_mlequ_preprocess(double* c, size_t idem, size_t n, + double* b, size_t m) { + size_t j = 0; + std::vector irowUsed(n, 0); - for (j = 0; j < n; j++) { - int numNonzero = 0; - size_t inonzero = npos; - for (size_t i = 0; i < n; i++) { - if (c[i + j * idem] != 0.0) { - numNonzero++; - inonzero = i; - } - } - if (numNonzero == 1) { - if (inonzero != j) { - if (irowUsed[inonzero] == 0) { - vcsUtil_swapRows(c, idem, n, b, m, j, inonzero); + for (j = 0; j < n; j++) { + int numNonzero = 0; + size_t inonzero = npos; + for (size_t i = 0; i < n; i++) { + if (c[i + j * idem] != 0.0) { + numNonzero++; + inonzero = i; + } + } + if (numNonzero == 1) { + if (inonzero != j) { + if (irowUsed[inonzero] == 0) { + vcsUtil_swapRows(c, idem, n, b, m, j, inonzero); #ifdef DEBUG_HKM - // mlequ_matrixDump(c, idem, n); + // mlequ_matrixDump(c, idem, n); #endif - } - } - irowUsed[j] = 1; - } - } + } + } + irowUsed[j] = 1; + } + } - for (j = 0; j < n; j++) { - if (c[j + j * idem] == 0.0) { - int numNonzero = 0; - size_t inonzero = npos; - for (size_t i = 0; i < n; i++) { - if (! irowUsed[i]) { - if (c[i + j * idem] != 0.0) { - if ((c[i + i * idem] == 0.0) || (c[j + i * idem] != 0.0)) { - numNonzero++; - inonzero = i; - } - } - } - } - if (numNonzero == 1) { - if (inonzero != j) { - if (irowUsed[inonzero] == 0) { - vcsUtil_swapRows(c, idem, n, b, m, j, inonzero); + for (j = 0; j < n; j++) { + if (c[j + j * idem] == 0.0) { + int numNonzero = 0; + size_t inonzero = npos; + for (size_t i = 0; i < n; i++) { + if (!irowUsed[i]) { + if (c[i + j * idem] != 0.0) { + if ((c[i + i * idem] == 0.0) + || (c[j + i * idem] != 0.0)) { + numNonzero++; + inonzero = i; + } + } + } + } + if (numNonzero == 1) { + if (inonzero != j) { + if (irowUsed[inonzero] == 0) { + vcsUtil_swapRows(c, idem, n, b, m, j, inonzero); #ifdef DEBUG_HKM - // mlequ_matrixDump(c, idem, n); + // mlequ_matrixDump(c, idem, n); #endif - } - } - irowUsed[j] = 1; - } - } - } + } + } + irowUsed[j] = 1; + } + } + } - for (j = 0; j < n; j++) { - if (c[j + j * idem] == 0.0) { - int numNonzero = 0; - size_t inonzero = npos; - for (size_t i = 0; i < n; i++) { - if (! irowUsed[i]) { - if (c[i + j * idem] != 0.0) { - if ((c[i + i * idem] == 0.0) || (c[j + i * idem] != 0.0)) { - numNonzero++; - inonzero = i; - } - } - } - } - if (inonzero != npos) { - if (inonzero != j) { - if (irowUsed[inonzero] == 0) { - vcsUtil_swapRows(c, idem, n, b, m, j, inonzero); + for (j = 0; j < n; j++) { + if (c[j + j * idem] == 0.0) { + int numNonzero = 0; + size_t inonzero = npos; + for (size_t i = 0; i < n; i++) { + if (!irowUsed[i]) { + if (c[i + j * idem] != 0.0) { + if ((c[i + i * idem] == 0.0) + || (c[j + i * idem] != 0.0)) { + numNonzero++; + inonzero = i; + } + } + } + } + if (inonzero != npos) { + if (inonzero != j) { + if (irowUsed[inonzero] == 0) { + vcsUtil_swapRows(c, idem, n, b, m, j, inonzero); #ifdef DEBUG_HKM - // mlequ_matrixDump(c, idem, n); + // mlequ_matrixDump(c, idem, n); #endif - } - } - } - } - } + } + } + } + } + } } //==================================================================================================================== // Invert an n x n matrix and solve m rhs's @@ -424,115 +419,113 @@ static void vcsUtil_mlequ_preprocess(double* c, size_t idem, size_t n, double* b * (each column is a new rhs) * @param m number of rhs's */ -int vcsUtil_mlequ(double* c, size_t idem, size_t n, double* b, size_t m) -{ - size_t k; +int vcsUtil_mlequ(double* c, size_t idem, size_t n, double* b, size_t m) { + size_t k; #ifdef DEBUG_HKM - // mlequ_matrixDump(c, idem, n); + // mlequ_matrixDump(c, idem, n); #endif - vcsUtil_mlequ_preprocess(c, idem, n, b, m); + vcsUtil_mlequ_preprocess(c, idem, n, b, m); #ifdef DEBUG_HKM - // mlequ_matrixDump(c, idem, n); - static int s_numCalls = 0; - s_numCalls++; + // mlequ_matrixDump(c, idem, n); + static int s_numCalls = 0; + s_numCalls++; #endif - double R; - if (n > idem || n <= 0) { - plogf("vcsUtil_mlequ ERROR: badly dimensioned matrix: %d %d\n", n, idem); - return 1; - } + double R; + if (n > idem || n <= 0) { + plogf("vcsUtil_mlequ ERROR: badly dimensioned matrix: %d %d\n", n, idem); + return 1; + } #ifdef DEBUG_HKM - int dmatrix = 0; - for (size_t i = 0; i < n; ++i) { - bool notFound = true; - for (size_t j = 0; j < n; ++j) { - if (c[i + j * idem] != 0.0) { - notFound = false; - } - } - if (notFound) { - printf(" vcsUtil_mlequ ERROR(): row %d is identically zero\n", i); - } - } - for (size_t j = 0; j < n; ++j) { - bool notFound = true; - for (size_t i = 0; i < n; ++i) { - if (c[i + j * idem] != 0.0) { - notFound = false; - } - } - if (notFound) { - printf(" vcsUtil_mlequ ERROR(): column %d is identically zero\n", j); - } - } - // if (s_numCalls >= 32) { - // printf("vcsUtil_mlequ: we are here\n"); - // dmatrix = 1; - // } + int dmatrix = 0; + for (size_t i = 0; i < n; ++i) { + bool notFound = true; + for (size_t j = 0; j < n; ++j) { + if (c[i + j * idem] != 0.0) { + notFound = false; + } + } + if (notFound) { + printf(" vcsUtil_mlequ ERROR(): row %d is identically zero\n", i); + } + } + for (size_t j = 0; j < n; ++j) { + bool notFound = true; + for (size_t i = 0; i < n; ++i) { + if (c[i + j * idem] != 0.0) { + notFound = false; + } + } + if (notFound) { + printf(" vcsUtil_mlequ ERROR(): column %d is identically zero\n", j); + } + } + // if (s_numCalls >= 32) { + // printf("vcsUtil_mlequ: we are here\n"); + // dmatrix = 1; + // } - if (dmatrix) { - mlequ_matrixDump(c, idem, n); - } + if (dmatrix) { + mlequ_matrixDump(c, idem, n); + } #endif - /* - * Loop over the rows - * -> At the end of each loop, the only nonzero entry in the column - * will be on the diagonal. We can therfore just invert the - * diagonal at the end of the program to solve the equation system. - */ - for (size_t i = 0; i < n; ++i) { - if (c[i + i * idem] == 0.0) { - /* - * Do a simple form of row pivoting to find a non-zero pivot - */ - for (k = i + 1; k < n; ++k) { - if (c[k + i * idem] != 0.0) { - goto FOUND_PIVOT; - } - } - plogf("vcsUtil_mlequ ERROR: Encountered a zero column: %d\n", i); + /* + * Loop over the rows + * -> At the end of each loop, the only nonzero entry in the column + * will be on the diagonal. We can therfore just invert the + * diagonal at the end of the program to solve the equation system. + */ + for (size_t i = 0; i < n; ++i) { + if (c[i + i * idem] == 0.0) { + /* + * Do a simple form of row pivoting to find a non-zero pivot + */ + for (k = i + 1; k < n; ++k) { + if (c[k + i * idem] != 0.0) { + goto FOUND_PIVOT; + } + } + plogf("vcsUtil_mlequ ERROR: Encountered a zero column: %d\n", i); #ifdef DEBUG_HKM - plogf(" call # %d\n", s_numCalls); + plogf(" call # %d\n", s_numCalls); #endif #ifdef DEBUG_HKM - mlequ_matrixDump(c, idem, n); + mlequ_matrixDump(c, idem, n); #endif - return 1; -FOUND_PIVOT: - ; - for (size_t j = 0; j < n; ++j) { - c[i + j * idem] += c[k + j * idem]; - } - for (size_t j = 0; j < m; ++j) { - b[i + j * idem] += b[k + j * idem]; - } - } + return 1; + FOUND_PIVOT: ; + for (size_t j = 0; j < n; ++j) { + c[i + j * idem] += c[k + j * idem]; + } + for (size_t j = 0; j < m; ++j) { + b[i + j * idem] += b[k + j * idem]; + } + } - for (size_t l = 0; l < n; ++l) { - if (l != i && c[l + i * idem] != 0.0) { - R = c[l + i * idem] / c[i + i * idem]; - c[l + i * idem] = 0.0; - for (size_t j = i+1; j < n; ++j) { - c[l + j * idem] -= c[i + j * idem] * R; - } - for (size_t j = 0; j < m; ++j) { - b[l + j * idem] -= b[i + j * idem] * R; - } - } - } - } - /* - * The negative in the last expression is due to the form of B upon - * input - */ - for (size_t i = 0; i < n; ++i) { - for (size_t j = 0; j < m; ++j) { - b[i + j * idem] = -b[i + j * idem] / c[i + i*idem]; - } - } - return 0; + for (size_t l = 0; l < n; ++l) { + if (l != i && c[l + i * idem] != 0.0) { + R = c[l + i * idem] / c[i + i * idem]; + c[l + i * idem] = 0.0; + for (size_t j = i + 1; j < n; ++j) { + c[l + j * idem] -= c[i + j * idem] * R; + } + for (size_t j = 0; j < m; ++j) { + b[l + j * idem] -= b[i + j * idem] * R; + } + } + } + } + /* + * The negative in the last expression is due to the form of B upon + * input + */ + for (size_t i = 0; i < n; ++i) { + for (size_t j = 0; j < m; ++j) { + b[i + j * idem] = -b[i + j * idem] / c[i + i * idem]; + } + } + return 0; } //==================================================================================================================== // Linear equation solution by Gauss-Jordan elimination for multiple rhs vectors @@ -564,104 +557,102 @@ FOUND_PIVOT: * (each column is a new rhs) * @param m number of rhs's */ -int vcsUtil_gaussj(double* c, size_t idem, size_t n, double* b, size_t m) -{ +int vcsUtil_gaussj(double* c, size_t idem, size_t n, double* b, size_t m) { - size_t i, j, k, l, ll; - size_t irow = npos; - size_t icol = npos; - bool needInverse = false; - double pivinv; + size_t i, j, k, l, ll; + size_t irow = npos; + size_t icol = npos; + bool needInverse = false; + double pivinv; #ifdef DEBUG_HKM - static int s_numCalls = 0; - s_numCalls++; + static int s_numCalls = 0; + s_numCalls++; #endif #ifdef DEBUG_HKM - // mlequ_matrixDump(c, idem, n); + // mlequ_matrixDump(c, idem, n); #endif - /* - * Preprocess the problem - */ - vcsUtil_mlequ_preprocess(c, idem, n, b, m); + /* + * Preprocess the problem + */ + vcsUtil_mlequ_preprocess(c, idem, n, b, m); #ifdef DEBUG_HKM - // mlequ_matrixDump(c, idem, n); + // mlequ_matrixDump(c, idem, n); #endif - std::vector indxc(n); - std::vector indxr(n); - std::vector ipiv(n, 0); - doublereal big = 0.0; - /* - * This is the main loop over the columns to be reduced. - */ - for (i = 0; i < n; i++) { - big = 0.0; - for (j = 0; j < n; j++) { - if (ipiv[j] != 1) { - for (k = 0; k < n; k++) { - if (ipiv[k] == 0) { - if (fabs(c[j + idem * k]) >= big) { - big = fabs(c[j + idem * k]); - irow = j; - icol = k; - } - } - } - } - } - ++(ipiv[icol]); - if (irow != icol) { - vcsUtil_swapRows(c, idem, n, b, m, irow, icol); - } - indxr[i] = irow; - indxc[i] = icol; - if (c[icol + idem * icol] == 0.0) { - plogf("vcsUtil_gaussj ERROR: Encountered a zero column: %d\n", i); - return 1; - } - pivinv = 1.0 / c[icol + idem * icol]; - c[icol + idem * icol] = 1.0; - for (l = 0; l < n; l++) { - c[icol + idem * l] *= pivinv; - } - for (l = 0; l < m; l++) { - b[icol + idem * l] *= pivinv; - } - for (ll = 0; ll < n; ll++) { - if (ll != icol) { - double dum = c[ll + idem * icol]; - c[ll + idem * icol] = 0; - for (l = 0; l < n; l++) { - c[ll + idem * l] -= c[icol + idem * l] * dum; - } - for (l = 0; l < m; l++) { - b[ll + idem * l] -= b[icol + idem * l] * dum; - } - } - } - } - if (needInverse) { - for (l = n-1; l != npos; l--) { - if (indxr[l] != indxc[l]) { - for (k = 0; k < n; k++) { - std::swap(c[k + idem * indxr[l]], c[k + idem * indxr[l]]); - } - } - } - } + std::vector indxc(n); + std::vector indxr(n); + std::vector ipiv(n, 0); + doublereal big = 0.0; + /* + * This is the main loop over the columns to be reduced. + */ + for (i = 0; i < n; i++) { + big = 0.0; + for (j = 0; j < n; j++) { + if (ipiv[j] != 1) { + for (k = 0; k < n; k++) { + if (ipiv[k] == 0) { + if (fabs(c[j + idem * k]) >= big) { + big = fabs(c[j + idem * k]); + irow = j; + icol = k; + } + } + } + } + } + ++(ipiv[icol]); + if (irow != icol) { + vcsUtil_swapRows(c, idem, n, b, m, irow, icol); + } + indxr[i] = irow; + indxc[i] = icol; + if (c[icol + idem * icol] == 0.0) { + plogf("vcsUtil_gaussj ERROR: Encountered a zero column: %d\n", i); + return 1; + } + pivinv = 1.0 / c[icol + idem * icol]; + c[icol + idem * icol] = 1.0; + for (l = 0; l < n; l++) { + c[icol + idem * l] *= pivinv; + } + for (l = 0; l < m; l++) { + b[icol + idem * l] *= pivinv; + } + for (ll = 0; ll < n; ll++) { + if (ll != icol) { + double dum = c[ll + idem * icol]; + c[ll + idem * icol] = 0; + for (l = 0; l < n; l++) { + c[ll + idem * l] -= c[icol + idem * l] * dum; + } + for (l = 0; l < m; l++) { + b[ll + idem * l] -= b[icol + idem * l] * dum; + } + } + } + } + if (needInverse) { + for (l = n - 1; l != npos; l--) { + if (indxr[l] != indxc[l]) { + for (k = 0; k < n; k++) { + std::swap(c[k + idem * indxr[l]], c[k + idem * indxr[l]]); + } + } + } + } - - /* - * The negative in the last expression is due to the form of B upon - * input - */ - for (i = 0; i < n; ++i) { - for (j = 0; j < m; ++j) { - b[i + j * idem] = -b[i + j * idem]; - } - } - return 0; + /* + * The negative in the last expression is due to the form of B upon + * input + */ + for (i = 0; i < n; ++i) { + for (j = 0; j < m; ++j) { + b[i + j * idem] = -b[i + j * idem]; + } + } + return 0; } //==================================================================================================================== @@ -674,36 +665,33 @@ int vcsUtil_gaussj(double* c, size_t idem, size_t n, double* b, size_t m) * - VCS_UNITS_KELVIN: 1.0 K-1 * - VCS_UNITS_MKS: joules kmol-1 K-1 = kg m2 s-2 kmol-1 K-1 */ -double vcsUtil_gasConstant(int mu_units) -{ - double r; - switch (mu_units) { - case VCS_UNITS_KCALMOL: - r = Cantera::GasConst_cal_mol_K * 1e-3; - break; - case VCS_UNITS_UNITLESS: - r = 1.0; - break; - case VCS_UNITS_KJMOL: - r = Cantera::GasConstant * 1e-6; - break; - case VCS_UNITS_KELVIN: - r = 1.0; - break; - case VCS_UNITS_MKS: - /* joules / kg-mol K = kg m2 / s2 kg-mol K */ - r = Cantera::GasConstant; - break; - default: - plogf("vcs_gasConstant error: uknown units: %d\n", - mu_units); - exit(EXIT_FAILURE); - } - return r; +double vcsUtil_gasConstant(int mu_units) { + double r; + switch (mu_units) { + case VCS_UNITS_KCALMOL: + r = Cantera::GasConst_cal_mol_K * 1e-3; + break; + case VCS_UNITS_UNITLESS: + r = 1.0; + break; + case VCS_UNITS_KJMOL: + r = Cantera::GasConstant * 1e-6; + break; + case VCS_UNITS_KELVIN: + r = 1.0; + break; + case VCS_UNITS_MKS: + /* joules / kg-mol K = kg m2 / s2 kg-mol K */ + r = Cantera::GasConstant; + break; + default: + plogf("vcs_gasConstant error: uknown units: %d\n", + mu_units); + exit(EXIT_FAILURE); + } + return r; } - - void vcs_print_line(const char* string, int num) /************************************************************************** @@ -714,86 +702,83 @@ void vcs_print_line(const char* string, int num) * ***************************************************************************/ { - if (string) { - for (int j = 0; j < num; j++) { - plogf("%s", string); - } - } - plogendl(); + if (string) { + for (int j = 0; j < num; j++) { + plogf("%s", string); + } + } + plogendl(); } - -const char* vcs_speciesType_string(int speciesStatus, int length) -{ - const char* sss; - switch (speciesStatus) { - case VCS_SPECIES_COMPONENT: - sss = "Component Species"; - break; - case VCS_SPECIES_MAJOR: - sss ="Major Species"; - break; - case VCS_SPECIES_MINOR: - sss ="Minor Species"; - break; - case VCS_SPECIES_ZEROEDPHASE: - if (length < 48) { - sss = "Set Zeroed-Phase"; - } else { - sss ="Purposely Zeroed-Phase Species (not in problem)"; - } - break; - case VCS_SPECIES_ZEROEDMS: - if (length < 23) { - sss = "Zeroed-MS Phase"; - } else { - sss ="Zeroed-MS Phase Species"; - } - break; - case VCS_SPECIES_ZEROEDSS: - if (length < 23) { - sss = "Zeroed-SS Phase"; - } else { - sss ="Zeroed-SS Phase Species"; - } - break; - case VCS_SPECIES_DELETED: - if (length < 22) { - sss = "Deleted Species"; - } else if (length < 40) { - sss = "Deleted-Small Species"; - } else { - sss ="Deleted-Small Species in a MS phase"; - } - break; - case VCS_SPECIES_ACTIVEBUTZERO: - if (length < 47) { - sss = "Tmp Zeroed in MS"; - } else { - sss ="Zeroed Species in an active MS phase (tmp)"; - } - break; - case VCS_SPECIES_STOICHZERO: - if (length < 56) { - sss = "Stoich Zeroed in MS"; - } else { - sss ="Zeroed Species in an active MS phase (Stoich Constraint)"; - } - break; - case VCS_SPECIES_INTERFACIALVOLTAGE: - if (length < 29) { - sss = "InterfaceVoltage"; - } else { - sss ="InterfaceVoltage Species"; - } - break; - default: - sss = "unknown species type"; - } - return sss; +const char* vcs_speciesType_string(int speciesStatus, int length) { + const char* sss; + switch (speciesStatus) { + case VCS_SPECIES_COMPONENT: + sss = "Component Species"; + break; + case VCS_SPECIES_MAJOR: + sss = "Major Species"; + break; + case VCS_SPECIES_MINOR: + sss = "Minor Species"; + break; + case VCS_SPECIES_ZEROEDPHASE: + if (length < 48) { + sss = "Set Zeroed-Phase"; + } else { + sss = "Purposely Zeroed-Phase Species (not in problem)"; + } + break; + case VCS_SPECIES_ZEROEDMS: + if (length < 23) { + sss = "Zeroed-MS Phase"; + } else { + sss = "Zeroed-MS Phase Species"; + } + break; + case VCS_SPECIES_ZEROEDSS: + if (length < 23) { + sss = "Zeroed-SS Phase"; + } else { + sss = "Zeroed-SS Phase Species"; + } + break; + case VCS_SPECIES_DELETED: + if (length < 22) { + sss = "Deleted Species"; + } else if (length < 40) { + sss = "Deleted-Small Species"; + } else { + sss = "Deleted-Small Species in a MS phase"; + } + break; + case VCS_SPECIES_ACTIVEBUTZERO: + if (length < 47) { + sss = "Tmp Zeroed in MS"; + } else { + sss = "Zeroed Species in an active MS phase (tmp)"; + } + break; + case VCS_SPECIES_STOICHZERO: + if (length < 56) { + sss = "Stoich Zeroed in MS"; + } else { + sss = "Zeroed Species in an active MS phase (Stoich Constraint)"; + } + break; + case VCS_SPECIES_INTERFACIALVOLTAGE: + if (length < 29) { + sss = "InterfaceVoltage"; + } else { + sss = "InterfaceVoltage Species"; + } + break; + default: + sss = "unknown species type"; + } + return sss; } - /************************************************************************ **/ void vcs_print_stringTrunc(const char* str, size_t space, int alignment) @@ -812,33 +797,33 @@ void vcs_print_stringTrunc(const char* str, size_t space, int alignment) * 2 left aligned ***********************************************************************/ { - size_t i, ls=0, rs=0; - size_t len = strlen(str); - if ((len) >= space) { - for (i = 0; i < space; i++) { - plogf("%c", str[i]); - } - } else { - if (alignment == 1) { - ls = space - len; - } else if (alignment == 2) { - rs = space - len; - } else { - ls = (space - len) / 2; - rs = space - len - ls; - } - if (ls != 0) { - for (i = 0; i < ls; i++) { - plogf(" "); - } - } - plogf("%s", str); - if (rs != 0) { - for (i = 0; i < rs; i++) { - plogf(" "); - } - } - } + size_t i, ls = 0, rs = 0; + size_t len = strlen(str); + if ((len) >= space) { + for (i = 0; i < space; i++) { + plogf("%c", str[i]); + } + } else { + if (alignment == 1) { + ls = space - len; + } else if (alignment == 2) { + rs = space - len; + } else { + ls = (space - len) / 2; + rs = space - len - ls; + } + if (ls != 0) { + for (i = 0; i < ls; i++) { + plogf(" "); + } + } + plogf("%s", str); + if (rs != 0) { + for (i = 0; i < rs; i++) { + plogf(" "); + } + } + } } /*****************************************************************************/ @@ -855,12 +840,81 @@ bool vcs_doubleEqual(double d1, double d2) * accuracy. *************************************************************************/ { - double denom = fabs(d1) + fabs(d2) + 1.0; - double fac = fabs(d1 - d2) / denom; - if (fac > 1.0E-10) { - return false; - } - return true; + double denom = fabs(d1) + fabs(d2) + 1.0; + double fac = fabs(d1 - d2) / denom; + if (fac > 1.0E-10) { + return false; + } + return true; } +//===================================================================================================================== +// Sorts a vector of ints in place from lowest to the highest values +/* + * The vector is returned sorted from lowest to highest. + * + * @param x Reference to a vector of ints. + */ +void vcs_heapsort(std::vector & x) { + int n = x.size(); + if (n < 2) + return; + doublereal rra; + int ll = n / 2; + int iret = n - 1; + + while (1 > 0) { + if (ll > 0) { + ll--; + rra = x[ll]; + } else { + rra = x[iret]; + x[iret] = x[0]; + iret--; + if (iret == 0) { + x[0] = rra; + return; + } + } + + int i = ll; + int j = ll + ll + 1; + + while (j <= iret) { + if (j < iret) { + if (x[j] < x[j + 1]) + j++; + } + if (rra < x[j]) { + x[i] = x[j]; + i = j; + j = j + j + 1; + } else { + j = iret + 1; + } + } + x[i] = rra; + } +} +//===================================================================================================================== +// Sorts a vector of ints and eliminates duplicates from the resulting list +/* + * @param xOrderedUnique Ordered vector of unique ints that were part of the original list + * @param x Reference to a constant vector of ints. + */ +void vcs_orderedUnique(std::vector & xOrderedUnique, const std::vector & x) +{ + std::vector xordered(x); + vcs_heapsort(xordered); + int lastV = x[0] - 1; + xOrderedUnique.clear(); + for (int i = 0; i < (int) xordered.size(); i++) { + if (lastV != xordered[i]) { + xOrderedUnique.push_back(xordered[i]); + lastV = xordered[i]; + } + } +} +//===================================================================================================================== + } diff --git a/src/kinetics/InterfaceKinetics.cpp b/src/kinetics/InterfaceKinetics.cpp index c191fe837..77dbad560 100644 --- a/src/kinetics/InterfaceKinetics.cpp +++ b/src/kinetics/InterfaceKinetics.cpp @@ -304,6 +304,8 @@ void InterfaceKinetics::getActivityConcentrations(doublereal* const conc) * Update the equilibrium constants in molar units for all * reversible reactions. Irreversible reactions have their * equilibrium constant set to zero. + * For reactions involving charged species the equilibrium + * constant is adjusted according to the electrostatic potential. */ void InterfaceKinetics::updateKc() { @@ -311,7 +313,10 @@ void InterfaceKinetics::updateKc() //static vector_fp mu(nTotalSpecies()); if (m_nrev > 0) { - + /* + * Get the vector of standard state electrochemical potentials for species in the Interfacial + * kinetics object and store it in m_mu0[] + */ size_t nsp, ik = 0; doublereal rt = GasConstant*thermo(0).temperature(); doublereal rrt = 1.0 / rt; @@ -327,8 +332,7 @@ void InterfaceKinetics::updateKc() } // compute Delta mu^0 for all reversible reactions - m_rxnstoich.getRevReactionDelta(m_ii, DATA_PTR(m_mu0), - DATA_PTR(m_rkcn)); + m_rxnstoich.getRevReactionDelta(m_ii, DATA_PTR(m_mu0), DATA_PTR(m_rkcn)); for (size_t i = 0; i < m_nrev; i++) { size_t irxn = m_revindex[i]; @@ -336,6 +340,7 @@ void InterfaceKinetics::updateKc() throw CanteraError("InterfaceKinetics", "illegal value: irxn = "+int2str(irxn)); } + // WARNING this may overflow HKM m_rkcn[irxn] = exp(m_rkcn[irxn]*rrt); } for (size_t i = 0; i != m_nirrev; ++i) { @@ -680,7 +685,7 @@ void InterfaceKinetics::updateROP() if (m_rxnPhaseIsReactant[j][rp]) { if (! m_phaseExists[rp]) { m_ropnet[j] = 0.0; - m_ropr[j] = m_ropr[j] = 0.0; + m_ropr[j] = m_ropf[j] = 0.0; } } } @@ -753,7 +758,7 @@ InterfaceKinetics::adjustRatesForIntermediatePhases() } #endif //================================================================================================= -//================================================================================================= + /* * * getDeltaGibbs(): @@ -1335,7 +1340,7 @@ solvePseudoSteadyStateProblem(int ifuncOverride, doublereal timeScaleOverride) } //================================================================================================ -void InterfaceKinetics::setPhaseExistence(const size_t iphase, const bool exists) +void InterfaceKinetics::setPhaseExistence(const size_t iphase, const int exists) { if (iphase >= m_thermo.size()) { throw CanteraError("InterfaceKinetics:setPhaseExistence", "out of bounds"); @@ -1343,6 +1348,9 @@ void InterfaceKinetics::setPhaseExistence(const size_t iphase, const bool exists if (exists) { if (!m_phaseExists[iphase]) { m_phaseExistsCheck--; + if (m_phaseExistsCheck < 0) { + m_phaseExistsCheck = 0; + } m_phaseExists[iphase] = true; } m_phaseIsStable[iphase] = true; @@ -1353,6 +1361,7 @@ void InterfaceKinetics::setPhaseExistence(const size_t iphase, const bool exists } m_phaseIsStable[iphase] = false; } + } //================================================================================================ // Gets the phase existence int for the ith phase