From 511a39fc7393b7b1f9401d9b7b336d517c577948 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Sun, 3 Apr 2016 21:04:09 -0400 Subject: [PATCH] [Equil] Standardize VCS solver to always use MKS units --- include/cantera/equil/vcs_VolPhase.h | 36 +----- include/cantera/equil/vcs_defs.h | 35 ----- include/cantera/equil/vcs_internal.h | 11 -- include/cantera/equil/vcs_prob.h | 29 +---- include/cantera/equil/vcs_solve.h | 57 +-------- include/cantera/equil/vcs_species_thermo.h | 10 +- src/equil/vcs_MultiPhaseEquil.cpp | 46 +------ src/equil/vcs_TP.cpp | 12 -- src/equil/vcs_VolPhase.cpp | 4 - src/equil/vcs_nondim.cpp | 142 ++++++--------------- src/equil/vcs_prob.cpp | 22 +--- src/equil/vcs_report.cpp | 5 +- src/equil/vcs_solve.cpp | 9 +- src/equil/vcs_solve_TP.cpp | 18 +-- src/equil/vcs_species_thermo.cpp | 25 +--- src/equil/vcs_util.cpp | 19 --- 16 files changed, 63 insertions(+), 417 deletions(-) diff --git a/include/cantera/equil/vcs_VolPhase.h b/include/cantera/equil/vcs_VolPhase.h index 355477362..775d0f1bc 100644 --- a/include/cantera/equil/vcs_VolPhase.h +++ b/include/cantera/equil/vcs_VolPhase.h @@ -292,9 +292,6 @@ public: const ThermoPhase* ptrThermoPhase() const; //! Return the total moles in the phase - /*! - * Units -> depends on VCS_UnitsFormat variable. Cantera -> J/kmol - */ double totalMoles() const; //! Returns the mole fraction of the kspec species @@ -586,31 +583,6 @@ public: */ size_t ChargeNeutralityElement; - //! Units for the chemical potential data, pressure data, volume, - //! and species amounts - /*! - * All internally stored quantities will have these units. Also, printed - * quantities will display in these units. Input quantities are expected in - * these units. - * - * | | | Chem_Pot | Pres | vol | moles| - * |---|--------------------|-------------------------|------|------|------| - * |-1 | VCS_UNITS_KCALMOL | kcal/gmol | Pa | m**3 | kmol | - * | 0 | VCS_UNITS_UNITLESS | MU / RT -> no units | Pa | m**3 | kmol | - * | 1 | VCS_UNITS_KJMOL | kJ / gmol | Pa | m**3 | kmol | - * | 2 | VCS_UNITS_KELVIN | KELVIN -> MU / R | Pa | m**3 | kmol | - * | 3 | VCS_UNITS_MKS | Joules / Kmol (Cantera) | Pa | m**3 | kmol | - * - * see vcs_defs.h for more information. - * - * Currently, this value should be the same as the owning VCS_PROB or - * VCS_SOLVE object. There is no code for handling anything else atm. - * - * (This variable is needed for the vcsc code, where it is not equal to - * VCS_UNITS_MKS). - */ - int p_VCS_UnitsFormat; - //! Convention for the activity formulation /*! * * 0 = molar based activities (default) @@ -773,9 +745,7 @@ private: /*! * Note, This is the chemical potential derived strictly from the polynomial * in temperature. Pressure effects have to be added in to get to the - * standard state. - * - * Units -> depends on VCS_UnitsFormat variable. Cantera -> J/kmol + * standard state. Units are J/kmol. */ mutable vector_fp SS0ChemicalPotential; @@ -783,9 +753,7 @@ private: //! current Temperature and pressure. /*! * Note, This is the chemical potential at unit activity. Thus, we can call - * it the standard state chemical potential as well. - * - * Units -> depends on VCS_UnitsFormat variable. Cantera -> J/kmol. + * it the standard state chemical potential as well. Units are J/kmol. */ mutable vector_fp StarChemicalPotential; diff --git a/include/cantera/equil/vcs_defs.h b/include/cantera/equil/vcs_defs.h index 0768eb0c8..e70208430 100644 --- a/include/cantera/equil/vcs_defs.h +++ b/include/cantera/equil/vcs_defs.h @@ -227,41 +227,6 @@ namespace Cantera //@} -/*! - * @name Units for the chemical potential data and pressure variables - * - * @verbatim - Chem_Pot Pres vol moles - ------------------------------------------------- - VCS_UNITS_KCALMOL = kcal/mol Pa m**3 kmol - VCS_UNITS_UNITLESS = MU / RT -> no units Pa m**3 kmol - VCS_UNITS_KJMOL = kJ / mol Pa m**3 kmol - VCS_UNITS_KELVIN = KELVIN -> MU / R Pa m**3 kmol - VCS_UNITS_MKS = Joules / Kmol (Cantera) Pa m**3 kmol - - Energy: - VCS_UNITS_KCALMOL = kcal/mol - VCS_UNITS_UNITLESS = MU / RT -> no units - VCS_UNITS_KJMOL = kJ / mol - VCS_UNITS_KELVIN = KELVIN -> MU / R - VCS_UNITS_MKS = J / kmol - - Pressure: (Pref and Pres) - VCS_UNITS_KCALMOL = Pa - VCS_UNITS_UNITLESS = Pa - VCS_UNITS_KJMOL = Pa - VCS_UNITS_KELVIN = Pa - VCS_UNITS_MKS = Pa = kg / m s2 - @endverbatim - * @{ - */ -#define VCS_UNITS_KCALMOL -1 -#define VCS_UNITS_UNITLESS 0 -#define VCS_UNITS_KJMOL 1 -#define VCS_UNITS_KELVIN 2 -#define VCS_UNITS_MKS 3 -//@} - /*! * @name Types of Element Constraint Equations * diff --git a/include/cantera/equil/vcs_internal.h b/include/cantera/equil/vcs_internal.h index b3f649336..c7d30be0a 100644 --- a/include/cantera/equil/vcs_internal.h +++ b/include/cantera/equil/vcs_internal.h @@ -86,17 +86,6 @@ public: double T_Time_vcs; }; -//! Returns the value of the gas constant in the units specified by parameter -/*! - * @param mu_units Specifies the units. - * - VCS_UNITS_KCALMOL: kcal gmol-1 K-1 - * - VCS_UNITS_UNITLESS: 1.0 K-1 - * - VCS_UNITS_KJMOL: kJ gmol-1 K-1 - * - VCS_UNITS_KELVIN: 1.0 K-1 - * - VCS_UNITS_MKS: joules kmol-1 K-1 = kg m2 s-2 kmol-1 K-1 - */ -double vcsUtil_gasConstant(int mu_units); - //! Definition of the function pointer for the root finder /*! * see vcsUtil_root1d for a definition of how to use this. diff --git a/include/cantera/equil/vcs_prob.h b/include/cantera/equil/vcs_prob.h index cb6dff107..8768dfdb5 100644 --- a/include/cantera/equil/vcs_prob.h +++ b/include/cantera/equil/vcs_prob.h @@ -48,7 +48,7 @@ public: size_t NPHASE0; //! Vector of chemical potentials of the species. This is a calculated - //! output quantity. length = number of species. units = m_VCS_UnitsFormat + //! output quantity. length = number of species. vector_fp m_gibbsSpecies; //! Total number of moles of the kth species. @@ -59,8 +59,6 @@ public: * * On output, this contains the solution for the total number of moles of * the kth species. - * - * units = m_VCS_UnitsFormat */ vector_fp w; @@ -71,7 +69,7 @@ public: //! Element abundances for jth element /*! * This is input from the input file and is considered a constant from - * thereon within the vcs_solve_TP(). units = m_VCS_UnitsFormat + * thereon within the vcs_solve_TP(). */ vector_fp gai; @@ -101,15 +99,10 @@ public: double T; //! Pressure - /*! - * units given by m_VCS_UnitsFormat - * -> are now PA - */ double PresPA; //! Volume of the entire system /*! - * units given by m_VCS_UnitsFormat * Note, this is an output variable atm */ double Vol; @@ -121,24 +114,6 @@ public: */ vector_fp VolPM; - //! Units for the chemical potential data, pressure data, volume, and - //! species amounts - /*! - * All internally stored quantities will have these units. Also, printed - * quantities will display in these units. - * - * | | | Chem_Pot | Pres | vol | moles| - * |---|----------------------|-------------------------|------|-------|------| - * |-1 | `VCS_UNITS_KCALMOL` | kcal/mol | atm | cm**3 | gmol | - * | 0 | `VCS_UNITS_UNITLESS` | MU / RT -> no units | atm | cm**3 | gmol | - * | 1 | `VCS_UNITS_KJMOL` | kJ / mol | atm | cm**3 | gmol | - * | 2 | `VCS_UNITS_KELVIN` | KELVIN -> MU / R | atm | cm**3 | gmol | - * | 3 | `VCS_UNITS_MKS` | Joules / Kmol (Cantera) | Pa | m**3 | kmol | - * - * see vcs_defs.h for more information - */ - int m_VCS_UnitsFormat; - //! Specification of the initial estimate method /*! * * 0: user estimate diff --git a/include/cantera/equil/vcs_solve.h b/include/cantera/equil/vcs_solve.h index 7ec9d8c24..483abdd7c 100644 --- a/include/cantera/equil/vcs_solve.h +++ b/include/cantera/equil/vcs_solve.h @@ -646,7 +646,7 @@ public: * 0 -> Dont print any intermediate results * @param maxit Maximum number of iterations for the algorithm * @param T Value of the Temperature (Kelvin) - * @param pres Value of the Pressure (units given by m_VCS_UnitsFormat variable + * @param pres Value of the Pressure * @return Returns an integer representing the success of the algorithm * * 0 = Equilibrium Achieved * * 1 = Range space error encountered. The element abundance criteria are @@ -893,26 +893,6 @@ public: */ int vcs_rearrange(); - //! Returns the multiplier for electric charge terms - /* - * This is basically equal to F/RT - * - * @param mu_units integer representing the dimensional units system - * @param TKelvin double Temperature in Kelvin - * @returns the value of F/RT - */ - double vcs_nondim_Farad(int mu_units, double TKelvin) const; - - //! Returns the multiplier for the nondimensionalization of the equations - /*! - * This is basically equal to RT - * - * @param mu_units integer representing the dimensional units system - * @param TKelvin double Temperature in Kelvin - * @returns the value of RT - */ - double vcs_nondimMult_TP(int mu_units, double TKelvin) const; - //! Nondimensionalize the problem data /*! * Nondimensionalize the free energies using the divisor, R * T @@ -938,14 +918,6 @@ public: */ void vcs_redim_TP(); - //! Print the string representing the Chemical potential units - /*! - * This gets printed using plogf() - * - * @param unitsFormat Integer representing the units system - */ - void vcs_printChemPotUnits(int unitsFormat) const; - //! Computes the current elemental abundances vector /*! * Computes the elemental abundances vector, m_elemAbundances[], and stores @@ -1636,18 +1608,7 @@ public: //! Temperature (Kelvin) double m_temperature; - //! Pressure (units are determined by m_VCS_UnitsFormat - /*! - * | Values | units | - * | ------ | ----- - * | -1: | atm | - * | 0: | atm | - * | 1: | atm | - * | 2: | atm | - * | 3: | Pa | - * - * Units being changed to Pa - */ + //! Pressure double m_pressurePA; //! Total kmoles of inert to add to each phase @@ -1869,7 +1830,7 @@ public: //! Partial molar volumes of the species /*! - * units = mks (m^3/kmol) -determined by m_VCS_UnitsFormat + * units = mks (m^3/kmol) * Length = number of species */ vector_fp m_PMVolumeSpecies; @@ -1901,18 +1862,6 @@ public: */ int m_timing_print_lvl; - //! Units for the chemical potential data - /*! - * | Value | chemical potential units | pressure units | - * | ----- | ------------------------ | -------------- | - * | -1 | kcal/mol | Pa | - * | 0 | MU/RT | Pa | - * | 1 | kJ/mol | Pa | - * | 2 | Kelvin | Pa | - * | 3 | J / kmol | Pa | - */ - int m_VCS_UnitsFormat; - friend class vcs_phaseStabilitySolve; }; diff --git a/include/cantera/equil/vcs_species_thermo.h b/include/cantera/equil/vcs_species_thermo.h index 3df449153..424d534b2 100644 --- a/include/cantera/equil/vcs_species_thermo.h +++ b/include/cantera/equil/vcs_species_thermo.h @@ -30,8 +30,7 @@ class vcs_VolPhase; /*! * Identifies the thermo model for the species. This structure is shared by * volumetric and surface species. However, each will have its own types of - * thermodynamic models. These quantities all have appropriate units. The units - * are specified by VCS_UnitsFormat. + * thermodynamic models. These quantities all have appropriate units. */ class VCS_SPECIES_THERMO { @@ -85,8 +84,6 @@ public: //! parameter that is used in the VCS_SSVOL_CONSTANT model. double SSStar_Vol0; - int m_VCS_UnitsFormat; - VCS_SPECIES_THERMO(size_t indexPhase, size_t indexSpeciesPhase); virtual ~VCS_SPECIES_THERMO() {} @@ -102,7 +99,7 @@ public: * * @param kspec species global index * @param TKelvin Temperature in Kelvin - * @param pres pressure is given in units specified by if__ variable. + * @param pres pressure in Pa * @return standard state free energy in units of Kelvin. */ virtual double GStar_R_calc(size_t kspec, double TKelvin, double pres); @@ -121,8 +118,7 @@ public: * This function calculates the standard state molar volume for species, * kspec, at the temperature TKelvin and pressure, Pres, * - * @return standard state volume in cm**3 per mol. - * (if__=3) m**3 / kmol + * @return standard state volume in m**3 / kmol */ virtual double VolStar_calc(size_t kglob, double TKelvin, double Pres); diff --git a/src/equil/vcs_MultiPhaseEquil.cpp b/src/equil/vcs_MultiPhaseEquil.cpp index 9b5d0da89..222451a6f 100644 --- a/src/equil/vcs_MultiPhaseEquil.cpp +++ b/src/equil/vcs_MultiPhaseEquil.cpp @@ -514,24 +514,8 @@ int vcs_MultiPhaseEquil::equilibrate_TP(int estimateEquil, plogf("\n"); plogf("----------------------------------------" "---------------------\n"); - plogf(" Name Mole_Number"); - if (m_vprob.m_VCS_UnitsFormat == VCS_UNITS_MKS) { - plogf("(kmol)"); - } else { - plogf("(gmol)"); - } - plogf(" Mole_Fraction Chem_Potential"); - if (m_vprob.m_VCS_UnitsFormat == VCS_UNITS_KCALMOL) { - plogf(" (kcal/mol)\n"); - } else if (m_vprob.m_VCS_UnitsFormat == VCS_UNITS_UNITLESS) { - plogf(" (Dimensionless)\n"); - } else if (m_vprob.m_VCS_UnitsFormat == VCS_UNITS_KJMOL) { - plogf(" (kJ/mol)\n"); - } else if (m_vprob.m_VCS_UnitsFormat == VCS_UNITS_KELVIN) { - plogf(" (Kelvin)\n"); - } else if (m_vprob.m_VCS_UnitsFormat == VCS_UNITS_MKS) { - plogf(" (J/kmol)\n"); - } + plogf(" Name Mole_Number(kmol)"); + plogf(" Mole_Fraction Chem_Potential (J/kmol)\n"); plogf("--------------------------------------------------" "-----------\n"); for (size_t i = 0; i < m_vprob.nspecies; i++) { @@ -713,7 +697,6 @@ int vcs_Cantera_to_vprob(MultiPhase* mphase, VCS_PROB* vprob) vprob->nspecies = totNumSpecies; vprob->ne = 0; vprob->NPhase = totNumPhases; - vprob->m_VCS_UnitsFormat = VCS_UNITS_MKS; // Set the initial estimate to a machine generated estimate for now // We will work out the details later. vprob->iest = -1; @@ -752,7 +735,6 @@ int vcs_Cantera_to_vprob(MultiPhase* mphase, VCS_PROB* vprob) VolPhase->m_gasPhase = gasPhase; // Tell the vcs_VolPhase pointer about cantera - VolPhase->p_VCS_UnitsFormat = vprob->m_VCS_UnitsFormat; VolPhase->setPtrThermoPhase(tPhase); VolPhase->setTotalMoles(0.0); @@ -874,7 +856,6 @@ int vcs_Cantera_to_vprob(MultiPhase* mphase, VCS_PROB* vprob) // Transfer the thermo specification of the species // vprob->SpeciesThermo[] - ts_ptr->m_VCS_UnitsFormat = VolPhase->p_VCS_UnitsFormat; // Add lookback connectivity into the thermo object first ts_ptr->IndexPhase = iphase; @@ -945,11 +926,10 @@ int vcs_Cantera_to_vprob(MultiPhase* mphase, VCS_PROB* vprob) // Now, calculate a sample naught Gibbs free energy calculation // at the specified temperature. - double R = vcsUtil_gasConstant(vprob->m_VCS_UnitsFormat); 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)/ R; + ts_ptr->SS0_feSave = VolPhase->G0_calc_one(k)/ GasConstant; ts_ptr->SS0_TSave = vprob->T; } } @@ -1214,24 +1194,8 @@ int vcs_MultiPhaseEquil::determine_PhaseStability(int iph, double& funcStab, int plogf("\n"); plogf("----------------------------------------" "---------------------\n"); - plogf(" Name Mole_Number"); - if (m_vprob.m_VCS_UnitsFormat == VCS_UNITS_MKS) { - plogf("(kmol)"); - } else { - plogf("(gmol)"); - } - plogf(" Mole_Fraction Chem_Potential"); - if (m_vprob.m_VCS_UnitsFormat == VCS_UNITS_KCALMOL) { - plogf(" (kcal/mol)\n"); - } else if (m_vprob.m_VCS_UnitsFormat == VCS_UNITS_UNITLESS) { - plogf(" (Dimensionless)\n"); - } else if (m_vprob.m_VCS_UnitsFormat == VCS_UNITS_KJMOL) { - plogf(" (kJ/mol)\n"); - } else if (m_vprob.m_VCS_UnitsFormat == VCS_UNITS_KELVIN) { - plogf(" (Kelvin)\n"); - } else if (m_vprob.m_VCS_UnitsFormat == VCS_UNITS_MKS) { - plogf(" (J/kmol)\n"); - } + plogf(" Name Mole_Number(kmol)"); + plogf(" Mole_Fraction Chem_Potential (J/kmol)\n"); plogf("-------------------------------------------------------------\n"); for (size_t i = 0; i < m_vprob.nspecies; i++) { plogf("%-12s", m_vprob.SpName[i]); diff --git a/src/equil/vcs_TP.cpp b/src/equil/vcs_TP.cpp index 23b31209a..15afa5150 100644 --- a/src/equil/vcs_TP.cpp +++ b/src/equil/vcs_TP.cpp @@ -45,24 +45,12 @@ int VCS_SOLVE::vcs_TP(int ipr, int ip1, int maxit, double T_arg, double pres_arg int VCS_SOLVE::vcs_evalSS_TP(int ipr, int ip1, double Temp, double pres) { - // We need to special case VCS_UNITS_UNITLESS, here. cpc_ts_GStar_calc() - // returns units of Kelvin. Also, the temperature comes into play in - // calculating the ideal equation of state contributions, and other - // equations of state also. Therefore, we will emulate the VCS_UNITS_KELVIN - // case, here by changing the initial Gibbs free energy units to Kelvin - // before feeding them to the cpc_ts_GStar_calc() routine. Then, we will - // revert them back to unitless at the end of this routine. for (size_t iph = 0; iph < m_numPhases; iph++) { vcs_VolPhase* vph = m_VolPhaseList[iph]; vph->setState_TP(m_temperature, m_pressurePA); vph->sendToVCS_GStar(&m_SSfeSpecies[0]); } - if (m_VCS_UnitsFormat == VCS_UNITS_UNITLESS) { - for (size_t i = 0; i < m_numSpeciesTot; ++i) { - m_SSfeSpecies[i] /= Temp; - } - } return VCS_SUCCESS; } diff --git a/src/equil/vcs_VolPhase.cpp b/src/equil/vcs_VolPhase.cpp index 3e74f1a85..7ed499638 100644 --- a/src/equil/vcs_VolPhase.cpp +++ b/src/equil/vcs_VolPhase.cpp @@ -26,7 +26,6 @@ vcs_VolPhase::vcs_VolPhase(VCS_SOLVE* owningSolverObject) : m_gasPhase(false), m_eqnState(VCS_EOS_CONSTANT), ChargeNeutralityElement(npos), - p_VCS_UnitsFormat(VCS_UNITS_MKS), p_activityConvention(0), m_numElemConstraints(0), m_elemGlobalIndex(0), @@ -68,7 +67,6 @@ vcs_VolPhase::vcs_VolPhase(const vcs_VolPhase& b) : m_gasPhase(b.m_gasPhase), m_eqnState(b.m_eqnState), ChargeNeutralityElement(b.ChargeNeutralityElement), - p_VCS_UnitsFormat(b.p_VCS_UnitsFormat), p_activityConvention(b.p_activityConvention), m_numElemConstraints(b.m_numElemConstraints), m_numSpecies(b.m_numSpecies), @@ -113,7 +111,6 @@ vcs_VolPhase& vcs_VolPhase::operator=(const vcs_VolPhase& b) m_gasPhase = b.m_gasPhase; m_eqnState = b.m_eqnState; ChargeNeutralityElement = b.ChargeNeutralityElement; - p_VCS_UnitsFormat = b.p_VCS_UnitsFormat; p_activityConvention= b.p_activityConvention; m_numSpecies = b.m_numSpecies; m_numElemConstraints = b.m_numElemConstraints; @@ -715,7 +712,6 @@ void vcs_VolPhase::setPtrThermoPhase(ThermoPhase* tp_ptr) Temp_ = TP_ptr->temperature(); Pres_ = TP_ptr->pressure(); setState_TP(Temp_, Pres_); - p_VCS_UnitsFormat = VCS_UNITS_MKS; m_phi = TP_ptr->electricPotential(); size_t nsp = TP_ptr->nSpecies(); size_t nelem = TP_ptr->nElements(); diff --git a/src/equil/vcs_nondim.cpp b/src/equil/vcs_nondim.cpp index b1d386ab1..4c4d750d9 100644 --- a/src/equil/vcs_nondim.cpp +++ b/src/equil/vcs_nondim.cpp @@ -15,51 +15,12 @@ namespace Cantera { -double VCS_SOLVE::vcs_nondim_Farad(int mu_units, double TKelvin) const -{ - if (TKelvin <= 0.0) { - TKelvin = 293.15; - } - switch (mu_units) { - case VCS_UNITS_MKS: - case VCS_UNITS_KJMOL: - case VCS_UNITS_KCALMOL: - return ElectronCharge * Avogadro / (TKelvin * GasConstant); - case VCS_UNITS_UNITLESS: - return ElectronCharge * Avogadro; - case VCS_UNITS_KELVIN: - return ElectronCharge * Avogadro/ TKelvin; - default: - throw CanteraError("vcs_nondim_Farad", "unknown units: {}", mu_units); - } -} - -double VCS_SOLVE::vcs_nondimMult_TP(int mu_units, double TKelvin) const -{ - if (TKelvin <= 0.0) { - TKelvin = 293.15; - } - switch (mu_units) { - case VCS_UNITS_KCALMOL: - return TKelvin * GasConst_cal_mol_K * 1e-3; - case VCS_UNITS_UNITLESS: - return 1.0; - case VCS_UNITS_KJMOL: - return TKelvin * GasConstant * 1e-6; - case VCS_UNITS_KELVIN: - return TKelvin; - case VCS_UNITS_MKS: - return TKelvin * GasConstant; - default: - throw CanteraError("vcs_nondimMult_TP", "unknown units: {}", mu_units); - } -} void VCS_SOLVE::vcs_nondim_TP() { if (m_unitsState == VCS_DIMENSIONAL_G) { m_unitsState = VCS_NONDIMENSIONAL_G; - double tf = 1.0 / vcs_nondimMult_TP(m_VCS_UnitsFormat, m_temperature); + double tf = 1.0 / (GasConstant * m_temperature); for (size_t i = 0; i < m_numSpeciesTot; ++i) { // Modify the standard state and total chemical potential data, // FF(I), to make it dimensionless, i.e., mu / RT. Thus, we may @@ -70,7 +31,7 @@ void VCS_SOLVE::vcs_nondim_TP() m_feSpecies_old[i] *= tf; } - m_Faraday_dim = vcs_nondim_Farad(m_VCS_UnitsFormat, m_temperature); + m_Faraday_dim = ElectronCharge * Avogadro / (m_temperature * GasConstant); // Scale the total moles if necessary: First find out the total moles double tmole_orig = vcs_tmoles(); @@ -104,26 +65,24 @@ void VCS_SOLVE::vcs_nondim_TP() } if (m_totalMoleScale != 1.0) { - if (m_VCS_UnitsFormat == VCS_UNITS_MKS) { - if (m_debug_print_lvl >= 2) { - plogf(" --- vcs_nondim_TP() called: USING A MOLE SCALE OF %g until further notice", m_totalMoleScale); - plogendl(); - } - for (size_t i = 0; i < m_numSpeciesTot; ++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) { - m_elemAbundancesGoal[i] *= (1.0 / m_totalMoleScale); + if (m_debug_print_lvl >= 2) { + plogf(" --- vcs_nondim_TP() called: USING A MOLE SCALE OF %g until further notice", m_totalMoleScale); + plogendl(); + } + for (size_t i = 0; i < m_numSpeciesTot; ++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) { + m_elemAbundancesGoal[i] *= (1.0 / m_totalMoleScale); + } - for (size_t iph = 0; iph < m_numPhases; iph++) { - TPhInertMoles[iph] *= (1.0 / m_totalMoleScale); - if (TPhInertMoles[iph] != 0.0) { - vcs_VolPhase* vphase = m_VolPhaseList[iph]; - vphase->setTotalMolesInert(TPhInertMoles[iph]); - } + for (size_t iph = 0; iph < m_numPhases; iph++) { + TPhInertMoles[iph] *= (1.0 / m_totalMoleScale); + if (TPhInertMoles[iph] != 0.0) { + vcs_VolPhase* vphase = m_VolPhaseList[iph]; + vphase->setTotalMolesInert(TPhInertMoles[iph]); } } vcs_tmoles(); @@ -135,7 +94,7 @@ void VCS_SOLVE::vcs_redim_TP() { if (m_unitsState != VCS_DIMENSIONAL_G) { m_unitsState = VCS_DIMENSIONAL_G; - double tf = vcs_nondimMult_TP(m_VCS_UnitsFormat, m_temperature); + double tf = m_temperature * GasConstant; for (size_t i = 0; i < m_numSpeciesTot; ++i) { // Modify the standard state and total chemical potential data, @@ -148,52 +107,27 @@ void VCS_SOLVE::vcs_redim_TP() m_Faraday_dim *= tf; } if (m_totalMoleScale != 1.0) { - if (m_VCS_UnitsFormat == VCS_UNITS_MKS) { - if (m_debug_print_lvl >= 2) { - plogf(" --- vcs_redim_TP() called: getting rid of mole scale of %g", m_totalMoleScale); - plogendl(); - } - for (size_t i = 0; i < m_numSpeciesTot; ++i) { - if (m_speciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - m_molNumSpecies_old[i] *= m_totalMoleScale; - } - } - for (size_t i = 0; i < m_numElemConstraints; ++i) { - m_elemAbundancesGoal[i] *= m_totalMoleScale; - } - - for (size_t iph = 0; iph < m_numPhases; iph++) { - TPhInertMoles[iph] *= m_totalMoleScale; - if (TPhInertMoles[iph] != 0.0) { - vcs_VolPhase* vphase = m_VolPhaseList[iph]; - vphase->setTotalMolesInert(TPhInertMoles[iph]); - } - } - vcs_tmoles(); + if (m_debug_print_lvl >= 2) { + plogf(" --- vcs_redim_TP() called: getting rid of mole scale of %g", m_totalMoleScale); + plogendl(); + } + for (size_t i = 0; i < m_numSpeciesTot; ++i) { + if (m_speciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { + m_molNumSpecies_old[i] *= m_totalMoleScale; + } + } + for (size_t i = 0; i < m_numElemConstraints; ++i) { + m_elemAbundancesGoal[i] *= m_totalMoleScale; } - } -} -void VCS_SOLVE::vcs_printChemPotUnits(int unitsFormat) const -{ - switch (unitsFormat) { - case VCS_UNITS_KCALMOL: - plogf("kcal/gmol"); - break; - case VCS_UNITS_UNITLESS: - plogf("dimensionless"); - break; - case VCS_UNITS_KJMOL: - plogf("kJ/gmol"); - break; - case VCS_UNITS_KELVIN: - plogf("Kelvin"); - break; - case VCS_UNITS_MKS: - plogf("J/kmol"); - break; - default: - throw CanteraError("VCS_SOLVE::vcs_printChemPotUnits", "unknown units!"); + for (size_t iph = 0; iph < m_numPhases; iph++) { + TPhInertMoles[iph] *= m_totalMoleScale; + if (TPhInertMoles[iph] != 0.0) { + vcs_VolPhase* vphase = m_VolPhaseList[iph]; + vphase->setTotalMolesInert(TPhInertMoles[iph]); + } + } + vcs_tmoles(); } } diff --git a/src/equil/vcs_prob.cpp b/src/equil/vcs_prob.cpp index 815b38e66..8d5652628 100644 --- a/src/equil/vcs_prob.cpp +++ b/src/equil/vcs_prob.cpp @@ -34,8 +34,6 @@ VCS_PROB::VCS_PROB(size_t nsp, size_t nel, size_t nph) : T(298.15), PresPA(1.0), Vol(0.0), - // Set the units for the chemical potential data to be unitless - m_VCS_UnitsFormat(VCS_UNITS_UNITLESS), // The default is to not expect an initial estimate of the species // concentrations iest(-1), @@ -224,30 +222,14 @@ void VCS_PROB::prob_report(int print_lvl) plogf("\nElemental Abundances: "); plogf(" Target_kmol ElemType ElActive\n"); - double fac = 1.0; - if (m_VCS_UnitsFormat == VCS_UNITS_MKS) { - fac = 1.0; - } for (size_t i = 0; i < ne; ++i) { writeline(' ', 26, false); plogf("%-2.2s", ElName[i]); - plogf("%20.12E ", fac * gai[i]); + plogf("%20.12E ", gai[i]); plogf("%3d %3d\n", m_elType[i], ElActive[i]); } - plogf("\nChemical Potentials: "); - if (m_VCS_UnitsFormat == VCS_UNITS_UNITLESS) { - plogf("(unitless)"); - } else if (m_VCS_UnitsFormat == VCS_UNITS_KCALMOL) { - plogf("(kcal/gmol)"); - } else if (m_VCS_UnitsFormat == VCS_UNITS_KJMOL) { - plogf("(kJ/gmol)"); - } else if (m_VCS_UnitsFormat == VCS_UNITS_KELVIN) { - plogf("(Kelvin)"); - } else if (m_VCS_UnitsFormat == VCS_UNITS_MKS) { - plogf("(J/kmol)"); - } - plogf("\n"); + plogf("\nChemical Potentials: (J/kmol)\n"); plogf(" Species (phase) " " SS0ChemPot StarChemPot\n"); for (size_t iphase = 0; iphase < NPhase; iphase++) { diff --git a/src/equil/vcs_report.cpp b/src/equil/vcs_report.cpp index 10b51fbc4..d8d8e7124 100644 --- a/src/equil/vcs_report.cpp +++ b/src/equil/vcs_report.cpp @@ -252,10 +252,7 @@ int VCS_SOLVE::vcs_report(int iconv) writeline('-', 93, true, true); plogf("Chemical Potentials of the Species: (dimensionless)\n"); - double rt = vcs_nondimMult_TP(m_VCS_UnitsFormat, m_temperature); - plogf("\t\t(RT = %g ", rt); - vcs_printChemPotUnits(m_VCS_UnitsFormat); - plogf(")\n"); + plogf("\t\t(RT = %g J/kmol)\n", GasConstant * m_temperature); plogf(" Name TKMoles StandStateChemPot " " ln(AC) ln(X_i) | F z_i phi | ChemPot | (-lnMnaught)"); plogf("| (MolNum ChemPot)|"); diff --git a/src/equil/vcs_solve.cpp b/src/equil/vcs_solve.cpp index 8ef199853..ed2d50f6c 100644 --- a/src/equil/vcs_solve.cpp +++ b/src/equil/vcs_solve.cpp @@ -49,8 +49,7 @@ VCS_SOLVE::VCS_SOLVE() : m_Faraday_dim(ElectronCharge * Avogadro), m_VCount(0), m_debug_print_lvl(0), - m_timing_print_lvl(1), - m_VCS_UnitsFormat(VCS_UNITS_UNITLESS) + m_timing_print_lvl(1) { } @@ -88,8 +87,6 @@ void VCS_SOLVE::vcs_initSizes(const size_t nspecies0, const size_t nelements, " Number of species is nonpositive\n"); } - m_VCS_UnitsFormat = VCS_UNITS_UNITLESS; - /* * We will initialize sc[] to note the fact that it needs to be * filled with meaningful information. @@ -464,9 +461,6 @@ int VCS_SOLVE::vcs_prob_specifyFully(const VCS_PROB* pub) TPhInertMoles[iph] = Vphase->totalMolesInert(); } - // if__ : Copy over the units for the chemical potential - m_VCS_UnitsFormat = pub->m_VCS_UnitsFormat; - // tolerance requirements -> copy them over here and later m_tolmaj = pub->tolmaj; m_tolmin = pub->tolmin; @@ -629,7 +623,6 @@ int VCS_SOLVE::vcs_prob_specify(const VCS_PROB* pub) m_temperature = pub->T; m_pressurePA = pub->PresPA; - m_VCS_UnitsFormat = pub->m_VCS_UnitsFormat; m_doEstimateEquil = pub->iest; m_totalVol = pub->Vol; m_tolmaj = pub->tolmaj; diff --git a/src/equil/vcs_solve_TP.cpp b/src/equil/vcs_solve_TP.cpp index 503f4f3cf..59f797d84 100644 --- a/src/equil/vcs_solve_TP.cpp +++ b/src/equil/vcs_solve_TP.cpp @@ -111,21 +111,7 @@ int VCS_SOLVE::vcs_solve_TP(int print_lvl, int printDetails, int maxit) if (m_doEstimateEquil == 0) { plogf("\n USER ESTIMATE OF EQUILIBRIUM\n"); } - if (m_VCS_UnitsFormat == VCS_UNITS_KCALMOL) { - plogf(" Stan. Chem. Pot. in kcal/mole\n"); - } - if (m_VCS_UnitsFormat == VCS_UNITS_UNITLESS) { - plogf(" Stan. Chem. Pot. is MU/RT\n"); - } - if (m_VCS_UnitsFormat == VCS_UNITS_KJMOL) { - plogf(" Stan. Chem. Pot. in KJ/mole\n"); - } - if (m_VCS_UnitsFormat == VCS_UNITS_KELVIN) { - plogf(" Stan. Chem. Pot. in Kelvin\n"); - } - if (m_VCS_UnitsFormat == VCS_UNITS_MKS) { - plogf(" Stan. Chem. Pot. in J/kmol\n"); - } + plogf(" Stan. Chem. Pot. in J/kmol\n"); plogf("\n SPECIES FORMULA VECTOR "); writeline(' ', 41, false); plogf(" STAN_CHEM_POT EQUILIBRIUM_EST. Species_Type\n\n"); @@ -134,7 +120,7 @@ int VCS_SOLVE::vcs_solve_TP(int print_lvl, int printDetails, int maxit) plogf("%-4.4s ", m_elementName[i]); } plogf(" PhaseID\n"); - double RT = vcs_nondimMult_TP(m_VCS_UnitsFormat, m_temperature); + double RT = GasConstant * m_temperature; for (size_t i = 0; i < m_numSpeciesTot; ++i) { plogf(" %-18.18s", m_speciesName[i]); for (size_t j = 0; j < m_numElemConstraints; ++j) { diff --git a/src/equil/vcs_species_thermo.cpp b/src/equil/vcs_species_thermo.cpp index 5e8c12ab6..86d7d7889 100644 --- a/src/equil/vcs_species_thermo.cpp +++ b/src/equil/vcs_species_thermo.cpp @@ -33,8 +33,7 @@ VCS_SPECIES_THERMO::VCS_SPECIES_THERMO(size_t indexPhase, SS0_Pref(1.01325E5), SSStar_Model(VCS_SSSTAR_CONSTANT), SSStar_Vol_Model(VCS_SSVOL_IDEALGAS), - SSStar_Vol0(-1.0), - m_VCS_UnitsFormat(VCS_UNITS_UNITLESS) + SSStar_Vol0(-1.0) { SS0_Pref = 1.01325E5; } @@ -53,8 +52,7 @@ VCS_SPECIES_THERMO::VCS_SPECIES_THERMO(const VCS_SPECIES_THERMO& b) : SS0_Pref(b.SS0_Pref), SSStar_Model(b.SSStar_Model), SSStar_Vol_Model(b.SSStar_Vol_Model), - SSStar_Vol0(b.SSStar_Vol0), - m_VCS_UnitsFormat(b.m_VCS_UnitsFormat) + SSStar_Vol0(b.SSStar_Vol0) { } @@ -76,7 +74,6 @@ VCS_SPECIES_THERMO::operator=(const VCS_SPECIES_THERMO& b) SSStar_Model = b.SSStar_Model; SSStar_Vol_Model = b.SSStar_Vol_Model; SSStar_Vol0 = b.SSStar_Vol0; - m_VCS_UnitsFormat = b.m_VCS_UnitsFormat; } return *this; } @@ -90,23 +87,14 @@ double VCS_SPECIES_THERMO::GStar_R_calc(size_t kglob, double TKelvin, double pres) { double fe = G0_R_calc(kglob, TKelvin); - if (m_VCS_UnitsFormat != VCS_UNITS_MKS) { - throw CanteraError("VCS_SPECIES_THERMO::GStar_R_calc", - "Possible inconsistency"); - } OwningPhase->setState_TP(TKelvin, pres); fe = OwningPhase->GStar_calc_one(IndexSpeciesPhase); - double R = vcsUtil_gasConstant(m_VCS_UnitsFormat); - return fe / R; + return fe / GasConstant; } double VCS_SPECIES_THERMO::VolStar_calc(size_t kglob, double TKelvin, double presPA) { - if (m_VCS_UnitsFormat != VCS_UNITS_MKS) { - throw CanteraError("VCS_SPECIES_THERMO::VolStar_calc", - "Possible inconsistency"); - } OwningPhase->setState_TP(TKelvin, presPA); return OwningPhase->VolStar_calc_one(IndexSpeciesPhase); } @@ -119,14 +107,9 @@ double VCS_SPECIES_THERMO::G0_R_calc(size_t kglob, double TKelvin) if (TKelvin == SS0_TSave) { return SS0_feSave; } - if (m_VCS_UnitsFormat != VCS_UNITS_MKS) { - throw CanteraError("VCS_SPECIES_THERMO::G0_R_calc", - "Possible inconsistency"); - } OwningPhase->setState_T(TKelvin); double fe = OwningPhase->G0_calc_one(IndexSpeciesPhase); - double R = vcsUtil_gasConstant(m_VCS_UnitsFormat); - fe /= R; + fe /= GasConstant; SS0_feSave = fe; SS0_TSave = TKelvin; return fe; diff --git a/src/equil/vcs_util.cpp b/src/equil/vcs_util.cpp index 719f94f01..7b555bb6c 100644 --- a/src/equil/vcs_util.cpp +++ b/src/equil/vcs_util.cpp @@ -57,25 +57,6 @@ size_t vcs_optMax(const double* x, const double* xSize, size_t j, size_t n) return largest; } -double vcsUtil_gasConstant(int mu_units) -{ - switch (mu_units) { - case VCS_UNITS_KCALMOL: - return GasConst_cal_mol_K * 1e-3; - case VCS_UNITS_UNITLESS: - return 1.0; - case VCS_UNITS_KJMOL: - return GasConstant * 1e-6; - case VCS_UNITS_KELVIN: - return 1.0; - case VCS_UNITS_MKS: - // joules / kg-mol K = kg m2 / s2 kg-mol K - return GasConstant; - default: - throw CanteraError("vcsUtil_gasConstant", "uknown units: {}", mu_units); - } -} - const char* vcs_speciesType_string(int speciesStatus, int length) { switch (speciesStatus) {