From 178b452a49b431c39368baee05ac13b19d4f54c4 Mon Sep 17 00:00:00 2001
From: Ray Speth <speth@mit.edu>
Date: Thu, 8 Sep 2016 22:32:25 -0400
Subject: [PATCH] [Transport] Automatically reinitialize after change in number
 of species

This makes it possible for GasTransport objects to work reasonably well with
ThermoPhase objects that have had species added to them. The main changes
required are to make the initialization functions idempotent.
---
 .../cython/cantera/test/test_transport.py     | 45 +++++++++++++++++++
 src/transport/GasTransport.cpp                | 13 ++++++
 src/transport/MixTransport.cpp                |  2 +-
 src/transport/MultiTransport.cpp              |  2 +-
 4 files changed, 60 insertions(+), 2 deletions(-)

diff --git a/interfaces/cython/cantera/test/test_transport.py b/interfaces/cython/cantera/test/test_transport.py
index 17fb99db2..c033ca3c1 100644
--- a/interfaces/cython/cantera/test/test_transport.py
+++ b/interfaces/cython/cantera/test/test_transport.py
@@ -43,6 +43,51 @@ class TestTransport(utilities.CanteraTest):
         self.assertTrue(all(self.phase.multi_diff_coeffs.flat >= 0.0))
         self.assertTrue(all(self.phase.thermal_diff_coeffs.flat != 0.0))
 
+    def test_add_species_mix(self):
+        S = {s.name: s for s in ct.Species.listFromFile('gri30.xml')}
+
+        base = ['H', 'H2', 'OH', 'O2', 'AR']
+        extra = ['H2O', 'CH4']
+
+        state = 500, 2e5, 'H2:0.4, O2:0.29, CH4:0.01, H2O:0.3'
+
+        gas1 = ct.Solution(thermo='IdealGas', species=[S[s] for s in base+extra])
+        gas1.transport_model = 'Mix'
+        gas1.TPX = state
+
+        gas2 = ct.Solution(thermo='IdealGas', species=[S[s] for s in base])
+        gas2.transport_model = 'Mix'
+        for s in extra:
+            gas2.add_species(S[s])
+        gas2.TPX = state
+
+        self.assertNear(gas1.viscosity, gas2.viscosity)
+        self.assertNear(gas1.thermal_conductivity, gas2.thermal_conductivity)
+        self.assertArrayNear(gas1.binary_diff_coeffs, gas2.binary_diff_coeffs)
+        self.assertArrayNear(gas1.mix_diff_coeffs, gas2.mix_diff_coeffs)
+
+    def test_add_species_multi(self):
+        S = {s.name: s for s in ct.Species.listFromFile('gri30.xml')}
+
+        base = ['H', 'H2', 'OH', 'O2', 'AR', 'N2']
+        extra = ['H2O', 'CH4']
+
+        state = 500, 2e5, 'H2:0.3, O2:0.28, CH4:0.02, H2O:0.3, N2:0.1'
+
+        gas1 = ct.Solution(thermo='IdealGas', species=[S[s] for s in base+extra])
+        gas1.transport_model = 'Multi'
+        gas1.TPX = state
+
+        gas2 = ct.Solution(thermo='IdealGas', species=[S[s] for s in base])
+        gas2.transport_model = 'Multi'
+        for s in extra:
+            gas2.add_species(S[s])
+        gas2.TPX = state
+
+        self.assertNear(gas1.thermal_conductivity, gas2.thermal_conductivity)
+        self.assertArrayNear(gas1.multi_diff_coeffs, gas2.multi_diff_coeffs)
+
+
 class TestTransportGeometryFlags(utilities.CanteraTest):
     phase_data = """
 units(length="cm", time="s", quantity="mol", act_energy="cal/mol")
diff --git a/src/transport/GasTransport.cpp b/src/transport/GasTransport.cpp
index 9b2a42585..5535498a0 100644
--- a/src/transport/GasTransport.cpp
+++ b/src/transport/GasTransport.cpp
@@ -101,6 +101,11 @@ GasTransport& GasTransport::operator=(const GasTransport& right)
 
 void GasTransport::update_T()
 {
+    if (m_thermo->nSpecies() != m_nsp) {
+        // Rebuild data structures if number of species has changed
+        init(m_thermo, m_mode, m_log_level);
+    }
+
     double T = m_thermo->temperature();
     if (T == m_temp) {
         return;
@@ -517,6 +522,10 @@ void GasTransport::fitCollisionIntegrals(MMCollisionInt& integrals)
         }
     }
     vector_fp fitlist;
+    m_omega22_poly.clear();
+    m_astar_poly.clear();
+    m_bstar_poly.clear();
+    m_cstar_poly.clear();
     for (size_t i = 0; i < m_nsp; i++) {
         for (size_t j = i; j < m_nsp; j++) {
             // Chemkin fits only delta* = 0
@@ -557,6 +566,9 @@ void GasTransport::fitProperties(MMCollisionInt& integrals)
     vector_fp tlog(np), spvisc(np), spcond(np);
     vector_fp w(np), w2(np);
 
+    m_visccoeffs.clear();
+    m_condcoeffs.clear();
+
     // generate array of log(t) values
     for (size_t n = 0; n < np; n++) {
         double t = m_thermo->minTemp() + dt*n;
@@ -715,6 +727,7 @@ void GasTransport::fitProperties(MMCollisionInt& integrals)
 
     mxerr = 0.0, mxrelerr = 0.0;
     vector_fp diff(np + 1);
+    m_diffcoeffs.clear();
     for (size_t k = 0; k < m_nsp; k++) {
         for (size_t j = k; j < m_nsp; j++) {
             for (size_t n = 0; n < np; n++) {
diff --git a/src/transport/MixTransport.cpp b/src/transport/MixTransport.cpp
index c57b76433..b84aca997 100644
--- a/src/transport/MixTransport.cpp
+++ b/src/transport/MixTransport.cpp
@@ -123,7 +123,7 @@ void MixTransport::getSpeciesFluxes(size_t ndim, const doublereal* const grad_T,
 void MixTransport::update_T()
 {
     doublereal t = m_thermo->temperature();
-    if (t == m_temp) {
+    if (t == m_temp && m_nsp == m_thermo->nSpecies()) {
         return;
     }
     if (t < 0.0) {
diff --git a/src/transport/MultiTransport.cpp b/src/transport/MultiTransport.cpp
index bbc922342..7be5c765f 100644
--- a/src/transport/MultiTransport.cpp
+++ b/src/transport/MultiTransport.cpp
@@ -388,7 +388,7 @@ void MultiTransport::getMultiDiffCoeffs(const size_t ld, doublereal* const d)
 
 void MultiTransport::update_T()
 {
-    if (m_temp == m_thermo->temperature()) {
+    if (m_temp == m_thermo->temperature() && m_nsp == m_thermo->nSpecies()) {
         return;
     }
     GasTransport::update_T();