[1D] Add infrastructure for analytic Jacobian columns

Adds a per-domain `jacobian_mode` flag ("finite-difference" or
"analytic"), two virtual hooks (`hasAnalyticJacobian(j,n)` and
`evalJacobianAnalytic(x, jac)`), and the column-skipping logic in
`OneDim::evalJacobian` that calls `evalJacobianAnalytic` after the FD
loop. No domain claims any columns yet, so this is a behavioral no-op.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This commit is contained in:
Ray Speth
2026-07-03 09:54:37 -05:00
committed by Ingmar Schoegl
co-authored by Claude Opus 4.8
parent 5b503cd7eb
commit 431d7eaf08
5 changed files with 88 additions and 0 deletions
+39
View File
@@ -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<const double> x, SystemJacobian& jac) {}
/**
* Set grid refinement criteria. @see Refiner::setCriteria.
* @since New in %Cantera 3.2
@@ -715,6 +753,7 @@ protected:
vector<string> 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<Solution> m_solution;
+2
View File
@@ -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":
+14
View File
@@ -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):
+12
View File
@@ -286,7 +286,14 @@ void OneDim::evalJacobian(span<const double> 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<const double> 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);
+21
View File
@@ -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)