diff --git a/include/cantera/oneD/Domain1D.h b/include/cantera/oneD/Domain1D.h index c3c5a5bd4..e67fb8b55 100644 --- a/include/cantera/oneD/Domain1D.h +++ b/include/cantera/oneD/Domain1D.h @@ -20,6 +20,7 @@ class Kinetics; class Transport; class Solution; class SolutionArray; +class SystemJacobian; /** * Base class for one-dimensional domains. @@ -259,6 +260,43 @@ public: return m_min[n]; } + //! Set the method used to evaluate this domain's Jacobian columns. + //! @param mode Either `"finite-difference"` (default) or `"analytic"`. + //! In `"analytic"` mode, derived classes that implement analytic + //! Jacobian elements (see hasAnalyticJacobian()) compute them directly + //! instead of by finite-differencing the residual. + //! @since New in %Cantera 4.0. + void setJacobianMode(const string& mode) { + if (mode != "finite-difference" && mode != "analytic") { + throw CanteraError("Domain1D::setJacobianMode", + "Unknown Jacobian mode '{}'", mode); + } + m_jacobianMode = mode; + } + + //! Get the method used to evaluate this domain's Jacobian columns. + //! @since New in %Cantera 4.0. + const string& jacobianMode() const { + return m_jacobianMode; + } + + //! Returns `true` if this domain computes the Jacobian column for + //! component `n` at (domain-local) grid point `j` analytically, in which + //! case the finite-difference evaluation of that column is skipped and the + //! domain must provide all entries of the column in evalJacobianAnalytic(). + //! @since New in %Cantera 4.0. + virtual bool hasAnalyticJacobian(size_t j, size_t n) const { + return false; + } + + //! Add this domain's analytic Jacobian entries (for columns claimed by + //! hasAnalyticJacobian()) to `jac` via SystemJacobian::setValue(), using + //! global row/column indices (domain-local index + loc()). Entries must be + //! the steady-state Jacobian; transient diagonal terms are handled + //! separately. Base class implementation does nothing. + //! @since New in %Cantera 4.0. + virtual void evalJacobianAnalytic(span x, SystemJacobian& jac) {} + /** * Set grid refinement criteria. @see Refiner::setCriteria. * @since New in %Cantera 3.2 @@ -715,6 +753,7 @@ protected: vector m_name; //!< Names of solution components int m_bw = -1; //!< See bandwidth() bool m_force_full_update = false; //!< see forceFullUpdate() + string m_jacobianMode = "finite-difference"; //!< see setJacobianMode() //! Composite thermo/kinetics/transport handler shared_ptr m_solution; diff --git a/interfaces/cython/cantera/_onedim.pxd b/interfaces/cython/cantera/_onedim.pxd index 279454ffa..044ff75ca 100644 --- a/interfaces/cython/cantera/_onedim.pxd +++ b/interfaces/cython/cantera/_onedim.pxd @@ -55,6 +55,8 @@ cdef extern from "cantera/oneD/Domain1D.h": string domainType() shared_ptr[CxxSolutionArray] toArray(cbool) except +translate_exception void fromArray(shared_ptr[CxxSolutionArray]) except +translate_exception + void setJacobianMode(string&) except +translate_exception + string jacobianMode() cdef extern from "cantera/oneD/Boundary1D.h": diff --git a/interfaces/cython/cantera/_onedim.pyx b/interfaces/cython/cantera/_onedim.pyx index c0a72f7ff..095222aae 100644 --- a/interfaces/cython/cantera/_onedim.pyx +++ b/interfaces/cython/cantera/_onedim.pyx @@ -434,6 +434,20 @@ cdef class Domain1D: else: return self.domain.transient_atol(self.component_index(component)) + property jacobian_mode: + """ + Method used to evaluate this domain's Jacobian columns: + ``'finite-difference'`` (default) or ``'analytic'``. In ``'analytic'`` + mode, domains that support it compute some Jacobian columns + analytically instead of by finite differences. + + .. versionadded:: 4.0 + """ + def __get__(self): + return pystr(self.domain.jacobianMode()) + def __set__(self, mode): + self.domain.setJacobianMode(stringify(mode)) + property name: """ The name / id of this domain """ def __get__(self): diff --git a/src/oneD/OneDim.cpp b/src/oneD/OneDim.cpp index 04638339e..131beb344 100644 --- a/src/oneD/OneDim.cpp +++ b/src/oneD/OneDim.cpp @@ -286,7 +286,14 @@ void OneDim::evalJacobian(span x0) size_t ipt = 0; for (size_t j = 0; j < points(); j++) { size_t nv = nVars(j); + Domain1D* dom = pointDomain(ipt); + size_t jLocal = j - dom->firstPoint(); for (size_t n = 0; n < nv; n++) { + if (dom->hasAnalyticJacobian(jLocal, n)) { + // skip FD perturbation; analytic fill handled below + ipt++; + continue; + } // perturb x(n); preserve sign(x(n)) double xsave = x0[ipt]; double dx = fabs(xsave) * m_jacobianRelPerturb + m_jacobianAbsPerturb; @@ -317,6 +324,11 @@ void OneDim::evalJacobian(span x0) } } + // analytic contributions for claimed columns + for (auto& dom : m_dom) { + dom->evalJacobianAnalytic(x0, *m_jac); + } + m_jac->updateElapsed(double(clock() - t0) / CLOCKS_PER_SEC); m_jac->incrementEvals(); m_jac->setAge(0); diff --git a/test/python/test_onedim.py b/test/python/test_onedim.py index c3c81477f..8a135fb69 100644 --- a/test/python/test_onedim.py +++ b/test/python/test_onedim.py @@ -2361,3 +2361,24 @@ class TestEvalJacobian: # update. Differences are bounded by the neglected dD/dY terms (<1%). scale = np.abs(fd_col).max() assert np.abs(jac_col - fd_col).max() < 2e-2 * scale + + +class TestJacobianMode: + def test_mode_roundtrip(self): + gas = ct.Solution("h2o2.yaml") + flame = ct.FreeFlow(gas) + assert flame.jacobian_mode == "finite-difference" + flame.jacobian_mode = "analytic" + assert flame.jacobian_mode == "analytic" + with pytest.raises(ct.CanteraError, match="Unknown Jacobian mode"): + flame.jacobian_mode = "automagic" + + def test_analytic_mode_unclaimed_is_identical(self): + # Until Flow1D implements claims, analytic mode must not change anything. + # NOTE: once Flow1D claims Y-columns (a later task) this becomes a + # closeness test; it is replaced there by TestAnalyticVsFD. + gas, sim = make_flame() + J_fd = get_jacobian(sim) + sim.flame.jacobian_mode = "analytic" + J_an = get_jacobian(sim) + assert np.array_equal(J_fd, J_an)