High: [Crystallizers] Concentration Jacobian diagonal has wrong sign for crystallization sink

[bernalde]

Maintainer triage (June 29, 2026)

• Status: Verified bug against the current master source in this fork.
• Severity: High.
• Area labels: area:crystallizers.
• Tackle batch: P4 - high-impact numerical/model correctness.

The source review confirmed that the cited code path and failure mode are still present. The original audit details are preserved below for reproduction notes and suggested fixes.

---

Original audit details

Part of the model-correctness audit — see #17.

Severity: high | Category: sign-convention | Location: PharmaPy/Crystallizers.py lines 1401-1410

What's wrong

In jac_states, the concentration RHS is dcomp_i/dt = -transf/vol * (kron_jtg_i - conc_i/rho_l) (material_balances line 1513, with transf == tr from line 1401). Differentiating w.r.t. conc_j gives two terms: dfconc_dconc_ij = -(1/vol)[ (d tr/d conc_j)(kron_i - conc_i/rho_l) - trdelta_ij/rho_l ]. The second (diagonal) term must enter the bracket with a MINUS sign because d(kron_i - conc_i/rho_l)/d conc_j = -delta_ij/rho_l. The code instead builds second_conc = tr/rho_lnp.eye(...) and adds it: dfconc_dconc = -1/vol_liq*(dtr_dconc_tgfirst_conc + second_conc). This makes the diagonal contribution -tr/(volrho_l) instead of the correct +tr/(vol*rho_l), i.e. the sign of the entire self-derivative term is inverted for every species (including the target).

Why it's incorrect

Hand differentiation of dcomp_i = -tr/vol*(kron_i - conc_i/rho_l): the chain-rule on the bracket yields -tr/vol * (-delta_ij/rho_l) = +tr/(vol*rho_l)delta_ij. The code produces -tr/(volrho_l)*delta_ij. The off-diagonal/outer-product (first_conc) term is correct and the analogous dfconc_dvol (line 1421) and dfconc_dmu2 (line 1378) are correct, which confirms the leading -1/vol sign convention; only second_conc has the wrong relative sign. This directly corrupts the species-concentration block of the sensitivity RHS (rhs_sensitivity line 697 uses jac_states_vals) and any Newton/Jacobian-based integration that consumes the analytical jac_states.

Suggested regression test

Finite-difference cross-check: set method='moments' on a BatchCryst, call jac_states(time, states, params, return_only=False) and compare against numerical_jac_central(lambda st: unit_model(time, st, params), states). The num_distr:-1 x num_distr:-1 sub-block diagonal entries will have opposite sign to the FD reference; fixing requires second_conc = -tr/rho_l*np.eye(len(w_conc)) (or subtracting it in the bracket).

---

Verification: 2/2 adversarial reviewers confirmed (refuted: 0).

[bernalde]

Independent Codex verification: CONFIRMED

severity_assessment: high
evidence: The RHS at PharmaPy/Crystallizers.py:1512-1513 is -tr/vol * (kron - conc/rho_l). Differentiating the self term gives +tr/(vol*rho_l), but PharmaPy/Crystallizers.py:1404-1408 adds second_conc inside the leading negative, producing the opposite sign; a scalar finite difference confirmed +0.0004 versus code -0.0004.
fix_note: Use second_conc = -tr/rho_l * np.eye(len(w_conc)).

[parkyr]

Reproduction confirmed

Confirmed against current master: the concentration Jacobian self term has the sign described in the existing verification comment.

Environment:

• Base commit: 1a3ab3155216f07e9a18c60687a7b6b1c5959202
• Test command: /home/yirangpark/repos/secquoia/PharmaPy/.venv/bin/python -m pytest tests/test_batch_cryst_concentration_jacobian.py -q

Failing test to add at tests/test_batch_cryst_concentration_jacobian.py:

import sys
from types import ModuleType

import numpy as np
import pytest


pytestmark = pytest.mark.unit


def _stub_assimulo_modules(monkeypatch):
    assimulo = ModuleType("assimulo")

    solvers = ModuleType("assimulo.solvers")
    solvers.CVode = object

    problem = ModuleType("assimulo.problem")
    problem.Explicit_Problem = object

    monkeypatch.setitem(sys.modules, "assimulo", assimulo)
    monkeypatch.setitem(sys.modules, "assimulo.solvers", solvers)
    monkeypatch.setitem(sys.modules, "assimulo.problem", problem)


def _import_batch_cryst(monkeypatch):
    try:
        from PharmaPy.Crystallizers import BatchCryst
    except ModuleNotFoundError as exc:
        if exc.name != "assimulo":
            raise
        _stub_assimulo_modules(monkeypatch)

        from PharmaPy.Crystallizers import BatchCryst

    return BatchCryst


class _Liquid:
    def getDensity(self, temp=None):
        return 1000.0


class _Solid:
    kv = 2.0

    def getDensity(self, temp=None):
        return 1500.0


class _Kinetics:
    prim_nucl = 0.0
    sec_nucl = 0.0
    growth = 2.0e12
    params = {
        "growth": [0.0, 0.0, 1.0],
        "nucl_prim": [0.0, 0.0, 1.0],
        "nucl_sec": [0.0, 0.0, 0.0, 0.0],
    }

    def get_solubility(self, temp, conc):
        return 0.25


def test_batch_cryst_concentration_jacobian_diagonal_matches_material_balance(monkeypatch):
    BatchCryst = _import_batch_cryst(monkeypatch)

    crystallizer = BatchCryst.__new__(BatchCryst)
    crystallizer.num_distr = 4
    crystallizer.num_species = 2
    crystallizer.target_ind = 0
    crystallizer.kron_jtg = np.array([1.0, 0.0])
    crystallizer.controls = {
        "temp": {
            "fun": lambda time: 298.15,
            "args": (),
            "kwargs": {},
        }
    }
    crystallizer.Liquid_1 = _Liquid()
    crystallizer.Solid_1 = _Solid()
    crystallizer._Kinetics = _Kinetics()

    moments = np.array([1.0, 2.0, 4.0, 8.0])
    mass_conc = np.array([0.55, 0.20])
    vol_liq = 2.0
    states = np.concatenate((moments, mass_conc, [vol_liq]))

    jacobian = crystallizer.jac_states(
        time=0.0,
        states=states,
        params=None,
        return_only=False,
    )

    rho_l = crystallizer.Liquid_1.getDensity(temp=298.15)
    rho_c = crystallizer.Solid_1.getDensity(temp=298.15)
    tr = (
        3
        * crystallizer.Solid_1.kv
        * crystallizer.Kinetics.growth
        * moments[2]
        * rho_c
        * (1e-6) ** 3
    )
    dtr_dconc_tg = crystallizer.Kinetics.params["growth"][-1] * tr / (
        mass_conc[crystallizer.target_ind] - crystallizer.Kinetics.get_solubility(298.15, mass_conc)
    )
    first_conc = np.outer(
        crystallizer.kron_jtg - mass_conc / rho_l,
        crystallizer.kron_jtg,
    )
    expected_conc_block = -1 / vol_liq * (
        dtr_dconc_tg * first_conc - tr / rho_l * np.eye(len(mass_conc))
    )

    np.testing.assert_allclose(
        jacobian[crystallizer.num_distr:-1, crystallizer.num_distr:-1],
        expected_conc_block,
        rtol=1e-12,
        atol=1e-12,
    )

Red output:

E   AssertionError:
E   Not equal to tolerance rtol=1e-12, atol=1e-12
E
E   Mismatched elements: 2 / 4 (50%)
E   Max absolute difference: 7.2e-05
E   Max relative difference: 2.
E    x: array([[-1.1997e-01, -0.0000e+00],
E          [ 2.4000e-05, -3.6000e-05]])
E    y: array([[-1.19898e-01, -0.00000e+00],
E          [ 2.40000e-05,  3.60000e-05]])

Observed vs expected:

Quantity	Observed	Expected
Concentration Jacobian block [1, 1]	-3.60000e-05	+3.60000e-05
Concentration Jacobian block [0, 0]	-1.19970e-01	-1.19898e-01

Acceptance criteria:

• BatchCryst.jac_states(..., return_only=False) uses the differentiated material-balance sign for the concentration self-derivative term.
• The regression test above fails on this base commit and passes after the concentration-block sign is corrected.
• Existing non-Assimulo pytest tests continue to pass.