[Reactor] Automatically clone Solution objects

Author: speth

include/cantera/zeroD/ReactorBase.h

@@ -84,6 +84,9 @@ class ReactorBase
     //!              ReactorBase with Solution object.
     void setSolution(shared_ptr<Solution> sol);
 
+    //! Access the Solution object used to represent the contents of this reactor.
+    shared_ptr<Solution> solution() { return m_solution; }
+
     //! @name Methods to set up a simulation
     //! @{
 

interfaces/cython/cantera/reactor.pxd

@@ -38,6 +38,7 @@ cdef extern from "cantera/zerodim.h" namespace "Cantera":
         CxxReactorBase() except +translate_exception
         string type()
         void setSolution(shared_ptr[CxxSolution]) except +translate_exception
+        shared_ptr[CxxSolution] solution()
         void restoreState() except +translate_exception
         void syncState() except +translate_exception
         double volume()

interfaces/cython/cantera/reactor.pyx

@@ -28,8 +28,7 @@ cdef class ReactorBase:
         self._outlets = []
         self._walls = []
         self._surfaces = []
-        if isinstance(contents, _SolutionBase):
-            self._contents = contents
+        self._contents = _wrap_Solution(self.rbase.solution())
 
         if volume is not None:
             self.volume = volume

samples/cxx/kinetics1/kinetics1.cpp

@@ -45,7 +45,7 @@ int kinetics1(int np, void* p)
     // create a 2D array to hold the output variables,
     // and store the values for the initial state
     Array2D states(nsp+4, 1);
-    saveSoln(0, 0.0, *(sol->thermo()), states);
+    saveSoln(0, 0.0, *r->solution()->thermo(), states);
 
     // create a container object for reactor to run the simulation
     ReactorNet sim(r);
@@ -56,13 +56,13 @@ int kinetics1(int np, void* p)
         double tm = i*dt;
         sim.advance(tm);
         cout << "time = " << tm << " s" << endl;
-        saveSoln(tm, *(sol->thermo()), states);
+        saveSoln(tm, *r->solution()->thermo(), states);
     }
     clock_t t1 = clock(); // save end time
 
 
     // make a CSV output file
-    writeCsv("kin1.csv", *sol->thermo(), states);
+    writeCsv("kin1.csv", *r->solution()->thermo(), states);
 
     // print final temperature and timing data
     double tmm = 1.0*(t1 - t0)/CLOCKS_PER_SEC;

samples/cxx/openmp_ignition/openmp_ignition.cpp

@@ -27,15 +27,14 @@ void run()
     // Containers for Cantera objects to be used in different. Each thread needs
     // to have its own set of linked Cantera objects. Multiple threads accessing
     // the same objects at the same time will cause errors.
-    vector<shared_ptr<Solution>> sols;
     vector<shared_ptr<IdealGasConstPressureReactor>> reactors;
     vector<unique_ptr<ReactorNet>> nets;
 
     // Create and link the Cantera objects for each thread. This step should be
-    // done in serial
+    // done in serial. Only one initial Solution object is needed because a copy will
+    // be created for each new Reactor automatically.
+    auto sol = newSolution("gri30.yaml", "gri30", "none");
     for (int i = 0; i < nThreads; i++) {
-        auto sol = newSolution("gri30.yaml", "gri30", "none");
-        sols.emplace_back(sol);
         reactors.emplace_back(new IdealGasConstPressureReactor(sol));
         nets.emplace_back(new ReactorNet(reactors.back()));
     }
@@ -61,7 +60,7 @@ void run()
     for (int i = 0; i < nPoints; i++) {
         // Get the Cantera objects that were initialized for this thread
         size_t j = omp_get_thread_num();
-        auto gas = sols[j]->thermo();
+        auto gas = reactors[j]->solution()->thermo();
         Reactor& reactor = *reactors[j];
         ReactorNet& net = *nets[j];
 

src/zeroD/Reactor.cpp

@@ -25,7 +25,7 @@ namespace Cantera
 Reactor::Reactor(shared_ptr<Solution> sol, const string& name)
     : ReactorBase(sol, name)
 {
-    m_kin = sol->kinetics().get();
+    m_kin = m_solution->kinetics().get();
     if (m_kin->nReactions() == 0) {
         setChemistry(false);
     } else {

src/zeroD/ReactorBase.cpp

@@ -24,9 +24,9 @@ ReactorBase::ReactorBase(shared_ptr<Solution> sol, const string& name)
         throw CanteraError("ReactorBase::ReactorBase",
                            "Missing or incomplete Solution object.");
     }
-    m_solution = sol;
+    m_solution = sol->clone({}, true, false);
     m_solution->thermo()->addSpeciesLock();
-    setThermo(*sol->thermo());
+    setThermo(*m_solution->thermo());
 }
 
 ReactorBase::~ReactorBase()
@@ -63,9 +63,9 @@ void ReactorBase::setSolution(shared_ptr<Solution> sol)
     }
     m_solution = sol;
     m_solution->thermo()->addSpeciesLock();
-    setThermo(*sol->thermo());
+    setThermo(*m_solution->thermo());
     try {
-        setKinetics(*sol->kinetics());
+        setKinetics(*m_solution->kinetics());
     } catch (NotImplementedError&) {
         // kinetics not used (example: Reservoir)
     }

src/zeroD/ReactorSurface.cpp

@@ -16,12 +16,35 @@ namespace Cantera
 ReactorSurface::ReactorSurface(shared_ptr<Solution> soln,
                                const vector<shared_ptr<ReactorBase>>& reactors,
                                const string& name)
-    : ReactorSurface(soln, name)
+    : ReactorBase(name)
 {
-    for (auto& r : reactors) {
-        r->addSurface(this);
-        m_reactors.push_back(r.get());
+    vector<shared_ptr<Solution>> adjacent;
+    for (auto R : reactors) {
+        adjacent.push_back(R->solution());
+        m_reactors.push_back(R.get());
+        R->addSurface(this);
     }
+    m_solution = soln->clone(adjacent, true, false);
+    m_solution->thermo()->addSpeciesLock();
+    setThermo(*m_solution->thermo());
+    if (!std::dynamic_pointer_cast<SurfPhase>(soln->thermo())) {
+        throw CanteraError("ReactorSurface::ReactorSurface",
+            "Solution object must have a SurfPhase object as the thermo manager.");
+    }
+
+    if (!soln->kinetics() ) {
+        throw CanteraError("ReactorSurface::ReactorSurface",
+            "Solution object must have kinetics manager.");
+    } else if (!std::dynamic_pointer_cast<InterfaceKinetics>(soln->kinetics())) {
+        throw CanteraError("ReactorSurface::ReactorSurface",
+            "Kinetics manager must be an InterfaceKinetics object.");
+    }
+    // todo: move all member variables to use shared pointers after Cantera 3.2
+    m_kinetics = m_solution->kinetics().get();
+    m_thermo = m_solution->thermo().get();
+    m_surf = dynamic_cast<SurfPhase*>(m_thermo);
+    m_cov.resize(m_surf->nSpecies());
+    m_surf->getCoverages(m_cov.data());
 }
 
 ReactorSurface::ReactorSurface(shared_ptr<Solution> sol, const string& name)

test/python/test_reactor.py

@@ -259,7 +259,7 @@ def test_equalize_pressure(self):
         self.net.advance(1.0)
 
         assert self.net.time == approx(1.0)
-        assert self.gas1.P == approx(self.gas2.P)
+        assert self.r1.thermo.P == approx(self.r2.thermo.P)
         assert self.r1.T != approx(self.r2.T)
 
     def test_tolerances(self, rtol_lim=1e-10, atol_lim=1e-20):
@@ -1203,7 +1203,7 @@ def test_nonreacting(self):
         mdot_o = 5.0
         T0 = 900.0
         self.setup_reactor(T0, 10*ct.one_atm, mdot_f, mdot_o)
-        self.gas.set_multiplier(0.0)
+        self.combustor.thermo.set_multiplier(0.0)
         t,T = self.integrate(100.0)
 
         for i in range(len(t)):
@@ -1406,7 +1406,7 @@ def test_with_surface_reactions(self):
         self.net1.rtol = self.net2.rtol = 1e-9
         self.integrate(surf=True)
 
-    def test_surf_install_deprecated(self):
+    def test_surf_install_deprecated(self, allow_deprecated):
         self.create_reactors(add_surf=True, use_surf_install=True)
         self.net1.atol = self.net2.atol = 1e-18
         self.net1.rtol = self.net2.rtol = 1e-9
@@ -1642,9 +1642,10 @@ def update(self, gas):
 
 @ct.extension(name="fail-rate", data=FailRateData)
 class FailRate(ct.ExtensibleRate):
-    def __init__(self, *args, recoverable, **kwargs):
+    def __init__(self, *args, recoverable=None, **kwargs):
         super().__init__(*args, **kwargs)
-        self.recoverable = recoverable
+        if recoverable is not None:
+            self.recoverable = recoverable
         self.count = 0
 
     def eval(self, data):
@@ -1654,6 +1655,12 @@ def eval(self, data):
                 raise ValueError("spam")
         return 0.0
 
+    def get_parameters(self, params):
+        params["recoverable"] = self.recoverable
+
+    def set_parameters(self, params, rate_coeff_units):
+        self.recoverable = params["recoverable"]
+
 
 class TestFlowReactor:
     gas_def = """
@@ -1859,7 +1866,7 @@ def test_recoverable_integrator_errors(self):
             sim.step()
 
         # At least some "recoverable" errors occurred
-        assert fail.rate.count > 0
+        assert r.thermo.reaction(gas.n_reactions - 1).rate.count > 0
 
     def test_max_steps(self):
         surf, gas = self.import_phases()
@@ -1982,8 +1989,8 @@ def test_reinitialization(self):
         cov1 = rsurf.kinetics.coverages
 
         # Reset the reactor to the same initial state
-        gas.TPX = 1700, 4000, 'NH3:1.0, SiF4:0.4'
-        surf.TP = gas.TP
+        r.thermo.TPX = 1700, 4000, 'NH3:1.0, SiF4:0.4'
+        surf.TP = 1700, 4000
         r.mass_flow_rate = 0.01
         r.syncState()
 
@@ -2354,7 +2361,7 @@ def integrate(r, net):
 
             rA, rB, surfX, surfY, net = setup(order)
             for (obj,k) in [(surfY,2), (surfX,2), (rB,18),
-                            (surfX,0), (rB,2)]:
+                            (surfY,0), (rB,2)]:
                 obj.add_sensitivity_reaction(k)
 
             integrate(rB, net)
@@ -2963,7 +2970,7 @@ def after_eval(self,t,LHS,RHS):
         # compare heat added (add_heat) to the equivalent energy contained by the solid
         # and gaseous mass in the reactor
         r1_heat = (mass_lump * cp_lump * (r1.thermo.T - 500) +
-                   mass_gas * (self.gas.enthalpy_mass - gas_initial_enthalpy))
+                   mass_gas * (r1.thermo.enthalpy_mass - gas_initial_enthalpy))
         add_heat = Q * time
         assert add_heat == approx(r1_heat, abs=1e-5)
 
@@ -3051,9 +3058,9 @@ def replace_update_surface_state(self, y):
         Hweight = ct.Element("H").weight
         total_sites = rsurf.area * surf.site_density
         def masses():
-            mass_H = (gas.elemental_mass_fraction("H") * r1.mass +
+            mass_H = (r1.thermo.elemental_mass_fraction("H") * r1.mass +
                       total_sites * r1.surfaces[0].kinetics["H(s)"].X * Hweight)
-            mass_O = gas.elemental_mass_fraction("O") * r1.mass
+            mass_O = r1.thermo.elemental_mass_fraction("O") * r1.mass
             return mass_H, mass_O
 
         net.step()
@@ -3263,11 +3270,11 @@ def test_jacobian(self):
 
         upstream = ct.Reservoir(gas)
         gas.equilibrate("HP")
-        gas.set_multiplier(0.0)
         downstream = ct.Reservoir(gas)
         V0 = 1e-3
         mdot = 120
         r = ct.MoleReactor(gas, volume=V0)
+        r.thermo.set_multiplier(0.0)
         inlet = ct.MassFlowController(upstream, r, mdot=mdot)
         ct.MassFlowController(r, downstream, mdot=mdot)
         net = ct.ReactorNet([r])