From bda0e9c113a239267d4e52a87d8e1692172d4f79 Mon Sep 17 00:00:00 2001 From: Anthony Walker Date: Tue, 10 Oct 2023 10:47:23 -0700 Subject: [PATCH 1/9] 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/preconditioners.pxd | 2 + interfaces/cython/cantera/preconditioners.pyx | 10 +++ 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 | 86 ++++++++++++++++++- 17 files changed, 250 insertions(+), 44 deletions(-) diff --git a/include/cantera/base/Delegator.h b/include/cantera/base/Delegator.h index ea02f18a12..032ae58bab 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 9a0301cec6..3d661d5112 100644 --- a/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h +++ b/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h @@ -43,7 +43,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 14e3527f9b..057192151e 100644 --- a/include/cantera/zeroD/IdealGasMoleReactor.h +++ b/include/cantera/zeroD/IdealGasMoleReactor.h @@ -38,12 +38,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 526668745c..f16e6465b0 100644 --- a/include/cantera/zeroD/Reactor.h +++ b/include/cantera/zeroD/Reactor.h @@ -187,7 +187,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. @@ -196,9 +196,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 1a30e79f71..1e3c3ce1c1 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 934a36ed0b..2607969e93 100644 --- a/interfaces/cython/cantera/_utils.pxd +++ b/interfaces/cython/cantera/_utils.pxd @@ -8,6 +8,7 @@ from libcpp.unordered_map cimport unordered_map 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" @@ -115,3 +116,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 1eb3c11c2b..58244c4351 100644 --- a/interfaces/cython/cantera/_utils.pyx +++ b/interfaces/cython/cantera/_utils.pyx @@ -5,6 +5,7 @@ import sys import os import warnings from cpython.ref cimport PyObject +from cython.operator cimport dereference import numbers import importlib.metadata from collections import namedtuple @@ -525,6 +526,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 721ddd736a..09843a3faa 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 fe2dc16cd5..a38e5f0144 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/preconditioners.pxd b/interfaces/cython/cantera/preconditioners.pxd index ebe80f6a02..c8e3414c72 100644 --- a/interfaces/cython/cantera/preconditioners.pxd +++ b/interfaces/cython/cantera/preconditioners.pxd @@ -12,6 +12,8 @@ cdef extern from "cantera/numerics/PreconditionerBase.h" namespace "Cantera": CxxPreconditionerBase() string preconditionerSide() void setPreconditionerSide(string) except +translate_exception + double gamma() + void setGamma(double) except +translate_exception cdef extern from "cantera/numerics/AdaptivePreconditioner.h" namespace "Cantera": cdef cppclass CxxAdaptivePreconditioner "Cantera::AdaptivePreconditioner" \ diff --git a/interfaces/cython/cantera/preconditioners.pyx b/interfaces/cython/cantera/preconditioners.pyx index ebe5bcc300..d489f2c29b 100644 --- a/interfaces/cython/cantera/preconditioners.pyx +++ b/interfaces/cython/cantera/preconditioners.pyx @@ -26,6 +26,16 @@ cdef class PreconditionerBase: def __set__(self, side): self.pbase.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 AdaptivePreconditioner(PreconditionerBase): precon_type = "Adaptive" precon_linear_solver_type = "GMRES" diff --git a/interfaces/cython/cantera/reactor.pxd b/interfaces/cython/cantera/reactor.pxd index c06a9e0b8c..765839f060 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 .preconditioners cimport * @@ -59,6 +60,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 addSurface(CxxReactorSurface*) void setAdvanceLimit(string&, double) except +translate_exception void addSensitivityReaction(size_t) except +translate_exception diff --git a/interfaces/cython/cantera/reactor.pyx b/interfaces/cython/cantera/reactor.pyx index 39a05f19d4..5354508a72 100644 --- a/interfaces/cython/cantera/reactor.pyx +++ b/interfaces/cython/cantera/reactor.pyx @@ -654,7 +654,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 757a6cb95b..0f512ef4f6 100644 --- a/src/zeroD/IdealGasConstPressureMoleReactor.cpp +++ b/src/zeroD/IdealGasConstPressureMoleReactor.cpp @@ -119,14 +119,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(); @@ -172,7 +171,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()); } } @@ -201,7 +200,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 @@ -226,14 +225,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 5586a7cc80..40f566621d 100644 --- a/src/zeroD/IdealGasMoleReactor.cpp +++ b/src/zeroD/IdealGasMoleReactor.cpp @@ -153,14 +153,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(); @@ -185,7 +183,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()); } } @@ -213,7 +211,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 @@ -242,14 +240,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 3d1c4f2635..7f197e56ad 100644 --- a/src/zeroD/ReactorNet.cpp +++ b/src/zeroD/ReactorNet.cpp @@ -578,18 +578,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 15de152b7a..43b67e6182 100644 --- a/test/python/test_reactor.py +++ b/test/python/test_reactor.py @@ -3133,6 +3133,86 @@ def deltaC(): vdot=r2.expansion_rate) # Regression test values - assert r1.thermo.P == approx(151561.15, rel=1e-6) - assert r1.thermo["H2"].Y[0] == approx(0.13765976, rel=1e-6) - assert r2.thermo["O2"].Y[0] == approx(0.94617029, rel=1e-6) + self.assertNear(r1.thermo.P, 151561.15, rtol=1e-6) + self.assertNear(r1.thermo["H2"].Y[0], 0.13765976, rtol=1e-6) + self.assertNear(r2.thermo["O2"].Y[0], 0.94617029, rtol=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.thermo.P - self.neighbor.thermo.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 From c1cb104bdfd4005935b221438e9486ad54def1ce Mon Sep 17 00:00:00 2001 From: Anthony Walker Date: Wed, 11 Oct 2023 20:45:42 -0700 Subject: [PATCH 2/9] Complete finite difference jacobian --- include/cantera/zeroD/ReactorNet.h | 38 ++++++++++++++ src/zeroD/ReactorNet.cpp | 84 +++++++++++++++++++++++++----- test/python/test_reactor.py | 11 ++++ 3 files changed, 120 insertions(+), 13 deletions(-) diff --git a/include/cantera/zeroD/ReactorNet.h b/include/cantera/zeroD/ReactorNet.h index f1ffe5b249..181283c00b 100644 --- a/include/cantera/zeroD/ReactorNet.h +++ b/include/cantera/zeroD/ReactorNet.h @@ -299,7 +299,45 @@ class ReactorNet : public FuncEval //! @param settings the settings map propagated to all reactors and kinetics objects virtual void setDerivativeSettings(AnyMap& settings); + //! 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: + //! 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 7f197e56ad..80e0129498 100644 --- a/src/zeroD/ReactorNet.cpp +++ b/src/zeroD/ReactorNet.cpp @@ -584,20 +584,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->setup(); @@ -626,4 +616,72 @@ 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; +} + + } diff --git a/test/python/test_reactor.py b/test/python/test_reactor.py index 43b67e6182..2795ada8df 100644 --- a/test/python/test_reactor.py +++ b/test/python/test_reactor.py @@ -1437,6 +1437,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 net1 = ct.ReactorNet() From 1c0dbee2dabd4ef27b848b320c5c67196a683704 Mon Sep 17 00:00:00 2001 From: Anthony Walker Date: Tue, 17 Oct 2023 10:33:12 -0700 Subject: [PATCH 3/9] Build jacobian structure defined for walls and flow devices --- include/cantera/zeroD/FlowDevice.h | 40 ++++++++++++ .../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, 506 insertions(+), 10 deletions(-) diff --git a/include/cantera/zeroD/FlowDevice.h b/include/cantera/zeroD/FlowDevice.h index 00e622b2ba..16451366e1 100644 --- a/include/cantera/zeroD/FlowDevice.h +++ b/include/cantera/zeroD/FlowDevice.h @@ -9,6 +9,7 @@ #include "cantera/base/ct_defs.h" #include "cantera/base/global.h" #include "cantera/base/ctexceptions.h" +#include "cantera/numerics/eigen_sparse.h" namespace Cantera { @@ -132,9 +133,48 @@ class FlowDevice 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; diff --git a/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h b/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h index 3d661d5112..a77a491913 100644 --- a/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h +++ b/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h @@ -47,6 +47,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; }; + protected: void setThermo(ThermoPhase& thermo) override; diff --git a/include/cantera/zeroD/IdealGasMoleReactor.h b/include/cantera/zeroD/IdealGasMoleReactor.h index 057192151e..dca17d00d4 100644 --- a/include/cantera/zeroD/IdealGasMoleReactor.h +++ b/include/cantera/zeroD/IdealGasMoleReactor.h @@ -43,6 +43,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: void setThermo(ThermoPhase& thermo) override; diff --git a/include/cantera/zeroD/MoleReactor.h b/include/cantera/zeroD/MoleReactor.h index 8f4f423088..4b2c51760c 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 { @@ -38,6 +39,8 @@ class MoleReactor : public Reactor string componentName(size_t k) 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 @@ -60,6 +63,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 f16e6465b0..1a62cc148f 100644 --- a/include/cantera/zeroD/Reactor.h +++ b/include/cantera/zeroD/Reactor.h @@ -219,7 +219,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. //! @@ -314,6 +334,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 c048471276..89a60f5883 100644 --- a/include/cantera/zeroD/ReactorBase.h +++ b/include/cantera/zeroD/ReactorBase.h @@ -266,6 +266,58 @@ class ReactorBase //! Set the ReactorNet that this reactor belongs to. void setNetwork(ReactorNet* net); + /*! 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: //! Specify the mixture contained in the reactor. Note that a pointer to //! this substance is stored, and as the integration proceeds, the state of @@ -282,6 +334,12 @@ class ReactorBase //! 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 = 1.0; //!< Current volume of the reactor [m^3] double m_enthalpy = 0.0; //!< Current specific enthalpy of the reactor [J/kg] diff --git a/include/cantera/zeroD/ReactorNet.h b/include/cantera/zeroD/ReactorNet.h index 181283c00b..7e817c7efb 100644 --- a/include/cantera/zeroD/ReactorNet.h +++ b/include/cantera/zeroD/ReactorNet.h @@ -236,6 +236,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'`. @@ -400,6 +414,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 14b9dbd034..2368c1adee 100644 --- a/include/cantera/zeroD/Wall.h +++ b/include/cantera/zeroD/Wall.h @@ -8,6 +8,7 @@ #include "cantera/base/ctexceptions.h" #include "cantera/zeroD/ReactorBase.h" +#include "cantera/numerics/eigen_sparse.h" namespace Cantera { @@ -104,6 +105,46 @@ class WallBase 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: string m_name; //!< Wall name. bool m_defaultNameSet = false; //!< `true` if default name has been previously set. @@ -115,6 +156,10 @@ class WallBase 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. @@ -219,6 +264,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 765839f060..c25443a8da 100644 --- a/interfaces/cython/cantera/reactor.pxd +++ b/interfaces/cython/cantera/reactor.pxd @@ -209,6 +209,8 @@ cdef extern from "cantera/zerodim.h" namespace "Cantera": void setPreconditioner(shared_ptr[CxxPreconditionerBase] 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 5354508a72..9bfae66e5c 100644 --- a/interfaces/cython/cantera/reactor.pyx +++ b/interfaces/cython/cantera/reactor.pyx @@ -2023,6 +2023,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 0f512ef4f6..f02c99f249 100644 --- a/src/zeroD/IdealGasConstPressureMoleReactor.cpp +++ b/src/zeroD/IdealGasConstPressureMoleReactor.cpp @@ -203,6 +203,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); @@ -228,7 +229,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 @@ -266,4 +273,23 @@ string IdealGasConstPressureMoleReactor::componentName(size_t k) { "Index is out of bounds."); } +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 40f566621d..98728778b9 100644 --- a/src/zeroD/IdealGasMoleReactor.cpp +++ b/src/zeroD/IdealGasMoleReactor.cpp @@ -214,6 +214,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); @@ -243,7 +244,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 8777010d94..7583c2cdd5 100644 --- a/src/zeroD/Reactor.cpp +++ b/src/zeroD/Reactor.cpp @@ -41,6 +41,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::setKinetics(Kinetics& kin) @@ -596,4 +607,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 80e0129498..28d8ea96c7 100644 --- a/src/zeroD/ReactorNet.cpp +++ b/src/zeroD/ReactorNet.cpp @@ -535,6 +535,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 @@ -622,6 +632,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 @@ -639,6 +669,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 2d1d61a819..3bf2cbc6f9 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 { @@ -88,4 +90,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->contents().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->contents().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->contents().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 2795ada8df..381a08584d 100644 --- a/test/python/test_reactor.py +++ b/test/python/test_reactor.py @@ -1432,21 +1432,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 446e3f1c069257d7e86641012bd1d7724d33e8ec Mon Sep 17 00:00:00 2001 From: Anthony Walker Date: Tue, 17 Oct 2023 23:03:54 -0700 Subject: [PATCH 4/9] 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 | 8 +++++--- 3 files changed, 11 insertions(+), 3 deletions(-) diff --git a/interfaces/cython/cantera/_utils.pxd b/interfaces/cython/cantera/_utils.pxd index 2607969e93..c7cd1d2fcf 100644 --- a/interfaces/cython/cantera/_utils.pxd +++ b/interfaces/cython/cantera/_utils.pxd @@ -119,3 +119,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 58244c4351..9f66f49480 100644 --- a/interfaces/cython/cantera/_utils.pyx +++ b/interfaces/cython/cantera/_utils.pyx @@ -545,6 +545,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 a38e5f0144..46aa03f9af 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,10 +247,12 @@ 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] funcInfo.setExceptionType(exc_type) From 00e517449ae66644b0625cc5bbd646e3d51ef50c Mon Sep 17 00:00:00 2001 From: Anthony Walker Date: Tue, 24 Oct 2023 11:52:47 -0400 Subject: [PATCH 5/9] 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 aa34f3363e..941bdaf958 100644 --- a/test/zeroD/test_zeroD.cpp +++ b/test/zeroD/test_zeroD.cpp @@ -208,6 +208,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 ec29e54d9ba0a589267f537ef1520be9d2562223 Mon Sep 17 00:00:00 2001 From: Anthony Walker Date: Tue, 24 Oct 2023 12:11:50 -0400 Subject: [PATCH 6/9] 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 16451366e1..45c6319d77 100644 --- a/include/cantera/zeroD/FlowDevice.h +++ b/include/cantera/zeroD/FlowDevice.h @@ -133,40 +133,47 @@ class FlowDevice 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 1a62cc148f..3acb946cef 100644 --- a/include/cantera/zeroD/Reactor.h +++ b/include/cantera/zeroD/Reactor.h @@ -206,8 +206,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. @@ -231,7 +230,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 89a60f5883..97051d2a8f 100644 --- a/include/cantera/zeroD/ReactorBase.h +++ b/include/cantera/zeroD/ReactorBase.h @@ -266,48 +266,50 @@ class ReactorBase //! Set the ReactorNet that this reactor belongs to. void setNetwork(ReactorNet* net); - /*! 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 7e817c7efb..013b237fd2 100644 --- a/include/cantera/zeroD/ReactorNet.h +++ b/include/cantera/zeroD/ReactorNet.h @@ -341,8 +341,7 @@ class ReactorNet : public FuncEval protected: //! 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 2368c1adee..a3937df1ca 100644 --- a/include/cantera/zeroD/Wall.h +++ b/include/cantera/zeroD/Wall.h @@ -105,44 +105,48 @@ class WallBase 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 28d8ea96c7..5ed00469d2 100644 --- a/src/zeroD/ReactorNet.cpp +++ b/src/zeroD/ReactorNet.cpp @@ -593,10 +593,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 e1e3cab1b6f20cf579ff74152696aed38f3a71a6 Mon Sep 17 00:00:00 2001 From: Anthony Walker Date: Tue, 24 Oct 2023 15:14:56 -0400 Subject: [PATCH 7/9] 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 -- include/cantera/zeroD/ReactorBase.h | 26 ------- 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 ++++++++++ 11 files changed, 111 insertions(+), 77 deletions(-) diff --git a/include/cantera/zeroD/FlowDevice.h b/include/cantera/zeroD/FlowDevice.h index 45c6319d77..01404d48ea 100644 --- a/include/cantera/zeroD/FlowDevice.h +++ b/include/cantera/zeroD/FlowDevice.h @@ -158,7 +158,9 @@ class FlowDevice //! @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 a77a491913..d54060060c 100644 --- a/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h +++ b/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h @@ -47,12 +47,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 dca17d00d4..fe245d8407 100644 --- a/include/cantera/zeroD/IdealGasMoleReactor.h +++ b/include/cantera/zeroD/IdealGasMoleReactor.h @@ -43,12 +43,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/include/cantera/zeroD/ReactorBase.h b/include/cantera/zeroD/ReactorBase.h index 97051d2a8f..64c4b4da3b 100644 --- a/include/cantera/zeroD/ReactorBase.h +++ b/include/cantera/zeroD/ReactorBase.h @@ -266,32 +266,6 @@ class ReactorBase //! Set the ReactorNet that this reactor belongs to. void setNetwork(ReactorNet* net); - //! 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 diff --git a/interfaces/cython/cantera/_utils.pxd b/interfaces/cython/cantera/_utils.pxd index c7cd1d2fcf..02f152de62 100644 --- a/interfaces/cython/cantera/_utils.pxd +++ b/interfaces/cython/cantera/_utils.pxd @@ -117,7 +117,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 9f66f49480..a8e165b99f 100644 --- a/interfaces/cython/cantera/_utils.pyx +++ b/interfaces/cython/cantera/_utils.pyx @@ -526,41 +526,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 46aa03f9af..9ecaf95ebe 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 7583c2cdd5..e93cc5e035 100644 --- a/src/zeroD/Reactor.cpp +++ b/src/zeroD/Reactor.cpp @@ -623,7 +623,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 5ed00469d2..edd81f2a5a 100644 --- a/src/zeroD/ReactorNet.cpp +++ b/src/zeroD/ReactorNet.cpp @@ -681,7 +681,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 381a08584d..034b4f47ca 100644 --- a/test/python/test_reactor.py +++ b/test/python/test_reactor.py @@ -183,6 +183,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.thermo.X[k1H2] > 0.3 or self.r2.thermo.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.gas1.species_names + self.gas2.species_names) + r1_species_start = self.r1.component_index(self.gas1.species_name(0)) + r2_species_start = self.r2.component_index(self.gas2.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() @@ -1424,7 +1479,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() @@ -1432,6 +1487,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 941bdaf958..be0658c4ac 100644 --- a/test/zeroD/test_zeroD.cpp +++ b/test/zeroD/test_zeroD.cpp @@ -265,6 +265,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); @@ -289,6 +295,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 7794f7eac11adedb3de09f68e4757be6c2aad777 Mon Sep 17 00:00:00 2001 From: Anthony Walker Date: Mon, 1 Jan 2024 17:32:43 -0500 Subject: [PATCH 8/9] 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 1e3c3ce1c1..50a07dcc8f 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 c25443a8da..6523dd96b3 100644 --- a/interfaces/cython/cantera/reactor.pxd +++ b/interfaces/cython/cantera/reactor.pxd @@ -222,6 +222,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 9bfae66e5c..7d9acfbadc 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 @@ -714,6 +716,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 034b4f47ca..2e9052e38e 100644 --- a/test/python/test_reactor.py +++ b/test/python/test_reactor.py @@ -3340,3 +3340,26 @@ def replace_build_jacobian(self, jac_vector): # test that jacobian wall element is hard coded value 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()) From 3b892b644abcac93c29e6036f5002380b4a0793a Mon Sep 17 00:00:00 2001 From: Anthony Walker Date: Sun, 3 Mar 2024 10:00:32 -0500 Subject: [PATCH 9/9] Working on PR updates, have to fix a few more bugs in regards to equations, handling volumes, etc. --- doc/sphinx/develop/index.md | 23 +++--- include/cantera/zeroD/FlowDevice.h | 26 ++----- .../zeroD/IdealGasConstPressureMoleReactor.h | 4 +- include/cantera/zeroD/IdealGasMoleReactor.h | 4 +- include/cantera/zeroD/Reactor.h | 23 +----- include/cantera/zeroD/ReactorBase.h | 42 ++++++----- 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, 162 insertions(+), 261 deletions(-) diff --git a/doc/sphinx/develop/index.md b/doc/sphinx/develop/index.md index e38e040c55..d6fa04700f 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 @@ -47,13 +48,13 @@ reactor-integration ## Adding New Features to Cantera -- [](CONTRIBUTING) -- [](style-guidelines) -- [](vscode-tips) -- [](writing-tests) -- [](running-tests) -- [](writing-examples) -- [](doc-formatting) +- [](CONTRIBUTING) +- [](style-guidelines) +- [](vscode-tips) +- [](writing-tests) +- [](running-tests) +- [](writing-examples) +- [](doc-formatting) ```{toctree} :caption: Adding New Features to Cantera diff --git a/include/cantera/zeroD/FlowDevice.h b/include/cantera/zeroD/FlowDevice.h index 01404d48ea..a12b6cb050 100644 --- a/include/cantera/zeroD/FlowDevice.h +++ b/include/cantera/zeroD/FlowDevice.h @@ -141,9 +141,10 @@ class FlowDevice //! @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"); } @@ -155,29 +156,12 @@ class FlowDevice //! @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 d54060060c..db186fc7ba 100644 --- a/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h +++ b/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h @@ -47,9 +47,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 fe245d8407..c5de108132 100644 --- a/include/cantera/zeroD/IdealGasMoleReactor.h +++ b/include/cantera/zeroD/IdealGasMoleReactor.h @@ -43,9 +43,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 3acb946cef..7185395259 100644 --- a/include/cantera/zeroD/Reactor.h +++ b/include/cantera/zeroD/Reactor.h @@ -70,12 +70,6 @@ class Reactor : public ReactorBase m_chem = cflag; } - //! Returns `true` if changes in the reactor composition due to chemical reactions - //! are enabled. - bool chemistryEnabled() const { - return m_chem; - } - void setEnergy(int eflag=1) override { if (eflag > 0) { m_energy = true; @@ -84,11 +78,6 @@ class Reactor : public ReactorBase } } - //! Returns `true` if solution of the energy equation is enabled. - bool energyEnabled() const { - return m_energy; - } - //! Number of equations (state variables) for this reactor size_t neq() { if (!m_nv) { @@ -195,15 +184,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 @@ -321,8 +302,6 @@ class Reactor : public ReactorBase vector m_wdot; //!< Species net molar production rates vector m_uk; //!< Species molar internal energies - bool m_chem = false; - bool m_energy = true; size_t m_nv = 0; size_t m_nv_surf; //!!< Number of variables associated with reactor surfaces diff --git a/include/cantera/zeroD/ReactorBase.h b/include/cantera/zeroD/ReactorBase.h index 64c4b4da3b..4956264b1f 100644 --- a/include/cantera/zeroD/ReactorBase.h +++ b/include/cantera/zeroD/ReactorBase.h @@ -266,26 +266,15 @@ class ReactorBase //! Set the ReactorNet that this reactor belongs to. void setNetwork(ReactorNet* net); - //! 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 @@ -294,6 +283,17 @@ class ReactorBase //! Return the offset between species and state variables virtual size_t speciesOffset() const { return m_sidx; }; + //! Returns `true` if solution of the energy equation is enabled. + virtual bool energyEnabled() const { + return m_energy; + } + + //! Returns `true` if changes in the reactor composition due to chemical reactions + //! are enabled. + bool chemistryEnabled() const { + return m_chem; + } + protected: //! Specify the mixture contained in the reactor. Note that a pointer to //! this substance is stored, and as the integration proceeds, the state of @@ -338,6 +338,12 @@ class ReactorBase //! Composite thermo/kinetics/transport handler shared_ptr m_solution; + + //! A bool that enables the energy equation + bool m_energy = true; + + //! A bool that enables the chemical kinetics equations + bool m_chem = false; }; } diff --git a/include/cantera/zeroD/ReactorDelegator.h b/include/cantera/zeroD/ReactorDelegator.h index 50a07dcc8f..ce826f7201 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 013b237fd2..9bbc387c23 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" @@ -239,16 +241,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 @@ -417,6 +410,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 a3937df1ca..7d5840f406 100644 --- a/include/cantera/zeroD/Wall.h +++ b/include/cantera/zeroD/Wall.h @@ -108,47 +108,27 @@ class WallBase //! 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: string m_name; //!< Wall name. bool m_defaultNameSet = false; //!< `true` if default name has been previously set. @@ -160,10 +140,6 @@ class WallBase 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. @@ -268,9 +244,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 9ecaf95ebe..075f6e762e 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 f02c99f249..b80cb01fe5 100644 --- a/src/zeroD/IdealGasConstPressureMoleReactor.cpp +++ b/src/zeroD/IdealGasConstPressureMoleReactor.cpp @@ -273,23 +273,11 @@ string IdealGasConstPressureMoleReactor::componentName(size_t k) { "Index is out of bounds."); } -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 98728778b9..452129d755 100644 --- a/src/zeroD/IdealGasMoleReactor.cpp +++ b/src/zeroD/IdealGasMoleReactor.cpp @@ -244,32 +244,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 e93cc5e035..9ecaa56c0f 100644 --- a/src/zeroD/Reactor.cpp +++ b/src/zeroD/Reactor.cpp @@ -607,7 +607,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++) { @@ -629,4 +629,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 edd81f2a5a..e53c0a40d5 100644 --- a/src/zeroD/ReactorNet.cpp +++ b/src/zeroD/ReactorNet.cpp @@ -109,7 +109,26 @@ void ReactorNet::initialize() "FlowReactors must be used alone."); } } - + // 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); @@ -529,6 +548,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 @@ -629,29 +658,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 @@ -674,24 +683,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); } } } @@ -703,7 +700,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 3bf2cbc6f9..5cbc2200c5 100644 --- a/src/zeroD/Wall.cpp +++ b/src/zeroD/Wall.cpp @@ -91,7 +91,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; @@ -99,12 +100,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); } } @@ -115,38 +117,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->contents().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->contents().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 2e9052e38e..ba1e333c12 100644 --- a/test/python/test_reactor.py +++ b/test/python/test_reactor.py @@ -1499,51 +1499,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 @@ -1551,14 +1544,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 be0658c4ac..56c67264a6 100644 --- a/test/zeroD/test_zeroD.cpp +++ b/test/zeroD/test_zeroD.cpp @@ -214,13 +214,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; @@ -240,11 +240,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()); @@ -255,22 +255,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); @@ -295,16 +289,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; @@ -319,10 +313,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); }