From 8aa93625f4cf8ba21a0b54de269152d2221c4882 Mon Sep 17 00:00:00 2001 From: Anthony Walker Date: Tue, 10 Oct 2023 10:47:23 -0700 Subject: [PATCH 01/13] Extend jacobian building to Extensible Reactor system --- include/cantera/base/Delegator.h | 28 +++++++ .../zeroD/IdealGasConstPressureMoleReactor.h | 2 +- include/cantera/zeroD/IdealGasMoleReactor.h | 8 +- include/cantera/zeroD/Reactor.h | 26 +++++- include/cantera/zeroD/ReactorDelegator.h | 8 ++ interfaces/cython/cantera/_utils.pxd | 4 + interfaces/cython/cantera/_utils.pyx | 32 ++++++++ interfaces/cython/cantera/delegator.pxd | 14 +++- interfaces/cython/cantera/delegator.pyx | 17 +++- interfaces/cython/cantera/jacobians.pxd | 2 + interfaces/cython/cantera/jacobians.pyx | 9 +++ interfaces/cython/cantera/reactor.pxd | 2 + interfaces/cython/cantera/reactor.pyx | 3 +- .../IdealGasConstPressureMoleReactor.cpp | 15 ++-- src/zeroD/IdealGasMoleReactor.cpp | 14 +--- src/zeroD/ReactorNet.cpp | 23 ++++-- test/python/test_reactor.py | 80 +++++++++++++++++++ 17 files changed, 246 insertions(+), 41 deletions(-) diff --git a/include/cantera/base/Delegator.h b/include/cantera/base/Delegator.h index 3a470e5e11b..12122f51931 100644 --- a/include/cantera/base/Delegator.h +++ b/include/cantera/base/Delegator.h @@ -10,6 +10,7 @@ #include "cantera/base/Units.h" #include "cantera/base/ctexceptions.h" #include "cantera/base/ExtensionManager.h" +#include "cantera/numerics/eigen_sparse.h" #include #include @@ -226,6 +227,20 @@ class Delegator *m_funcs_v_d_dp_dp[name] = makeDelegate(func, when, *m_funcs_v_d_dp_dp[name]); } + //! Set delegates for member functions with the signature + //! `void(vector>&)` + void setDelegate(const string& name, + const function>&)>& func, + const string& when) + { + if (!m_funcs_v_vet.count(name)) { + throw NotImplementedError("Delegator::setDelegate", + "for function '{}' with signature " + "'void(vector>&)'.", name); + } + *m_funcs_v_vet[name] = makeDelegate(func, when, *m_funcs_v_vet[name]); + } + //! Set delegates for member functions with the signature //! `void(double*, double*, double*)` void setDelegate( @@ -386,6 +401,16 @@ class Delegator m_funcs_v_dp_dp_dp[name] = ⌖ } + //! Install a function with the signature `void(vector>&)` + //! as being delegatable + void install(const string& name, + function>&)>& target, + const function>&)>& base) + { + target = base; + m_funcs_v_vet[name] = ⌖ + } + //! Install a function with the signature `double(void*)` as being delegatable void install(const string& name, function& target, const function& func) @@ -514,6 +539,7 @@ class Delegator //! - `sz` for `size_t` //! - `AM` for `AnyMap` //! - `US` for `UnitStack` + //! - `VET` for `vector>` //! - prefix `c` for `const` arguments //! - suffix `r` for reference arguments //! - suffix `p` for pointer arguments @@ -532,6 +558,7 @@ class Delegator function, double, double*, double*)>*> m_funcs_v_d_dp_dp; map, double*, double*, double*)>*> m_funcs_v_dp_dp_dp; + map>&)>*> m_funcs_v_vet; // Delegates with a return value map> m_base_d_vp; @@ -542,6 +569,7 @@ class Delegator map> m_base_sz_csr; map*> m_funcs_sz_csr; + //! @} //! Handles to wrappers for the delegated object in external language interfaces. diff --git a/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h b/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h index 27de6d57c46..7313d0912b3 100644 --- a/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h +++ b/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h @@ -39,7 +39,7 @@ class IdealGasConstPressureMoleReactor : public ConstPressureMoleReactor //! Neglects derivatives with respect to mole fractions that would generate a //! fully-dense Jacobian. Currently also neglects terms related to interactions //! between reactors, for example via inlets and outlets. - Eigen::SparseMatrix jacobian() override; + void buildJacobian(vector>& jacVector) override; bool preconditionerSupported() const override { return true; }; diff --git a/include/cantera/zeroD/IdealGasMoleReactor.h b/include/cantera/zeroD/IdealGasMoleReactor.h index 9a2fceece4f..9fe2a1ead44 100644 --- a/include/cantera/zeroD/IdealGasMoleReactor.h +++ b/include/cantera/zeroD/IdealGasMoleReactor.h @@ -40,12 +40,8 @@ class IdealGasMoleReactor : public MoleReactor void updateState(double* y) override; - //! Calculate an approximate Jacobian to accelerate preconditioned solvers - - //! Neglects derivatives with respect to mole fractions that would generate a - //! fully-dense Jacobian. Currently, also neglects terms related to interactions - //! between reactors, for example via inlets and outlets. - Eigen::SparseMatrix jacobian() override; + //! Calculate the Jacobian to accelerate preconditioned solvers + void buildJacobian(vector>& jacVector) override; bool preconditionerSupported() const override {return true;}; diff --git a/include/cantera/zeroD/Reactor.h b/include/cantera/zeroD/Reactor.h index 8dcaf5853d2..03b1d509523 100644 --- a/include/cantera/zeroD/Reactor.h +++ b/include/cantera/zeroD/Reactor.h @@ -199,7 +199,7 @@ class Reactor : public ReactorBase //! @param limit value for step size limit void setAdvanceLimit(const string& nm, const double limit); - //! Calculate the Jacobian of a specific Reactor specialization. + //! A wrapper for the Jacobian function to return the Eigen::SparseMatrix //! @warning Depending on the particular implementation, this may return an //! approximate Jacobian intended only for use in forming a preconditioner for //! iterative solvers. @@ -208,9 +208,31 @@ class Reactor : public ReactorBase //! @warning This method is an experimental part of the %Cantera //! API and may be changed or removed without notice. virtual Eigen::SparseMatrix jacobian() { - throw NotImplementedError("Reactor::jacobian"); + m_jac_trips.clear(); + // Add before, during, after evals + buildJacobian(m_jac_trips); + // construct jacobian from vector + Eigen::SparseMatrix jac(m_nv, m_nv); + jac.setFromTriplets(m_jac_trips.begin(), m_jac_trips.end()); + return jac; + } + + //! Calculate the Jacobian of a specific Reactor specialization. + //! @param jac_vector vector where jacobian triplets are added + //! @param offset offset added to the row and col indices of the elements + //! @warning Depending on the particular implementation, this may return an + //! approximate Jacobian intended only for use in forming a preconditioner for + //! iterative solvers. + //! @ingroup derivGroup + //! + //! @warning This method is an experimental part of the %Cantera + //! API and may be changed or removed without notice. + virtual void buildJacobian(vector>& jacVector) { + throw NotImplementedError(type() + "::buildJacobian"); } + // virtual void jacobian() + //! Calculate the reactor-specific Jacobian using a finite difference method. //! //! This method is used only for informational purposes. Jacobian calculations diff --git a/include/cantera/zeroD/ReactorDelegator.h b/include/cantera/zeroD/ReactorDelegator.h index f2556a78717..1790331c100 100644 --- a/include/cantera/zeroD/ReactorDelegator.h +++ b/include/cantera/zeroD/ReactorDelegator.h @@ -10,6 +10,7 @@ #include "cantera/base/Delegator.h" #include "cantera/zeroD/ReactorSurface.h" #include "cantera/thermo/SurfPhase.h" +#include "cantera/numerics/eigen_sparse.h" namespace Cantera { @@ -94,6 +95,8 @@ class ReactorDelegator : public Delegator, public R, public ReactorAccessor [this](const string& nm) { return R::componentIndex(nm); }); install("speciesIndex", m_speciesIndex, [this](const string& nm) { return R::speciesIndex(nm); }); + install("buildJacobian", m_build_jacobian, + [this](vector>& jv) { R::buildJacobian(jv); }); } // Overrides of Reactor methods @@ -160,6 +163,10 @@ class ReactorDelegator : public Delegator, public R, public ReactorAccessor return m_speciesIndex(nm); } + void buildJacobian(vector>& jacVector) override { + m_build_jacobian(jacVector); + } + // Public access to protected Reactor variables needed by derived classes void setNEq(size_t n) override { @@ -204,6 +211,7 @@ class ReactorDelegator : public Delegator, public R, public ReactorAccessor function m_componentName; function m_componentIndex; function m_speciesIndex; + function>&)> m_build_jacobian; }; } diff --git a/interfaces/cython/cantera/_utils.pxd b/interfaces/cython/cantera/_utils.pxd index ade43f1e140..77f98c958b9 100644 --- a/interfaces/cython/cantera/_utils.pxd +++ b/interfaces/cython/cantera/_utils.pxd @@ -9,6 +9,7 @@ from libcpp.memory cimport unique_ptr, make_unique from .ctcxx cimport * from .units cimport UnitSystem, CxxUnits +from .delegator cimport CxxEigenTriplet cdef extern from "cantera/base/AnyMap.h" namespace "Cantera": cdef cppclass CxxAnyValue "Cantera::AnyValue" @@ -118,3 +119,6 @@ cdef anymap_to_py(CxxAnyMap& m) cdef CxxAnyValue python_to_anyvalue(item, name=*) except * cdef anyvalue_to_python(string name, CxxAnyValue& v) + +cdef vector[CxxEigenTriplet] python_to_triplets(triplets) except * +cdef triplets_to_python(vector[CxxEigenTriplet]& triplets) diff --git a/interfaces/cython/cantera/_utils.pyx b/interfaces/cython/cantera/_utils.pyx index ef423ba45bb..956558788a7 100644 --- a/interfaces/cython/cantera/_utils.pyx +++ b/interfaces/cython/cantera/_utils.pyx @@ -5,6 +5,7 @@ import os import warnings from cpython.ref cimport PyObject from libcpp.utility cimport move +from cython.operator cimport dereference import numbers as _numbers import importlib as _importlib from collections import namedtuple as _namedtuple @@ -527,6 +528,37 @@ cdef vector[vector[string]] list2_string_to_anyvalue(data): v[i][j] = stringify(jtem) return v +cdef vector[CxxEigenTriplet] python_to_triplets(triplets): + # check that triplet dimensions are met + trip_shape = np.shape(triplets) + assert len(trip_shape) == 2 + assert trip_shape[1] == 3 + cdef vector[CxxEigenTriplet] trips + # c++ variables + cdef size_t row + cdef size_t col + cdef double val + cdef CxxEigenTriplet* trip_ptr + for r, c, v in triplets: + row = r + col = c + val = v + trip_ptr = new CxxEigenTriplet(row, col, val) + trips.push_back(dereference(trip_ptr)) + return trips + +cdef triplets_to_python(vector[CxxEigenTriplet]& triplets): + values = [] + cdef size_t row + cdef size_t col + cdef double val + for t in triplets: + row = t.row() + col = t.col() + val = t.value() + values.append((row, col, val)) + return values + def _py_to_any_to_py(dd): # used for internal testing purposes only cdef string name = stringify("test") diff --git a/interfaces/cython/cantera/delegator.pxd b/interfaces/cython/cantera/delegator.pxd index 721ddd736a3..09843a3faa0 100644 --- a/interfaces/cython/cantera/delegator.pxd +++ b/interfaces/cython/cantera/delegator.pxd @@ -7,6 +7,15 @@ from .ctcxx cimport * from .func1 cimport * from .units cimport CxxUnitStack +# from .kinetics cimport * + +cdef extern from "cantera/numerics/eigen_sparse.h" namespace "Eigen": + cdef cppclass CxxEigenTriplet "Eigen::Triplet": + CxxEigenTriplet() + CxxEigenTriplet(size_t, size_t, double) + size_t row() + size_t col() + size_t value() cdef extern from "" namespace "std" nogil: cdef cppclass size_array1 "std::array": @@ -52,7 +61,7 @@ cdef extern from "cantera/base/Delegator.h" namespace "Cantera": void setDelegate(string&, function[int(double&, void*)], string&) except +translate_exception void setDelegate(string&, function[int(string&, size_t)], string&) except +translate_exception void setDelegate(string&, function[int(size_t&, string&)], string&) except +translate_exception - + void setDelegate(string&, function[void(vector[CxxEigenTriplet]&)], string&) except +translate_exception cdef extern from "cantera/cython/funcWrapper.h": # pyOverride is actually a templated function, but we have to specify the individual @@ -77,6 +86,9 @@ cdef extern from "cantera/cython/funcWrapper.h": cdef function[int(string&, size_t)] pyOverride(PyObject*, int(PyFuncInfo&, string&, size_t)) cdef function[int(size_t&, const string&)] pyOverride( PyObject*, int(PyFuncInfo&, size_t&, const string&)) + cdef function[void(vector[CxxEigenTriplet]&)] pyOverride( + PyObject*, void(PyFuncInfo&, vector[CxxEigenTriplet]&)) + cdef extern from "cantera/base/ExtensionManager.h" namespace "Cantera": cdef cppclass CxxExtensionManager "Cantera::ExtensionManager": diff --git a/interfaces/cython/cantera/delegator.pyx b/interfaces/cython/cantera/delegator.pyx index c5b01bfd531..f209b18e275 100644 --- a/interfaces/cython/cantera/delegator.pyx +++ b/interfaces/cython/cantera/delegator.pyx @@ -9,7 +9,7 @@ from libc.stdlib cimport malloc from libc.string cimport strcpy from ._utils import CanteraError -from ._utils cimport stringify, pystr, anymap_to_py, py_to_anymap +from ._utils cimport stringify, pystr, anymap_to_py, py_to_anymap, triplets_to_python, python_to_triplets from .units cimport Units, UnitStack # from .reaction import ExtensibleRate, ExtensibleRateData from .reaction cimport (ExtensibleRate, ExtensibleRateData, CxxReaction, @@ -244,6 +244,18 @@ cdef void callback_v_d_dp_dp(PyFuncInfo& funcInfo, size_array2 sizes, double arg funcInfo.setExceptionType(exc_type) funcInfo.setExceptionValue(exc_value) +# Wrapper for void(vector&) +cdef void callback_v_vet(PyFuncInfo& funcInfo, vector[CxxEigenTriplet]& jac_vector) noexcept: + try: + python_trips = triplets_to_python(jac_vector) + # convert vector to python object + (funcInfo.func())(python_trips) + jac_vector = python_to_triplets(python_trips) + except BaseException as e: + exc_type, exc_value = sys.exc_info()[:2] + funcInfo.setExceptionType(exc_type) + funcInfo.setExceptionValue(exc_value) + cdef int assign_delegates(obj, CxxDelegator* delegator) except -1: """ Use methods defined in the Python class ``obj`` as delegates for the C++ @@ -305,7 +317,6 @@ cdef int assign_delegates(obj, CxxDelegator* delegator) except -1: if when is None: continue - cxx_name = stringify(options[0]) callback = options[1].replace(' ', '') @@ -366,6 +377,8 @@ cdef int assign_delegates(obj, CxxDelegator* delegator) except -1: elif callback == 'void(double,double*,double*)': delegator.setDelegate(cxx_name, pyOverride(method, callback_v_d_dp_dp), cxx_when) + elif callback == 'void(vector[CxxEigenTriplet]&)': + delegator.setDelegate(cxx_name, pyOverride(method, callback_v_vet), cxx_when) else: raise ValueError("Don't know how to set delegates for functions " f"with signature '{callback}'") diff --git a/interfaces/cython/cantera/jacobians.pxd b/interfaces/cython/cantera/jacobians.pxd index a5c6c45b18e..3a5cd3032b3 100644 --- a/interfaces/cython/cantera/jacobians.pxd +++ b/interfaces/cython/cantera/jacobians.pxd @@ -12,6 +12,8 @@ cdef extern from "cantera/numerics/SystemJacobian.h" namespace "Cantera": CxxSystemJacobian() string preconditionerSide() string type() + double gamma() + void setGamma(double) except +translate_exception void setPreconditionerSide(string) except +translate_exception cdef extern from "cantera/numerics/EigenSparseJacobian.h" namespace "Cantera": diff --git a/interfaces/cython/cantera/jacobians.pyx b/interfaces/cython/cantera/jacobians.pyx index bc951af06f1..432d2489b09 100644 --- a/interfaces/cython/cantera/jacobians.pyx +++ b/interfaces/cython/cantera/jacobians.pyx @@ -63,6 +63,15 @@ cdef class SystemJacobian: def __set__(self, side): self._base.get().setPreconditionerSide(stringify(side)) + property gamma: + """ Get/Set the value of gamma used in the expression P = (I - gamma * J). + """ + def __get__(self): + return self.pbase.get().gamma() + + def __set__(self, value): + self.pbase.get().setGamma(value) + cdef class EigenSparseJacobian(SystemJacobian): """ diff --git a/interfaces/cython/cantera/reactor.pxd b/interfaces/cython/cantera/reactor.pxd index 2db654e3365..6c19b9f4603 100644 --- a/interfaces/cython/cantera/reactor.pxd +++ b/interfaces/cython/cantera/reactor.pxd @@ -5,6 +5,7 @@ #distutils: language = c++ from .ctcxx cimport * +from .delegator cimport CxxEigenTriplet from .kinetics cimport * from .func1 cimport * from .jacobians cimport * @@ -58,6 +59,7 @@ cdef extern from "cantera/zerodim.h" namespace "Cantera": void getState(double*) except +translate_exception CxxSparseMatrix jacobian() except +translate_exception CxxSparseMatrix finiteDifferenceJacobian() except +translate_exception + void buildJacobian(vector[CxxEigenTriplet]&) except +translate_exception void setAdvanceLimit(string&, double) except +translate_exception void addSensitivitySpeciesEnthalpy(size_t) except +translate_exception diff --git a/interfaces/cython/cantera/reactor.pyx b/interfaces/cython/cantera/reactor.pyx index 921cbeebb0b..66e83868b28 100644 --- a/interfaces/cython/cantera/reactor.pyx +++ b/interfaces/cython/cantera/reactor.pyx @@ -645,7 +645,8 @@ cdef class ExtensibleReactor(Reactor): 'eval_surfaces': ('evalSurfaces', 'void(double*,double*,double*)'), 'component_name': ('componentName', 'string(size_t)'), 'component_index': ('componentIndex', 'size_t(string)'), - 'species_index': ('speciesIndex', 'size_t(string)') + 'species_index': ('speciesIndex', 'size_t(string)'), + 'build_jacobian': ('buildJacobian', 'void(vector[CxxEigenTriplet]&)') } def __cinit__(self, *args, **kwargs): diff --git a/src/zeroD/IdealGasConstPressureMoleReactor.cpp b/src/zeroD/IdealGasConstPressureMoleReactor.cpp index 648dea27489..42538b3df67 100644 --- a/src/zeroD/IdealGasConstPressureMoleReactor.cpp +++ b/src/zeroD/IdealGasConstPressureMoleReactor.cpp @@ -114,14 +114,13 @@ void IdealGasConstPressureMoleReactor::eval(double time, double* LHS, double* RH } } -Eigen::SparseMatrix IdealGasConstPressureMoleReactor::jacobian() +void IdealGasConstPressureMoleReactor::buildJacobian( + vector>& jacVector) { if (m_nv == 0) { throw CanteraError("IdealGasConstPressureMoleReactor::jacobian", "Reactor must be initialized first."); } - // clear former jacobian elements - m_jac_trips.clear(); // dnk_dnj represents d(dot(n_k)) / d (n_j) but is first assigned as // d (dot(omega)) / d c_j, it is later transformed appropriately. Eigen::SparseMatrix dnk_dnj = m_kin->netProductionRates_ddCi(); @@ -169,7 +168,7 @@ Eigen::SparseMatrix IdealGasConstPressureMoleReactor::jacobian() if (static_cast(it.row()) < m_nsp) { it.valueRef() = it.value() + netProductionRates[it.row()] * molarVol; } - m_jac_trips.emplace_back(static_cast(it.row() + m_sidx), + jacVector.emplace_back(static_cast(it.row() + m_sidx), static_cast(it.col() + m_sidx), it.value()); } } @@ -198,7 +197,7 @@ Eigen::SparseMatrix IdealGasConstPressureMoleReactor::jacobian() for (size_t j = 0; j < m_nv; j++) { double ydotPerturbed = rhsPerturbed[j] / lhsPerturbed[j]; double ydotCurrent = rhsCurrent[j] / lhsCurrent[j]; - m_jac_trips.emplace_back(static_cast(j), 0, + jacVector.emplace_back(static_cast(j), 0, (ydotPerturbed - ydotCurrent) / deltaTemp); } // d T_dot/dnj @@ -223,14 +222,10 @@ Eigen::SparseMatrix IdealGasConstPressureMoleReactor::jacobian() Eigen::VectorXd hk_dnkdnj_sums = dnk_dnj.transpose() * enthalpy; // Add derivatives to jac by spanning columns for (size_t j = 0; j < ssize; j++) { - m_jac_trips.emplace_back(0, static_cast(j + m_sidx), + jacVector.emplace_back(0, static_cast(j + m_sidx), (specificHeat[j] * qdot - NCp * hk_dnkdnj_sums[j]) * denom); } } - // convert triplets to sparse matrix - Eigen::SparseMatrix jac(m_nv, m_nv); - jac.setFromTriplets(m_jac_trips.begin(), m_jac_trips.end()); - return jac; } size_t IdealGasConstPressureMoleReactor::componentIndex(const string& nm) const diff --git a/src/zeroD/IdealGasMoleReactor.cpp b/src/zeroD/IdealGasMoleReactor.cpp index d7d39a7140c..9292e8d4afe 100644 --- a/src/zeroD/IdealGasMoleReactor.cpp +++ b/src/zeroD/IdealGasMoleReactor.cpp @@ -183,14 +183,12 @@ void IdealGasMoleReactor::eval(double time, double* LHS, double* RHS) } } -Eigen::SparseMatrix IdealGasMoleReactor::jacobian() +void IdealGasMoleReactor::buildJacobian(vector>& jacVector) { if (m_nv == 0) { throw CanteraError("IdealGasMoleReactor::jacobian", "Reactor must be initialized first."); } - // clear former jacobian elements - m_jac_trips.clear(); // dnk_dnj represents d(dot(n_k)) / d (n_j) but is first assigned as // d (dot(omega)) / d c_j, it is later transformed appropriately. Eigen::SparseMatrix dnk_dnj = m_kin->netProductionRates_ddCi(); @@ -215,7 +213,7 @@ Eigen::SparseMatrix IdealGasMoleReactor::jacobian() // as it substantially reduces matrix sparsity for (int k = 0; k < dnk_dnj.outerSize(); k++) { for (Eigen::SparseMatrix::InnerIterator it(dnk_dnj, k); it; ++it) { - m_jac_trips.emplace_back(static_cast(it.row() + m_sidx), + jacVector.emplace_back(static_cast(it.row() + m_sidx), static_cast(it.col() + m_sidx), it.value()); } } @@ -243,7 +241,7 @@ Eigen::SparseMatrix IdealGasMoleReactor::jacobian() for (size_t j = 0; j < m_nv; j++) { double ydotPerturbed = rhsPerturbed[j] / lhsPerturbed[j]; double ydotCurrent = rhsCurrent[j] / lhsCurrent[j]; - m_jac_trips.emplace_back(static_cast(j), 0, + jacVector.emplace_back(static_cast(j), 0, (ydotPerturbed - ydotCurrent) / deltaTemp); } // d T_dot/dnj @@ -272,14 +270,10 @@ Eigen::SparseMatrix IdealGasMoleReactor::jacobian() Eigen::VectorXd uk_dnkdnj_sums = dnk_dnj.transpose() * internal_energy; // add derivatives to jacobian for (size_t j = 0; j < ssize; j++) { - m_jac_trips.emplace_back(0, static_cast(j + m_sidx), + jacVector.emplace_back(0, static_cast(j + m_sidx), (specificHeat[j] * qdot - NCv * uk_dnkdnj_sums[j]) * denom); } } - // convert triplets to sparse matrix - Eigen::SparseMatrix jac(m_nv, m_nv); - jac.setFromTriplets(m_jac_trips.begin(), m_jac_trips.end()); - return jac; } } diff --git a/src/zeroD/ReactorNet.cpp b/src/zeroD/ReactorNet.cpp index 28642a0cc70..1a4ea0fe906 100644 --- a/src/zeroD/ReactorNet.cpp +++ b/src/zeroD/ReactorNet.cpp @@ -769,18 +769,25 @@ void ReactorNet::preconditionerSetup(double t, double* y, double gamma) vector yCopy(m_nv); // Get state of reactor getState(yCopy.data()); - // transform state based on preconditioner rules + // Transform state based on preconditioner rules precon->stateAdjustment(yCopy); - // update network with adjusted state + // Update network with adjusted state updateState(yCopy.data()); + // Create jacobian triplet vector + vector> jac_vector; + vector jstarts; // Get jacobians and give elements to preconditioners + jstarts.push_back(jac_vector.size()); for (size_t i = 0; i < m_reactors.size(); i++) { - Eigen::SparseMatrix rJac = m_reactors[i]->jacobian(); - for (int k=0; k::InnerIterator it(rJac, k); it; ++it) { - precon->setValue(it.row() + m_start[i], it.col() + m_start[i], - it.value()); - } + m_reactors[i]->buildJacobian(jac_vector); + jstarts.push_back(jac_vector.size()); + } + // Add to preconditioner with offset + for (size_t i=0; i < m_reactors.size(); i++) { + for (size_t j = jstarts[i]; j < jstarts[i+1]; j++) { + auto it = jac_vector[j]; + precon->setValue(it.row() + m_start[i], it.col() + m_start[i], + it.value()); } } // post reactor setup operations diff --git a/test/python/test_reactor.py b/test/python/test_reactor.py index 736e95ee8a9..59881785c32 100644 --- a/test/python/test_reactor.py +++ b/test/python/test_reactor.py @@ -3357,6 +3357,86 @@ def deltaC(): assert r1.phase["H2"].Y[0] == approx(0.13765976, rel=1e-6) assert r2.phase["O2"].Y[0] == approx(0.94617029, rel=1e-6) + def test_after_jacobian(self): + class AfterJacobianReactor(ct.ExtensibleIdealGasMoleReactor): + def __init__(self, *args, neighbor, **kwargs): + super().__init__(*args, **kwargs) + self.v_wall = 0 + self.k_wall = 1e-5 + self.neighbor = neighbor + + def after_initialize(self, t0): + self.n_vars += 1 + self.i_wall = self.n_vars - 1 + + def after_get_state(self, y): + y[self.i_wall] = self.v_wall + + def after_update_state(self, y): + self.v_wall = y[self.i_wall] + self.walls[0].velocity = self.v_wall + + def after_eval(self, t, LHS, RHS): + # Extra equation is d(v_wall)/dt = k * delta P + a = self.k_wall * (self.phase.P - self.neighbor.phase.P) + RHS[self.i_wall] = a + + def before_component_index(self, name): + if name == 'v_wall': + return self.i_wall + + def before_component_name(self, i): + if i == self.i_wall: + return 'v_wall' + + def after_build_jacobian(self, jac_vector): + jac_vector.append((self.i_wall, self.i_wall, 1e20)) + + self.gas.TP = 300, ct.one_atm + res = ct.Reservoir(self.gas) + self.gas.TP = 300, 2 * ct.one_atm + r = AfterJacobianReactor(self.gas, neighbor=res) + w = ct.Wall(r, res) + net = ct.ReactorNet([r]) + precon = ct.AdaptivePreconditioner() + net.preconditioner = precon + net.step() + # test that jacobian wall element is hard coded value + jac = r.jacobian + assert jac[r.i_wall, r.i_wall] == 1e20 + pmat = precon.matrix + assert pmat[r.i_wall, r.i_wall] == (1 - precon.gamma * 1e20) + + def test_before_jacobian(self): + class BeforeJacobianReactor(ct.ExtensibleIdealGasMoleReactor): + + def before_build_jacobian(self, jac_vector): + jac_vector.append((0, 0, 1e10)) + + self.gas.TP = 300, ct.one_atm + r = BeforeJacobianReactor(self.gas) + net = ct.ReactorNet([r]) + net.preconditioner = ct.AdaptivePreconditioner() + net.step() + # test that jacobian wall element is hard coded value + jac = r.jacobian + assert jac[0, 0] == 1e10 + + def test_replace_jacobian(self): + class ReplaceJacobianReactor(ct.ExtensibleIdealGasMoleReactor): + + def replace_build_jacobian(self, jac_vector): + jac_vector.append((0, 0, 0)) + + self.gas.TP = 300, ct.one_atm + r = ReplaceJacobianReactor(self.gas) + net = ct.ReactorNet([r]) + net.preconditioner = ct.AdaptivePreconditioner() + net.step() + # test that jacobian wall element is hard coded value + jac = r.jacobian + assert np.sum(jac) == 0 + class TestSteadySolver: @pytest.mark.parametrize("reactor_class", From 8363ec441921d5a8c9c1521f5a06266260e1999c Mon Sep 17 00:00:00 2001 From: Anthony Walker Date: Wed, 11 Oct 2023 20:45:42 -0700 Subject: [PATCH 02/13] Complete finite difference jacobian --- include/cantera/zeroD/ReactorNet.h | 37 +++++++++++++ src/zeroD/ReactorNet.cpp | 84 +++++++++++++++++++++++++----- test/python/test_reactor.py | 11 ++++ 3 files changed, 119 insertions(+), 13 deletions(-) diff --git a/include/cantera/zeroD/ReactorNet.h b/include/cantera/zeroD/ReactorNet.h index 2a885a3406c..1925a9eb1c3 100644 --- a/include/cantera/zeroD/ReactorNet.h +++ b/include/cantera/zeroD/ReactorNet.h @@ -368,12 +368,49 @@ class ReactorNet : public FuncEval //! `1 - max_i(|y[i]-y_base[i]| / limit[i])` so a zero indicates a component has //! reached its limit; otherwise `gout[0]` is positive. void evalRootFunctions(double t, const double* y, double* gout) override; + //! Calculate the system Jacobian using a finite difference method. + //! + //! This method is used only for informational purposes. Jacobian calculations + //! for the full reactor system are handled internally by CVODES. + //! + //! @warning This method is an experimental part of the %Cantera + //! API and may be changed or removed without notice. + virtual Eigen::SparseMatrix finiteDifferenceJacobian(); + + //! A wrapper method to calculate the system jacobian as Eigen::SparseMatrix + //! @warning Depending on the particular implementation, this may return an + //! approximate Jacobian intended only for use in forming a preconditioner for + //! iterative solvers. + //! + //! @warning This method is an experimental part of the %Cantera + //! API and may be changed or removed without notice. + virtual Eigen::SparseMatrix jacobian() { + vector> jac_trips; + // Add before, during, after evals + buildJacobian(jac_trips); + // construct jacobian from vector + Eigen::SparseMatrix jac(m_nv, m_nv); + jac.setFromTriplets(jac_trips.begin(), jac_trips.end()); + return jac; + } protected: //! Add the reactor *r* to this reactor network. //! @since Changed in %Cantera 3.2. Previous version used a reference. void addReactor(shared_ptr reactor); + //! Calculate the Jacobian of the entire reactor network. + //! @param jac_vector vector where jacobian triplets are added + //! @param offset offset added to the row and col indices of the elements + //! @warning Depending on the particular implementation, this may return an + //! approximate Jacobian intended only for use in forming a preconditioner for + //! iterative solvers. + //! @ingroup derivGroup + //! + //! @warning This method is an experimental part of the %Cantera + //! API and may be changed or removed without notice. + virtual void buildJacobian(vector>& jacVector); + //! Check that preconditioning is supported by all reactors in the network virtual void checkPreconditionerSupported() const; diff --git a/src/zeroD/ReactorNet.cpp b/src/zeroD/ReactorNet.cpp index 1a4ea0fe906..90f0699f11b 100644 --- a/src/zeroD/ReactorNet.cpp +++ b/src/zeroD/ReactorNet.cpp @@ -775,20 +775,10 @@ void ReactorNet::preconditionerSetup(double t, double* y, double gamma) updateState(yCopy.data()); // Create jacobian triplet vector vector> jac_vector; - vector jstarts; - // Get jacobians and give elements to preconditioners - jstarts.push_back(jac_vector.size()); - for (size_t i = 0; i < m_reactors.size(); i++) { - m_reactors[i]->buildJacobian(jac_vector); - jstarts.push_back(jac_vector.size()); - } + buildJacobian(jac_vector); // Add to preconditioner with offset - for (size_t i=0; i < m_reactors.size(); i++) { - for (size_t j = jstarts[i]; j < jstarts[i+1]; j++) { - auto it = jac_vector[j]; - precon->setValue(it.row() + m_start[i], it.col() + m_start[i], - it.value()); - } + for (auto it : jac_vector) { + precon->setValue(it.row(), it.col(), it.value()); } // post reactor setup operations precon->updatePreconditioner(); @@ -817,6 +807,74 @@ void ReactorNet::checkPreconditionerSupported() const { } } +void ReactorNet::buildJacobian(vector>& jacVector) +{ + // network must be initialized for the jacobian + if (! m_init) { + initialize(); + } + // Create jacobian triplet vector + vector jstarts; + // Get jacobians and give elements to preconditioners + jstarts.push_back(jacVector.size()); + for (size_t i = 0; i < m_reactors.size(); i++) { + m_reactors[i]->buildJacobian(jacVector); + jstarts.push_back(jacVector.size()); + } + // Add to preconditioner with offset + for (size_t i=0; i < m_reactors.size(); i++) { + for (size_t j = jstarts[i]; j < jstarts[i+1]; j++) { + auto it = jacVector[j]; + auto newTrip = Eigen::Triplet(it.row() + m_start[i], it.col() + + m_start[i], it.value()); + jacVector[j] = newTrip; + } + } +} + +Eigen::SparseMatrix ReactorNet::finiteDifferenceJacobian() +{ + // network must be initialized for the jacobian + if (! m_init) { + initialize(); + } + + // clear former jacobian elements + vector> jac_trips; + + // Get the current state + Eigen::ArrayXd yCurrent(m_nv); + getState(yCurrent.data()); + + Eigen::ArrayXd yPerturbed = yCurrent; + Eigen::ArrayXd ydotCurrent(m_nv), ydotPerturbed(m_nv); + + eval(m_time, yCurrent.data(), ydotCurrent.data(), m_sens_params.data()); + double rel_perturb = std::sqrt(std::numeric_limits::epsilon()); + + for (size_t j = 0; j < m_nv; j++) { + yPerturbed = yCurrent; + double delta_y = std::max(std::abs(yCurrent[j]), 1000 * m_atols) * rel_perturb; + yPerturbed[j] += delta_y; + ydotPerturbed = 0; + eval(m_time, yPerturbed.data(), ydotPerturbed.data(), m_sens_params.data()); + // d ydot_i/dy_j + for (size_t i = 0; i < m_nv; i++) { + if (ydotCurrent[i] != ydotPerturbed[i]) { + jac_trips.emplace_back( + static_cast(i), static_cast(j), + (ydotPerturbed[i] - ydotCurrent[i]) / delta_y); + } + } + } + updateState(yCurrent.data()); + + Eigen::SparseMatrix jac(m_nv, m_nv); + jac.setFromTriplets(jac_trips.begin(), jac_trips.end()); + return jac; +} + + SteadyReactorSolver::SteadyReactorSolver(ReactorNet* net, double* x0) : m_net(net) { diff --git a/test/python/test_reactor.py b/test/python/test_reactor.py index 59881785c32..8204fb56bf2 100644 --- a/test/python/test_reactor.py +++ b/test/python/test_reactor.py @@ -1543,6 +1543,17 @@ class TestIdealGasMoleReactor(TestMoleReactor): reactorClass = ct.IdealGasMoleReactor test_preconditioner_unsupported = None + @pytest.mark.diagnose + def test_heat_transfer_network(self): + # Result should be the same if (m * cp) / (U * A) is held constant + self.make_reactors(T1=300, T2=1000) + self.add_wall(U=200, A=1.0) + self.net.preconditioner = ct.AdaptivePreconditioner() + print(self.net.finite_difference_jacobian) + # self.net.advance(1.0) + # T1a = self.r1.T + # T2a = self.r2.T + def test_adaptive_precon_integration(self): # Network one with non-mole reactor T0 = 900 From c90fd415be816f8b481b4a164cd6dcd9d5ae9a1c Mon Sep 17 00:00:00 2001 From: Anthony Walker Date: Tue, 17 Oct 2023 10:33:12 -0700 Subject: [PATCH 03/13] Build jacobian structure defined for walls and flow devices --- include/cantera/zeroD/FlowDevice.h | 43 +++++++++++++ .../zeroD/IdealGasConstPressureMoleReactor.h | 12 ++++ include/cantera/zeroD/IdealGasMoleReactor.h | 12 ++++ include/cantera/zeroD/MoleReactor.h | 6 ++ include/cantera/zeroD/Reactor.h | 26 +++++++- include/cantera/zeroD/ReactorBase.h | 58 +++++++++++++++++ include/cantera/zeroD/ReactorNet.h | 18 ++++++ include/cantera/zeroD/Wall.h | 49 ++++++++++++++ interfaces/cython/cantera/reactor.pxd | 2 + interfaces/cython/cantera/reactor.pyx | 25 ++++++++ .../IdealGasConstPressureMoleReactor.cpp | 26 ++++++++ src/zeroD/IdealGasMoleReactor.cpp | 26 ++++++++ src/zeroD/Reactor.cpp | 33 ++++++++++ src/zeroD/ReactorNet.cpp | 55 ++++++++++++++++ src/zeroD/Wall.cpp | 64 +++++++++++++++++++ test/python/test_reactor.py | 64 ++++++++++++++++--- 16 files changed, 509 insertions(+), 10 deletions(-) diff --git a/include/cantera/zeroD/FlowDevice.h b/include/cantera/zeroD/FlowDevice.h index 65630bc8c0a..b9b5b34ad5d 100644 --- a/include/cantera/zeroD/FlowDevice.h +++ b/include/cantera/zeroD/FlowDevice.h @@ -10,6 +10,7 @@ #include "cantera/base/global.h" #include "cantera/base/ctexceptions.h" #include "ConnectorNode.h" +#include "cantera/numerics/eigen_sparse.h" namespace Cantera { @@ -140,7 +141,49 @@ class FlowDevice : public ConnectorNode m_time = time; } + /*! Build the Jacobian terms specific to the flow device for the given connected reactor. + * @param r a pointer to the calling reactor + * @param jacVector a vector of triplets to be added to the jacobian for the reactor + * @warning This function is an experimental part of the %Cantera API and may be changed + * or removed without notice. + * @since New in %Cantera 3.0. + */ + virtual void buildReactorJacobian(ReactorBase* r, vector>& jacVector) { + throw NotImplementedError(type() + "::buildReactorJacobian"); + } + + /*! Build the Jacobian terms specific to the flow device for the network. These terms + * will be adjusted to the networks indexing system outside of the reactor. + * @param jacVector a vector of triplets to be added to the jacobian for the reactor + * @warning This function is an experimental part of the %Cantera API and may be changed + * or removed without notice. + * @since New in %Cantera 3.0. + */ + virtual void buildNetworkJacobian(vector>& jacVector) { + throw NotImplementedError(type() + "::buildNetworkJacobian"); + } + + /*! Specify the jacobian terms have been calculated and should not be recalculated. + * @warning This function is an experimental part of the %Cantera API and may be changed + * or removed without notice. + * @since New in %Cantera 3.0. + */ + void jacobianCalculated() { m_jac_calculated = true; }; + + /*! Specify that jacobian terms have not been calculated and should be recalculated. + * @warning This function is an experimental part of the %Cantera API and may be changed + * or removed without notice. + * @since New in %Cantera 3.0. + */ + void jacobianNotCalculated() { m_jac_calculated = false; }; + protected: + string m_name; //!< Flow device name. + bool m_defaultNameSet = false; //!< `true` if default name has been previously set. + //! a variable to switch on and off so calculations are not doubled by the calling + //! reactor or network + bool m_jac_calculated = false; + double m_mdot = Undef; //! Function set by setPressureFunction; used by updateMassFlowRate diff --git a/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h b/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h index 7313d0912b3..2bb92fa35f1 100644 --- a/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h +++ b/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h @@ -43,6 +43,18 @@ class IdealGasConstPressureMoleReactor : public ConstPressureMoleReactor bool preconditionerSupported() const override { return true; }; + double temperatureDerivative() override { return 1; }; + + double temperatureRadiationDerivative() override { + return 4 * std::pow(temperature(), 3); + } + + double moleDerivative(size_t index) override; + + double moleRadiationDerivative(size_t index) override; + + size_t speciesOffset() const override { return m_sidx; }; + size_t componentIndex(const string& nm) const override; string componentName(size_t k) override; double upperBound(size_t k) const override; diff --git a/include/cantera/zeroD/IdealGasMoleReactor.h b/include/cantera/zeroD/IdealGasMoleReactor.h index 9fe2a1ead44..cd30ddf02d4 100644 --- a/include/cantera/zeroD/IdealGasMoleReactor.h +++ b/include/cantera/zeroD/IdealGasMoleReactor.h @@ -45,6 +45,18 @@ class IdealGasMoleReactor : public MoleReactor bool preconditionerSupported() const override {return true;}; + double temperatureDerivative() override { return 1; }; + + double temperatureRadiationDerivative() override { + return 4 * std::pow(temperature(), 3); + } + + double moleDerivative(size_t index) override; + + double moleRadiationDerivative(size_t index) override; + + size_t speciesOffset() const override { return m_sidx; }; + protected: vector m_uk; //!< Species molar internal energies }; diff --git a/include/cantera/zeroD/MoleReactor.h b/include/cantera/zeroD/MoleReactor.h index e15e07695d2..723d4b29df4 100644 --- a/include/cantera/zeroD/MoleReactor.h +++ b/include/cantera/zeroD/MoleReactor.h @@ -7,6 +7,7 @@ #define CT_MOLEREACTOR_H #include "Reactor.h" +#include "cantera/zeroD/Wall.h" namespace Cantera { @@ -40,6 +41,8 @@ class MoleReactor : public Reactor double lowerBound(size_t k) const override; void resetBadValues(double* y) override; + size_t energyIndex() const override { return m_eidx; }; + protected: //! For each surface in the reactor, update vector of triplets with all relevant //! surface jacobian derivatives of species with respect to species @@ -62,6 +65,9 @@ class MoleReactor : public Reactor //! const value for the species start index const size_t m_sidx = 2; + + //! index of state variable associated with energy + const size_t m_eidx = 0; }; } diff --git a/include/cantera/zeroD/Reactor.h b/include/cantera/zeroD/Reactor.h index 03b1d509523..20bc24c1232 100644 --- a/include/cantera/zeroD/Reactor.h +++ b/include/cantera/zeroD/Reactor.h @@ -231,7 +231,27 @@ class Reactor : public ReactorBase throw NotImplementedError(type() + "::buildJacobian"); } - // virtual void jacobian() + //! Calculate the Jacobian of a Reactor specialization for wall contributions. + //! @param jacVector vector where jacobian triplets are added + //! @warning Depending on the particular implementation, this may return an + //! approximate Jacobian intended only for use in forming a preconditioner for + //! iterative solvers. + //! @ingroup derivGroup + //! + //! @warning This method is an experimental part of the %Cantera + //! API and may be changed or removed without notice. + virtual void buildWallJacobian(vector>& jacVector); + + //! Calculate flow contributions to the Jacobian of a Reactor specialization. + //! @param jac_vector vector where jacobian triplets are added + //! @warning Depending on the particular implementation, this may return an + //! approximate Jacobian intended only for use in forming a preconditioner for + //! iterative solvers. + //! @ingroup derivGroup + //! + //! @warning This method is an experimental part of the %Cantera + //! API and may be changed or removed without notice. + virtual void buildFlowJacobian(vector>& jacVector); //! Calculate the reactor-specific Jacobian using a finite difference method. //! @@ -320,6 +340,10 @@ class Reactor : public ReactorBase //! Vector of triplets representing the jacobian vector> m_jac_trips; + //! Boolean to skip walls in jacobian + bool m_jac_skip_walls = false; + //! Boolean to skip flow devices in jacobian + bool m_jac_skip_flow_devices = false; }; } diff --git a/include/cantera/zeroD/ReactorBase.h b/include/cantera/zeroD/ReactorBase.h index 0fed7946eb5..8f5a17ae870 100644 --- a/include/cantera/zeroD/ReactorBase.h +++ b/include/cantera/zeroD/ReactorBase.h @@ -302,12 +302,70 @@ class ReactorBase return m_sensParams.size(); } + /*! Calculate the derivative of temperature with respect to the temperature in the + * heat transfer equation based on the reactor specific equation of state. + * This function should also transform the state of the derivative to that + * appropriate for the jacobian's state/ + * @warning This function is an experimental part of the %Cantera API and may be changed + * or removed without notice. + * @since New in %Cantera 3.0. + */ + virtual double temperatureDerivative() { + throw NotImplementedError("Reactor::temperatureDerivative"); + } + + /*! Calculate the derivative of temperature with respect to the temperature in the + * heat transfer radiation equation based on the reactor specific equation of state. + * This function should also transform the state of the derivative to that + * appropriate for the jacobian's state/ + * @warning This function is an experimental part of the %Cantera API and may be changed + * or removed without notice. + * @since New in %Cantera 3.0. + */ + virtual double temperatureRadiationDerivative() { + throw NotImplementedError("Reactor::temperatureRadiationDerivative"); + } + + /*! Calculate the derivative of T with respect to the ith species in the heat + * transfer equation based on the reactor specific equation of state. + * @param index index of the species the derivative is with respect too + * @warning This function is an experimental part of the %Cantera API and may be changed + * or removed without notice. + * @since New in %Cantera 3.0. + */ + virtual double moleDerivative(size_t index) { + throw NotImplementedError("Reactor::moleDerivative"); + } + + /*! Calculate the derivative of T with respect to the ith species in the heat + * transfer radiation equation based on the reactor specific equation of state. + * @param index index of the species the derivative is with respect too + * @warning This function is an experimental part of the %Cantera API and may be changed + * or removed without notice. + * @since New in %Cantera 3.0. + */ + virtual double moleRadiationDerivative(size_t index) { + throw NotImplementedError("Reactor::moleRadiationDerivative"); + } + + //! Return the index associated with energy of the system + virtual size_t energyIndex() const { return m_eidx; }; + + //! Return the offset between species and state variables + virtual size_t speciesOffset() const { return m_sidx; }; + protected: explicit ReactorBase(const string& name="(none)"); //! Number of homogeneous species in the mixture size_t m_nsp = 0; + //! species offset in the state vector + const size_t m_sidx = 3; + + //! index of state variable associated with energy + const size_t m_eidx = 1; + ThermoPhase* m_thermo = nullptr; double m_vol = 0.0; //!< Current volume of the reactor [m^3] double m_mass = 0.0; //!< Current mass of the reactor [kg] diff --git a/include/cantera/zeroD/ReactorNet.h b/include/cantera/zeroD/ReactorNet.h index 1925a9eb1c3..24daa83bf75 100644 --- a/include/cantera/zeroD/ReactorNet.h +++ b/include/cantera/zeroD/ReactorNet.h @@ -281,6 +281,20 @@ class ReactorNet : public FuncEval //! reactor network. size_t globalComponentIndex(const string& component, size_t reactor=0); + //! Return the index corresponding to the start of the reactor specific state + //! vector in the reactor with index *reactor* in the global state vector for the + //! reactor network. + size_t globalStartIndex(Reactor* curr_reactor) { + for (size_t i = 0; i < m_reactors.size(); i++) { + if (curr_reactor == m_reactors[i]) { + return m_start[i]; + } + } + throw CanteraError("ReactorNet::globalStartIndex: ", + curr_reactor->name(), " not found in network."); + return npos; + } + //! Return the name of the i-th component of the global state vector. The //! name returned includes both the name of the reactor and the specific //! component, for example `'reactor1: CH4'`. @@ -481,6 +495,10 @@ class ReactorNet : public FuncEval //! "left hand side" of each governing equation vector m_LHS; vector m_RHS; + + //! derivative settings + bool m_jac_skip_walls = false; + bool m_jac_skip_flow_devices = false; }; diff --git a/include/cantera/zeroD/Wall.h b/include/cantera/zeroD/Wall.h index 8d2531badae..3d2c46d2301 100644 --- a/include/cantera/zeroD/Wall.h +++ b/include/cantera/zeroD/Wall.h @@ -9,6 +9,7 @@ #include "cantera/base/ctexceptions.h" #include "cantera/zeroD/ReactorBase.h" #include "ConnectorNode.h" +#include "cantera/numerics/eigen_sparse.h" namespace Cantera { @@ -85,6 +86,46 @@ class WallBase : public ConnectorNode m_time = time; } + /*! Build the Jacobian terms specific to the flow device for the given connected reactor. + * @param r a pointer to the calling reactor + * @param jacVector a vector of triplets to be added to the jacobian for the reactor + * @warning This function is an experimental part of the %Cantera API and may be + * changed + * or removed without notice. + * @since New in %Cantera 3.0. + */ + virtual void buildReactorJacobian(ReactorBase* r, vector>& jacVector) { + throw NotImplementedError("WallBase::buildReactorJacobian"); + } + + /*! Build the Jacobian terms specific to the flow device for the network. These terms + * will be adjusted to the networks indexing system outside of the reactor. + * @param jacVector a vector of triplets to be added to the jacobian for the reactor + * @warning This function is an experimental part of the %Cantera API and may be + * changed + * or removed without notice. + * @since New in %Cantera 3.0. + */ + virtual void buildNetworkJacobian(vector>& jacVector) { + throw NotImplementedError("WallBase::buildNetworkJacobian"); + } + + /*! Specify the jacobian terms have been calculated and should not be recalculated. + * @warning This function is an experimental part of the %Cantera API and may be + * changed + * or removed without notice. + * @since New in %Cantera 3.0. + */ + void jacobianCalculated() { m_jac_calculated = true; }; + + /*! Specify that jacobian terms have not been calculated and should be recalculated. + * @warning This function is an experimental part of the %Cantera API and may be + * changed + * or removed without notice. + * @since New in %Cantera 3.0. + */ + void jacobianNotCalculated() { m_jac_calculated = false; }; + protected: ReactorBase* m_left = nullptr; ReactorBase* m_right = nullptr; @@ -93,6 +134,10 @@ class WallBase : public ConnectorNode double m_time = 0.0; double m_area = 1.0; + + //! a variable to switch on and off so calculations are not doubled by the calling + //! reactor or network + bool m_jac_calculated = false; }; //! Represents a wall between between two ReactorBase objects. @@ -200,6 +245,10 @@ class Wall : public WallBase return m_k; } + virtual void buildReactorJacobian(ReactorBase* r, vector>& jacVector) override; + + virtual void buildNetworkJacobian(vector>& jacVector) override; + protected: //! expansion rate coefficient diff --git a/interfaces/cython/cantera/reactor.pxd b/interfaces/cython/cantera/reactor.pxd index 6c19b9f4603..5a1ba565fcc 100644 --- a/interfaces/cython/cantera/reactor.pxd +++ b/interfaces/cython/cantera/reactor.pxd @@ -207,6 +207,8 @@ cdef extern from "cantera/zerodim.h" namespace "Cantera": void setPreconditioner(shared_ptr[CxxSystemJacobian] preconditioner) void setDerivativeSettings(CxxAnyMap&) CxxAnyMap solverStats() except +translate_exception + CxxSparseMatrix jacobian() except +translate_exception + CxxSparseMatrix finiteDifferenceJacobian() except +translate_exception cdef extern from "cantera/zeroD/ReactorDelegator.h" namespace "Cantera": cdef cppclass CxxReactorAccessor "Cantera::ReactorAccessor": diff --git a/interfaces/cython/cantera/reactor.pyx b/interfaces/cython/cantera/reactor.pyx index 66e83868b28..80c707ffda0 100644 --- a/interfaces/cython/cantera/reactor.pyx +++ b/interfaces/cython/cantera/reactor.pyx @@ -2082,6 +2082,31 @@ cdef class ReactorNet: def __set__(self, settings): self.net.setDerivativeSettings(py_to_anymap(settings)) + property jacobian: + """ + Get the system Jacobian or an approximation thereof. + + **Warning**: Depending on the particular implementation, this may return an + approximate Jacobian intended only for use in forming a preconditioner for + iterative solvers, excluding terms that would generate a fully-dense Jacobian. + + **Warning**: This method is an experimental part of the Cantera API and may be + changed or removed without notice. + """ + def __get__(self): + return get_from_sparse(self.net.jacobian(), self.n_vars, self.n_vars) + + property finite_difference_jacobian: + """ + Get the system Jacobian, calculated using a finite difference method. + + **Warning:** this property is an experimental part of the Cantera API and + may be changed or removed without notice. + """ + def __get__(self): + return get_from_sparse(self.net.finiteDifferenceJacobian(), + self.n_vars, self.n_vars) + def draw(self, *, graph_attr=None, node_attr=None, edge_attr=None, heat_flow_attr=None, mass_flow_attr=None, moving_wall_edge_attr=None, surface_edge_attr=None, show_wall_velocity=True, print_state=False, diff --git a/src/zeroD/IdealGasConstPressureMoleReactor.cpp b/src/zeroD/IdealGasConstPressureMoleReactor.cpp index 42538b3df67..0b5b24f9684 100644 --- a/src/zeroD/IdealGasConstPressureMoleReactor.cpp +++ b/src/zeroD/IdealGasConstPressureMoleReactor.cpp @@ -200,6 +200,7 @@ void IdealGasConstPressureMoleReactor::buildJacobian( jacVector.emplace_back(static_cast(j), 0, (ydotPerturbed - ydotCurrent) / deltaTemp); } + // d T_dot/dnj // allocate vectors for whole system Eigen::VectorXd enthalpy = Eigen::VectorXd::Zero(ssize); @@ -225,7 +226,13 @@ void IdealGasConstPressureMoleReactor::buildJacobian( jacVector.emplace_back(0, static_cast(j + m_sidx), (specificHeat[j] * qdot - NCp * hk_dnkdnj_sums[j]) * denom); } + + // build wall jacobian + buildWallJacobian(jacVector); } + + // build flow jacobian + buildFlowJacobian(jacVector); } size_t IdealGasConstPressureMoleReactor::componentIndex(const string& nm) const @@ -282,4 +289,23 @@ double IdealGasConstPressureMoleReactor::lowerBound(size_t k) const { } } +double IdealGasConstPressureMoleReactor::moleDerivative(size_t index) +{ + // derivative of temperature transformed by ideal gas law + vector moles(m_nsp); + getMoles(moles.data()); + double dTdni = pressure() * m_vol / GasConstant / std::accumulate(moles.begin(), moles.end(), 0.0); + return dTdni; +} + +double IdealGasConstPressureMoleReactor::moleRadiationDerivative(size_t index) +{ + // derivative of temperature transformed by ideal gas law + vector moles(m_nsp); + getMoles(moles.data()); + double dT4dni = std::pow(pressure() * m_vol / GasConstant, 4); + dT4dni *= std::pow(1 / std::accumulate(moles.begin(), moles.end(), 0.0), 5); + return dT4dni; +} + } diff --git a/src/zeroD/IdealGasMoleReactor.cpp b/src/zeroD/IdealGasMoleReactor.cpp index 9292e8d4afe..ec3e25d39e3 100644 --- a/src/zeroD/IdealGasMoleReactor.cpp +++ b/src/zeroD/IdealGasMoleReactor.cpp @@ -244,6 +244,7 @@ void IdealGasMoleReactor::buildJacobian(vector>& jacVecto jacVector.emplace_back(static_cast(j), 0, (ydotPerturbed - ydotCurrent) / deltaTemp); } + // d T_dot/dnj Eigen::VectorXd netProductionRates = Eigen::VectorXd::Zero(ssize); Eigen::VectorXd internal_energy = Eigen::VectorXd::Zero(ssize); @@ -273,7 +274,32 @@ void IdealGasMoleReactor::buildJacobian(vector>& jacVecto jacVector.emplace_back(0, static_cast(j + m_sidx), (specificHeat[j] * qdot - NCv * uk_dnkdnj_sums[j]) * denom); } + + // build wall jacobian + buildWallJacobian(jacVector); } + + // build flow jacobian + buildFlowJacobian(jacVector); +} + +double IdealGasMoleReactor::moleDerivative(size_t index) +{ + // derivative of temperature transformed by ideal gas law + vector moles(m_nsp); + getMoles(moles.data()); + double dTdni = pressure() * m_vol / GasConstant / std::accumulate(moles.begin(), moles.end(), 0.0); + return dTdni; +} + +double IdealGasMoleReactor::moleRadiationDerivative(size_t index) +{ + // derivative of temperature transformed by ideal gas law + vector moles(m_nsp); + getMoles(moles.data()); + double dT4dni = std::pow(pressure() * m_vol / GasConstant, 4); + dT4dni *= std::pow(1 / std::accumulate(moles.begin(), moles.end(), 0.0), 5); + return dT4dni; } } diff --git a/src/zeroD/Reactor.cpp b/src/zeroD/Reactor.cpp index 58b7ee125c4..2b6e78d0f5c 100644 --- a/src/zeroD/Reactor.cpp +++ b/src/zeroD/Reactor.cpp @@ -42,6 +42,17 @@ void Reactor::setDerivativeSettings(AnyMap& settings) for (auto S : m_surfaces) { S->kinetics()->setDerivativeSettings(settings); } + + // set reactor settings + bool force = settings.empty(); + if (force || settings.hasKey("skip-walls")) { + m_jac_skip_walls = settings.getBool("skip-walls", + false); + } + if (force || settings.hasKey("skip-flow-devices")) { + m_jac_skip_flow_devices = settings.getBool("skip-flow-devices", + false); + } } void Reactor::getState(double* y) @@ -631,4 +642,26 @@ void Reactor::setAdvanceLimit(const string& nm, const double limit) } } +void Reactor:: buildWallJacobian(vector>& jacVector) +{ + if (!m_jac_skip_walls) { + for (size_t i = 0; i < m_wall.size(); i++) { + m_wall[i]->buildReactorJacobian(this, jacVector); + } + } +} + +void Reactor::buildFlowJacobian(vector>& jacVector) +{ + if (!m_jac_skip_flow_devices) { + for (size_t i = 0; i < m_outlet.size(); i++) { + m_outlet[i]->buildReactorJacobian(this, jacVector); + } + + for (size_t i = 0; i buildReactorJacobian(this, jacVector); + } + } +} + } diff --git a/src/zeroD/ReactorNet.cpp b/src/zeroD/ReactorNet.cpp index 90f0699f11b..e6768f2d2b6 100644 --- a/src/zeroD/ReactorNet.cpp +++ b/src/zeroD/ReactorNet.cpp @@ -726,6 +726,16 @@ void ReactorNet::setDerivativeSettings(AnyMap& settings) for (size_t i = 0; i < m_reactors.size(); i++) { m_reactors[i]->setDerivativeSettings(settings); } + // set network settings + bool force = settings.empty(); + if (force || settings.hasKey("skip-walls")) { + m_jac_skip_walls = settings.getBool("skip-walls", + false); + } + if (force || settings.hasKey("skip-flow-devices")) { + m_jac_skip_flow_devices = settings.getBool("skip-flow-devices", + false); + } } AnyMap ReactorNet::solverStats() const @@ -813,6 +823,26 @@ void ReactorNet::buildJacobian(vector>& jacVector) if (! m_init) { initialize(); } + // loop through and set connectors not found + for (auto r : m_reactors) { + // walls + if (!m_jac_skip_walls) { + for (size_t i = 0; i < r->nWalls(); i++) { + r->wall(i).jacobianNotCalculated(); + } + } + // flow devices + if (!m_jac_skip_flow_devices) { + // outlets + for (size_t i = 0; i < r->nOutlets(); i++) { + r->outlet(i).jacobianNotCalculated(); + } + // inlets + for (size_t i = 0; i < r->nInlets(); i++) { + r->inlet(i).jacobianNotCalculated(); + } + } + } // Create jacobian triplet vector vector jstarts; // Get jacobians and give elements to preconditioners @@ -830,6 +860,31 @@ void ReactorNet::buildJacobian(vector>& jacVector) jacVector[j] = newTrip; } } + + // loop through all connections and then set them found so calculations are not + // repeated + for (auto r : m_reactors) { + // walls + if (!m_jac_skip_walls) { + for (size_t i = 0; i < r->nWalls(); i++) { + r->wall(i).buildNetworkJacobian(jacVector); + r->wall(i).jacobianCalculated(); + } + } + // flow devices + if (m_jac_skip_flow_devices) { + // outlets + for (size_t i = 0; i < r->nOutlets(); i++) { + r->outlet(i).buildNetworkJacobian(jacVector); + r->outlet(i).jacobianCalculated(); + } + // inlets + for (size_t i = 0; i < r->nInlets(); i++) { + r->inlet(i).buildNetworkJacobian(jacVector); + r->inlet(i).jacobianCalculated(); + } + } + } } Eigen::SparseMatrix ReactorNet::finiteDifferenceJacobian() diff --git a/src/zeroD/Wall.cpp b/src/zeroD/Wall.cpp index 19c2006f6d6..b02e8b10c35 100644 --- a/src/zeroD/Wall.cpp +++ b/src/zeroD/Wall.cpp @@ -6,6 +6,8 @@ #include "cantera/base/stringUtils.h" #include "cantera/numerics/Func1.h" #include "cantera/zeroD/Wall.h" +#include "cantera/thermo/ThermoPhase.h" +#include "cantera/zeroD/ReactorNet.h" namespace Cantera { @@ -75,4 +77,66 @@ double Wall::heatRate() return q1; } + +void Wall::buildReactorJacobian(ReactorBase* r, vector>& jacVector) +{ + // get derivative of heat transfer for both reactors + vector> network; + size_t nsp = r->phase()->thermo()->nSpecies(); + size_t sidx = r->speciesOffset(); + size_t eidx = r->energyIndex(); + // define a scalar for direction based on left and right + double scalar = (r == m_left) ? 1.0 : -1.0; + // elements within the current reactor + // find dQdni for the current reactor w.r.t current reactor + for (size_t i = sidx; i < nsp + sidx; i++) { + double dQdni = m_rrth * m_area * scalar * r->moleDerivative(i); + dQdni += m_emiss * m_area * scalar * r->moleRadiationDerivative(i); + jacVector.emplace_back(eidx, i, dQdni); + } +} + +void Wall::buildNetworkJacobian(vector>& jacVector) +{ + // No interdependent terms for reservoirs + if (m_right->type() == "Reservoir" || m_left->type() == "Reservoir") { + return; + } + // return if the jacobian has been calculated + if (m_jac_calculated) { + return; + } + // get derivatives for inter-dependent reactor terms + //variables for the right side + vector> network; + size_t r_nsp = m_right->phase()->thermo()->nSpecies(); + size_t r_sidx = m_right->speciesOffset(); + size_t r_net = m_right->network().globalStartIndex((Reactor* ) m_right); + size_t r_eidx = m_right->energyIndex(); + + // variables for the left side + size_t l_nsp = m_left->phase()->thermo()->nSpecies(); + size_t l_sidx = m_left->speciesOffset(); + size_t l_net = m_left->network().globalStartIndex((Reactor* ) m_left); + size_t l_eidx = m_left->energyIndex(); + + if (((Reactor* ) m_right)->energyEnabled()) { + // find dQdni for the right reactor w.r.t left reactor + for (size_t i = l_sidx; i < l_sidx + l_nsp; i++) { + double dQdni = m_rrth * m_area * m_left->moleDerivative(i); + dQdni += m_emiss * m_area * m_left->moleRadiationDerivative(i); + jacVector.emplace_back(r_eidx + r_net, i + l_net, dQdni); + } + } + + if (((Reactor* ) m_left)->energyEnabled()) { + // find dQdni for the left reactor w.r.t right reactor + for (size_t i = r_sidx; i < r_sidx + r_nsp; i++) { + double dQdni = - m_rrth * m_area * m_right->moleDerivative(i); + dQdni -= m_emiss * m_area * m_right->moleRadiationDerivative(i); + jacVector.emplace_back(l_eidx + l_net, i + r_net, dQdni); + } + } +} + } diff --git a/test/python/test_reactor.py b/test/python/test_reactor.py index 8204fb56bf2..2ac1f4bd17e 100644 --- a/test/python/test_reactor.py +++ b/test/python/test_reactor.py @@ -1538,21 +1538,67 @@ def test_get_solver_type(self): self.net2.initialize() assert self.net2.linear_solver_type == "GMRES" + def test_heat_transfer_network(self): + # create first reactor + gas1 = ct.Solution("h2o2.yaml", "ohmech") + gas1.TPX = 600, ct.one_atm, "O2:1.0" + r1 = self.reactorClass(gas1) + r1.chemistry_enabled = False + + # create second reactor + gas2 = ct.Solution("h2o2.yaml", "ohmech") + gas2.TPX = 300, ct.one_atm, "O2:1.0" + r2 = self.reactorClass(gas2) + r2.chemistry_enabled = False + + # create wall + U = 2.0 + A = 3.0 + w = ct.Wall(r1, r2, U=U, A=A) + net = ct.ReactorNet([r1, r2]) + jac = net.jacobian + + # check for values + for i in range(gas1.n_species): + assert np.isclose(jac[0, i + 1 + r1.n_vars], - U * A * gas2.T) + assert np.isclose(jac[r1.n_vars, i + 1], U * A * gas1.T) + if (i + 1 != 2): + assert np.isclose(jac[0, i + 1], U * A * gas1.T, rtol=1e-2) + assert np.isclose(jac[r1.n_vars, i + 1 + r1.n_vars], - U * A * gas2.T, rtol=1e-2) + class TestIdealGasMoleReactor(TestMoleReactor): reactorClass = ct.IdealGasMoleReactor test_preconditioner_unsupported = None - @pytest.mark.diagnose def test_heat_transfer_network(self): - # Result should be the same if (m * cp) / (U * A) is held constant - self.make_reactors(T1=300, T2=1000) - self.add_wall(U=200, A=1.0) - self.net.preconditioner = ct.AdaptivePreconditioner() - print(self.net.finite_difference_jacobian) - # self.net.advance(1.0) - # T1a = self.r1.T - # T2a = self.r2.T + # create first reactor + gas1 = ct.Solution("h2o2.yaml", "ohmech") + gas1.TPX = 600, ct.one_atm, "O2:1.0" + r1 = self.reactorClass(gas1) + r1.chemistry_enabled = False + + # create second reactor + gas2 = ct.Solution("h2o2.yaml", "ohmech") + gas2.TPX = 300, ct.one_atm, "O2:1.0" + r2 = self.reactorClass(gas2) + r2.chemistry_enabled = False + + # create wall + U = 2.0 + A = 3.0 + w = ct.Wall(r1, r2, U=U, A=A) + net = ct.ReactorNet([r1, r2]) + jac = net.jacobian + + # check for values + for i in range(gas1.n_species): + assert np.isclose(jac[0, i + 2 + r1.n_vars], - U * A * gas2.T) + assert np.isclose(jac[r1.n_vars, i + 2], U * A * gas1.T) + if (i + 2 != 3): + assert np.isclose(jac[0, i + 2], U * A * gas1.T, rtol=1e-2) + assert np.isclose(jac[r1.n_vars, i + 2 + r1.n_vars], - U * A * gas2.T, rtol=1e-2) + def test_adaptive_precon_integration(self): # Network one with non-mole reactor From 17323ab9231123717aa6ccdbaf50c20a55344102 Mon Sep 17 00:00:00 2001 From: Anthony Walker Date: Tue, 17 Oct 2023 23:03:54 -0700 Subject: [PATCH 04/13] Update Extensible Jacobian interface to not reconstruct entire triplet vector --- interfaces/cython/cantera/_utils.pxd | 2 ++ interfaces/cython/cantera/_utils.pyx | 4 ++++ interfaces/cython/cantera/delegator.pyx | 10 ++++++---- 3 files changed, 12 insertions(+), 4 deletions(-) diff --git a/interfaces/cython/cantera/_utils.pxd b/interfaces/cython/cantera/_utils.pxd index 77f98c958b9..25707835407 100644 --- a/interfaces/cython/cantera/_utils.pxd +++ b/interfaces/cython/cantera/_utils.pxd @@ -122,3 +122,5 @@ cdef anyvalue_to_python(string name, CxxAnyValue& v) cdef vector[CxxEigenTriplet] python_to_triplets(triplets) except * cdef triplets_to_python(vector[CxxEigenTriplet]& triplets) + +cdef CxxEigenTriplet get_triplet(row, col, val) except * diff --git a/interfaces/cython/cantera/_utils.pyx b/interfaces/cython/cantera/_utils.pyx index 956558788a7..62216e67e88 100644 --- a/interfaces/cython/cantera/_utils.pyx +++ b/interfaces/cython/cantera/_utils.pyx @@ -547,6 +547,10 @@ cdef vector[CxxEigenTriplet] python_to_triplets(triplets): trips.push_back(dereference(trip_ptr)) return trips +cdef CxxEigenTriplet get_triplet(row, col, val): + cdef CxxEigenTriplet* trip_ptr = new CxxEigenTriplet(row, col, val) + return dereference(trip_ptr) + cdef triplets_to_python(vector[CxxEigenTriplet]& triplets): values = [] cdef size_t row diff --git a/interfaces/cython/cantera/delegator.pyx b/interfaces/cython/cantera/delegator.pyx index f209b18e275..0dab6e43621 100644 --- a/interfaces/cython/cantera/delegator.pyx +++ b/interfaces/cython/cantera/delegator.pyx @@ -9,7 +9,7 @@ from libc.stdlib cimport malloc from libc.string cimport strcpy from ._utils import CanteraError -from ._utils cimport stringify, pystr, anymap_to_py, py_to_anymap, triplets_to_python, python_to_triplets +from ._utils cimport stringify, pystr, anymap_to_py, py_to_anymap, triplets_to_python, python_to_triplets, get_triplet from .units cimport Units, UnitStack # from .reaction import ExtensibleRate, ExtensibleRateData from .reaction cimport (ExtensibleRate, ExtensibleRateData, CxxReaction, @@ -247,12 +247,14 @@ cdef void callback_v_d_dp_dp(PyFuncInfo& funcInfo, size_array2 sizes, double arg # Wrapper for void(vector&) cdef void callback_v_vet(PyFuncInfo& funcInfo, vector[CxxEigenTriplet]& jac_vector) noexcept: try: - python_trips = triplets_to_python(jac_vector) + python_trips = [] # convert vector to python object (funcInfo.func())(python_trips) - jac_vector = python_to_triplets(python_trips) + # add the triplets to the jacobian vector + for r, c, v in python_trips: + jac_vector.push_back(get_triplet(r, c, v)) except BaseException as e: - exc_type, exc_value = sys.exc_info()[:2] + exc_type, exc_value = _sys.exc_info()[:2] funcInfo.setExceptionType(exc_type) funcInfo.setExceptionValue(exc_value) From 9b5a46c03ecedbab9ff3fc9ae84bd19999963d47 Mon Sep 17 00:00:00 2001 From: Anthony Walker Date: Tue, 24 Oct 2023 11:52:47 -0400 Subject: [PATCH 05/13] Cpp Wall Jacobian Test --- test/zeroD/test_zeroD.cpp | 81 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 81 insertions(+) diff --git a/test/zeroD/test_zeroD.cpp b/test/zeroD/test_zeroD.cpp index 2dc6f020b24..343aaa3a86d 100644 --- a/test/zeroD/test_zeroD.cpp +++ b/test/zeroD/test_zeroD.cpp @@ -300,6 +300,87 @@ TEST(AdaptivePreconditionerTests, test_precon_solver_stats) EXPECT_GE(stats["nonlinear_conv_fails"].asInt(), 0); } +TEST(JacobianTests, test_wall_jacobian_build) +{ + // create first reactor + auto sol1 = newSolution("h2o2.yaml"); + sol1->thermo()->setState_TPY(1000.0, OneAtm, " O2:1.0"); + IdealGasMoleReactor reactor1; + reactor1.insert(sol1); + reactor1.setInitialVolume(1.0); + // create second reactor + auto sol2 = newSolution("h2o2.yaml"); + sol2->thermo()->setState_TPY(900.0, OneAtm, " O2:1.0"); + IdealGasConstPressureMoleReactor reactor2; + reactor2.insert(sol2); + reactor2.setInitialVolume(1.0); + // create the wall + Wall w; + w.install(reactor1, reactor2); + w.setArea(2.0); + w.setHeatTransferCoeff(3.0); + // setup reactor network and integrate + ReactorNet network; + network.addReactor(reactor1); + network.addReactor(reactor2); + network.initialize(); + // create jacobian the size of network + Eigen::SparseMatrix wallJacMat; + wallJacMat.resize(network.neq(), network.neq()); + // manually get wall jacobian elements + vector> wallJac; + // build jac for reactor 1 wall only + w.buildReactorJacobian(&reactor1, wallJac); + wallJacMat.setFromTriplets(wallJac.begin(), wallJac.end()); + // check appropriate values + // double tol = 1e-8; + for (int k = 0; k < wallJacMat.outerSize(); k++) { + for (Eigen::SparseMatrix::InnerIterator it(wallJacMat, k); it; ++it) { + EXPECT_DOUBLE_EQ(it.value(), 6000.0); + EXPECT_EQ(it.row(), 0); // check that it is the first row + EXPECT_GE(it.col(), reactor1.speciesOffset()); + EXPECT_LT(it.col(), reactor1.neq()); + } + } + // build jac for reactor 2 wall only + wallJac.clear(); + w.buildReactorJacobian(&reactor2, wallJac); + wallJacMat.setZero(); + wallJacMat.setFromTriplets(wallJac.begin(), wallJac.end()); + // check appropriate values + // double tol = 1e-8; + for (int k = 0; k < wallJacMat.outerSize(); k++) { + for (Eigen::SparseMatrix::InnerIterator it(wallJacMat, k); it; ++it) { + EXPECT_DOUBLE_EQ(it.value(), -5400.0); + EXPECT_EQ(it.row(), 0); // check that it is the first row + EXPECT_GE(it.col(), reactor2.speciesOffset()); + EXPECT_LT(it.col(), reactor2.neq()); + } + } + // build jac for network terms + wallJac.clear(); + w.buildNetworkJacobian(wallJac); + wallJacMat.setZero(); + wallJacMat.setFromTriplets(wallJac.begin(), wallJac.end()); + // check appropriate values + // double tol = 1e-8; + for (int k = 0; k < wallJacMat.outerSize(); k++) { + for (Eigen::SparseMatrix::InnerIterator it(wallJacMat, k); it; ++it) { + if (it.value() < 0) { + EXPECT_DOUBLE_EQ(it.value(), -5400.0); + EXPECT_EQ(it.row(), 0); // check that it is the first row + EXPECT_GE(it.col(), reactor1.neq() + reactor2.speciesOffset()); + EXPECT_LT(it.col(), reactor1.neq() + reactor2.neq()); + } else { + EXPECT_DOUBLE_EQ(it.value(), 6000.0); + EXPECT_EQ(it.row(), reactor1.neq()); // check that it is the first row + EXPECT_GE(it.col(), reactor2.speciesOffset()); + EXPECT_LT(it.col(), reactor1.neq()); + } + } + } +} + int main(int argc, char** argv) { printf("Running main() from test_zeroD.cpp\n"); From 87d6ce180e9c900464beee4d39e201cbe7d3e950 Mon Sep 17 00:00:00 2001 From: Anthony Walker Date: Tue, 24 Oct 2023 12:11:50 -0400 Subject: [PATCH 06/13] Doxygen fixes --- include/cantera/zeroD/FlowDevice.h | 55 ++++++++++++++----------- include/cantera/zeroD/Reactor.h | 5 +-- include/cantera/zeroD/ReactorBase.h | 62 +++++++++++++++-------------- include/cantera/zeroD/ReactorNet.h | 3 +- include/cantera/zeroD/Wall.h | 60 +++++++++++++++------------- src/zeroD/ReactorNet.cpp | 6 +-- 6 files changed, 101 insertions(+), 90 deletions(-) diff --git a/include/cantera/zeroD/FlowDevice.h b/include/cantera/zeroD/FlowDevice.h index b9b5b34ad5d..8662e40cf04 100644 --- a/include/cantera/zeroD/FlowDevice.h +++ b/include/cantera/zeroD/FlowDevice.h @@ -141,40 +141,47 @@ class FlowDevice : public ConnectorNode m_time = time; } - /*! Build the Jacobian terms specific to the flow device for the given connected reactor. - * @param r a pointer to the calling reactor - * @param jacVector a vector of triplets to be added to the jacobian for the reactor - * @warning This function is an experimental part of the %Cantera API and may be changed - * or removed without notice. - * @since New in %Cantera 3.0. - */ + //! Build the Jacobian terms specific to the flow device for the given connected + //! reactor. + //! @param r a pointer to the calling reactor + //! @param jacVector a vector of triplets to be added to the jacobian for the + //! reactor + //! @warning This function is an experimental part of the %Cantera API and may be + //! changed + //! or removed without notice. + //! @since New in %Cantera 3.0. + //! virtual void buildReactorJacobian(ReactorBase* r, vector>& jacVector) { throw NotImplementedError(type() + "::buildReactorJacobian"); } - /*! Build the Jacobian terms specific to the flow device for the network. These terms - * will be adjusted to the networks indexing system outside of the reactor. - * @param jacVector a vector of triplets to be added to the jacobian for the reactor - * @warning This function is an experimental part of the %Cantera API and may be changed - * or removed without notice. - * @since New in %Cantera 3.0. - */ + //! Build the Jacobian terms specific to the flow device for the network. These + //! terms + //! will be adjusted to the networks indexing system outside of the reactor. + //! @param jacVector a vector of triplets to be added to the jacobian for the + //! reactor + //! @warning This function is an experimental part of the %Cantera API and may be + //! changed + //! or removed without notice. + //! @since New in %Cantera 3.0. + //! virtual void buildNetworkJacobian(vector>& jacVector) { throw NotImplementedError(type() + "::buildNetworkJacobian"); } - /*! Specify the jacobian terms have been calculated and should not be recalculated. - * @warning This function is an experimental part of the %Cantera API and may be changed - * or removed without notice. - * @since New in %Cantera 3.0. - */ + //! Specify the jacobian terms have been calculated and should not be recalculated. + //! @warning This function is an experimental part of the %Cantera API and may be + //! changed + //! or removed without notice. + //! @since New in %Cantera 3.0. + //! void jacobianCalculated() { m_jac_calculated = true; }; - /*! Specify that jacobian terms have not been calculated and should be recalculated. - * @warning This function is an experimental part of the %Cantera API and may be changed - * or removed without notice. - * @since New in %Cantera 3.0. - */ + //! Specify that jacobian terms have not been calculated and should be recalculated. + //! @warning This function is an experimental part of the %Cantera API and may be changed + //! or removed without notice. + //! @since New in %Cantera 3.0. + //! void jacobianNotCalculated() { m_jac_calculated = false; }; protected: diff --git a/include/cantera/zeroD/Reactor.h b/include/cantera/zeroD/Reactor.h index 20bc24c1232..0ffecc16d67 100644 --- a/include/cantera/zeroD/Reactor.h +++ b/include/cantera/zeroD/Reactor.h @@ -218,8 +218,7 @@ class Reactor : public ReactorBase } //! Calculate the Jacobian of a specific Reactor specialization. - //! @param jac_vector vector where jacobian triplets are added - //! @param offset offset added to the row and col indices of the elements + //! @param jacVector vector where jacobian triplets are added //! @warning Depending on the particular implementation, this may return an //! approximate Jacobian intended only for use in forming a preconditioner for //! iterative solvers. @@ -243,7 +242,7 @@ class Reactor : public ReactorBase virtual void buildWallJacobian(vector>& jacVector); //! Calculate flow contributions to the Jacobian of a Reactor specialization. - //! @param jac_vector vector where jacobian triplets are added + //! @param jacVector vector where jacobian triplets are added //! @warning Depending on the particular implementation, this may return an //! approximate Jacobian intended only for use in forming a preconditioner for //! iterative solvers. diff --git a/include/cantera/zeroD/ReactorBase.h b/include/cantera/zeroD/ReactorBase.h index 8f5a17ae870..9d70d634acf 100644 --- a/include/cantera/zeroD/ReactorBase.h +++ b/include/cantera/zeroD/ReactorBase.h @@ -302,48 +302,50 @@ class ReactorBase return m_sensParams.size(); } - /*! Calculate the derivative of temperature with respect to the temperature in the - * heat transfer equation based on the reactor specific equation of state. - * This function should also transform the state of the derivative to that - * appropriate for the jacobian's state/ - * @warning This function is an experimental part of the %Cantera API and may be changed - * or removed without notice. - * @since New in %Cantera 3.0. - */ + //! Calculate the derivative of temperature with respect to the temperature in the + //! heat transfer equation based on the reactor specific equation of state. + //! This function should also transform the state of the derivative to that + //! appropriate for the jacobian's state/ + //! @warning This function is an experimental part of the %Cantera API and may be changed + //! or removed without notice. + //! @since New in %Cantera 3.0. + //! virtual double temperatureDerivative() { throw NotImplementedError("Reactor::temperatureDerivative"); } - /*! Calculate the derivative of temperature with respect to the temperature in the - * heat transfer radiation equation based on the reactor specific equation of state. - * This function should also transform the state of the derivative to that - * appropriate for the jacobian's state/ - * @warning This function is an experimental part of the %Cantera API and may be changed - * or removed without notice. - * @since New in %Cantera 3.0. - */ + //! Calculate the derivative of temperature with respect to the temperature in the + //! heat transfer radiation equation based on the reactor specific equation of + //! state. + //! This function should also transform the state of the derivative to that + //! appropriate for the jacobian's state/ + //! @warning This function is an experimental part of the %Cantera API and may be + //! changed + //! or removed without notice. + //! @since New in %Cantera 3.0. + //! virtual double temperatureRadiationDerivative() { throw NotImplementedError("Reactor::temperatureRadiationDerivative"); } - /*! Calculate the derivative of T with respect to the ith species in the heat - * transfer equation based on the reactor specific equation of state. - * @param index index of the species the derivative is with respect too - * @warning This function is an experimental part of the %Cantera API and may be changed - * or removed without notice. - * @since New in %Cantera 3.0. - */ + //! Calculate the derivative of T with respect to the ith species in the heat + //! transfer equation based on the reactor specific equation of state. + //! @param index index of the species the derivative is with respect too + //! @warning This function is an experimental part of the %Cantera API and may be changed + //! or removed without notice. + //! @since New in %Cantera 3.0. + //! virtual double moleDerivative(size_t index) { throw NotImplementedError("Reactor::moleDerivative"); } - /*! Calculate the derivative of T with respect to the ith species in the heat - * transfer radiation equation based on the reactor specific equation of state. - * @param index index of the species the derivative is with respect too - * @warning This function is an experimental part of the %Cantera API and may be changed - * or removed without notice. - * @since New in %Cantera 3.0. - */ + //! Calculate the derivative of T with respect to the ith species in the heat + //! transfer radiation equation based on the reactor specific equation of state. + //! @param index index of the species the derivative is with respect too + //! @warning This function is an experimental part of the %Cantera API and may be changed + //! or removed without notice. + //! @since New in %Cantera 3.0. + //! virtual double moleRadiationDerivative(size_t index) { throw NotImplementedError("Reactor::moleRadiationDerivative"); } diff --git a/include/cantera/zeroD/ReactorNet.h b/include/cantera/zeroD/ReactorNet.h index 24daa83bf75..55901d2e6eb 100644 --- a/include/cantera/zeroD/ReactorNet.h +++ b/include/cantera/zeroD/ReactorNet.h @@ -414,8 +414,7 @@ class ReactorNet : public FuncEval void addReactor(shared_ptr reactor); //! Calculate the Jacobian of the entire reactor network. - //! @param jac_vector vector where jacobian triplets are added - //! @param offset offset added to the row and col indices of the elements + //! @param jacVector vector where jacobian triplets are added //! @warning Depending on the particular implementation, this may return an //! approximate Jacobian intended only for use in forming a preconditioner for //! iterative solvers. diff --git a/include/cantera/zeroD/Wall.h b/include/cantera/zeroD/Wall.h index 3d2c46d2301..9fb79322452 100644 --- a/include/cantera/zeroD/Wall.h +++ b/include/cantera/zeroD/Wall.h @@ -86,44 +86,48 @@ class WallBase : public ConnectorNode m_time = time; } - /*! Build the Jacobian terms specific to the flow device for the given connected reactor. - * @param r a pointer to the calling reactor - * @param jacVector a vector of triplets to be added to the jacobian for the reactor - * @warning This function is an experimental part of the %Cantera API and may be - * changed - * or removed without notice. - * @since New in %Cantera 3.0. - */ + //! Build the Jacobian terms specific to the flow device for the given connected + //! reactor. + //! @param r a pointer to the calling reactor + //! @param jacVector a vector of triplets to be added to the jacobian for the + //! reactor + //! @warning This function is an experimental part of the %Cantera API and may be + //! changed + //! or removed without notice. + //! @since New in %Cantera 3.0. + //! virtual void buildReactorJacobian(ReactorBase* r, vector>& jacVector) { throw NotImplementedError("WallBase::buildReactorJacobian"); } - /*! Build the Jacobian terms specific to the flow device for the network. These terms - * will be adjusted to the networks indexing system outside of the reactor. - * @param jacVector a vector of triplets to be added to the jacobian for the reactor - * @warning This function is an experimental part of the %Cantera API and may be - * changed - * or removed without notice. - * @since New in %Cantera 3.0. - */ + //! Build the Jacobian terms specific to the flow device for the network. These + //! terms + //! will be adjusted to the networks indexing system outside of the reactor. + //! @param jacVector a vector of triplets to be added to the jacobian for the + //! reactor + //! @warning This function is an experimental part of the %Cantera API and may be + //! changed + //! or removed without notice. + //! @since New in %Cantera 3.0. + //! virtual void buildNetworkJacobian(vector>& jacVector) { throw NotImplementedError("WallBase::buildNetworkJacobian"); } - /*! Specify the jacobian terms have been calculated and should not be recalculated. - * @warning This function is an experimental part of the %Cantera API and may be - * changed - * or removed without notice. - * @since New in %Cantera 3.0. - */ + //! Specify the jacobian terms have been calculated and should not be recalculated. + //! @warning This function is an experimental part of the %Cantera API and may be + //! changed + //! or removed without notice. + //! @since New in %Cantera 3.0. + //! void jacobianCalculated() { m_jac_calculated = true; }; - /*! Specify that jacobian terms have not been calculated and should be recalculated. - * @warning This function is an experimental part of the %Cantera API and may be - * changed - * or removed without notice. - * @since New in %Cantera 3.0. - */ + //! Specify that jacobian terms have not been calculated and should be recalculated. + //! @warning This function is an experimental part of the %Cantera API and may be + //! changed + //! or removed without notice. + //! @since New in %Cantera 3.0. + //! void jacobianNotCalculated() { m_jac_calculated = false; }; protected: diff --git a/src/zeroD/ReactorNet.cpp b/src/zeroD/ReactorNet.cpp index e6768f2d2b6..543979bca9c 100644 --- a/src/zeroD/ReactorNet.cpp +++ b/src/zeroD/ReactorNet.cpp @@ -784,10 +784,10 @@ void ReactorNet::preconditionerSetup(double t, double* y, double gamma) // Update network with adjusted state updateState(yCopy.data()); // Create jacobian triplet vector - vector> jac_vector; - buildJacobian(jac_vector); + vector> jacVector; + buildJacobian(jacVector); // Add to preconditioner with offset - for (auto it : jac_vector) { + for (auto it : jacVector) { precon->setValue(it.row(), it.col(), it.value()); } // post reactor setup operations From 9d6077b185a3955c491a49c10c43a6f274426fb6 Mon Sep 17 00:00:00 2001 From: Anthony Walker Date: Tue, 24 Oct 2023 15:14:56 -0400 Subject: [PATCH 07/13] Add tests, fix bugs, remove added unused code in cython and cpp --- include/cantera/zeroD/FlowDevice.h | 4 +- .../zeroD/IdealGasConstPressureMoleReactor.h | 6 -- include/cantera/zeroD/IdealGasMoleReactor.h | 6 -- interfaces/cython/cantera/_utils.pxd | 3 - interfaces/cython/cantera/_utils.pyx | 31 --------- interfaces/cython/cantera/delegator.pyx | 2 +- src/zeroD/Reactor.cpp | 2 +- src/zeroD/ReactorNet.cpp | 2 +- test/python/test_reactor.py | 69 ++++++++++++++++++- test/zeroD/test_zeroD.cpp | 37 ++++++++++ 10 files changed, 111 insertions(+), 51 deletions(-) diff --git a/include/cantera/zeroD/FlowDevice.h b/include/cantera/zeroD/FlowDevice.h index 8662e40cf04..79a07aee193 100644 --- a/include/cantera/zeroD/FlowDevice.h +++ b/include/cantera/zeroD/FlowDevice.h @@ -166,7 +166,9 @@ class FlowDevice : public ConnectorNode //! @since New in %Cantera 3.0. //! virtual void buildNetworkJacobian(vector>& jacVector) { - throw NotImplementedError(type() + "::buildNetworkJacobian"); + if (!m_jac_calculated) { + throw NotImplementedError(type() + "::buildNetworkJacobian"); + } } //! Specify the jacobian terms have been calculated and should not be recalculated. diff --git a/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h b/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h index 2bb92fa35f1..1ad790f9586 100644 --- a/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h +++ b/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h @@ -43,12 +43,6 @@ class IdealGasConstPressureMoleReactor : public ConstPressureMoleReactor bool preconditionerSupported() const override { return true; }; - double temperatureDerivative() override { return 1; }; - - double temperatureRadiationDerivative() override { - return 4 * std::pow(temperature(), 3); - } - double moleDerivative(size_t index) override; double moleRadiationDerivative(size_t index) override; diff --git a/include/cantera/zeroD/IdealGasMoleReactor.h b/include/cantera/zeroD/IdealGasMoleReactor.h index cd30ddf02d4..ed1a0fb905b 100644 --- a/include/cantera/zeroD/IdealGasMoleReactor.h +++ b/include/cantera/zeroD/IdealGasMoleReactor.h @@ -45,12 +45,6 @@ class IdealGasMoleReactor : public MoleReactor bool preconditionerSupported() const override {return true;}; - double temperatureDerivative() override { return 1; }; - - double temperatureRadiationDerivative() override { - return 4 * std::pow(temperature(), 3); - } - double moleDerivative(size_t index) override; double moleRadiationDerivative(size_t index) override; diff --git a/interfaces/cython/cantera/_utils.pxd b/interfaces/cython/cantera/_utils.pxd index 25707835407..3949a296d40 100644 --- a/interfaces/cython/cantera/_utils.pxd +++ b/interfaces/cython/cantera/_utils.pxd @@ -120,7 +120,4 @@ cdef anymap_to_py(CxxAnyMap& m) cdef CxxAnyValue python_to_anyvalue(item, name=*) except * cdef anyvalue_to_python(string name, CxxAnyValue& v) -cdef vector[CxxEigenTriplet] python_to_triplets(triplets) except * -cdef triplets_to_python(vector[CxxEigenTriplet]& triplets) - cdef CxxEigenTriplet get_triplet(row, col, val) except * diff --git a/interfaces/cython/cantera/_utils.pyx b/interfaces/cython/cantera/_utils.pyx index 62216e67e88..d9d9ce50dd6 100644 --- a/interfaces/cython/cantera/_utils.pyx +++ b/interfaces/cython/cantera/_utils.pyx @@ -528,41 +528,10 @@ cdef vector[vector[string]] list2_string_to_anyvalue(data): v[i][j] = stringify(jtem) return v -cdef vector[CxxEigenTriplet] python_to_triplets(triplets): - # check that triplet dimensions are met - trip_shape = np.shape(triplets) - assert len(trip_shape) == 2 - assert trip_shape[1] == 3 - cdef vector[CxxEigenTriplet] trips - # c++ variables - cdef size_t row - cdef size_t col - cdef double val - cdef CxxEigenTriplet* trip_ptr - for r, c, v in triplets: - row = r - col = c - val = v - trip_ptr = new CxxEigenTriplet(row, col, val) - trips.push_back(dereference(trip_ptr)) - return trips - cdef CxxEigenTriplet get_triplet(row, col, val): cdef CxxEigenTriplet* trip_ptr = new CxxEigenTriplet(row, col, val) return dereference(trip_ptr) -cdef triplets_to_python(vector[CxxEigenTriplet]& triplets): - values = [] - cdef size_t row - cdef size_t col - cdef double val - for t in triplets: - row = t.row() - col = t.col() - val = t.value() - values.append((row, col, val)) - return values - def _py_to_any_to_py(dd): # used for internal testing purposes only cdef string name = stringify("test") diff --git a/interfaces/cython/cantera/delegator.pyx b/interfaces/cython/cantera/delegator.pyx index 0dab6e43621..2d6d2471c05 100644 --- a/interfaces/cython/cantera/delegator.pyx +++ b/interfaces/cython/cantera/delegator.pyx @@ -9,7 +9,7 @@ from libc.stdlib cimport malloc from libc.string cimport strcpy from ._utils import CanteraError -from ._utils cimport stringify, pystr, anymap_to_py, py_to_anymap, triplets_to_python, python_to_triplets, get_triplet +from ._utils cimport stringify, pystr, anymap_to_py, py_to_anymap, get_triplet from .units cimport Units, UnitStack # from .reaction import ExtensibleRate, ExtensibleRateData from .reaction cimport (ExtensibleRate, ExtensibleRateData, CxxReaction, diff --git a/src/zeroD/Reactor.cpp b/src/zeroD/Reactor.cpp index 2b6e78d0f5c..812f82bfd69 100644 --- a/src/zeroD/Reactor.cpp +++ b/src/zeroD/Reactor.cpp @@ -658,7 +658,7 @@ void Reactor::buildFlowJacobian(vector>& jacVector) m_outlet[i]->buildReactorJacobian(this, jacVector); } - for (size_t i = 0; i buildReactorJacobian(this, jacVector); } } diff --git a/src/zeroD/ReactorNet.cpp b/src/zeroD/ReactorNet.cpp index 543979bca9c..cf8fe0c3797 100644 --- a/src/zeroD/ReactorNet.cpp +++ b/src/zeroD/ReactorNet.cpp @@ -872,7 +872,7 @@ void ReactorNet::buildJacobian(vector>& jacVector) } } // flow devices - if (m_jac_skip_flow_devices) { + if (!m_jac_skip_flow_devices) { // outlets for (size_t i = 0; i < r->nOutlets(); i++) { r->outlet(i).buildNetworkJacobian(jacVector); diff --git a/test/python/test_reactor.py b/test/python/test_reactor.py index 2ac1f4bd17e..2545f9018b3 100644 --- a/test/python/test_reactor.py +++ b/test/python/test_reactor.py @@ -170,6 +170,61 @@ def test_finite_difference_jacobian(self): if name in constant: assert all(J[i, species_start:] == 0), (i, name) + def test_network_finite_difference_jacobian(self): + self.make_reactors(T1=900, P1=101325, X1="H2:0.4, O2:0.4, N2:0.2") + k1H2 = self.gas1.species_index("H2") + k2H2 = self.gas1.species_index("H2") + while self.r1.phase.X[k1H2] > 0.3 or self.r2.phase.X[k2H2] > 0.3: + self.net.step() + + J = self.net.finite_difference_jacobian + assert J.shape == (self.net.n_vars, self.net.n_vars) + + # state variables that should be constant, depending on reactor type + constant = {"mass", "volume", "int_energy", "enthalpy", "pressure"} + variable = {"temperature"} + for i in range(3): + name = self.r1.component_name(i) + if name in constant: + assert all(J[i,:] == 0), (i, name) + elif name in variable: + assert any(J[i,:] != 0) + # check in second reactor + name = self.r2.component_name(i) + if name in constant: + assert all(J[i + self.r1.n_vars,:] == 0), (i, name) + elif name in variable: + assert any(J[i + self.r1.n_vars,:] != 0) + + # Disabling energy equation should zero these terms + self.r1.energy_enabled = False + self.r2.energy_enabled = False + J = self.net.finite_difference_jacobian + for i in range(3): + name = self.r1.component_name(i) + if name == "temperature": + assert all(J[i,:] == 0) + name = self.r2.component_name(i) + if name == "temperature": + assert all(J[i + self.r1.n_vars,:] == 0) + + # Disabling species equations should zero these terms + self.r1.energy_enabled = True + self.r1.chemistry_enabled = False + self.r2.energy_enabled = True + self.r2.chemistry_enabled = False + J = self.net.finite_difference_jacobian + constant = set(self.r1.phase.species_names + self.r2.phase.species_names) + r1_species_start = self.r1.component_index(self.r1.phase.species_name(0)) + r2_species_start = self.r2.component_index(self.r2.phase.species_name(0)) + for i in range(self.r1.n_vars): + name = self.r1.component_name(i) + if name in constant: + assert all(J[i, r1_species_start:] == 0), (i, name) + name = self.r2.component_name(i) + if name in constant: + assert all(J[i + self.r1.n_vars, (r2_species_start + self.r1.n_vars):] == 0), (i, name) + def test_timestepping(self): self.make_reactors() @@ -1530,7 +1585,7 @@ def create_reactors(self, **kwargs): self.precon = ct.AdaptivePreconditioner() self.net2.preconditioner = self.precon self.net2.derivative_settings = {"skip-third-bodies":True, "skip-falloff":True, - "skip-coverage-dependence":True} + "skip-coverage-dependence":True, "skip-flow-devices": True} def test_get_solver_type(self): self.create_reactors() @@ -1538,6 +1593,18 @@ def test_get_solver_type(self): self.net2.initialize() assert self.net2.linear_solver_type == "GMRES" + def test_mass_flow_jacobian(self): + self.create_reactors(add_mdot=True) + # reset derivative settings + self.net2.derivative_settings = {"skip-third-bodies":True, "skip-falloff":True, + "skip-coverage-dependence":True, "skip-flow-devices": False} + + with pytest.raises(NotImplementedError, match="MassFlowController::buildReactorJacobian"): + J = self.net2.jacobian + + with pytest.raises(NotImplementedError, match="MassFlowController::buildReactorJacobian"): + J = self.r2.jacobian + def test_heat_transfer_network(self): # create first reactor gas1 = ct.Solution("h2o2.yaml", "ohmech") diff --git a/test/zeroD/test_zeroD.cpp b/test/zeroD/test_zeroD.cpp index 343aaa3a86d..56d9f226095 100644 --- a/test/zeroD/test_zeroD.cpp +++ b/test/zeroD/test_zeroD.cpp @@ -357,6 +357,12 @@ TEST(JacobianTests, test_wall_jacobian_build) EXPECT_LT(it.col(), reactor2.neq()); } } + // check that switch works + wallJac.clear(); + w.jacobianCalculated(); + w.buildNetworkJacobian(wallJac); + EXPECT_EQ(wallJac.size(), 0); + w.jacobianNotCalculated(); // build jac for network terms wallJac.clear(); w.buildNetworkJacobian(wallJac); @@ -381,6 +387,37 @@ TEST(JacobianTests, test_wall_jacobian_build) } } +TEST(JacobianTests, test_flow_jacobian_build) +{ + // create reservoir reactor + auto sol = newSolution("h2o2.yaml"); + sol->thermo()->setState_TPY(1000.0, OneAtm, "O2:1.0"); + Reservoir res; + res.insert(sol); + // create reactor + IdealGasConstPressureMoleReactor reactor; + reactor.insert(sol); + reactor.setInitialVolume(1.0); + // create the flow device + MassFlowController mfc; + mfc.install(res, reactor); + mfc.setMassFlowCoeff(1.0); + // setup reactor network and integrate + ReactorNet network; + network.addReactor(reactor); + network.initialize(); + // manually get wall jacobian elements + vector> flowJac; + // expect errors from building jacobians + EXPECT_THROW(mfc.buildReactorJacobian(&reactor, flowJac), NotImplementedError); + // check the jacobian calculated flag and throw/catch errors accordingly + mfc.jacobianCalculated(); + mfc.buildNetworkJacobian(flowJac); + EXPECT_EQ(flowJac.size(), 0); + mfc.jacobianNotCalculated(); + EXPECT_THROW(mfc.buildNetworkJacobian(flowJac), NotImplementedError); +} + int main(int argc, char** argv) { printf("Running main() from test_zeroD.cpp\n"); From 17857dfe3f3e925a0cf6437c1e7970011a684e84 Mon Sep 17 00:00:00 2001 From: Anthony Walker Date: Mon, 1 Jan 2024 17:32:43 -0500 Subject: [PATCH 08/13] Add a default evaluation function for reactors It can sometimes be useful to add on to the default evaluation in such a way that the standard before or after delegations cannot handle. A default evaluation function allows for this. --- include/cantera/zeroD/ReactorDelegator.h | 8 ++++++++ interfaces/cython/cantera/reactor.pxd | 1 + interfaces/cython/cantera/reactor.pyx | 12 ++++++++++++ test/python/test_reactor.py | 23 +++++++++++++++++++++++ 4 files changed, 44 insertions(+) diff --git a/include/cantera/zeroD/ReactorDelegator.h b/include/cantera/zeroD/ReactorDelegator.h index 1790331c100..25069427b49 100644 --- a/include/cantera/zeroD/ReactorDelegator.h +++ b/include/cantera/zeroD/ReactorDelegator.h @@ -52,6 +52,10 @@ class ReactorAccessor //! Set the state of the thermo object for surface *n* to correspond to the //! state of that surface virtual void restoreSurfaceState(size_t n) = 0; + + //! Public access to the default evaluation function so it can be used in + //! replace functions + virtual void defaultEval(double t, double* LHS, double* RHS) = 0; }; //! Delegate methods of the Reactor class to external functions @@ -167,6 +171,10 @@ class ReactorDelegator : public Delegator, public R, public ReactorAccessor m_build_jacobian(jacVector); } + void defaultEval(double t, double* LHS, double* RHS) override { + R::eval(t, LHS, RHS); + } + // Public access to protected Reactor variables needed by derived classes void setNEq(size_t n) override { diff --git a/interfaces/cython/cantera/reactor.pxd b/interfaces/cython/cantera/reactor.pxd index 5a1ba565fcc..d22a72a293f 100644 --- a/interfaces/cython/cantera/reactor.pxd +++ b/interfaces/cython/cantera/reactor.pxd @@ -220,6 +220,7 @@ cdef extern from "cantera/zeroD/ReactorDelegator.h" namespace "Cantera": void setHeatRate(double) void restoreThermoState() except +translate_exception void restoreSurfaceState(size_t) except +translate_exception + void defaultEval(double time, double* LHS, double* RHS) ctypedef CxxReactorAccessor* CxxReactorAccessorPtr diff --git a/interfaces/cython/cantera/reactor.pyx b/interfaces/cython/cantera/reactor.pyx index 80c707ffda0..1b7a87fb0da 100644 --- a/interfaces/cython/cantera/reactor.pyx +++ b/interfaces/cython/cantera/reactor.pyx @@ -2,6 +2,8 @@ # at https://cantera.org/license.txt for license and copyright information. import warnings +import numpy as np +from collections import defaultdict as _defaultdict import numbers as _numbers from cython.operator cimport dereference as deref import numpy as np @@ -705,6 +707,16 @@ cdef class ExtensibleReactor(Reactor): """ self.accessor.restoreSurfaceState(n) + def default_eval(self, time, LHS, RHS): + """ + Evaluation of the base reactors `eval` function to be used in `replace` + functions and maintain original functionality. + """ + assert len(LHS) == self.n_vars and len(RHS) == self.n_vars + cdef np.ndarray[np.double_t, ndim=1, mode="c"] rhs = np.frombuffer(RHS) + cdef np.ndarray[np.double_t, ndim=1, mode="c"] lhs = np.frombuffer(LHS) + self.accessor.defaultEval(time, &lhs[0], &rhs[0]) + cdef class ExtensibleIdealGasReactor(ExtensibleReactor): """ diff --git a/test/python/test_reactor.py b/test/python/test_reactor.py index 2545f9018b3..e17da1beb2e 100644 --- a/test/python/test_reactor.py +++ b/test/python/test_reactor.py @@ -3561,6 +3561,29 @@ def replace_build_jacobian(self, jac_vector): jac = r.jacobian assert np.sum(jac) == 0 + def test_replace_with_default_eval(self): + class ReplaceEvalReactor(ct.ExtensibleIdealGasConstPressureMoleReactor): + + def replace_eval(self, t, LHS, RHS): + self.default_eval(t, LHS, RHS) + + # setup thermo object + gas = ct.Solution("h2o2.yaml", "ohmech") + gas.set_equivalence_ratio(0.5, "H2:1.0", "O2:1.0") + gas.equilibrate("HP") + # replacement reactor + r = ReplaceEvalReactor(gas) + r.volume = 1.0 + # default reactor + rstd = ct.IdealGasConstPressureMoleReactor(gas) + rstd.volume = r.volume + # network of both reactors + net = ct.ReactorNet([r, rstd]) + net.preconditioner = ct.AdaptivePreconditioner() + net.advance_to_steady_state() + # reactors should have the same solution because the default is used + self.assertArrayNear(r.get_state(), rstd.get_state()) + class TestSteadySolver: @pytest.mark.parametrize("reactor_class", From 8dfd7edf07153604a8a0f1b28bcf1af2f33dfa32 Mon Sep 17 00:00:00 2001 From: Anthony Walker Date: Sun, 3 Mar 2024 10:00:32 -0500 Subject: [PATCH 09/13] Working on PR updates, have to fix a few more bugs in regards to equations, handling volumes, etc. --- doc/sphinx/develop/index.md | 9 ++- include/cantera/zeroD/FlowDevice.h | 26 ++----- .../zeroD/IdealGasConstPressureMoleReactor.h | 4 +- include/cantera/zeroD/IdealGasMoleReactor.h | 4 +- include/cantera/zeroD/Reactor.h | 10 +-- include/cantera/zeroD/ReactorBase.h | 51 ++----------- include/cantera/zeroD/ReactorDelegator.h | 7 +- include/cantera/zeroD/ReactorNet.h | 17 ++--- include/cantera/zeroD/Wall.h | 50 ++++--------- interfaces/cython/cantera/delegator.pyx | 1 + .../IdealGasConstPressureMoleReactor.cpp | 18 +---- src/zeroD/IdealGasMoleReactor.cpp | 22 +----- src/zeroD/Reactor.cpp | 11 ++- src/zeroD/ReactorNet.cpp | 75 +++++++++---------- src/zeroD/Wall.cpp | 32 ++++---- test/python/test_reactor.py | 36 ++++----- test/zeroD/test_zeroD.cpp | 32 +++----- 17 files changed, 138 insertions(+), 267 deletions(-) diff --git a/doc/sphinx/develop/index.md b/doc/sphinx/develop/index.md index 8b19f523289..2cf48501711 100644 --- a/doc/sphinx/develop/index.md +++ b/doc/sphinx/develop/index.md @@ -1,6 +1,7 @@ # Develop (sec-compiling)= + ## Compiling Cantera from Source If you're interested in contributing new features to Cantera, or you want to try the @@ -12,9 +13,9 @@ Cantera](compiling/configure-build) on your computer. The following additional references may also be useful: -- [](compiling/dependencies.md) -- [](compiling/config-options) -- [](compiling/special-cases) +- [](compiling/dependencies.md) +- [](compiling/config-options) +- [](compiling/special-cases) ```{toctree} :caption: Compiling Cantera from Source @@ -35,7 +36,7 @@ compiling/special-cases This section is a work in progress. ``` -- [](reactor-integration) +- [](reactor-integration) ```{toctree} :caption: How Cantera Works diff --git a/include/cantera/zeroD/FlowDevice.h b/include/cantera/zeroD/FlowDevice.h index 79a07aee193..c567d25e0d7 100644 --- a/include/cantera/zeroD/FlowDevice.h +++ b/include/cantera/zeroD/FlowDevice.h @@ -149,9 +149,10 @@ class FlowDevice : public ConnectorNode //! @warning This function is an experimental part of the %Cantera API and may be //! changed //! or removed without notice. - //! @since New in %Cantera 3.0. + //! @since New in %Cantera 3.1. //! - virtual void buildReactorJacobian(ReactorBase* r, vector>& jacVector) { + virtual void buildReactorJacobian(ReactorBase* r, + vector>& jacVector) { throw NotImplementedError(type() + "::buildReactorJacobian"); } @@ -163,29 +164,12 @@ class FlowDevice : public ConnectorNode //! @warning This function is an experimental part of the %Cantera API and may be //! changed //! or removed without notice. - //! @since New in %Cantera 3.0. + //! @since New in %Cantera 3.1. //! virtual void buildNetworkJacobian(vector>& jacVector) { - if (!m_jac_calculated) { - throw NotImplementedError(type() + "::buildNetworkJacobian"); - } + throw NotImplementedError(type() + "::buildNetworkJacobian"); } - //! Specify the jacobian terms have been calculated and should not be recalculated. - //! @warning This function is an experimental part of the %Cantera API and may be - //! changed - //! or removed without notice. - //! @since New in %Cantera 3.0. - //! - void jacobianCalculated() { m_jac_calculated = true; }; - - //! Specify that jacobian terms have not been calculated and should be recalculated. - //! @warning This function is an experimental part of the %Cantera API and may be changed - //! or removed without notice. - //! @since New in %Cantera 3.0. - //! - void jacobianNotCalculated() { m_jac_calculated = false; }; - protected: string m_name; //!< Flow device name. bool m_defaultNameSet = false; //!< `true` if default name has been previously set. diff --git a/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h b/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h index 1ad790f9586..62f75206a38 100644 --- a/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h +++ b/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h @@ -43,9 +43,7 @@ class IdealGasConstPressureMoleReactor : public ConstPressureMoleReactor bool preconditionerSupported() const override { return true; }; - double moleDerivative(size_t index) override; - - double moleRadiationDerivative(size_t index) override; + double temperature_ddni(size_t index) override; size_t speciesOffset() const override { return m_sidx; }; diff --git a/include/cantera/zeroD/IdealGasMoleReactor.h b/include/cantera/zeroD/IdealGasMoleReactor.h index ed1a0fb905b..48b6a0ccd66 100644 --- a/include/cantera/zeroD/IdealGasMoleReactor.h +++ b/include/cantera/zeroD/IdealGasMoleReactor.h @@ -45,9 +45,7 @@ class IdealGasMoleReactor : public MoleReactor bool preconditionerSupported() const override {return true;}; - double moleDerivative(size_t index) override; - - double moleRadiationDerivative(size_t index) override; + double temperature_ddni(size_t index) override; size_t speciesOffset() const override { return m_sidx; }; diff --git a/include/cantera/zeroD/Reactor.h b/include/cantera/zeroD/Reactor.h index 0ffecc16d67..14e9124e99e 100644 --- a/include/cantera/zeroD/Reactor.h +++ b/include/cantera/zeroD/Reactor.h @@ -207,15 +207,7 @@ class Reactor : public ReactorBase //! //! @warning This method is an experimental part of the %Cantera //! API and may be changed or removed without notice. - virtual Eigen::SparseMatrix jacobian() { - m_jac_trips.clear(); - // Add before, during, after evals - buildJacobian(m_jac_trips); - // construct jacobian from vector - Eigen::SparseMatrix jac(m_nv, m_nv); - jac.setFromTriplets(m_jac_trips.begin(), m_jac_trips.end()); - return jac; - } + virtual Eigen::SparseMatrix jacobian(); //! Calculate the Jacobian of a specific Reactor specialization. //! @param jacVector vector where jacobian triplets are added diff --git a/include/cantera/zeroD/ReactorBase.h b/include/cantera/zeroD/ReactorBase.h index 9d70d634acf..43cc1bbd9d6 100644 --- a/include/cantera/zeroD/ReactorBase.h +++ b/include/cantera/zeroD/ReactorBase.h @@ -302,52 +302,15 @@ class ReactorBase return m_sensParams.size(); } - //! Calculate the derivative of temperature with respect to the temperature in the - //! heat transfer equation based on the reactor specific equation of state. - //! This function should also transform the state of the derivative to that - //! appropriate for the jacobian's state/ - //! @warning This function is an experimental part of the %Cantera API and may be changed - //! or removed without notice. - //! @since New in %Cantera 3.0. - //! - virtual double temperatureDerivative() { - throw NotImplementedError("Reactor::temperatureDerivative"); - } - - //! Calculate the derivative of temperature with respect to the temperature in the - //! heat transfer radiation equation based on the reactor specific equation of - //! state. - //! This function should also transform the state of the derivative to that - //! appropriate for the jacobian's state/ - //! @warning This function is an experimental part of the %Cantera API and may be - //! changed - //! or removed without notice. - //! @since New in %Cantera 3.0. - //! - virtual double temperatureRadiationDerivative() { - throw NotImplementedError("Reactor::temperatureRadiationDerivative"); - } - - //! Calculate the derivative of T with respect to the ith species in the heat - //! transfer equation based on the reactor specific equation of state. + //! Calculate the derivative of T with respect to the ith species in the energy + //! conservation equation based on the reactor specific equation of state. //! @param index index of the species the derivative is with respect too - //! @warning This function is an experimental part of the %Cantera API and may be changed - //! or removed without notice. - //! @since New in %Cantera 3.0. - //! - virtual double moleDerivative(size_t index) { - throw NotImplementedError("Reactor::moleDerivative"); - } - - //! Calculate the derivative of T with respect to the ith species in the heat - //! transfer radiation equation based on the reactor specific equation of state. - //! @param index index of the species the derivative is with respect too - //! @warning This function is an experimental part of the %Cantera API and may be changed - //! or removed without notice. - //! @since New in %Cantera 3.0. + //! @warning This function is an experimental part of the %Cantera API and may + //! be changed or removed without notice. + //! @since New in %Cantera 3.1. //! - virtual double moleRadiationDerivative(size_t index) { - throw NotImplementedError("Reactor::moleRadiationDerivative"); + virtual double temperature_ddni(size_t index) { + throw NotImplementedError("Reactor::temperature_ddni"); } //! Return the index associated with energy of the system diff --git a/include/cantera/zeroD/ReactorDelegator.h b/include/cantera/zeroD/ReactorDelegator.h index 25069427b49..5860e6a7992 100644 --- a/include/cantera/zeroD/ReactorDelegator.h +++ b/include/cantera/zeroD/ReactorDelegator.h @@ -74,7 +74,8 @@ class ReactorDelegator : public Delegator, public R, public ReactorAccessor install("updateState", m_updateState, [this](std::array sizes, double* y) { R::updateState(y); }); install("updateSurfaceState", m_updateSurfaceState, - [this](std::array sizes, double* y) { R::updateSurfaceState(y); }); + [this](std::array sizes, double* y) + { R::updateSurfaceState(y); }); install("getSurfaceInitialConditions", m_getSurfaceInitialConditions, [this](std::array sizes, double* y) { R::getSurfaceInitialConditions(y); @@ -89,8 +90,8 @@ class ReactorDelegator : public Delegator, public R, public ReactorAccessor ); install("evalWalls", m_evalWalls, [this](double t) { R::evalWalls(t); }); install("evalSurfaces", m_evalSurfaces, - [this](std::array sizes, double* LHS, double* RHS, double* sdot) { - R::evalSurfaces(LHS, RHS, sdot); + [this](std::array sizes, double* LHS, double* RHS, double* sd) { + R::evalSurfaces(LHS, RHS, sd); } ); install("componentName", m_componentName, diff --git a/include/cantera/zeroD/ReactorNet.h b/include/cantera/zeroD/ReactorNet.h index 55901d2e6eb..631e80c35d1 100644 --- a/include/cantera/zeroD/ReactorNet.h +++ b/include/cantera/zeroD/ReactorNet.h @@ -7,6 +7,8 @@ #define CT_REACTORNET_H #include "Reactor.h" +#include "Wall.h" +#include "FlowDevice.h" #include "cantera/numerics/FuncEval.h" #include "cantera/numerics/SteadyStateSystem.h" @@ -284,16 +286,7 @@ class ReactorNet : public FuncEval //! Return the index corresponding to the start of the reactor specific state //! vector in the reactor with index *reactor* in the global state vector for the //! reactor network. - size_t globalStartIndex(Reactor* curr_reactor) { - for (size_t i = 0; i < m_reactors.size(); i++) { - if (curr_reactor == m_reactors[i]) { - return m_start[i]; - } - } - throw CanteraError("ReactorNet::globalStartIndex: ", - curr_reactor->name(), " not found in network."); - return npos; - } + size_t globalStartIndex(ReactorBase* curr_reactor); //! Return the name of the i-th component of the global state vector. The //! name returned includes both the name of the reactor and the specific @@ -498,6 +491,10 @@ class ReactorNet : public FuncEval //! derivative settings bool m_jac_skip_walls = false; bool m_jac_skip_flow_devices = false; + //! set to store walls for Jacobian calculation + set m_walls; + //! set to store flow devices for Jacobian calculation + set m_flow_devices; }; diff --git a/include/cantera/zeroD/Wall.h b/include/cantera/zeroD/Wall.h index 9fb79322452..9530ab6329d 100644 --- a/include/cantera/zeroD/Wall.h +++ b/include/cantera/zeroD/Wall.h @@ -89,47 +89,27 @@ class WallBase : public ConnectorNode //! Build the Jacobian terms specific to the flow device for the given connected //! reactor. //! @param r a pointer to the calling reactor - //! @param jacVector a vector of triplets to be added to the jacobian for the - //! reactor + //! @param jacVector a vector of triplets to be added to the reactor Jacobian //! @warning This function is an experimental part of the %Cantera API and may be - //! changed - //! or removed without notice. - //! @since New in %Cantera 3.0. + //! changed or removed without notice. + //! @since New in %Cantera 3.1. //! - virtual void buildReactorJacobian(ReactorBase* r, vector>& jacVector) { + virtual void buildReactorJacobian(ReactorBase* r, + vector>& jacVector) { throw NotImplementedError("WallBase::buildReactorJacobian"); } - //! Build the Jacobian terms specific to the flow device for the network. These - //! terms - //! will be adjusted to the networks indexing system outside of the reactor. - //! @param jacVector a vector of triplets to be added to the jacobian for the - //! reactor + //! Build the Jacobian cross-reactor terms specific to the flow device for the + //! network. + //! @param jacVector a vector of triplets to be added to the network Jacobian //! @warning This function is an experimental part of the %Cantera API and may be - //! changed - //! or removed without notice. - //! @since New in %Cantera 3.0. + //! changed or removed without notice. + //! @since New in %Cantera 3.1. //! virtual void buildNetworkJacobian(vector>& jacVector) { throw NotImplementedError("WallBase::buildNetworkJacobian"); } - //! Specify the jacobian terms have been calculated and should not be recalculated. - //! @warning This function is an experimental part of the %Cantera API and may be - //! changed - //! or removed without notice. - //! @since New in %Cantera 3.0. - //! - void jacobianCalculated() { m_jac_calculated = true; }; - - //! Specify that jacobian terms have not been calculated and should be recalculated. - //! @warning This function is an experimental part of the %Cantera API and may be - //! changed - //! or removed without notice. - //! @since New in %Cantera 3.0. - //! - void jacobianNotCalculated() { m_jac_calculated = false; }; - protected: ReactorBase* m_left = nullptr; ReactorBase* m_right = nullptr; @@ -138,10 +118,6 @@ class WallBase : public ConnectorNode double m_time = 0.0; double m_area = 1.0; - - //! a variable to switch on and off so calculations are not doubled by the calling - //! reactor or network - bool m_jac_calculated = false; }; //! Represents a wall between between two ReactorBase objects. @@ -249,9 +225,11 @@ class Wall : public WallBase return m_k; } - virtual void buildReactorJacobian(ReactorBase* r, vector>& jacVector) override; + void buildReactorJacobian(ReactorBase* r, + vector>& jacVector) override; - virtual void buildNetworkJacobian(vector>& jacVector) override; + void buildNetworkJacobian(vector>& jacVector) + override; protected: diff --git a/interfaces/cython/cantera/delegator.pyx b/interfaces/cython/cantera/delegator.pyx index 2d6d2471c05..2166b245f74 100644 --- a/interfaces/cython/cantera/delegator.pyx +++ b/interfaces/cython/cantera/delegator.pyx @@ -319,6 +319,7 @@ cdef int assign_delegates(obj, CxxDelegator* delegator) except -1: if when is None: continue + cxx_name = stringify(options[0]) callback = options[1].replace(' ', '') diff --git a/src/zeroD/IdealGasConstPressureMoleReactor.cpp b/src/zeroD/IdealGasConstPressureMoleReactor.cpp index 0b5b24f9684..93f7caf57ef 100644 --- a/src/zeroD/IdealGasConstPressureMoleReactor.cpp +++ b/src/zeroD/IdealGasConstPressureMoleReactor.cpp @@ -289,23 +289,11 @@ double IdealGasConstPressureMoleReactor::lowerBound(size_t k) const { } } -double IdealGasConstPressureMoleReactor::moleDerivative(size_t index) +double IdealGasConstPressureMoleReactor::temperature_ddni(size_t index) { // derivative of temperature transformed by ideal gas law - vector moles(m_nsp); - getMoles(moles.data()); - double dTdni = pressure() * m_vol / GasConstant / std::accumulate(moles.begin(), moles.end(), 0.0); - return dTdni; -} - -double IdealGasConstPressureMoleReactor::moleRadiationDerivative(size_t index) -{ - // derivative of temperature transformed by ideal gas law - vector moles(m_nsp); - getMoles(moles.data()); - double dT4dni = std::pow(pressure() * m_vol / GasConstant, 4); - dT4dni *= std::pow(1 / std::accumulate(moles.begin(), moles.end(), 0.0), 5); - return dT4dni; + double n_total = m_mass / m_thermo->meanMolecularWeight(); + return pressure() * m_vol / GasConstant / n_total; } } diff --git a/src/zeroD/IdealGasMoleReactor.cpp b/src/zeroD/IdealGasMoleReactor.cpp index ec3e25d39e3..0b9e0ef78b1 100644 --- a/src/zeroD/IdealGasMoleReactor.cpp +++ b/src/zeroD/IdealGasMoleReactor.cpp @@ -274,32 +274,16 @@ void IdealGasMoleReactor::buildJacobian(vector>& jacVecto jacVector.emplace_back(0, static_cast(j + m_sidx), (specificHeat[j] * qdot - NCv * uk_dnkdnj_sums[j]) * denom); } - - // build wall jacobian buildWallJacobian(jacVector); } - - // build flow jacobian buildFlowJacobian(jacVector); } -double IdealGasMoleReactor::moleDerivative(size_t index) -{ - // derivative of temperature transformed by ideal gas law - vector moles(m_nsp); - getMoles(moles.data()); - double dTdni = pressure() * m_vol / GasConstant / std::accumulate(moles.begin(), moles.end(), 0.0); - return dTdni; -} - -double IdealGasMoleReactor::moleRadiationDerivative(size_t index) +double IdealGasMoleReactor::temperature_ddni(size_t index) { // derivative of temperature transformed by ideal gas law - vector moles(m_nsp); - getMoles(moles.data()); - double dT4dni = std::pow(pressure() * m_vol / GasConstant, 4); - dT4dni *= std::pow(1 / std::accumulate(moles.begin(), moles.end(), 0.0), 5); - return dT4dni; + double n_total = m_mass / m_thermo->meanMolecularWeight(); + return pressure() * m_vol / GasConstant / n_total; } } diff --git a/src/zeroD/Reactor.cpp b/src/zeroD/Reactor.cpp index 812f82bfd69..1e3d337f3d8 100644 --- a/src/zeroD/Reactor.cpp +++ b/src/zeroD/Reactor.cpp @@ -642,7 +642,7 @@ void Reactor::setAdvanceLimit(const string& nm, const double limit) } } -void Reactor:: buildWallJacobian(vector>& jacVector) +void Reactor::buildWallJacobian(vector>& jacVector) { if (!m_jac_skip_walls) { for (size_t i = 0; i < m_wall.size(); i++) { @@ -664,4 +664,13 @@ void Reactor::buildFlowJacobian(vector>& jacVector) } } +Eigen::SparseMatrix Reactor::jacobian() { + m_jac_trips.clear(); + // Add before, during, after evals + buildJacobian(m_jac_trips); + // construct jacobian from vector + Eigen::SparseMatrix jac(m_nv, m_nv); + jac.setFromTriplets(m_jac_trips.begin(), m_jac_trips.end()); + return jac; + } } diff --git a/src/zeroD/ReactorNet.cpp b/src/zeroD/ReactorNet.cpp index cf8fe0c3797..efea6ff1ee6 100644 --- a/src/zeroD/ReactorNet.cpp +++ b/src/zeroD/ReactorNet.cpp @@ -159,7 +159,26 @@ void ReactorNet::initialize() " when creating the Reactor or ReactorSurface objects.", shared); } } - + // Create walls and flow devices sets + for (auto r : m_reactors) { + // walls + if (!m_jac_skip_walls) { + for (size_t i = 0; i < r->nWalls(); i++) { + m_walls.insert(&(r->wall(i))); + } + } + // flow devices + if (!m_jac_skip_flow_devices) { + // outlets + for (size_t i = 0; i < r->nOutlets(); i++) { + m_flow_devices.insert(&(r->outlet(i))); + } + // inlets + for (size_t i = 0; i < r->nInlets(); i++) { + m_flow_devices.insert(&(r->inlet(i))); + } + } + } m_ydot.resize(m_nv,0.0); m_yest.resize(m_nv,0.0); m_advancelimits.resize(m_nv,-1.0); @@ -720,6 +739,16 @@ size_t ReactorNet::registerSensitivityParameter( return m_sens_params.size() - 1; } +size_t ReactorNet::globalStartIndex(ReactorBase* curr_reactor) { + for (size_t i = 0; i < m_reactors.size(); i++) { + if (curr_reactor == m_reactors[i]) { + return m_start[i]; + } + } + throw CanteraError("ReactorNet::globalStartIndex: ", + curr_reactor->name(), " not found in network."); + } + void ReactorNet::setDerivativeSettings(AnyMap& settings) { // Apply given settings to all reactors @@ -820,29 +849,9 @@ void ReactorNet::checkPreconditionerSupported() const { void ReactorNet::buildJacobian(vector>& jacVector) { // network must be initialized for the jacobian - if (! m_init) { + if (!m_init) { initialize(); } - // loop through and set connectors not found - for (auto r : m_reactors) { - // walls - if (!m_jac_skip_walls) { - for (size_t i = 0; i < r->nWalls(); i++) { - r->wall(i).jacobianNotCalculated(); - } - } - // flow devices - if (!m_jac_skip_flow_devices) { - // outlets - for (size_t i = 0; i < r->nOutlets(); i++) { - r->outlet(i).jacobianNotCalculated(); - } - // inlets - for (size_t i = 0; i < r->nInlets(); i++) { - r->inlet(i).jacobianNotCalculated(); - } - } - } // Create jacobian triplet vector vector jstarts; // Get jacobians and give elements to preconditioners @@ -865,24 +874,12 @@ void ReactorNet::buildJacobian(vector>& jacVector) // repeated for (auto r : m_reactors) { // walls - if (!m_jac_skip_walls) { - for (size_t i = 0; i < r->nWalls(); i++) { - r->wall(i).buildNetworkJacobian(jacVector); - r->wall(i).jacobianCalculated(); - } + for (const auto& wall : m_walls) { + wall->buildNetworkJacobian(jacVector); } // flow devices - if (!m_jac_skip_flow_devices) { - // outlets - for (size_t i = 0; i < r->nOutlets(); i++) { - r->outlet(i).buildNetworkJacobian(jacVector); - r->outlet(i).jacobianCalculated(); - } - // inlets - for (size_t i = 0; i < r->nInlets(); i++) { - r->inlet(i).buildNetworkJacobian(jacVector); - r->inlet(i).jacobianCalculated(); - } + for (const auto& flow_device : m_flow_devices) { + flow_device->buildNetworkJacobian(jacVector); } } } @@ -894,7 +891,7 @@ Eigen::SparseMatrix ReactorNet::finiteDifferenceJacobian() initialize(); } - // clear former jacobian elements + // allocate jacobian triplet vector vector> jac_trips; // Get the current state diff --git a/src/zeroD/Wall.cpp b/src/zeroD/Wall.cpp index b02e8b10c35..e511d776949 100644 --- a/src/zeroD/Wall.cpp +++ b/src/zeroD/Wall.cpp @@ -78,7 +78,8 @@ double Wall::heatRate() } -void Wall::buildReactorJacobian(ReactorBase* r, vector>& jacVector) +void Wall::buildReactorJacobian(ReactorBase* r, + vector>& jacVector) { // get derivative of heat transfer for both reactors vector> network; @@ -86,12 +87,13 @@ void Wall::buildReactorJacobian(ReactorBase* r, vector>& size_t sidx = r->speciesOffset(); size_t eidx = r->energyIndex(); // define a scalar for direction based on left and right - double scalar = (r == m_left) ? 1.0 : -1.0; + double direction = (r == m_left) ? 1.0 : -1.0; // elements within the current reactor // find dQdni for the current reactor w.r.t current reactor for (size_t i = sidx; i < nsp + sidx; i++) { - double dQdni = m_rrth * m_area * scalar * r->moleDerivative(i); - dQdni += m_emiss * m_area * scalar * r->moleRadiationDerivative(i); + double dQdni = m_rrth * m_area * direction * r->temperature_ddni(i); + dQdni += m_emiss * m_area * direction * r->temperature_ddni(i) * 4 + * pow(r->temperature(), 3); jacVector.emplace_back(eidx, i, dQdni); } } @@ -102,38 +104,36 @@ void Wall::buildNetworkJacobian(vector>& jacVector) if (m_right->type() == "Reservoir" || m_left->type() == "Reservoir") { return; } - // return if the jacobian has been calculated - if (m_jac_calculated) { - return; - } // get derivatives for inter-dependent reactor terms //variables for the right side vector> network; size_t r_nsp = m_right->phase()->thermo()->nSpecies(); size_t r_sidx = m_right->speciesOffset(); - size_t r_net = m_right->network().globalStartIndex((Reactor* ) m_right); + size_t r_net = m_right->network().globalStartIndex(m_right); size_t r_eidx = m_right->energyIndex(); // variables for the left side size_t l_nsp = m_left->phase()->thermo()->nSpecies(); size_t l_sidx = m_left->speciesOffset(); - size_t l_net = m_left->network().globalStartIndex((Reactor* ) m_left); + size_t l_net = m_left->network().globalStartIndex(m_left); size_t l_eidx = m_left->energyIndex(); - if (((Reactor* ) m_right)->energyEnabled()) { + if (m_right->energyEnabled()) { // find dQdni for the right reactor w.r.t left reactor for (size_t i = l_sidx; i < l_sidx + l_nsp; i++) { - double dQdni = m_rrth * m_area * m_left->moleDerivative(i); - dQdni += m_emiss * m_area * m_left->moleRadiationDerivative(i); + double dQdni = m_rrth * m_area * m_left->temperature_ddni(i); + dQdni += m_emiss * m_area * m_left->temperature_ddni(i) * 4 + * pow(m_left->temperature(), 3); jacVector.emplace_back(r_eidx + r_net, i + l_net, dQdni); } } - if (((Reactor* ) m_left)->energyEnabled()) { + if (m_left->energyEnabled()) { // find dQdni for the left reactor w.r.t right reactor for (size_t i = r_sidx; i < r_sidx + r_nsp; i++) { - double dQdni = - m_rrth * m_area * m_right->moleDerivative(i); - dQdni -= m_emiss * m_area * m_right->moleRadiationDerivative(i); + double dQdni = - m_rrth * m_area * m_right->temperature_ddni(i); + dQdni -= m_emiss * m_area * m_right->temperature_ddni(i) * 4 + * pow(m_right->temperature(), 3); jacVector.emplace_back(l_eidx + l_net, i + r_net, dQdni); } } diff --git a/test/python/test_reactor.py b/test/python/test_reactor.py index e17da1beb2e..aafd1ab279e 100644 --- a/test/python/test_reactor.py +++ b/test/python/test_reactor.py @@ -1605,51 +1605,44 @@ def test_mass_flow_jacobian(self): with pytest.raises(NotImplementedError, match="MassFlowController::buildReactorJacobian"): J = self.r2.jacobian + @pytest.mark.xfail def test_heat_transfer_network(self): # create first reactor gas1 = ct.Solution("h2o2.yaml", "ohmech") gas1.TPX = 600, ct.one_atm, "O2:1.0" r1 = self.reactorClass(gas1) - r1.chemistry_enabled = False # create second reactor gas2 = ct.Solution("h2o2.yaml", "ohmech") gas2.TPX = 300, ct.one_atm, "O2:1.0" - r2 = self.reactorClass(gas2) - r2.chemistry_enabled = False + r2 = ct.reactorClass(gas2) # create wall - U = 2.0 + U = 2. A = 3.0 w = ct.Wall(r1, r2, U=U, A=A) - net = ct.ReactorNet([r1, r2]) + net = ct.ReactorNet([r1,]) jac = net.jacobian - - # check for values - for i in range(gas1.n_species): - assert np.isclose(jac[0, i + 1 + r1.n_vars], - U * A * gas2.T) - assert np.isclose(jac[r1.n_vars, i + 1], U * A * gas1.T) - if (i + 1 != 2): - assert np.isclose(jac[0, i + 1], U * A * gas1.T, rtol=1e-2) - assert np.isclose(jac[r1.n_vars, i + 1 + r1.n_vars], - U * A * gas2.T, rtol=1e-2) - + fd_jac = net.finite_difference_jacobian + for i in range(jac.shape[0]): + for j in range(jac.shape[1]): + assert np.isclose(jac[i, j], fd_jac[i, j]) class TestIdealGasMoleReactor(TestMoleReactor): reactorClass = ct.IdealGasMoleReactor test_preconditioner_unsupported = None + @pytest.mark.xfail def test_heat_transfer_network(self): # create first reactor gas1 = ct.Solution("h2o2.yaml", "ohmech") gas1.TPX = 600, ct.one_atm, "O2:1.0" r1 = self.reactorClass(gas1) - r1.chemistry_enabled = False # create second reactor gas2 = ct.Solution("h2o2.yaml", "ohmech") gas2.TPX = 300, ct.one_atm, "O2:1.0" r2 = self.reactorClass(gas2) - r2.chemistry_enabled = False # create wall U = 2.0 @@ -1657,14 +1650,11 @@ def test_heat_transfer_network(self): w = ct.Wall(r1, r2, U=U, A=A) net = ct.ReactorNet([r1, r2]) jac = net.jacobian - + fd_jac = net.finite_difference_jacobian # check for values - for i in range(gas1.n_species): - assert np.isclose(jac[0, i + 2 + r1.n_vars], - U * A * gas2.T) - assert np.isclose(jac[r1.n_vars, i + 2], U * A * gas1.T) - if (i + 2 != 3): - assert np.isclose(jac[0, i + 2], U * A * gas1.T, rtol=1e-2) - assert np.isclose(jac[r1.n_vars, i + 2 + r1.n_vars], - U * A * gas2.T, rtol=1e-2) + for i in range(jac.shape[0]): + for j in range(jac.shape[1]): + assert np.isclose(jac[i, j], fd_jac[i, j]) def test_adaptive_precon_integration(self): diff --git a/test/zeroD/test_zeroD.cpp b/test/zeroD/test_zeroD.cpp index 56d9f226095..0885e2818fc 100644 --- a/test/zeroD/test_zeroD.cpp +++ b/test/zeroD/test_zeroD.cpp @@ -306,13 +306,13 @@ TEST(JacobianTests, test_wall_jacobian_build) auto sol1 = newSolution("h2o2.yaml"); sol1->thermo()->setState_TPY(1000.0, OneAtm, " O2:1.0"); IdealGasMoleReactor reactor1; - reactor1.insert(sol1); + reactor1.setSolution(sol1); reactor1.setInitialVolume(1.0); // create second reactor auto sol2 = newSolution("h2o2.yaml"); sol2->thermo()->setState_TPY(900.0, OneAtm, " O2:1.0"); IdealGasConstPressureMoleReactor reactor2; - reactor2.insert(sol2); + reactor2.setSolution(sol2); reactor2.setInitialVolume(1.0); // create the wall Wall w; @@ -332,11 +332,11 @@ TEST(JacobianTests, test_wall_jacobian_build) // build jac for reactor 1 wall only w.buildReactorJacobian(&reactor1, wallJac); wallJacMat.setFromTriplets(wallJac.begin(), wallJac.end()); - // check appropriate values - // double tol = 1e-8; + // check that wall jacobian forms correct value + double v1 = sol1->thermo()->temperature() * w.area() * w.getHeatTransferCoeff(); for (int k = 0; k < wallJacMat.outerSize(); k++) { for (Eigen::SparseMatrix::InnerIterator it(wallJacMat, k); it; ++it) { - EXPECT_DOUBLE_EQ(it.value(), 6000.0); + EXPECT_DOUBLE_EQ(it.value(), v1); EXPECT_EQ(it.row(), 0); // check that it is the first row EXPECT_GE(it.col(), reactor1.speciesOffset()); EXPECT_LT(it.col(), reactor1.neq()); @@ -347,22 +347,16 @@ TEST(JacobianTests, test_wall_jacobian_build) w.buildReactorJacobian(&reactor2, wallJac); wallJacMat.setZero(); wallJacMat.setFromTriplets(wallJac.begin(), wallJac.end()); - // check appropriate values - // double tol = 1e-8; + // check that wall jacobian forms correct value + double v2 = sol2->thermo()->temperature() * w.area() * w.getHeatTransferCoeff(); for (int k = 0; k < wallJacMat.outerSize(); k++) { for (Eigen::SparseMatrix::InnerIterator it(wallJacMat, k); it; ++it) { - EXPECT_DOUBLE_EQ(it.value(), -5400.0); + EXPECT_DOUBLE_EQ(it.value(), -v2); EXPECT_EQ(it.row(), 0); // check that it is the first row EXPECT_GE(it.col(), reactor2.speciesOffset()); EXPECT_LT(it.col(), reactor2.neq()); } } - // check that switch works - wallJac.clear(); - w.jacobianCalculated(); - w.buildNetworkJacobian(wallJac); - EXPECT_EQ(wallJac.size(), 0); - w.jacobianNotCalculated(); // build jac for network terms wallJac.clear(); w.buildNetworkJacobian(wallJac); @@ -387,16 +381,16 @@ TEST(JacobianTests, test_wall_jacobian_build) } } -TEST(JacobianTests, test_flow_jacobian_build) +TEST(JacobianTests, test_flow_jacobian_not_implemented) { // create reservoir reactor auto sol = newSolution("h2o2.yaml"); sol->thermo()->setState_TPY(1000.0, OneAtm, "O2:1.0"); Reservoir res; - res.insert(sol); + res.setSolution(sol); // create reactor IdealGasConstPressureMoleReactor reactor; - reactor.insert(sol); + reactor.setSolution(sol); reactor.setInitialVolume(1.0); // create the flow device MassFlowController mfc; @@ -411,10 +405,6 @@ TEST(JacobianTests, test_flow_jacobian_build) // expect errors from building jacobians EXPECT_THROW(mfc.buildReactorJacobian(&reactor, flowJac), NotImplementedError); // check the jacobian calculated flag and throw/catch errors accordingly - mfc.jacobianCalculated(); - mfc.buildNetworkJacobian(flowJac); - EXPECT_EQ(flowJac.size(), 0); - mfc.jacobianNotCalculated(); EXPECT_THROW(mfc.buildNetworkJacobian(flowJac), NotImplementedError); } From 18e709fa0e8f4a3a0c71c43d33b6738ea85fef21 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Wed, 3 Sep 2025 20:16:56 -0400 Subject: [PATCH 10/13] [Test] Replace deprecated syntax in reactor tests --- test/python/test_reactor.py | 2 +- test/zeroD/test_zeroD.cpp | 56 ++++++++++++++++--------------------- 2 files changed, 25 insertions(+), 33 deletions(-) diff --git a/test/python/test_reactor.py b/test/python/test_reactor.py index aafd1ab279e..0ec99270432 100644 --- a/test/python/test_reactor.py +++ b/test/python/test_reactor.py @@ -3572,7 +3572,7 @@ def replace_eval(self, t, LHS, RHS): net.preconditioner = ct.AdaptivePreconditioner() net.advance_to_steady_state() # reactors should have the same solution because the default is used - self.assertArrayNear(r.get_state(), rstd.get_state()) + assert r.get_state() == approx(rstd.get_state()) class TestSteadySolver: diff --git a/test/zeroD/test_zeroD.cpp b/test/zeroD/test_zeroD.cpp index 0885e2818fc..614f6de565a 100644 --- a/test/zeroD/test_zeroD.cpp +++ b/test/zeroD/test_zeroD.cpp @@ -305,24 +305,20 @@ TEST(JacobianTests, test_wall_jacobian_build) // create first reactor auto sol1 = newSolution("h2o2.yaml"); sol1->thermo()->setState_TPY(1000.0, OneAtm, " O2:1.0"); - IdealGasMoleReactor reactor1; - reactor1.setSolution(sol1); - reactor1.setInitialVolume(1.0); + auto reactor1 = make_shared(sol1); + reactor1->setInitialVolume(1.0); // create second reactor auto sol2 = newSolution("h2o2.yaml"); sol2->thermo()->setState_TPY(900.0, OneAtm, " O2:1.0"); - IdealGasConstPressureMoleReactor reactor2; - reactor2.setSolution(sol2); - reactor2.setInitialVolume(1.0); + auto reactor2 = make_shared(sol2); + reactor2->setInitialVolume(1.0); // create the wall - Wall w; - w.install(reactor1, reactor2); + Wall w(reactor1, reactor2); w.setArea(2.0); w.setHeatTransferCoeff(3.0); // setup reactor network and integrate - ReactorNet network; - network.addReactor(reactor1); - network.addReactor(reactor2); + vector> reactors = {reactor1, reactor2}; + ReactorNet network(reactors); network.initialize(); // create jacobian the size of network Eigen::SparseMatrix wallJacMat; @@ -330,7 +326,7 @@ TEST(JacobianTests, test_wall_jacobian_build) // manually get wall jacobian elements vector> wallJac; // build jac for reactor 1 wall only - w.buildReactorJacobian(&reactor1, wallJac); + w.buildReactorJacobian(reactor1.get(), wallJac); wallJacMat.setFromTriplets(wallJac.begin(), wallJac.end()); // check that wall jacobian forms correct value double v1 = sol1->thermo()->temperature() * w.area() * w.getHeatTransferCoeff(); @@ -338,13 +334,13 @@ TEST(JacobianTests, test_wall_jacobian_build) for (Eigen::SparseMatrix::InnerIterator it(wallJacMat, k); it; ++it) { EXPECT_DOUBLE_EQ(it.value(), v1); EXPECT_EQ(it.row(), 0); // check that it is the first row - EXPECT_GE(it.col(), reactor1.speciesOffset()); - EXPECT_LT(it.col(), reactor1.neq()); + EXPECT_GE(it.col(), reactor1->speciesOffset()); + EXPECT_LT(it.col(), reactor1->neq()); } } // build jac for reactor 2 wall only wallJac.clear(); - w.buildReactorJacobian(&reactor2, wallJac); + w.buildReactorJacobian(reactor2.get(), wallJac); wallJacMat.setZero(); wallJacMat.setFromTriplets(wallJac.begin(), wallJac.end()); // check that wall jacobian forms correct value @@ -353,8 +349,8 @@ TEST(JacobianTests, test_wall_jacobian_build) for (Eigen::SparseMatrix::InnerIterator it(wallJacMat, k); it; ++it) { EXPECT_DOUBLE_EQ(it.value(), -v2); EXPECT_EQ(it.row(), 0); // check that it is the first row - EXPECT_GE(it.col(), reactor2.speciesOffset()); - EXPECT_LT(it.col(), reactor2.neq()); + EXPECT_GE(it.col(), reactor2->speciesOffset()); + EXPECT_LT(it.col(), reactor2->neq()); } } // build jac for network terms @@ -369,13 +365,13 @@ TEST(JacobianTests, test_wall_jacobian_build) if (it.value() < 0) { EXPECT_DOUBLE_EQ(it.value(), -5400.0); EXPECT_EQ(it.row(), 0); // check that it is the first row - EXPECT_GE(it.col(), reactor1.neq() + reactor2.speciesOffset()); - EXPECT_LT(it.col(), reactor1.neq() + reactor2.neq()); + EXPECT_GE(it.col(), reactor1->neq() + reactor2->speciesOffset()); + EXPECT_LT(it.col(), reactor1->neq() + reactor2->neq()); } else { EXPECT_DOUBLE_EQ(it.value(), 6000.0); - EXPECT_EQ(it.row(), reactor1.neq()); // check that it is the first row - EXPECT_GE(it.col(), reactor2.speciesOffset()); - EXPECT_LT(it.col(), reactor1.neq()); + EXPECT_EQ(it.row(), reactor1->neq()); // check that it is the first row + EXPECT_GE(it.col(), reactor2->speciesOffset()); + EXPECT_LT(it.col(), reactor1->neq()); } } } @@ -386,24 +382,20 @@ TEST(JacobianTests, test_flow_jacobian_not_implemented) // create reservoir reactor auto sol = newSolution("h2o2.yaml"); sol->thermo()->setState_TPY(1000.0, OneAtm, "O2:1.0"); - Reservoir res; - res.setSolution(sol); + auto res = make_shared(sol); // create reactor - IdealGasConstPressureMoleReactor reactor; - reactor.setSolution(sol); - reactor.setInitialVolume(1.0); + auto reactor = make_shared(sol); + reactor->setInitialVolume(1.0); // create the flow device - MassFlowController mfc; - mfc.install(res, reactor); + MassFlowController mfc(res, reactor); mfc.setMassFlowCoeff(1.0); // setup reactor network and integrate - ReactorNet network; - network.addReactor(reactor); + ReactorNet network(reactor); network.initialize(); // manually get wall jacobian elements vector> flowJac; // expect errors from building jacobians - EXPECT_THROW(mfc.buildReactorJacobian(&reactor, flowJac), NotImplementedError); + EXPECT_THROW(mfc.buildReactorJacobian(reactor.get(), flowJac), NotImplementedError); // check the jacobian calculated flag and throw/catch errors accordingly EXPECT_THROW(mfc.buildNetworkJacobian(flowJac), NotImplementedError); } From c769cd8f62aa255a3133a9842aee2b2df843dbc9 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Thu, 4 Sep 2025 13:41:19 -0400 Subject: [PATCH 11/13] [Python] Fix derived pointer initialization for Jacobians --- interfaces/cython/cantera/jacobians.pxd | 4 ---- interfaces/cython/cantera/jacobians.pyx | 15 +++++++++------ 2 files changed, 9 insertions(+), 10 deletions(-) diff --git a/interfaces/cython/cantera/jacobians.pxd b/interfaces/cython/cantera/jacobians.pxd index 3a5cd3032b3..3c0593b9452 100644 --- a/interfaces/cython/cantera/jacobians.pxd +++ b/interfaces/cython/cantera/jacobians.pxd @@ -46,20 +46,16 @@ cdef extern from "cantera/numerics/SystemJacobianFactory.h" namespace "Cantera": cdef class SystemJacobian: @staticmethod cdef wrap(shared_ptr[CxxSystemJacobian]) - cdef set_cxx_object(self) cdef shared_ptr[CxxSystemJacobian] _base cdef class EigenSparseJacobian(SystemJacobian): - cdef set_cxx_object(self) cdef CxxEigenSparseJacobian* sparse_jac cdef class EigenSparseDirectJacobian(EigenSparseJacobian): pass cdef class AdaptivePreconditioner(EigenSparseJacobian): - cdef set_cxx_object(self) cdef CxxAdaptivePreconditioner* adaptive cdef class BandedJacobian(SystemJacobian): - cdef set_cxx_object(self) cdef CxxMultiJac* band_jac diff --git a/interfaces/cython/cantera/jacobians.pyx b/interfaces/cython/cantera/jacobians.pyx index 432d2489b09..b6eaa3db74f 100644 --- a/interfaces/cython/cantera/jacobians.pyx +++ b/interfaces/cython/cantera/jacobians.pyx @@ -48,7 +48,7 @@ cdef class SystemJacobian: jac.set_cxx_object() return jac - cdef set_cxx_object(self): + def set_cxx_object(self): pass property side: @@ -67,10 +67,10 @@ cdef class SystemJacobian: """ Get/Set the value of gamma used in the expression P = (I - gamma * J). """ def __get__(self): - return self.pbase.get().gamma() + return self._base.get().gamma() def __set__(self, value): - self.pbase.get().setGamma(value) + self._base.get().setGamma(value) cdef class EigenSparseJacobian(SystemJacobian): @@ -81,7 +81,8 @@ cdef class EigenSparseJacobian(SystemJacobian): _type = "eigen-sparse" - cdef set_cxx_object(self): + def set_cxx_object(self): + super().set_cxx_object() self.sparse_jac = self._base.get() def print_contents(self): @@ -117,7 +118,8 @@ cdef class AdaptivePreconditioner(EigenSparseJacobian): _type = "Adaptive" linear_solver_type = "GMRES" - cdef set_cxx_object(self): + def set_cxx_object(self): + super().set_cxx_object() self.adaptive = self._base.get() property threshold: @@ -188,5 +190,6 @@ cdef class BandedJacobian(SystemJacobian): _type = "banded-direct" linear_solver_type = "direct" - cdef set_cxx_object(self): + def set_cxx_object(self): + super().set_cxx_object() self.band_jac = self._base.get() From 960b8d626567119e4a7ca2a13fe2d2741d234641 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Sun, 7 Sep 2025 13:08:31 -0400 Subject: [PATCH 12/13] [Reactor] Fix species counting in case of user-defined solution components --- src/zeroD/IdealGasConstPressureMoleReactor.cpp | 5 ++++- src/zeroD/IdealGasMoleReactor.cpp | 5 ++++- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/src/zeroD/IdealGasConstPressureMoleReactor.cpp b/src/zeroD/IdealGasConstPressureMoleReactor.cpp index 93f7caf57ef..84354c0b4ef 100644 --- a/src/zeroD/IdealGasConstPressureMoleReactor.cpp +++ b/src/zeroD/IdealGasConstPressureMoleReactor.cpp @@ -125,7 +125,10 @@ void IdealGasConstPressureMoleReactor::buildJacobian( // d (dot(omega)) / d c_j, it is later transformed appropriately. Eigen::SparseMatrix dnk_dnj = m_kin->netProductionRates_ddCi(); // species size that accounts for surface species - size_t ssize = m_nv - m_sidx; + size_t ssize = m_nsp; + for (auto surf : m_surfaces) { + ssize += surf->thermo()->nSpecies(); + } // map derivatives from the surface chemistry jacobian // to the reactor jacobian if (!m_surfaces.empty()) { diff --git a/src/zeroD/IdealGasMoleReactor.cpp b/src/zeroD/IdealGasMoleReactor.cpp index 0b9e0ef78b1..26f7026fb5a 100644 --- a/src/zeroD/IdealGasMoleReactor.cpp +++ b/src/zeroD/IdealGasMoleReactor.cpp @@ -193,7 +193,10 @@ void IdealGasMoleReactor::buildJacobian(vector>& jacVecto // d (dot(omega)) / d c_j, it is later transformed appropriately. Eigen::SparseMatrix dnk_dnj = m_kin->netProductionRates_ddCi(); // species size that accounts for surface species - size_t ssize = m_nv - m_sidx; + size_t ssize = m_nsp; + for (auto surf : m_surfaces) { + ssize += surf->thermo()->nSpecies(); + } // map derivatives from the surface chemistry jacobian // to the reactor jacobian if (!m_surfaces.empty()) { From e899b210edd14267031e28dd7110369fc7344bb1 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Fri, 12 Dec 2025 12:43:27 -0500 Subject: [PATCH 13/13] [Reactor] Fix duplicate calculation of connector Jacobians --- src/zeroD/ReactorNet.cpp | 14 +++++--------- 1 file changed, 5 insertions(+), 9 deletions(-) diff --git a/src/zeroD/ReactorNet.cpp b/src/zeroD/ReactorNet.cpp index efea6ff1ee6..61539238aa9 100644 --- a/src/zeroD/ReactorNet.cpp +++ b/src/zeroD/ReactorNet.cpp @@ -872,15 +872,11 @@ void ReactorNet::buildJacobian(vector>& jacVector) // loop through all connections and then set them found so calculations are not // repeated - for (auto r : m_reactors) { - // walls - for (const auto& wall : m_walls) { - wall->buildNetworkJacobian(jacVector); - } - // flow devices - for (const auto& flow_device : m_flow_devices) { - flow_device->buildNetworkJacobian(jacVector); - } + for (const auto& wall : m_walls) { + wall->buildNetworkJacobian(jacVector); + } + for (const auto& flow_device : m_flow_devices) { + flow_device->buildNetworkJacobian(jacVector); } }