[Thermo] Add property to call species charges in phase

AUTHORS

@@ -24,6 +24,7 @@ John Hewson (@jchewson), Sandia National Laboratory
 Trevor Hickey (@luxe)
 Yuanjie Jiang
 Benjamin Kee (@lionkey)
+Daniel Korff (@korffdm), Colorado School of Mines
 Jon Kristofer
 Samesh Lakothia (@sameshl)
 Kyle Linevitch, Jr. (@KyleLinevitchJr)

include/cantera/clib/ct.h

@@ -40,6 +40,7 @@ extern "C" {
     CANTERA_CAPI int thermo_setMassFractionsByName(int n, const char* y);
     CANTERA_CAPI int thermo_getAtomicWeights(int n, size_t lenm, double* atw);
     CANTERA_CAPI int thermo_getMolecularWeights(int n, size_t lenm, double* mw);
+    CANTERA_CAPI int thermo_getSpeciesCharges(int n, size_t lenm, double* sc);
     CANTERA_CAPI int thermo_getElementName(int n, size_t k, size_t lennm, char* nm);
     CANTERA_CAPI int thermo_getSpeciesName(int n, size_t m, size_t lennm, char* nm);
     CANTERA_CAPI int thermo_getName(int n, size_t lennm, char* nm);

include/cantera/cython/wrappers.h

@@ -73,6 +73,7 @@ THERMO_1D(getConcentrations)
 THERMO_1D(setConcentrations)
 
 THERMO_1D(getMolecularWeights)
+THERMO_1D(getSpeciesCharges)
 THERMO_1D(getChemPotentials)
 THERMO_1D(getElectrochemPotentials)
 THERMO_1D(getPartialMolarEnthalpies)

include/cantera/thermo/Phase.h

@@ -471,6 +471,10 @@ class Phase
     //! units = kg / kmol
     const vector_fp& molecularWeights() const;
 
+    //! Copy the vector of species charges into array charges.
+    //!     @param charges Output array of species charges (elem. charge)
+    const vector_fp& getSpeciesCharges(double* charges) const;
+
     /// @name Composition
     //@{
 

interfaces/cython/cantera/_cantera.pxd

@@ -906,6 +906,7 @@ cdef extern from "cantera/cython/wrappers.h":
 
     # other ThermoPhase methods
     cdef void thermo_getMolecularWeights(CxxThermoPhase*, double*) except +translate_exception
+    cdef void thermo_getSpeciesCharges(CxxThermoPhase*, double*) except +translate_exception
 
     # Kinetics per-reaction properties
     cdef void kin_getFwdRatesOfProgress(CxxKinetics*, double*) except +translate_exception

interfaces/cython/cantera/composite.py

@@ -483,7 +483,7 @@ class SolutionArray:
         'element_name', 'element_names', 'atomic_weight', 'atomic_weights',
         'n_species', 'species_name', 'species_names', 'species_index',
         'species', 'n_atoms', 'molecular_weights', 'min_temp', 'max_temp',
-        'reference_pressure',
+        'reference_pressure', 'species_charges',
         # From Kinetics
         'n_total_species', 'n_reactions', 'n_phases', 'reaction_phase_index',
         'kinetics_species_index', 'reaction', 'reactions', 'modify_reaction',

interfaces/cython/cantera/test/test_thermo.py

@@ -134,6 +134,15 @@ def test_weights(self):
                 test_weight += aw * self.phase.n_atoms(i,j)
             self.assertNear(test_weight, mw)
 
+    def test_charges(self):
+        gas = ct.Solution('ch4_ion.yaml')
+        charges = gas.species_charges
+        test = {'E': -1., 'N2': 0., 'H3O+': 1.}
+        for species, charge in test.items():
+            self.assertIn(species, gas.species_names)
+            index = gas.species_index(species)
+            self.assertEqual(charges[index], charge)
+
     def test_setComposition(self):
         X = np.zeros(self.phase.n_species)
         X[2] = 1.0
@@ -1657,7 +1666,7 @@ def test_extra(self):
         states = ct.SolutionArray(self.gas, 10, extra=extra)
         keys = list(states._extra.keys())
         self.assertEqual(keys[0], 'grid')
-        
+
         with self.assertRaises(ValueError):
             states = ct.SolutionArray(self.gas, extra=['creation_rates'])
 

interfaces/cython/cantera/thermo.pyx

@@ -660,6 +660,11 @@ cdef class ThermoPhase(_SolutionBase):
         def __get__(self):
             return self._getArray1(thermo_getMolecularWeights)
 
+    property species_charges:
+        """Array of species charges [-]."""
+        def __get__(self):
+            return self._getArray1(thermo_getSpeciesCharges)
+
     property mean_molecular_weight:
         """The mean molecular weight (molar mass) [kg/kmol]."""
         def __get__(self):

interfaces/matlab/toolbox/@ThermoPhase/speciesCharges.m

@@ -0,0 +1,12 @@
+function ch = speciesCharges(tp)
+% SPECIESCHARGES  Get the elementary charge of the species.
+% x = speciesCharges(tp)
+%
+% :param tp:
+%     Instance of class :mat:func:`ThermoPhase` (or another
+%     object that derives from ThermoPhase)
+% :return:
+%     Vector of species charges. Units: elem. charge
+%
+
+ch = phase_get(tp.tp_id, 23);

src/clib/ct.cpp

@@ -284,6 +284,18 @@ extern "C" {
         }
     }
 
+    int thermo_getSpeciesCharges(int n, size_t lenm, double* sc)
+    {
+        try {
+            ThermoPhase& p = ThermoCabinet::item(n);
+            p.checkSpeciesArraySize(lenm);
+            const vector_fp& ch = p.getSpeciesCharges(sc);
+            return 0;
+        } catch (...) {
+            return handleAllExceptions(-1, ERR);
+        }
+    }
+
     int thermo_getName(int n, size_t lennm, char* nm)
     {
         try {

src/matlab/phasemethods.cpp

@@ -196,6 +196,9 @@ void phasemethods(int nlhs, mxArray* plhs[],
         case 22:
             iok = thermo_getMolecularWeights(ph,nsp, &x[0]);
             break;
+        case 23:
+            iok = thermo_getSpeciesCharges(ph, nsp, &x[0]);
+            break;
         default:
             mexErrMsgTxt("Unknown job number");
         }

src/thermo/Phase.cpp

@@ -559,6 +559,13 @@ const vector_fp& Phase::molecularWeights() const
     return m_molwts;
 }
 
+const vector_fp& Phase::getSpeciesCharges(double* charges) const
+{
+    const vector_fp& ch = m_speciesCharge;
+    copy(ch.begin(), ch.end(), charges);
+    return m_speciesCharge;
+}
+
 compositionMap Phase::getMoleFractionsByName(double threshold) const
 {
     compositionMap comp;