diff --git a/Cantera/src/equil/vcs_DoubleStarStar.cpp b/Cantera/src/equil/vcs_DoubleStarStar.cpp index 3d2034d59..4a95e111b 100644 --- a/Cantera/src/equil/vcs_DoubleStarStar.cpp +++ b/Cantera/src/equil/vcs_DoubleStarStar.cpp @@ -126,6 +126,10 @@ double * const * const DoubleStarStar::baseDataAddr() { return (double * const * const) &(m_colAddr[0]); } +double const * const * const DoubleStarStar::constBaseDataAddr() const { + return (double const * const * const) &(m_colAddr[0]); +} + // Number of rows int DoubleStarStar::nRows() const { return m_nrows; diff --git a/Cantera/src/equil/vcs_DoubleStarStar.h b/Cantera/src/equil/vcs_DoubleStarStar.h index bb4f55545..79db395ea 100644 --- a/Cantera/src/equil/vcs_DoubleStarStar.h +++ b/Cantera/src/equil/vcs_DoubleStarStar.h @@ -94,6 +94,14 @@ public: */ double * const * const baseDataAddr(); + //! Returns a const double ** pointer to the base address + /*! + * This is the second way to get to the data + * This returns a double ** which can later be used in + * Dmatrix[icol][irow] notation to get to the data + */ + double const * const * const constBaseDataAddr() const; + //! Number of rows int nRows() const; diff --git a/Cantera/src/equil/vcs_MultiPhaseEquil.cpp b/Cantera/src/equil/vcs_MultiPhaseEquil.cpp index cae3f8ad0..0fc2a75c3 100644 --- a/Cantera/src/equil/vcs_MultiPhaseEquil.cpp +++ b/Cantera/src/equil/vcs_MultiPhaseEquil.cpp @@ -987,7 +987,7 @@ namespace VCSnonideal { * Query the ThermoPhase object to find out what convention * it uses for the specification of activity and Standard State. */ - VolPhase->m_activityConvention = tPhase->activityConvention(); + VolPhase->p_activityConvention = tPhase->activityConvention(); /* * Assign the value of eqn of state * -> Handle conflicts here. @@ -1092,7 +1092,7 @@ namespace VCSnonideal { /* * Transfer the type of unknown */ - vprob->SpeciesUnknownType[kT] = VolPhase->m_speciesUnknownType[k]; + vprob->SpeciesUnknownType[kT] = VolPhase->speciesUnknownType(k); /* * Transfer the species information from the * volPhase structure to the VPROB structure diff --git a/Cantera/src/equil/vcs_VolPhase.cpp b/Cantera/src/equil/vcs_VolPhase.cpp index 9de383f69..b21a5574d 100644 --- a/Cantera/src/equil/vcs_VolPhase.cpp +++ b/Cantera/src/equil/vcs_VolPhase.cpp @@ -37,18 +37,17 @@ namespace VCSnonideal { m_singleSpecies(true), m_gasPhase(false), m_eqnState(VCS_EOS_CONSTANT), - m_numElemConstraints(0), ChargeNeutralityElement(-1), + p_VCS_UnitsFormat(VCS_UNITS_MKS), + p_activityConvention(0), + m_numElemConstraints(0), m_elemGlobalIndex(0), NVolSpecies(0), m_totalMolesInert(0.0), - m_activityConvention(0), m_isIdealSoln(false), m_existence(0), m_MFStartIndex(0), IndSpecies(0), - //IndSpeciesContig(true), - p_VCS_UnitsFormat(VCS_UNITS_MKS), m_useCanteraCalls(false), TP_ptr(0), v_totalMoles(0.0), @@ -99,15 +98,15 @@ namespace VCSnonideal { m_singleSpecies(b.m_singleSpecies), m_gasPhase(b.m_gasPhase), m_eqnState(b.m_eqnState), - m_numElemConstraints(b.m_numElemConstraints), ChargeNeutralityElement(b.ChargeNeutralityElement), + p_VCS_UnitsFormat(b.p_VCS_UnitsFormat), + p_activityConvention(b.p_activityConvention), + m_numElemConstraints(b.m_numElemConstraints), NVolSpecies(b.NVolSpecies), m_totalMolesInert(b.m_totalMolesInert), - m_activityConvention(b.m_activityConvention), m_isIdealSoln(b.m_isIdealSoln), m_existence(b.m_existence), m_MFStartIndex(b.m_MFStartIndex), - p_VCS_UnitsFormat(b.p_VCS_UnitsFormat), m_useCanteraCalls(b.m_useCanteraCalls), TP_ptr(b.TP_ptr), v_totalMoles(b.v_totalMoles), @@ -163,13 +162,13 @@ namespace VCSnonideal { m_elementNames[e] = b.m_elementNames[e]; } - ElActive = b.ElActive; + m_elementActive = b.m_elementActive; m_elementType = b.m_elementType; - FormulaMatrix.resize(m_numElemConstraints, NVolSpecies, 0.0); + m_formulaMatrix.resize(m_numElemConstraints, NVolSpecies, 0.0); for (int e = 0; e < m_numElemConstraints; e++) { for (int k = 0; k < NVolSpecies; k++) { - FormulaMatrix[e][k] = b.FormulaMatrix[e][k]; + m_formulaMatrix[e][k] = b.m_formulaMatrix[e][k]; } } @@ -178,7 +177,7 @@ namespace VCSnonideal { NVolSpecies = b.NVolSpecies; PhaseName = b.PhaseName; m_totalMolesInert = b.m_totalMolesInert; - m_activityConvention = b.m_activityConvention; + p_activityConvention= b.p_activityConvention; m_isIdealSoln = b.m_isIdealSoln; m_existence = b.m_existence; m_MFStartIndex = b.m_MFStartIndex; @@ -339,9 +338,9 @@ namespace VCSnonideal { m_elementNames.resize(numElemConstraints); - ElActive.resize(numElemConstraints+1, 1); + m_elementActive.resize(numElemConstraints+1, 1); m_elementType.resize(numElemConstraints, VCS_ELEM_TYPE_ABSPOS); - FormulaMatrix.resize(numElemConstraints, NVolSpecies, 0.0); + m_formulaMatrix.resize(numElemConstraints, NVolSpecies, 0.0); m_elementNames.resize(numElemConstraints, ""); m_elemGlobalIndex.resize(numElemConstraints, -1); @@ -1414,7 +1413,7 @@ namespace VCSnonideal { ename = tPhase->elementName(eT); if (ename == "E") { eFound = eT; - ElActive[eT] = 0; + m_elementActive[eT] = 0; m_elementType[eT] = VCS_ELEM_TYPE_ELECTRONCHARGE; } } @@ -1430,7 +1429,7 @@ namespace VCSnonideal { if (eFound == -2) { eFound = ne; m_elementType[ne] = VCS_ELEM_TYPE_ELECTRONCHARGE; - ElActive[ne] = 0; + m_elementActive[ne] = 0; std::string ename = "E"; m_elementNames[ne] = ename; ne++; @@ -1439,7 +1438,7 @@ namespace VCSnonideal { } - FormulaMatrix.resize(ne, ns, 0.0); + m_formulaMatrix.resize(ne, ns, 0.0); m_speciesUnknownType.resize(ns, VCS_SPECIES_TYPE_MOLNUM); @@ -1466,7 +1465,7 @@ namespace VCSnonideal { m_elementNames[e] = ename; } - double * const * const fm = FormulaMatrix.baseDataAddr(); + double * const * const fm = m_formulaMatrix.baseDataAddr(); for (k = 0; k < ns; k++) { e = 0; for (eT = 0; eT < nebase; eT++) { @@ -1516,5 +1515,18 @@ namespace VCSnonideal { void vcs_VolPhase::setElementType(const int e, const int eType) { m_elementType[e] = eType; } + + double const * const * const vcs_VolPhase::getFormulaMatrix() const { + double const * const * const fm = m_formulaMatrix.constBaseDataAddr(); + return fm; + } + + int vcs_VolPhase::speciesUnknownType(const int k) const { + return m_speciesUnknownType[k]; + } + + int vcs_VolPhase::elementActive(const int e) const { + return m_elementActive[e]; + } } diff --git a/Cantera/src/equil/vcs_VolPhase.h b/Cantera/src/equil/vcs_VolPhase.h index 4477dae8e..ee31f70a8 100644 --- a/Cantera/src/equil/vcs_VolPhase.h +++ b/Cantera/src/equil/vcs_VolPhase.h @@ -511,6 +511,29 @@ namespace VCSnonideal { */ int transferElementsFM(const Cantera::ThermoPhase * const tPhase); + //! Get a constant form of the Species Formula Matrix + /*! + * Returns a double ** pointer such that + * + * fm[e][f] is the formula matrix entry for element e for species k + */ + double const * const * const getFormulaMatrix() const; + + //! Returns the type of the species unknown + /*! + * @param k species index + * + * returns the SpeciesUnknownType[k] = type of species + * Normal -> VCS_SPECIES_TYPE_MOLUNK + * ( unknown is the mole number in the phase) + * metal electron -> VCS_SPECIES_INTERFACIALVOLTAGE + * ( unknown is the interfacial voltage (volts) + */ + int speciesUnknownType(const int k) const; + + + int elementActive(const int e) const; + private: //! Evaluate the activity coefficients at the current conditions @@ -580,6 +603,7 @@ namespace VCSnonideal { */ void _updateMoleFractionDependencies(); + /************************************************************************* * MEMBER DATA * @@ -627,14 +651,6 @@ namespace VCSnonideal { */ int m_eqnState; - private: - //! Number of element constraints within the problem - /*! - * This is usually equal to the number of elements. - */ - int m_numElemConstraints; - - public: //! This is the element number for the charge neutrality //! condition of the phase /*! @@ -643,19 +659,58 @@ namespace VCSnonideal { */ int ChargeNeutralityElement; + //! Units for the chemical potential data, pressure data, volume, + //! and species amounts + /*! + * All internally storred quantities will have these units. Also, printed + * quantitities will display in these units. Input quantities are expected + * in these units. + * + * Chem_Pot Pres vol moles + * ---------------------------------------------------------------------- + * -1 VCS_UNITS_KCALMOL = kcal/gmol Pa m**3 kmol + * 0 VCS_UNITS_UNITLESS = MU / RT -> no units Pa m**3 kmol + * 1 VCS_UNITS_KJMOL = kJ / gmol Pa m**3 kmol + * 2 VCS_UNITS_KELVIN = KELVIN -> MU / R Pa m**3 kmol + * 3 VCS_UNITS_MKS = Joules / Kmol (Cantera) Pa m**3 kmol + * ---------------------------------------------------------------------- + * + * see vcs_defs.h for more information. + * + * Currently, this value should be the same as the owning VCS_PROB or + * VCS_SOLVE object. There is no code for handling anything else atm. + * + * (This variable is needed for the vcsc code, where it is not equal + * to VCS_UNITS_MKS). + */ + int p_VCS_UnitsFormat; + + //! Convention for the activity formulation + /*! + * 0 = molar based activities (default) + * 1 = Molality based activities + * mu = mu_0 + ln a_molality + * standard state is based on unity molality + */ + int p_activityConvention; + private: + //! Number of element constraints within the problem + /*! + * This is usually equal to the number of elements. + */ + int m_numElemConstraints; + //! vector of strings containing the element constraint names /*! * Length = nElemConstraints */ std::vector m_elementNames; - public: //! boolean indicating whether an element constraint is active //! for the current problem - std::vector ElActive; + std::vector m_elementActive; - private: //! Type of the element constraint /*! * m_elType[j] = type of the element @@ -669,7 +724,6 @@ namespace VCSnonideal { */ std::vector m_elementType; - public: //! Formula Matrix for the phase /*! * FormulaMatrix[j][kspec] @@ -677,7 +731,7 @@ namespace VCSnonideal { * Number of elements, j, * in the kspec species */ - DoubleStarStar FormulaMatrix; + DoubleStarStar m_formulaMatrix; //! Type of the species unknown /*! @@ -689,8 +743,7 @@ namespace VCSnonideal { */ std::vector m_speciesUnknownType; - private: - //! Index of the element number in the global list of elements + //! Index of the element number in the global list of elements //! storred in VCS_PROB or VCS_SOLVE std::vector m_elemGlobalIndex; @@ -705,22 +758,11 @@ namespace VCSnonideal { //! Total moles of inert in the phase double m_totalMolesInert; - public: - //! Convention for the activity formulation - /*! - * 0 = molar based activities (default) - * 1 = Molality based activities - * mu = mu_0 + ln a_molality - * standard state is based on unity molality - */ - int m_activityConvention; - //! Boolean indicating whether the phase is an ideal solution //! and therefore it's molar-based activity coefficients are //! uniformly equal to one. bool m_isIdealSoln; - private: //! Current state of existence: /*! * 0 : Doesn't exist currently @@ -755,34 +797,6 @@ namespace VCSnonideal { */ std::vector ListSpeciesPtr; - public: - //! Units for the chemical potential data, pressure data, volume, - //! and species amounts - /*! - * All internally storred quantities will have these units. Also, printed - * quantitities will display in these units. Input quantities are expected - * in these units. - * - * Chem_Pot Pres vol moles - * ---------------------------------------------------------------------- - * -1 VCS_UNITS_KCALMOL = kcal/gmol Pa m**3 kmol - * 0 VCS_UNITS_UNITLESS = MU / RT -> no units Pa m**3 kmol - * 1 VCS_UNITS_KJMOL = kJ / gmol Pa m**3 kmol - * 2 VCS_UNITS_KELVIN = KELVIN -> MU / R Pa m**3 kmol - * 3 VCS_UNITS_MKS = Joules / Kmol (Cantera) Pa m**3 kmol - * ---------------------------------------------------------------------- - * - * see vcs_defs.h for more information. - * - * Currently, this value should be the same as the owning VCS_PROB or - * VCS_SOLVE object. There is no code for handling anything else atm. - * - * (This variable is needed for the vcsc code, where it is not equal - * to VCS_UNITS_MKS). - */ - int p_VCS_UnitsFormat; - - private: //! If this is true, then calculations are actually performed within //! Cantera bool m_useCanteraCalls; diff --git a/Cantera/src/equil/vcs_prob.cpp b/Cantera/src/equil/vcs_prob.cpp index 3e8952bed..b418ac837 100644 --- a/Cantera/src/equil/vcs_prob.cpp +++ b/Cantera/src/equil/vcs_prob.cpp @@ -384,7 +384,7 @@ namespace VCSnonideal { } if (foundPos == -1) { int elType = volPhase->elementType(eVP); - int elactive = volPhase->ElActive[eVP]; + int elactive = volPhase->elementActive(eVP); e = addElement(enVP.c_str(), elType, elactive); volPhase->setElemGlobalIndex(eVP, e); } @@ -442,7 +442,7 @@ namespace VCSnonideal { plogf("Shouldn't be here\n"); exit(-1); } - double *const *const fm = volPhase->FormulaMatrix.baseDataAddr(); + double const *const *const fm = volPhase->getFormulaMatrix(); for (eVP = 0; eVP < volPhase->nElemConstraints(); eVP++) { e = volPhase->elemGlobalIndex(eVP); #ifdef DEBUG_MODE diff --git a/Cantera/src/equil/vcs_solve.cpp b/Cantera/src/equil/vcs_solve.cpp index 3ce0a2ee6..d7474ac3a 100644 --- a/Cantera/src/equil/vcs_solve.cpp +++ b/Cantera/src/equil/vcs_solve.cpp @@ -736,8 +736,8 @@ namespace VCSnonideal { */ for (iph = 0; iph < nph; iph++) { Vphase = m_VolPhaseList[iph]; - m_phaseActConvention[iph] = Vphase->m_activityConvention; - if (Vphase->m_activityConvention != 0) { + m_phaseActConvention[iph] = Vphase->p_activityConvention; + if (Vphase->p_activityConvention != 0) { /* * We assume here that species 0 is the solvent. * The solvent isn't on a unity activity basis @@ -751,7 +751,7 @@ namespace VCSnonideal { double mnaught = m_wtSpecies[iSolvent] / 1000.; for (int k = 1; k < Vphase->NVolSpecies; k++) { int kspec = Vphase->spGlobalIndexVCS(k); - m_actConventionSpecies[kspec] = Vphase->m_activityConvention; + m_actConventionSpecies[kspec] = Vphase->p_activityConvention; m_lnMnaughtSpecies[kspec] = log(mnaught); } } @@ -983,7 +983,7 @@ namespace VCSnonideal { pub->mf[kT], vPhase->molefraction(k)); exit(-1); } - if (pubPhase->m_speciesUnknownType[k] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { + if (pubPhase->speciesUnknownType(k) != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { sumMoles += pub->w[kT]; } }