From 71c85078ae04ceb13046bd3719a3b645715bcc23 Mon Sep 17 00:00:00 2001
From: Ingmar Schoegl <ischoegl@lsu.edu>
Date: Sat, 22 Jun 2024 11:23:23 +0200
Subject: [PATCH 1/8] [oneD] Rename StFlow to Flow1D

---
 .../reference/onedim/governing-equations.md   |   6 +-
 include/cantera/clib/ctonedim.h               |  14 +-
 include/cantera/oneD/Boundary1D.h             |  10 +-
 include/cantera/oneD/{StFlow.h => Flow1D.h}   |  16 +-
 include/cantera/oneD/IonFlow.h                |   4 +-
 include/cantera/onedim.h                      |   2 +-
 interfaces/cython/cantera/_onedim.pxd         |   6 +-
 interfaces/cython/cantera/_onedim.pyx         |   2 +-
 .../matlab_experimental/OneDim/Domain1D.m     |   6 +-
 .../matlab_experimental/OneDim/Flow1D.m       |  10 +-
 samples/cxx/flamespeed/flamespeed.cpp         |   4 +-
 src/clib/ctonedim.cpp                         |  30 +--
 src/oneD/Boundary1D.cpp                       |   8 +-
 src/oneD/DomainFactory.cpp                    |  16 +-
 src/oneD/{StFlow.cpp => Flow1D.cpp}           | 180 +++++++++---------
 src/oneD/IonFlow.cpp                          |  16 +-
 src/oneD/Sim1D.cpp                            |  12 +-
 src/oneD/refine.cpp                           |   4 +-
 test/clib/test_ctonedim.cpp                   |  12 +-
 test/oneD/test_oneD.cpp                       |   8 +-
 20 files changed, 183 insertions(+), 183 deletions(-)
 rename include/cantera/oneD/{StFlow.h => Flow1D.h} (98%)
 rename src/oneD/{StFlow.cpp => Flow1D.cpp} (88%)

diff --git a/doc/sphinx/reference/onedim/governing-equations.md b/doc/sphinx/reference/onedim/governing-equations.md
index 15d6adf975c..b6a04be8345 100644
--- a/doc/sphinx/reference/onedim/governing-equations.md
+++ b/doc/sphinx/reference/onedim/governing-equations.md
@@ -7,7 +7,7 @@ solution to reduce the three-dimensional governing equations to a single dimensi
 ## Axisymmetric Stagnation Flow
 
 The governing equations for a steady axisymmetric stagnation flow follow those derived
-in Section 7.2 of {cite:t}`kee2017` and are implemented by class {ct}`StFlow`.
+in Section 7.2 of {cite:t}`kee2017` and are implemented by class {ct}`Flow1D`.
 
 *Continuity*:
 
@@ -99,7 +99,7 @@ $$
 where $D_{ki}$ is the multicomponent diffusion coefficient and $D_k^T$ is the Soret
 diffusion coefficient. Inclusion of the Soret calculation must be explicitly enabled
 when setting up the simulation, on top of specifying a multicomponent transport model,
-for example by using the {ct}`StFlow::enableSoret` method (C++) or setting the
+for example by using the {ct}`Flow1D::enableSoret` method (C++) or setting the
 {py:attr}`~cantera.FlameBase.soret_enabled` property (Python).
 
 ## Boundary Conditions
@@ -295,4 +295,4 @@ $$
 
 $$
 \dot{m}_\t{in,right} = - \rho(z_{\t{in,right}}) U_o(z_{\t{in,right}})
-$$
\ No newline at end of file
+$$
diff --git a/include/cantera/clib/ctonedim.h b/include/cantera/clib/ctonedim.h
index 438aa7df64d..e1267daf760 100644
--- a/include/cantera/clib/ctonedim.h
+++ b/include/cantera/clib/ctonedim.h
@@ -59,14 +59,14 @@ extern "C" {
 
     CANTERA_CAPI int inlet_setSpreadRate(int i, double v);
 
-    CANTERA_CAPI int stflow_new(int iph, int ikin, int itr, int itype);
-    CANTERA_CAPI int stflow_setTransport(int i, int itr);
-    CANTERA_CAPI int stflow_enableSoret(int i, int iSoret);
-    CANTERA_CAPI int stflow_setPressure(int i, double p);
-    CANTERA_CAPI double stflow_pressure(int i);
-    CANTERA_CAPI int stflow_setFixedTempProfile(int i, size_t n, const double* pos,
+    CANTERA_CAPI int flow1D_new(int iph, int ikin, int itr, int itype);
+    CANTERA_CAPI int flow1D_setTransport(int i, int itr);
+    CANTERA_CAPI int flow1D_enableSoret(int i, int iSoret);
+    CANTERA_CAPI int flow1D_setPressure(int i, double p);
+    CANTERA_CAPI double flow1D_pressure(int i);
+    CANTERA_CAPI int flow1D_setFixedTempProfile(int i, size_t n, const double* pos,
             size_t m, const double* temp);
-    CANTERA_CAPI int stflow_solveEnergyEqn(int i, int flag);
+    CANTERA_CAPI int flow1D_solveEnergyEqn(int i, int flag);
 
     CANTERA_CAPI int sim1D_new(size_t nd, const int* domains);
     CANTERA_CAPI int sim1D_del(int i);
diff --git a/include/cantera/oneD/Boundary1D.h b/include/cantera/oneD/Boundary1D.h
index da39a495baf..3c285bb3741 100644
--- a/include/cantera/oneD/Boundary1D.h
+++ b/include/cantera/oneD/Boundary1D.h
@@ -13,7 +13,7 @@
 #include "Domain1D.h"
 #include "cantera/thermo/SurfPhase.h"
 #include "cantera/kinetics/InterfaceKinetics.h"
-#include "StFlow.h"
+#include "Flow1D.h"
 
 namespace Cantera
 {
@@ -109,8 +109,8 @@ class Boundary1D : public Domain1D
 protected:
     void _init(size_t n);
 
-    StFlow* m_flow_left = nullptr;
-    StFlow* m_flow_right = nullptr;
+    Flow1D* m_flow_left = nullptr;
+    Flow1D* m_flow_right = nullptr;
     size_t m_ilr = 0;
     size_t m_left_nv = 0;
     size_t m_right_nv = 0;
@@ -176,7 +176,7 @@ class Inlet1D : public Boundary1D
     size_t m_nsp = 0;
     vector<double> m_yin;
     string m_xstr;
-    StFlow* m_flow = nullptr;
+    Flow1D* m_flow = nullptr;
 };
 
 /**
@@ -293,7 +293,7 @@ class OutletRes1D : public Boundary1D
     size_t m_nsp = 0;
     vector<double> m_yres;
     string m_xstr;
-    StFlow* m_flow = nullptr;
+    Flow1D* m_flow = nullptr;
 };
 
 /**
diff --git a/include/cantera/oneD/StFlow.h b/include/cantera/oneD/Flow1D.h
similarity index 98%
rename from include/cantera/oneD/StFlow.h
rename to include/cantera/oneD/Flow1D.h
index 8f96b13b0a6..a4bc7ccd526 100644
--- a/include/cantera/oneD/StFlow.h
+++ b/include/cantera/oneD/Flow1D.h
@@ -1,10 +1,10 @@
-//! @file StFlow.h
+//! @file Flow1D.h
 
 // This file is part of Cantera. See License.txt in the top-level directory or
 // at https://cantera.org/license.txt for license and copyright information.
 
-#ifndef CT_STFLOW_H
-#define CT_STFLOW_H
+#ifndef CT_FLOW1D_H
+#define CT_FLOW1D_H
 
 #include "Domain1D.h"
 #include "cantera/base/Array.h"
@@ -42,7 +42,7 @@ class Transport;
  *  similarity solution for chemically-reacting, axisymmetric flows.
  *  @ingroup flowGroup
  */
-class StFlow : public Domain1D
+class Flow1D : public Domain1D
 {
 public:
     //--------------------------------
@@ -54,19 +54,19 @@ class StFlow : public Domain1D
     //!     to evaluate all thermodynamic, kinetic, and transport properties.
     //! @param nsp Number of species.
     //! @param points Initial number of grid points
-    StFlow(ThermoPhase* ph = 0, size_t nsp = 1, size_t points = 1);
+    Flow1D(ThermoPhase* ph = 0, size_t nsp = 1, size_t points = 1);
 
     //! Delegating constructor
-    StFlow(shared_ptr<ThermoPhase> th, size_t nsp = 1, size_t points = 1);
+    Flow1D(shared_ptr<ThermoPhase> th, size_t nsp = 1, size_t points = 1);
 
     //! Create a new flow domain.
     //! @param sol  Solution object used to evaluate all thermodynamic, kinetic, and
     //!     transport properties
     //! @param id  name of flow domain
     //! @param points  initial number of grid points
-    StFlow(shared_ptr<Solution> sol, const string& id="", size_t points=1);
+    Flow1D(shared_ptr<Solution> sol, const string& id="", size_t points=1);
 
-    ~StFlow();
+    ~Flow1D();
 
     string domainType() const override;
 
diff --git a/include/cantera/oneD/IonFlow.h b/include/cantera/oneD/IonFlow.h
index 484bdb51080..5392d9616dc 100644
--- a/include/cantera/oneD/IonFlow.h
+++ b/include/cantera/oneD/IonFlow.h
@@ -6,7 +6,7 @@
 #ifndef CT_IONFLOW_H
 #define CT_IONFLOW_H
 
-#include "cantera/oneD/StFlow.h"
+#include "cantera/oneD/Flow1D.h"
 
 namespace Cantera
 {
@@ -25,7 +25,7 @@ namespace Cantera
  *
  * @ingroup flowGroup
  */
-class IonFlow : public StFlow
+class IonFlow : public Flow1D
 {
 public:
     IonFlow(ThermoPhase* ph = 0, size_t nsp = 1, size_t points = 1);
diff --git a/include/cantera/onedim.h b/include/cantera/onedim.h
index a74c0c3eda4..eef1579b0c8 100644
--- a/include/cantera/onedim.h
+++ b/include/cantera/onedim.h
@@ -17,7 +17,7 @@
 #include "oneD/Sim1D.h"
 #include "oneD/Domain1D.h"
 #include "oneD/Boundary1D.h"
-#include "oneD/StFlow.h"
+#include "oneD/Flow1D.h"
 #include "oneD/refine.h"
 
 #endif
diff --git a/interfaces/cython/cantera/_onedim.pxd b/interfaces/cython/cantera/_onedim.pxd
index cf2800bf30f..3863e088096 100644
--- a/interfaces/cython/cantera/_onedim.pxd
+++ b/interfaces/cython/cantera/_onedim.pxd
@@ -64,8 +64,8 @@ cdef extern from "cantera/oneD/Boundary1D.h":
         cbool coverageEnabled()
 
 
-cdef extern from "cantera/oneD/StFlow.h":
-    cdef cppclass CxxStFlow "Cantera::StFlow" (CxxDomain1D):
+cdef extern from "cantera/oneD/Flow1D.h":
+    cdef cppclass CxxFlow1D "Cantera::Flow1D" (CxxDomain1D):
         void setTransportModel(const string&) except +translate_exception
         string type()
         string transportModel()
@@ -173,7 +173,7 @@ cdef class ReactingSurface1D(Boundary1D):
     cdef public Kinetics surface
 
 cdef class FlowBase(Domain1D):
-    cdef CxxStFlow* flow
+    cdef CxxFlow1D* flow
 
 cdef class Sim1D:
     cdef CxxSim1D* sim
diff --git a/interfaces/cython/cantera/_onedim.pyx b/interfaces/cython/cantera/_onedim.pyx
index 4da961939b7..e7335ab0a1a 100644
--- a/interfaces/cython/cantera/_onedim.pyx
+++ b/interfaces/cython/cantera/_onedim.pyx
@@ -449,7 +449,7 @@ cdef class FlowBase(Domain1D):
         self._domain = CxxNewDomain1D(
             stringify(self._domain_type), phase._base, stringify(name))
         self.domain = self._domain.get()
-        self.flow = <CxxStFlow*>self.domain
+        self.flow = <CxxFlow1D*>self.domain
 
     def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
diff --git a/interfaces/matlab_experimental/OneDim/Domain1D.m b/interfaces/matlab_experimental/OneDim/Domain1D.m
index 0a992756ca5..658e8bb5a8b 100644
--- a/interfaces/matlab_experimental/OneDim/Domain1D.m
+++ b/interfaces/matlab_experimental/OneDim/Domain1D.m
@@ -108,12 +108,12 @@ function delete(d)
 
         function set.energyEnabled(d, flag)
             d.energyEnabled = flag;
-            ctFunc('stflow_solveEnergyEqn', d.domainID, int8(flag));
+            ctFunc('flow1D_solveEnergyEqn', d.domainID, int8(flag));
         end
 
         function set.soretEnabled(d, flag)
             d.soretEnabled = flag;
-            ctFunc('stflow_enableSoret', d.domainID, int8(flag));
+            ctFunc('flow1D_enableSoret', d.domainID, int8(flag));
         end
 
         %% Domain Get Methods
@@ -359,7 +359,7 @@ function setTransientTolerances(d, component, rtol, atol)
         end
 
         function set.transport(d, itr)
-            ctFunc('stflow_setTransport', d.domainID, itr);
+            ctFunc('flow1D_setTransport', d.domainID, itr);
         end
 
         function setupGrid(d, grid)
diff --git a/interfaces/matlab_experimental/OneDim/Flow1D.m b/interfaces/matlab_experimental/OneDim/Flow1D.m
index 253699ed815..1da4c7fc1fd 100644
--- a/interfaces/matlab_experimental/OneDim/Flow1D.m
+++ b/interfaces/matlab_experimental/OneDim/Flow1D.m
@@ -31,11 +31,11 @@
         %% Flow1D Class Methods
 
         function pressure = get.P(d)
-            pressure = ctFunc('stflow_pressure', d.domainID);
+            pressure = ctFunc('flow1D_pressure', d.domainID);
         end
 
         function set.P(d, p)
-            ctFunc('stflow_setPressure', d.domainID, p);
+            ctFunc('flow1D_setPressure', d.domainID, p);
         end
 
         function setFixedTempProfile(d, profile)
@@ -57,11 +57,11 @@ function setFixedTempProfile(d, profile)
 
             if sz(1) == 2
                 l = length(profile(1, :));
-                ctFunc('stflow_setFixedTempProfile', d.domainID, ...
+                ctFunc('flow1D_setFixedTempProfile', d.domainID, ...
                         l, profile(1, :), l, profile(2, :));
             elseif sz(2) == 2
                 l = length(profile(:, 1));
-                ctFunc('stflow_setFixedTempProfile', d.domainID, ...
+                ctFunc('flow1D_setFixedTempProfile', d.domainID, ...
                         l, profile(:, 1), l, profile(:, 2));
             else
                 error('Wrong temperature profile array shape.');
@@ -71,4 +71,4 @@ function setFixedTempProfile(d, profile)
 
     end
 
-end
\ No newline at end of file
+end
diff --git a/samples/cxx/flamespeed/flamespeed.cpp b/samples/cxx/flamespeed/flamespeed.cpp
index 09e7749ee05..2859038d613 100644
--- a/samples/cxx/flamespeed/flamespeed.cpp
+++ b/samples/cxx/flamespeed/flamespeed.cpp
@@ -8,7 +8,7 @@
  *
  *     flamespeed [equivalence_ratio] [refine_grid] [loglevel]
  *
- * Requires: cantera >= 3.0
+ * Requires: cantera >= 3.1
  *
  * .. tags:: C++, combustion, 1D flow, premixed flame, flame speed, saving output
  */
@@ -56,7 +56,7 @@ int flamespeed(double phi, bool refine_grid, int loglevel)
 
         //-------- step 1: create the flow -------------
 
-        auto flow = newDomain<StFlow>("gas-flow", sol, "flow");
+        auto flow = newDomain<Flow1D>("gas-flow", sol, "flow");
         flow->setFreeFlow();
 
         // create an initial grid
diff --git a/src/clib/ctonedim.cpp b/src/clib/ctonedim.cpp
index 5bc04656108..355bd15080f 100644
--- a/src/clib/ctonedim.cpp
+++ b/src/clib/ctonedim.cpp
@@ -400,11 +400,11 @@ extern "C" {
 
     //------------------ stagnation flow domains --------------------
 
-    int stflow_new(int iph, int ikin, int itr, int itype)
+    int flow1D_new(int iph, int ikin, int itr, int itype)
     {
         try {
             auto ph = ThermoCabinet::at(iph);
-            auto x = make_shared<StFlow>(ph, ph->nSpecies(), 2);
+            auto x = make_shared<Flow1D>(ph, ph->nSpecies(), 2);
             if (itype == 1) {
                 x->setAxisymmetricFlow();
             } else if (itype == 2) {
@@ -420,47 +420,47 @@ extern "C" {
         }
     }
 
-    int stflow_setTransport(int i, int itr)
+    int flow1D_setTransport(int i, int itr)
     {
         try {
-            DomainCabinet::get<StFlow>(i).setTransport(TransportCabinet::at(itr));
+            DomainCabinet::get<Flow1D>(i).setTransport(TransportCabinet::at(itr));
             return 0;
         } catch (...) {
             return handleAllExceptions(-1, ERR);
         }
     }
 
-    int stflow_enableSoret(int i, int iSoret)
+    int flow1D_enableSoret(int i, int iSoret)
     {
         try {
             bool withSoret = (iSoret > 0);
-            DomainCabinet::get<StFlow>(i).enableSoret(withSoret);
+            DomainCabinet::get<Flow1D>(i).enableSoret(withSoret);
             return 0;
         } catch (...) {
             return handleAllExceptions(-1, ERR);
         }
     }
 
-    int stflow_setPressure(int i, double p)
+    int flow1D_setPressure(int i, double p)
     {
         try {
-            DomainCabinet::get<StFlow>(i).setPressure(p);
+            DomainCabinet::get<Flow1D>(i).setPressure(p);
             return 0;
         } catch (...) {
             return handleAllExceptions(-1, ERR);
         }
     }
 
-    double stflow_pressure(int i)
+    double flow1D_pressure(int i)
     {
         try {
-            return DomainCabinet::get<StFlow>(i).pressure();
+            return DomainCabinet::get<Flow1D>(i).pressure();
         } catch (...) {
             return handleAllExceptions(DERR, DERR);
         }
     }
 
-    int stflow_setFixedTempProfile(int i, size_t n, const double* pos,
+    int flow1D_setFixedTempProfile(int i, size_t n, const double* pos,
                                    size_t m, const double* temp)
     {
         try {
@@ -469,20 +469,20 @@ extern "C" {
                 vpos[j] = pos[j];
                 vtemp[j] = temp[j];
             }
-            DomainCabinet::get<StFlow>(i).setFixedTempProfile(vpos, vtemp);
+            DomainCabinet::get<Flow1D>(i).setFixedTempProfile(vpos, vtemp);
             return 0;
         } catch (...) {
             return handleAllExceptions(-1, ERR);
         }
     }
 
-    int stflow_solveEnergyEqn(int i, int flag)
+    int flow1D_solveEnergyEqn(int i, int flag)
     {
         try {
             if (flag > 0) {
-                DomainCabinet::get<StFlow>(i).solveEnergyEqn(npos);
+                DomainCabinet::get<Flow1D>(i).solveEnergyEqn(npos);
             } else {
-                DomainCabinet::get<StFlow>(i).fixTemperature(npos);
+                DomainCabinet::get<Flow1D>(i).fixTemperature(npos);
             }
             return 0;
         } catch (...) {
diff --git a/src/oneD/Boundary1D.cpp b/src/oneD/Boundary1D.cpp
index 5b8cc7d8185..f698e88d14e 100644
--- a/src/oneD/Boundary1D.cpp
+++ b/src/oneD/Boundary1D.cpp
@@ -6,7 +6,7 @@
 #include "cantera/base/SolutionArray.h"
 #include "cantera/oneD/Boundary1D.h"
 #include "cantera/oneD/OneDim.h"
-#include "cantera/oneD/StFlow.h"
+#include "cantera/oneD/Flow1D.h"
 
 using namespace std;
 
@@ -42,7 +42,7 @@ void Boundary1D::_init(size_t n)
             }
             m_left_loc = container().start(m_index-1);
             m_left_points = r.nPoints();
-            m_flow_left = dynamic_cast<StFlow*>(&r);
+            m_flow_left = dynamic_cast<Flow1D*>(&r);
             if (m_flow_left != nullptr) {
                 m_phase_left = &m_flow_left->phase();
             }
@@ -64,7 +64,7 @@ void Boundary1D::_init(size_t n)
                 m_right_nsp = 0;
             }
             m_right_loc = container().start(m_index+1);
-            m_flow_right = dynamic_cast<StFlow*>(&r);
+            m_flow_right = dynamic_cast<Flow1D*>(&r);
             if (m_flow_right != nullptr) {
                 m_phase_right = &m_flow_right->phase();
             }
@@ -212,7 +212,7 @@ void Inlet1D::eval(size_t jg, double* xg, double* rg,
                 // When using two-point control, the mass flow rate at the left inlet is
                 // not specified. Instead, the mass flow rate is dictated by the
                 // velocity at the left inlet, which comes from the U variable. The
-                // default boundary condition specified in the StFlow.cpp file already
+                // default boundary condition specified in the Flow1D.cpp file already
                 // handles this case. We only need to update the stored value of m_mdot
                 // so that other equations that use the quantity are consistent.
                 m_mdot = m_flow->density(0)*xb[c_offset_U];
diff --git a/src/oneD/DomainFactory.cpp b/src/oneD/DomainFactory.cpp
index e11957e2259..593c4ddc05d 100644
--- a/src/oneD/DomainFactory.cpp
+++ b/src/oneD/DomainFactory.cpp
@@ -5,7 +5,7 @@
 
 #include "cantera/oneD/DomainFactory.h"
 #include "cantera/oneD/Boundary1D.h"
-#include "cantera/oneD/StFlow.h"
+#include "cantera/oneD/Flow1D.h"
 #include "cantera/oneD/IonFlow.h"
 #include "cantera/transport/Transport.h"
 
@@ -39,37 +39,37 @@ DomainFactory::DomainFactory()
         return new ReactingSurf1D(solution, id);
     });
     reg("gas-flow", [](shared_ptr<Solution> solution, const string& id) {
-        return new StFlow(solution, id);
+        return new Flow1D(solution, id);
     });
     reg("ion-flow", [](shared_ptr<Solution> solution, const string& id) {
         return new IonFlow(solution, id);
     });
     reg("free-flow", [](shared_ptr<Solution> solution, const string& id) {
-        StFlow* ret;
+        Flow1D* ret;
         if (solution->transport()->transportModel() == "ionized-gas") {
             ret = new IonFlow(solution, id);
         } else {
-            ret = new StFlow(solution, id);
+            ret = new Flow1D(solution, id);
         }
         ret->setFreeFlow();
         return ret;
     });
     reg("axisymmetric-flow", [](shared_ptr<Solution> solution, const string& id) {
-        StFlow* ret;
+        Flow1D* ret;
         if (solution->transport()->transportModel() == "ionized-gas") {
             ret = new IonFlow(solution, id);
         } else {
-            ret = new StFlow(solution, id);
+            ret = new Flow1D(solution, id);
         }
         ret->setAxisymmetricFlow();
         return ret;
     });
     reg("unstrained-flow", [](shared_ptr<Solution> solution, const string& id) {
-        StFlow* ret;
+        Flow1D* ret;
         if (solution->transport()->transportModel() == "ionized-gas") {
             ret = new IonFlow(solution, id);
         } else {
-            ret = new StFlow(solution, id);
+            ret = new Flow1D(solution, id);
         }
         ret->setUnstrainedFlow();
         return ret;
diff --git a/src/oneD/StFlow.cpp b/src/oneD/Flow1D.cpp
similarity index 88%
rename from src/oneD/StFlow.cpp
rename to src/oneD/Flow1D.cpp
index 0ab9cd50e44..c28d5a41e5a 100644
--- a/src/oneD/StFlow.cpp
+++ b/src/oneD/Flow1D.cpp
@@ -1,10 +1,10 @@
-//! @file StFlow.cpp
+//! @file Flow1D.cpp
 
 // This file is part of Cantera. See License.txt in the top-level directory or
 // at https://cantera.org/license.txt for license and copyright information.
 
 #include "cantera/base/SolutionArray.h"
-#include "cantera/oneD/StFlow.h"
+#include "cantera/oneD/Flow1D.h"
 #include "cantera/oneD/refine.h"
 #include "cantera/transport/Transport.h"
 #include "cantera/transport/TransportFactory.h"
@@ -16,7 +16,7 @@ using namespace std;
 namespace Cantera
 {
 
-StFlow::StFlow(ThermoPhase* ph, size_t nsp, size_t points) :
+Flow1D::Flow1D(ThermoPhase* ph, size_t nsp, size_t points) :
     Domain1D(nsp+c_offset_Y, points),
     m_nsp(nsp)
 {
@@ -90,23 +90,23 @@ StFlow::StFlow(ThermoPhase* ph, size_t nsp, size_t points) :
     m_kRadiating[1] = m_thermo->speciesIndex("H2O");
 }
 
-StFlow::StFlow(shared_ptr<ThermoPhase> th, size_t nsp, size_t points)
-    : StFlow(th.get(), nsp, points)
+Flow1D::Flow1D(shared_ptr<ThermoPhase> th, size_t nsp, size_t points)
+    : Flow1D(th.get(), nsp, points)
 {
     auto sol = Solution::create();
     sol->setThermo(th);
     setSolution(sol);
 }
 
-StFlow::StFlow(shared_ptr<Solution> sol, const string& id, size_t points)
-    : StFlow(sol->thermo().get(), sol->thermo()->nSpecies(), points)
+Flow1D::Flow1D(shared_ptr<Solution> sol, const string& id, size_t points)
+    : Flow1D(sol->thermo().get(), sol->thermo()->nSpecies(), points)
 {
     setSolution(sol);
     m_id = id;
     m_kin = m_solution->kinetics().get();
     m_trans = m_solution->transport().get();
     if (m_trans->transportModel() == "none") {
-        throw CanteraError("StFlow::StFlow",
+        throw CanteraError("Flow1D::Flow1D",
             "An appropriate transport model\nshould be set when instantiating the "
             "Solution ('gas') object.");
     }
@@ -116,14 +116,14 @@ StFlow::StFlow(shared_ptr<Solution> sol, const string& id, size_t points)
     });
 }
 
-StFlow::~StFlow()
+Flow1D::~Flow1D()
 {
     if (m_solution) {
         m_solution->removeChangedCallback(this);
     }
 }
 
-string StFlow::domainType() const {
+string Flow1D::domainType() const {
     if (m_isFree) {
         return "free-flow";
     }
@@ -133,20 +133,20 @@ string StFlow::domainType() const {
     return "unstrained-flow";
 }
 
-void StFlow::setKinetics(shared_ptr<Kinetics> kin)
+void Flow1D::setKinetics(shared_ptr<Kinetics> kin)
 {
     m_kin = kin.get();
     m_solution->setKinetics(kin);
 }
 
-void StFlow::setTransport(shared_ptr<Transport> trans)
+void Flow1D::setTransport(shared_ptr<Transport> trans)
 {
     if (!trans) {
-        throw CanteraError("StFlow::setTransport", "Unable to set empty transport.");
+        throw CanteraError("Flow1D::setTransport", "Unable to set empty transport.");
     }
     m_trans = trans.get();
     if (m_trans->transportModel() == "none") {
-        throw CanteraError("StFlow::setTransport", "Invalid Transport model 'none'.");
+        throw CanteraError("Flow1D::setTransport", "Invalid Transport model 'none'.");
     }
     m_do_multicomponent = (m_trans->transportModel() == "multicomponent" ||
         m_trans->transportModel() == "multicomponent-CK");
@@ -159,7 +159,7 @@ void StFlow::setTransport(shared_ptr<Transport> trans)
     m_solution->setTransport(trans);
 }
 
-void StFlow::resize(size_t ncomponents, size_t points)
+void Flow1D::resize(size_t ncomponents, size_t points)
 {
     Domain1D::resize(ncomponents, points);
     m_rho.resize(m_points, 0.0);
@@ -185,14 +185,14 @@ void StFlow::resize(size_t ncomponents, size_t points)
     m_z.resize(m_points);
 }
 
-void StFlow::setupGrid(size_t n, const double* z)
+void Flow1D::setupGrid(size_t n, const double* z)
 {
     resize(m_nv, n);
 
     m_z[0] = z[0];
     for (size_t j = 1; j < m_points; j++) {
         if (z[j] <= z[j-1]) {
-            throw CanteraError("StFlow::setupGrid",
+            throw CanteraError("Flow1D::setupGrid",
                                "grid points must be monotonically increasing");
         }
         m_z[j] = z[j];
@@ -200,7 +200,7 @@ void StFlow::setupGrid(size_t n, const double* z)
     }
 }
 
-void StFlow::resetBadValues(double* xg)
+void Flow1D::resetBadValues(double* xg)
 {
     double* x = xg + loc();
     for (size_t j = 0; j < m_points; j++) {
@@ -210,26 +210,26 @@ void StFlow::resetBadValues(double* xg)
     }
 }
 
-void StFlow::setTransportModel(const string& trans)
+void Flow1D::setTransportModel(const string& trans)
 {
     m_solution->setTransportModel(trans);
 }
 
-string StFlow::transportModel() const {
+string Flow1D::transportModel() const {
     return m_trans->transportModel();
 }
 
-void StFlow::setFluxGradientBasis(ThermoBasis fluxGradientBasis) {
+void Flow1D::setFluxGradientBasis(ThermoBasis fluxGradientBasis) {
     m_fluxGradientBasis = fluxGradientBasis;
     if (transportModel() != "mixture-averaged-CK" &&
         transportModel() != "mixture-averaged") {
-        warn_user("StFlow::setFluxGradientBasis",
+        warn_user("Flow1D::setFluxGradientBasis",
                   "Setting fluxGradientBasis only affects "
                   "the mixture-averaged diffusion model.");
     }
 }
 
-void StFlow::_getInitialSoln(double* x)
+void Flow1D::_getInitialSoln(double* x)
 {
     for (size_t j = 0; j < m_points; j++) {
         T(x,j) = m_thermo->temperature();
@@ -238,7 +238,7 @@ void StFlow::_getInitialSoln(double* x)
     }
 }
 
-void StFlow::setGas(const double* x, size_t j)
+void Flow1D::setGas(const double* x, size_t j)
 {
     m_thermo->setTemperature(T(x,j));
     const double* yy = x + m_nv*j + c_offset_Y;
@@ -246,7 +246,7 @@ void StFlow::setGas(const double* x, size_t j)
     m_thermo->setPressure(m_press);
 }
 
-void StFlow::setGasAtMidpoint(const double* x, size_t j)
+void Flow1D::setGasAtMidpoint(const double* x, size_t j)
 {
     m_thermo->setTemperature(0.5*(T(x,j)+T(x,j+1)));
     const double* yyj = x + m_nv*j + c_offset_Y;
@@ -258,10 +258,10 @@ void StFlow::setGasAtMidpoint(const double* x, size_t j)
     m_thermo->setPressure(m_press);
 }
 
-void StFlow::_finalize(const double* x)
+void Flow1D::_finalize(const double* x)
 {
     if (!m_do_multicomponent && m_do_soret) {
-        throw CanteraError("StFlow::_finalize",
+        throw CanteraError("Flow1D::_finalize",
             "Thermal diffusion (the Soret effect) is enabled, and requires "
             "using a multicomponent transport model.");
     }
@@ -305,7 +305,7 @@ void StFlow::_finalize(const double* x)
     }
 }
 
-void StFlow::eval(size_t jGlobal, double* xGlobal, double* rsdGlobal,
+void Flow1D::eval(size_t jGlobal, double* xGlobal, double* rsdGlobal,
                   integer* diagGlobal, double rdt)
 {
     // If evaluating a Jacobian, and the global point is outside the domain of
@@ -344,7 +344,7 @@ void StFlow::eval(size_t jGlobal, double* xGlobal, double* rsdGlobal,
     evalSpecies(x, rsd, diag, rdt, jmin, jmax);
 }
 
-void StFlow::updateProperties(size_t jg, double* x, size_t jmin, size_t jmax)
+void Flow1D::updateProperties(size_t jg, double* x, size_t jmin, size_t jmax)
 {
     // properties are computed for grid points from j0 to j1
     size_t j0 = std::max<size_t>(jmin, 1) - 1;
@@ -368,7 +368,7 @@ void StFlow::updateProperties(size_t jg, double* x, size_t jmin, size_t jmax)
     updateDiffFluxes(x, j0, j1);
 }
 
-void StFlow::computeRadiation(double* x, size_t jmin, size_t jmax)
+void Flow1D::computeRadiation(double* x, size_t jmin, size_t jmax)
 {
     // Variable definitions for the Planck absorption coefficient and the
     // radiation calculation:
@@ -415,7 +415,7 @@ void StFlow::computeRadiation(double* x, size_t jmin, size_t jmax)
     }
 }
 
-void StFlow::evalContinuity(double* x, double* rsd, int* diag,
+void Flow1D::evalContinuity(double* x, double* rsd, int* diag,
                             double rdt, size_t jmin, size_t jmax)
 {
     // The left boundary has the same form for all cases.
@@ -474,7 +474,7 @@ void StFlow::evalContinuity(double* x, double* rsd, int* diag,
     }
 }
 
-void StFlow::evalMomentum(double* x, double* rsd, int* diag,
+void Flow1D::evalMomentum(double* x, double* rsd, int* diag,
                           double rdt, size_t jmin, size_t jmax)
 {
     if (!m_usesLambda) { //disable this equation
@@ -506,7 +506,7 @@ void StFlow::evalMomentum(double* x, double* rsd, int* diag,
     }
 }
 
-void StFlow::evalLambda(double* x, double* rsd, int* diag,
+void Flow1D::evalLambda(double* x, double* rsd, int* diag,
                         double rdt, size_t jmin, size_t jmax)
 {
     if (!m_usesLambda) { // disable this equation
@@ -549,7 +549,7 @@ void StFlow::evalLambda(double* x, double* rsd, int* diag,
     }
 }
 
-void StFlow::evalEnergy(double* x, double* rsd, int* diag,
+void Flow1D::evalEnergy(double* x, double* rsd, int* diag,
                         double rdt, size_t jmin, size_t jmax)
 {
     if (jmin == 0) { // left boundary
@@ -587,7 +587,7 @@ void StFlow::evalEnergy(double* x, double* rsd, int* diag,
     }
 }
 
-void StFlow::evalUo(double* x, double* rsd, int* diag,
+void Flow1D::evalUo(double* x, double* rsd, int* diag,
                     double rdt, size_t jmin, size_t jmax)
 {
     if (!m_twoPointControl) { // disable this equation
@@ -629,7 +629,7 @@ void StFlow::evalUo(double* x, double* rsd, int* diag,
     }
 }
 
-void StFlow::evalSpecies(double* x, double* rsd, int* diag,
+void Flow1D::evalSpecies(double* x, double* rsd, int* diag,
                          double rdt, size_t jmin, size_t jmax)
 {
     if (jmin == 0) { // left boundary
@@ -668,7 +668,7 @@ void StFlow::evalSpecies(double* x, double* rsd, int* diag,
     }
 }
 
-void StFlow::evalElectricField(double* x, double* rsd, int* diag,
+void Flow1D::evalElectricField(double* x, double* rsd, int* diag,
                                double rdt, size_t jmin, size_t jmax)
 {
     for (size_t j = jmin; j <= jmax; j++) {
@@ -677,7 +677,7 @@ void StFlow::evalElectricField(double* x, double* rsd, int* diag,
     }
 }
 
-void StFlow::updateTransport(double* x, size_t j0, size_t j1)
+void Flow1D::updateTransport(double* x, size_t j0, size_t j1)
 {
      if (m_do_multicomponent) {
         for (size_t j = j0; j < j1; j++) {
@@ -725,7 +725,7 @@ void StFlow::updateTransport(double* x, size_t j0, size_t j1)
     }
 }
 
-void StFlow::show(const double* x)
+void Flow1D::show(const double* x)
 {
     writelog("    Pressure:  {:10.4g} Pa\n", m_press);
 
@@ -742,7 +742,7 @@ void StFlow::show(const double* x)
     }
 }
 
-void StFlow::updateDiffFluxes(const double* x, size_t j0, size_t j1)
+void Flow1D::updateDiffFluxes(const double* x, size_t j0, size_t j1)
 {
     if (m_do_multicomponent) {
         for (size_t j = j0; j < j1; j++) {
@@ -788,7 +788,7 @@ void StFlow::updateDiffFluxes(const double* x, size_t j0, size_t j1)
     }
 }
 
-string StFlow::componentName(size_t n) const
+string Flow1D::componentName(size_t n) const
 {
     switch (n) {
     case c_offset_U:
@@ -812,7 +812,7 @@ string StFlow::componentName(size_t n) const
     }
 }
 
-size_t StFlow::componentIndex(const string& name) const
+size_t Flow1D::componentIndex(const string& name) const
 {
     if (name=="velocity") {
         return c_offset_U;
@@ -832,12 +832,12 @@ size_t StFlow::componentIndex(const string& name) const
                 return n;
             }
         }
-        throw CanteraError("StFlow1D::componentIndex",
+        throw CanteraError("Flow1D1D::componentIndex",
                            "no component named " + name);
     }
 }
 
-bool StFlow::componentActive(size_t n) const
+bool Flow1D::componentActive(size_t n) const
 {
     switch (n) {
     case c_offset_V: // spread_rate
@@ -853,7 +853,7 @@ bool StFlow::componentActive(size_t n) const
     }
 }
 
-AnyMap StFlow::getMeta() const
+AnyMap Flow1D::getMeta() const
 {
     AnyMap state = Domain1D::getMeta();
     state["transport-model"] = m_trans->transportModel();
@@ -913,7 +913,7 @@ AnyMap StFlow::getMeta() const
     return state;
 }
 
-shared_ptr<SolutionArray> StFlow::asArray(const double* soln) const
+shared_ptr<SolutionArray> Flow1D::asArray(const double* soln) const
 {
     auto arr = SolutionArray::create(
         m_solution, static_cast<int>(nPoints()), getMeta());
@@ -947,7 +947,7 @@ shared_ptr<SolutionArray> StFlow::asArray(const double* soln) const
     return arr;
 }
 
-void StFlow::fromArray(SolutionArray& arr, double* soln)
+void Flow1D::fromArray(SolutionArray& arr, double* soln)
 {
     Domain1D::setMeta(arr.meta());
     arr.setLoc(0);
@@ -968,7 +968,7 @@ void StFlow::fromArray(SolutionArray& arr, double* soln)
                 soln[index(i,j)] = data[j];
             }
         } else {
-            warn_user("StFlow::fromArray", "Saved state does not contain values for "
+            warn_user("Flow1D::fromArray", "Saved state does not contain values for "
                 "component '{}' in domain '{}'.", name, id());
         }
     }
@@ -977,7 +977,7 @@ void StFlow::fromArray(SolutionArray& arr, double* soln)
     setMeta(arr.meta());
 }
 
-void StFlow::setMeta(const AnyMap& state)
+void Flow1D::setMeta(const AnyMap& state)
 {
     if (state.hasKey("energy-enabled")) {
         const AnyValue& ee = state["energy-enabled"];
@@ -1047,13 +1047,13 @@ void StFlow::setMeta(const AnyMap& state)
             m_tLeft = cm["left-temperature"].asDouble();
             m_tRight = cm["right-temperature"].asDouble();
         } else {
-            warn_user("StFlow::setMeta", "Unknown continuation method '{}'.",
+            warn_user("Flow1D::setMeta", "Unknown continuation method '{}'.",
                 cm["type"].asString());
         }
     }
 }
 
-void StFlow::solveEnergyEqn(size_t j)
+void Flow1D::solveEnergyEqn(size_t j)
 {
     bool changed = false;
     if (j == npos) {
@@ -1077,43 +1077,43 @@ void StFlow::solveEnergyEqn(size_t j)
     }
 }
 
-size_t StFlow::getSolvingStage() const
+size_t Flow1D::getSolvingStage() const
 {
-    throw NotImplementedError("StFlow::getSolvingStage",
+    throw NotImplementedError("Flow1D::getSolvingStage",
         "Not used by '{}' objects.", type());
 }
 
-void StFlow::setSolvingStage(const size_t stage)
+void Flow1D::setSolvingStage(const size_t stage)
 {
-    throw NotImplementedError("StFlow::setSolvingStage",
+    throw NotImplementedError("Flow1D::setSolvingStage",
         "Not used by '{}' objects.", type());
 }
 
-void StFlow::solveElectricField(size_t j)
+void Flow1D::solveElectricField(size_t j)
 {
-    throw NotImplementedError("StFlow::solveElectricField",
+    throw NotImplementedError("Flow1D::solveElectricField",
         "Not used by '{}' objects.", type());
 }
 
-void StFlow::fixElectricField(size_t j)
+void Flow1D::fixElectricField(size_t j)
 {
-    throw NotImplementedError("StFlow::fixElectricField",
+    throw NotImplementedError("Flow1D::fixElectricField",
         "Not used by '{}' objects.", type());
 }
 
-bool StFlow::doElectricField(size_t j) const
+bool Flow1D::doElectricField(size_t j) const
 {
-    throw NotImplementedError("StFlow::doElectricField",
+    throw NotImplementedError("Flow1D::doElectricField",
         "Not used by '{}' objects.", type());
 }
 
-void StFlow::setBoundaryEmissivities(double e_left, double e_right)
+void Flow1D::setBoundaryEmissivities(double e_left, double e_right)
 {
     if (e_left < 0 || e_left > 1) {
-        throw CanteraError("StFlow::setBoundaryEmissivities",
+        throw CanteraError("Flow1D::setBoundaryEmissivities",
             "The left boundary emissivity must be between 0.0 and 1.0!");
     } else if (e_right < 0 || e_right > 1) {
-        throw CanteraError("StFlow::setBoundaryEmissivities",
+        throw CanteraError("Flow1D::setBoundaryEmissivities",
             "The right boundary emissivity must be between 0.0 and 1.0!");
     } else {
         m_epsilon_left = e_left;
@@ -1121,7 +1121,7 @@ void StFlow::setBoundaryEmissivities(double e_left, double e_right)
     }
 }
 
-void StFlow::fixTemperature(size_t j)
+void Flow1D::fixTemperature(size_t j)
 {
     bool changed = false;
     if (j == npos) {
@@ -1145,7 +1145,7 @@ void StFlow::fixTemperature(size_t j)
     }
 }
 
-void StFlow::grad_hk(const double* x, size_t j)
+void Flow1D::grad_hk(const double* x, size_t j)
 {
     for(size_t k = 0; k < m_nsp; k++) {
         if (u(x, j) > 0.0) {
@@ -1158,122 +1158,122 @@ void StFlow::grad_hk(const double* x, size_t j)
 }
 
 // Two-point control functions
-double StFlow::leftControlPointTemperature() const
+double Flow1D::leftControlPointTemperature() const
 {
     if (m_twoPointControl) {
         if (m_zLeft != Undef) {
             return m_tLeft;
         } else {
-            throw CanteraError("StFlow::leftControlPointTemperature",
+            throw CanteraError("Flow1D::leftControlPointTemperature",
                 "Invalid operation: left control point location is not set");
         }
     } else {
-        throw CanteraError("StFlow::leftControlPointTemperature",
+        throw CanteraError("Flow1D::leftControlPointTemperature",
                 "Invalid operation: two-point control is not enabled.");
     }
 }
 
-double StFlow::leftControlPointCoordinate() const
+double Flow1D::leftControlPointCoordinate() const
 {
     if (m_twoPointControl) {
         if (m_zLeft != Undef) {
             return m_zLeft;
         } else {
-            throw CanteraError("StFlow::leftControlPointCoordinate",
+            throw CanteraError("Flow1D::leftControlPointCoordinate",
                 "Invalid operation: left control point location is not set");
         }
     } else {
-        throw CanteraError("StFlow::leftControlPointCoordinate",
+        throw CanteraError("Flow1D::leftControlPointCoordinate",
                 "Invalid operation: two-point control is not enabled.");
     }
 }
 
-void StFlow::setLeftControlPointTemperature(double temperature)
+void Flow1D::setLeftControlPointTemperature(double temperature)
 {
     if (m_twoPointControl) {
         if (m_zLeft != Undef) {
             m_tLeft = temperature;
         } else {
-            throw CanteraError("StFlow::setLeftControlPointTemperature",
+            throw CanteraError("Flow1D::setLeftControlPointTemperature",
                 "Invalid operation: left control point location is not set");
         }
     } else {
-        throw CanteraError("StFlow::setLeftControlPointTemperature",
+        throw CanteraError("Flow1D::setLeftControlPointTemperature",
                 "Invalid operation: two-point control is not enabled.");
     }
 }
 
-void StFlow::setLeftControlPointCoordinate(double z_left)
+void Flow1D::setLeftControlPointCoordinate(double z_left)
 {
     if (m_twoPointControl) {
             m_zLeft = z_left;
     } else {
-        throw CanteraError("StFlow::setLeftControlPointCoordinate",
+        throw CanteraError("Flow1D::setLeftControlPointCoordinate",
                 "Invalid operation: two-point control is not enabled.");
     }
 }
 
-double StFlow::rightControlPointTemperature() const
+double Flow1D::rightControlPointTemperature() const
 {
     if (m_twoPointControl) {
         if (m_zRight != Undef) {
             return m_tRight;
         } else {
-            throw CanteraError("StFlow::rightControlPointTemperature",
+            throw CanteraError("Flow1D::rightControlPointTemperature",
                 "Invalid operation: right control point location is not set");
         }
     } else {
-        throw CanteraError("StFlow::rightControlPointTemperature",
+        throw CanteraError("Flow1D::rightControlPointTemperature",
                 "Invalid operation: two-point control is not enabled.");
     }
 }
 
-double StFlow::rightControlPointCoordinate() const
+double Flow1D::rightControlPointCoordinate() const
 {
     if (m_twoPointControl) {
         if (m_zRight != Undef) {
             return m_zRight;
         } else {
-            throw CanteraError("StFlow::rightControlPointCoordinate",
+            throw CanteraError("Flow1D::rightControlPointCoordinate",
                 "Invalid operation: right control point location is not set");
         }
     } else {
-        throw CanteraError("StFlow::rightControlPointCoordinate",
+        throw CanteraError("Flow1D::rightControlPointCoordinate",
                 "Invalid operation: two-point control is not enabled.");
     }
 }
 
-void StFlow::setRightControlPointTemperature(double temperature)
+void Flow1D::setRightControlPointTemperature(double temperature)
 {
     if (m_twoPointControl) {
         if (m_zRight != Undef) {
             m_tRight = temperature;
         } else {
-            throw CanteraError("StFlow::setRightControlPointTemperature",
+            throw CanteraError("Flow1D::setRightControlPointTemperature",
                 "Invalid operation: right control point location is not set");
         }
     } else {
-        throw CanteraError("StFlow::setRightControlPointTemperature",
+        throw CanteraError("Flow1D::setRightControlPointTemperature",
                 "Invalid operation: two-point control is not enabled.");
     }
 }
 
-void StFlow::setRightControlPointCoordinate(double z_right)
+void Flow1D::setRightControlPointCoordinate(double z_right)
 {
     if (m_twoPointControl) {
             m_zRight = z_right;
     } else {
-        throw CanteraError("StFlow::setRightControlPointCoordinate",
+        throw CanteraError("Flow1D::setRightControlPointCoordinate",
                 "Invalid operation: two-point control is not enabled.");
     }
 }
 
-void StFlow::enableTwoPointControl(bool twoPointControl)
+void Flow1D::enableTwoPointControl(bool twoPointControl)
 {
     if (isStrained()) {
         m_twoPointControl = twoPointControl;
     } else {
-        throw CanteraError("StFlow::enableTwoPointControl",
+        throw CanteraError("Flow1D::enableTwoPointControl",
             "Invalid operation: two-point control can only be used"
             "with axisymmetric flames.");
     }
diff --git a/src/oneD/IonFlow.cpp b/src/oneD/IonFlow.cpp
index c1363266d14..9630bd6f19c 100644
--- a/src/oneD/IonFlow.cpp
+++ b/src/oneD/IonFlow.cpp
@@ -4,7 +4,7 @@
 // at https://cantera.org/license.txt for license and copyright information.
 
 #include "cantera/oneD/IonFlow.h"
-#include "cantera/oneD/StFlow.h"
+#include "cantera/oneD/Flow1D.h"
 #include "cantera/oneD/refine.h"
 #include "cantera/transport/Transport.h"
 #include "cantera/numerics/funcs.h"
@@ -16,7 +16,7 @@ namespace Cantera
 {
 
 IonFlow::IonFlow(ThermoPhase* ph, size_t nsp, size_t points) :
-    StFlow(ph, nsp, points)
+    Flow1D(ph, nsp, points)
 {
     // make a local copy of species charge
     for (size_t k = 0; k < m_nsp; k++) {
@@ -82,7 +82,7 @@ string IonFlow::domainType() const {
 }
 
 void IonFlow::resize(size_t components, size_t points){
-    StFlow::resize(components, points);
+    Flow1D::resize(components, points);
     m_mobility.resize(m_nsp*m_points);
     m_do_species.resize(m_nsp,true);
     m_do_electric_field.resize(m_points,false);
@@ -93,13 +93,13 @@ bool IonFlow::componentActive(size_t n) const
     if (n == c_offset_E) {
         return true;
     } else {
-        return StFlow::componentActive(n);
+        return Flow1D::componentActive(n);
     }
 }
 
 void IonFlow::updateTransport(double* x, size_t j0, size_t j1)
 {
-    StFlow::updateTransport(x,j0,j1);
+    Flow1D::updateTransport(x,j0,j1);
     for (size_t j = j0; j < j1; j++) {
         setGasAtMidpoint(x,j);
         m_trans->getMobilities(&m_mobility[j*m_nsp]);
@@ -202,7 +202,7 @@ void IonFlow::setSolvingStage(const size_t stage)
 void IonFlow::evalElectricField(double* x, double* rsd, int* diag,
                                 double rdt, size_t jmin, size_t jmax)
 {
-    StFlow::evalElectricField(x, rsd, diag, rdt, jmin, jmax);
+    Flow1D::evalElectricField(x, rsd, diag, rdt, jmin, jmax);
     if (m_stage != 2) {
         return;
     }
@@ -227,7 +227,7 @@ void IonFlow::evalElectricField(double* x, double* rsd, int* diag,
 void IonFlow::evalSpecies(double* x, double* rsd, int* diag,
                           double rdt, size_t jmin, size_t jmax)
 {
-    StFlow::evalSpecies(x, rsd, diag, rdt, jmin, jmax);
+    Flow1D::evalSpecies(x, rsd, diag, rdt, jmin, jmax);
     if (m_stage != 2) {
         return;
     }
@@ -311,7 +311,7 @@ void IonFlow::setElectronTransport(vector<double>& tfix, vector<double>& diff_e,
 
 void IonFlow::_finalize(const double* x)
 {
-    StFlow::_finalize(x);
+    Flow1D::_finalize(x);
 
     bool p = m_do_electric_field[0];
     if (p) {
diff --git a/src/oneD/Sim1D.cpp b/src/oneD/Sim1D.cpp
index 555001dde59..3d03688afee 100644
--- a/src/oneD/Sim1D.cpp
+++ b/src/oneD/Sim1D.cpp
@@ -7,7 +7,7 @@
 
 #include "cantera/oneD/Sim1D.h"
 #include "cantera/oneD/MultiJac.h"
-#include "cantera/oneD/StFlow.h"
+#include "cantera/oneD/Flow1D.h"
 #include "cantera/oneD/MultiNewton.h"
 #include "cantera/oneD/refine.h"
 #include "cantera/numerics/funcs.h"
@@ -622,7 +622,7 @@ int Sim1D::setFixedTemperature(double t)
         size_t mfixed = npos;
 
         // loop over current grid to determine where new point is needed
-        StFlow* d_free = dynamic_cast<StFlow*>(&domain(n));
+        Flow1D* d_free = dynamic_cast<Flow1D*>(&domain(n));
         size_t npnow = d.nPoints();
         size_t nstart = znew.size();
         if (d_free && d_free->isFree()) {
@@ -702,7 +702,7 @@ double Sim1D::fixedTemperature()
 {
     double t_fixed = std::numeric_limits<double>::quiet_NaN();
     for (size_t n = 0; n < nDomains(); n++) {
-        StFlow* d = dynamic_cast<StFlow*>(&domain(n));
+        Flow1D* d = dynamic_cast<Flow1D*>(&domain(n));
         if (d && d->isFree() && d->m_tfixed > 0) {
             t_fixed = d->m_tfixed;
             break;
@@ -715,7 +715,7 @@ double Sim1D::fixedTemperatureLocation()
 {
     double z_fixed = std::numeric_limits<double>::quiet_NaN();
     for (size_t n = 0; n < nDomains(); n++) {
-        StFlow* d = dynamic_cast<StFlow*>(&domain(n));
+        Flow1D* d = dynamic_cast<Flow1D*>(&domain(n));
         if (d && d->isFree() && d->m_tfixed > 0) {
             z_fixed = d->m_zfixed;
             break;
@@ -735,7 +735,7 @@ void Sim1D::setLeftControlPoint(double temperature)
             continue;
         }
 
-        StFlow& d_axis = dynamic_cast<StFlow&>(domain(n));
+        Flow1D& d_axis = dynamic_cast<Flow1D&>(domain(n));
         size_t np = d_axis.nPoints();
 
         // Skip if two-point control is not enabled
@@ -786,7 +786,7 @@ void Sim1D::setRightControlPoint(double temperature)
             continue;
         }
 
-        StFlow& d_axis = dynamic_cast<StFlow&>(domain(n));
+        Flow1D& d_axis = dynamic_cast<Flow1D&>(domain(n));
         size_t np = d_axis.nPoints();
 
         // Skip if two-point control is not enabled
diff --git a/src/oneD/refine.cpp b/src/oneD/refine.cpp
index 719d26c47b2..5364e7e25e3 100644
--- a/src/oneD/refine.cpp
+++ b/src/oneD/refine.cpp
@@ -4,7 +4,7 @@
 // at https://cantera.org/license.txt for license and copyright information.
 
 #include "cantera/oneD/refine.h"
-#include "cantera/oneD/StFlow.h"
+#include "cantera/oneD/Flow1D.h"
 #include "cantera/base/global.h"
 
 using namespace std;
@@ -144,7 +144,7 @@ int Refiner::analyze(size_t n, const double* z, const double* x)
         }
     }
 
-    StFlow* fflame = dynamic_cast<StFlow*>(m_domain);
+    Flow1D* fflame = dynamic_cast<Flow1D*>(m_domain);
 
     // Refine based on properties of the grid itself
     for (size_t j = 1; j < n-1; j++) {
diff --git a/test/clib/test_ctonedim.cpp b/test/clib/test_ctonedim.cpp
index 5a7a178b7be..feff0abda88 100644
--- a/test/clib/test_ctonedim.cpp
+++ b/test/clib/test_ctonedim.cpp
@@ -26,7 +26,7 @@ TEST(ctonedim, freeflow)
     int flow = domain_new("free-flow", sol, "flow");
     ASSERT_GE(flow, 0);
     domain_setID(flow, "flow");
-    ASSERT_NEAR(stflow_pressure(flow), P, 1e-5);
+    ASSERT_NEAR(flow1D_pressure(flow), P, 1e-5);
 
     int buflen = domain_type3(flow, 0, 0);
     char* buf = new char[buflen];
@@ -56,10 +56,10 @@ TEST(ctonedim, freeflow_from_parts)
     thermo_setPressure(ph, P);
 
     int itype = 2; // free flow
-    int flow = stflow_new(ph, kin, tr, itype);
+    int flow = flow1D_new(ph, kin, tr, itype);
     ASSERT_GE(flow, 0);
     domain_setID(flow, "flow");
-    ASSERT_NEAR(stflow_pressure(flow), P, 1e-5);
+    ASSERT_NEAR(flow1D_pressure(flow), P, 1e-5);
 
     int buflen = domain_type3(flow, 0, 0);
     char* buf = new char[buflen];
@@ -124,7 +124,7 @@ TEST(ctonedim, catcomb_stack)
 
     // flow
     int itype = 1; // free flow
-    int flow = stflow_new(gas, gas_kin, trans, itype);
+    int flow = flow1D_new(gas, gas_kin, trans, itype);
     domain_setID(flow, "flow");
 
     // reacting surface
@@ -173,7 +173,7 @@ TEST(ctonedim, freeflame_from_parts)
 
     // flow
     int itype = 2; // free flow
-    int flow = stflow_new(ph, kin, tr, itype);
+    int flow = flow1D_new(ph, kin, tr, itype);
     domain_setID(flow, "flow");
 
     // grid
@@ -232,7 +232,7 @@ TEST(ctonedim, freeflame_from_parts)
     sim1D_setFixedTemperature(flame, 0.85 * T + .15 * Tad);
 
     // solve and save
-    stflow_solveEnergyEqn(flow, 1);
+    flow1D_solveEnergyEqn(flow, 1);
     bool refine_grid = false;
     int loglevel = 0;
     sim1D_solve(flame, loglevel, refine_grid);
diff --git a/test/oneD/test_oneD.cpp b/test/oneD/test_oneD.cpp
index 47de95d47b7..b765c4f2b6d 100644
--- a/test/oneD/test_oneD.cpp
+++ b/test/oneD/test_oneD.cpp
@@ -35,7 +35,7 @@ TEST(onedim, freeflame)
     double Tad = gas->temperature();
 
     // flow
-    auto flow = newDomain<StFlow>("free-flow", sol, "flow");
+    auto flow = newDomain<Flow1D>("free-flow", sol, "flow");
 
     // grid
     int nz = 21;
@@ -106,11 +106,11 @@ TEST(onedim, flame_types)
 {
     auto sol = newSolution("h2o2.yaml", "ohmech", "mixture-averaged");
 
-    auto free = newDomain<StFlow>("free-flow", sol, "flow");
+    auto free = newDomain<Flow1D>("free-flow", sol, "flow");
     ASSERT_EQ(free->type(), "free-flow");
-    auto symm = newDomain<StFlow>("axisymmetric-flow", sol, "flow");
+    auto symm = newDomain<Flow1D>("axisymmetric-flow", sol, "flow");
     ASSERT_EQ(symm->type(), "axisymmetric-flow");
-    auto burner = newDomain<StFlow>("unstrained-flow", sol, "flow");
+    auto burner = newDomain<Flow1D>("unstrained-flow", sol, "flow");
     ASSERT_EQ(burner->type(), "unstrained-flow");
 }
 

From 54c2e975e8f5d9857b268599d0d5fa92fe8afcfa Mon Sep 17 00:00:00 2001
From: Ingmar Schoegl <ischoegl@lsu.edu>
Date: Mon, 11 Sep 2023 19:01:38 -0500
Subject: [PATCH 2/8] [oneD] Remove Flow1D::wdot

---
 include/cantera/oneD/Flow1D.h | 4 ----
 src/oneD/Flow1D.cpp           | 4 ++--
 2 files changed, 2 insertions(+), 6 deletions(-)

diff --git a/include/cantera/oneD/Flow1D.h b/include/cantera/oneD/Flow1D.h
index a4bc7ccd526..27ccb8bb44d 100644
--- a/include/cantera/oneD/Flow1D.h
+++ b/include/cantera/oneD/Flow1D.h
@@ -396,10 +396,6 @@ class Flow1D : public Domain1D
     AnyMap getMeta() const override;
     void setMeta(const AnyMap& state) override;
 
-    double wdot(size_t k, size_t j) const {
-        return m_wdot(k,j);
-    }
-
     //! Update the properties (thermo, transport, and diffusion flux).
     //! This function is called in eval after the points which need
     //! to be updated are defined.
diff --git a/src/oneD/Flow1D.cpp b/src/oneD/Flow1D.cpp
index c28d5a41e5a..b12bfc92ccc 100644
--- a/src/oneD/Flow1D.cpp
+++ b/src/oneD/Flow1D.cpp
@@ -569,7 +569,7 @@ void Flow1D::evalEnergy(double* x, double* rsd, int* diag,
             double sum = 0.0;
             for (size_t k = 0; k < m_nsp; k++) {
                 double flxk = 0.5*(m_flux(k,j-1) + m_flux(k,j));
-                sum += wdot(k,j)*m_hk(k,j);
+                sum += m_wdot(k,j)*m_hk(k,j);
                 sum += flxk * m_dhk_dz(k,j) / m_wt[k];
             }
 
@@ -660,7 +660,7 @@ void Flow1D::evalSpecies(double* x, double* rsd, int* diag,
         for (size_t k = 0; k < m_nsp; k++) {
             double convec = rho_u(x,j)*dYdz(x,k,j);
             double diffus = 2.0*(m_flux(k,j) - m_flux(k,j-1)) / (z(j+1) - z(j-1));
-            rsd[index(c_offset_Y + k, j)] = (m_wt[k]*(wdot(k,j))
+            rsd[index(c_offset_Y + k, j)] = (m_wt[k]*m_wdot(k,j)
                                              - convec - diffus)/m_rho[j]
                                             - rdt*(Y(x,k,j) - Y_prev(k,j));
             diag[index(c_offset_Y + k, j)] = 1;

From 475cb81cedfd3823e336445aff25a931462cc923 Mon Sep 17 00:00:00 2001
From: Ingmar Schoegl <ischoegl@lsu.edu>
Date: Mon, 11 Sep 2023 19:48:56 -0500
Subject: [PATCH 3/8] [oneD] Reintroduce StFlow

---
 include/cantera/oneD/StFlow.h |  69 ++++++++
 src/oneD/DomainFactory.cpp    |   4 +
 src/oneD/StFlow.cpp           | 315 ++++++++++++++++++++++++++++++++++
 3 files changed, 388 insertions(+)
 create mode 100644 include/cantera/oneD/StFlow.h
 create mode 100644 src/oneD/StFlow.cpp

diff --git a/include/cantera/oneD/StFlow.h b/include/cantera/oneD/StFlow.h
new file mode 100644
index 00000000000..07a8df317f4
--- /dev/null
+++ b/include/cantera/oneD/StFlow.h
@@ -0,0 +1,69 @@
+//! @file StFlow.h
+
+// This file is part of Cantera. See License.txt in the top-level directory or
+// at https://cantera.org/license.txt for license and copyright information.
+
+#ifndef CT_STFLOW_H
+#define CT_STFLOW_H
+
+#include "Flow1D.h"
+
+namespace Cantera
+{
+
+/**
+ *  This class represents 1D flow domains that satisfy the one-dimensional
+ *  similarity solution for chemically-reacting, axisymmetric flows.
+ *
+ *  @deprecated To be removed after Cantera 3.1; replaced by Flow1D.
+ *  @ingroup flowGroup
+ */
+class StFlow : public Flow1D
+{
+public:
+    //! Create a new flow domain.
+    //! @param ph Object representing the gas phase. This object will be used
+    //!     to evaluate all thermodynamic, kinetic, and transport properties.
+    //! @param nsp Number of species.
+    //! @param points Initial number of grid points
+    StFlow(ThermoPhase* ph = 0, size_t nsp = 1, size_t points = 1);
+
+    //! Delegating constructor
+    StFlow(shared_ptr<ThermoPhase> th, size_t nsp = 1, size_t points = 1);
+
+    //! Create a new flow domain.
+    //! @param sol  Solution object used to evaluate all thermodynamic, kinetic, and
+    //!     transport properties
+    //! @param id  name of flow domain
+    //! @param points  initial number of grid points
+    StFlow(shared_ptr<Solution> sol, const string& id="", size_t points=1);
+
+    void eval(size_t j, double* x, double* r, integer* mask, double rdt) override;
+
+    //! Evaluate all residual components at the right boundary.
+    virtual void evalRightBoundary(double* x, double* res, int* diag, double rdt);
+
+    //! Evaluate the residual corresponding to the continuity equation at all
+    //! interior grid points.
+    virtual void evalContinuity(size_t j, double* x, double* r, int* diag, double rdt);
+
+protected:
+    double wdot(size_t k, size_t j) const {
+        return m_wdot(k,j);
+    }
+
+    //! Write the net production rates at point `j` into array `m_wdot`
+    void getWdot(double* x, size_t j) {
+        setGas(x,j);
+        m_kin->getNetProductionRates(&m_wdot(0,j));
+    }
+
+    //! Evaluate the residual function. This function is called in eval
+    //! after updateProperties is called.
+    virtual void evalResidual(double* x, double* rsd, int* diag,
+                              double rdt, size_t jmin, size_t jmax);
+};
+
+}
+
+#endif
diff --git a/src/oneD/DomainFactory.cpp b/src/oneD/DomainFactory.cpp
index 593c4ddc05d..de164e3a32b 100644
--- a/src/oneD/DomainFactory.cpp
+++ b/src/oneD/DomainFactory.cpp
@@ -7,6 +7,7 @@
 #include "cantera/oneD/Boundary1D.h"
 #include "cantera/oneD/Flow1D.h"
 #include "cantera/oneD/IonFlow.h"
+#include "cantera/oneD/StFlow.h"
 #include "cantera/transport/Transport.h"
 
 namespace Cantera
@@ -41,6 +42,9 @@ DomainFactory::DomainFactory()
     reg("gas-flow", [](shared_ptr<Solution> solution, const string& id) {
         return new Flow1D(solution, id);
     });
+    reg("legacy-flow", [](shared_ptr<Solution> solution, const string& id) {
+        return new StFlow(solution, id);
+    });
     reg("ion-flow", [](shared_ptr<Solution> solution, const string& id) {
         return new IonFlow(solution, id);
     });
diff --git a/src/oneD/StFlow.cpp b/src/oneD/StFlow.cpp
new file mode 100644
index 00000000000..f6631441aa5
--- /dev/null
+++ b/src/oneD/StFlow.cpp
@@ -0,0 +1,315 @@
+//! @file StFlow.cpp
+
+// This file is part of Cantera. See License.txt in the top-level directory or
+// at https://cantera.org/license.txt for license and copyright information.
+
+#include "cantera/oneD/StFlow.h"
+#include "cantera/base/global.h"
+
+using namespace std;
+
+namespace Cantera
+{
+
+StFlow::StFlow(ThermoPhase* ph, size_t nsp, size_t points) :
+    Flow1D(ph, nsp, points)
+{
+    warn_deprecated("StFlow::StFlow",
+        "To be removed after Cantera 3.1. Class replaced by Flow1D.");
+}
+
+StFlow::StFlow(shared_ptr<ThermoPhase> th, size_t nsp, size_t points)
+    : Flow1D(th, nsp, points)
+{
+    warn_deprecated("StFlow::StFlow",
+        "To be removed after Cantera 3.1. Class replaced by Flow1D.");
+}
+
+StFlow::StFlow(shared_ptr<Solution> sol, const string& id, size_t points)
+    : Flow1D(sol, id, points)
+{
+    warn_deprecated("StFlow::StFlow",
+        "To be removed after Cantera 3.1. Class replaced by Flow1D.");
+}
+
+void StFlow::eval(size_t jg, double* xg, double* rg, integer* diagg, double rdt)
+{
+    // if evaluating a Jacobian, and the global point is outside the domain of
+    // influence for this domain, then skip evaluating the residual
+    if (jg != npos && (jg + 1 < firstPoint() || jg > lastPoint() + 1)) {
+        return;
+    }
+
+    // start of local part of global arrays
+    double* x = xg + loc();
+    double* rsd = rg + loc();
+    integer* diag = diagg + loc();
+
+    size_t jmin, jmax;
+    if (jg == npos) { // evaluate all points
+        jmin = 0;
+        jmax = m_points - 1;
+    } else { // evaluate points for Jacobian
+        size_t jpt = (jg == 0) ? 0 : jg - firstPoint();
+        jmin = std::max<size_t>(jpt, 1) - 1;
+        jmax = std::min(jpt+1,m_points-1);
+    }
+
+    updateProperties(jg, x, jmin, jmax);
+    evalResidual(x, rsd, diag, rdt, jmin, jmax);
+    evalUo(x, rsd, diag, rdt, jmin, jmax);
+}
+
+void StFlow::evalResidual(double* x, double* rsd, int* diag,
+                          double rdt, size_t jmin, size_t jmax)
+{
+    //----------------------------------------------------
+    // evaluate the residual equations at all required
+    // grid points
+    //----------------------------------------------------
+
+    // calculation of qdotRadiation (see docstring of enableRadiation)
+    if (m_do_radiation) {
+        // variable definitions for the Planck absorption coefficient and the
+        // radiation calculation:
+        double k_P_ref = 1.0*OneAtm;
+
+        // polynomial coefficients:
+        const double c_H2O[6] = {-0.23093, -1.12390, 9.41530, -2.99880,
+                                     0.51382, -1.86840e-5};
+        const double c_CO2[6] = {18.741, -121.310, 273.500, -194.050,
+                                     56.310, -5.8169};
+
+        // calculation of the two boundary values
+        double boundary_Rad_left = m_epsilon_left * StefanBoltz * pow(T(x, 0), 4);
+        double boundary_Rad_right = m_epsilon_right * StefanBoltz * pow(T(x, m_points - 1), 4);
+
+        // loop over all grid points
+        for (size_t j = jmin; j < jmax; j++) {
+            // helping variable for the calculation
+            double radiative_heat_loss = 0;
+
+            // calculation of the mean Planck absorption coefficient
+            double k_P = 0;
+            // absorption coefficient for H2O
+            if (m_kRadiating[1] != npos) {
+                double k_P_H2O = 0;
+                for (size_t n = 0; n <= 5; n++) {
+                    k_P_H2O += c_H2O[n] * pow(1000 / T(x, j), (double) n);
+                }
+                k_P_H2O /= k_P_ref;
+                k_P += m_press * X(x, m_kRadiating[1], j) * k_P_H2O;
+            }
+            // absorption coefficient for CO2
+            if (m_kRadiating[0] != npos) {
+                double k_P_CO2 = 0;
+                for (size_t n = 0; n <= 5; n++) {
+                    k_P_CO2 += c_CO2[n] * pow(1000 / T(x, j), (double) n);
+                }
+                k_P_CO2 /= k_P_ref;
+                k_P += m_press * X(x, m_kRadiating[0], j) * k_P_CO2;
+            }
+
+            // calculation of the radiative heat loss term
+            radiative_heat_loss = 2 * k_P *(2 * StefanBoltz * pow(T(x, j), 4)
+            - boundary_Rad_left - boundary_Rad_right);
+
+            // set the radiative heat loss vector
+            m_qdotRadiation[j] = radiative_heat_loss;
+        }
+    }
+
+    for (size_t j = jmin; j <= jmax; j++) {
+        //----------------------------------------------
+        //         left boundary
+        //----------------------------------------------
+
+        if (j == 0) {
+            // these may be modified by a boundary object
+
+            // Continuity. This propagates information right-to-left, since
+            // rho_u at point 0 is dependent on rho_u at point 1, but not on
+            // mdot from the inlet.
+            rsd[index(c_offset_U,0)] =
+                -(rho_u(x,1) - rho_u(x,0))/m_dz[0]
+                -(density(1)*V(x,1) + density(0)*V(x,0));
+
+            // the inlet (or other) object connected to this one will modify
+            // these equations by subtracting its values for V, T, and mdot. As
+            // a result, these residual equations will force the solution
+            // variables to the values for the boundary object
+            rsd[index(c_offset_V,0)] = V(x,0);
+            rsd[index(c_offset_T,0)] = T(x,0);
+            if (m_usesLambda) {
+                rsd[index(c_offset_L, 0)] = -rho_u(x, 0);
+            } else {
+                rsd[index(c_offset_L, 0)] = lambda(x, 0);
+                diag[index(c_offset_L, 0)] = 0;
+            }
+
+            // The default boundary condition for species is zero flux. However,
+            // the boundary object may modify this.
+            double sum = 0.0;
+            for (size_t k = 0; k < m_nsp; k++) {
+                sum += Y(x,k,0);
+                rsd[index(c_offset_Y + k, 0)] =
+                    -(m_flux(k,0) + rho_u(x,0)* Y(x,k,0));
+            }
+            rsd[index(c_offset_Y + leftExcessSpecies(), 0)] = 1.0 - sum;
+
+            // set residual of poisson's equ to zero
+            rsd[index(c_offset_E, 0)] = x[index(c_offset_E, j)];
+        } else if (j == m_points - 1) {
+            evalRightBoundary(x, rsd, diag, rdt);
+        } else { // interior points
+            evalContinuity(j, x, rsd, diag, rdt);
+            // set residual of poisson's equ to zero
+            rsd[index(c_offset_E, j)] = x[index(c_offset_E, j)];
+
+            //------------------------------------------------
+            //    Radial momentum equation
+            //
+            //    \rho dV/dt + \rho u dV/dz + \rho V^2
+            //       = d(\mu dV/dz)/dz - lambda
+            //-------------------------------------------------
+            if (m_usesLambda) {
+                rsd[index(c_offset_V,j)] =
+                    (shear(x, j) - lambda(x, j) - rho_u(x, j) * dVdz(x, j)
+                    - m_rho[j] * V(x, j) * V(x, j)) / m_rho[j]
+                    - rdt * (V(x, j) - V_prev(j));
+                diag[index(c_offset_V, j)] = 1;
+            } else {
+                rsd[index(c_offset_V, j)] = V(x, j);
+                diag[index(c_offset_V, j)] = 0;
+            }
+
+            //-------------------------------------------------
+            //    Species equations
+            //
+            //   \rho dY_k/dt + \rho u dY_k/dz + dJ_k/dz
+            //   = M_k\omega_k
+            //-------------------------------------------------
+            getWdot(x,j);
+            for (size_t k = 0; k < m_nsp; k++) {
+                double convec = rho_u(x,j)*dYdz(x,k,j);
+                double diffus = 2.0*(m_flux(k,j) - m_flux(k,j-1))
+                                / (z(j+1) - z(j-1));
+                rsd[index(c_offset_Y + k, j)]
+                = (m_wt[k]*(wdot(k,j))
+                   - convec - diffus)/m_rho[j]
+                  - rdt*(Y(x,k,j) - Y_prev(k,j));
+                diag[index(c_offset_Y + k, j)] = 1;
+            }
+
+            //-----------------------------------------------
+            //    energy equation
+            //
+            //    \rho c_p dT/dt + \rho c_p u dT/dz
+            //    = d(k dT/dz)/dz
+            //      - sum_k(\omega_k h_k_ref)
+            //      - sum_k(J_k c_p_k / M_k) dT/dz
+            //-----------------------------------------------
+            if (m_do_energy[j]) {
+
+                setGas(x,j);
+                double dtdzj = dTdz(x,j);
+                double sum = 0.0;
+
+                grad_hk(x, j);
+                for (size_t k = 0; k < m_nsp; k++) {
+                    double flxk = 0.5*(m_flux(k,j-1) + m_flux(k,j));
+                    sum += wdot(k,j)*m_hk(k,j);
+                    sum += flxk * m_dhk_dz(k,j) / m_wt[k];
+                }
+
+                rsd[index(c_offset_T, j)] = - m_cp[j]*rho_u(x,j)*dtdzj
+                                            - divHeatFlux(x,j) - sum;
+                rsd[index(c_offset_T, j)] /= (m_rho[j]*m_cp[j]);
+                rsd[index(c_offset_T, j)] -= rdt*(T(x,j) - T_prev(j));
+                rsd[index(c_offset_T, j)] -= (m_qdotRadiation[j] / (m_rho[j] * m_cp[j]));
+                diag[index(c_offset_T, j)] = 1;
+            } else {
+                // residual equations if the energy equation is disabled
+                rsd[index(c_offset_T, j)] = T(x,j) - T_fixed(j);
+                diag[index(c_offset_T, j)] = 0;
+            }
+
+            if (m_usesLambda) {
+                rsd[index(c_offset_L, j)] = lambda(x, j) - lambda(x, j - 1);
+            } else {
+                rsd[index(c_offset_L, j)] = lambda(x, j);
+            }
+            diag[index(c_offset_L, j)] = 0;
+        }
+    }
+}
+
+void StFlow::evalRightBoundary(double* x, double* rsd, int* diag, double rdt)
+{
+    size_t j = m_points - 1;
+
+    // the boundary object connected to the right of this one may modify or
+    // replace these equations. The default boundary conditions are zero u, V,
+    // and T, and zero diffusive flux for all species.
+
+    rsd[index(c_offset_V,j)] = V(x,j);
+    diag[index(c_offset_V,j)] = 0;
+    double sum = 0.0;
+    // set residual of poisson's equ to zero
+    rsd[index(c_offset_E, j)] = x[index(c_offset_E, j)];
+    for (size_t k = 0; k < m_nsp; k++) {
+        sum += Y(x,k,j);
+        rsd[index(k+c_offset_Y,j)] = m_flux(k,j-1) + rho_u(x,j)*Y(x,k,j);
+    }
+    rsd[index(c_offset_Y + rightExcessSpecies(), j)] = 1.0 - sum;
+    diag[index(c_offset_Y + rightExcessSpecies(), j)] = 0;
+    if (m_usesLambda) {
+        rsd[index(c_offset_U, j)] = rho_u(x, j);
+    } else {
+        rsd[index(c_offset_U, j)] = rho_u(x, j) - rho_u(x, j-1);
+    }
+
+    rsd[index(c_offset_L, j)] = lambda(x, j) - lambda(x, j-1);
+    diag[index(c_offset_L, j)] = 0;
+    rsd[index(c_offset_T, j)] = T(x, j);
+}
+
+void StFlow::evalContinuity(size_t j, double* x, double* rsd, int* diag, double rdt)
+{
+    //algebraic constraint
+    diag[index(c_offset_U, j)] = 0;
+    //----------------------------------------------
+    //    Continuity equation
+    //
+    //    d(\rho u)/dz + 2\rho V = 0
+    //----------------------------------------------
+    if (m_usesLambda) {
+        // Note that this propagates the mass flow rate information to the left
+        // (j+1 -> j) from the value specified at the right boundary. The
+        // lambda information propagates in the opposite direction.
+        rsd[index(c_offset_U,j)] =
+            -(rho_u(x,j+1) - rho_u(x,j))/m_dz[j]
+            -(density(j+1)*V(x,j+1) + density(j)*V(x,j));
+    } else if (m_isFree) {
+        // terms involving V are zero as V=0 by definition
+        if (grid(j) > m_zfixed) {
+            rsd[index(c_offset_U,j)] =
+                - (rho_u(x,j) - rho_u(x,j-1))/m_dz[j-1];
+        } else if (grid(j) == m_zfixed) {
+            if (m_do_energy[j]) {
+                rsd[index(c_offset_U,j)] = (T(x,j) - m_tfixed);
+            } else {
+                rsd[index(c_offset_U,j)] = (rho_u(x,j)
+                                            - m_rho[0]*0.3); // why 0.3?
+            }
+        } else if (grid(j) < m_zfixed) {
+            rsd[index(c_offset_U,j)] =
+                - (rho_u(x,j+1) - rho_u(x,j))/m_dz[j];
+        }
+    } else {
+        // unstrained with fixed mass flow rate
+        rsd[index(c_offset_U, j)] = rho_u(x, j) - rho_u(x, j - 1);
+    }
+}
+
+} // namespace

From 51265a44b159c380a85f3633f8cd0c75c9277b61 Mon Sep 17 00:00:00 2001
From: Ingmar Schoegl <ischoegl@lsu.edu>
Date: Mon, 11 Sep 2023 19:49:33 -0500
Subject: [PATCH 4/8] [unittests] Cover legacy StFlow class

---
 test/oneD/test_oneD.cpp | 93 +++++++++++++++++++++++++++++++++++++++++
 1 file changed, 93 insertions(+)

diff --git a/test/oneD/test_oneD.cpp b/test/oneD/test_oneD.cpp
index b765c4f2b6d..0e5e1d72ebd 100644
--- a/test/oneD/test_oneD.cpp
+++ b/test/oneD/test_oneD.cpp
@@ -5,6 +5,7 @@
 #include "cantera/core.h"
 #include "cantera/onedim.h"
 #include "cantera/oneD/DomainFactory.h"
+#include "cantera/oneD/StFlow.h"
 #include "cantera/oneD/IonFlow.h"
 
 using namespace Cantera;
@@ -102,6 +103,98 @@ TEST(onedim, freeflame)
     }
 }
 
+TEST(onedim, legacy)
+{
+    //! @todo: Remove after Cantera 3.1
+    auto sol = newSolution("h2o2.yaml", "ohmech", "mixture-averaged");
+    auto gas = sol->thermo();
+    size_t nsp = gas->nSpecies();
+
+    // reactants
+    double uin = .3;
+    double T = 300;
+    double P = 101325;
+    string X = "H2:0.65, O2:0.5, AR:2";
+    gas->setState_TPX(T, P, X);
+    double rho_in = gas->density();
+    vector<double> yin(nsp);
+    gas->getMassFractions(&yin[0]);
+
+    // product estimate
+    gas->equilibrate("HP");
+    vector<double> yout(nsp);
+    gas->getMassFractions(&yout[0]);
+    double rho_out = gas->density();
+    double Tad = gas->temperature();
+
+    // flow
+    suppress_deprecation_warnings();
+    auto flow = newDomain<StFlow>("legacy-flow", sol);
+    make_deprecation_warnings_fatal();
+    flow->setID("flow");
+    flow->setFreeFlow();
+
+    // grid
+    int nz = 21;
+    double lz = 0.02;
+    vector<double> z(nz);
+    double dz = lz;
+    dz /= (double)(nz - 1);
+    for (int iz = 0; iz < nz; iz++) {
+        z[iz] = iz * dz;
+    }
+    flow->setupGrid(nz, &z[0]);
+
+    // inlet
+    auto inlet = newDomain<Inlet1D>("inlet", sol);
+    inlet->setMoleFractions(X);
+    inlet->setMdot(uin * rho_in);
+    inlet->setTemperature(T);
+
+    // outlet
+    auto outlet = newDomain<Outlet1D>("outlet", sol);
+    double uout = inlet->mdot() / rho_out;
+
+    // set up simulation
+    vector<shared_ptr<Domain1D>> domains { inlet, flow, outlet };
+    Sim1D flame(domains);
+    int dom = static_cast<int>(flame.domainIndex("flow"));
+    ASSERT_EQ(dom, 1);
+
+    // set up initial guess
+    vector<double> locs{0.0, 0.3, 0.7, 1.0};
+    vector<double> value{uin, uin, uout, uout};
+    flame.setInitialGuess("velocity", locs, value);
+    value = {T, T, Tad, Tad};
+    flame.setInitialGuess("T", locs, value);
+    for (size_t i = 0; i < nsp; i++) {
+        value = {yin[i], yin[i], yout[i], yout[i]};
+        flame.setInitialGuess(gas->speciesName(i), locs, value);
+    }
+
+    // simulation settings
+    double ratio = 15.0;
+    double slope = 0.3;
+    double curve = 0.5;
+    flame.setRefineCriteria(dom, ratio, slope, curve);
+    flame.setFixedTemperature(0.85 * T + .15 * Tad);
+
+    // solve
+    flow->solveEnergyEqn();
+    bool refine_grid = false;
+    int loglevel = 0;
+    flame.solve(loglevel, refine_grid);
+
+    ASSERT_EQ(flow->nPoints(), static_cast<size_t>(nz + 1));
+    size_t comp = flow->componentIndex("T");
+    double Tprev = flame.value(dom, comp, 0);
+    for (size_t n = 0; n < flow->nPoints(); n++) {
+        T = flame.value(dom, comp, n);
+        ASSERT_GE(T, Tprev);
+        Tprev = T;
+    }
+}
+
 TEST(onedim, flame_types)
 {
     auto sol = newSolution("h2o2.yaml", "ohmech", "mixture-averaged");

From 6f4066aca0da469336d22b165b6939076665ab78 Mon Sep 17 00:00:00 2001
From: Ingmar Schoegl <ischoegl@lsu.edu>
Date: Mon, 11 Sep 2023 20:15:30 -0500
Subject: [PATCH 5/8] [oneD] Prevent compiler warning

---
 include/cantera/oneD/Flow1D.h | 6 ++++++
 include/cantera/oneD/StFlow.h | 4 +---
 src/oneD/Flow1D.cpp           | 6 ++++++
 3 files changed, 13 insertions(+), 3 deletions(-)

diff --git a/include/cantera/oneD/Flow1D.h b/include/cantera/oneD/Flow1D.h
index 27ccb8bb44d..840bcf2d1d6 100644
--- a/include/cantera/oneD/Flow1D.h
+++ b/include/cantera/oneD/Flow1D.h
@@ -542,6 +542,12 @@ class Flow1D : public Domain1D
     virtual void evalElectricField(double* x, double* rsd, int* diag,
                                    double rdt, size_t jmin, size_t jmax);
 
+    //! Alternate version of evalContinuity with legacy signature.
+    //! Implemented by StFlow; included here to prevent compiler warnings about shadowed
+    //! virtual functions.
+    //! @deprecated To be removed after Cantera 3.1.
+    virtual void evalContinuity(size_t j, double* x, double* r, int* diag, double rdt);
+
     /**
      * Evaluate the oxidizer axial velocity equation residual.
      *
diff --git a/include/cantera/oneD/StFlow.h b/include/cantera/oneD/StFlow.h
index 07a8df317f4..40e826a7fcb 100644
--- a/include/cantera/oneD/StFlow.h
+++ b/include/cantera/oneD/StFlow.h
@@ -43,9 +43,7 @@ class StFlow : public Flow1D
     //! Evaluate all residual components at the right boundary.
     virtual void evalRightBoundary(double* x, double* res, int* diag, double rdt);
 
-    //! Evaluate the residual corresponding to the continuity equation at all
-    //! interior grid points.
-    virtual void evalContinuity(size_t j, double* x, double* r, int* diag, double rdt);
+    void evalContinuity(size_t j, double* x, double* r, int* diag, double rdt) override;
 
 protected:
     double wdot(size_t k, size_t j) const {
diff --git a/src/oneD/Flow1D.cpp b/src/oneD/Flow1D.cpp
index b12bfc92ccc..f1e0808e2b5 100644
--- a/src/oneD/Flow1D.cpp
+++ b/src/oneD/Flow1D.cpp
@@ -677,6 +677,12 @@ void Flow1D::evalElectricField(double* x, double* rsd, int* diag,
     }
 }
 
+void Flow1D::evalContinuity(size_t j, double* x, double* rsd, int* diag, double rdt)
+{
+    throw CanteraError("Flow1D::evalContinuity",
+        "Overloaded by StFlow; to be removed after Cantera 3.1");
+}
+
 void Flow1D::updateTransport(double* x, size_t j0, size_t j1)
 {
      if (m_do_multicomponent) {

From 0fe746e965bcde659264221cb86e716c2cd87d82 Mon Sep 17 00:00:00 2001
From: Ingmar Schoegl <ischoegl@lsu.edu>
Date: Mon, 11 Sep 2023 20:23:12 -0500
Subject: [PATCH 6/8] [oneD] Add named doxygen section for governing equations

Also, add missing member variable documentation.
---
 include/cantera/oneD/Flow1D.h | 16 +++++++++++-----
 1 file changed, 11 insertions(+), 5 deletions(-)

diff --git a/include/cantera/oneD/Flow1D.h b/include/cantera/oneD/Flow1D.h
index 840bcf2d1d6..fa6f6046c85 100644
--- a/include/cantera/oneD/Flow1D.h
+++ b/include/cantera/oneD/Flow1D.h
@@ -418,6 +418,10 @@ class Flow1D : public Domain1D
      */
     void computeRadiation(double* x, size_t jmin, size_t jmax);
 
+    //! @name Governing Equations
+    //! Methods are used to evaluate residuals for each of the governing equations.
+    //! @{
+
     /**
      * Evaluate the continuity equation residual.
      *
@@ -542,6 +546,8 @@ class Flow1D : public Domain1D
     virtual void evalElectricField(double* x, double* rsd, int* diag,
                                    double rdt, size_t jmin, size_t jmax);
 
+    //! @} End of Governing Equations
+
     //! Alternate version of evalContinuity with legacy signature.
     //! Implemented by StFlow; included here to prevent compiler warnings about shadowed
     //! virtual functions.
@@ -649,7 +655,7 @@ class Flow1D : public Domain1D
     }
     //! @}
 
-    //! @name convective spatial derivatives.
+    //! @name Convective spatial derivatives
     //!
     //! These use upwind differencing, assuming u(z) is negative
     //! @{
@@ -701,12 +707,12 @@ class Flow1D : public Domain1D
     vector<double> m_dz;
 
     // mixture thermo properties
-    vector<double> m_rho;
-    vector<double> m_wtm;
+    vector<double> m_rho;  //!< Vector of size #m_nsp to cache densities
+    vector<double> m_wtm;  //!< Vector of size #m_nsp to cache mean molecular weights
 
     // species thermo properties
     vector<double> m_wt;
-    vector<double> m_cp;
+    vector<double> m_cp;  //!< Vector of size #m_nsp to cache specific heat capacities
 
     // transport properties
     vector<double> m_visc;
@@ -724,7 +730,7 @@ class Flow1D : public Domain1D
     //! Array of size #m_nsp by #m_points-1 for saving enthalpy fluxes
     Array2D m_dhk_dz;
 
-    // production rates
+    //! Array of size #m_nsp by #m_points for saving species production rates
     Array2D m_wdot;
 
     size_t m_nsp; //!< Number of species in the mechanism

From 1d8f1c682e8cbec96cb6953818a72df60641d8d1 Mon Sep 17 00:00:00 2001
From: Ingmar Schoegl <ischoegl@lsu.edu>
Date: Sat, 22 Jun 2024 11:31:40 +0200
Subject: [PATCH 7/8] [oneD] Reorder content for named doxygen section

---
 include/cantera/oneD/Flow1D.h |  76 +++++++-------
 src/oneD/Flow1D.cpp           | 188 +++++++++++++++++-----------------
 2 files changed, 135 insertions(+), 129 deletions(-)

diff --git a/include/cantera/oneD/Flow1D.h b/include/cantera/oneD/Flow1D.h
index fa6f6046c85..aa48218b8e8 100644
--- a/include/cantera/oneD/Flow1D.h
+++ b/include/cantera/oneD/Flow1D.h
@@ -172,7 +172,6 @@ class Flow1D : public Domain1D
     //! Returns true if the specified component is an active part of the solver state
     virtual bool componentActive(size_t n) const;
 
-    //! Print the solution.
     void show(const double* x) override;
 
     shared_ptr<SolutionArray> asArray(const double* soln) const override;
@@ -396,6 +395,43 @@ class Flow1D : public Domain1D
     AnyMap getMeta() const override;
     void setMeta(const AnyMap& state) override;
 
+    //! @name Updates of cached properties
+    //! These methods are called by eval() to update cached properties and data that are
+    //! used for the evaluation of the governing equations.
+    //! @{
+
+    /**
+     * Update the thermodynamic properties from point j0 to point j1
+     * (inclusive), based on solution x.
+     *
+     * The gas state is set to be consistent with the solution at the
+     * points from j0 to j1.
+     *
+     * Properties that are computed and cached are:
+     * * #m_rho (density)
+     * * #m_wtm (mean molecular weight)
+     * * #m_cp (specific heat capacity)
+     * * #m_hk (species specific enthalpies)
+     * * #m_wdot (species production rates)
+     */
+    void updateThermo(const double* x, size_t j0, size_t j1) {
+        for (size_t j = j0; j <= j1; j++) {
+            setGas(x,j);
+            m_rho[j] = m_thermo->density();
+            m_wtm[j] = m_thermo->meanMolecularWeight();
+            m_cp[j] = m_thermo->cp_mass();
+            m_thermo->getPartialMolarEnthalpies(&m_hk(0, j));
+            m_kin->getNetProductionRates(&m_wdot(0, j));
+        }
+    }
+
+    //! Update the transport properties at grid points in the range from `j0`
+    //! to `j1`, based on solution `x`.
+    virtual void updateTransport(double* x, size_t j0, size_t j1);
+
+    //! Update the diffusive mass fluxes.
+    virtual void updateDiffFluxes(const double* x, size_t j0, size_t j1);
+
     //! Update the properties (thermo, transport, and diffusion flux).
     //! This function is called in eval after the points which need
     //! to be updated are defined.
@@ -418,9 +454,11 @@ class Flow1D : public Domain1D
      */
     void computeRadiation(double* x, size_t jmin, size_t jmax);
 
+    //! @}
+
     //! @name Governing Equations
-    //! Methods are used to evaluate residuals for each of the governing equations.
-    //! @{
+    //! Methods called by eval() to calculate residuals for individual governing
+    //! equations.
 
     /**
      * Evaluate the continuity equation residual.
@@ -572,31 +610,6 @@ class Flow1D : public Domain1D
     virtual void evalUo(double* x, double* rsd, int* diag,
                         double rdt, size_t jmin, size_t jmax);
 
-    /**
-     * Update the thermodynamic properties from point j0 to point j1
-     * (inclusive), based on solution x.
-     *
-     * The gas state is set to be consistent with the solution at the
-     * points from j0 to j1.
-     *
-     * Properties that are computed and cached are:
-     * * #m_rho (density)
-     * * #m_wtm (mean molecular weight)
-     * * #m_cp (specific heat capacity)
-     * * #m_hk (species specific enthalpies)
-     * * #m_wdot (species production rates)
-     */
-    void updateThermo(const double* x, size_t j0, size_t j1) {
-        for (size_t j = j0; j <= j1; j++) {
-            setGas(x,j);
-            m_rho[j] = m_thermo->density();
-            m_wtm[j] = m_thermo->meanMolecularWeight();
-            m_cp[j] = m_thermo->cp_mass();
-            m_thermo->getPartialMolarEnthalpies(&m_hk(0, j));
-            m_kin->getNetProductionRates(&m_wdot(0, j));
-        }
-    }
-
     //! @name Solution components
     //! @{
 
@@ -691,9 +704,6 @@ class Flow1D : public Domain1D
         return m*m_nsp*m_nsp + m_nsp*j + k;
     }
 
-    //! Update the diffusive mass fluxes.
-    virtual void updateDiffFluxes(const double* x, size_t j0, size_t j1);
-
     //! Get the gradient of species specific molar enthalpies
     virtual void grad_hk(const double* x, size_t j);
 
@@ -775,10 +785,6 @@ class Flow1D : public Domain1D
     bool m_isFree;
     bool m_usesLambda;
 
-    //! Update the transport properties at grid points in the range from `j0`
-    //! to `j1`, based on solution `x`.
-    virtual void updateTransport(double* x, size_t j0, size_t j1);
-
     //! Flag for activating two-point flame control
     bool m_twoPointControl = false;
 
diff --git a/src/oneD/Flow1D.cpp b/src/oneD/Flow1D.cpp
index f1e0808e2b5..cb82b164d13 100644
--- a/src/oneD/Flow1D.cpp
+++ b/src/oneD/Flow1D.cpp
@@ -368,6 +368,100 @@ void Flow1D::updateProperties(size_t jg, double* x, size_t jmin, size_t jmax)
     updateDiffFluxes(x, j0, j1);
 }
 
+void Flow1D::updateTransport(double* x, size_t j0, size_t j1)
+{
+     if (m_do_multicomponent) {
+        for (size_t j = j0; j < j1; j++) {
+            setGasAtMidpoint(x,j);
+            double wtm = m_thermo->meanMolecularWeight();
+            double rho = m_thermo->density();
+            m_visc[j] = (m_dovisc ? m_trans->viscosity() : 0.0);
+            m_trans->getMultiDiffCoeffs(m_nsp, &m_multidiff[mindex(0,0,j)]);
+
+            // Use m_diff as storage for the factor outside the summation
+            for (size_t k = 0; k < m_nsp; k++) {
+                m_diff[k+j*m_nsp] = m_wt[k] * rho / (wtm*wtm);
+            }
+
+            m_tcon[j] = m_trans->thermalConductivity();
+            if (m_do_soret) {
+                m_trans->getThermalDiffCoeffs(m_dthermal.ptrColumn(0) + j*m_nsp);
+            }
+        }
+    } else { // mixture averaged transport
+        for (size_t j = j0; j < j1; j++) {
+            setGasAtMidpoint(x,j);
+            m_visc[j] = (m_dovisc ? m_trans->viscosity() : 0.0);
+
+            if (m_fluxGradientBasis == ThermoBasis::molar) {
+                m_trans->getMixDiffCoeffs(&m_diff[j*m_nsp]);
+            } else {
+                m_trans->getMixDiffCoeffsMass(&m_diff[j*m_nsp]);
+            }
+
+            double rho = m_thermo->density();
+
+            if (m_fluxGradientBasis == ThermoBasis::molar) {
+                double wtm = m_thermo->meanMolecularWeight();
+                for (size_t k=0; k < m_nsp; k++) {
+                    m_diff[k+j*m_nsp] *= m_wt[k] * rho / wtm;
+                }
+            } else {
+                for (size_t k=0; k < m_nsp; k++) {
+                    m_diff[k+j*m_nsp] *= rho;
+                }
+            }
+            m_tcon[j] = m_trans->thermalConductivity();
+        }
+    }
+}
+
+void Flow1D::updateDiffFluxes(const double* x, size_t j0, size_t j1)
+{
+    if (m_do_multicomponent) {
+        for (size_t j = j0; j < j1; j++) {
+            double dz = z(j+1) - z(j);
+            for (size_t k = 0; k < m_nsp; k++) {
+                double sum = 0.0;
+                for (size_t m = 0; m < m_nsp; m++) {
+                    sum += m_wt[m] * m_multidiff[mindex(k,m,j)] * (X(x,m,j+1)-X(x,m,j));
+                }
+                m_flux(k,j) = sum * m_diff[k+j*m_nsp] / dz;
+            }
+        }
+    } else {
+        for (size_t j = j0; j < j1; j++) {
+            double sum = 0.0;
+            double dz = z(j+1) - z(j);
+            if (m_fluxGradientBasis == ThermoBasis::molar) {
+                for (size_t k = 0; k < m_nsp; k++) {
+                    m_flux(k,j) = m_diff[k+m_nsp*j] * (X(x,k,j) - X(x,k,j+1))/dz;
+                    sum -= m_flux(k,j);
+                }
+            } else {
+                for (size_t k = 0; k < m_nsp; k++) {
+                    m_flux(k,j) = m_diff[k+m_nsp*j] * (Y(x,k,j) - Y(x,k,j+1))/dz;
+                    sum -= m_flux(k,j);
+                }
+            }
+            // correction flux to insure that \sum_k Y_k V_k = 0.
+            for (size_t k = 0; k < m_nsp; k++) {
+                m_flux(k,j) += sum*Y(x,k,j);
+            }
+        }
+    }
+
+    if (m_do_soret) {
+        for (size_t m = j0; m < j1; m++) {
+            double gradlogT = 2.0 * (T(x,m+1) - T(x,m)) /
+                              ((T(x,m+1) + T(x,m)) * (z(m+1) - z(m)));
+            for (size_t k = 0; k < m_nsp; k++) {
+                m_flux(k,m) -= m_dthermal(k,m)*gradlogT;
+            }
+        }
+    }
+}
+
 void Flow1D::computeRadiation(double* x, size_t jmin, size_t jmax)
 {
     // Variable definitions for the Planck absorption coefficient and the
@@ -683,54 +777,6 @@ void Flow1D::evalContinuity(size_t j, double* x, double* rsd, int* diag, double
         "Overloaded by StFlow; to be removed after Cantera 3.1");
 }
 
-void Flow1D::updateTransport(double* x, size_t j0, size_t j1)
-{
-     if (m_do_multicomponent) {
-        for (size_t j = j0; j < j1; j++) {
-            setGasAtMidpoint(x,j);
-            double wtm = m_thermo->meanMolecularWeight();
-            double rho = m_thermo->density();
-            m_visc[j] = (m_dovisc ? m_trans->viscosity() : 0.0);
-            m_trans->getMultiDiffCoeffs(m_nsp, &m_multidiff[mindex(0,0,j)]);
-
-            // Use m_diff as storage for the factor outside the summation
-            for (size_t k = 0; k < m_nsp; k++) {
-                m_diff[k+j*m_nsp] = m_wt[k] * rho / (wtm*wtm);
-            }
-
-            m_tcon[j] = m_trans->thermalConductivity();
-            if (m_do_soret) {
-                m_trans->getThermalDiffCoeffs(m_dthermal.ptrColumn(0) + j*m_nsp);
-            }
-        }
-    } else { // mixture averaged transport
-        for (size_t j = j0; j < j1; j++) {
-            setGasAtMidpoint(x,j);
-            m_visc[j] = (m_dovisc ? m_trans->viscosity() : 0.0);
-
-            if (m_fluxGradientBasis == ThermoBasis::molar) {
-                m_trans->getMixDiffCoeffs(&m_diff[j*m_nsp]);
-            } else {
-                m_trans->getMixDiffCoeffsMass(&m_diff[j*m_nsp]);
-            }
-
-            double rho = m_thermo->density();
-
-            if (m_fluxGradientBasis == ThermoBasis::molar) {
-                double wtm = m_thermo->meanMolecularWeight();
-                for (size_t k=0; k < m_nsp; k++) {
-                    m_diff[k+j*m_nsp] *= m_wt[k] * rho / wtm;
-                }
-            } else {
-                for (size_t k=0; k < m_nsp; k++) {
-                    m_diff[k+j*m_nsp] *= rho;
-                }
-            }
-            m_tcon[j] = m_trans->thermalConductivity();
-        }
-    }
-}
-
 void Flow1D::show(const double* x)
 {
     writelog("    Pressure:  {:10.4g} Pa\n", m_press);
@@ -748,52 +794,6 @@ void Flow1D::show(const double* x)
     }
 }
 
-void Flow1D::updateDiffFluxes(const double* x, size_t j0, size_t j1)
-{
-    if (m_do_multicomponent) {
-        for (size_t j = j0; j < j1; j++) {
-            double dz = z(j+1) - z(j);
-            for (size_t k = 0; k < m_nsp; k++) {
-                double sum = 0.0;
-                for (size_t m = 0; m < m_nsp; m++) {
-                    sum += m_wt[m] * m_multidiff[mindex(k,m,j)] * (X(x,m,j+1)-X(x,m,j));
-                }
-                m_flux(k,j) = sum * m_diff[k+j*m_nsp] / dz;
-            }
-        }
-    } else {
-        for (size_t j = j0; j < j1; j++) {
-            double sum = 0.0;
-            double dz = z(j+1) - z(j);
-            if (m_fluxGradientBasis == ThermoBasis::molar) {
-                for (size_t k = 0; k < m_nsp; k++) {
-                    m_flux(k,j) = m_diff[k+m_nsp*j] * (X(x,k,j) - X(x,k,j+1))/dz;
-                    sum -= m_flux(k,j);
-                }
-            } else {
-                for (size_t k = 0; k < m_nsp; k++) {
-                    m_flux(k,j) = m_diff[k+m_nsp*j] * (Y(x,k,j) - Y(x,k,j+1))/dz;
-                    sum -= m_flux(k,j);
-                }
-            }
-            // correction flux to insure that \sum_k Y_k V_k = 0.
-            for (size_t k = 0; k < m_nsp; k++) {
-                m_flux(k,j) += sum*Y(x,k,j);
-            }
-        }
-    }
-
-    if (m_do_soret) {
-        for (size_t m = j0; m < j1; m++) {
-            double gradlogT = 2.0 * (T(x,m+1) - T(x,m)) /
-                              ((T(x,m+1) + T(x,m)) * (z(m+1) - z(m)));
-            for (size_t k = 0; k < m_nsp; k++) {
-                m_flux(k,m) -= m_dthermal(k,m)*gradlogT;
-            }
-        }
-    }
-}
-
 string Flow1D::componentName(size_t n) const
 {
     switch (n) {

From 0e7e224273d2074784e67d669d668085cc55cc11 Mon Sep 17 00:00:00 2001
From: Ingmar Schoegl <ischoegl@lsu.edu>
Date: Sun, 23 Jun 2024 08:51:22 +0200
Subject: [PATCH 8/8] [clib] Retain stflow_ routines but throw error

---
 include/cantera/clib/ctonedim.h | 10 +++++
 src/clib/ctonedim.cpp           | 73 ++++++++++++++++++++++++++++++++-
 test/clib/test_ctonedim.cpp     | 22 ++++++++++
 test/clib/utils.h               |  0
 4 files changed, 104 insertions(+), 1 deletion(-)
 create mode 100644 test/clib/utils.h

diff --git a/include/cantera/clib/ctonedim.h b/include/cantera/clib/ctonedim.h
index e1267daf760..b5f3023a5e1 100644
--- a/include/cantera/clib/ctonedim.h
+++ b/include/cantera/clib/ctonedim.h
@@ -68,6 +68,16 @@ extern "C" {
             size_t m, const double* temp);
     CANTERA_CAPI int flow1D_solveEnergyEqn(int i, int flag);
 
+    //! @todo: Remove all functions with `stflow` prefix after Cantera 3.1
+    CANTERA_CAPI int stflow_new(int iph, int ikin, int itr, int itype);
+    CANTERA_CAPI int stflow_setTransport(int i, int itr);
+    CANTERA_CAPI int stflow_enableSoret(int i, int iSoret);
+    CANTERA_CAPI int stflow_setPressure(int i, double p);
+    CANTERA_CAPI double stflow_pressure(int i);
+    CANTERA_CAPI int stflow_setFixedTempProfile(int i, size_t n, const double* pos,
+            size_t m, const double* temp);
+    CANTERA_CAPI int stflow_solveEnergyEqn(int i, int flag);
+
     CANTERA_CAPI int sim1D_new(size_t nd, const int* domains);
     CANTERA_CAPI int sim1D_del(int i);
     CANTERA_CAPI int sim1D_setValue(int i, int dom, int comp, int localPoint, double value);
diff --git a/src/clib/ctonedim.cpp b/src/clib/ctonedim.cpp
index 355bd15080f..4896bb9781f 100644
--- a/src/clib/ctonedim.cpp
+++ b/src/clib/ctonedim.cpp
@@ -398,7 +398,7 @@ extern "C" {
         }
     }
 
-    //------------------ stagnation flow domains --------------------
+    //------------------ flow domains --------------------
 
     int flow1D_new(int iph, int ikin, int itr, int itype)
     {
@@ -490,6 +490,77 @@ extern "C" {
         }
     }
 
+    int stflow_new(int iph, int ikin, int itr, int itype)
+    {
+        try {
+            throw NotImplementedError("stflow_new",
+                "Function replaced by 'flow1D_new'.");
+        } catch (...) {
+            return handleAllExceptions(-1, ERR);
+        }
+    }
+
+    int stflow_setTransport(int i, int itr)
+    {
+        try {
+            throw NotImplementedError("stflow_setTransport",
+                "Function replaced by 'flow1D_setTransport'.");
+        } catch (...) {
+            return handleAllExceptions(-1, ERR);
+        }
+    }
+
+    int stflow_enableSoret(int i, int iSoret)
+    {
+        try {
+            throw NotImplementedError("stflow_enableSoret",
+                "Function replaced by 'flow1D_enableSoret'.");
+        } catch (...) {
+            return handleAllExceptions(-1, ERR);
+        }
+    }
+
+    int stflow_setPressure(int i, double p)
+    {
+        try {
+            throw NotImplementedError("stflow_setPressure",
+                "Function replaced by 'flow1D_setPressure'.");
+        } catch (...) {
+            return handleAllExceptions(-1, ERR);
+        }
+    }
+
+    double stflow_pressure(int i)
+    {
+        try {
+            throw NotImplementedError("stflow_pressure",
+                "Function replaced by 'flow1D_pressure'.");
+        } catch (...) {
+            return handleAllExceptions(DERR, DERR);
+        }
+    }
+
+    int stflow_setFixedTempProfile(int i, size_t n, const double* pos,
+                                   size_t m, const double* temp)
+    {
+        try {
+            throw NotImplementedError("stflow_setFixedTempProfile",
+                "Function replaced by 'flow1D_setFixedTempProfile'.");
+        } catch (...) {
+            return handleAllExceptions(-1, ERR);
+        }
+    }
+
+    int stflow_solveEnergyEqn(int i, int flag)
+    {
+        try {
+            throw NotImplementedError("stflow_solveEnergyEqn",
+                "Function replaced by 'flow1D_solveEnergyEqn'.");
+        } catch (...) {
+            return handleAllExceptions(-1, ERR);
+        }
+    }
+
     //------------------- Sim1D --------------------------------------
 
     int sim1D_new(size_t nd, const int* domains)
diff --git a/test/clib/test_ctonedim.cpp b/test/clib/test_ctonedim.cpp
index feff0abda88..cf4752b20d3 100644
--- a/test/clib/test_ctonedim.cpp
+++ b/test/clib/test_ctonedim.cpp
@@ -254,3 +254,25 @@ TEST(ctonedim, freeflame_from_parts)
         Tprev = T;
     }
 }
+
+TEST(ctonedim, stflow_tests)
+{
+    //! @todo: To be removed after Cantera 3.1
+    ct_resetStorage();
+    auto gas = newThermo("h2o2.yaml", "ohmech");
+
+    int sol = soln_newSolution("h2o2.yaml", "ohmech", "default");
+    int ph = soln_thermo(sol);
+    int kin = soln_kinetics(sol);
+    int tr = soln_transport(sol);
+
+    // spot check some errors
+    int itype = 2; // free flow
+    int ret = stflow_new(ph, kin, tr, itype);
+    ASSERT_EQ(ret, -1);  // -1 is an error code
+
+    int flow = flow1D_new(ph, kin, tr, itype);
+    ASSERT_EQ(stflow_setTransport(flow, tr), -1);
+    ASSERT_EQ(stflow_pressure(flow), DERR);  // DERR is an error code
+    ASSERT_EQ(stflow_setPressure(flow, OneAtm), -1);
+}
diff --git a/test/clib/utils.h b/test/clib/utils.h
new file mode 100644
index 00000000000..e69de29bb2d