diff --git a/include/cantera/equil/vcs_DoubleStarStar.h b/include/cantera/equil/vcs_DoubleStarStar.h deleted file mode 100644 index b7d1f661c..000000000 --- a/include/cantera/equil/vcs_DoubleStarStar.h +++ /dev/null @@ -1,114 +0,0 @@ -/** - * @file vcs_DoubleStarStar.h - * - * Header file for class DoubleStarStar - */ -#ifndef VCS_DOUBLESTARSTAR_H -#define VCS_DOUBLESTARSTAR_H - -#include - -namespace VCSnonideal -{ - -using std::size_t; - -//! A class for 2D double arrays stored in column-major -//! (Fortran-compatible) form. -/*! - * In this form, the data entry for an `n` row, `m` colum matrix is index = - * `i + (n-1) * j` where `Matrix[j][i]` references the element in row `i`, - * column `j`. - * - * The way this is instantiated is via the constructor, - * DoubleStarStar Dmatrix(mcol, mrow)`. - * - * The way this is referenced is via the notation: `Dmatrix[icol][irow]`. - */ -class DoubleStarStar -{ -public: - - //! Default constructor. Create an empty array. - DoubleStarStar(); - - //! Constructor. - /*! - * Create an `nrow` by `mcol` double array, and initialize all elements - * to `v`. - * - * @param mcol Number of columns - * @param nrow Number of rows - * @param v value used to initialize elements - */ - DoubleStarStar(size_t mcol, size_t nrow, double v = 0.0); - - DoubleStarStar(const DoubleStarStar& y); - DoubleStarStar& operator=(const DoubleStarStar& y); - - //! Resize the array, and fill the new entries with `v` - /*! - * @param mcol This is the number of columns in the new matrix - * @param nrow This is the number of rows - * @param v Default fill value -> defaults to zero. - */ - void resize(size_t mcol, size_t nrow, double v = 0.0); - - //! Pointer to the top of the column - /*! - * @param jcol This is the jth column - * - * @return returns the pointer to the top of the jth column - */ - double* operator[](size_t jcol); - - //! Returns a const Pointer to the top of the jth column - /*! - * @param jcol This is the jth column - * - * @return returns the pointer to the top of the jth column - */ - const double* operator[](size_t jcol) const; - - //! Returns a `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* 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* constBaseDataAddr() const; - - //! Number of rows - size_t nRows() const; - - //! Number of columns - size_t nColumns() const; - -private: - //! Storage area - std::vector m_data; - - //! Vector of addresses for the top of the columns - /*! - * Length = mcol - */ - std::vector m_colAddr; - - //! number of rows - size_t m_nrows; - - //! number of columns - size_t m_ncols; -}; - -} - -#endif diff --git a/include/cantera/equil/vcs_VolPhase.h b/include/cantera/equil/vcs_VolPhase.h index 147d78efe..68b6221fb 100644 --- a/include/cantera/equil/vcs_VolPhase.h +++ b/include/cantera/equil/vcs_VolPhase.h @@ -11,8 +11,8 @@ #ifndef VCS_VOLPHASE_H #define VCS_VOLPHASE_H -#include "cantera/equil/vcs_DoubleStarStar.h" #include "cantera/equil/vcs_SpeciesProperties.h" +#include "cantera/base/Array.h" // Forward reference for ThermoPhase object within the Cantera namespace namespace Cantera @@ -288,12 +288,12 @@ public: * * @param LnAcJac_VCS jacobian parameter * The Jacobians are actually d( lnActCoeff) / d (MolNumber); - * dLnActCoeffdMolNumber[j][k] + * dLnActCoeffdMolNumber(k,j) * * j = id of the species mole number * k = id of the species activity coefficient */ - void sendToVCS_LnActCoeffJac(double* const* const LnACJac_VCS); + void sendToVCS_LnActCoeffJac(Cantera::Array2D& LnACJac_VCS); //! Set the pointer for Cantera's ThermoPhase parameter /*! @@ -502,7 +502,7 @@ public: * Returns a `double**` pointer such that `fm[e][f]` is the formula * matrix entry for element `e` for species `k` */ - double const* const* getFormulaMatrix() const; + const Cantera::Array2D& getFormulaMatrix() const; //! Returns the type of the species unknown /*! @@ -680,10 +680,10 @@ private: //! Formula Matrix for the phase /*! - * FormulaMatrix[j][kspec] = Formula Matrix for the species + * FormulaMatrix(kspec,j) = Formula Matrix for the species * Number of elements, j, in the kspec species */ - DoubleStarStar m_formulaMatrix; + Cantera::Array2D m_formulaMatrix; //! Type of the species unknown /*! @@ -838,11 +838,11 @@ private: //! Vector of the derivatives of the ln activity coefficient wrt to the //! current mole number multiplied by the current phase moles /*! - * np_dLnActCoeffdMolNumber[j][k]; + * np_dLnActCoeffdMolNumber(k,j); * - j = id of the species mole number * - k = id of the species activity coefficient */ - mutable DoubleStarStar np_dLnActCoeffdMolNumber; + mutable Cantera::Array2D np_dLnActCoeffdMolNumber; //! Status /*! diff --git a/include/cantera/equil/vcs_prob.h b/include/cantera/equil/vcs_prob.h index a6bf17b06..8456bdb0a 100644 --- a/include/cantera/equil/vcs_prob.h +++ b/include/cantera/equil/vcs_prob.h @@ -11,7 +11,7 @@ #ifndef _VCS_PROB_H #define _VCS_PROB_H -#include "vcs_DoubleStarStar.h" +#include "cantera/base/Array.h" #include "vcs_IntStarStar.h" #include "cantera/equil/vcs_defs.h" #include @@ -88,9 +88,9 @@ public: //! Formula Matrix for the problem /*! - * FormulaMatrix[j][kspec] = Number of elements, j, in the kspec species + * FormulaMatrix(kspec,j) = Number of elements, j, in the kspec species */ - DoubleStarStar FormulaMatrix; + Cantera::Array2D FormulaMatrix; //! Specifies the species unknown type /*! diff --git a/include/cantera/equil/vcs_solve.h b/include/cantera/equil/vcs_solve.h index c221a4614..f48fc9843 100644 --- a/include/cantera/equil/vcs_solve.h +++ b/include/cantera/equil/vcs_solve.h @@ -25,9 +25,9 @@ #include "cantera/base/ct_defs.h" #include "cantera/equil/vcs_defs.h" -#include "cantera/equil/vcs_DoubleStarStar.h" #include "cantera/equil/vcs_IntStarStar.h" #include "cantera/equil/vcs_internal.h" +#include "cantera/base/Array.h" namespace VCSnonideal { @@ -163,13 +163,13 @@ public: * Rearranges the solution data to put the component data at the * front of the species list. * - * Then, calculates m_stoichCoeffRxnMatrix[irxn][jcomp] the formation + * Then, calculates m_stoichCoeffRxnMatrix(jcomp,irxn) the formation * reactions for all noncomponent species in the mechanism. Also * calculates DNG(I) and DNL(I), the net mole change for each formation * reaction. Also, initializes IR(I) to the default state. * - * @param[in] doJustComponents If true, the m_stoichCoeffRxnMatrix[][] and - * m_deltaMolNumPhase[] are not calculated. + * @param[in] doJustComponents If true, the m_stoichCoeffRxnMatrix and + * m_deltaMolNumPhase are not calculated. * * @param[in] aw Vector of mole fractions which will be used to construct an * optimal basis from. @@ -193,14 +193,11 @@ public: * calculates the #m_numComponents species. It switches their positions * in the species vector so that they occupy the first #m_numComponents * spots in the species vector. - * - #m_stoichCoeffRxnMatrix[irxn][jcomp] Stoichiometric coefficient + * - #m_stoichCoeffRxnMatrix(jcomp,irxn) Stoichiometric coefficient * matrix for the reaction mechanism expressed in Reduced Canonical * Form. jcomp refers to the component number, and irxn refers to the * irxn_th non-component species. - * - #m_deltaMolNumPhase[irxn]: Change in the number of total number of - * moles of species in all phases due to the noncomponent formation - * reaction, irxn. - * - #m_deltaMolNumPhase[irxn][iphase]: Change in the number of moles in + * - #m_deltaMolNumPhase(iphase,irxn): Change in the number of moles in * phase, iphase, due to the noncomponent formation reaction, irxn. * - #m_phaseParticipation[irxn]: This is 1 if the phase, iphase, * participates in the formation reaction, irxn, and zero otherwise. @@ -547,7 +544,7 @@ public: /*! * Formation reactions are * reactions which create each noncomponent species from the component - * species. m_stoichCoeffRxnMatrix[irxn][jcomp] are the stoichiometric + * species. m_stoichCoeffRxnMatrix(jcomp,irxn) are the stoichiometric * coefficients for these reactions. A stoichiometric coefficient of * one is assumed for species irxn in this reaction. * @@ -1331,20 +1328,6 @@ private: */ bool vcs_globStepDamp(); - //! Switch rows and columns of a square matrix - /*! - * Switches the row and column of a matrix. So that after - * - * J[k1][j] = J_old[k2][j] and J[j][k1] = J_old[j][k2] - * J[k2][j] = J_old[k1][j] and J[j][k2] = J_old[j][k1] - * - * @param Jac Double pointer to the Jacobian - * @param k1 first row/column value to be switched - * @param k2 second row/column value to be switched - */ - void vcs_switch2D(double* const* const Jac, - const size_t k1, const size_t k2) const; - //! Calculate the norm of a deltaGibbs free energy vector /*! * Positive DG for species which don't exist are ignored. @@ -1545,11 +1528,11 @@ public: //! Formula matrix for the problem /*! - * FormulaMatrix[j][kspec] = Number of elements, j, in the kspec species + * FormulaMatrix(kspec,j) = Number of elements, j, in the kspec species * * Both element and species indices are swapped. */ - DoubleStarStar m_formulaMatrix; + Cantera::Array2D m_formulaMatrix; //! Stoichiometric coefficient matrix for the reaction mechanism expressed in Reduced Canonical Form. /*! @@ -1559,15 +1542,15 @@ public: * * NOTE: kspec = irxn + m_numComponents * - * m_stoichCoeffRxnMatrix[irxn][j] : + * m_stoichCoeffRxnMatrix(j,irxn) : * j refers to the component number, and irxn refers to the irxn_th non-component species. * The stoichiometric coefficients multilplied by the Formula coefficients of the * component species add up to the negative value of the number of elements in * the species kspec. * - * length = [nspecies0][nelements0] + * size = nelements0 x nspecies0 */ - DoubleStarStar m_stoichCoeffRxnMatrix; + Cantera::Array2D m_stoichCoeffRxnMatrix; //! Absolute size of the stoichiometric coefficients /*! @@ -1645,9 +1628,9 @@ public: //! Change in the number of moles of phase, iphase, due to the //! noncomponent formation reaction, irxn, for species, k: /*! - * m_deltaMolNumPhase[irxn][iphase] = k = nc + irxn + * m_deltaMolNumPhase(iphase,irxn) = k = nc + irxn */ - DoubleStarStar m_deltaMolNumPhase; + Cantera::Array2D m_deltaMolNumPhase; //! This is 1 if the phase, iphase, participates in the formation reaction //! irxn, and zero otherwise. PhaseParticipation[irxn][iphase] @@ -1944,12 +1927,12 @@ public: //! Change in the log of the activity coefficient with respect to the mole number //! multiplied by the phase mole number /*! - * length = [nspecies][nspecies] + * size = nspecies x nspecies * * This is a temporary array that gets regenerated every time it's * needed. It is not swapped wrt species. */ - DoubleStarStar m_np_dLnActCoeffdMolNum; + Cantera::Array2D m_np_dLnActCoeffdMolNum; //! Molecular weight of each species /*! diff --git a/src/equil/vcs_DoubleStarStar.cpp b/src/equil/vcs_DoubleStarStar.cpp deleted file mode 100644 index ead6bb44f..000000000 --- a/src/equil/vcs_DoubleStarStar.cpp +++ /dev/null @@ -1,138 +0,0 @@ -/** - * @file vcs_DoubleStarStar.cpp - * - * Implementation file for class DoubleStarStar - */ -#include "cantera/equil/vcs_DoubleStarStar.h" - -namespace VCSnonideal -{ - -DoubleStarStar::DoubleStarStar() : - m_nrows(0), - m_ncols(0) -{ - m_data.clear(); - m_colAddr.clear(); -} - -DoubleStarStar::DoubleStarStar(size_t m, size_t n, double v) : - m_nrows(n), - m_ncols(m) -{ - m_data.resize(n*m); - std::fill(m_data.begin(), m_data.end(), v); - m_colAddr.resize(m); - for (size_t jcol = 0; jcol < m_ncols; jcol++) { - m_colAddr[jcol] = &(m_data[jcol*m_nrows]); - } -} - -DoubleStarStar::DoubleStarStar(const DoubleStarStar& y) -{ - m_nrows = y.m_nrows; - m_ncols = y.m_ncols; - m_data.resize(m_nrows*m_ncols); - m_data = y.m_data; - m_colAddr.resize(m_ncols); - if (!m_data.empty()) { - for (size_t jcol = 0; jcol < m_ncols; jcol++) { - m_colAddr[jcol] = &(m_data[jcol*m_nrows]); - } - } -} - -DoubleStarStar& DoubleStarStar::operator=(const DoubleStarStar& y) -{ - if (&y == this) { - return *this; - } - m_nrows = y.m_nrows; - m_ncols = y.m_ncols; - m_data.resize(m_nrows*m_ncols); - m_data = y.m_data; - m_colAddr.resize(m_ncols); - if (!m_data.empty()) { - for (size_t jcol = 0; jcol < m_ncols; jcol++) { - m_colAddr[jcol] = &(m_data[jcol*m_nrows]); - } - } - return *this; -} - -void DoubleStarStar::resize(size_t m, size_t n, double v) -{ - std::vector old_data; - bool doCopy = false; - if (m_nrows > 0 && m_ncols > 0) { - if (m_nrows != n) { - doCopy = true; - old_data = m_data; - } - } - m_data.resize(n*m, v); - if (doCopy) { - if (n >= m_nrows && m >= m_ncols) { - for (size_t jcol = 0; jcol < m_ncols; jcol++) { - for (size_t irow = 0; irow < m_nrows; irow++) { - m_data[jcol*n + irow] = old_data[jcol*m_nrows + irow]; - } - for (size_t irow = m_nrows; irow < n; irow++) { - m_data[jcol*n + irow] = v; - } - } - for (size_t jcol = m_ncols; jcol < m; jcol++) { - for (size_t irow = 0; irow < n; irow++) { - m_data[jcol*n + irow] = v; - } - } - } else { - std::fill(m_data.begin(), m_data.end(), v); - for (size_t jcol = 0; jcol < m_ncols; jcol++) { - for (size_t irow = 0; irow < m_nrows; irow++) { - m_data[jcol*n + irow] = old_data[jcol*m_nrows + irow]; - } - } - } - } - m_nrows = n; - m_ncols = m; - m_colAddr.resize(m_ncols); - if (!m_data.empty()) { - for (size_t jcol = 0; jcol < m_ncols; jcol++) { - m_colAddr[jcol] = &(m_data[jcol*m_nrows]); - } - } -} - -double* DoubleStarStar::operator[](size_t jcol) -{ - return m_colAddr[jcol]; -} - -const double* DoubleStarStar::operator[](size_t jcol) const -{ - return (const double*) m_colAddr[jcol]; -} - -double* const* DoubleStarStar::baseDataAddr() -{ - return (double* const*) &(m_colAddr[0]); -} - -double const* const* DoubleStarStar::constBaseDataAddr() const -{ - return (double const* const*) &(m_colAddr[0]); -} - -size_t DoubleStarStar::nRows() const -{ - return m_nrows; -} - -size_t DoubleStarStar::nColumns() const -{ - return m_ncols; -} - -} diff --git a/src/equil/vcs_MultiPhaseEquil.cpp b/src/equil/vcs_MultiPhaseEquil.cpp index 0d4330710..152dec054 100644 --- a/src/equil/vcs_MultiPhaseEquil.cpp +++ b/src/equil/vcs_MultiPhaseEquil.cpp @@ -987,7 +987,7 @@ int vcs_Cantera_to_vprob(Cantera::MultiPhase* mphase, sProp->WtSpecies = tPhase->molecularWeight(k); sProp->FormulaMatrixCol.resize(vprob->ne, 0.0); for (size_t e = 0; e < vprob->ne; e++) { - sProp->FormulaMatrixCol[e] = vprob->FormulaMatrix[e][kT]; + sProp->FormulaMatrixCol[e] = vprob->FormulaMatrix(kT,e); } sProp->Charge = tPhase->charge(k); sProp->SurfaceSpecies = false; @@ -1291,8 +1291,6 @@ void vcs_MultiPhaseEquil::getStoichVector(size_t rxn, Cantera::vector_fp& nu) nu[i] = 0.0; } size_t nc = numComponents(); - // scMatrix [nrxn][ncomp] - const DoubleStarStar& scMatrix = m_vsolve.m_stoichCoeffRxnMatrix; const std::vector& indSpecies = m_vsolve.m_speciesMapIndex; if (rxn > nsp - nc) { return; @@ -1301,7 +1299,7 @@ void vcs_MultiPhaseEquil::getStoichVector(size_t rxn, Cantera::vector_fp& nu) nu[j] = 1.0; for (size_t kc = 0; kc < nc; kc++) { j = indSpecies[kc]; - nu[j] = scMatrix[rxn][kc]; + nu[j] = m_vsolve.m_stoichCoeffRxnMatrix(kc,rxn); } } diff --git a/src/equil/vcs_VolPhase.cpp b/src/equil/vcs_VolPhase.cpp index 312608906..15aeec0a7 100644 --- a/src/equil/vcs_VolPhase.cpp +++ b/src/equil/vcs_VolPhase.cpp @@ -131,12 +131,7 @@ vcs_VolPhase& vcs_VolPhase::operator=(const vcs_VolPhase& b) } m_elementActive = b.m_elementActive; m_elementType = b.m_elementType; - m_formulaMatrix.resize(m_numElemConstraints, m_numSpecies, 0.0); - for (size_t e = 0; e < m_numElemConstraints; e++) { - for (size_t k = 0; k < m_numSpecies; k++) { - m_formulaMatrix[e][k] = b.m_formulaMatrix[e][k]; - } - } + m_formulaMatrix = b.m_formulaMatrix; m_speciesUnknownType = b.m_speciesUnknownType; m_elemGlobalIndex = b.m_elemGlobalIndex; PhaseName = b.PhaseName; @@ -299,7 +294,7 @@ void vcs_VolPhase::elemResize(const size_t numElemConstraints) m_elementActive.resize(numElemConstraints+1, 1); m_elementType.resize(numElemConstraints, VCS_ELEM_TYPE_ABSPOS); - m_formulaMatrix.resize(numElemConstraints, m_numSpecies, 0.0); + m_formulaMatrix.resize(m_numSpecies, numElemConstraints, 0.0); m_elementNames.resize(numElemConstraints, ""); m_elemGlobalIndex.resize(numElemConstraints, npos); @@ -750,10 +745,10 @@ void vcs_VolPhase::_updateLnActCoeffJac() if (!TP_ptr) { return; } - TP_ptr->getdlnActCoeffdlnN(m_numSpecies, &np_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 = phaseTotalMoles * Xmol_[j]; - double* const np_lnActCoeffCol = np_dLnActCoeffdMolNumber[j]; + double* const np_lnActCoeffCol = np_dLnActCoeffdMolNumber.ptrColumn(j); if (moles_j_base < 1.0E-200) { moles_j_base = 1.0E-7 * moles_j_base + 1.0E-13 * phaseTotalMoles + 1.0E-150; } @@ -812,7 +807,7 @@ void vcs_VolPhase::_updateLnActCoeffJac() _updateActCoeff(); } -void vcs_VolPhase::sendToVCS_LnActCoeffJac(double* const* const np_LnACJac_VCS) +void vcs_VolPhase::sendToVCS_LnActCoeffJac(Cantera::Array2D& np_LnACJac_VCS) { /* * update the Ln Act Coeff jacobian entries with respect to the @@ -826,11 +821,9 @@ void vcs_VolPhase::sendToVCS_LnActCoeffJac(double* const* const np_LnACJac_VCS) */ for (size_t j = 0; j < m_numSpecies; j++) { size_t jglob = IndSpecies[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]; - np_lnACJacVCS_col[kglob] = np_lnACJac_col[k]; + np_LnACJac_VCS(kglob,jglob) = np_dLnActCoeffdMolNumber(k,j); } } } @@ -1217,7 +1210,7 @@ size_t vcs_VolPhase::transferElementsFM(const Cantera::ThermoPhase* const tPhase } - m_formulaMatrix.resize(ne, ns, 0.0); + m_formulaMatrix.resize(ns, ne, 0.0); m_speciesUnknownType.resize(ns, VCS_SPECIES_TYPE_MOLNUM); @@ -1241,21 +1234,20 @@ size_t vcs_VolPhase::transferElementsFM(const Cantera::ThermoPhase* const tPhase m_elementNames[e] = "cn_" + pname; } - double* const* const fm = m_formulaMatrix.baseDataAddr(); for (size_t k = 0; k < ns; k++) { e = 0; for (size_t eT = 0; eT < nebase; eT++) { - fm[e][k] = tPhase->nAtoms(k, eT); + m_formulaMatrix(k,e) = tPhase->nAtoms(k, eT); e++; } if (eFound != npos) { - fm[eFound][k] = - tPhase->charge(k); + m_formulaMatrix(k,eFound) = - tPhase->charge(k); } } if (cne) { for (size_t k = 0; k < ns; k++) { - fm[ChargeNeutralityElement][k] = tPhase->charge(k); + m_formulaMatrix(k,ChargeNeutralityElement) = tPhase->charge(k); } } @@ -1285,9 +1277,9 @@ void vcs_VolPhase::setElementType(const size_t e, const int eType) m_elementType[e] = eType; } -double const* const* vcs_VolPhase::getFormulaMatrix() const +const Cantera::Array2D& vcs_VolPhase::getFormulaMatrix() const { - return m_formulaMatrix.constBaseDataAddr(); + return m_formulaMatrix; } int vcs_VolPhase::speciesUnknownType(const size_t k) const diff --git a/src/equil/vcs_elem.cpp b/src/equil/vcs_elem.cpp index 4146eb298..cd9935d3c 100644 --- a/src/equil/vcs_elem.cpp +++ b/src/equil/vcs_elem.cpp @@ -21,7 +21,7 @@ void VCS_SOLVE::vcs_elab() m_elemAbundances[j] = 0.0; for (size_t i = 0; i < m_numSpeciesTot; ++i) { if (m_speciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - m_elemAbundances[j] += m_formulaMatrix[j][i] * m_molNumSpecies_old[i]; + m_elemAbundances[j] += m_formulaMatrix(i,j) * m_molNumSpecies_old[i]; } } } @@ -57,7 +57,7 @@ bool VCS_SOLVE::vcs_elabcheck(int ibound) */ bool multisign = false; for (size_t kspec = 0; kspec < m_numSpeciesTot; kspec++) { - double eval = m_formulaMatrix[i][kspec]; + double eval = m_formulaMatrix(kspec,i); if (eval < 0.0) { multisign = true; } @@ -98,7 +98,7 @@ void VCS_SOLVE::vcs_elabPhase(size_t iphase, double* const elemAbundPhase) for (size_t i = 0; i < m_numSpeciesTot; ++i) { if (m_speciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { if (m_phaseID[i] == iphase) { - elemAbundPhase[j] += m_formulaMatrix[j][i] * m_molNumSpecies_old[i]; + elemAbundPhase[j] += m_formulaMatrix(i,j) * m_molNumSpecies_old[i]; } } } @@ -142,7 +142,7 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) bool multisign = false; for (size_t kspec = 0; kspec < m_numSpeciesTot; kspec++) { if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - double eval = m_formulaMatrix[i][kspec]; + double eval = m_formulaMatrix(kspec,i); if (eval < 0.0) { multisign = true; } @@ -155,7 +155,7 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) if (numNonZero < 2) { for (size_t kspec = 0; kspec < m_numSpeciesTot; kspec++) { if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - double eval = m_formulaMatrix[i][kspec]; + double eval = m_formulaMatrix(kspec,i); if (eval > 0.0) { m_molNumSpecies_old[kspec] = m_elemAbundancesGoal[i] / eval; changed = true; @@ -167,7 +167,7 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) size_t compID = npos; for (size_t kspec = 0; kspec < m_numComponents; kspec++) { if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - double eval = m_formulaMatrix[i][kspec]; + double eval = m_formulaMatrix(kspec,i); if (eval > 0.0) { compID = kspec; numCompNonZero++; @@ -178,10 +178,10 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) double diff = m_elemAbundancesGoal[i]; for (size_t kspec = m_numComponents; kspec < m_numSpeciesTot; kspec++) { if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - double eval = m_formulaMatrix[i][kspec]; + double eval = m_formulaMatrix(kspec,i); diff -= eval * m_molNumSpecies_old[kspec]; } - m_molNumSpecies_old[compID] = std::max(0.0,diff/m_formulaMatrix[i][compID]); + m_molNumSpecies_old[compID] = std::max(0.0,diff/m_formulaMatrix(compID,i)); changed = true; } } @@ -207,7 +207,7 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) if (elType == VCS_ELEM_TYPE_ABSPOS) { for (size_t kspec = 0; kspec < m_numSpeciesTot; kspec++) { if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - double atomComp = m_formulaMatrix[i][kspec]; + double atomComp = m_formulaMatrix(kspec,i); if (atomComp > 0.0) { double maxPermissible = m_elemAbundancesGoal[i] / atomComp; if (m_molNumSpecies_old[kspec] > maxPermissible) { @@ -255,7 +255,7 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) retn = 1; } for (size_t j = 0; j < m_numComponents; ++j) { - aa[j + i*m_numElemConstraints] = - m_formulaMatrix[j][i]; + aa[j + i*m_numElemConstraints] = - m_formulaMatrix(i,j); } } int info; @@ -338,7 +338,7 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) double saveDir = 0.0; bool goodSpec = true; for (size_t i = 0; i < m_numComponents; ++i) { - double dir = m_formulaMatrix[i][kspec] * (m_elemAbundancesGoal[i] - m_elemAbundances[i]); + double dir = m_formulaMatrix(kspec,i) * (m_elemAbundancesGoal[i] - m_elemAbundances[i]); if (fabs(dir) > 1.0E-10) { if (dir > 0.0) { if (saveDir < 0.0) { @@ -353,7 +353,7 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) } saveDir = dir; } else { - if (m_formulaMatrix[i][kspec] != 0.) { + if (m_formulaMatrix(kspec,i) != 0.) { goodSpec = false; break; } @@ -363,8 +363,8 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) int its = 0; double xx = 0.0; for (size_t i = 0; i < m_numComponents; ++i) { - if (m_formulaMatrix[i][kspec] != 0.0) { - xx += (m_elemAbundancesGoal[i] - m_elemAbundances[i]) / m_formulaMatrix[i][kspec]; + if (m_formulaMatrix(kspec,i) != 0.0) { + xx += (m_elemAbundancesGoal[i] - m_elemAbundances[i]) / m_formulaMatrix(kspec,i); its++; } } @@ -397,8 +397,8 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) (m_elType[i] == VCS_ELEM_TYPE_ABSPOS && m_elemAbundancesGoal[i] == 0.0)) { for (size_t kspec = 0; kspec < m_numSpeciesRdc; kspec++) { if (m_elemAbundances[i] > 0.0) { - if (m_formulaMatrix[i][kspec] < 0.0) { - m_molNumSpecies_old[kspec] -= m_elemAbundances[i] / m_formulaMatrix[i][kspec] ; + if (m_formulaMatrix(kspec,i) < 0.0) { + m_molNumSpecies_old[kspec] -= m_elemAbundances[i] / m_formulaMatrix(kspec,i) ; if (m_molNumSpecies_old[kspec] < 0.0) { m_molNumSpecies_old[kspec] = 0.0; } @@ -407,8 +407,8 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) } } if (m_elemAbundances[i] < 0.0) { - if (m_formulaMatrix[i][kspec] > 0.0) { - m_molNumSpecies_old[kspec] -= m_elemAbundances[i] / m_formulaMatrix[i][kspec]; + if (m_formulaMatrix(kspec,i) > 0.0) { + m_molNumSpecies_old[kspec] -= m_elemAbundances[i] / m_formulaMatrix(kspec,i); if (m_molNumSpecies_old[kspec] < 0.0) { m_molNumSpecies_old[kspec] = 0.0; } @@ -435,13 +435,13 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) bool useZeroed = true; for (size_t kspec = 0; kspec < m_numSpeciesRdc; kspec++) { if (dev < 0.0) { - if (m_formulaMatrix[i][kspec] < 0.0) { + if (m_formulaMatrix(kspec,i) < 0.0) { if (m_molNumSpecies_old[kspec] > 0.0) { useZeroed = false; } } } else { - if (m_formulaMatrix[i][kspec] > 0.0) { + if (m_formulaMatrix(kspec,i) > 0.0) { if (m_molNumSpecies_old[kspec] > 0.0) { useZeroed = false; } @@ -451,8 +451,8 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) for (size_t kspec = 0; kspec < m_numSpeciesRdc; kspec++) { if (m_molNumSpecies_old[kspec] > 0.0 || useZeroed) { if (dev < 0.0) { - if (m_formulaMatrix[i][kspec] < 0.0) { - double delta = dev / m_formulaMatrix[i][kspec] ; + if (m_formulaMatrix(kspec,i) < 0.0) { + double delta = dev / m_formulaMatrix(kspec,i) ; m_molNumSpecies_old[kspec] += delta; if (m_molNumSpecies_old[kspec] < 0.0) { m_molNumSpecies_old[kspec] = 0.0; @@ -462,8 +462,8 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) } } if (dev > 0.0) { - if (m_formulaMatrix[i][kspec] > 0.0) { - double delta = dev / m_formulaMatrix[i][kspec] ; + if (m_formulaMatrix(kspec,i) > 0.0) { + double delta = dev / m_formulaMatrix(kspec,i) ; m_molNumSpecies_old[kspec] += delta; if (m_molNumSpecies_old[kspec] < 0.0) { m_molNumSpecies_old[kspec] = 0.0; diff --git a/src/equil/vcs_elem_rearrange.cpp b/src/equil/vcs_elem_rearrange.cpp index 88b31c30b..3ad9d84c7 100644 --- a/src/equil/vcs_elem_rearrange.cpp +++ b/src/equil/vcs_elem_rearrange.cpp @@ -102,7 +102,7 @@ int VCS_SOLVE::vcs_elem_rearrange(double* const aw, double* const sa, * from the current component. */ for (size_t j = 0; j < ncomponents; ++j) { - sm[j + jr*ncomponents] = m_formulaMatrix[k][j]; + sm[j + jr*ncomponents] = m_formulaMatrix(j,k); } if (jl > 0) { /* @@ -204,7 +204,7 @@ void VCS_SOLVE::vcs_switch_elem_pos(size_t ipos, size_t jpos) std::swap(m_elType[ipos], m_elType[jpos]); std::swap(m_elementActive[ipos], m_elementActive[jpos]); for (size_t j = 0; j < m_numSpeciesTot; ++j) { - std::swap(m_formulaMatrix[ipos][j], m_formulaMatrix[jpos][j]); + std::swap(m_formulaMatrix(j,ipos), m_formulaMatrix(j,jpos)); } std::swap(m_elementName[ipos], m_elementName[jpos]); } diff --git a/src/equil/vcs_inest.cpp b/src/equil/vcs_inest.cpp index 67df2f3ab..32dae0db9 100644 --- a/src/equil/vcs_inest.cpp +++ b/src/equil/vcs_inest.cpp @@ -51,7 +51,7 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm, if (m_elementActive[j]) { double tmp = 0.0; for (size_t kspec = 0; kspec < nspecies; ++kspec) { - tmp += m_formulaMatrix[j][kspec] * m_molNumSpecies_old[kspec]; + tmp += m_formulaMatrix(kspec,j) * m_molNumSpecies_old[kspec]; } plogf("%s ", pprefix); plogf(" %-9.9s", (m_elementName[j]).c_str()); @@ -184,11 +184,11 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm, * exp(-m_deltaGRxn_new[irxn]); for (size_t k = 0; k < m_numComponents; ++k) { - m_deltaMolNumSpecies[k] += m_stoichCoeffRxnMatrix[irxn][k] * m_deltaMolNumSpecies[kspec]; + m_deltaMolNumSpecies[k] += m_stoichCoeffRxnMatrix(k,irxn) * m_deltaMolNumSpecies[kspec]; } for (iph = 0; iph < m_numPhases; iph++) { - m_deltaPhaseMoles[iph] += m_deltaMolNumPhase[irxn][iph] * m_deltaMolNumSpecies[kspec]; + m_deltaPhaseMoles[iph] += m_deltaMolNumPhase(iph,irxn) * m_deltaMolNumSpecies[kspec]; } } } diff --git a/src/equil/vcs_phaseStability.cpp b/src/equil/vcs_phaseStability.cpp index ffee1388b..f9d0a48b8 100644 --- a/src/equil/vcs_phaseStability.cpp +++ b/src/equil/vcs_phaseStability.cpp @@ -66,7 +66,7 @@ bool VCS_SOLVE::vcs_popPhasePossible(const size_t iphasePop) const */ for (size_t j = 0; j < m_numComponents; ++j) { if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) { - double stoicC = m_stoichCoeffRxnMatrix[irxn][j]; + double stoicC = m_stoichCoeffRxnMatrix(j,irxn); if (stoicC != 0.0) { double negChangeComp = - stoicC; if (negChangeComp > 0.0) { @@ -91,11 +91,11 @@ bool VCS_SOLVE::vcs_popPhasePossible(const size_t iphasePop) const for (size_t jrxn = 0; jrxn < m_numRxnRdc; jrxn++) { bool foundJrxn = false; // First, if the component is a product of the reaction - if (m_stoichCoeffRxnMatrix[jrxn][kspec] > 0.0) { + if (m_stoichCoeffRxnMatrix(kspec,jrxn) > 0.0) { foundJrxn = true; // We can do the reaction if all other reactant components have positive mole fractions for (size_t kcomp = 0; kcomp < m_numComponents; kcomp++) { - if (m_stoichCoeffRxnMatrix[jrxn][kcomp] < 0.0) { + if (m_stoichCoeffRxnMatrix(kcomp,jrxn) < 0.0) { if (m_molNumSpecies_old[kcomp] <= VCS_DELETE_ELEMENTABS_CUTOFF*0.5) { foundJrxn = false; } @@ -106,7 +106,7 @@ bool VCS_SOLVE::vcs_popPhasePossible(const size_t iphasePop) const } } // Second we are here if the component is a reactant in the reaction, and the reaction goes backwards. - else if (m_stoichCoeffRxnMatrix[jrxn][kspec] < 0.0) { + else if (m_stoichCoeffRxnMatrix(kspec,jrxn) < 0.0) { foundJrxn = true; size_t jspec = jrxn + m_numComponents; if (m_molNumSpecies_old[jspec] <= VCS_DELETE_ELEMENTABS_CUTOFF*0.5) { @@ -115,7 +115,7 @@ bool VCS_SOLVE::vcs_popPhasePossible(const size_t iphasePop) const } // We can do the backwards reaction if all of the product components species are positive for (size_t kcomp = 0; kcomp < m_numComponents; kcomp++) { - if (m_stoichCoeffRxnMatrix[jrxn][kcomp] > 0.0) { + if (m_stoichCoeffRxnMatrix(kcomp,jrxn) > 0.0) { if (m_molNumSpecies_old[kcomp] <= VCS_DELETE_ELEMENTABS_CUTOFF*0.5) { foundJrxn = false; } @@ -162,7 +162,7 @@ int VCS_SOLVE::vcs_phasePopDeterminePossibleList() vcs_VolPhase* Vphase = m_VolPhaseList[iph]; int existence = Vphase->exists(); if (existence < 0) { - if (m_stoichCoeffRxnMatrix[irxn][j] > 0.0) { + if (m_stoichCoeffRxnMatrix(j,irxn) > 0.0) { if (std::find(jList.begin(), jList.end(), iph) != jList.end()) { jList.push_back(iph); } @@ -198,13 +198,13 @@ int VCS_SOLVE::vcs_phasePopDeterminePossibleList() for (size_t j = 0; j < m_numComponents; j++) { if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) { if (m_molNumSpecies_old[j] <= 0.0) { - if (m_stoichCoeffRxnMatrix[irxn][j] < 0.0) { + if (m_stoichCoeffRxnMatrix(j,irxn) < 0.0) { bool foundPos = false; for (size_t kk = 0; kk < nsp; kk++) { size_t kkspec = Vphase->spGlobalIndexVCS(kk); if (kkspec >= m_numComponents) { size_t iirxn = kkspec - m_numComponents; - if (m_stoichCoeffRxnMatrix[iirxn][j] > 0.0) { + if (m_stoichCoeffRxnMatrix(j,iirxn) > 0.0) { foundPos = true; } } @@ -393,11 +393,10 @@ int VCS_SOLVE::vcs_popPhaseRxnStepSizes(const size_t iphasePop) // Section for a single-species phase if (Vphase->m_singleSpecies) { double s = 0.0; - double* dnPhase_irxn = m_deltaMolNumPhase[irxn]; for (size_t 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]; + s += SQUARE(m_stoichCoeffRxnMatrix(j,irxn)) / m_molNumSpecies_old[j]; } } } @@ -405,7 +404,7 @@ int VCS_SOLVE::vcs_popPhaseRxnStepSizes(const size_t iphasePop) Vphase = m_VolPhaseList[j]; if (! Vphase->m_singleSpecies) { if (m_tPhaseMoles_old[j] > 0.0) { - s -= SQUARE(dnPhase_irxn[j]) / m_tPhaseMoles_old[j]; + s -= SQUARE(m_deltaMolNumPhase(j,irxn)) / m_tPhaseMoles_old[j]; } } } @@ -423,7 +422,7 @@ int VCS_SOLVE::vcs_popPhaseRxnStepSizes(const size_t iphasePop) * section to do damping of the m_deltaMolNumSpecies[] */ for (size_t j = 0; j < m_numComponents; ++j) { - double stoicC = m_stoichCoeffRxnMatrix[irxn][j]; + double stoicC = m_stoichCoeffRxnMatrix(j,irxn); if (stoicC != 0.0) { if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) { double negChangeComp = - stoicC * m_deltaMolNumSpecies[kspec]; @@ -468,7 +467,7 @@ int VCS_SOLVE::vcs_popPhaseRxnStepSizes(const size_t iphasePop) if (kspec >= m_numComponents) { irxn = kspec - m_numComponents; for (size_t j = 0; j < m_numComponents; ++j) { - double stoicC = m_stoichCoeffRxnMatrix[irxn][j]; + double stoicC = m_stoichCoeffRxnMatrix(j,irxn); if (stoicC != 0.0) { if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) { molNumSpecies_tmp[j] += stoicC * delmol; @@ -628,7 +627,7 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph) size_t kspec = Vphase->spGlobalIndexVCS(k); if (kspec >= m_numComponents) { size_t irxn = kspec - m_numComponents; - fracDelta_old[kc] += m_stoichCoeffRxnMatrix[irxn][kc_spec] * fracDelta_old[k]; + fracDelta_old[kc] += m_stoichCoeffRxnMatrix(kc_spec,irxn) * fracDelta_old[k]; } } } @@ -687,10 +686,9 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph) if (i == 0) { m_deltaGRxn_Deficient[irxn] = m_deltaGRxn_old[irxn]; } - double* dtmp_ptr = m_stoichCoeffRxnMatrix[irxn]; - if (dtmp_ptr[kc_spec] != 0.0) { + if (m_stoichCoeffRxnMatrix(kc_spec,irxn) != 0.0) { m_deltaGRxn_Deficient[irxn] += - dtmp_ptr[kc_spec] * (m_feSpecies_Deficient[kc_spec]- m_feSpecies_old[kc_spec]); + m_stoichCoeffRxnMatrix(kc_spec,irxn) * (m_feSpecies_Deficient[kc_spec]- m_feSpecies_old[kc_spec]); } } @@ -743,7 +741,7 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph) size_t kspec = Vphase->spGlobalIndexVCS(k); if (kspec >= m_numComponents) { size_t irxn = kspec - m_numComponents; - fracDelta_raw[kc] += m_stoichCoeffRxnMatrix[irxn][kc_spec] * fracDelta_raw[k]; + fracDelta_raw[kc] += m_stoichCoeffRxnMatrix(kc_spec,irxn) * fracDelta_raw[k]; } } } diff --git a/src/equil/vcs_prep.cpp b/src/equil/vcs_prep.cpp index dd6493c17..6e81adc2a 100644 --- a/src/equil/vcs_prep.cpp +++ b/src/equil/vcs_prep.cpp @@ -213,7 +213,7 @@ int VCS_SOLVE::vcs_prep() vcs_vdzero(m_feSpecies_old, m_numSpeciesTot); vcs_vdzero(m_feSpecies_new, m_numSpeciesTot); vcs_vdzero(m_molNumSpecies_new, m_numSpeciesTot); - vcs_dzero(&(m_deltaMolNumPhase[0][0]), m_numSpeciesTot * m_numPhases); + vcs_dzero(&(m_deltaMolNumPhase(0,0)), m_numSpeciesTot * m_numPhases); vcs_izero(&(m_phaseParticipation[0][0]), m_numSpeciesTot * m_numPhases); vcs_dzero(VCS_DATA_PTR(m_deltaPhaseMoles), m_numPhases); vcs_dzero(VCS_DATA_PTR(m_tPhaseMoles_new), m_numPhases); diff --git a/src/equil/vcs_prob.cpp b/src/equil/vcs_prob.cpp index 30194c229..d362195ca 100644 --- a/src/equil/vcs_prob.cpp +++ b/src/equil/vcs_prob.cpp @@ -73,7 +73,7 @@ VCS_PROB::VCS_PROB(size_t nsp, size_t nel, size_t nph) : w.resize(nspecies, 0.0); mf.resize(nspecies, 0.0); gai.resize(ne, 0.0); - FormulaMatrix.resize(ne, nspecies, 0.0); + FormulaMatrix.resize(nspecies, ne, 0.0); SpeciesUnknownType.resize(nspecies, VCS_SPECIES_TYPE_MOLNUM); VolPM.resize(nspecies, 0.0); PhaseID.resize(nspecies, npos); @@ -123,7 +123,7 @@ void VCS_PROB::resizeSpecies(size_t nsp, int force) m_gibbsSpecies.resize(nsp, 0.0); w.resize(nsp, 0.0); mf.resize(nsp, 0.0); - FormulaMatrix.resize(NE0, nsp, 0.0); + FormulaMatrix.resize(nsp, NE0, 0.0); SpeciesUnknownType.resize(nsp, VCS_SPECIES_TYPE_MOLNUM); VolPM.resize(nsp, 0.0); PhaseID.resize(nsp, 0); @@ -143,7 +143,7 @@ void VCS_PROB::resizeElements(size_t nel, int force) { if (force || nel > NE0) { gai.resize(nel, 0.0); - FormulaMatrix.resize(nel, NSPECIES0, 0.0); + FormulaMatrix.resize(NSPECIES0, nel, 0.0); ElName.resize(nel, ""); m_elType.resize(nel, VCS_ELEM_TYPE_ABSPOS); ElActive.resize(nel, 1); @@ -157,13 +157,12 @@ void VCS_PROB::resizeElements(size_t nel, int force) void VCS_PROB::set_gai() { double* ElemAbund = VCS_DATA_PTR(gai); - double* const* const fm = FormulaMatrix.baseDataAddr(); vcs_dzero(ElemAbund, ne); for (size_t j = 0; j < ne; j++) { for (size_t kspec = 0; kspec < nspecies; kspec++) { if (SpeciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - ElemAbund[j] += fm[j][kspec] * w[kspec]; + ElemAbund[j] += FormulaMatrix(kspec,j) * w[kspec]; } } } @@ -345,13 +344,13 @@ size_t VCS_PROB::addOnePhaseSpecies(vcs_VolPhase* volPhase, size_t k, size_t kT) plogf("Shouldn't be here\n"); exit(EXIT_FAILURE); } - double const* const* const fm = volPhase->getFormulaMatrix(); + const Cantera::Array2D& fm = volPhase->getFormulaMatrix(); for (size_t eVP = 0; eVP < volPhase->nElemConstraints(); eVP++) { size_t e = volPhase->elemGlobalIndex(eVP); if (DEBUG_MODE_ENABLED && e == npos) { exit(EXIT_FAILURE); } - FormulaMatrix[e][kT] = fm[eVP][k]; + FormulaMatrix(kT,e) = fm(k,eVP); } /* * Tell the phase object about the current position of the diff --git a/src/equil/vcs_report.cpp b/src/equil/vcs_report.cpp index 0a81b68e6..4931f630b 100644 --- a/src/equil/vcs_report.cpp +++ b/src/equil/vcs_report.cpp @@ -185,7 +185,7 @@ int VCS_SOLVE::vcs_report(int iconv) plogf("%-10.10s", m_speciesName[kspec].c_str()); plogf("|%10.3g |", m_molNumSpecies_old[kspec]*molScale); for (size_t j = 0; j < m_numComponents; j++) { - plogf(" %6.2f", m_stoichCoeffRxnMatrix[irxn][j]); + plogf(" %6.2f", m_stoichCoeffRxnMatrix(j,irxn)); } plogf(" |%10.3g |", m_deltaGRxn_new[irxn]); plogf("\n"); diff --git a/src/equil/vcs_rxnadj.cpp b/src/equil/vcs_rxnadj.cpp index 82095ff83..27271b628 100644 --- a/src/equil/vcs_rxnadj.cpp +++ b/src/equil/vcs_rxnadj.cpp @@ -68,9 +68,6 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) sprintf(ANOTE, "ZeroedPhase: Phase is artificially zeroed"); } } else if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - - double* dnPhase_irxn = m_deltaMolNumPhase[irxn]; - if (m_molNumSpecies_old[kspec] == 0.0 && (!m_SSPhase[kspec])) { /********************************************************************/ /******* MULTISPECIES PHASE WITH total moles equal to zero *********/ @@ -179,7 +176,7 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) for (size_t 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]; + s += SQUARE(m_stoichCoeffRxnMatrix(j,irxn)) / m_molNumSpecies_old[j]; } } } @@ -187,7 +184,7 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) vcs_VolPhase* Vphase = m_VolPhaseList[j]; if (!Vphase->m_singleSpecies) { if (m_tPhaseMoles_old[j] > 0.0) { - s -= SQUARE(dnPhase_irxn[j]) / m_tPhaseMoles_old[j]; + s -= SQUARE(m_deltaMolNumPhase(j,irxn)) / m_tPhaseMoles_old[j]; } } } @@ -209,7 +206,7 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) m_deltaMolNumSpecies[kspec] = -m_deltaGRxn_new[irxn] / s; // New section to do damping of the m_deltaMolNumSpecies[] for (size_t j = 0; j < m_numComponents; ++j) { - double stoicC = m_stoichCoeffRxnMatrix[irxn][j]; + double stoicC = m_stoichCoeffRxnMatrix(j,irxn); if (stoicC != 0.0) { double negChangeComp = -stoicC * m_deltaMolNumSpecies[kspec]; if (negChangeComp > m_molNumSpecies_old[j]) { @@ -259,8 +256,8 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) dss = m_molNumSpecies_old[kspec]; k = kspec; for (size_t j = 0; j < m_numComponents; ++j) { - if (m_stoichCoeffRxnMatrix[irxn][j] > 0.0) { - double xx = m_molNumSpecies_old[j] / m_stoichCoeffRxnMatrix[irxn][j]; + if (m_stoichCoeffRxnMatrix(j,irxn) > 0.0) { + double xx = m_molNumSpecies_old[j] / m_stoichCoeffRxnMatrix(j,irxn); if (xx < dss) { dss = xx; k = j; @@ -271,8 +268,8 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) } else { dss = 1.0e10; for (size_t j = 0; j < m_numComponents; ++j) { - if (m_stoichCoeffRxnMatrix[irxn][j] < 0.0) { - double xx = -m_molNumSpecies_old[j] / m_stoichCoeffRxnMatrix[irxn][j]; + if (m_stoichCoeffRxnMatrix(j,irxn) < 0.0) { + double xx = -m_molNumSpecies_old[j] / m_stoichCoeffRxnMatrix(j,irxn); if (xx < dss) { dss = xx; k = j; @@ -318,7 +315,7 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) } m_deltaMolNumSpecies[kspec] = dss; for (size_t j = 0; j < m_numComponents; ++j) { - m_deltaMolNumSpecies[j] = dss * m_stoichCoeffRxnMatrix[irxn][j]; + m_deltaMolNumSpecies[j] = dss * m_stoichCoeffRxnMatrix(j,irxn); } iphDel = m_phaseID[k]; @@ -398,8 +395,6 @@ int VCS_SOLVE::vcs_rxn_adj_cg() } size_t kspec = m_indexRxnToSpecies[irxn]; - double* dnPhase_irxn = m_deltaMolNumPhase[irxn]; - if (m_molNumSpecies_old[kspec] == 0.0 && (!m_SSPhase[kspec])) { /* *******************************************************************/ /* **** MULTISPECIES PHASE WITH total moles equal to zero ************/ @@ -467,13 +462,13 @@ int VCS_SOLVE::vcs_rxn_adj_cg() } for (size_t j = 0; j < m_numComponents; ++j) { if (!m_SSPhase[j]) { - s += SQUARE(m_stoichCoeffRxnMatrix[irxn][j]) / m_molNumSpecies_old[j]; + s += SQUARE(m_stoichCoeffRxnMatrix(j,irxn)) / m_molNumSpecies_old[j]; } } for (size_t j = 0; j < m_numPhases; j++) { if (!(m_VolPhaseList[j])->m_singleSpecies) { if (m_tPhaseMoles_old[j] > 0.0) { - s -= SQUARE(dnPhase_irxn[j]) / m_tPhaseMoles_old[j]; + s -= SQUARE(m_deltaMolNumPhase(j,irxn)) / m_tPhaseMoles_old[j]; } } } @@ -497,8 +492,8 @@ int VCS_SOLVE::vcs_rxn_adj_cg() dss = m_molNumSpecies_old[kspec]; k = kspec; for (size_t j = 0; j < m_numComponents; ++j) { - if (m_stoichCoeffRxnMatrix[irxn][j] > 0.0) { - double xx = m_molNumSpecies_old[j] / m_stoichCoeffRxnMatrix[irxn][j]; + if (m_stoichCoeffRxnMatrix(j,irxn) > 0.0) { + double xx = m_molNumSpecies_old[j] / m_stoichCoeffRxnMatrix(j,irxn); if (xx < dss) { dss = xx; k = j; @@ -509,8 +504,8 @@ int VCS_SOLVE::vcs_rxn_adj_cg() } else { dss = 1.0e10; for (size_t j = 0; j < m_numComponents; ++j) { - if (m_stoichCoeffRxnMatrix[irxn][j] < 0.0) { - double xx = -m_molNumSpecies_old[j] / m_stoichCoeffRxnMatrix[irxn][j]; + if (m_stoichCoeffRxnMatrix(j,irxn) < 0.0) { + double xx = -m_molNumSpecies_old[j] / m_stoichCoeffRxnMatrix(j,irxn); if (xx < dss) { dss = xx; k = j; @@ -530,8 +525,8 @@ int VCS_SOLVE::vcs_rxn_adj_cg() m_molNumSpecies_old[kspec] += dss; m_tPhaseMoles_old[m_phaseID[kspec]] += dss; for (size_t 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_molNumSpecies_old[j] += dss * m_stoichCoeffRxnMatrix(j,irxn); + m_tPhaseMoles_old[m_phaseID[j]] += dss * m_stoichCoeffRxnMatrix(j,irxn); } m_molNumSpecies_old[k] = 0.0; m_tPhaseMoles_old[m_phaseID[k]] = 0.0; @@ -607,11 +602,11 @@ double VCS_SOLVE::vcs_Hessian_actCoeff_diag(size_t irxn) if (np_kspec < 1.0E-13) { np_kspec = 1.0E-13; } - double* sc_irxn = m_stoichCoeffRxnMatrix[irxn]; + double* sc_irxn = m_stoichCoeffRxnMatrix.ptrColumn(irxn); /* * First the diagonal term of the Jacobian */ - double s = m_np_dLnActCoeffdMolNum[kspec][kspec] / np_kspec; + double 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. @@ -622,12 +617,12 @@ double VCS_SOLVE::vcs_Hessian_actCoeff_diag(size_t irxn) if (m_phaseID[k] == m_phaseID[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; + s += sc_irxn[k] * sc_irxn[l] * m_np_dLnActCoeffdMolNum(l,k) / np; } } } if (kph == m_phaseID[l]) { - s += sc_irxn[l] * (m_np_dLnActCoeffdMolNum[kspec][l] + m_np_dLnActCoeffdMolNum[l][kspec]) / np_kspec; + s += sc_irxn[l] * (m_np_dLnActCoeffdMolNum(l,kspec) + m_np_dLnActCoeffdMolNum(kspec,l)) / np_kspec; } } @@ -655,7 +650,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_np_dLnActCoeffdMolNum.baseDataAddr()); + Vphase->sendToVCS_LnActCoeffJac(m_np_dLnActCoeffdMolNum); } } } @@ -671,9 +666,8 @@ double VCS_SOLVE::deltaG_Recalc_Rxn(const int stateCalc, const size_t irxn, cons } } double deltaG = mu_i[kspec]; - double* sc_irxn = m_stoichCoeffRxnMatrix[irxn]; for (size_t k = 0; k < m_numComponents; k++) { - deltaG += sc_irxn[k] * mu_i[k]; + deltaG += m_stoichCoeffRxnMatrix(k,irxn) * mu_i[k]; } return deltaG; } @@ -684,7 +678,7 @@ double VCS_SOLVE::vcs_line_search(const size_t irxn, const double dx_orig, char* size_t kspec = m_indexRxnToSpecies[irxn]; const int MAXITS = 10; double dx = dx_orig; - double* sc_irxn = m_stoichCoeffRxnMatrix[irxn]; + double* sc_irxn = m_stoichCoeffRxnMatrix.ptrColumn(irxn); double* molNumBase = VCS_DATA_PTR(m_molNumSpecies_old); double* acBase = VCS_DATA_PTR(m_actCoeffSpecies_old); double* ac = VCS_DATA_PTR(m_actCoeffSpecies_new); diff --git a/src/equil/vcs_setMolesLinProg.cpp b/src/equil/vcs_setMolesLinProg.cpp index dd95b5afb..161a0bfd7 100644 --- a/src/equil/vcs_setMolesLinProg.cpp +++ b/src/equil/vcs_setMolesLinProg.cpp @@ -116,7 +116,7 @@ int VCS_SOLVE::vcs_setMolesLinProg() ik = m_numComponents + irxn; dg_rt = m_SSfeSpecies[ik]; dxi_min = 1.0e10; - const double* sc_irxn = m_stoichCoeffRxnMatrix[irxn]; + const double* sc_irxn = m_stoichCoeffRxnMatrix.ptrColumn(irxn); for (size_t jcomp = 0; jcomp < m_numElemConstraints; jcomp++) { dg_rt += m_SSfeSpecies[jcomp] * sc_irxn[jcomp]; } diff --git a/src/equil/vcs_solve.cpp b/src/equil/vcs_solve.cpp index d2b0c5078..cdeadfdf3 100644 --- a/src/equil/vcs_solve.cpp +++ b/src/equil/vcs_solve.cpp @@ -97,7 +97,7 @@ void VCS_SOLVE::vcs_initSizes(const size_t nspecies0, const size_t nelements, * We will initialize sc[] to note the fact that it needs to be * filled with meaningful information. */ - m_stoichCoeffRxnMatrix.resize(nspecies0, nelements, 0.0); + m_stoichCoeffRxnMatrix.resize(nelements, nspecies0, 0.0); m_scSize.resize(nspecies0, 0.0); m_spSize.resize(nspecies0, 1.0); @@ -108,7 +108,7 @@ void VCS_SOLVE::vcs_initSizes(const size_t nspecies0, const size_t nelements, m_speciesUnknownType.resize(nspecies0, VCS_SPECIES_TYPE_MOLNUM); - m_deltaMolNumPhase.resize(nspecies0, nphase0, 0.0); + m_deltaMolNumPhase.resize(nphase0, nspecies0, 0.0); m_phaseParticipation.resize(nspecies0, nphase0, 0); m_phasePhi.resize(nphase0, 0.0); @@ -130,7 +130,7 @@ void VCS_SOLVE::vcs_initSizes(const size_t nspecies0, const size_t nelements, m_TmpPhase.resize(nphase0, 0.0); m_TmpPhase2.resize(nphase0, 0.0); - m_formulaMatrix.resize(nelements, nspecies0); + m_formulaMatrix.resize(nspecies0, nelements); TPhInertMoles.resize(nphase0, 0.0); @@ -445,10 +445,10 @@ int VCS_SOLVE::vcs_prob_specifyFully(const VCS_PROB* pub) for (size_t i = 0; i < nspecies; i++) { bool nonzero = false; for (size_t j = 0; j < nelements; j++) { - if (pub->FormulaMatrix[j][i] != 0.0) { + if (pub->FormulaMatrix(i,j) != 0.0) { nonzero = true; } - m_formulaMatrix[j][i] = pub->FormulaMatrix[j][i]; + m_formulaMatrix(i,j) = pub->FormulaMatrix(i,j); } if (!nonzero) { plogf("vcs_prob_specifyFully:: species %d %s has a zero formula matrix!\n", i, @@ -522,7 +522,7 @@ int VCS_SOLVE::vcs_prob_specifyFully(const VCS_PROB* pub) for (size_t kspec = 0; kspec < nspecies; kspec++) { if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { sum += m_molNumSpecies_old[kspec]; - m_elemAbundancesGoal[j] += m_formulaMatrix[j][kspec] * m_molNumSpecies_old[kspec]; + m_elemAbundancesGoal[j] += m_formulaMatrix(kspec,j) * m_molNumSpecies_old[kspec]; } } if (pub->m_elType[j] == VCS_ELEM_TYPE_LATTICERATIO) { diff --git a/src/equil/vcs_solve_TP.cpp b/src/equil/vcs_solve_TP.cpp index b9d66c2ff..a1281177d 100644 --- a/src/equil/vcs_solve_TP.cpp +++ b/src/equil/vcs_solve_TP.cpp @@ -168,7 +168,7 @@ int VCS_SOLVE::vcs_solve_TP(int print_lvl, int printDetails, int maxit) for (size_t i = 0; i < m_numSpeciesTot; ++i) { plogf(" %-18.18s", m_speciesName[i].c_str()); for (size_t j = 0; j < m_numElemConstraints; ++j) { - plogf("% -7.3g ", m_formulaMatrix[j][i]); + plogf("% -7.3g ", m_formulaMatrix(i,j)); } plogf(" %3d ", m_phaseID[i]); writeline(' ', std::max(55-int(m_numElemConstraints)*8, 0), false); @@ -603,7 +603,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, } for (size_t irxn = 0; irxn < m_numRxnRdc; irxn++) { size_t kspec = m_indexRxnToSpecies[irxn]; - double* sc_irxn = m_stoichCoeffRxnMatrix[irxn]; + double* sc_irxn = m_stoichCoeffRxnMatrix.ptrColumn(irxn); size_t iph = m_phaseID[kspec]; vcs_VolPhase* Vphase = m_VolPhaseList[iph]; if (DEBUG_MODE_ENABLED) { @@ -662,7 +662,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, for (size_t j = 0; j < m_numElemConstraints; ++j) { int elType = m_elType[j]; if (elType == VCS_ELEM_TYPE_ABSPOS) { - double atomComp = m_formulaMatrix[j][kspec]; + double atomComp = m_formulaMatrix(kspec,j); if (atomComp > 0.0) { double maxPermissible = m_elemAbundancesGoal[j] / atomComp; if (maxPermissible < VCS_DELETE_MINORSPECIES_CUTOFF) { @@ -998,9 +998,8 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, * Calculate the tentative change in the total number of * moles in all of the phases */ - double* dnPhase_irxn = m_deltaMolNumPhase[irxn]; for (iph = 0; iph < m_numPhases; iph++) { - m_deltaPhaseMoles[iph] += dx * dnPhase_irxn[iph]; + m_deltaPhaseMoles[iph] += dx * m_deltaMolNumPhase(iph,irxn); } } @@ -1397,27 +1396,27 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, (m_molNumSpecies_old[j] * m_spSize[j] * 1.01); } } - if (doSwap && m_stoichCoeffRxnMatrix[i][j] != 0.0) { + if (doSwap && m_stoichCoeffRxnMatrix(j,i) != 0.0) { if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- Get a new basis because "); plogf("%s", m_speciesName[l].c_str()); plogf(" is better than comp "); plogf("%s", m_speciesName[j].c_str()); plogf(" and share nonzero stoic: %-9.1f", - m_stoichCoeffRxnMatrix[i][j]); + m_stoichCoeffRxnMatrix(j,i)); plogendl(); } forceComponentCalc = 1; return; } #ifdef DEBUG_NOT - if (m_speciesStatus[l] == VCS_SPECIES_ZEROEDMS && m_molNumSpecies_old[j] == 0.0 && m_stoichCoeffRxnMatrix[i][j] != 0.0 && dg[i] < 0.0) { + if (m_speciesStatus[l] == VCS_SPECIES_ZEROEDMS && m_molNumSpecies_old[j] == 0.0 && m_stoichCoeffRxnMatrix(j,i) != 0.0 && dg[i] < 0.0) { if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- Get a new basis because %s", m_speciesName[l].c_str()); plogf(" has dg < 0.0 and comp %s has zero mole num", m_speciesName[j].c_str()); plogf(" and share nonzero stoic: %-9.1f", - m_stoichCoeffRxnMatrix[i][j]); + m_stoichCoeffRxnMatrix(j,i)); plogendl(); } return; @@ -1754,7 +1753,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_np_dLnActCoeffdMolNum[kspec][kspec] / (m_tPhaseMoles_old[iph]); + s = m_np_dLnActCoeffdMolNum(kspec,kspec) / (m_tPhaseMoles_old[iph]); /* * We fit it to a power law approximation of the activity coefficient * @@ -1846,7 +1845,7 @@ int VCS_SOLVE::delta_species(const size_t kspec, double* const delta_ptr) * one would work, */ double dx = delta; - double* sc_irxn = m_stoichCoeffRxnMatrix[irxn]; + double* sc_irxn = m_stoichCoeffRxnMatrix.ptrColumn(irxn); for (size_t j = 0; j < m_numComponents; ++j) { if (m_molNumSpecies_old[j] > 0.0) { double tmp = sc_irxn[j] * dx; @@ -2116,8 +2115,8 @@ bool VCS_SOLVE::vcs_delete_multiphase(const size_t iph) for (size_t kspec = m_numComponents; kspec < m_numSpeciesRdc; ++kspec) { size_t irxn = kspec - m_numComponents; if (m_phaseID[kspec] != iph) { - if (m_stoichCoeffRxnMatrix[irxn][kcomp] != 0.0) { - dxWant = -m_molNumSpecies_old[kcomp] / m_stoichCoeffRxnMatrix[irxn][kcomp]; + if (m_stoichCoeffRxnMatrix(kcomp,irxn) != 0.0) { + dxWant = -m_molNumSpecies_old[kcomp] / m_stoichCoeffRxnMatrix(kcomp,irxn); if (dxWant + m_molNumSpecies_old[kspec] < 0.0) { dxPerm = -m_molNumSpecies_old[kspec]; } @@ -2126,9 +2125,9 @@ bool VCS_SOLVE::vcs_delete_multiphase(const size_t iph) if (m_phaseID[jcomp] == iph) { dxPerm = 0.0; } else { - dj = dxWant * m_stoichCoeffRxnMatrix[irxn][jcomp]; + dj = dxWant * m_stoichCoeffRxnMatrix(jcomp,irxn); if (dj + m_molNumSpecies_old[kcomp] < 0.0) { - dxPerm2 = -m_molNumSpecies_old[kcomp] / m_stoichCoeffRxnMatrix[irxn][jcomp]; + dxPerm2 = -m_molNumSpecies_old[kcomp] / m_stoichCoeffRxnMatrix(jcomp,irxn); } if (fabs(dxPerm2) < fabs(dxPerm)) { dxPerm = dxPerm2; @@ -2556,7 +2555,7 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[], bool lindep; size_t juse = npos; size_t jlose = npos; - double* dptr, *scrxn_ptr; + double* scrxn_ptr; Cantera::clockWC tickTock; if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" "); @@ -2589,7 +2588,7 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[], vcs_print_stringTrunc(m_speciesName[k].c_str(), 11, 1); plogf(" | "); for (size_t j = 0; j < m_numElemConstraints; j++) { - plogf(" %8.2g", m_formulaMatrix[j][k]); + plogf(" %8.2g", m_formulaMatrix(k,j)); } plogf("\n"); } @@ -2694,7 +2693,7 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[], for (size_t j = 0; j < m_numElemConstraints; ++j) { if (m_elementActive[j]) { if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) { - double nu = m_formulaMatrix[j][kspec]; + double nu = m_formulaMatrix(kspec,j); if (nu != 0.0) { nonZeroesKspec++; maxConcPossKspec = std::min(m_elemAbundancesGoal[j] / nu, maxConcPossKspec); @@ -2774,7 +2773,7 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[], */ size_t jl = jr; for (size_t j = 0; j < m_numElemConstraints; ++j) { - sm[j + jr*m_numElemConstraints] = m_formulaMatrix[j][k]; + sm[j + jr*m_numElemConstraints] = m_formulaMatrix(k,j); } if (jl > 0) { /* @@ -2893,24 +2892,24 @@ L_END_LOOP: */ for (size_t j = 0; j < ncTrial; ++j) { for (size_t i = 0; i < ncTrial; ++i) { - sm[i + j*m_numElemConstraints] = m_formulaMatrix[i][j]; + sm[i + j*m_numElemConstraints] = m_formulaMatrix(j,i); } } for (size_t i = 0; i < m_numRxnTot; ++i) { k = m_indexRxnToSpecies[i]; for (size_t j = 0; j < ncTrial; ++j) { - m_stoichCoeffRxnMatrix[i][j] = - m_formulaMatrix[j][k]; + m_stoichCoeffRxnMatrix(j,i) = - m_formulaMatrix(k,j); } } // Solve the linear system to calculate the reaction matrix, - // m_stoichCoeffRxnMatrix[][]. + // m_stoichCoeffRxnMatrix. ct_dgetrf(ncTrial, ncTrial, sm, m_numElemConstraints, &ipiv[0], info); if (info) { plogf("vcs_solve_TP ERROR: Error factorizing stoichiometric coefficient matrix\n"); return VCS_FAILED_CONVERGENCE; } ct_dgetrs(ctlapack::NoTranspose, ncTrial, m_numRxnTot, sm, m_numElemConstraints, - &ipiv[0], m_stoichCoeffRxnMatrix[0], m_numElemConstraints, info); + &ipiv[0], m_stoichCoeffRxnMatrix.ptrColumn(0), m_numElemConstraints, info); /* * NOW, if we have interfacial voltage unknowns, what we did @@ -2939,9 +2938,9 @@ L_END_LOOP: for (size_t j = 0; j < ncTrial; ++j) { for (size_t i = 0; i < ncTrial; ++i) { if (i == jlose) { - sm[i + j*m_numElemConstraints] = m_formulaMatrix[juse][j]; + sm[i + j*m_numElemConstraints] = m_formulaMatrix(j,juse); } else { - sm[i + j*m_numElemConstraints] = m_formulaMatrix[i][j]; + sm[i + j*m_numElemConstraints] = m_formulaMatrix(j,i); } } } @@ -2949,9 +2948,9 @@ L_END_LOOP: k = m_indexRxnToSpecies[i]; for (size_t j = 0; j < ncTrial; ++j) { if (j == jlose) { - aw[j] = - m_formulaMatrix[juse][k]; + aw[j] = - m_formulaMatrix(k,juse); } else { - aw[j] = - m_formulaMatrix[j][k]; + aw[j] = - m_formulaMatrix(k,j); } } } @@ -2965,7 +2964,7 @@ L_END_LOOP: &ipiv[0], aw, m_numElemConstraints, info); size_t i = k - ncTrial; for (size_t j = 0; j < ncTrial; j++) { - m_stoichCoeffRxnMatrix[i][j] = aw[j]; + m_stoichCoeffRxnMatrix(j,i) = aw[j]; } } } @@ -2977,7 +2976,7 @@ L_END_LOOP: for (size_t i = 0; i < m_numRxnTot; i++) { double szTmp = 0.0; for (size_t j = 0; j < ncTrial; j++) { - szTmp += fabs(m_stoichCoeffRxnMatrix[i][j]); + szTmp += fabs(m_stoichCoeffRxnMatrix(j,i)); } m_scSize[i] = szTmp; } @@ -3009,7 +3008,7 @@ L_END_LOOP: plogf("|% -10.3E|", m_molNumSpecies_old[m_indexRxnToSpecies[i]]); } for (size_t j = 0; j < ncTrial; j++) { - plogf(" %+7.3f", m_stoichCoeffRxnMatrix[i][j]); + plogf(" %+7.3f", m_stoichCoeffRxnMatrix(j,i)); } plogf("\n"); } @@ -3027,17 +3026,17 @@ L_END_LOOP: double sum; for (size_t j = 0; j < ncTrial; ++j) { if (j == jlose) { - sum = m_formulaMatrix[juse][k]; + sum = m_formulaMatrix(k,juse); for (size_t n = 0; n < ncTrial; n++) { - double numElements = m_formulaMatrix[juse][n]; - double coeff = m_stoichCoeffRxnMatrix[i][n]; + double numElements = m_formulaMatrix(n,juse); + double coeff = m_stoichCoeffRxnMatrix(n,i); sum += coeff * numElements; } } else { - sum = m_formulaMatrix[j][k]; + sum = m_formulaMatrix(k,j); for (size_t n = 0; n < ncTrial; n++) { - double numElements = m_formulaMatrix[j][n]; - double coeff = m_stoichCoeffRxnMatrix[i][n]; + double numElements = m_formulaMatrix(n,j); + double coeff = m_stoichCoeffRxnMatrix(n,i); sum += coeff * numElements; } } @@ -3077,27 +3076,26 @@ L_END_LOOP: /* * Zero out the change of Phase Moles array */ - vcs_dzero(m_deltaMolNumPhase[0], (NSPECIES0)*(NPHASE0)); + vcs_dzero(&m_deltaMolNumPhase(0,0), (NSPECIES0)*(NPHASE0)); vcs_izero(m_phaseParticipation[0], (NSPECIES0)*(NPHASE0)); /* * Loop over each reaction, creating the change in Phase Moles - * array, m_deltaMolNumPhase[irxn][iphase], + * array, m_deltaMolNumPhase(iphase,irxn), * and the phase participation array, PhaseParticipation[irxn][iphase] */ for (size_t irxn = 0; irxn < m_numRxnTot; ++irxn) { - scrxn_ptr = m_stoichCoeffRxnMatrix[irxn]; - dptr = m_deltaMolNumPhase[irxn]; + scrxn_ptr = m_stoichCoeffRxnMatrix.ptrColumn(irxn); size_t kspec = m_indexRxnToSpecies[irxn]; size_t iph = m_phaseID[kspec]; int* pp_ptr = m_phaseParticipation[irxn]; - dptr[iph] = 1.0; + m_deltaMolNumPhase(iph,irxn) = 1.0; pp_ptr[iph]++; for (size_t j = 0; j < ncTrial; ++j) { iph = m_phaseID[j]; if (fabs(scrxn_ptr[j]) <= 1.0e-6) { scrxn_ptr[j] = 0.0; } else { - dptr[iph] += scrxn_ptr[j]; + m_deltaMolNumPhase(iph,irxn) += scrxn_ptr[j]; pp_ptr[iph]++; } } @@ -3199,7 +3197,7 @@ int VCS_SOLVE::vcs_species_type(const size_t kspec) const for (size_t j = 0; j < m_numElemConstraints; ++j) { int elType = m_elType[j]; if (elType == VCS_ELEM_TYPE_ABSPOS) { - double atomComp = m_formulaMatrix[j][kspec]; + double atomComp = m_formulaMatrix(kspec,j); if (atomComp > 0.0) { double maxPermissible = m_elemAbundancesGoal[j] / atomComp; if (maxPermissible < VCS_DELETE_MINORSPECIES_CUTOFF) { @@ -3229,7 +3227,7 @@ int VCS_SOLVE::vcs_species_type(const size_t kspec) const */ if (irxn >= 0) { for (size_t j = 0; j < m_numComponents; ++j) { - double stoicC = m_stoichCoeffRxnMatrix[irxn][j]; + double stoicC = m_stoichCoeffRxnMatrix(j,irxn); if (stoicC != 0.0) { double negChangeComp = - stoicC; if (negChangeComp > 0.0) { @@ -3348,7 +3346,7 @@ int VCS_SOLVE::vcs_species_type(const size_t kspec) const double szAdj = m_scSize[irxn] * std::sqrt((double)m_numRxnTot); for (size_t k = 0; k < m_numComponents; ++k) { if (!(m_SSPhase[k])) { - if (m_stoichCoeffRxnMatrix[irxn][k] != 0.0) { + if (m_stoichCoeffRxnMatrix(k,irxn) != 0.0) { if (m_molNumSpecies_old[kspec] * szAdj >= m_molNumSpecies_old[k] * 0.01) { return VCS_SPECIES_MAJOR; } @@ -3758,7 +3756,7 @@ void VCS_SOLVE::prneav() const for (size_t j = 0; j < m_numElemConstraints; ++j) { for (size_t i = 0; i < m_numSpeciesTot; ++i) { if (m_speciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - eav[j] += m_formulaMatrix[j][i] * m_molNumSpecies_old[i]; + eav[j] += m_formulaMatrix(i,j) * m_molNumSpecies_old[i]; } } } @@ -3950,20 +3948,6 @@ bool VCS_SOLVE::vcs_evaluate_speciesType() return (m_numRxnMinorZeroed >= m_numRxnRdc); } -void VCS_SOLVE::vcs_switch2D(double* const* const Jac, - const size_t k1, const size_t k2) const -{ - if (k1 == k2) { - return; - } - for (size_t i = 0; i < m_numSpeciesTot; i++) { - std::swap(Jac[k1][i], Jac[k2][i]); - } - for (size_t i = 0; i < m_numSpeciesTot; i++) { - std::swap(Jac[i][k1], Jac[i][k2]); - } -} - void VCS_SOLVE::vcs_deltag(const int l, const bool doDeleted, const int vcsState, const bool alterZeroedPhases) { @@ -4011,7 +3995,7 @@ void VCS_SOLVE::vcs_deltag(const int l, const bool doDeleted, if (m_speciesStatus[kspec] != VCS_SPECIES_MINOR) { icase = 0; deltaGRxn[irxn] = feSpecies[m_indexRxnToSpecies[irxn]]; - double* dtmp_ptr = m_stoichCoeffRxnMatrix[irxn]; + double* dtmp_ptr = m_stoichCoeffRxnMatrix.ptrColumn(irxn); for (kspec = 0; kspec < m_numComponents; ++kspec) { deltaGRxn[irxn] += dtmp_ptr[kspec] * feSpecies[kspec]; if (molNumSpecies[kspec] < VCS_DELETE_MINORSPECIES_CUTOFF && dtmp_ptr[kspec] < 0.0) { @@ -4030,7 +4014,7 @@ void VCS_SOLVE::vcs_deltag(const int l, const bool doDeleted, for (size_t irxn = 0; irxn < irxnl; ++irxn) { icase = 0; deltaGRxn[irxn] = feSpecies[m_indexRxnToSpecies[irxn]]; - double* dtmp_ptr = m_stoichCoeffRxnMatrix[irxn]; + double* dtmp_ptr = m_stoichCoeffRxnMatrix.ptrColumn(irxn); for (size_t kspec = 0; kspec < m_numComponents; ++kspec) { deltaGRxn[irxn] += dtmp_ptr[kspec] * feSpecies[kspec]; if (molNumSpecies[kspec] < VCS_DELETE_MINORSPECIES_CUTOFF && @@ -4051,7 +4035,7 @@ void VCS_SOLVE::vcs_deltag(const int l, const bool doDeleted, if (m_speciesStatus[kspec] <= VCS_SPECIES_MINOR) { icase = 0; deltaGRxn[irxn] = feSpecies[m_indexRxnToSpecies[irxn]]; - double* dtmp_ptr = m_stoichCoeffRxnMatrix[irxn]; + double* dtmp_ptr = m_stoichCoeffRxnMatrix.ptrColumn(irxn); for (kspec = 0; kspec < m_numComponents; ++kspec) { deltaGRxn[irxn] += dtmp_ptr[kspec] * feSpecies[kspec]; if (m_molNumSpecies_old[kspec] < VCS_DELETE_MINORSPECIES_CUTOFF && @@ -4207,7 +4191,7 @@ void VCS_SOLVE::vcs_printDeltaG(const int stateCalc) plogf("|%10.3g|", m_molNumSpecies_old[m_indexRxnToSpecies[i]]); } for (size_t j = 0; j < m_numComponents; j++) { - plogf(" %6.2f", m_stoichCoeffRxnMatrix[i][j]); + plogf(" %6.2f", m_stoichCoeffRxnMatrix(j,i)); } plogf("\n"); } @@ -4342,9 +4326,8 @@ void VCS_SOLVE::vcs_deltag_Phase(const size_t iphase, const bool doDeleted, if (kspec >= m_numComponents) { size_t irxn = kspec - m_numComponents; deltaGRxn[irxn] = feSpecies[kspec]; - double* dtmp_ptr = m_stoichCoeffRxnMatrix[irxn]; for (size_t kcomp = 0; kcomp < m_numComponents; ++kcomp) { - deltaGRxn[irxn] += dtmp_ptr[kcomp] * feSpecies[kcomp]; + deltaGRxn[irxn] += m_stoichCoeffRxnMatrix(kcomp,irxn) * feSpecies[kcomp]; } } } @@ -4362,9 +4345,8 @@ void VCS_SOLVE::vcs_deltag_Phase(const size_t iphase, const bool doDeleted, zeroedPhase = false; } deltaGRxn[irxn] = feSpecies[kspec]; - double* dtmp_ptr = m_stoichCoeffRxnMatrix[irxn]; for (size_t kcomp = 0; kcomp < m_numComponents; ++kcomp) { - deltaGRxn[irxn] += dtmp_ptr[kcomp] * feSpecies[kcomp]; + deltaGRxn[irxn] += m_stoichCoeffRxnMatrix(kcomp,irxn) * feSpecies[kcomp]; } } } @@ -4491,10 +4473,15 @@ void VCS_SOLVE::vcs_switch_pos(const bool ifunc, const size_t k1, const size_t k std::swap(m_PMVolumeSpecies[k1], m_PMVolumeSpecies[k2]); for (size_t j = 0; j < m_numElemConstraints; ++j) { - std::swap(m_formulaMatrix[j][k1], m_formulaMatrix[j][k2]); + std::swap(m_formulaMatrix(k1,j), m_formulaMatrix(k2,j)); } - if (m_useActCoeffJac) { - vcs_switch2D(m_np_dLnActCoeffdMolNum.baseDataAddr(), k1, k2); + if (m_useActCoeffJac && k1 != k2) { + for (size_t i = 0; i < m_numSpeciesTot; i++) { + std::swap(m_np_dLnActCoeffdMolNum(k1,i), m_np_dLnActCoeffdMolNum(k2,i)); + } + for (size_t i = 0; i < m_numSpeciesTot; i++) { + std::swap(m_np_dLnActCoeffdMolNum(i,k1), m_np_dLnActCoeffdMolNum(i,k2)); + } } std::swap(m_speciesStatus[k1], m_speciesStatus[k2]); /* @@ -4514,11 +4501,11 @@ void VCS_SOLVE::vcs_switch_pos(const bool ifunc, const size_t k1, const size_t k i1 , i2); } for (size_t j = 0; j < m_numComponents; ++j) { - std::swap(m_stoichCoeffRxnMatrix[i1][j], m_stoichCoeffRxnMatrix[i2][j]); + std::swap(m_stoichCoeffRxnMatrix(j,i1), m_stoichCoeffRxnMatrix(j,i2)); } std::swap(m_scSize[i1], m_scSize[i2]); for (size_t iph = 0; iph < m_numPhases; iph++) { - std::swap(m_deltaMolNumPhase[i1][iph], m_deltaMolNumPhase[i2][iph]); + std::swap(m_deltaMolNumPhase(iph,i1), m_deltaMolNumPhase(iph,i2)); std::swap(m_phaseParticipation[i1][iph], m_phaseParticipation[i2][iph]); } @@ -4574,7 +4561,7 @@ double VCS_SOLVE::vcs_birthGuess(const int kspec) * back from this routine. This evaluation should * be respected. */ - double* sc_irxn = m_stoichCoeffRxnMatrix[irxn]; + double* sc_irxn = m_stoichCoeffRxnMatrix.ptrColumn(irxn); for (size_t j = 0; j < m_numComponents; ++j) { // Only loop over element constraints that involve positive def. constraints if (m_speciesUnknownType[j] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {