diff --git a/include/cantera/equil/vcs_MultiPhaseEquil.h b/include/cantera/equil/vcs_MultiPhaseEquil.h index 5f147f05f..658646c47 100644 --- a/include/cantera/equil/vcs_MultiPhaseEquil.h +++ b/include/cantera/equil/vcs_MultiPhaseEquil.h @@ -15,26 +15,6 @@ namespace Cantera { -//! Translate a MultiPhase object into a VCS_SOLVE problem definition object -/*! - * @param mphase MultiPhase object that is the source for all of the information - * @param vsolve VCS_SOLVE problem definition that gets all of the information - * - * Note, both objects share the underlying ThermoPhase objects. So, neither can - * be const objects. - */ -int vcs_Cantera_to_vprob(MultiPhase* mphase, VCS_SOLVE* vsolve); - -//! Translate a MultiPhase information into a VCS_SOLVE problem definition object -/*! - * This version updates the problem statement information only. All species and - * phase definitions remain the same. - * - * @param mphase MultiPhase object that is the source for all of the information - * @param vsolve VCS_SOLVE problem definition that gets all of the information - */ -int vcs_Cantera_update_vprob(MultiPhase* mphase, VCS_SOLVE* vsolve); - //! %Cantera's Interface to the Multiphase chemical equilibrium solver. /*! * Class vcs_MultiPhaseEquil is designed to be used to set a mixture containing @@ -243,10 +223,6 @@ public: */ void reportCSV(const std::string& reportFile); - // Friend functions - friend int vcs_Cantera_to_vprob(MultiPhase* mphase, VCS_SOLVE* vsolve); - friend int vcs_Cantera_update_vprob(MultiPhase* mphase, VCS_SOLVE* vsolve); - protected: //! Vector that takes into account of the current sorting of the species /*! diff --git a/include/cantera/equil/vcs_defs.h b/include/cantera/equil/vcs_defs.h index 5d84223f0..5ea310c4a 100644 --- a/include/cantera/equil/vcs_defs.h +++ b/include/cantera/equil/vcs_defs.h @@ -25,15 +25,6 @@ namespace Cantera #define VCS_MP_FAIL -6 //@} -/*! - * @name Type of the underlying equilibrium solve - * @{ - */ - -//! Current, it is always done holding T and P constant. -#define VCS_PROBTYPE_TP 0 -//@} - /*! * @name Sizes of Phases and Cutoff Mole Numbers * diff --git a/include/cantera/equil/vcs_solve.h b/include/cantera/equil/vcs_solve.h index 0f1f13ac0..3ef15ea00 100644 --- a/include/cantera/equil/vcs_solve.h +++ b/include/cantera/equil/vcs_solve.h @@ -46,20 +46,9 @@ class VCS_SOLVE public: //! Initialize the sizes within the VCS_SOLVE object /*! - * This resizes all of the internal arrays within the object. This routine - * operates in two modes. If all of the parameters are the same as it - * currently exists in the object, nothing is done by this routine; a quick - * exit is carried out and all of the data in the object persists. - * - * If any of the parameters are different than currently exists in the - * object, then all of the data in the object must be redone. It may not be - * zeroed, but it must be redone. - * - * @param nspecies0 Number of species within the object - * @param nelements Number of element constraints within the problem - * @param nphase0 Number of phases defined within the problem. + * This resizes all of the internal arrays within the object. */ - VCS_SOLVE(size_t nspecies0, size_t nelements, size_t nphase0); + VCS_SOLVE(MultiPhase* mphase, int printLvl=0); ~VCS_SOLVE(); @@ -67,16 +56,6 @@ public: /*! * This is the main interface routine to the equilibrium solver * - * @param ifunc Determines the operation to be done: Valid values: - * 0 -> Solve a new problem by initializing structures first. An - * initial estimate may or may not have been already - * determined. This is indicated in the VCS_SOLVE structure. - * 1 -> The problem has already been initialized and set up. We - * call this routine to resolve it using the problem - * statement and solution estimate contained in the VCS_SOLVE - * structure. - * 2 -> Don't solve a problem. Destroy all the private - * structures. * @param ipr Printing of results * ipr = 1 -> Print problem statement and final results to * standard output @@ -87,7 +66,7 @@ public: * @return nonzero value: failure to solve the problem at hand. zero : * success */ - int vcs(int ifunc, int ipr, int ip1, int maxit); + int vcs(int ipr, int ip1, int maxit); //! Main routine that solves for equilibrium at constant T and P using a //! variant of the VCS method @@ -530,10 +509,6 @@ public: //! consistent with what's needed for solution. It is called one time for //! each new problem structure definition. /*! - * This routine is always followed by vcs_prep(). Therefore, tasks that need - * to be done for every call to vcsc() should be placed in vcs_prep() and - * not in this routine. - * * The problem structure refers to: * * - the number and identity of the species. @@ -552,42 +527,18 @@ public: * done so. During this process the order of the species is changed * in the private data structure. All references to the species * properties must employ the ind[] index vector. + * 4. Initialization of arrays to zero. + * 5. Check to see if the problem is well posed (If all the element + * abundances are zero, the algorithm will fail) * * @param printLvl Print level of the routine - * @return VCS_SUCCESS = everything went OK - */ - int vcs_prep_oneTime(int printLvl); - - //! Prepare the object for solution - /*! - * This routine is mostly concerned with changing the private data to be - * consistent with that needed for solution. It is called for every - * invocation of the vcs_solve() except for the cleanup invocation. - * - * Tasks: - * 1. Initialization of arrays to zero. - * * @return * VCS_SUCCESS = everything went OK; * VCS_PUB_BAD = There is an irreconcilable difference in the * public data structure from when the problem was * initially set up. */ - int vcs_prep(); - - //! In this routine, we check for things that will cause the algorithm - //! to fail. - /*! - * We check to see if the problem is well posed. If it is not, we return - * false and print out error conditions. - * - * Current there is one condition. If all the element abundances are zero, - * the algorithm will fail. - * - * @return If true, the problem is well-posed. If false, the problem - * is not well posed. - */ - bool vcs_wellPosed(); + int vcs_prep(int printLvl); //! Rearrange the constraint equations represented by the Formula //! Matrix so that the operational ones are in the front @@ -817,12 +768,6 @@ public: //! Fully specify the problem to be solved int vcs_prob_specifyFully(); - //! Specify the problem to be solved, incrementally - /*! - * It's assumed we are solving the same problem. - */ - int vcs_prob_specify(); - private: //! Zero out the concentration of a species. /*! @@ -1063,20 +1008,6 @@ private: public: //! @{ Variables moved from VCS_PROB - //! Problem type. I.e., the identity of what is held constant. Currently, T - //! and P are held constant, and this input is ignored - int prob_type; - - //! Total number of species in the problems - size_t nspecies; - - //! Number of element constraints in the equilibrium problem - size_t ne; - //! Number of element constraints used to size data structures - //! involving elements - size_t NE0; - //! Number of phases in the problem - size_t NPhase; //! Vector of chemical potentials of the species. This is a calculated //! output quantity. length = number of species. vector_fp m_gibbsSpecies; @@ -1119,42 +1050,6 @@ public: */ vector_int SpeciesUnknownType; - //! Temperature (Kelvin) - /*! - * Specification of the temperature for the equilibrium problem - */ - double T; - - //! Pressure - double PresPA; - - //! Volume of the entire system - /*! - * Note, this is an output variable atm - */ - double Vol; - - //! Partial Molar Volumes of species - /*! - * This is a calculated vector, calculated from w[]. - * length number of species. - */ - vector_fp VolPM; - - //! Specification of the initial estimate method - /*! - * * 0: user estimate - * * 1: user estimate if satisifies elements - * * -1: machine estimate - */ - int iest; - - //! Tolerance requirement for major species - double tolmaj; - - //! Tolerance requirement for minor species - double tolmin; - //! Mapping between the species and the phases std::vector PhaseID; @@ -1177,20 +1072,11 @@ public: //! Array of phase structures std::vector VPhaseList; - // String containing the title of the run - std::string Title; - //! Vector of pointers to thermo structures which identify the model and //! parameters for evaluating the thermodynamic functions for that //! particular species std::vector SpeciesThermo; - //! Number of iterations. This is an output variable - int m_Iterations; - - //! Number of basis optimizations used. This is an output variable. - int m_NumBasisOptimizations; - //! Print level for print routines int m_printLvl; @@ -1258,20 +1144,11 @@ public: void reportCSV(const std::string& reportFile); - //! value of the number of species used to size data structures - size_t NSPECIES0; - - //! value of the number of phases used to size data structures - size_t NPHASE0; - //! Total number of species in the problems - size_t m_numSpeciesTot; + size_t m_nsp; //! Number of element constraints in the problem - /*! - * This is typically equal to the number of elements in the problem - */ - size_t m_numElemConstraints; + size_t m_nelem; //! Number of components calculated for the problem size_t m_numComponents; @@ -1620,9 +1497,6 @@ public: //! Array of Phase Structures. Length = number of phases. std::vector m_VolPhaseList; - //! String containing the title of the run - std::string m_title; - //! This specifies the current state of units for the Gibbs free energy //! properties in the program. /*! diff --git a/src/equil/vcs_MultiPhaseEquil.cpp b/src/equil/vcs_MultiPhaseEquil.cpp index dbadb8887..b8466d1d7 100644 --- a/src/equil/vcs_MultiPhaseEquil.cpp +++ b/src/equil/vcs_MultiPhaseEquil.cpp @@ -22,20 +22,10 @@ using namespace std; namespace Cantera { vcs_MultiPhaseEquil::vcs_MultiPhaseEquil(MultiPhase* mix, int printLvl) : - m_mix(0), + m_mix(mix), m_printLvl(printLvl), - m_vsolve(mix->nSpecies(), mix->nElements(), mix->nPhases()) + m_vsolve(mix, printLvl) { - m_mix = mix; - m_vsolve.m_printLvl = m_printLvl; - m_vsolve.m_mix = m_mix; - - // Work out the details of the VCS_SOLVE construction and Transfer the - // current problem to the VCS_SOLVE object - int res = vcs_Cantera_to_vprob(mix, &m_vsolve); - if (res != 0) { - plogf("problems\n"); - } } int vcs_MultiPhaseEquil::equilibrate_TV(int XY, doublereal xtarget, @@ -430,20 +420,7 @@ int vcs_MultiPhaseEquil::equilibrate_TP(int estimateEquil, clockWC tickTock; m_printLvl = printLvl; m_vsolve.m_printLvl = printLvl; - - // Extract the current state information from the MultiPhase object and - // Transfer it to VCS_PROB object. - int res = vcs_Cantera_update_vprob(m_mix, &m_vsolve); - if (res != 0) { - plogf("problems\n"); - } - - // Set the estimation technique - if (estimateEquil) { - m_vsolve.iest = estimateEquil; - } else { - m_vsolve.iest = 0; - } + m_vsolve.m_doEstimateEquil = estimateEquil; // Check obvious bounds on the temperature and pressure NOTE, we may want to // do more here with the real bounds given by the ThermoPhase objects. @@ -456,10 +433,6 @@ int vcs_MultiPhaseEquil::equilibrate_TP(int estimateEquil, "Pressure less than zero on input"); } - // Print out the problem specification from the point of - // view of the vprob object. - m_vsolve.prob_report(m_printLvl); - //! Call the thermo Program int ip1 = m_printLvl; int ipr = std::max(0, m_printLvl-1); @@ -468,7 +441,7 @@ int vcs_MultiPhaseEquil::equilibrate_TP(int estimateEquil, } else { ip1 = 0; } - int iSuccess = m_vsolve.vcs(0, ipr, ip1, maxit); + int iSuccess = m_vsolve.vcs(ipr, ip1, maxit); // Transfer the information back to the MultiPhase object. Note we don't // just call setMoles, because some multispecies solution phases may be @@ -477,7 +450,7 @@ int vcs_MultiPhaseEquil::equilibrate_TP(int estimateEquil, // about likely reemergent states. m_mix->uploadMoleFractionsFromPhases(); size_t kGlob = 0; - for (size_t ip = 0; ip < m_vsolve.NPhase; ip++) { + for (size_t ip = 0; ip < m_vsolve.m_numPhases; ip++) { double phaseMole = 0.0; ThermoPhase& tref = m_mix->phase(ip); for (size_t k = 0; k < tref.nSpecies(); k++, kGlob++) { @@ -493,8 +466,8 @@ int vcs_MultiPhaseEquil::equilibrate_TP(int estimateEquil, plogf("\nVCS FAILED TO CONVERGE!\n"); } plogf("\n"); - plogf("Temperature = %g Kelvin\n", m_vsolve.T); - plogf("Pressure = %g Pa\n", m_vsolve.PresPA); + plogf("Temperature = %g Kelvin\n", m_vsolve.m_temperature); + plogf("Pressure = %g Pa\n", m_vsolve.m_pressurePA); plogf("\n"); plogf("----------------------------------------" "---------------------\n"); @@ -502,7 +475,7 @@ int vcs_MultiPhaseEquil::equilibrate_TP(int estimateEquil, plogf(" Mole_Fraction Chem_Potential (J/kmol)\n"); plogf("--------------------------------------------------" "-----------\n"); - for (size_t i = 0; i < m_vsolve.nspecies; i++) { + for (size_t i = 0; i < m_mix->nSpecies(); i++) { plogf("%-12s", m_mix->speciesName(i)); if (m_vsolve.SpeciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { plogf(" %15.3e %15.3e ", 0.0, m_vsolve.mf[i]); @@ -533,7 +506,7 @@ int vcs_MultiPhaseEquil::equilibrate_TP(int estimateEquil, void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile) { - size_t nphase = m_vsolve.NPhase; + size_t nphase = m_vsolve.m_numPhases; FILE* FP = fopen(reportFile.c_str(), "w"); if (!FP) { @@ -572,8 +545,8 @@ void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile) fprintf(FP,"Temperature = %11.5g kelvin\n", m_mix->temperature()); fprintf(FP,"Pressure = %11.5g Pascal\n", m_mix->pressure()); fprintf(FP,"Total Volume = %11.5g m**3\n", vol); - fprintf(FP,"Number Basis optimizations = %d\n", m_vsolve.m_NumBasisOptimizations); - fprintf(FP,"Number VCS iterations = %d\n", m_vsolve.m_Iterations); + fprintf(FP,"Number Basis optimizations = %d\n", m_vsolve.m_VCount->Basis_Opts); + fprintf(FP,"Number VCS iterations = %d\n", m_vsolve.m_VCount->Its); for (size_t iphase = 0; iphase < nphase; iphase++) { size_t istart = m_mix->speciesIndex(0, iphase); @@ -666,399 +639,4 @@ void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile) fclose(FP); } -// HKM -> Work on transferring the current value of the voltages into the -// equilibrium problem. -int vcs_Cantera_to_vprob(MultiPhase* mphase, VCS_SOLVE* vsolve) -{ - VCS_SPECIES_THERMO* ts_ptr = 0; - - // Calculate the total number of species and phases in the problem - size_t totNumPhases = mphase->nPhases(); - size_t totNumSpecies = mphase->nSpecies(); - - // Problem type has yet to be worked out. - vsolve->prob_type = 0; - vsolve->nspecies = totNumSpecies; - vsolve->ne = 0; - vsolve->NPhase = totNumPhases; - // Set the initial estimate to a machine generated estimate for now - // We will work out the details later. - vsolve->iest = -1; - vsolve->T = mphase->temperature(); - vsolve->PresPA = mphase->pressure(); - vsolve->Vol = mphase->volume(); - vsolve->Title = "MultiPhase Object"; - - int printLvl = vsolve->m_printLvl; - - // Loop over the phases, transferring pertinent information - int kT = 0; - for (size_t iphase = 0; iphase < totNumPhases; iphase++) { - // Get the ThermoPhase object - assume volume phase - ThermoPhase* tPhase = &mphase->phase(iphase); - size_t nelem = tPhase->nElements(); - - // Query Cantera for the equation of state type of the current phase. - std::string eos = tPhase->type(); - bool gasPhase = (eos == "IdealGas"); - - // Find out the number of species in the phase - size_t nSpPhase = tPhase->nSpecies(); - // Find out the name of the phase - string phaseName = tPhase->name(); - - // Call the basic vcs_VolPhase creation routine. - // Properties set here: - // ->PhaseNum = phase number in the thermo problem - // ->GasPhase = Boolean indicating whether it is a gas phase - // ->NumSpecies = number of species in the phase - // ->TMolesInert = Inerts in the phase = 0.0 for cantera - // ->PhaseName = Name of the phase - vcs_VolPhase* VolPhase = vsolve->VPhaseList[iphase]; - VolPhase->resize(iphase, nSpPhase, nelem, phaseName.c_str(), 0.0); - VolPhase->m_gasPhase = gasPhase; - - // Tell the vcs_VolPhase pointer about cantera - VolPhase->setPtrThermoPhase(tPhase); - VolPhase->setTotalMoles(0.0); - - // Set the electric potential of the volume phase from the - // ThermoPhase object's value. - VolPhase->setElectricPotential(tPhase->electricPotential()); - - // Query the ThermoPhase object to find out what convention - // it uses for the specification of activity and Standard State. - VolPhase->p_activityConvention = tPhase->activityConvention(); - - // Assign the value of eqn of state. Handle conflicts here. - if (eos == "IdealGas") { - VolPhase->m_eqnState = VCS_EOS_IDEAL_GAS; - } else if (eos == "ConstDensity") { - VolPhase->m_eqnState = VCS_EOS_CONSTANT; - } else if (eos == "StoichSubstance") { - VolPhase->m_eqnState = VCS_EOS_STOICH_SUB; - } else if (eos == "IdealSolidSoln") { - VolPhase->m_eqnState = VCS_EOS_IDEAL_SOLN; - } else if (eos == "Surf" || eos == "Edge") { - throw CanteraError("VCSnonideal", - "Surface/edge phase not handled yet."); - } else { - if (printLvl > 1) { - writelog("Unknown Cantera EOS to VCSnonideal: '{}'\n", eos); - } - VolPhase->m_eqnState = VCS_EOS_UNK_CANTERA; - } - - // Transfer all of the element information from the ThermoPhase object - // to the vcs_VolPhase object. Also decide whether we need a new charge - // neutrality element in the phase to enforce a charge neutrality - // constraint. We also decide whether this is a single species phase - // with the voltage being the independent variable setting the chemical - // potential of the electrons. - VolPhase->transferElementsFM(tPhase); - - // Combine the element information in the vcs_VolPhase - // object into the vprob object. - vsolve->addPhaseElements(VolPhase); - VolPhase->setState_TP(vsolve->T, vsolve->PresPA); - vector_fp muPhase(tPhase->nSpecies(),0.0); - tPhase->getChemPotentials(&muPhase[0]); - double tMoles = 0.0; - - // Loop through each species in the current phase - for (size_t k = 0; k < nSpPhase; k++) { - // Obtain the molecular weight of the species from the - // ThermoPhase object - vsolve->WtSpecies[kT] = tPhase->molecularWeight(k); - - // Obtain the charges of the species from the ThermoPhase object - vsolve->Charge[kT] = tPhase->charge(k); - - // Set the phaseid of the species - vsolve->PhaseID[kT] = iphase; - - // Transfer the type of unknown - vsolve->SpeciesUnknownType[kT] = VolPhase->speciesUnknownType(k); - if (vsolve->SpeciesUnknownType[kT] == VCS_SPECIES_TYPE_MOLNUM) { - // Set the initial number of kmoles of the species - // and the mole fraction vector - vsolve->w[kT] = mphase->speciesMoles(kT); - tMoles += vsolve->w[kT]; - vsolve->mf[kT] = mphase->moleFraction(kT); - } else if (vsolve->SpeciesUnknownType[kT] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - vsolve->w[kT] = tPhase->electricPotential(); - vsolve->mf[kT] = mphase->moleFraction(kT); - } else { - throw CanteraError(" vcs_Cantera_to_vsolve() ERROR", - "Unknown species type: {}", vsolve->SpeciesUnknownType[kT]); - } - - // transfer chemical potential vector - vsolve->m_gibbsSpecies[kT] = muPhase[k]; - - // Transfer the species information from the - // volPhase structure to the VPROB structure - // This includes: - // FormulaMatrix[][] - // VolPhase->IndSpecies[] - vsolve->addOnePhaseSpecies(VolPhase, k, kT); - - // Get a pointer to the thermo object - ts_ptr = vsolve->SpeciesThermo[kT]; - - // Fill in the vcs_SpeciesProperty structure - vcs_SpeciesProperties* sProp = VolPhase->speciesProperty(k); - sProp->NumElements = vsolve->ne; - sProp->SpName = mphase->speciesName(kT); - sProp->SpeciesThermo = ts_ptr; - sProp->WtSpecies = tPhase->molecularWeight(k); - sProp->FormulaMatrixCol.resize(vsolve->ne, 0.0); - for (size_t e = 0; e < vsolve->ne; e++) { - sProp->FormulaMatrixCol[e] = vsolve->FormulaMatrix(kT,e); - } - sProp->Charge = tPhase->charge(k); - sProp->SurfaceSpecies = false; - sProp->VolPM = 0.0; - - // Transfer the thermo specification of the species - // vsolve->SpeciesThermo[] - - // Add lookback connectivity into the thermo object first - ts_ptr->IndexPhase = iphase; - ts_ptr->IndexSpeciesPhase = k; - ts_ptr->OwningPhase = VolPhase; - - // get a reference to the Cantera species thermo. - MultiSpeciesThermo& sp = tPhase->speciesThermo(); - - int spType = sp.reportType(k); - if (spType == SIMPLE) { - double c[4]; - double minTemp, maxTemp, refPressure; - sp.reportParams(k, spType, c, minTemp, maxTemp, refPressure); - ts_ptr->SS0_Model = VCS_SS0_CONSTANT; - ts_ptr->SS0_T0 = c[0]; - ts_ptr->SS0_H0 = c[1]; - ts_ptr->SS0_S0 = c[2]; - ts_ptr->SS0_Cp0 = c[3]; - if (gasPhase) { - ts_ptr->SSStar_Model = VCS_SSSTAR_IDEAL_GAS; - ts_ptr->SSStar_Vol_Model = VCS_SSVOL_IDEALGAS; - } else { - ts_ptr->SSStar_Model = VCS_SSSTAR_CONSTANT; - ts_ptr->SSStar_Vol_Model = VCS_SSVOL_CONSTANT; - } - } else { - if (vsolve->m_printLvl > 2) { - plogf("vcs_Cantera_convert: Species Type %d not known \n", - spType); - } - ts_ptr->SS0_Model = VCS_SS0_NOTHANDLED; - ts_ptr->SSStar_Model = VCS_SSSTAR_NOTHANDLED; - } - - // Transfer the Volume Information -> NEEDS WORK - if (gasPhase) { - ts_ptr->SSStar_Vol_Model = VCS_SSVOL_IDEALGAS; - ts_ptr->SSStar_Vol0 = 82.05 * 273.15 / 1.0; - } else { - vector_fp phaseTermCoeff(nSpPhase, 0.0); - int nCoeff; - tPhase->getParameters(nCoeff, &phaseTermCoeff[0]); - ts_ptr->SSStar_Vol_Model = VCS_SSVOL_CONSTANT; - ts_ptr->SSStar_Vol0 = phaseTermCoeff[k]; - } - kT++; - } - - // Now go back through the species in the phase and assign a valid mole - // fraction to all phases, even if the initial estimate of the total - // number of moles is zero. - if (tMoles > 0.0) { - for (size_t k = 0; k < nSpPhase; k++) { - size_t kTa = VolPhase->spGlobalIndexVCS(k); - vsolve->mf[kTa] = vsolve->w[kTa] / tMoles; - } - } else { - // Perhaps, we could do a more sophisticated treatment below. - // But, will start with this. - for (size_t k = 0; k < nSpPhase; k++) { - size_t kTa = VolPhase->spGlobalIndexVCS(k); - vsolve->mf[kTa]= 1.0 / (double) nSpPhase; - } - } - - VolPhase->setMolesFromVCS(VCS_STATECALC_OLD, &vsolve->w[0]); - - // Now, calculate a sample naught Gibbs free energy calculation - // at the specified temperature. - for (size_t k = 0; k < nSpPhase; k++) { - vcs_SpeciesProperties* sProp = VolPhase->speciesProperty(k); - ts_ptr = sProp->SpeciesThermo; - ts_ptr->SS0_feSave = VolPhase->G0_calc_one(k)/ GasConstant; - ts_ptr->SS0_TSave = vsolve->T; - } - } - - // Transfer initial element abundances to the vprob object. - // We have to find the mapping index from one to the other - vsolve->gai.resize(vsolve->ne, 0.0); - vsolve->set_gai(); - - // Printout the species information: PhaseID's and mole nums - if (vsolve->m_printLvl > 1) { - writeline('=', 80, true, true); - writeline('=', 16, false); - plogf(" Cantera_to_vprob: START OF PROBLEM STATEMENT "); - writeline('=', 20); - writeline('=', 80); - plogf(" Phase IDs of species\n"); - plogf(" species phaseID phaseName "); - plogf(" Initial_Estimated_kMols\n"); - for (size_t i = 0; i < vsolve->nspecies; i++) { - size_t iphase = vsolve->PhaseID[i]; - - vcs_VolPhase* VolPhase = vsolve->VPhaseList[iphase]; - plogf("%16s %5d %16s", mphase->speciesName(i).c_str(), iphase, - VolPhase->PhaseName.c_str()); - if (vsolve->SpeciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - plogf(" Volts = %-10.5g\n", vsolve->w[i]); - } else { - plogf(" %-10.5g\n", vsolve->w[i]); - } - } - - // Printout of the Phase structure information - writeline('-', 80, true, true); - plogf(" Information about phases\n"); - plogf(" PhaseName PhaseNum SingSpec GasPhase EqnState NumSpec"); - plogf(" TMolesInert Tmoles(kmol)\n"); - - for (size_t iphase = 0; iphase < vsolve->NPhase; iphase++) { - vcs_VolPhase* VolPhase = vsolve->VPhaseList[iphase]; - plogf("%16s %5d %5d %8d %16s %8d %16e ", VolPhase->PhaseName.c_str(), - VolPhase->VP_ID_, VolPhase->m_singleSpecies, - VolPhase->m_gasPhase, VolPhase->eos_name(), - VolPhase->nSpecies(), VolPhase->totalMolesInert()); - plogf("%16e\n", VolPhase->totalMoles()); - } - - writeline('=', 80, true, true); - writeline('=', 16, false); - plogf(" Cantera_to_vprob: END OF PROBLEM STATEMENT "); - writeline('=', 20); - writeline('=', 80); - plogf("\n"); - } - return VCS_SUCCESS; -} - -int vcs_Cantera_update_vprob(MultiPhase* mphase, VCS_SOLVE* vsolve) -{ - size_t totNumPhases = mphase->nPhases(); - size_t kT = 0; - vector_fp tmpMoles; - // Problem type has yet to be worked out. - vsolve->prob_type = 0; - // Whether we have an estimate or not gets overwritten on - // the call to the equilibrium solver. - vsolve->iest = -1; - vsolve->T = mphase->temperature(); - vsolve->PresPA = mphase->pressure(); - vsolve->Vol = mphase->volume(); - - for (size_t iphase = 0; iphase < totNumPhases; iphase++) { - ThermoPhase* tPhase = &mphase->phase(iphase); - vcs_VolPhase* volPhase = vsolve->VPhaseList[iphase]; - - // Set the electric potential of the volume phase from the - // ThermoPhase object's value. - volPhase->setElectricPotential(tPhase->electricPotential()); - - volPhase->setState_TP(vsolve->T, vsolve->PresPA); - vector_fp muPhase(tPhase->nSpecies(),0.0); - tPhase->getChemPotentials(&muPhase[0]); - - // Loop through each species in the current phase - size_t nSpPhase = tPhase->nSpecies(); - tmpMoles.resize(nSpPhase); - for (size_t k = 0; k < nSpPhase; k++) { - tmpMoles[k] = mphase->speciesMoles(kT); - vsolve->w[kT] = mphase->speciesMoles(kT); - vsolve->mf[kT] = mphase->moleFraction(kT); - - // transfer chemical potential vector - vsolve->m_gibbsSpecies[kT] = muPhase[k]; - - kT++; - } - if (volPhase->phiVarIndex() != npos) { - size_t kphi = volPhase->phiVarIndex(); - size_t kglob = volPhase->spGlobalIndexVCS(kphi); - vsolve->w[kglob] = tPhase->electricPotential(); - } - volPhase->setMolesFromVCS(VCS_STATECALC_OLD, &vsolve->w[0]); - if ((nSpPhase == 1) && (volPhase->phiVarIndex() == 0)) { - volPhase->setExistence(VCS_PHASE_EXIST_ALWAYS); - } else if (volPhase->totalMoles() > 0.0) { - volPhase->setExistence(VCS_PHASE_EXIST_YES); - } else { - volPhase->setExistence(VCS_PHASE_EXIST_NO); - } - } - - // Transfer initial element abundances to the vprob object. Put them in the - // front of the object. There may be more constraints than there are - // elements. But, we know the element abundances are in the front of the - // vector. - vsolve->set_gai(); - - // Printout the species information: PhaseID's and mole nums - if (vsolve->m_printLvl > 1) { - writeline('=', 80, true, true); - writeline('=', 20, false); - plogf(" Cantera_to_vprob: START OF PROBLEM STATEMENT "); - writeline('=', 20); - writeline('=', 80); - plogf("\n"); - plogf(" Phase IDs of species\n"); - plogf(" species phaseID phaseName "); - plogf(" Initial_Estimated_kMols\n"); - for (size_t i = 0; i < vsolve->nspecies; i++) { - size_t iphase = vsolve->PhaseID[i]; - vcs_VolPhase* VolPhase = vsolve->VPhaseList[iphase]; - plogf("%16s %5d %16s", mphase->speciesName(i).c_str(), iphase, - VolPhase->PhaseName.c_str()); - if (vsolve->SpeciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - plogf(" Volts = %-10.5g\n", vsolve->w[i]); - } else { - plogf(" %-10.5g\n", vsolve->w[i]); - } - } - - // Printout of the Phase structure information - writeline('-', 80, true, true); - plogf(" Information about phases\n"); - plogf(" PhaseName PhaseNum SingSpec GasPhase EqnState NumSpec"); - plogf(" TMolesInert Tmoles(kmol)\n"); - - for (size_t iphase = 0; iphase < vsolve->NPhase; iphase++) { - vcs_VolPhase* VolPhase = vsolve->VPhaseList[iphase]; - plogf("%16s %5d %5d %8d %16s %8d %16e ", VolPhase->PhaseName.c_str(), - VolPhase->VP_ID_, VolPhase->m_singleSpecies, - VolPhase->m_gasPhase, VolPhase->eos_name(), - VolPhase->nSpecies(), VolPhase->totalMolesInert()); - plogf("%16e\n", VolPhase->totalMoles()); - } - - writeline('=', 80, true, true); - writeline('=', 20, false); - plogf(" Cantera_to_vprob: END OF PROBLEM STATEMENT "); - writeline('=', 20); - writeline('=', 80); - plogf("\n"); - } - return VCS_SUCCESS; -} - } diff --git a/src/equil/vcs_TP.cpp b/src/equil/vcs_TP.cpp index c25eb1a41..5c5f29fdd 100644 --- a/src/equil/vcs_TP.cpp +++ b/src/equil/vcs_TP.cpp @@ -60,7 +60,7 @@ int VCS_SOLVE::vcs_evalSS_TP(int ipr, int ip1, double Temp, double pres) void VCS_SOLVE::vcs_fePrep_TP() { - for (size_t i = 0; i < m_numSpeciesTot; ++i) { + for (size_t i = 0; i < m_nsp; ++i) { // For single species phases, initialize the chemical potential with the // value of the standard state chemical potential. This value doesn't // change during the calculation diff --git a/src/equil/vcs_elem.cpp b/src/equil/vcs_elem.cpp index 4154641d1..0cd8c10e4 100644 --- a/src/equil/vcs_elem.cpp +++ b/src/equil/vcs_elem.cpp @@ -16,9 +16,9 @@ namespace Cantera { void VCS_SOLVE::vcs_elab() { - for (size_t j = 0; j < m_numElemConstraints; ++j) { + for (size_t j = 0; j < m_nelem; ++j) { m_elemAbundances[j] = 0.0; - for (size_t i = 0; i < m_numSpeciesTot; ++i) { + for (size_t i = 0; i < m_nsp; ++i) { if (m_speciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { m_elemAbundances[j] += m_formulaMatrix(i,j) * m_molNumSpecies_old[i]; } @@ -30,7 +30,7 @@ bool VCS_SOLVE::vcs_elabcheck(int ibound) { size_t top = m_numComponents; if (ibound) { - top = m_numElemConstraints; + top = m_nelem; } for (size_t i = 0; i < top; ++i) { @@ -50,7 +50,7 @@ bool VCS_SOLVE::vcs_elabcheck(int ibound) // addition of terms. We are limited to 13 digits of finite // arithmetic accuracy. bool multisign = false; - for (size_t kspec = 0; kspec < m_numSpeciesTot; kspec++) { + for (size_t kspec = 0; kspec < m_nsp; kspec++) { double eval = m_formulaMatrix(kspec,i); if (eval < 0.0) { multisign = true; @@ -84,9 +84,9 @@ bool VCS_SOLVE::vcs_elabcheck(int ibound) void VCS_SOLVE::vcs_elabPhase(size_t iphase, double* const elemAbundPhase) { - for (size_t j = 0; j < m_numElemConstraints; ++j) { + for (size_t j = 0; j < m_nelem; ++j) { elemAbundPhase[j] = 0.0; - for (size_t i = 0; i < m_numSpeciesTot; ++i) { + for (size_t i = 0; i < m_nsp; ++i) { if (m_speciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE && m_phaseID[i] == iphase) { elemAbundPhase[j] += m_formulaMatrix(i,j) * m_molNumSpecies_old[i]; } @@ -101,29 +101,29 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) vector_fp ga_save(m_elemAbundances); if (m_debug_print_lvl >= 2) { plogf(" --- vcsc_elcorr: Element abundances correction routine"); - if (m_numElemConstraints != m_numComponents) { - plogf(" (m_numComponents != m_numElemConstraints)"); + if (m_nelem != m_numComponents) { + plogf(" (m_numComponents != m_nelem)"); } plogf("\n"); } - for (size_t i = 0; i < m_numElemConstraints; ++i) { + for (size_t i = 0; i < m_nelem; ++i) { x[i] = m_elemAbundances[i] - m_elemAbundancesGoal[i]; } double l2before = 0.0; - for (size_t i = 0; i < m_numElemConstraints; ++i) { + for (size_t i = 0; i < m_nelem; ++i) { l2before += x[i] * x[i]; } - l2before = sqrt(l2before/m_numElemConstraints); + l2before = sqrt(l2before/m_nelem); // Special section to take out single species, single component, // moles. These are species which have non-zero entries in the // formula matrix, and no other species have zero values either. bool changed = false; - for (size_t i = 0; i < m_numElemConstraints; ++i) { + for (size_t i = 0; i < m_nelem; ++i) { int numNonZero = 0; bool multisign = false; - for (size_t kspec = 0; kspec < m_numSpeciesTot; kspec++) { + for (size_t kspec = 0; kspec < m_nsp; kspec++) { if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { double eval = m_formulaMatrix(kspec,i); if (eval < 0.0) { @@ -136,7 +136,7 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) } if (!multisign) { if (numNonZero < 2) { - for (size_t kspec = 0; kspec < m_numSpeciesTot; kspec++) { + for (size_t kspec = 0; kspec < m_nsp; kspec++) { if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { double eval = m_formulaMatrix(kspec,i); if (eval > 0.0) { @@ -159,7 +159,7 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) } if (numCompNonZero == 1) { double diff = m_elemAbundancesGoal[i]; - for (size_t kspec = m_numComponents; kspec < m_numSpeciesTot; kspec++) { + for (size_t kspec = m_numComponents; kspec < m_nsp; kspec++) { if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { double eval = m_formulaMatrix(kspec,i); diff -= eval * m_molNumSpecies_old[kspec]; @@ -183,10 +183,10 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) // Note, also we can do this over ne, the number of elements, not just the // number of components. changed = false; - for (size_t i = 0; i < m_numElemConstraints; ++i) { + for (size_t i = 0; i < m_nelem; ++i) { int elType = m_elType[i]; if (elType == VCS_ELEM_TYPE_ABSPOS) { - for (size_t kspec = 0; kspec < m_numSpeciesTot; kspec++) { + for (size_t kspec = 0; kspec < m_nsp; kspec++) { if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { double atomComp = m_formulaMatrix(kspec,i); if (atomComp > 0.0) { @@ -238,7 +238,7 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) } } - solve(A, x, 1, m_numElemConstraints); + solve(A, x, 1, m_nelem); // Now apply the new direction without creating negative species. double par = 0.5; @@ -288,7 +288,7 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) // Try some ad hoc procedures for fixing the problem if (retn >= 2) { // First find a species whose adjustment is a win-win situation. - for (size_t kspec = 0; kspec < m_numSpeciesTot; kspec++) { + for (size_t kspec = 0; kspec < m_nsp; kspec++) { if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { continue; } @@ -348,7 +348,7 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) goto L_CLEANUP; } - for (size_t i = 0; i < m_numElemConstraints; ++i) { + for (size_t i = 0; i < m_nelem; ++i) { if (m_elType[i] == VCS_ELEM_TYPE_CHARGENEUTRALITY || (m_elType[i] == VCS_ELEM_TYPE_ABSPOS && m_elemAbundancesGoal[i] == 0.0)) { for (size_t kspec = 0; kspec < m_numSpeciesRdc; kspec++) { @@ -374,7 +374,7 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) // For electron charges element types, we try positive deltas in the species // concentrations to match the desired electron charge exactly. - for (size_t i = 0; i < m_numElemConstraints; ++i) { + for (size_t i = 0; i < m_nelem; ++i) { double dev = m_elemAbundancesGoal[i] - m_elemAbundances[i]; if (m_elType[i] == VCS_ELEM_TYPE_ELECTRONCHARGE && (fabs(dev) > 1.0E-300)) { bool useZeroed = true; @@ -418,14 +418,14 @@ L_CLEANUP: ; vcs_tmoles(); double l2after = 0.0; - for (size_t i = 0; i < m_numElemConstraints; ++i) { + for (size_t i = 0; i < m_nelem; ++i) { l2after += pow(m_elemAbundances[i] - m_elemAbundancesGoal[i], 2); } - l2after = sqrt(l2after/m_numElemConstraints); + l2after = sqrt(l2after/m_nelem); if (m_debug_print_lvl >= 2) { plogf(" --- Elem_Abund: Correct Initial " " Final\n"); - for (size_t i = 0; i < m_numElemConstraints; ++i) { + for (size_t i = 0; i < m_nelem; ++i) { plogf(" --- "); plogf("%-2.2s", m_elementName[i]); plogf(" %20.12E %20.12E %20.12E\n", m_elemAbundancesGoal[i], ga_save[i], m_elemAbundances[i]); diff --git a/src/equil/vcs_elem_rearrange.cpp b/src/equil/vcs_elem_rearrange.cpp index dce1dee98..67130609a 100644 --- a/src/equil/vcs_elem_rearrange.cpp +++ b/src/equil/vcs_elem_rearrange.cpp @@ -33,7 +33,7 @@ int VCS_SOLVE::vcs_elem_rearrange(double* const aw, double* const sa, double test = -1.0E10; while (lindep) { lindep = false; - for (size_t i = 0; i < m_numElemConstraints; ++i) { + for (size_t i = 0; i < m_nelem; ++i) { test -= 1.0; aw[i] = m_elemAbundancesGoal[i]; if (test == aw[i]) { @@ -53,14 +53,14 @@ int VCS_SOLVE::vcs_elem_rearrange(double* const aw, double* const sa, while (true) { // Search the remaining part of the mole fraction vector, AW, for // the largest remaining species. Return its identity in K. - k = m_numElemConstraints; - for (size_t ielem = jr; ielem < m_numElemConstraints; ielem++) { + k = m_nelem; + for (size_t ielem = jr; ielem < m_nelem; ielem++) { if (m_elementActive[ielem] && aw[ielem] != test) { k = ielem; break; } } - if (k == m_numElemConstraints) { + if (k == m_nelem) { throw CanteraError("vcs_elem_rearrange", "Shouldn't be here. Algorithm misfired."); } @@ -136,7 +136,7 @@ void VCS_SOLVE::vcs_switch_elem_pos(size_t ipos, size_t jpos) if (ipos == jpos) { return; } - AssertThrowMsg(ipos < m_numElemConstraints && jpos < m_numElemConstraints, + AssertThrowMsg(ipos < m_nelem && jpos < m_nelem, "vcs_switch_elem_pos", "inappropriate args: {} {}", ipos, jpos); @@ -158,7 +158,7 @@ void VCS_SOLVE::vcs_switch_elem_pos(size_t ipos, size_t jpos) std::swap(m_elementMapIndex[ipos], m_elementMapIndex[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) { + for (size_t j = 0; j < m_nsp; ++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 1be483f86..cdbfce978 100644 --- a/src/equil/vcs_inest.cpp +++ b/src/equil/vcs_inest.cpp @@ -19,7 +19,6 @@ static char pprefix[20] = " --- vcs_inest: "; void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm, double* const ss, double test) { - size_t nspecies = m_numSpeciesTot; size_t nrxn = m_numRxnTot; // CALL ROUTINE TO SOLVE MAX(CC*molNum) SUCH THAT AX*molNum = BB AND @@ -30,7 +29,7 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm, plogf("%s Mole Numbers returned from linear programming (vcs_inest initial guess):\n", pprefix); plogf("%s SPECIES MOLE_NUMBER -SS_ChemPotential\n", pprefix); - for (size_t kspec = 0; kspec < nspecies; ++kspec) { + for (size_t kspec = 0; kspec < m_nsp; ++kspec) { plogf("%s ", pprefix); plogf("%-12.12s", m_speciesName[kspec]); plogf(" %15.5g %12.3g\n", m_molNumSpecies_old[kspec], -m_SSfeSpecies[kspec]); @@ -38,10 +37,10 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm, plogf("%s Element Abundance Agreement returned from linear " "programming (vcs_inest initial guess):\n", pprefix); plogf("%s Element Goal Actual\n", pprefix); - for (size_t j = 0; j < m_numElemConstraints; j++) { + for (size_t j = 0; j < m_nelem; j++) { if (m_elementActive[j]) { double tmp = 0.0; - for (size_t kspec = 0; kspec < nspecies; ++kspec) { + for (size_t kspec = 0; kspec < m_nsp; ++kspec) { tmp += m_formulaMatrix(kspec,j) * m_molNumSpecies_old[kspec]; } plogf("%s ", pprefix); @@ -55,7 +54,7 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm, // Make sure all species have positive definite mole numbers Set voltages to // zero for now, until we figure out what to do m_deltaMolNumSpecies.assign(m_deltaMolNumSpecies.size(), 0.0); - for (size_t kspec = 0; kspec < nspecies; ++kspec) { + for (size_t kspec = 0; kspec < m_nsp; ++kspec) { if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { if (m_molNumSpecies_old[kspec] <= 0.0) { // HKM Should eventually include logic here for non SS phases @@ -113,7 +112,7 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm, } vcs_deltag(0, true, VCS_STATECALC_NEW); if (m_debug_print_lvl >= 2) { - for (size_t kspec = 0; kspec < nspecies; ++kspec) { + for (size_t kspec = 0; kspec < m_nsp; ++kspec) { plogf("%s", pprefix); plogf("%-12.12s", m_speciesName[kspec]); if (kspec < m_numComponents) { @@ -168,7 +167,7 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm, } } if (m_debug_print_lvl >= 2) { - for (size_t kspec = 0; kspec < nspecies; ++kspec) { + for (size_t kspec = 0; kspec < m_nsp; ++kspec) { if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { plogf("%sdirection (", pprefix); plogf("%-12.12s", m_speciesName[kspec]); @@ -212,7 +211,7 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm, m_deltaMolNumSpecies[kspec] = 0.0; } } - for (size_t kspec = m_numComponents; kspec < nspecies; ++kspec) { + for (size_t kspec = m_numComponents; kspec < m_nsp; ++kspec) { if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE && m_deltaMolNumSpecies[kspec] != 0.0) { m_molNumSpecies_old[kspec] = m_deltaMolNumSpecies[kspec] * par; @@ -228,9 +227,9 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm, // CONVERGENCE FORCING SECTION vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD); - vcs_dfe(VCS_STATECALC_OLD, 0, 0, nspecies); + vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_nsp); double s = 0.0; - for (size_t kspec = 0; kspec < nspecies; ++kspec) { + for (size_t kspec = 0; kspec < m_nsp; ++kspec) { s += m_deltaMolNumSpecies[kspec] * m_feSpecies_old[kspec]; } if (s == 0.0) { @@ -269,7 +268,7 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm, plogf("%s Final Mole Numbers produced by inest:\n", pprefix); plogf("%s SPECIES MOLE_NUMBER\n", pprefix); - for (size_t kspec = 0; kspec < nspecies; ++kspec) { + for (size_t kspec = 0; kspec < m_nsp; ++kspec) { plogf("%s %-12.12s %g\n", pprefix, m_speciesName[kspec], m_molNumSpecies_old[kspec]); } @@ -298,10 +297,10 @@ int VCS_SOLVE::vcs_inest_TP() } // temporary space for usage in this routine and in subroutines - vector_fp sm(m_numElemConstraints*m_numElemConstraints, 0.0); - vector_fp ss(m_numElemConstraints, 0.0); - vector_fp sa(m_numElemConstraints, 0.0); - vector_fp aw(m_numSpeciesTot+ m_numElemConstraints, 0.0); + vector_fp sm(m_nelem*m_nelem, 0.0); + vector_fp ss(m_nelem, 0.0); + vector_fp sa(m_nelem, 0.0); + vector_fp aw(m_nsp + m_nelem, 0.0); // Go get the estimate of the solution if (m_debug_print_lvl >= 2) { @@ -343,7 +342,7 @@ int VCS_SOLVE::vcs_inest_TP() plogf("%sElement Abundances RANGE ERROR\n", pprefix); plogf("%s - Initial guess satisfies NC=%d element abundances, " "BUT not NE=%d element abundances\n", pprefix, - m_numComponents, m_numElemConstraints); + m_numComponents, m_nelem); } } } @@ -355,7 +354,7 @@ int VCS_SOLVE::vcs_inest_TP() plogf("%sElement Abundances RANGE ERROR\n", pprefix); plogf("%s - Initial guess satisfies NC=%d element abundances, " "BUT not NE=%d element abundances\n", pprefix, - m_numComponents, m_numElemConstraints); + m_numComponents, m_nelem); } } } diff --git a/src/equil/vcs_nondim.cpp b/src/equil/vcs_nondim.cpp index 933d18c33..dbe60871f 100644 --- a/src/equil/vcs_nondim.cpp +++ b/src/equil/vcs_nondim.cpp @@ -19,7 +19,7 @@ void VCS_SOLVE::vcs_nondim_TP() if (m_unitsState == VCS_DIMENSIONAL_G) { m_unitsState = VCS_NONDIMENSIONAL_G; double tf = 1.0 / (GasConstant * m_temperature); - for (size_t i = 0; i < m_numSpeciesTot; ++i) { + for (size_t i = 0; i < m_nsp; ++i) { // Modify the standard state and total chemical potential data, // FF(I), to make it dimensionless, i.e., mu / RT. Thus, we may // divide it by the temperature. @@ -37,7 +37,7 @@ void VCS_SOLVE::vcs_nondim_TP() // Then add in the total moles of elements that are goals. Either one or // the other is specified here. double esum = 0.0; - for (size_t i = 0; i < m_numElemConstraints; ++i) { + for (size_t i = 0; i < m_nelem; ++i) { if (m_elType[i] == VCS_ELEM_TYPE_ABSPOS) { esum += fabs(m_elemAbundancesGoal[i]); } @@ -66,12 +66,12 @@ void VCS_SOLVE::vcs_nondim_TP() if (m_debug_print_lvl >= 2) { plogf(" --- vcs_nondim_TP() called: USING A MOLE SCALE OF %g until further notice\n", m_totalMoleScale); } - for (size_t i = 0; i < m_numSpeciesTot; ++i) { + for (size_t i = 0; i < m_nsp; ++i) { if (m_speciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { m_molNumSpecies_old[i] *= (1.0 / m_totalMoleScale); } } - for (size_t i = 0; i < m_numElemConstraints; ++i) { + for (size_t i = 0; i < m_nelem; ++i) { m_elemAbundancesGoal[i] *= (1.0 / m_totalMoleScale); } @@ -92,7 +92,7 @@ void VCS_SOLVE::vcs_redim_TP() if (m_unitsState != VCS_DIMENSIONAL_G) { m_unitsState = VCS_DIMENSIONAL_G; double tf = m_temperature * GasConstant; - for (size_t i = 0; i < m_numSpeciesTot; ++i) { + for (size_t i = 0; i < m_nsp; ++i) { // Modify the standard state and total chemical potential data, // FF(I), to make it have units, i.e. mu = RT * mu_star @@ -107,12 +107,12 @@ void VCS_SOLVE::vcs_redim_TP() if (m_debug_print_lvl >= 2) { plogf(" --- vcs_redim_TP() called: getting rid of mole scale of %g\n", m_totalMoleScale); } - for (size_t i = 0; i < m_numSpeciesTot; ++i) { + for (size_t i = 0; i < m_nsp; ++i) { if (m_speciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { m_molNumSpecies_old[i] *= m_totalMoleScale; } } - for (size_t i = 0; i < m_numElemConstraints; ++i) { + for (size_t i = 0; i < m_nelem; ++i) { m_elemAbundancesGoal[i] *= m_totalMoleScale; } diff --git a/src/equil/vcs_prep.cpp b/src/equil/vcs_prep.cpp index 9c93bb63c..4484a4027 100644 --- a/src/equil/vcs_prep.cpp +++ b/src/equil/vcs_prep.cpp @@ -15,7 +15,7 @@ namespace Cantera void VCS_SOLVE::vcs_SSPhase() { vector_int numPhSpecies(m_numPhases, 0); - for (size_t kspec = 0; kspec < m_numSpeciesTot; ++kspec) { + for (size_t kspec = 0; kspec < m_nsp; ++kspec) { numPhSpecies[m_phaseID[kspec]]++; } @@ -36,7 +36,7 @@ void VCS_SOLVE::vcs_SSPhase() // Fill in some useful arrays here that have to do with the static // information concerning the phase ID of species. SSPhase = Boolean // indicating whether a species is in a single species phase or not. - for (size_t kspec = 0; kspec < m_numSpeciesTot; kspec++) { + for (size_t kspec = 0; kspec < m_nsp; kspec++) { size_t iph = m_phaseID[kspec]; vcs_VolPhase* Vphase = m_VolPhaseList[iph]; if (Vphase->m_singleSpecies) { @@ -47,7 +47,7 @@ void VCS_SOLVE::vcs_SSPhase() } } -int VCS_SOLVE::vcs_prep_oneTime(int printLvl) +int VCS_SOLVE::vcs_prep(int printLvl) { int retn = VCS_SUCCESS; m_debug_print_lvl = printLvl; @@ -58,18 +58,18 @@ int VCS_SOLVE::vcs_prep_oneTime(int printLvl) // Set an initial estimate for the number of noncomponent species equal to // nspecies - nelements. This may be changed below - if (m_numElemConstraints > m_numSpeciesTot) { + if (m_nelem > m_nsp) { m_numRxnTot = 0; } else { - m_numRxnTot = m_numSpeciesTot - m_numElemConstraints; + m_numRxnTot = m_nsp - m_nelem; } m_numRxnRdc = m_numRxnTot; - m_numSpeciesRdc = m_numSpeciesTot; + m_numSpeciesRdc = m_nsp; for (size_t i = 0; i < m_numRxnRdc; ++i) { - m_indexRxnToSpecies[i] = m_numElemConstraints + i; + m_indexRxnToSpecies[i] = m_nelem + i; } - for (size_t kspec = 0; kspec < m_numSpeciesTot; ++kspec) { + for (size_t kspec = 0; kspec < m_nsp; ++kspec) { size_t pID = m_phaseID[kspec]; size_t spPhIndex = m_speciesLocalPhaseIndex[kspec]; vcs_VolPhase* vPhase = m_VolPhaseList[pID]; @@ -104,7 +104,7 @@ int VCS_SOLVE::vcs_prep_oneTime(int printLvl) bool modifiedSoln = false; if (m_doEstimateEquil < 0) { double sum = 0.0; - for (size_t kspec = 0; kspec < m_numSpeciesTot; ++kspec) { + for (size_t kspec = 0; kspec < m_nsp; ++kspec) { if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_MOLNUM) { sum += fabs(m_molNumSpecies_old[kspec]); } @@ -113,7 +113,7 @@ int VCS_SOLVE::vcs_prep_oneTime(int printLvl) modifiedSoln = true; double pres = (m_pressurePA <= 0.0) ? 1.01325E5 : m_pressurePA; retn = vcs_evalSS_TP(0, 0, m_temperature, pres); - for (size_t kspec = 0; kspec < m_numSpeciesTot; ++kspec) { + for (size_t kspec = 0; kspec < m_nsp; ++kspec) { if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_MOLNUM) { m_molNumSpecies_old[kspec] = - m_SSfeSpecies[kspec]; } else { @@ -126,15 +126,15 @@ int VCS_SOLVE::vcs_prep_oneTime(int printLvl) // NC = number of components is in the vcs.h common block. This call to // BASOPT doesn't calculate the stoichiometric reaction matrix. - vector_fp awSpace(m_numSpeciesTot + (m_numElemConstraints + 2)*(m_numElemConstraints), 0.0); + vector_fp awSpace(m_nsp + (m_nelem + 2)*(m_nelem), 0.0); double* aw = &awSpace[0]; if (aw == NULL) { plogf("vcs_prep_oneTime: failed to get memory: global bailout\n"); return VCS_NOMEMORY; } - double* sa = aw + m_numSpeciesTot; - double* sm = sa + m_numElemConstraints; - double* ss = sm + (m_numElemConstraints)*(m_numElemConstraints); + double* sa = aw + m_nsp; + double* sm = sa + m_nelem; + double* ss = sm + m_nelem * m_nelem; bool conv; retn = vcs_basopt(true, aw, sa, sm, ss, test, &conv); if (retn != VCS_SUCCESS) { @@ -145,8 +145,8 @@ int VCS_SOLVE::vcs_prep_oneTime(int printLvl) return retn; } - if (m_numSpeciesTot >= m_numComponents) { - m_numRxnTot = m_numRxnRdc = m_numSpeciesTot - m_numComponents; + if (m_nsp >= m_numComponents) { + m_numRxnTot = m_numRxnRdc = m_nsp - m_numComponents; for (size_t i = 0; i < m_numRxnRdc; ++i) { m_indexRxnToSpecies[i] = m_numComponents + i; } @@ -155,11 +155,11 @@ int VCS_SOLVE::vcs_prep_oneTime(int printLvl) } // The elements might need to be rearranged. - awSpace.resize(m_numElemConstraints + (m_numElemConstraints + 2)*(m_numElemConstraints), 0.0); + awSpace.resize(m_nelem + (m_nelem + 2)*m_nelem, 0.0); aw = &awSpace[0]; - sa = aw + m_numElemConstraints; - sm = sa + m_numElemConstraints; - ss = sm + (m_numElemConstraints)*(m_numElemConstraints); + sa = aw + m_nelem; + sm = sa + m_nelem; + ss = sm + m_nelem * m_nelem; retn = vcs_elem_rearrange(aw, sa, sm, ss); if (retn != VCS_SUCCESS) { plogf("vcs_prep_oneTime:"); @@ -172,15 +172,11 @@ int VCS_SOLVE::vcs_prep_oneTime(int printLvl) // If we mucked up the solution unknowns because they were all // zero to start with, set them back to zero here if (modifiedSoln) { - for (size_t kspec = 0; kspec < m_numSpeciesTot; ++kspec) { + for (size_t kspec = 0; kspec < m_nsp; ++kspec) { m_molNumSpecies_old[kspec] = 0.0; } } - return VCS_SUCCESS; -} -int VCS_SOLVE::vcs_prep() -{ // Initialize various arrays in the data to zero m_feSpecies_old.assign(m_feSpecies_old.size(), 0.0); m_feSpecies_new.assign(m_feSpecies_new.size(), 0.0); @@ -192,20 +188,18 @@ int VCS_SOLVE::vcs_prep() // Calculate the total number of moles in all phases. vcs_tmoles(); - return VCS_SUCCESS; -} -bool VCS_SOLVE::vcs_wellPosed() -{ + // Check to see if the current problem is well posed. double sum = 0.0; - for (size_t e = 0; e < ne; e++) { + for (size_t e = 0; e < m_nelem; e++) { sum += m_mix->elementMoles(e); } if (sum < 1.0E-20) { - plogf("vcs_wellPosed: Element abundance is close to zero\n"); - return false; + // Check to see if the current problem is well posed. + plogf("vcs has determined the problem is not well posed: Bailing\n"); + return VCS_PUB_BAD; } - return true; + return VCS_SUCCESS; } } diff --git a/src/equil/vcs_prob.cpp b/src/equil/vcs_prob.cpp index 3ceb13708..82cd8d39f 100644 --- a/src/equil/vcs_prob.cpp +++ b/src/equil/vcs_prob.cpp @@ -25,8 +25,8 @@ namespace Cantera void VCS_SOLVE::set_gai() { gai.assign(gai.size(), 0.0); - for (size_t j = 0; j < ne; j++) { - for (size_t kspec = 0; kspec < nspecies; kspec++) { + for (size_t j = 0; j < m_nelem; j++) { + for (size_t kspec = 0; kspec < m_nsp; kspec++) { if (SpeciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { gai[j] += FormulaMatrix(kspec,j) * w[kspec]; } @@ -46,24 +46,20 @@ void VCS_SOLVE::prob_report(int print_lvl) writeline('=', 31); writeline('=', 80); plogf("\n"); - if (prob_type == 0) { - plogf("\tSolve a constant T, P problem:\n"); - plogf("\t\tT = %g K\n", T); - double pres_atm = PresPA / 1.01325E5; + plogf("\tSolve a constant T, P problem:\n"); + plogf("\t\tT = %g K\n", m_temperature); + double pres_atm = m_pressurePA / 1.01325E5; - plogf("\t\tPres = %g atm\n", pres_atm); - } else { - throw CanteraError("VCS_PROB::prob_report", "Unknown problem type"); - } + plogf("\t\tPres = %g atm\n", pres_atm); plogf("\n"); plogf(" Phase IDs of species\n"); plogf(" species phaseID phaseName "); plogf(" Initial_Estimated_Moles Species_Type\n"); - for (size_t i = 0; i < nspecies; i++) { + for (size_t i = 0; i < m_nsp; i++) { vcs_VolPhase* Vphase = VPhaseList[PhaseID[i]]; plogf("%16s %5d %16s", m_mix->speciesName(i), PhaseID[i], Vphase->PhaseName); - if (iest >= 0) { + if (m_doEstimateEquil >= 0) { plogf(" %-10.5g", w[i]); } else { plogf(" N/A"); @@ -85,13 +81,13 @@ void VCS_SOLVE::prob_report(int print_lvl) " EqnState NumSpec"); plogf(" TMolesInert TKmoles\n"); - for (size_t iphase = 0; iphase < NPhase; iphase++) { + for (size_t iphase = 0; iphase < m_numPhases; iphase++) { vcs_VolPhase* Vphase = VPhaseList[iphase]; plogf("%16s %5d %5d %8d ", Vphase->PhaseName, Vphase->VP_ID_, Vphase->m_singleSpecies, Vphase->m_gasPhase); plogf("%16s %8d %16e ", Vphase->eos_name(), Vphase->nSpecies(), Vphase->totalMolesInert()); - if (iest >= 0) { + if (m_doEstimateEquil >= 0) { plogf("%16e\n", Vphase->totalMoles()); } else { plogf(" N/A\n"); @@ -100,7 +96,7 @@ void VCS_SOLVE::prob_report(int print_lvl) plogf("\nElemental Abundances: "); plogf(" Target_kmol ElemType ElActive\n"); - for (size_t i = 0; i < ne; ++i) { + for (size_t i = 0; i < m_nelem; ++i) { writeline(' ', 26, false); plogf("%-2.2s", m_elementName[i]); plogf("%20.12E ", gai[i]); @@ -110,9 +106,9 @@ void VCS_SOLVE::prob_report(int print_lvl) plogf("\nChemical Potentials: (J/kmol)\n"); plogf(" Species (phase) " " SS0ChemPot StarChemPot\n"); - for (size_t iphase = 0; iphase < NPhase; iphase++) { + for (size_t iphase = 0; iphase < m_numPhases; iphase++) { vcs_VolPhase* Vphase = VPhaseList[iphase]; - Vphase->setState_TP(T, PresPA); + Vphase->setState_TP(m_temperature, m_pressurePA); for (size_t kindex = 0; kindex < Vphase->nSpecies(); kindex++) { size_t kglob = Vphase->spGlobalIndexVCS(kindex); plogf("%16s ", m_mix->speciesName(kglob)); @@ -145,7 +141,7 @@ void VCS_SOLVE::addPhaseElements(vcs_VolPhase* volPhase) // Search for matches with the existing elements. If found, then fill in // the entry in the global mapping array. - for (size_t e = 0; e < ne; e++) { + for (size_t e = 0; e < m_nelem; e++) { std::string en = m_elementName[e]; if (!strcmp(enVP.c_str(), en.c_str())) { volPhase->setElemGlobalIndex(eVP, e); @@ -167,14 +163,13 @@ size_t VCS_SOLVE::addElement(const char* elNameNew, int elType, int elactive) throw CanteraError("VCS_SOLVE::addElement", "error: element must have a name"); } - ne++; - m_numElemConstraints++; + m_nelem++; m_numComponents++; gai.push_back(0.0); - FormulaMatrix.resize(NSPECIES0, ne, 0.0); - m_formulaMatrix.resize(NSPECIES0, ne); - m_stoichCoeffRxnMatrix.resize(ne, NSPECIES0, 0.0); + FormulaMatrix.resize(m_nsp, m_nelem, 0.0); + m_formulaMatrix.resize(m_nsp, m_nelem); + m_stoichCoeffRxnMatrix.resize(m_nelem, m_nsp, 0.0); m_elType.push_back(elType); ElActive.push_back(elactive); m_elementActive.push_back(elactive); @@ -182,13 +177,12 @@ size_t VCS_SOLVE::addElement(const char* elNameNew, int elType, int elactive) m_elemAbundancesGoal.push_back(0.0); m_elementMapIndex.push_back(0); m_elementName.push_back(elNameNew); - NE0 = ne; - return ne - 1; + return m_nelem - 1; } size_t VCS_SOLVE::addOnePhaseSpecies(vcs_VolPhase* volPhase, size_t k, size_t kT) { - if (kT > nspecies) { + if (kT > m_nsp) { // Need to expand the number of species here throw CanteraError("VCS_PROB::addOnePhaseSpecies", "Shouldn't be here"); } diff --git a/src/equil/vcs_report.cpp b/src/equil/vcs_report.cpp index 95a38a0e8..acc35cbfb 100644 --- a/src/equil/vcs_report.cpp +++ b/src/equil/vcs_report.cpp @@ -12,13 +12,12 @@ namespace Cantera int VCS_SOLVE::vcs_report(int iconv) { bool printActualMoles = true, inertYes = false; - size_t nspecies = m_numSpeciesTot; char originalUnitsState = m_unitsState; - std::vector sortindex(nspecies,0); - vector_fp xy(nspecies,0.0); + std::vector sortindex(m_nsp, 0); + vector_fp xy(m_nsp, 0.0); // SORT DEPENDENT SPECIES IN DECREASING ORDER OF MOLES - for (size_t i = 0; i < nspecies; ++i) { + for (size_t i = 0; i < m_nsp; ++i) { sortindex[i] = i; xy[i] = m_molNumSpecies_old[i]; } @@ -45,7 +44,7 @@ int VCS_SOLVE::vcs_report(int iconv) } vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD); - vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_numSpeciesTot); + vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_nsp); // PRINT OUT RESULTS plogf("\n\n\n\n"); @@ -117,9 +116,9 @@ int VCS_SOLVE::vcs_report(int iconv) TPhInertMoles[i] / m_tPhaseMoles_old[i], 0.0); } } - if (m_numSpeciesRdc != nspecies) { + if (m_numSpeciesRdc != m_nsp) { plogf("\n SPECIES WITH LESS THAN 1.0E-32 KMOLES:\n\n"); - for (size_t kspec = m_numSpeciesRdc; kspec < nspecies; ++kspec) { + for (size_t kspec = m_numSpeciesRdc; kspec < m_nsp; ++kspec) { plogf(" %-12.12s", m_speciesName[kspec]); // Note m_deltaGRxn_new[] stores in kspec slot not irxn slot, after solve plogf(" %14.7E %14.7E %12.4E", @@ -171,30 +170,30 @@ int VCS_SOLVE::vcs_report(int iconv) plogf("\n"); // TABLE OF PHASE INFORMATION - vector_fp gaPhase(m_numElemConstraints, 0.0); - vector_fp gaTPhase(m_numElemConstraints, 0.0); + vector_fp gaPhase(m_nelem, 0.0); + vector_fp gaTPhase(m_nelem, 0.0); double totalMoles = 0.0; double gibbsPhase = 0.0; double gibbsTotal = 0.0; plogf("\n\n"); plogf("\n"); - writeline('-', m_numElemConstraints*10 + 58); + writeline('-', m_nelem*10 + 58); plogf(" | ElementID |"); - for (size_t j = 0; j < m_numElemConstraints; j++) { + for (size_t j = 0; j < m_nelem; j++) { plogf(" %3d", j); } plogf(" | |\n"); plogf(" | Element |"); - for (size_t j = 0; j < m_numElemConstraints; j++) { + for (size_t j = 0; j < m_nelem; j++) { plogf(" %10.10s", m_elementName[j]); } plogf(" | |\n"); plogf(" PhaseName |KMolTarget |"); - for (size_t j = 0; j < m_numElemConstraints; j++) { + for (size_t j = 0; j < m_nelem; j++) { plogf(" %10.3g", m_elemAbundancesGoal[j]); } plogf(" | Gibbs Total |\n"); - writeline('-', m_numElemConstraints*10 + 58); + writeline('-', m_nelem*10 + 58); for (size_t iphase = 0; iphase < m_numPhases; iphase++) { plogf(" %3d ", iphase); vcs_VolPhase* VPhase = m_VolPhaseList[iphase]; @@ -206,7 +205,7 @@ int VCS_SOLVE::vcs_report(int iconv) throw CanteraError("VCS_SOLVE::vcs_report", "we have a problem"); } vcs_elabPhase(iphase, &gaPhase[0]); - for (size_t j = 0; j < m_numElemConstraints; j++) { + for (size_t j = 0; j < m_nelem; j++) { plogf(" %10.3g", gaPhase[j]); gaTPhase[j] += gaPhase[j]; } @@ -215,14 +214,14 @@ int VCS_SOLVE::vcs_report(int iconv) gibbsTotal += gibbsPhase; plogf(" | %18.11E |\n", gibbsPhase); } - writeline('-', m_numElemConstraints*10 + 58); + writeline('-', m_nelem*10 + 58); plogf(" TOTAL |%10.3e |", totalMoles); - for (size_t j = 0; j < m_numElemConstraints; j++) { + for (size_t j = 0; j < m_nelem; j++) { plogf(" %10.3g", gaTPhase[j]); } plogf(" | %18.11E |\n", gibbsTotal); - writeline('-', m_numElemConstraints*10 + 58); + writeline('-', m_nelem*10 + 58); plogf("\n"); // GLOBAL SATISFACTION INFORMATION @@ -238,7 +237,7 @@ int VCS_SOLVE::vcs_report(int iconv) plogf("\nElemental Abundances (kmol): "); plogf(" Actual Target Type ElActive\n"); - for (size_t i = 0; i < m_numElemConstraints; ++i) { + for (size_t i = 0; i < m_nelem; ++i) { writeline(' ', 26, false); plogf("%-2.2s", m_elementName[i]); plogf("%20.12E %20.12E", m_elemAbundances[i]*molScale, m_elemAbundancesGoal[i]*molScale); @@ -255,7 +254,7 @@ int VCS_SOLVE::vcs_report(int iconv) " ln(AC) ln(X_i) | F z_i phi | ChemPot | (-lnMnaught)"); plogf("| (MolNum ChemPot)|"); writeline('-', 147, true, true); - for (size_t i = 0; i < nspecies; ++i) { + for (size_t i = 0; i < m_nsp; ++i) { size_t j = sortindex[i]; size_t pid = m_phaseID[j]; plogf(" %-12.12s", m_speciesName[j]); diff --git a/src/equil/vcs_rxnadj.cpp b/src/equil/vcs_rxnadj.cpp index 97245a94a..2d3fd8fed 100644 --- a/src/equil/vcs_rxnadj.cpp +++ b/src/equil/vcs_rxnadj.cpp @@ -242,7 +242,7 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) } // Delete the single species phase - for (size_t j = 0; j < m_numSpeciesTot; j++) { + for (size_t j = 0; j < m_nsp; j++) { m_deltaMolNumSpecies[j] = 0.0; } m_deltaMolNumSpecies[kspec] = dss; diff --git a/src/equil/vcs_setMolesLinProg.cpp b/src/equil/vcs_setMolesLinProg.cpp index eac8c4603..1e0d713a5 100644 --- a/src/equil/vcs_setMolesLinProg.cpp +++ b/src/equil/vcs_setMolesLinProg.cpp @@ -40,13 +40,13 @@ int VCS_SOLVE::vcs_setMolesLinProg() int iter = 0; bool abundancesOK = true; bool usedZeroedSpecies; - vector_fp sm(m_numElemConstraints*m_numElemConstraints, 0.0); - vector_fp ss(m_numElemConstraints, 0.0); - vector_fp sa(m_numElemConstraints, 0.0); - vector_fp wx(m_numElemConstraints, 0.0); - vector_fp aw(m_numSpeciesTot, 0.0); + vector_fp sm(m_nelem * m_nelem, 0.0); + vector_fp ss(m_nelem, 0.0); + vector_fp sa(m_nelem, 0.0); + vector_fp wx(m_nelem, 0.0); + vector_fp aw(m_nsp, 0.0); - for (size_t ik = 0; ik < m_numSpeciesTot; ik++) { + for (size_t ik = 0; ik < m_nsp; ik++) { if (m_speciesUnknownType[ik] != VCS_SPECIES_INTERFACIALVOLTAGE) { m_molNumSpecies_old[ik] = max(0.0, m_molNumSpecies_old[ik]); } @@ -98,7 +98,7 @@ int VCS_SOLVE::vcs_setMolesLinProg() double dg_rt = m_SSfeSpecies[ik]; dxi_min = 1.0e10; const double* sc_irxn = m_stoichCoeffRxnMatrix.ptrColumn(irxn); - for (size_t jcomp = 0; jcomp < m_numElemConstraints; jcomp++) { + for (size_t jcomp = 0; jcomp < m_nelem; jcomp++) { dg_rt += m_SSfeSpecies[jcomp] * sc_irxn[jcomp]; } // fwd or rev direction. diff --git a/src/equil/vcs_solve.cpp b/src/equil/vcs_solve.cpp index 023425194..a295fc241 100644 --- a/src/equil/vcs_solve.cpp +++ b/src/equil/vcs_solve.cpp @@ -21,75 +21,58 @@ namespace Cantera int vcs_timing_print_lvl = 1; -VCS_SOLVE::VCS_SOLVE(size_t nspecies_, size_t nelements_, size_t nphase_) : - prob_type(VCS_PROBTYPE_TP), - nspecies(nspecies_), - ne(0), - NE0(0), - NPhase(nphase_), - T(298.15), - PresPA(1.0), - Vol(0.0), - // The default is to not expect an initial estimate of the species - // concentrations - iest(-1), - tolmaj(1.0E-8), - tolmin(1.0E-6), - m_Iterations(0), - m_NumBasisOptimizations(0), - m_printLvl(0), +VCS_SOLVE::VCS_SOLVE(MultiPhase* mphase, int printLvl) : + m_printLvl(printLvl), vcs_debug_print_lvl(0), - NSPECIES0(nspecies_), - NPHASE0(nphase_), - m_numSpeciesTot(nspecies_), - m_numElemConstraints(0), + m_mix(mphase), + m_nsp(mphase->nSpecies()), + m_nelem(0), m_numComponents(0), m_numRxnTot(0), m_numSpeciesRdc(0), m_numRxnRdc(0), m_numRxnMinorZeroed(0), - m_numPhases(0), - m_doEstimateEquil(0), + m_numPhases(mphase->nPhases()), + m_doEstimateEquil(-1), m_totalMolNum(0.0), - m_temperature(0.0), - m_pressurePA(0.0), - m_tolmaj(0.0), - m_tolmin(0.0), - m_tolmaj2(0.0), - m_tolmin2(0.0), + m_temperature(mphase->temperature()), + m_pressurePA(mphase->pressure()), + m_tolmaj(1.0E-8), + m_tolmin(1.0E-6), + m_tolmaj2(1.0E-10), + m_tolmin2(1.0E-8), m_unitsState(VCS_DIMENSIONAL_G), m_totalMoleScale(1.0), m_useActCoeffJac(0), - m_totalVol(0.0), + m_totalVol(mphase->volume()), m_Faraday_dim(ElectronCharge * Avogadro), m_VCount(0), m_debug_print_lvl(0), m_timing_print_lvl(1) { - m_gibbsSpecies.resize(nspecies, 0.0); - w.resize(nspecies, 0.0); - mf.resize(nspecies, 0.0); - SpeciesUnknownType.resize(nspecies, VCS_SPECIES_TYPE_MOLNUM); - VolPM.resize(nspecies, 0.0); - PhaseID.resize(nspecies, npos); - WtSpecies.resize(nspecies, 0.0); - Charge.resize(nspecies, 0.0); - SpeciesThermo.resize(nspecies,0); - for (size_t kspec = 0; kspec < nspecies; kspec++) { + m_gibbsSpecies.resize(m_nsp, 0.0); + w.resize(m_nsp, 0.0); + mf.resize(m_nsp, 0.0); + SpeciesUnknownType.resize(m_nsp, VCS_SPECIES_TYPE_MOLNUM); + PhaseID.resize(m_nsp, npos); + WtSpecies.resize(m_nsp, 0.0); + Charge.resize(m_nsp, 0.0); + SpeciesThermo.resize(m_nsp,0); + for (size_t kspec = 0; kspec < m_nsp; kspec++) { SpeciesThermo[kspec] = new VCS_SPECIES_THERMO(0, 0); } - VPhaseList.resize(NPhase, 0); - for (size_t iphase = 0; iphase < NPhase; iphase++) { + VPhaseList.resize(m_numPhases, 0); + for (size_t iphase = 0; iphase < m_numPhases; iphase++) { VPhaseList[iphase] = new vcs_VolPhase(); } string ser = "VCS_SOLVE: ERROR:\n\t"; - if (nspecies_ <= 0) { + if (m_nsp <= 0) { plogf("%s Number of species is nonpositive\n", ser); throw CanteraError("VCS_SOLVE()", ser + " Number of species is nonpositive\n"); } - if (nphase_ <= 0) { + if (m_numPhases <= 0) { plogf("%s Number of phases is nonpositive\n", ser); throw CanteraError("VCS_SOLVE()", ser + " Number of species is nonpositive\n"); @@ -99,68 +82,68 @@ VCS_SOLVE::VCS_SOLVE(size_t nspecies_, size_t nelements_, size_t nphase_) : * We will initialize sc[] to note the fact that it needs to be * filled with meaningful information. */ - m_scSize.resize(NSPECIES0, 0.0); - m_spSize.resize(NSPECIES0, 1.0); - m_SSfeSpecies.resize(NSPECIES0, 0.0); - m_feSpecies_new.resize(NSPECIES0, 0.0); - m_molNumSpecies_old.resize(NSPECIES0, 0.0); - m_speciesUnknownType.resize(NSPECIES0, VCS_SPECIES_TYPE_MOLNUM); - m_deltaMolNumPhase.resize(NPHASE0, NSPECIES0, 0.0); - m_phaseParticipation.resize(NPHASE0, NSPECIES0, 0); - m_phasePhi.resize(NPHASE0, 0.0); - m_molNumSpecies_new.resize(NSPECIES0, 0.0); - m_deltaGRxn_new.resize(NSPECIES0, 0.0); - m_deltaGRxn_old.resize(NSPECIES0, 0.0); - m_deltaGRxn_Deficient.resize(NSPECIES0, 0.0); - m_deltaGRxn_tmp.resize(NSPECIES0, 0.0); - m_deltaMolNumSpecies.resize(NSPECIES0, 0.0); - m_feSpecies_old.resize(NSPECIES0, 0.0); - m_tPhaseMoles_old.resize(NPHASE0, 0.0); - m_tPhaseMoles_new.resize(NPHASE0, 0.0); - m_deltaPhaseMoles.resize(NPHASE0, 0.0); - m_TmpPhase.resize(NPHASE0, 0.0); - m_TmpPhase2.resize(NPHASE0, 0.0); - TPhInertMoles.resize(NPHASE0, 0.0); + m_scSize.resize(m_nsp, 0.0); + m_spSize.resize(m_nsp, 1.0); + m_SSfeSpecies.resize(m_nsp, 0.0); + m_feSpecies_new.resize(m_nsp, 0.0); + m_molNumSpecies_old.resize(m_nsp, 0.0); + m_speciesUnknownType.resize(m_nsp, VCS_SPECIES_TYPE_MOLNUM); + m_deltaMolNumPhase.resize(m_numPhases, m_nsp, 0.0); + m_phaseParticipation.resize(m_numPhases, m_nsp, 0); + m_phasePhi.resize(m_numPhases, 0.0); + m_molNumSpecies_new.resize(m_nsp, 0.0); + m_deltaGRxn_new.resize(m_nsp, 0.0); + m_deltaGRxn_old.resize(m_nsp, 0.0); + m_deltaGRxn_Deficient.resize(m_nsp, 0.0); + m_deltaGRxn_tmp.resize(m_nsp, 0.0); + m_deltaMolNumSpecies.resize(m_nsp, 0.0); + m_feSpecies_old.resize(m_nsp, 0.0); + m_tPhaseMoles_old.resize(m_numPhases, 0.0); + m_tPhaseMoles_new.resize(m_numPhases, 0.0); + m_deltaPhaseMoles.resize(m_numPhases, 0.0); + m_TmpPhase.resize(m_numPhases, 0.0); + m_TmpPhase2.resize(m_numPhases, 0.0); + TPhInertMoles.resize(m_numPhases, 0.0); // ind[] is an index variable that keep track of solution vector rotations. - m_speciesMapIndex.resize(NSPECIES0, 0); - m_speciesLocalPhaseIndex.resize(NSPECIES0, 0); + m_speciesMapIndex.resize(m_nsp, 0); + m_speciesLocalPhaseIndex.resize(m_nsp, 0); // ir[] is an index vector that keeps track of the irxn to species mapping. // We can't fill it in until we know the number of c components in the // problem - m_indexRxnToSpecies.resize(NSPECIES0, 0); + m_indexRxnToSpecies.resize(m_nsp, 0); // Initialize all species to be major species - m_speciesStatus.resize(NSPECIES0, 1); + m_speciesStatus.resize(m_nsp, 1); - m_SSPhase.resize(2*NSPECIES0, 0); - m_phaseID.resize(NSPECIES0, 0); - m_speciesName.resize(NSPECIES0); + m_SSPhase.resize(2*m_nsp, 0); + m_phaseID.resize(m_nsp, 0); + m_speciesName.resize(m_nsp); // space for activity coefficients for all species. Set it equal to one. - m_actConventionSpecies.resize(NSPECIES0, 0); - m_phaseActConvention.resize(NPHASE0, 0); - m_lnMnaughtSpecies.resize(NSPECIES0, 0.0); - m_actCoeffSpecies_new.resize(NSPECIES0, 1.0); - m_actCoeffSpecies_old.resize(NSPECIES0, 1.0); - m_wtSpecies.resize(NSPECIES0, 0.0); - m_chargeSpecies.resize(NSPECIES0, 0.0); - m_speciesThermoList.resize(NSPECIES0, (VCS_SPECIES_THERMO*)0); + m_actConventionSpecies.resize(m_nsp, 0); + m_phaseActConvention.resize(m_numPhases, 0); + m_lnMnaughtSpecies.resize(m_nsp, 0.0); + m_actCoeffSpecies_new.resize(m_nsp, 1.0); + m_actCoeffSpecies_old.resize(m_nsp, 1.0); + m_wtSpecies.resize(m_nsp, 0.0); + m_chargeSpecies.resize(m_nsp, 0.0); + m_speciesThermoList.resize(m_nsp, (VCS_SPECIES_THERMO*)0); // Phase Info - m_VolPhaseList.resize(NPHASE0, 0); - for (size_t iph = 0; iph < NPHASE0; iph++) { + m_VolPhaseList.resize(m_numPhases, 0); + for (size_t iph = 0; iph < m_numPhases; iph++) { m_VolPhaseList[iph] = new vcs_VolPhase(this); } // For Future expansion m_useActCoeffJac = true; if (m_useActCoeffJac) { - m_np_dLnActCoeffdMolNum.resize(NSPECIES0, NSPECIES0, 0.0); + m_np_dLnActCoeffdMolNum.resize(m_nsp, m_nsp, 0.0); } - m_PMVolumeSpecies.resize(NSPECIES0, 0.0); + m_PMVolumeSpecies.resize(m_nsp, 0.0); // counters kept within vcs m_VCount = new VCS_COUNTERS(); @@ -170,7 +153,465 @@ VCS_SOLVE::VCS_SOLVE(size_t nspecies_, size_t nelements_, size_t nphase_) : m_timing_print_lvl = 0; } - return; + VCS_SPECIES_THERMO* ts_ptr = 0; + + // Loop over the phases, transferring pertinent information + int kT = 0; + for (size_t iphase = 0; iphase < m_numPhases; iphase++) { + // Get the ThermoPhase object - assume volume phase + ThermoPhase* tPhase = &mphase->phase(iphase); + size_t nelem = tPhase->nElements(); + + // Query Cantera for the equation of state type of the current phase. + std::string eos = tPhase->type(); + bool gasPhase = (eos == "IdealGas"); + + // Find out the number of species in the phase + size_t nSpPhase = tPhase->nSpecies(); + // Find out the name of the phase + string phaseName = tPhase->name(); + + // Call the basic vcs_VolPhase creation routine. + // Properties set here: + // ->PhaseNum = phase number in the thermo problem + // ->GasPhase = Boolean indicating whether it is a gas phase + // ->NumSpecies = number of species in the phase + // ->TMolesInert = Inerts in the phase = 0.0 for cantera + // ->PhaseName = Name of the phase + vcs_VolPhase* VolPhase = VPhaseList[iphase]; + VolPhase->resize(iphase, nSpPhase, nelem, phaseName.c_str(), 0.0); + VolPhase->m_gasPhase = gasPhase; + + // Tell the vcs_VolPhase pointer about cantera + VolPhase->setPtrThermoPhase(tPhase); + VolPhase->setTotalMoles(0.0); + + // Set the electric potential of the volume phase from the + // ThermoPhase object's value. + VolPhase->setElectricPotential(tPhase->electricPotential()); + + // Query the ThermoPhase object to find out what convention + // it uses for the specification of activity and Standard State. + VolPhase->p_activityConvention = tPhase->activityConvention(); + + // Assign the value of eqn of state. Handle conflicts here. + if (eos == "IdealGas") { + VolPhase->m_eqnState = VCS_EOS_IDEAL_GAS; + } else if (eos == "ConstDensity") { + VolPhase->m_eqnState = VCS_EOS_CONSTANT; + } else if (eos == "StoichSubstance") { + VolPhase->m_eqnState = VCS_EOS_STOICH_SUB; + } else if (eos == "IdealSolidSoln") { + VolPhase->m_eqnState = VCS_EOS_IDEAL_SOLN; + } else if (eos == "Surf" || eos == "Edge") { + throw CanteraError("VCSnonideal", + "Surface/edge phase not handled yet."); + } else { + if (m_printLvl > 1) { + writelog("Unknown Cantera EOS to VCSnonideal: '{}'\n", eos); + } + VolPhase->m_eqnState = VCS_EOS_UNK_CANTERA; + } + + // Transfer all of the element information from the ThermoPhase object + // to the vcs_VolPhase object. Also decide whether we need a new charge + // neutrality element in the phase to enforce a charge neutrality + // constraint. We also decide whether this is a single species phase + // with the voltage being the independent variable setting the chemical + // potential of the electrons. + VolPhase->transferElementsFM(tPhase); + + // Combine the element information in the vcs_VolPhase + // object into the vprob object. + addPhaseElements(VolPhase); + VolPhase->setState_TP(m_temperature, m_pressurePA); + vector_fp muPhase(tPhase->nSpecies(),0.0); + tPhase->getChemPotentials(&muPhase[0]); + double tMoles = 0.0; + + // Loop through each species in the current phase + for (size_t k = 0; k < nSpPhase; k++) { + // Obtain the molecular weight of the species from the + // ThermoPhase object + WtSpecies[kT] = tPhase->molecularWeight(k); + + // Obtain the charges of the species from the ThermoPhase object + Charge[kT] = tPhase->charge(k); + + // Set the phaseid of the species + PhaseID[kT] = iphase; + + // Transfer the type of unknown + SpeciesUnknownType[kT] = VolPhase->speciesUnknownType(k); + if (SpeciesUnknownType[kT] == VCS_SPECIES_TYPE_MOLNUM) { + // Set the initial number of kmoles of the species + // and the mole fraction vector + w[kT] = mphase->speciesMoles(kT); + tMoles += w[kT]; + mf[kT] = mphase->moleFraction(kT); + } else if (SpeciesUnknownType[kT] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { + w[kT] = tPhase->electricPotential(); + mf[kT] = mphase->moleFraction(kT); + } else { + throw CanteraError(" vcs_Cantera_to_vsolve() ERROR", + "Unknown species type: {}", SpeciesUnknownType[kT]); + } + + // transfer chemical potential vector + m_gibbsSpecies[kT] = muPhase[k]; + + // Transfer the species information from the + // volPhase structure to the VPROB structure + // This includes: + // FormulaMatrix[][] + // VolPhase->IndSpecies[] + addOnePhaseSpecies(VolPhase, k, kT); + + // Get a pointer to the thermo object + ts_ptr = SpeciesThermo[kT]; + + // Fill in the vcs_SpeciesProperty structure + vcs_SpeciesProperties* sProp = VolPhase->speciesProperty(k); + sProp->NumElements = m_nelem; + sProp->SpName = mphase->speciesName(kT); + sProp->SpeciesThermo = ts_ptr; + sProp->WtSpecies = tPhase->molecularWeight(k); + sProp->FormulaMatrixCol.resize(m_nelem, 0.0); + for (size_t e = 0; e < m_nelem; e++) { + sProp->FormulaMatrixCol[e] = FormulaMatrix(kT,e); + } + sProp->Charge = tPhase->charge(k); + sProp->SurfaceSpecies = false; + sProp->VolPM = 0.0; + + // Transfer the thermo specification of the species + // vsolve->SpeciesThermo[] + + // Add lookback connectivity into the thermo object first + ts_ptr->IndexPhase = iphase; + ts_ptr->IndexSpeciesPhase = k; + ts_ptr->OwningPhase = VolPhase; + + // get a reference to the Cantera species thermo. + MultiSpeciesThermo& sp = tPhase->speciesThermo(); + + int spType = sp.reportType(k); + if (spType == SIMPLE) { + double c[4]; + double minTemp, maxTemp, refPressure; + sp.reportParams(k, spType, c, minTemp, maxTemp, refPressure); + ts_ptr->SS0_Model = VCS_SS0_CONSTANT; + ts_ptr->SS0_T0 = c[0]; + ts_ptr->SS0_H0 = c[1]; + ts_ptr->SS0_S0 = c[2]; + ts_ptr->SS0_Cp0 = c[3]; + if (gasPhase) { + ts_ptr->SSStar_Model = VCS_SSSTAR_IDEAL_GAS; + ts_ptr->SSStar_Vol_Model = VCS_SSVOL_IDEALGAS; + } else { + ts_ptr->SSStar_Model = VCS_SSSTAR_CONSTANT; + ts_ptr->SSStar_Vol_Model = VCS_SSVOL_CONSTANT; + } + } else { + if (m_printLvl > 2) { + plogf("vcs_Cantera_convert: Species Type %d not known \n", + spType); + } + ts_ptr->SS0_Model = VCS_SS0_NOTHANDLED; + ts_ptr->SSStar_Model = VCS_SSSTAR_NOTHANDLED; + } + + // Transfer the Volume Information -> NEEDS WORK + if (gasPhase) { + ts_ptr->SSStar_Vol_Model = VCS_SSVOL_IDEALGAS; + ts_ptr->SSStar_Vol0 = 82.05 * 273.15 / 1.0; + } else { + vector_fp phaseTermCoeff(nSpPhase, 0.0); + int nCoeff; + tPhase->getParameters(nCoeff, &phaseTermCoeff[0]); + ts_ptr->SSStar_Vol_Model = VCS_SSVOL_CONSTANT; + ts_ptr->SSStar_Vol0 = phaseTermCoeff[k]; + } + kT++; + } + + // Now go back through the species in the phase and assign a valid mole + // fraction to all phases, even if the initial estimate of the total + // number of moles is zero. + if (tMoles > 0.0) { + for (size_t k = 0; k < nSpPhase; k++) { + size_t kTa = VolPhase->spGlobalIndexVCS(k); + mf[kTa] = w[kTa] / tMoles; + } + } else { + // Perhaps, we could do a more sophisticated treatment below. + // But, will start with this. + for (size_t k = 0; k < nSpPhase; k++) { + size_t kTa = VolPhase->spGlobalIndexVCS(k); + mf[kTa]= 1.0 / (double) nSpPhase; + } + } + + VolPhase->setMolesFromVCS(VCS_STATECALC_OLD, &w[0]); + + // Now, calculate a sample naught Gibbs free energy calculation + // at the specified temperature. + for (size_t k = 0; k < nSpPhase; k++) { + vcs_SpeciesProperties* sProp = VolPhase->speciesProperty(k); + ts_ptr = sProp->SpeciesThermo; + ts_ptr->SS0_feSave = VolPhase->G0_calc_one(k)/ GasConstant; + ts_ptr->SS0_TSave = m_temperature; + } + } + + // Transfer initial element abundances to the vprob object. + // We have to find the mapping index from one to the other + gai.resize(m_nelem, 0.0); + set_gai(); + + // Printout the species information: PhaseID's and mole nums + if (m_printLvl > 1) { + writeline('=', 80, true, true); + writeline('=', 16, false); + plogf(" Cantera_to_vprob: START OF PROBLEM STATEMENT "); + writeline('=', 20); + writeline('=', 80); + plogf(" Phase IDs of species\n"); + plogf(" species phaseID phaseName "); + plogf(" Initial_Estimated_kMols\n"); + for (size_t i = 0; i < m_nsp; i++) { + size_t iphase = PhaseID[i]; + + vcs_VolPhase* VolPhase = VPhaseList[iphase]; + plogf("%16s %5d %16s", mphase->speciesName(i).c_str(), iphase, + VolPhase->PhaseName.c_str()); + if (SpeciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { + plogf(" Volts = %-10.5g\n", w[i]); + } else { + plogf(" %-10.5g\n", w[i]); + } + } + + // Printout of the Phase structure information + writeline('-', 80, true, true); + plogf(" Information about phases\n"); + plogf(" PhaseName PhaseNum SingSpec GasPhase EqnState NumSpec"); + plogf(" TMolesInert Tmoles(kmol)\n"); + + for (size_t iphase = 0; iphase < m_numPhases; iphase++) { + vcs_VolPhase* VolPhase = VPhaseList[iphase]; + plogf("%16s %5d %5d %8d %16s %8d %16e ", VolPhase->PhaseName.c_str(), + VolPhase->VP_ID_, VolPhase->m_singleSpecies, + VolPhase->m_gasPhase, VolPhase->eos_name(), + VolPhase->nSpecies(), VolPhase->totalMolesInert()); + plogf("%16e\n", VolPhase->totalMoles()); + } + + writeline('=', 80, true, true); + writeline('=', 16, false); + plogf(" Cantera_to_vprob: END OF PROBLEM STATEMENT "); + writeline('=', 20); + writeline('=', 80); + plogf("\n"); + } + + // Copy over the species molecular weights + m_wtSpecies = WtSpecies; + + // Copy over the charges + m_chargeSpecies = Charge; + + // Copy the VCS_SPECIES_THERMO structures + for (size_t kspec = 0; kspec < m_nsp; kspec++) { + delete m_speciesThermoList[kspec]; + VCS_SPECIES_THERMO* spf = SpeciesThermo[kspec]; + m_speciesThermoList[kspec] = spf->duplMyselfAsVCS_SPECIES_THERMO(); + if (m_speciesThermoList[kspec] == NULL) { + throw CanteraError("VCS_SOLVE::VCS_SOLVE", + " duplMyselfAsVCS_SPECIES_THERMO returned an error!"); + } + } + + // Copy the species unknown type + m_speciesUnknownType = SpeciesUnknownType; + + // w[] -> Copy the equilibrium mole number estimate if it exists. + if (w.size() != 0) { + m_molNumSpecies_old = w; + } else { + m_doEstimateEquil = -1; + m_molNumSpecies_old.assign(m_molNumSpecies_old.size(), 0.0); + } + + // Formulate the Goal Element Abundance Vector + if (gai.size() != 0) { + for (size_t i = 0; i < m_nelem; i++) { + m_elemAbundancesGoal[i] = gai[i]; + if (m_elType[i] == VCS_ELEM_TYPE_LATTICERATIO && m_elemAbundancesGoal[i] < 1.0E-10) { + m_elemAbundancesGoal[i] = 0.0; + } + } + } else { + if (m_doEstimateEquil == 0) { + double sum = 0; + for (size_t j = 0; j < m_nelem; j++) { + m_elemAbundancesGoal[j] = 0.0; + for (size_t kspec = 0; kspec < m_nsp; kspec++) { + if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { + sum += m_molNumSpecies_old[kspec]; + m_elemAbundancesGoal[j] += m_formulaMatrix(kspec,j) * m_molNumSpecies_old[kspec]; + } + } + if (m_elType[j] == VCS_ELEM_TYPE_LATTICERATIO && m_elemAbundancesGoal[j] < 1.0E-10 * sum) { + m_elemAbundancesGoal[j] = 0.0; + } + } + } else { + throw CanteraError("VCS_SOLVE::VCS_SOLVE", + "Element Abundances, m_elemAbundancesGoal[], not specified"); + } + } + + // zero out values that will be filled in later + // + // TPhMoles[] -> Untouched here. These will be filled in vcs_prep.c + // TPhMoles1[] + // DelTPhMoles[] + + // TPhInertMoles[] -> must be copied over here + for (size_t iph = 0; iph < m_numPhases; iph++) { + vcs_VolPhase* Vphase = VPhaseList[iph]; + TPhInertMoles[iph] = Vphase->totalMolesInert(); + } + + // m_speciesIndexVector[] is an index variable that keep track of solution + // vector rotations. + for (size_t i = 0; i < m_nsp; i++) { + m_speciesMapIndex[i] = i; + } + + // IndEl[] is an index variable that keep track of element vector rotations. + for (size_t i = 0; i < m_nelem; i++) { + m_elementMapIndex[i] = i; + } + + // Define all species to be major species, initially. + for (size_t i = 0; i < m_nsp; i++) { + m_speciesStatus[i] = VCS_SPECIES_MAJOR; + } + + // PhaseID: Fill in the species to phase mapping. Check for bad values at + // the same time. + if (PhaseID.size() != 0) { + std::vector numPhSp(m_numPhases, 0); + for (size_t kspec = 0; kspec < m_nsp; kspec++) { + size_t iph = PhaseID[kspec]; + if (iph >= m_numPhases) { + throw CanteraError("VCS_SOLVE::VCS_SOLVE", + "Species to Phase Mapping, PhaseID, has a bad value\n" + "\tPhaseID[{}] = {}\n" + "Allowed values: 0 to {}", kspec, iph, m_numPhases - 1); + } + m_phaseID[kspec] = PhaseID[kspec]; + m_speciesLocalPhaseIndex[kspec] = numPhSp[iph]; + numPhSp[iph]++; + } + for (size_t iph = 0; iph < m_numPhases; iph++) { + vcs_VolPhase* Vphase = VPhaseList[iph]; + if (numPhSp[iph] != Vphase->nSpecies()) { + throw CanteraError("VCS_SOLVE::VCS_SOLVE", + "Number of species in phase {}, {}, doesn't match ({} != {}) [vphase = {}]", + ser, iph, Vphase->PhaseName, numPhSp[iph], Vphase->nSpecies(), (size_t) Vphase); + } + } + } else { + if (m_numPhases == 1) { + for (size_t kspec = 0; kspec < m_nsp; kspec++) { + m_phaseID[kspec] = 0; + m_speciesLocalPhaseIndex[kspec] = kspec; + } + } else { + throw CanteraError("VCS_SOLVE::VCS_SOLVE", + "Species to Phase Mapping, PhaseID, is not defined"); + } + } + + for (size_t i = 0; i < m_nelem; i++) { + if (m_elType[i] == VCS_ELEM_TYPE_CHARGENEUTRALITY) { + if (m_elemAbundancesGoal[i] != 0.0) { + if (fabs(m_elemAbundancesGoal[i]) > 1.0E-9) { + throw CanteraError("VCS_SOLVE::vcs_prob_specifyFully", + "Charge neutrality condition {} is signicantly " + "nonzero, {}. Giving up", + m_elementName[i], m_elemAbundancesGoal[i]); + } else { + if (m_debug_print_lvl >= 2) { + plogf("Charge neutrality condition %s not zero, %g. Setting it zero\n", + m_elementName[i], m_elemAbundancesGoal[i]); + } + m_elemAbundancesGoal[i] = 0.0; + } + + } + } + } + + // Copy over the species names + for (size_t i = 0; i < m_nsp; i++) { + m_speciesName[i] = m_mix->speciesName(i); + } + + // FormulaMatrix[] -> Copy the formula matrix over + for (size_t i = 0; i < m_nsp; i++) { + bool nonzero = false; + for (size_t j = 0; j < m_nelem; j++) { + if (FormulaMatrix(i,j) != 0.0) { + nonzero = true; + } + m_formulaMatrix(i,j) = FormulaMatrix(i,j); + } + if (!nonzero) { + throw CanteraError("VCS_SOLVE::VCS_SOLVE", + "species {} {} has a zero formula matrix!", i, m_speciesName[i]); + } + } + + // Copy over all of the phase information. Use the object's assignment + // operator + for (size_t iph = 0; iph < m_numPhases; iph++) { + *m_VolPhaseList[iph] = *VPhaseList[iph]; + + // Fix up the species thermo pointer in the vcs_SpeciesThermo object. It + // should point to the species thermo pointer in the private data space. + vcs_VolPhase* Vphase = m_VolPhaseList[iph]; + for (size_t k = 0; k < Vphase->nSpecies(); k++) { + vcs_SpeciesProperties* sProp = Vphase->speciesProperty(k); + size_t kT = Vphase->spGlobalIndexVCS(k); + sProp->SpeciesThermo = m_speciesThermoList[kT]; + } + } + + // Specify the Activity Convention information + for (size_t iph = 0; iph < m_numPhases; iph++) { + vcs_VolPhase* Vphase = m_VolPhaseList[iph]; + 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 The activity for the solvent assumes + // that the it goes to one as the species mole fraction goes to one; + // i.e., it's really on a molarity framework. So + // SpecLnMnaught[iSolvent] = 0.0, and the loop below starts at 1, + // not 0. + size_t iSolvent = Vphase->spGlobalIndexVCS(0); + double mnaught = m_wtSpecies[iSolvent] / 1000.; + for (size_t k = 1; k < Vphase->nSpecies(); k++) { + size_t kspec = Vphase->spGlobalIndexVCS(k); + m_actConventionSpecies[kspec] = Vphase->p_activityConvention; + m_lnMnaughtSpecies[kspec] = log(mnaught); + } + } + } + } VCS_SOLVE::~VCS_SOLVE() @@ -202,89 +643,55 @@ void VCS_SOLVE::vcs_delete_memory() delete m_VCount; m_VCount = 0; - NSPECIES0 = 0; - NPHASE0 = 0; - m_numElemConstraints = 0; + m_nsp = 0; + m_nelem = 0; m_numComponents = 0; } -int VCS_SOLVE::vcs(int ifunc, int ipr, int ip1, int maxit) +int VCS_SOLVE::vcs(int ipr, int ip1, int maxit) { int retn = 0, iconv = 0; clockWC tickTock; int iprintTime = std::max(ipr, ip1); - //iprintTime = std::min(iprintTime, m_timing_print_lvl); - // writelog("iprintTime = {}; m_timing_print_lvl = {}; ipr = {}; ip1 = {}\n", - // iprintTime, m_timing_print_lvl, ipr, ip1); - if (ifunc > 2) { - plogf("vcs: Unrecognized value of ifunc, %d: bailing!\n", - ifunc); - return VCS_PUB_BAD; + // This function is called to copy the public data and the current + // problem specification into the current object's data structure. + retn = vcs_prob_specifyFully(); + if (retn != 0) { + plogf("vcs_pub_to_priv returned a bad status, %d: bailing!\n", + retn); + return retn; } - if (ifunc == 0) { - // This function is called to copy the public data and the current - // problem specification into the current object's data structure. - retn = vcs_prob_specifyFully(); - if (retn != 0) { - plogf("vcs_pub_to_priv returned a bad status, %d: bailing!\n", - retn); - return retn; - } - // Prep the problem data - // - adjust the identity of any phases - // - determine the number of components in the problem - retn = vcs_prep_oneTime(ip1); - if (retn != 0) { - plogf("vcs_prep_oneTime returned a bad status, %d: bailing!\n", - retn); - return retn; - } + prob_report(m_printLvl); + + // Prep the problem data + // - adjust the identity of any phases + // - determine the number of components in the problem + retn = vcs_prep(ip1); + if (retn != 0) { + plogf("vcs_prep_oneTime returned a bad status, %d: bailing!\n", + retn); + return retn; } - if (ifunc == 1) { - // This function is called to copy the current problem into the current - // object's data structure. - retn = vcs_prob_specify(); - if (retn != 0) { - plogf("vcs_prob_specify returned a bad status, %d: bailing!\n", - retn); - return retn; - } + + // Once we have defined the global internal data structure defining the + // problem, then we go ahead and solve the problem. + // + // (right now, all we do is solve fixed T, P problems. Methods for other + // problem types will go in at this level. For example, solving for + // fixed T, V problems will involve a 2x2 Newton's method, using loops + // over vcs_TP() to calculate the residual and Jacobian) + iconv = vcs_TP(ipr, ip1, maxit, m_temperature, m_pressurePA); + + // If requested to print anything out, go ahead and do so; + if (ipr > 0) { + vcs_report(iconv); } - if (ifunc != 2) { - // Prep the problem data for this particular instantiation of - // the problem - retn = vcs_prep(); - if (retn != VCS_SUCCESS) { - plogf("vcs_prep returned a bad status, %d: bailing!\n", retn); - return retn; - } - // Check to see if the current problem is well posed. - if (!vcs_wellPosed()) { - plogf("vcs has determined the problem is not well posed: Bailing\n"); - return VCS_PUB_BAD; - } - - // Once we have defined the global internal data structure defining the - // problem, then we go ahead and solve the problem. - // - // (right now, all we do is solve fixed T, P problems. Methods for other - // problem types will go in at this level. For example, solving for - // fixed T, V problems will involve a 2x2 Newton's method, using loops - // over vcs_TP() to calculate the residual and Jacobian) - iconv = vcs_TP(ipr, ip1, maxit, T, PresPA); - - // If requested to print anything out, go ahead and do so; - if (ipr > 0) { - vcs_report(iconv); - } - - vcs_prob_update(); - } + vcs_prob_update(); // Report on the time if requested to do so double te = tickTock.secondsWC(); @@ -305,382 +712,112 @@ int VCS_SOLVE::vcs(int ifunc, int ipr, int ip1, int maxit) int VCS_SOLVE::vcs_prob_specifyFully() { - const char* ser = - "vcs_pub_to_priv ERROR :ill defined interface -> bailout:\n\t"; + size_t kT = 0; + // Whether we have an estimate or not gets overwritten on + // the call to the equilibrium solver. + m_temperature = m_mix->temperature(); + m_pressurePA = m_mix->pressure(); + m_totalVol = m_mix->volume(); - // First Check to see whether we have room for the current problem size - if (NSPECIES0 < nspecies) { - plogf("%sPrivate Data is dimensioned too small\n", ser); - return VCS_PUB_BAD; + vector invSpecies(m_nsp); + for (size_t k = 0; k < m_nsp; k++) { + invSpecies[m_speciesMapIndex[k]] = k; } - size_t nph = NPhase; - if (NPHASE0 < nph) { - plogf("%sPrivate Data is dimensioned too small\n", ser); - return VCS_PUB_BAD; + + for (size_t iphase = 0; iphase < m_numPhases; iphase++) { + ThermoPhase* tPhase = &m_mix->phase(iphase); + vcs_VolPhase* volPhase = VPhaseList[iphase]; + + volPhase->setState_TP(m_temperature, m_pressurePA); + vector_fp muPhase(tPhase->nSpecies(),0.0); + tPhase->getChemPotentials(&muPhase[0]); + + // Loop through each species in the current phase + size_t nSpPhase = tPhase->nSpecies(); + for (size_t k = 0; k < nSpPhase; k++) { + // transfer chemical potential vector + m_gibbsSpecies[invSpecies[kT]] = muPhase[k]; + + kT++; + } + if ((nSpPhase == 1) && (volPhase->phiVarIndex() == 0)) { + volPhase->setExistence(VCS_PHASE_EXIST_ALWAYS); + } else if (volPhase->totalMoles() > 0.0) { + volPhase->setExistence(VCS_PHASE_EXIST_YES); + } else { + volPhase->setExistence(VCS_PHASE_EXIST_NO); + } + } + + // Printout the species information: PhaseID's and mole nums + if (m_printLvl > 1) { + writeline('=', 80, true, true); + writeline('=', 20, false); + plogf(" Cantera_to_vprob: START OF PROBLEM STATEMENT "); + writeline('=', 20); + writeline('=', 80); + plogf("\n"); + plogf(" Phase IDs of species\n"); + plogf(" species phaseID phaseName "); + plogf(" Initial_Estimated_kMols\n"); + for (size_t i = 0; i < m_nsp; i++) { + size_t iphase = m_phaseID[i]; + vcs_VolPhase* VolPhase = m_VolPhaseList[iphase]; + plogf("%16s %5d %16s", m_speciesName[i].c_str(), iphase, + VolPhase->PhaseName.c_str()); + if (SpeciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { + plogf(" Volts = %-10.5g\n", w[i]); + } else { + plogf(" %-10.5g\n", w[i]); + } + } + + // Printout of the Phase structure information + writeline('-', 80, true, true); + plogf(" Information about phases\n"); + plogf(" PhaseName PhaseNum SingSpec GasPhase EqnState NumSpec"); + plogf(" TMolesInert Tmoles(kmol)\n"); + + for (size_t iphase = 0; iphase < m_numPhases; iphase++) { + vcs_VolPhase* VolPhase = VPhaseList[iphase]; + plogf("%16s %5d %5d %8d %16s %8d %16e ", VolPhase->PhaseName.c_str(), + VolPhase->VP_ID_, VolPhase->m_singleSpecies, + VolPhase->m_gasPhase, VolPhase->eos_name(), + VolPhase->nSpecies(), VolPhase->totalMolesInert()); + plogf("%16e\n", VolPhase->totalMoles()); + } + + writeline('=', 80, true, true); + writeline('=', 20, false); + plogf(" Cantera_to_vprob: END OF PROBLEM STATEMENT "); + writeline('=', 20); + writeline('=', 80); + plogf("\n"); } // OK, We have room. Now, transfer the integer numbers - m_numSpeciesTot = nspecies; - m_numSpeciesRdc = m_numSpeciesTot; + m_numSpeciesRdc = m_nsp; // m_numRxnTot = number of noncomponents, also equal to the number of // reactions. Note, it's possible that the number of elements is greater // than the number of species. In that case set the number of reactions to // zero. - if (ne > nspecies) { + if (m_nelem > m_nsp) { m_numRxnTot = 0; } else { - m_numRxnTot = nspecies - ne; + m_numRxnTot = m_nsp - m_nelem; } m_numRxnRdc = m_numRxnTot; // number of minor species rxn -> all species rxn are major at the start. m_numRxnMinorZeroed = 0; - // NPhase = number of phases - m_numPhases = nph; - m_debug_print_lvl = vcs_debug_print_lvl; - // FormulaMatrix[] -> Copy the formula matrix over - for (size_t i = 0; i < nspecies; i++) { - bool nonzero = false; - for (size_t j = 0; j < ne; j++) { - if (FormulaMatrix(i,j) != 0.0) { - nonzero = true; - } - m_formulaMatrix(i,j) = FormulaMatrix(i,j); - } - if (!nonzero) { - plogf("vcs_prob_specifyFully:: species %d %s has a zero formula matrix!\n", i, - m_speciesName[i]); - return VCS_PUB_BAD; - } - } - - // Copy over the species molecular weights - m_wtSpecies = WtSpecies; - - // Copy over the charges - m_chargeSpecies = Charge; - - // Copy the VCS_SPECIES_THERMO structures - for (size_t kspec = 0; kspec < nspecies; kspec++) { - delete m_speciesThermoList[kspec]; - VCS_SPECIES_THERMO* spf = SpeciesThermo[kspec]; - m_speciesThermoList[kspec] = spf->duplMyselfAsVCS_SPECIES_THERMO(); - if (m_speciesThermoList[kspec] == NULL) { - plogf(" duplMyselfAsVCS_SPECIES_THERMO returned an error!\n"); - return VCS_PUB_BAD; - } - } - - // Copy the species unknown type - m_speciesUnknownType = SpeciesUnknownType; - - // iest => Do we have an initial estimate of the species mole numbers ? - m_doEstimateEquil = iest; - - // w[] -> Copy the equilibrium mole number estimate if it exists. - if (w.size() != 0) { - m_molNumSpecies_old = w; - } else { - m_doEstimateEquil = -1; - m_molNumSpecies_old.assign(m_molNumSpecies_old.size(), 0.0); - } - - // Formulate the Goal Element Abundance Vector - if (gai.size() != 0) { - for (size_t i = 0; i < ne; i++) { - m_elemAbundancesGoal[i] = gai[i]; - if (m_elType[i] == VCS_ELEM_TYPE_LATTICERATIO && m_elemAbundancesGoal[i] < 1.0E-10) { - m_elemAbundancesGoal[i] = 0.0; - } - } - } else { - if (m_doEstimateEquil == 0) { - double sum = 0; - for (size_t j = 0; j < ne; j++) { - m_elemAbundancesGoal[j] = 0.0; - 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(kspec,j) * m_molNumSpecies_old[kspec]; - } - } - if (m_elType[j] == VCS_ELEM_TYPE_LATTICERATIO && m_elemAbundancesGoal[j] < 1.0E-10 * sum) { - m_elemAbundancesGoal[j] = 0.0; - } - } - } else { - plogf("%sElement Abundances, m_elemAbundancesGoal[], not specified\n", ser); - return VCS_PUB_BAD; - } - } - - // zero out values that will be filled in later - // - // TPhMoles[] -> Untouched here. These will be filled in vcs_prep.c - // TPhMoles1[] - // DelTPhMoles[] - // - // T, Pres, copy over here - if (T > 0.0) { - m_temperature = T; - } else { - m_temperature = 293.15; - } - if (PresPA > 0.0) { - m_pressurePA = PresPA; - } else { - m_pressurePA = OneAtm; - } - - // TPhInertMoles[] -> must be copied over here - for (size_t iph = 0; iph < nph; iph++) { - vcs_VolPhase* Vphase = VPhaseList[iph]; - TPhInertMoles[iph] = Vphase->totalMolesInert(); - } - - // tolerance requirements -> copy them over here and later - m_tolmaj = tolmaj; - m_tolmin = tolmin; - m_tolmaj2 = 0.01 * m_tolmaj; - m_tolmin2 = 0.01 * m_tolmin; - - // m_speciesIndexVector[] is an index variable that keep track of solution - // vector rotations. - for (size_t i = 0; i < nspecies; i++) { - m_speciesMapIndex[i] = i; - } - - // IndEl[] is an index variable that keep track of element vector rotations. - for (size_t i = 0; i < ne; i++) { - m_elementMapIndex[i] = i; - } - - // Define all species to be major species, initially. - for (size_t i = 0; i < nspecies; i++) { - m_speciesStatus[i] = VCS_SPECIES_MAJOR; - } - - // PhaseID: Fill in the species to phase mapping. Check for bad values at - // the same time. - if (PhaseID.size() != 0) { - std::vector numPhSp(nph, 0); - for (size_t kspec = 0; kspec < nspecies; kspec++) { - size_t iph = PhaseID[kspec]; - if (iph >= nph) { - plogf("%sSpecies to Phase Mapping, PhaseID, has a bad value\n", - ser); - plogf("\tPhaseID[%d] = %d\n", kspec, iph); - plogf("\tAllowed values: 0 to %d\n", nph - 1); - return VCS_PUB_BAD; - } - m_phaseID[kspec] = PhaseID[kspec]; - m_speciesLocalPhaseIndex[kspec] = numPhSp[iph]; - numPhSp[iph]++; - } - for (size_t iph = 0; iph < nph; iph++) { - vcs_VolPhase* Vphase = VPhaseList[iph]; - if (numPhSp[iph] != Vphase->nSpecies()) { - plogf("%sNumber of species in phase %d, %s, doesn't match (%d != %d) [vphase = %d]\n", - ser, iph, Vphase->PhaseName, numPhSp[iph], Vphase->nSpecies(), (size_t) Vphase); - return VCS_PUB_BAD; - } - } - } else { - if (m_numPhases == 1) { - for (size_t kspec = 0; kspec < nspecies; kspec++) { - m_phaseID[kspec] = 0; - m_speciesLocalPhaseIndex[kspec] = kspec; - } - } else { - plogf("%sSpecies to Phase Mapping, PhaseID, is not defined\n", ser); - return VCS_PUB_BAD; - } - } - - for (size_t i = 0; i < ne; i++) { - if (m_elType[i] == VCS_ELEM_TYPE_CHARGENEUTRALITY) { - if (m_elemAbundancesGoal[i] != 0.0) { - if (fabs(m_elemAbundancesGoal[i]) > 1.0E-9) { - throw CanteraError("VCS_SOLVE::vcs_prob_specifyFully", - "Charge neutrality condition {} is signicantly " - "nonzero, {}. Giving up", - m_elementName[i], m_elemAbundancesGoal[i]); - } else { - if (m_debug_print_lvl >= 2) { - plogf("Charge neutrality condition %s not zero, %g. Setting it zero\n", - m_elementName[i], m_elemAbundancesGoal[i]); - } - m_elemAbundancesGoal[i] = 0.0; - } - - } - } - } - - // Copy over the species names - for (size_t i = 0; i < nspecies; i++) { - m_speciesName[i] = m_mix->speciesName(i); - } - - // Copy over all of the phase information. Use the object's assignment - // operator - for (size_t iph = 0; iph < nph; iph++) { - *m_VolPhaseList[iph] = *VPhaseList[iph]; - - // Fix up the species thermo pointer in the vcs_SpeciesThermo object. It - // should point to the species thermo pointer in the private data space. - vcs_VolPhase* Vphase = m_VolPhaseList[iph]; - for (size_t k = 0; k < Vphase->nSpecies(); k++) { - vcs_SpeciesProperties* sProp = Vphase->speciesProperty(k); - size_t kT = Vphase->spGlobalIndexVCS(k); - sProp->SpeciesThermo = m_speciesThermoList[kT]; - } - } - - // Specify the Activity Convention information - for (size_t iph = 0; iph < nph; iph++) { - vcs_VolPhase* Vphase = m_VolPhaseList[iph]; - 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 The activity for the solvent assumes - // that the it goes to one as the species mole fraction goes to one; - // i.e., it's really on a molarity framework. So - // SpecLnMnaught[iSolvent] = 0.0, and the loop below starts at 1, - // not 0. - size_t iSolvent = Vphase->spGlobalIndexVCS(0); - double mnaught = m_wtSpecies[iSolvent] / 1000.; - for (size_t k = 1; k < Vphase->nSpecies(); k++) { - size_t kspec = Vphase->spGlobalIndexVCS(k); - m_actConventionSpecies[kspec] = Vphase->p_activityConvention; - m_lnMnaughtSpecies[kspec] = log(mnaught); - } - } - } - - // Copy the title info - if (Title.size() == 0) { - m_title = "Unspecified Problem Title"; - } else { - m_title = Title; - } - - // Copy the volume info - m_totalVol = Vol; - if (m_PMVolumeSpecies.size() != 0) { - m_PMVolumeSpecies = VolPM; - } - // Return the success flag return VCS_SUCCESS; } -int VCS_SOLVE::vcs_prob_specify() -{ - string yo("vcs_prob_specify ERROR: "); - int retn = VCS_SUCCESS; - - m_temperature = T; - m_pressurePA = PresPA; - m_doEstimateEquil = iest; - m_totalVol = Vol; - m_tolmaj = tolmaj; - m_tolmin = tolmin; - m_tolmaj2 = 0.01 * m_tolmaj; - m_tolmin2 = 0.01 * m_tolmin; - - for (size_t kspec = 0; kspec < m_numSpeciesTot; ++kspec) { - size_t k = m_speciesMapIndex[kspec]; - m_molNumSpecies_old[kspec] = w[k]; - m_molNumSpecies_new[kspec] = mf[k]; - m_feSpecies_old[kspec] = m_gibbsSpecies[k]; - } - - // Transfer the element abundance goals to the solve object - for (size_t i = 0; i < m_numElemConstraints; i++) { - size_t j = m_elementMapIndex[i]; - m_elemAbundancesGoal[i] = gai[j]; - } - - // Try to do the best job at guessing at the title - if (Title.size() == 0) { - if (m_title.size() == 0) { - m_title = "Unspecified Problem Title"; - } - } else { - m_title = Title; - } - - // Copy over the phase information. For each entry in the phase structure, - // determine if that entry can change from its initial value Either copy - // over the new value or create an error condition. - bool status_change = false; - for (size_t iph = 0; iph < m_numPhases; iph++) { - vcs_VolPhase* vPhase = m_VolPhaseList[iph]; - vcs_VolPhase* pub_phase_ptr = VPhaseList[iph]; - - if (vPhase->VP_ID_ != pub_phase_ptr->VP_ID_) { - plogf("%sPhase numbers have changed:%d %d\n", - yo, vPhase->VP_ID_, pub_phase_ptr->VP_ID_); - retn = VCS_PUB_BAD; - } - - if (vPhase->m_singleSpecies != pub_phase_ptr->m_singleSpecies) { - plogf("%sSingleSpecies value have changed:%d %d\n", - yo, vPhase->m_singleSpecies, pub_phase_ptr->m_singleSpecies); - retn = VCS_PUB_BAD; - } - - if (vPhase->m_gasPhase != pub_phase_ptr->m_gasPhase) { - plogf("%sGasPhase value have changed:%d %d\n", - yo, vPhase->m_gasPhase, pub_phase_ptr->m_gasPhase); - retn = VCS_PUB_BAD; - } - - vPhase->m_eqnState = pub_phase_ptr->m_eqnState; - - if (vPhase->nSpecies() != pub_phase_ptr->nSpecies()) { - plogf("%sNVolSpecies value have changed:%d %d\n", - yo, vPhase->nSpecies(), pub_phase_ptr->nSpecies()); - retn = VCS_PUB_BAD; - } - - if (vPhase->PhaseName != pub_phase_ptr->PhaseName) { - plogf("%sPhaseName value have changed:%s %s\n", - yo, vPhase->PhaseName, pub_phase_ptr->PhaseName); - retn = VCS_PUB_BAD; - } - - if (vPhase->totalMolesInert() != pub_phase_ptr->totalMolesInert()) { - status_change = true; - } - - // Copy over the number of inert moles if it has changed. - TPhInertMoles[iph] = pub_phase_ptr->totalMolesInert(); - vPhase->setTotalMolesInert(pub_phase_ptr->totalMolesInert()); - if (TPhInertMoles[iph] > 0.0) { - vPhase->setExistence(2); - vPhase->m_singleSpecies = false; - } - - // Copy over the interfacial potential - double phi = pub_phase_ptr->electricPotential(); - vPhase->setElectricPotential(phi); - } - - if (status_change) { - vcs_SSPhase(); - } - - // Calculate the total number of moles in all phases. - vcs_tmoles(); - return retn; -} - int VCS_SOLVE::vcs_prob_update() { size_t k1 = 0; @@ -688,10 +825,10 @@ int VCS_SOLVE::vcs_prob_update() m_totalVol = vcs_VolTotal(m_temperature, m_pressurePA, &m_molNumSpecies_old[0], &m_PMVolumeSpecies[0]); - for (size_t i = 0; i < m_numSpeciesTot; ++i) { + for (size_t i = 0; i < m_nsp; ++i) { // Find the index of I in the index vector, m_speciesIndexVector[]. Call // it K1 and continue. - for (size_t j = 0; j < m_numSpeciesTot; ++j) { + for (size_t j = 0; j < m_nsp; ++j) { k1 = j; if (m_speciesMapIndex[j] == i) { break; @@ -705,14 +842,10 @@ int VCS_SOLVE::vcs_prob_update() w[i] = 0.0; } m_gibbsSpecies[i] = m_feSpecies_old[k1]; - VolPM[i] = m_PMVolumeSpecies[k1]; } - T = m_temperature; - PresPA = m_pressurePA; - Vol = m_totalVol; size_t kT = 0; - for (size_t iph = 0; iph < NPhase; iph++) { + for (size_t iph = 0; iph < m_numPhases; iph++) { vcs_VolPhase* pubPhase = VPhaseList[iph]; vcs_VolPhase* vPhase = m_VolPhaseList[iph]; pubPhase->setTotalMolesInert(vPhase->totalMolesInert()); @@ -751,9 +884,6 @@ int VCS_SOLVE::vcs_prob_update() sumMoles, pubPhase->totalMoles()); } } - - m_Iterations = m_VCount->Its; - m_NumBasisOptimizations = m_VCount->Basis_Opts; return VCS_SUCCESS; } diff --git a/src/equil/vcs_solve_TP.cpp b/src/equil/vcs_solve_TP.cpp index fcf539673..80290a223 100644 --- a/src/equil/vcs_solve_TP.cpp +++ b/src/equil/vcs_solve_TP.cpp @@ -66,11 +66,11 @@ int VCS_SOLVE::vcs_solve_TP(int print_lvl, int printDetails, int maxit) clockWC ticktock; // temporary space for usage in this routine and in subroutines - m_sm.assign(m_numElemConstraints*m_numElemConstraints, 0.0); - m_ss.assign(m_numElemConstraints, 0.0); - m_sa.assign(m_numElemConstraints, 0.0); - m_aw.assign(m_numSpeciesTot, 0.0); - m_wx.assign(m_numElemConstraints, 0.0); + m_sm.assign(m_nelem * m_nelem, 0.0); + m_ss.assign(m_nelem, 0.0); + m_sa.assign(m_nelem, 0.0); + m_aw.assign(m_nsp, 0.0); + m_wx.assign(m_nelem, 0.0); int solveFail = false; @@ -80,8 +80,8 @@ int VCS_SOLVE::vcs_solve_TP(int print_lvl, int printDetails, int maxit) // Printout the initial conditions for problem if (print_lvl != 0) { plogf("VCS CALCULATION METHOD\n\n "); - plogf("%s\n", m_title); - plogf("\n\n%5d SPECIES\n%5d ELEMENTS\n", m_numSpeciesTot, m_numElemConstraints); + plogf("MultiPhase Object\n"); + plogf("\n\n%5d SPECIES\n%5d ELEMENTS\n", m_nsp, m_nelem); plogf("%5d COMPONENTS\n", m_numComponents); plogf("%5d PHASES\n", m_numPhases); plogf(" PRESSURE%22.8g %3s\n", m_pressurePA, "Pa "); @@ -95,7 +95,7 @@ int VCS_SOLVE::vcs_solve_TP(int print_lvl, int printDetails, int maxit) } plogf("\n ELEMENTAL ABUNDANCES CORRECT"); plogf(" FROM ESTIMATE Type\n\n"); - for (size_t i = 0; i < m_numElemConstraints; ++i) { + for (size_t i = 0; i < m_nelem; ++i) { writeline(' ', 26, false); plogf("%-2.2s", m_elementName[i]); plogf("%20.12E%20.12E %3d\n", m_elemAbundancesGoal[i], m_elemAbundances[i], @@ -114,18 +114,18 @@ int VCS_SOLVE::vcs_solve_TP(int print_lvl, int printDetails, int maxit) writeline(' ', 41, false); plogf(" STAN_CHEM_POT EQUILIBRIUM_EST. Species_Type\n\n"); writeline(' ', 20, false); - for (size_t i = 0; i < m_numElemConstraints; ++i) { + for (size_t i = 0; i < m_nelem; ++i) { plogf("%-4.4s ", m_elementName[i]); } plogf(" PhaseID\n"); double RT = GasConstant * m_temperature; - for (size_t i = 0; i < m_numSpeciesTot; ++i) { + for (size_t i = 0; i < m_nsp; ++i) { plogf(" %-18.18s", m_speciesName[i]); - for (size_t j = 0; j < m_numElemConstraints; ++j) { + for (size_t j = 0; j < m_nelem; ++j) { plogf("% -7.3g ", m_formulaMatrix(i,j)); } plogf(" %3d ", m_phaseID[i]); - writeline(' ', std::max(55-int(m_numElemConstraints)*8, 0), false); + writeline(' ', std::max(55-int(m_nelem)*8, 0), false); plogf("%12.5E %12.5E", RT * m_SSfeSpecies[i], m_molNumSpecies_old[i]); if (m_speciesUnknownType[i] == VCS_SPECIES_TYPE_MOLNUM) { plogf(" Mol_Num\n"); @@ -137,7 +137,7 @@ int VCS_SOLVE::vcs_solve_TP(int print_lvl, int printDetails, int maxit) } } - for (size_t i = 0; i < m_numSpeciesTot; ++i) { + for (size_t i = 0; i < m_nsp; ++i) { if (m_molNumSpecies_old[i] < 0.0) { plogf("On Input species %-12s has a negative MF, setting it small\n", m_speciesName[i]); @@ -250,7 +250,7 @@ int VCS_SOLVE::vcs_solve_TP(int print_lvl, int printDetails, int maxit) // Evaluate the final mole fractions storing them in wt[] m_molNumSpecies_new.assign(m_molNumSpecies_new.size(), 0.0); - for (size_t kspec = 0; kspec < m_numSpeciesTot; ++kspec) { + for (size_t kspec = 0; kspec < m_nsp; ++kspec) { if (m_SSPhase[kspec]) { m_molNumSpecies_new[kspec] = 1.0; } else { @@ -424,13 +424,13 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, // m_deltaMolNumSpecies[kspec] Delta in the Species Mole Numbers if (iphaseDelete != npos) { debuglog(" --- Main Loop Treatment -> Circumvented due to Phase Deletion\n", m_debug_print_lvl >= 2); - for (size_t k = 0; k < m_numSpeciesTot; k++) { + for (size_t k = 0; k < m_nsp; k++) { m_molNumSpecies_new[k] = m_molNumSpecies_old[k] + m_deltaMolNumSpecies[k]; size_t iph = m_phaseID[k]; m_tPhaseMoles_new[iph] += m_deltaMolNumSpecies[k]; } if (kspec >= m_numComponents) { - if (m_molNumSpecies_new[m_numSpeciesTot] != 0.0) { + if (m_molNumSpecies_new[m_nsp] != 0.0) { throw CanteraError("VCS_SOLVE::solve_tp_inner", "we shouldn't be here!"); } @@ -451,7 +451,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, // Calculate the new chemical potentials using the tentative solution // values. We only calculate a subset of these, because we have only // updated a subset of the W(). - vcs_dfe(VCS_STATECALC_NEW, 0, 0, m_numSpeciesTot); + vcs_dfe(VCS_STATECALC_NEW, 0, 0, m_nsp); // Evaluate DeltaG for all components if ITI=0, and for major components // only if ITI NE 0 @@ -515,7 +515,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, } } } else { - for (size_t j = 0; j < m_numElemConstraints; ++j) { + for (size_t j = 0; j < m_nelem; ++j) { int elType = m_elType[j]; if (elType == VCS_ELEM_TYPE_ABSPOS) { double atomComp = m_formulaMatrix(kspec,j); @@ -734,7 +734,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, ++m_numRxnMinorZeroed; allMinorZeroedSpecies = (m_numRxnMinorZeroed == m_numRxnRdc); - for (size_t kk = 0; kk < m_numSpeciesTot; kk++) { + for (size_t kk = 0; kk < m_nsp; kk++) { m_deltaMolNumSpecies[kk] = 0.0; m_molNumSpecies_new[kk] = m_molNumSpecies_old[kk]; } @@ -837,7 +837,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, plogf(" --- Reduction in step size due to component %s going negative = %11.3E\n", m_speciesName[ll], par); } - for (size_t i = 0; i < m_numSpeciesTot; ++i) { + for (size_t i = 0; i < m_nsp; ++i) { m_deltaMolNumSpecies[i] *= par; } for (size_t iph = 0; iph < m_numPhases; iph++) { @@ -847,14 +847,14 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, par = 1.0; } checkDelta1(&m_deltaMolNumSpecies[0], - &m_deltaPhaseMoles[0], m_numSpeciesTot); + &m_deltaPhaseMoles[0], m_nsp); // Now adjust the wt[kspec]'s so that the reflect the decrease in the // overall length of m_deltaMolNumSpecies[kspec] just calculated. At the // end of this section wt[], m_deltaMolNumSpecies[], tPhMoles, and // tPhMoles1 should all be consistent with a new estimate of the state // of the system. - for (size_t kspec = 0; kspec < m_numSpeciesTot; ++kspec) { + for (size_t kspec = 0; kspec < m_nsp; ++kspec) { m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + m_deltaMolNumSpecies[kspec]; if (m_molNumSpecies_new[kspec] < 0.0 && (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE)) { @@ -876,7 +876,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, // Calculate the new chemical potentials using the tentative solution // values. We only calculate a subset of these, because we have only // updated a subset of the W(). - vcs_dfe(VCS_STATECALC_NEW, 0, 0, m_numSpeciesTot); + vcs_dfe(VCS_STATECALC_NEW, 0, 0, m_nsp); // Evaluate DeltaG for all components if ITI=0, and for major components // only if ITI NE 0 @@ -962,7 +962,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, m_molNumSpecies_new[i], m_feSpecies_old[i], m_feSpecies_new[i], m_deltaGRxn_old[l1], m_deltaGRxn_new[l1]); } - for (size_t kspec = m_numSpeciesRdc; kspec < m_numSpeciesTot; ++kspec) { + for (size_t kspec = m_numSpeciesRdc; kspec < m_nsp; ++kspec) { size_t l1 = kspec - m_numComponents; plogf(" --- %-12.12s", m_speciesName[kspec]); plogf(" %2d %14.6E%14.6E%14.6E%14.6E%14.6E%14.6E\n", @@ -1280,7 +1280,7 @@ void VCS_SOLVE::solve_tp_equilib_check(bool& allMinorZeroedSpecies, // If we have deleted a species then we need to recheck the the deleted // species, before exiting - if (m_numSpeciesRdc != m_numSpeciesTot) { + if (m_numSpeciesRdc != m_nsp) { stage = RECHECK_DELETED; return; } @@ -1600,7 +1600,7 @@ void VCS_SOLVE::vcs_reinsert_deleted(size_t kspec) if (! m_SSPhase[kspec]) { if (Vphase->exists() == VCS_PHASE_EXIST_NO) { Vphase->setExistence(VCS_PHASE_EXIST_YES); - for (size_t k = 0; k < m_numSpeciesTot; k++) { + for (size_t k = 0; k < m_nsp; k++) { if (m_phaseID[k] == iph && m_speciesStatus[k] != VCS_SPECIES_DELETED) { m_speciesStatus[k] = VCS_SPECIES_MINOR; } @@ -1723,7 +1723,7 @@ bool VCS_SOLVE::vcs_delete_multiphase(const size_t iph) // reinstate the "major ones" in the future. Note, species in phases with // zero mole numbers are still considered active. Whether the phase pops // back into existence or not is checked as part of the main iteration loop. - for (size_t kspec = m_numSpeciesRdc; kspec < m_numSpeciesTot; ++kspec) { + for (size_t kspec = m_numSpeciesRdc; kspec < m_nsp; ++kspec) { if (m_phaseID[kspec] == iph) { m_molNumSpecies_old[kspec] = 0.0; m_molNumSpecies_new[kspec] = 0.0; @@ -1760,14 +1760,14 @@ int VCS_SOLVE::vcs_recheck_deleted() if (m_debug_print_lvl >= 2) { plogf(" --- Start rechecking deleted species in multispec phases\n"); } - if (m_numSpeciesRdc == m_numSpeciesTot) { + if (m_numSpeciesRdc == m_nsp) { return 0; } // Use the standard chemical potentials for the chemical potentials of // deleted species. Then, calculate Delta G for for formation reactions. // Note: fe[] here includes everything except for the ln(x[i]) term - for (size_t kspec = m_numSpeciesRdc; kspec < m_numSpeciesTot; ++kspec) { + for (size_t kspec = m_numSpeciesRdc; kspec < m_nsp; ++kspec) { size_t iph = m_phaseID[kspec]; m_feSpecies_new[kspec] = (m_SSfeSpecies[kspec] + log(m_actCoeffSpecies_old[kspec]) - m_lnMnaughtSpecies[kspec] @@ -1830,7 +1830,7 @@ int VCS_SOLVE::vcs_recheck_deleted() size_t VCS_SOLVE::vcs_add_all_deleted() { - if (m_numSpeciesRdc == m_numSpeciesTot) { + if (m_numSpeciesRdc == m_nsp) { return 0; } @@ -1841,7 +1841,7 @@ size_t VCS_SOLVE::vcs_add_all_deleted() // of the routine. m_molNumSpecies_new = m_molNumSpecies_old; for (int cits = 0; cits < 3; cits++) { - for (size_t kspec = m_numSpeciesRdc; kspec < m_numSpeciesTot; ++kspec) { + for (size_t kspec = m_numSpeciesRdc; kspec < m_nsp; ++kspec) { size_t iph = m_phaseID[kspec]; vcs_VolPhase* Vphase = m_VolPhaseList[iph]; if (m_molNumSpecies_new[kspec] == 0.0) { @@ -1898,7 +1898,7 @@ size_t VCS_SOLVE::vcs_add_all_deleted() } } vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD); - vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_numSpeciesTot); + vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_nsp); vcs_deltag(0, true, VCS_STATECALC_OLD); size_t retn = 0; @@ -2050,18 +2050,18 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[], plogf("\n"); plogf(" --- Formula Matrix used in BASOPT calculation\n"); plogf(" --- Active | "); - for (size_t j = 0; j < m_numElemConstraints; j++) { + for (size_t j = 0; j < m_nelem; j++) { plogf(" %1d ", m_elementActive[j]); } plogf("\n"); plogf(" --- Species | "); - for (size_t j = 0; j < m_numElemConstraints; j++) { + for (size_t j = 0; j < m_nelem; j++) { writelog(" {:>8.8s}", m_elementName[j]); } plogf("\n"); - for (k = 0; k < m_numSpeciesTot; k++) { + for (k = 0; k < m_nsp; k++) { writelog(" --- {:>11.11s} | ", m_speciesName[k]); - for (size_t j = 0; j < m_numElemConstraints; j++) { + for (size_t j = 0; j < m_nelem; j++) { plogf(" %8.2g", m_formulaMatrix(k,j)); } plogf("\n"); @@ -2073,17 +2073,17 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[], // Calculate the maximum value of the number of components possible. It's // equal to the minimum of the number of elements and the number of total // species. - size_t ncTrial = std::min(m_numElemConstraints, m_numSpeciesTot); + size_t ncTrial = std::min(m_nelem, m_nsp); m_numComponents = ncTrial; *usedZeroedSpecies = false; vector_int ipiv(ncTrial); // Use a temporary work array for the mole numbers, aw[] std::copy(m_molNumSpecies_old.begin(), - m_molNumSpecies_old.begin() + m_numSpeciesTot, aw); + m_molNumSpecies_old.begin() + m_nsp, aw); // Take out the Voltage unknowns from consideration - for (k = 0; k < m_numSpeciesTot; k++) { + for (k = 0; k < m_nsp; k++) { if (m_speciesUnknownType[k] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { aw[k] = test; } @@ -2101,7 +2101,7 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[], // the largest remaining species. Return its identity in K. The // first search criteria is always the largest positive magnitude of // the mole number. - k = vcs_basisOptMax(aw, jr, m_numSpeciesTot); + k = vcs_basisOptMax(aw, jr, m_nsp); // The fun really starts when you have run out of species that have // a significant concentration. It becomes extremely important to @@ -2147,11 +2147,11 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[], size_t kfound = npos; int minNonZeroes = 100000; int nonZeroesKspec = 0; - for (size_t kspec = ncTrial; kspec < m_numSpeciesTot; kspec++) { + for (size_t kspec = ncTrial; kspec < m_nsp; kspec++) { if (aw[kspec] >= 0.0 && m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { maxConcPossKspec = 1.0E10; nonZeroesKspec = 0; - for (size_t j = 0; j < m_numElemConstraints; ++j) { + for (size_t j = 0; j < m_nelem; ++j) { if (m_elementActive[j] && m_elType[j] == VCS_ELEM_TYPE_ABSPOS) { double nu = m_formulaMatrix(kspec,j); if (nu != 0.0) { @@ -2180,7 +2180,7 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[], if (kfound == npos) { double gmin = 0.0; kfound = k; - for (size_t kspec = ncTrial; kspec < m_numSpeciesTot; kspec++) { + for (size_t kspec = ncTrial; kspec < m_nsp; kspec++) { if (aw[kspec] >= 0.0) { size_t irxn = kspec - ncTrial; if (m_deltaGRxn_new[irxn] < gmin) { @@ -2197,18 +2197,18 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[], m_numComponents = jr; ncTrial = m_numComponents; size_t numPreDeleted = m_numRxnTot - m_numRxnRdc; - if (numPreDeleted != (m_numSpeciesTot - m_numSpeciesRdc)) { + if (numPreDeleted != (m_nsp - m_numSpeciesRdc)) { throw CanteraError("VCS_SOLVE::vcs_basopt", "we shouldn't be here"); } - m_numRxnTot = m_numSpeciesTot - ncTrial; + m_numRxnTot = m_nsp - ncTrial; m_numRxnRdc = m_numRxnTot - numPreDeleted; - m_numSpeciesRdc = m_numSpeciesTot - numPreDeleted; - for (size_t i = 0; i < m_numSpeciesTot; ++i) { + m_numSpeciesRdc = m_nsp - numPreDeleted; + for (size_t i = 0; i < m_nsp; ++i) { m_indexRxnToSpecies[i] = ncTrial + i; } if (m_debug_print_lvl >= 2) { plogf(" --- Total number of components found = %3d (ne = %d)\n ", - ncTrial, m_numElemConstraints); + ncTrial, m_nelem); } goto L_END_LOOP; } @@ -2222,8 +2222,8 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[], // Modified Gram-Schmidt Method, p. 202 Dalquist // QR factorization of a matrix without row pivoting. size_t jl = jr; - for (size_t j = 0; j < m_numElemConstraints; ++j) { - sm[j + jr*m_numElemConstraints] = m_formulaMatrix(k,j); + for (size_t j = 0; j < m_nelem; ++j) { + sm[j + jr*m_nelem] = m_formulaMatrix(k,j); } if (jl > 0) { // Compute the coefficients of JA column of the the upper @@ -2231,16 +2231,16 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[], // different than Dalquist) R_JA_JA = 1 for (size_t j = 0; j < jl; ++j) { ss[j] = 0.0; - for (size_t i = 0; i < m_numElemConstraints; ++i) { - ss[j] += sm[i + jr*m_numElemConstraints] * sm[i + j*m_numElemConstraints]; + for (size_t i = 0; i < m_nelem; ++i) { + ss[j] += sm[i + jr*m_nelem] * sm[i + j*m_nelem]; } ss[j] /= sa[j]; } // Now make the new column, (*,JR), orthogonal to the previous // columns for (size_t j = 0; j < jl; ++j) { - for (size_t i = 0; i < m_numElemConstraints; ++i) { - sm[i + jr*m_numElemConstraints] -= ss[j] * sm[i + j*m_numElemConstraints]; + for (size_t i = 0; i < m_nelem; ++i) { + sm[i + jr*m_nelem] -= ss[j] * sm[i + j*m_nelem]; } } } @@ -2248,8 +2248,8 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[], // Find the new length of the new column in Q. It will be used in // the denominator in future row calcs. sa[jr] = 0.0; - for (size_t ml = 0; ml < m_numElemConstraints; ++ml) { - sa[jr] += pow(sm[ml + jr*m_numElemConstraints], 2); + for (size_t ml = 0; ml < m_nelem; ++ml) { + sa[jr] += pow(sm[ml + jr*m_nelem], 2); } // IF NORM OF NEW ROW .LT. 1E-3 REJECT @@ -2337,23 +2337,23 @@ L_END_LOOP: } // Solve the linear system to calculate the reaction matrix, // m_stoichCoeffRxnMatrix. - solve(C, m_stoichCoeffRxnMatrix.ptrColumn(0), m_numRxnTot, m_numElemConstraints); + solve(C, m_stoichCoeffRxnMatrix.ptrColumn(0), m_numRxnTot, m_nelem); // NOW, if we have interfacial voltage unknowns, what we did was just wrong // -> hopefully it didn't blow up. Redo the problem. Search for inactive E juse = npos; jlose = npos; - for (size_t j = 0; j < m_numElemConstraints; j++) { + for (size_t j = 0; j < m_nelem; j++) { if (!m_elementActive[j] && !strcmp(m_elementName[j].c_str(), "E")) { juse = j; } } - for (size_t j = 0; j < m_numElemConstraints; j++) { + for (size_t j = 0; j < m_nelem; j++) { if (m_elementActive[j] && !strncmp((m_elementName[j]).c_str(), "cn_", 3)) { jlose = j; } } - for (k = 0; k < m_numSpeciesTot; k++) { + for (k = 0; k < m_nsp; k++) { if (m_speciesUnknownType[k] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { for (size_t j = 0; j < ncTrial; ++j) { for (size_t i = 0; i < ncTrial; ++i) { @@ -2375,7 +2375,7 @@ L_END_LOOP: } } - solve(C, aw, 1, m_numElemConstraints); + solve(C, aw, 1, m_nelem); size_t i = k - ncTrial; for (size_t j = 0; j < ncTrial; j++) { m_stoichCoeffRxnMatrix(j,i) = aw[j]; @@ -2580,7 +2580,7 @@ int VCS_SOLVE::vcs_species_type(const size_t kspec) const // see if the species has an element which is so low that species will // always be zero - for (size_t j = 0; j < m_numElemConstraints; ++j) { + for (size_t j = 0; j < m_nelem; ++j) { if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) { double atomComp = m_formulaMatrix(kspec,j); if (atomComp > 0.0) { @@ -2767,7 +2767,7 @@ void VCS_SOLVE::vcs_dfe(const int stateCalc, tlogMoles[iph] = tPhInertMoles[iph]; } - for (size_t kspec = 0; kspec < m_numSpeciesTot; kspec++) { + for (size_t kspec = 0; kspec < m_nsp; kspec++) { if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { size_t iph = m_phaseID[kspec]; tlogMoles[iph] += molNum[kspec]; @@ -2964,7 +2964,7 @@ double VCS_SOLVE::vcs_tmoles() for (size_t i = 0; i < m_numPhases; i++) { m_tPhaseMoles_old[i] = TPhInertMoles[i]; } - for (size_t i = 0; i < m_numSpeciesTot; i++) { + for (size_t i = 0; i < m_nsp; i++) { if (m_speciesUnknownType[i] == VCS_SPECIES_TYPE_MOLNUM) { m_tPhaseMoles_old[m_phaseID[i]] += m_molNumSpecies_old[i]; } @@ -2989,7 +2989,7 @@ void VCS_SOLVE::check_tmoles() const for (size_t i = 0; i < m_numPhases; i++) { double m_tPhaseMoles_old_a = TPhInertMoles[i]; - for (size_t k = 0; k < m_numSpeciesTot; k++) { + for (size_t k = 0; k < m_nsp; k++) { if (m_speciesUnknownType[k] == VCS_SPECIES_TYPE_MOLNUM && m_phaseID[k] == i) { m_tPhaseMoles_old_a += m_molNumSpecies_old[k]; } @@ -3031,7 +3031,7 @@ bool VCS_SOLVE::vcs_evaluate_speciesType() } else if (m_debug_print_lvl >= 5) { plogf(" --- Species Status decision is reevaluated\n"); } - for (size_t kspec = 0; kspec < m_numSpeciesTot; ++kspec) { + for (size_t kspec = 0; kspec < m_nsp; ++kspec) { m_speciesStatus[kspec] = vcs_species_type(kspec); if (m_debug_print_lvl >= 5) { plogf(" --- %-16s: ", m_speciesName[kspec]); @@ -3326,7 +3326,7 @@ void VCS_SOLVE::vcs_printDeltaG(const int stateCalc) writelog(" "); writeline('-', 132); - for (size_t kspec = 0; kspec < m_numSpeciesTot; kspec++) { + for (size_t kspec = 0; kspec < m_nsp; kspec++) { size_t irxn = npos; if (kspec >= m_numComponents) { irxn = kspec - m_numComponents; @@ -3402,7 +3402,7 @@ void VCS_SOLVE::vcs_switch_pos(const bool ifunc, const size_t k1, const size_t k if (k1 == k2) { return; } - if (k1 >= m_numSpeciesTot || k2 >= m_numSpeciesTot) { + if (k1 >= m_nsp || k2 >= m_nsp) { plogf("vcs_switch_pos: ifunc = 0: inappropriate args: %d %d\n", k1, k2); } @@ -3440,14 +3440,14 @@ void VCS_SOLVE::vcs_switch_pos(const bool ifunc, const size_t k1, const size_t k std::swap(m_speciesThermoList[k1], m_speciesThermoList[k2]); std::swap(m_PMVolumeSpecies[k1], m_PMVolumeSpecies[k2]); - for (size_t j = 0; j < m_numElemConstraints; ++j) { + for (size_t j = 0; j < m_nelem; ++j) { std::swap(m_formulaMatrix(k1,j), m_formulaMatrix(k2,j)); } if (m_useActCoeffJac && k1 != k2) { - for (size_t i = 0; i < m_numSpeciesTot; i++) { + for (size_t i = 0; i < m_nsp; i++) { std::swap(m_np_dLnActCoeffdMolNum(k1,i), m_np_dLnActCoeffdMolNum(k2,i)); } - for (size_t i = 0; i < m_numSpeciesTot; i++) { + for (size_t i = 0; i < m_nsp; i++) { std::swap(m_np_dLnActCoeffdMolNum(i,k1), m_np_dLnActCoeffdMolNum(i,k2)); } }