From 0bd88927f5e0703277390c94e3dafedd24391b3c Mon Sep 17 00:00:00 2001
From: Ray Speth <speth@mit.edu>
Date: Wed, 17 Aug 2016 19:58:24 -0400
Subject: [PATCH] [Equil] Fix ElemRearrange when unused elements are present

The solver could get stuck in an infinite loop when there were elements that
were not used in any species.
---
 .../cython/cantera/test/test_equilibrium.py   | 27 +++++++++++++++++++
 src/equil/BasisOptimize.cpp                   |  1 +
 2 files changed, 28 insertions(+)

diff --git a/interfaces/cython/cantera/test/test_equilibrium.py b/interfaces/cython/cantera/test/test_equilibrium.py
index 4623bb8ca..bc4433633 100644
--- a/interfaces/cython/cantera/test/test_equilibrium.py
+++ b/interfaces/cython/cantera/test/test_equilibrium.py
@@ -100,6 +100,33 @@ class MultiphaseEquilTest(EquilTestCases, utilities.CanteraTest):
         self.check(gas, CH4=0, O2=1, H2O=2, CO2=1)
 
 
+class EquilExtraElements(utilities.CanteraTest):
+    def setUp(self):
+        s = """ideal_gas(elements='H Ar C O Cl N',
+                         species='gri30: AR N2 CH4 O2 CO2 H2O CO H2 OH')"""
+        self.gas = ct.Solution(source=s)
+        self.gas.TP = 300, 101325
+        self.gas.set_equivalence_ratio(0.8, 'CH4', 'O2:1.0, N2:3.76')
+
+    def test_auto(self):
+        # Succeeds after falling back to VCS
+        self.gas.equilibrate('TP')
+        self.assertNear(self.gas['CH4'].X[0], 0.0)
+
+    @unittest.expectedFailure
+    def test_element_potential(self):
+        self.gas.equilibrate('TP', solver='element_potential')
+        self.assertNear(self.gas['CH4'].X[0], 0.0)
+
+    def test_gibbs(self):
+        self.gas.equilibrate('TP', solver='gibbs')
+        self.assertNear(self.gas['CH4'].X[0], 0.0)
+
+    def test_vcs(self):
+        self.gas.equilibrate('TP', solver='vcs')
+        self.assertNear(self.gas['CH4'].X[0], 0.0)
+
+
 class VCS_EquilTest(EquilTestCases, utilities.CanteraTest):
     def __init__(self, *args, **kwargs):
         EquilTestCases.__init__(self, 'vcs')
diff --git a/src/equil/BasisOptimize.cpp b/src/equil/BasisOptimize.cpp
index 0e90e0f78..04d0658ad 100644
--- a/src/equil/BasisOptimize.cpp
+++ b/src/equil/BasisOptimize.cpp
@@ -350,6 +350,7 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances,
             // Search the element vector. We first locate elements that are
             // present in any amount. Then, we locate elements that are not
             // present in any amount. Return its identity in K.
+            k = nelements;
             size_t kk;
             for (size_t ielem = jr; ielem < nelements; ielem++) {
                 kk = orderVectorElements[ielem];