Kinetics API

Chemical kinetics module: rate models, kinetic reactions, ODE problem setup, solvers, and calorimetry.

Rate constants and models

ChemistryLab.KINETICS_RATE_FACTORIES Constant

KINETICS_RATE_FACTORIES

Compiled ThermoFactory objects for each kinetic rate model. Populated by __init__() from KINETICS_RATE_MODELS.

Keys are model name symbols (e.g. :arrhenius). Values are ThermoFactory callables that return SymbolicFunc{1} instances.

Usage

factory = KINETICS_RATE_FACTORIES[:arrhenius]
k = factory(; k₀=1e-5, Ea=50000.0, T_ref=298.15, R_gas=R_GAS)
k(; T = 298.15)   # → 1e-5  (rate constant at reference temperature)

ChemistryLab.KINETICS_RATE_MODELS Constant

KINETICS_RATE_MODELS

Dictionary of raw kinetic rate-constant model expressions, analogous to THERMO_MODELS.

Each entry maps a model name (:arrhenius, …) to a Dict containing:

• :k — symbolic Expr for the rate constant as a function of variables.

• :vars — list of variable symbols (e.g. [:T]).

• :units — list of Symbol => Quantity pairs for parameters and variables.

• :output_unit — Quantity representing the output unit.

At package initialization, every entry is compiled into a ThermoFactory stored in KINETICS_RATE_FACTORIES.

Example

k_acid = KINETICS_RATE_FACTORIES[:arrhenius](;
    k₀    = 5.012e-1,   # mol/(m² s) at T_ref
    Ea    = 14400.0,    # J/mol
    T_ref = 298.15,     # K
)
k_acid(; T = 310.0)   # → Float64 rate constant

ChemistryLab.PK84_PARAMS_C2S Constant

PK84_PARAMS_C2S :: NamedTuple

Canonical Parrot & Killoh (1984) parameters for belite (C₂S = Ca₂SiO₄).

k₁ = 0.5 d⁻¹, n₁ = 1.0, k₂ = 0.006 d⁻¹, k₃ = 0.2 d⁻¹, n₃ = 5.0, Ea = 21 kJ/mol.

With n₁ = 1 the Avrami branch reduces to k₁(1-ξ), which never limits the rate: belite hydration is governed by the power law throughout.

See PK84_PARAMS_C3S.

ChemistryLab.PK84_PARAMS_C3A Constant

PK84_PARAMS_C3A :: NamedTuple

Canonical Parrot & Killoh (1984) parameters for tricalcium aluminate (C₃A = Ca₃Al₂O₆).

k₁ = 1.0 d⁻¹, n₁ = 0.85, k₂ = 0.04 d⁻¹, k₃ = 1.0 d⁻¹, n₃ = 3.2, Ea = 54 kJ/mol.

See PK84_PARAMS_C3S.

ChemistryLab.PK84_PARAMS_C3S Constant

PK84_PARAMS_C3S :: NamedTuple

Canonical Parrot & Killoh (1984) parameters for alite (C₃S = Ca₃SiO₅), valid for a Blaine fineness of 385 m²/kg and a reference temperature of 20 °C.

k₁ = 1.5 d⁻¹, n₁ = 0.7, k₂ = 0.05 d⁻¹, k₃ = 1.1 d⁻¹, n₃ = 3.3, Ea = 42 kJ/mol.

Pass to parrot_killoh_avrami, not to parrot_killoh — the two use different functional forms and their parameters are not transferable.

ChemistryLab.PK84_PARAMS_C4AF Constant

PK84_PARAMS_C4AF :: NamedTuple

Canonical Parrot & Killoh (1984) parameters for tetracalcium aluminoferrite (C₄AF = Ca₄Al₂Fe₂O₁₀).

k₁ = 0.37 d⁻¹, n₁ = 0.7, k₂ = 0.015 d⁻¹, k₃ = 0.4 d⁻¹, n₃ = 3.7, Ea = 32 kJ/mol.

See PK84_PARAMS_C3S.

ChemistryLab.PK_AVRAMI_SEED Constant

PK_AVRAMI_SEED

Lower bound imposed on the normalized degree of hydration inside the Avrami branch of parrot_killoh_avrami, so that the rate is strictly positive at α = 0 and the ODE leaves its degenerate initial point. See the discussion in parrot_killoh_avrami.

One visible consequence: a phase whose rate is governed by the Avrami branch near α = 0 (C₃S, C₃A, C₄AF) starts more slowly than one governed by the power law (C₂S), so alite only overtakes belite after a few minutes. The crossover falls inside the induction period, which this model does not describe anyway.

ChemistryLab.PK_PARAMS_C2S Constant

PK_PARAMS_C2S :: NamedTuple

Parrot & Killoh (1984) parameters for belite (C₂S = Ca₂SiO₄).

Values of the smoothed variant, of unestablished provenance (K₁=0.95, K₂=0.0005, K₃=0.0024 d⁻¹). Activation energy from Schindler & Folliard (2005). Reference temperature: 293.15 K (20 °C).

ChemistryLab.PK_PARAMS_C3A Constant

PK_PARAMS_C3A :: NamedTuple

Parrot & Killoh (1984) parameters for tricalcium aluminate (C₃A = Ca₃Al₂O₆) in the presence of sulfate (gypsum), corresponding to ettringite formation.

Values of the smoothed variant, of unestablished provenance (K₁=0.082, K₂=0.00024, K₃=0.0024 d⁻¹). Activation energy from Schindler & Folliard (2005). Reference temperature: 293.15 K (20 °C).

ChemistryLab.PK_PARAMS_C3S Constant

PK_PARAMS_C3S :: NamedTuple

Parrot & Killoh (1984) parameters for alite (C₃S = Ca₃SiO₅).

Values of the smoothed variant, of unestablished provenance (K₁=1.5, K₂=0.018, K₃=0.0024 d⁻¹). Activation energy from Schindler & Folliard (2005). Reference temperature: 293.15 K (20 °C).

Pass to parrot_killoh to build a KineticFunc:

pk = parrot_killoh(PK_PARAMS_C3S, "C3S")
# or with α_max limit (Powers 1948):
pk = parrot_killoh(PK_PARAMS_C3S, "C3S"; α_max = min(1.0, w_c / 0.42))

ChemistryLab.PK_PARAMS_C4AF Constant

PK_PARAMS_C4AF :: NamedTuple

Parrot & Killoh (1984) parameters for tetracalcium aluminoferrite (C₄AF = Ca₄Al₂Fe₂O₁₀).

Values of the smoothed variant, of unestablished provenance (K₁=0.165, K₂=0.0015, K₃=0.0024 d⁻¹). Activation energy from Schindler & Folliard (2005). Reference temperature: 293.15 K (20 °C).

ChemistryLab.WALLER_PARAMS_FLY_ASH Constant

WALLER_PARAMS_FLY_ASH :: NamedTuple

Waller (1999) parameters for class-F fly ash: τ = 80 d, n = 0.7, blaine_ref = 400 m²/kg, Ea = 83.14 kJ/mol.

Adjusted on SEM image analysis assuming a fly ash of 60 % pozzolanic activity.

Pass to waller.

ChemistryLab.WALLER_PARAMS_SILICA_FUME Constant

WALLER_PARAMS_SILICA_FUME :: NamedTuple

Waller (1999) parameters applied to silica fume — identical kinetics to WALLER_PARAMS_FLY_ASH, the higher reactivity being carried by the fineness. Pass blaine = 2000u"m^2/kg" to waller, the effective value recommended by Lavergne et al. (2018); the BET surface of silica fume (about 20 000 m²/kg) is not a Blaine fineness and must not be used here.

ChemistryLab.WALLER_PARAMS_SLAG Constant

WALLER_PARAMS_SLAG :: NamedTuple

Waller (1999) parameters for ground granulated blast-furnace slag: τ = 100 d, n = 0.7, blaine_ref = 400 m²/kg, Ea = 83.14 kJ/mol.

Slag is latent-hydraulic rather than pozzolanic; the longer characteristic time reflects its slower long-term reaction. Combine with an α_max below 1 (0.9 is customary) to account for the unreactive crystalline fraction.

ChemistryLab.KineticFunc Type

KineticFunc{F, R <: NamedTuple, Q}

Compiled kinetic rate function, analogous to NumericFunc for thermodynamics.

Calling convention (positional, not keyword):

kf(T, P, t, n, lna, n_initial) -> Real   # [mol/s]

where:

• T [K], P [Pa]: temperature and pressure (plain Real or ForwardDiff.Dual).

• t [s]: current time.

• n::StateView: moles of all species (named access: n["C3S"]).

• lna::StateView: log-activities of all species.

• n_initial::StateView: initial moles (always Float64).

• return: net dissolution rate [mol/s], positive = dissolution.

AD-compatible when the compiled closure is AD-compatible.

Examples

julia> idx = Dict("C3S" => 1);

julia> n_sv  = StateView([1.0], idx);

julia> lna_sv = StateView([0.0], idx);

julia> pk = parrot_killoh_avrami(PK84_PARAMS_C3S, "C3S");

julia> pk(293.15, 1e5, 0.0, n_sv, lna_sv, n_sv) > 0
true

ChemistryLab.PoreHumidity Type

PoreHumidity(retention, system; reference, T = temperature(reference))

The internal relative humidity of a sealed paste, computed from the composition it currently has.

Pass it as the humidity keyword of parrot_killoh_avrami or waller and the rate law stops reading a humidity imposed from outside and starts reading the one the material makes for itself. That is what closes the loop: hydration consumes water, the pore space empties, the humidity falls, and humidity_factor throttles the reaction — self-desiccation, which is what Powers' α_max = w/c / 0.42 describes empirically and what powers_alpha_max otherwise supplies as an input.

The humidity is the water activity the retention law returns at the current degree of saturation of the pore space, S = V_liquid / (V_ref − V_solid), both volumes recomputed from the composition the rate law is handed. V_ref is the fresh paste's total volume — the same reference the two-argument porosity uses, and the same sealed-curing convention: the volume the reactions destroy stays inside as empty porosity.

This is where the arrest comes from, and it is kinetic

The Kelvin term is far too small to arrest hydration thermodynamically. Measured on a CEM I paste, imposing a water activity anywhere from 0.95 down to 0.05 leaves the equilibrium assemblage unchanged: the shift is RT ln a_w = −553 J/mol of water at a_w = 0.80, worth about 1.8 kJ per mole of alite against a hydration Gibbs energy of order −100 kJ/mol. Nulling that would need a_w ≈ 5e-6, a Kelvin radius smaller than a water molecule.

A real paste stops at 75–80 % RH because transport and nucleation stop, not because the reaction has become unfavorable. So the humidity belongs in the rate law, through humidity_factor, and CapillaryWater is what makes the water activity of the equilibrium state mean the same thing.

Examples

h = PoreHumidity(retention, cs; reference = fresh)
rxn[:rate] = parrot_killoh_avrami(PK84_PARAMS_C3S, "C3S"; humidity = h)

See also: humidity_factor, WaterRetention, CapillaryWater, powers_alpha_max.

ChemistryLab.RateMechanism Type

struct RateMechanism{F<:AbstractFunc, T<:Real}

A single kinetic mechanism (acid/neutral/base/…) contributing to the overall mineral dissolution or precipitation rate.

The mechanism rate is:

r_mech = k(T) × [Π_catalysts aᵢ^nᵢ] × sign(1 - Ω) × |1 - Ω^p|^q

Fields

• k: rate constant as AbstractFunc (typically SymbolicFunc{1} from arrhenius_rate_constant). Called as k(; T=...).

• p: saturation exponent p in (1 - Ω^p)^q. Default 1.0.

• q: outer exponent q. Default 1.0.

• catalysts: vector of RateModelCatalyst (may be empty).

Examples

k_acid = arrhenius_rate_constant(5.012e-1, 14400.0)
mech   = RateMechanism(k_acid, 1.0, 1.0, [RateModelCatalyst("H+", 1.0)])

ChemistryLab.RateMechanism Method

RateMechanism(k::AbstractFunc, p::Real, q::Real) -> RateMechanism

Construct a RateMechanism with no catalyst contributions.

ChemistryLab.RateModelCatalyst Type

struct RateModelCatalyst{T<:Real}

Describes the contribution of a catalyst species to a reaction mechanism rate.

The catalyst multiplies the base rate by exp(n * ln aᵢ) = aᵢ^n, where aᵢ is the activity of the catalyst species.

Fields

• species: PHREEQC-format formula string of the catalyst species (e.g. "H+", "OH-").

• n: power exponent (dimensionless).

Examples

acid_catalyst   = RateModelCatalyst("H+",  0.5)    # ∝ a(H+)^0.5
base_catalyst   = RateModelCatalyst("OH-", 0.5)    # ∝ a(OH-)^0.5
co2_catalyst    = RateModelCatalyst("CO2", 1.0)    # ∝ a(CO2)

ChemistryLab.StateView Type

StateView{T, I <: AbstractDict}

Thin wrapper giving O(1) named access to a species data vector.

sv["C3S"] === sv.data[sv.index["C3S"]]

The index dict is built once at KineticsProblem construction; data is a plain vector (mutated in-place or re-wrapped each ODE step) — no dict allocation in the hot path.

Examples

julia> idx = Dict("Ca++" => 1, "C3S" => 2);

julia> sv = StateView([0.5, 1.0], idx);

julia> sv["C3S"]
1.0

julia> haskey(sv, "Ca++")
true

ChemistryLab.add_kinetics_rate_model Method

add_kinetics_rate_model(name::Symbol, dict_model::Dict)

Register a new kinetic rate-constant model in KINETICS_RATE_MODELS and compile it into KINETICS_RATE_FACTORIES.

dict_model must contain at minimum :k (expression), :vars (variable list), :units (parameter units), and :output_unit.

Example

add_kinetics_rate_model(:power_law, Dict(
    :k    => :(k₀ * (T / T_ref)^n),
    :vars => [:T],
    :units => [:T => u"K", :T_ref => u"K", :k₀ => u"mol/(m^2*s)", :n => u"1"],
    :output_unit => u"mol/(m^2*s)",
))

ChemistryLab.arrhenius_rate_constant Method

arrhenius_rate_constant(k₀, Ea; T_ref=298.15, R_gas=R_GAS) -> NumericFunc

Build a temperature-dependent Arrhenius rate constant as a NumericFunc:

k(T) = k₀ × exp(-Eₐ / R × (1/T - 1/T_ref))

The returned object is callable as k(; T=...) and fully AD-compatible (ForwardDiff-safe: the closure captures k₀, Ea, T_ref, R_gas directly, so dual numbers propagate correctly through all parameters).

Arithmetic between SymbolicFunc/NumericFunc objects is supported, so rate constants can be composed with activity or surface-area functions.

Arguments

• k₀: pre-exponential factor at T_ref. Plain Real → SI [mol/(m² s)]; Quantity → automatically converted (e.g. 5e-4u"mol/(m^2*s)").

• Ea: activation energy. Plain Real → SI [J/mol]; Quantity → converted (e.g. 62.0u"kJ/mol").

• T_ref: reference temperature. Plain Real → SI [K]; Quantity → converted (e.g. 298.15u"K"). Default 298.15.

• R_gas: gas constant [J/(mol K)] (plain Real only; default R_GAS, the CODATA value taken from DynamicQuantities.Constants).

Returns

A NumericFunc with variable T (in K) and refs = (T = T_ref * u"K",).

Examples

julia> k = arrhenius_rate_constant(5.0e-4, 62000.0);

julia> isapprox(k(; T = 298.15), 5.0e-4; rtol = 1e-10)
true

julia> k(; T = 350.0) > k(; T = 298.15)   # higher T → higher k
true

Unit-aware: k₀ in mmol/(m²·s), Ea in kJ/mol, T_ref in K — all converted to SI:

k = arrhenius_rate_constant(0.5u"mmol/(m^2*s)", 62.0u"kJ/mol"; T_ref = 298.15u"K")

AD-compatible through all parameters:

ForwardDiff.derivative(T  -> arrhenius_rate_constant(5e-4, 62000.0)(; T = T),  298.15)
ForwardDiff.derivative(Ea -> arrhenius_rate_constant(5e-4, Ea)(; T = 350.0),   62000.0)
ForwardDiff.derivative(k₀ -> arrhenius_rate_constant(k₀,   62000.0)(; T = 298.15), 5e-4)

ChemistryLab.blaine_factor Method

blaine_factor(blaine; blaine_ref = 385u"m^2/kg") -> Real

Fineness correction of the hydration rate: the Parrot & Killoh parameters were adjusted for a cement of Blaine fineness blaine_ref, and the rate scales as blaine / blaine_ref.

Both arguments accept a DynamicQuantities.Quantity or a plain Real in m²/kg. The default reference is 385 m²/kg for clinker phases; pass blaine_ref = 400u"m^2/kg" for the Waller kinetics of additions.

Examples

julia> blaine_factor(385u"m^2/kg")
1.0

julia> round(blaine_factor(462u"m^2/kg"); digits = 4)
1.2

ChemistryLab.humidity_factor Method

humidity_factor(h) -> Real

Reduction coefficient β_h applied to the hydration rate at internal relative humidity h ∈ [0, 1] (Parrot et al., as used by van Breugel):

Hydration is taken to stop below 80 % relative humidity, on thermodynamic grounds. The cut is a genuine discontinuity: the one-sided limit from above is β_h(0.80⁺) ≈ 0.0953 while the value at and below 0.80 is exactly 0. The jump is mild in practice, the rate having already fallen by an order of magnitude from β_h(0.99) ≈ 0.914.

Examples

julia> humidity_factor(0.75)
0.0

julia> humidity_factor(0.80)
0.0

julia> round(humidity_factor(0.801); digits = 4)
0.0968

ChemistryLab.parrot_killoh Method

parrot_killoh(params::NamedTuple, mineral_name::AbstractString; α_max=1.0) -> KineticFunc

Build a smoothed three-mechanism clinker hydration rate as a KineticFunc.

Deprecated, and no longer attributed to Parrott & Killoh

Use parrot_killoh_avrami with PK84_PARAMS_C3S and siblings instead. This function is kept so that existing scripts keep running, and it warns once per session.

Why the attribution is withdrawn. The formulas below are not those of Parrott & Killoh: the nucleation–growth term carries no Avrami logarithm, and K₃ — a shell-formation coefficient — sits in the diffusion expression where the canonical formulation uses K₂. Nor do the shipped parameters match any published set: N₁ = 3.3 is the canonical n₃, and the canonical k₃ = 1.1 has no counterpart at all. The primary source is a conference proceedings without a DOI (British Ceramic Proceedings 35, 41–53, 1984) that could not be consulted, so the attribution is retracted rather than repaired by an invented calibration.

Why it matters. With PK_PARAMS_* the diffusion branch takes over very early — measured at α/α_max = 0.003 for C₂S, 0.013 for C₃S, 0.057 for C₃A, and over the whole range for C₄AF — and governs throughout the interval a seven-day run traverses, because its prefactor 3K₃/N₃ = 0.0018 d⁻¹ is 28 times smaller than the canonical k₂ = 0.05 d⁻¹. The rate then integrates in closed form,

α(t) = α_max · [1 − (1 − √(2·K₃·t / N₃))³]

which is independent of the phase, because K₃ = 0.0024 d⁻¹ and N₃ = 4 are identical in all four parameter sets. Measured on a CEM I at w/c = 0.40 over seven days, C₃S, C₂S and C₃A all land on α ≈ 0.239, and C₄AF lower still at 0.193 — there its own nucleation-growth branch is slower than diffusion and limits instead. The weighted mean comes to 0.234 against the 0.61 the cement literature reports. The signature is unmistakable: K₁ spans a factor of 18 across the four phases and changes almost nothing.

params must be a NamedTuple with keys K₁, N₁, K₂, N₂, K₃, N₃, B, Ea, T_ref. All dimensional values accept plain Real (SI) or DynamicQuantities.Quantity.

mineral_name is the PHREEQC formula string (e.g. "C3S") used to look up the mineral moles in the n and n_initial StateViews.

Three competing mechanisms determine the rate (Parrot & Killoh 1984):

Mechanism	Formula
Nucleation–growth	r_NG = (K₁/N₁)(1-ξ)^N₁ / (1 + B·ξ^N₃)
Interaction	r_I = K₂(1-ξ)^N₂
Diffusion	r_D = 3K₃(1-ξ)^(2/3) / (N₃·(1-(1-ξ)^(1/3)))

The rate [mol/s] is n_initial × Aₜ × min(max(r_NG, r_I), r_D) where ξ = α / α_max is the normalized degree of hydration and Aₜ = exp(-Ea/R × (1/T - 1/T_ref)) is the Arrhenius factor.

α_max can be set to apply the Powers (1948) water/cement ratio limit: α_max = min(1.0, w_c / 0.42).

Returns

A KineticFunc — callable as pk(T, P, t, n::StateView, lna::StateView, n_initial::StateView) -> Real [mol/s]. AD-compatible (ForwardDiff-safe): no Float64 casts in the evaluation path.

Examples

Shown rather than run: the constructor warns by design, and a deprecation warning on stderr is not something a doctest should have to reproduce verbatim. The numerical behavior of this variant is pinned by the test suite instead — see the branch oracle in test/kinetics/test_rate_models.jl, which asserts that PK_PARAMS_* puts C3S, C2S and C3A on the diffusion branch within a few percent of hydration and that C4AF is limited by its own nucleation branch.

pk = parrot_killoh(PK_PARAMS_C3S, "C3S")
idx = Dict("C3S" => 1)
n0 = StateView([1.0], idx)
lna = StateView([0.0], idx)
pk(293.15, 1e5, 0.0, n0, lna, n0) > 0      # true

See also: PK_PARAMS_C3S, PK_PARAMS_C2S, PK_PARAMS_C3A, PK_PARAMS_C4AF.

ChemistryLab.parrot_killoh_avrami Method

parrot_killoh_avrami(params::NamedTuple, mineral_name::AbstractString;
                     α_max = 1.0, blaine = nothing, humidity = nothing) -> KineticFunc

Build the Parrot & Killoh (1984) clinker hydration rate in its canonical formulation, as reported by Lothenbach et al. (2008) and used by Lavergne et al. (2018).

params must be a NamedTuple with keys k₁, n₁, k₂, k₃, n₃, Ea, T_ref — see PK84_PARAMS_C3S and siblings. Dimensional values accept plain Real (SI) or DynamicQuantities.Quantity.

Three competing mechanisms limit the rate, and the slowest one wins:

Mechanism	Formula
Nucleation–growth (Avrami)	α̇₁ = (k₁/n₁)(1-ξ)(-ln(1-ξ))^(1-n₁)
Diffusion (Jander)	α̇₂ = k₂(1-ξ)^(2/3) / (1-(1-ξ)^(1/3))
Shell formation (power law)	α̇₃ = k₃(1-ξ)^n₃

so that α̇ = min(α̇₁, α̇₂, α̇₃), with ξ = α/α_max the normalized degree of hydration. The returned rate [mol/s] is n_initial × Aₜ × β_B × β_h × α̇, where Aₜ = exp(-Ea/R × (1/T - 1/T_ref)) is the Arrhenius factor, β_B the Blaine fineness factor (blaine_factor) and β_h the relative-humidity reduction (humidity_factor). Both default to 1 when their keyword is nothing.

Two Parrot–Killoh variants ship with ChemistryLab

parrot_killoh implements a different, smoothed variant (min(max(r_NG, r_I), r_D) with a B-damped nucleation term) together with the parameter set of PK_PARAMS_C3S and siblings. The two are not interchangeable: their parameters are not transferable, and only parrot_killoh_avrami with PK84_PARAMS_C3S reproduces the α(t) curves published in the cement literature cited above.

With the canonical parameters, C₂S has no nucleation–growth stage and C₃S has no diffusion-controlled stage — an artifact of the 1984 fit that the original authors acknowledged, and a convenient signature to check an implementation against.

Keyword arguments

• α_max: Powers (1948) water availability cap — see powers_alpha_max.

• blaine: Blaine fineness of the binder, as a Quantity or a plain Real in m²/kg. nothing (default) means no correction.

• humidity: internal relative humidity, either a constant in [0, 1] or a callable t -> h(t). nothing (default) means no correction.

Returns

A KineticFunc — callable as pk(T, P, t, n::StateView, lna::StateView, n_initial::StateView) -> Real [mol/s]. AD-compatible (ForwardDiff-safe): no Float64 casts in the evaluation path.

Examples

julia> pk = parrot_killoh_avrami(PK84_PARAMS_C3S, "C3S"; blaine = 380u"m^2/kg");

julia> idx = Dict("C3S" => 1);

julia> n0 = StateView([1.0], idx);

julia> lna = StateView([0.0], idx);

julia> pk(293.15, 1e5, 3600.0, StateView([0.9], idx), lna, n0) > 0
true

See also: PK84_PARAMS_C3S, waller, blaine_factor, humidity_factor, powers_alpha_max.

ChemistryLab.pore_saturation Method

pore_saturation(h::PoreHumidity, n) -> Real

Degree of saturation of the pore space at composition n, V_liquid / V_pore with V_pore = V_ref − V_solid. Clamped to [0, 1].

ChemistryLab.powers_alpha_max Method

powers_alpha_max(w_c; curing = :sealed) -> Real

Powers (1948) upper bound on the degree of hydration set by the availability of water, α_max = min(1, w/c / k): a paste below w/c = k cannot hydrate completely, k being 0.42 sealed or 0.36 water-cured, according to curing.

The 0.42 is not a stoichiometric demand, and reading it as one leads to the wrong conclusion about what a Gibbs minimization should return. It is about 0.23 g of non-evaporable water per gram of cement — the water written into the hydrate formulae, which is a mass balance — plus about 0.19 g of gel water held in the C-S-H gel pores, which is physically present and chemically unavailable. In a sealed paste hydration stops by self-desiccation with water still in the specimen, so this bound is a statement about transport and access, not about thermodynamics: an equilibrium calculation on the same mix consumes all the clinker well below 0.42, and only runs out of water near the stoichiometric demand.

The two curing conventions

curing = :sealed is a specimen that exchanges nothing with its surroundings — the convention of porosity and of the w/c example. curing = :saturated is a specimen kept under water after setting, free to draw in what the chemical shrinkage empties; the bound is then 0.36 and a mix that would arrest sealed can go on reacting. Neither is a property of the cement: they are two boundary conditions on the same paste, and which one applies is the caller's to state.

Where water is abundant — w/c above the coefficient, or a cure that keeps supplying it — the bound is 1 and this function stops doing anything, which is the correct answer rather than a degenerate case: nothing about water is then limiting the reaction.

Pass the result as the α_max keyword of parrot_killoh, parrot_killoh_avrami or waller.

It is a ceiling, not a schedule

α_max says how far the reaction can go, never how far it has got. At an early age the degree of reaction is set by the kinetics and is far below this bound; the bound binds only at long times, and only for the constituents whose kinetics would otherwise have taken them past it. For a constituent that reacts slowly — a slag, and a fly ash still more — the binding limit at 28 days is its own dissolution rate, not the water. Take the smaller of the two.

Examples

julia> powers_alpha_max(0.5)
1.0

julia> round(powers_alpha_max(0.32); digits = 4)
0.7619

julia> round(powers_alpha_max(0.32; curing = :saturated); digits = 4)
0.8889

ChemistryLab.saturation_ratio Method

saturation_ratio(stoich::AbstractVector, lna::AbstractVector,
                 ΔₐG⁰overRT::AbstractVector; ϵ=1e-16) -> Real

Compute the saturation ratio Ω = IAP / K for a kinetic reaction.

ln Ω = Σᵢ νᵢ ln aᵢ − ln K
     = Σᵢ νᵢ ln aᵢ + ΔᵣG⁰/(RT)   (note: ln K = −ΔᵣG⁰/RT = −Σᵢ νᵢ ΔₐG⁰ᵢ/RT)

where stoich[i] is the stoichiometric coefficient (positive for products, negative for reactants), lna[i] is the log-activity of species i, and ΔₐG⁰overRT[i] is the dimensionless standard Gibbs energy of formation ΔₐG⁰ᵢ / RT for species i.

Arguments

• stoich: stoichiometric coefficient vector for this reaction (length = number of species).

• lna: log-activity vector (same indexing as species in system).

• ΔₐG⁰overRT: dimensionless standard Gibbs energies ΔₐG⁰ᵢ/RT.

• ϵ: floor to avoid exp overflow when Ω → ∞.

Returns

Ω = exp(ln_IAP - ln_K) where ln_K = -ΔᵣG⁰/RT.

AD-compatible (ForwardDiff-safe).

ChemistryLab.waller Method

waller(params::NamedTuple, mineral_name::AbstractString;
       α_max = 1.0, blaine = nothing, humidity = nothing) -> KineticFunc

Build the Waller (1999) reaction rate of a pozzolanic or latent-hydraulic addition (fly ash, silica fume, ground granulated slag) as a KineticFunc.

The degree of reaction follows a sigmoid in log-time,

whose rate, written as a function of the current degree so that it composes with temperature, fineness and humidity corrections, is

params must be a NamedTuple with keys τ, n, Ea, T_ref and, optionally, blaine_ref — see WALLER_PARAMS_FLY_ASH.

Pozzolanic reactions are markedly more temperature-sensitive than the hydraulic reactions of clinker: the shipped activation energy is 83.14 kJ/mol against 21–54 kJ/mol for the clinker phases.

Silica fume is far finer than the cement (about 20 000 m²/kg by BET). Its reactivity is nonetheless represented here through the Blaine fineness, for which an effective 2000 m²/kg is the recommended default (WALLER_PARAMS_SILICA_FUME) — the two measurements probe different physical phenomena and are not interchangeable.

Keyword arguments

Identical to parrot_killoh_avrami. The Blaine correction is taken relative to params.blaine_ref (400 m²/kg for fly ash), not to the clinker reference of 385 m²/kg.

Examples

julia> fa = waller(WALLER_PARAMS_FLY_ASH, "FlyAsh");

julia> idx = Dict("FlyAsh" => 1);

julia> n0 = StateView([1.0], idx);

julia> lna = StateView([0.0], idx);

julia> fa(293.15, 1e5, 86400.0, StateView([0.9], idx), lna, n0) > 0
true

See also: WALLER_PARAMS_FLY_ASH, WALLER_PARAMS_SILICA_FUME, WALLER_PARAMS_SLAG, parrot_killoh_avrami.

Kinetic reactions and surface area

ChemistryLab.AbstractSurfaceModel Type

abstract type AbstractSurfaceModel end

Base type for models that compute the reactive surface area of a mineral phase.

Concrete subtypes must implement:

surface_area(model, n::Real, molar_mass::Real) -> Real

returning the reactive surface area in m².

All methods must be AD-compatible (no Float64 casts).

ChemistryLab.BETSurfaceArea Type

struct BETSurfaceArea{T<:Real} <: AbstractSurfaceModel

Reactive surface area that scales with the mineral mass, following a BET (Brunauer-Emmett-Teller) specific-surface-area measurement.

A = A_spec × n × M_mineral       [m²]

where A_spec [m²/kg] is the specific BET surface area, n is the current molar amount [mol], and M_mineral is the molar mass [kg/mol].

This is the standard approach in reactive-transport models (Palandri & Kharaka 2004).

Fields

• A_specific: specific BET surface area [m²/kg].

Examples

BETSurfaceArea(90.0)              # 90 m²/kg  (plain Real → SI)
BETSurfaceArea(0.09u"m^2/g")     # 0.09 m²/g → 90 m²/kg

ChemistryLab.BETSurfaceArea Method

BETSurfaceArea(A_specific) -> BETSurfaceArea

Construct a BETSurfaceArea. A_specific can be a plain Real (SI [m²/kg]) or a Quantity (automatically converted to m²/kg), e.g. 0.09u"m^2/g" → 90 m²/kg.

ChemistryLab.FixedSurfaceArea Type

struct FixedSurfaceArea{T<:Real} <: AbstractSurfaceModel

Constant reactive surface area, independent of mineral abundance.

Suitable for short simulations or when the surface area is externally controlled (e.g. from BET measurements on a fixed mass of powder).

Fields

• A: total reactive surface area [m²].

Examples

FixedSurfaceArea(0.5)            # 0.5 m²  (plain Real → SI)
FixedSurfaceArea(500.0u"cm^2")  # 500 cm² → 0.05 m²

ChemistryLab.FixedSurfaceArea Method

FixedSurfaceArea(A) -> FixedSurfaceArea

Construct a FixedSurfaceArea. A can be a plain Real (SI [m²]) or a Quantity (automatically converted to m²).

ChemistryLab.KineticReaction Type

struct KineticReaction{R<:AbstractReaction, F, H}

Associates a chemical Reaction with a compiled KineticFunc.

Following Leal et al. (2017), reactions — not individual species — carry kinetics. A single mineral can therefore appear as a reactant in multiple KineticReaction objects (e.g. C₃A → ettringite and C₃A → monosulphate for multi-pathway cement hydration). The ODE state is indexed by unique mineral species, and contributions from all reactions that consume the same mineral are accumulated.

Fields

• reaction: the underlying Reaction / CemReaction.

• rate_fn: a KineticFunc (or any callable matching the six-argument signature (T, P, t, n, lna, n_initial) -> Real) computing r [mol/s].

• idx_mineral: index of the primary (controlling) mineral species in the parent ChemicalSystem. Determined automatically as the first solid (AS_CRYSTAL) reactant.

• stoich: stoichiometric coefficient vector for all species in the system. Sign convention: positive for products, negative for reactants.

• heat_per_mol: enthalpy of reaction [J/mol], positive = exothermic (heat released). When nothing (default), the enthalpy is derived from the stoichiometric sum of species :ΔₐH⁰ values.

Constructors

From a species name (convenience, builds a minimal dissolution Reaction):

pk = parrot_killoh(PK_PARAMS_C3S, "C3S")
kr = KineticReaction(cs, "C3S", pk)
kr = KineticReaction(cs, "C3S", pk; heat_per_mol = 114_634.0)

From an explicit Reaction (multi-pathway):

pk_c3a = parrot_killoh(PK_PARAMS_C3A, "C3A")
kr_ett  = KineticReaction(cs, rxn_C3A_ettringite,   pk_c3a)
kr_mono = KineticReaction(cs, rxn_C3A_monosulphate, pk_c3a)

Reaction-centric (rate stored in rxn.properties[:rate]):

rxn[:rate] = parrot_killoh(PK_PARAMS_C3S, "C3S")
kr = KineticReaction(cs, rxn)

ChemistryLab.KineticReaction Method

KineticReaction(cs::ChemicalSystem, rxn::AbstractReaction) -> KineticReaction

Reaction-centric constructor: build a KineticReaction from a Reaction that carries its kinetics in reaction.properties.

Required property:

• rxn[:rate] — a KineticFunc or any callable matching (T, P, t, n, lna, n_initial) -> Real. Non-KineticFunc callables are wrapped automatically in a KineticFunc with empty refs.

Optional property:

• rxn[:heat_per_mol] — a Number giving the molar enthalpy [J/mol] for calorimetry.

Examples

pk = parrot_killoh(PK_PARAMS_C3S, "C3S")
rxn[:rate]         = pk
rxn[:heat_per_mol] = 114_634.0
kr = KineticReaction(cs, rxn)

# Build problem directly from a list of annotated Reaction objects:
kp = KineticsProblem(cs, [rxn_C3S, rxn_C3A, rxn_C2S], state0, tspan)

ChemistryLab.KineticReaction Method

KineticReaction(cs::ChemicalSystem, species_name::AbstractString, rate_fn;
                heat_per_mol=nothing) -> KineticReaction

Convenience constructor: look up species_name in cs and build a minimal dissolution Reaction (species as sole reactant, no products) automatically.

The default stoichiometry places -1.0 at the mineral index and 0.0 everywhere else.

ChemistryLab.KineticReaction Method

KineticReaction(rxn, rate_fn, idx_mineral, stoich; heat_per_mol=nothing)

Low-level constructor: explicit Reaction, rate callable, index, and stoichiometry.

ChemistryLab.KineticReaction Method

KineticReaction(cs::ChemicalSystem, rxn::AbstractReaction, rate_fn;
                heat_per_mol=nothing) -> KineticReaction

Construct a KineticReaction from an explicit Reaction object.

The controlling mineral index (idx_mineral) is determined automatically as the index of the first solid (AS_CRYSTAL) reactant found in rxn.reactants that is present in cs. The stoichiometric vector is derived from the reaction stoichiometry.

This constructor is the recommended entry point for multi-pathway kinetics:

pk_c3a = parrot_killoh(PK_PARAMS_C3A, "C3A")
kr_ett  = KineticReaction(cs, cs.dict_reactions["C3A_ettringite"],   pk_c3a)
kr_mono = KineticReaction(cs, cs.dict_reactions["C3A_monosulphate"], pk_c3a)
kp = KineticsProblem(cs, [kr_C3S, kr_ett, kr_mono], state0, tspan)

ChemistryLab.first_order_rate Method

first_order_rate(k, cs, rxn, surface_model; p=1.0, q=1.0, ϵ=1e-16) -> KineticFunc

Build a single-mechanism first-order TST rate as a KineticFunc.

r = A(n) × k(T) × sign(1 - Ω) × |1 - Ω^p|^q

This is a convenience wrapper around transition_state with one no-catalyst mechanism. Useful as a minimal test case or for empirical fits.

Arguments

• k: rate constant as an AbstractFunc (e.g. from arrhenius_rate_constant).

• cs, rxn, surface_model: same as transition_state.

• p, q: saturation exponents (defaults 1.0).

• ϵ: regularization floor (default 1e-16).

Examples

k = arrhenius_rate_constant(1e-7, 40000.0)
rf = first_order_rate(k, cs, rxn, BETSurfaceArea(90.0))
kr = KineticReaction(cs, rxn, rf)

ChemistryLab.molar_mass Method

molar_mass(kr::KineticReaction) -> Float64

Return the molar mass of the mineral species [kg/mol], used internally for BETSurfaceArea calculations.

Searches kr.reaction.reactants for a species with an :M property. Falls back to 0.1 kg/mol when :M is unavailable.

ChemistryLab.surface_area Method

surface_area(model::BETSurfaceArea, n::Real, molar_mass::Real) -> Real

Return A_specific × n × molar_mass [m²]. Clamps at zero to avoid negative surface areas when n → 0. AD-compatible.

ChemistryLab.surface_area Method

surface_area(model::FixedSurfaceArea, n::Real, molar_mass::Real) -> Real

Return the fixed surface area model.A [m²], independent of moles n. AD-compatible.

ChemistryLab.transition_state Method

transition_state(mechanisms, cs, rxn, surface_model; ϵ=1e-16) -> KineticFunc

Build a Transition-State Theory (TST) dissolution/precipitation rate function from a list of RateMechanism objects, returning a KineticFunc.

The compiled closure captures:

• the mineral name and molar mass (from rxn + cs)

• the surface model (FixedSurfaceArea or BETSurfaceArea)

• stoichiometry and ΔₐG⁰ callables for all aqueous species (T-dependent Ω)

The net rate [mol/s] is:

r = A(n) × Σ_m [ k_m(T) × Π_cat(aᵢ^nᵢ) × (1 - Ω^p) × |1 - Ω^p|^(q-1) ]

where Ω(T) = exp(Σ νᵢ ln aᵢ + Σ νᵢ ΔₐG°ᵢ(T)/(RT)) is re-evaluated at every ODE step — correct for variable-temperature semi-adiabatic calorimetry.

Arguments

• mechanisms: vector of RateMechanism (acid, neutral, base, …).

• cs: ChemicalSystem supplying ΔₐG⁰ callables for aqueous species.

• rxn: AbstractReaction defining stoichiometry and the mineral species.

• surface_model: AbstractSurfaceModel — captures area as a function of n.

• ϵ: regularization floor near Ω = 1 (default 1e-16).

Returns

A KineticFunc callable as f(T, P, t, n::StateView, lna::StateView, n_initial::StateView) -> Real [mol/s].

AD-compatible: all operations use generic Julia arithmetic; no Float64 casts.

References

• Palandri, J.L. & Kharaka, Y.K. (2004). USGS Open-File Report 2004-1068.

• Leal, A.M.M. et al. (2017). Pure Appl. Chem. 89, 597–643.

Kinetics problem

ChemistryLab.EQ_RESIDUAL_TOL Constant

EQ_RESIDUAL_TOL

Largest per-element balance violation, relative to that element's own total, that still counts as a converged speciation. Above it the point is not handed on as a warm start and the solve is retried from a budget-clipped guess.

ChemistryLab.RESTORE_MAXIT Constant

RESTORE_MAXIT

Alternating-projection sweeps allowed when restoring the feasibility of an in-run guess. Exposed because the right value is a trade: the projection converges linearly and a cement can need tens of thousands of sweeps, while this runs at every right-hand-side evaluation.

ChemistryLab._CONTINUATION_STEPS Constant

_CONTINUATION_STEPS

How many bisection steps a replay may take between the last certified instant and one it cannot certify directly. Each step halves the jump in the component totals, so eight of them reduce it by a factor 256.

ChemistryLab._EQ_GUESS_FLOOR Constant

_EQ_GUESS_FLOOR

Lower floor applied to the starting guess of the equilibrium sub-solve, chosen strictly above the 1e-16 lower bound that EquilibriumProblem imposes.

An interior-point method started on its own bound stalls short of its tolerance and returns MaxIters. Loosening that tolerance is not the fix: measured against Reaktoro on the calcite reference case, the worst species error grows from 4.3 % at tol = 1e-10 to 38 % at 1e-8 and 252 % at 1e-7. Moving the guess inside keeps the tight tolerance and removes the stalls.

ChemistryLab._RETRY_ABS_TOL Constant

_RETRY_ABS_TOL

Element-balance violation in moles above which a replayed speciation is solved a second time from a guess carrying no active set. Chosen well above machine precision and well below anything chemically meaningful.

ChemistryLab.KineticsProblem Type

struct KineticsProblem{CS, CAL, ES, AM}

Encapsulates a kinetics simulation following Leal et al. (2017).

The ODE state vector u is structured as:

• Without re-speciation: u = [nₖ₁, …, nₖ_K, ξ₁, …, ξ_M, [T | Q]]

• With re-speciation: u = [bₑ₁, …, bₑ_C, nₖ₁, …, nₖ_K, ξ₁, …, ξ_M, [T | Q]]

where bₑ are the element amounts in the equilibrium partition, nₖ the moles of kinetic species and ξ the extents of the kinetic reactions. The trailing slot is present only with a calorimeter: the temperature for SemiAdiabaticCalorimeter, the accumulated heat for IsothermalCalorimeter.

Fields

• system: ChemicalSystem.

• kinetic_reactions: vector of KineticReaction objects.

• initial_state: ChemicalState providing initial moles, T, P.

• tspan: (t_start, t_end) time interval [s].

• calorimeter: nothing, IsothermalCalorimeter, or SemiAdiabaticCalorimeter.

• activity_model: AbstractActivityModel for log-activities.

• equilibrium_solver: solver for re-speciation, or nothing.

• idx_kinetic: indices of kinetic species in system.species.

• idx_equilibrium: indices of equilibrium species.

• ν: stoichiometric matrix (M × N) = SM.N' restricted to kinetic reactions.

• νe, νk: partitions of ν for equilibrium / kinetic species.

• Ae: formula matrix restricted to equilibrium species (C × Nₑ).

See also: integrate, KineticsSolver.

ChemistryLab._abs_residual Method

_abs_residual(Ae, n_e, be) -> Float64

Largest element-balance violation in moles. Reported alongside the relative measure because it is the one a chemist can judge: 1e-10 mol is machine precision whatever the system, and 7e-2 mol is not.

ChemistryLab._budget_clip! Method

_budget_clip!(n_eq, Ae, be)

Clip a starting guess to the element budget: no species may exceed what the totals bₑ can supply, nⱼ ≤ minᵢ bᵢ/Aᵢⱼ over the rows it consumes.

This does not change the feasible set — it only moves the guess into it. It is what unblocks the OPC case: once the sulfate is spent the warm start still carried ettringite at the aluminum budget, three times the sulfur available, and the interior-point solve could not walk back from there.

ChemistryLab._build_kinetics_problem Method

KineticsProblem(cs, kinetic_reactions, initial_state, tspan; ...) -> KineticsProblem

Construct a KineticsProblem from an explicit list of reactions.

Each element of kinetic_reactions must be either a KineticReaction or a Reaction with a :rate entry in its properties. Reaction objects are automatically wrapped via KineticReaction(cs, rxn).

KineticsProblem(cs, initial_state, tspan; ...) -> KineticsProblem

Construct from a ChemicalSystem that has kinetic_species declared (reactions and rates auto-generated via the kinetic_species keyword).

Arguments

• cs: ChemicalSystem.

• kinetic_reactions: AbstractVector of KineticReaction or Reaction objects carrying a :rate property.

• initial_state: ChemicalState providing initial moles, T, P.

• tspan: (t0, tf) time interval. Plain Real → [s]; Quantity → converted.

• calorimeter: nothing (no thermal coupling), IsothermalCalorimeter, or SemiAdiabaticCalorimeter.

• activity_model: activity model for log-activity computation (default: dilute).

• equilibrium_solver: nothing (no re-speciation) or an EquilibriumSolver.

Examples

# From explicit reactions
rxn = Reaction(OrderedDict(sp("C3S") => 1.0, sp("H2O@") => 3.33),
               OrderedDict(sp("Jennite") => 0.167, sp("Portlandite") => 1.5))
rxn[:rate] = parrot_killoh(PK_PARAMS_C3S, "C3S"; α_max)
kp = KineticsProblem(cs, [rxn], state0, (0.0, 7 * 86400.0))

# From kinetic_species in ChemicalSystem
cs = ChemicalSystem(species, primaries;
    kinetic_species = Dict("C3S" => pk_C3S, "C2S" => pk_C2S))
kp = KineticsProblem(cs, state0, (0.0, 7 * 86400.0))

ChemistryLab._equilibrium_subsystem Method

_equilibrium_subsystem(system, idx_equilibrium) -> ChemicalSystem

The chemical system restricted to the equilibrium partition, as the partitioned formulation of (Leal et al., 2017) requires.

Its formula matrix is exactly system.CSM.A[:, idx_equilibrium], in the same species order, so the element amounts bₑ carried by the ODE state are handed to it unchanged. The primaries are those of the parent system that survive the restriction — the kinetic minerals never do, not being in the partition.

ChemistryLab._one_speciation Method

_one_speciation(p, guess, be) -> (ok, n_e, abs_res)

One equilibrium solve of φ(bₑ) from guess. Returns ok = false only if the solve threw; a solve that returns a poor composition still returns true, with its element-balance violation in moles, so the caller can compare attempts.

ChemistryLab._reconstruction_guess! Method

_reconstruction_guess!(buf, p, be) -> buf

The independent guess: the composition the specimen was cast with, carried onto the current element budget. It holds no active set at all, which is exactly what recommends it where the assemblage is switching.

ChemistryLab._respeciate_solve! Method

_respeciate_solve!(p, n_eq, be) -> Bool

Solve φ(bₑ) from the guess n_eq, write the result into p.n_full, and record the element-balance residual. Returns false if the solve threw.

Two guesses are available and the order between them is chosen, not fixed. The warm one is the previous speciation; the reconstruction is the cast composition carried onto the current budget. Whichever runs first, the other is tried when the first leaves too much matter unaccounted for, and the better of the two is kept — so the answer does not depend on the order, only the cost does.

WHY THE ORDER IS WORTH CHOOSING. The warm start carries the previous ACTIVE SET, and where the assemblage switches it is the wrong one: an interior-point method started inside a set of phases that no longer exists does not cross over, it exhausts its iterations. Measured on a six-hour hydration, 143 respeciations:

	first solve	then
the 98 calls the warm guess handled	1247 ms	—
the 45 calls where the assemblage switched	8087 ms, rejected	1467 ms from the reconstruction

Eight seconds spent learning that a guess is wrong, then one and a half to get the answer without it. That was 66 % of the whole integration.

An assemblage switch is not a single-step event — a phase takes several steps to exhaust — so the previous call's outcome says which guess to try first, and eq_switching carries it. This is information the problem already has, not a tuning parameter: no threshold is introduced, and the tolerance that decides "too much matter unaccounted for" is the same _RETRY_ABS_TOL as before.

ChemistryLab._restore_feasibility! Method

_restore_feasibility!(n_eq, Ae, be; maxit, tol) -> n_eq

Move a starting guess into {Aₑn = bₑ, n ≥ 0} by alternating projection: clamp to the box, then project onto the affine set through the 7×7 system AₑAₑᵀ.

The Gibbs minimization is posed with A n = b as a hard equality, so a guess that violates it starts the interior-point method outside its own feasible set. On a full OPC this was not a detail: the solve stopped on a point demanding 0.732 mol of sulfate against the 0.267 mol available — 174 % over — and no tolerance, iteration budget, barrier setting or bound changed it, because none of them addresses an infeasible start. A projected-gradient phase-1 proved the set is non-empty (residual 3·10⁻¹²); this is the cheap way to land in it.

Ending on the box clamp rather than the affine step matters: the affine projection alone leaves small negative amounts, which the barrier cannot accept.

maxit deserves attention. Alternating projection converges linearly, at a rate set by the angle between the box and the affine set, and on a cement that angle can be small: at the six-hour instant of an ordinary Portland cement — where the iron row carries 0.013 mol across thirteen species — 200 sweeps left a residual of 8.4e-1, 2000 left 6.7e-2, and 20 000 were needed to reach 6.7e-9.

The default is small on purpose, and measured: inside the ODE right-hand side this runs at every evaluation, and on a full OPC the worst in-run balance is 1.1 mol at 200 sweeps against 8.5 at 2000 and 41 at 100 000. A better guess producing a worse answer is the back-end's own unpredictability; until that is understood the in-run budget stays where it measures best. Note that the ranking depends on the back-end and should be re-measured if it changes.

A replay (speciated_states) runs a handful of times and buys accuracy instead, asking for a much larger budget.

ChemistryLab._row_residual Method

_row_residual(Ae, n_e, be) -> Float64

Largest element-balance violation |Aₑnₑ − bₑ| relative to that element's own budget. Unlike a single global scale it cannot hide a small element behind a large one, which is what let a 0.465 mol sulfur violation report as 1.4e-2.

ChemistryLab._with_equilibrium_solver Method

_with_equilibrium_solver(kp::KineticsProblem, es) -> KineticsProblem

Return kp with its equilibrium solver replaced by es, or kp itself when es is nothing. Used by integrate so that a solver passed on the KineticsSolver reaches the ODE; a solver already set on the problem wins when both are given, and the mismatch is reported.

ChemistryLab.build_kinetics_ode Method

build_kinetics_ode(kp::KineticsProblem) -> Function

Build the ODE right-hand-side f!(du, u, p, t) implementing Leal et al. (2017).

State layout:

• u[1:n_be] = bₑ (element amounts in equilibrium partition)

• u[n_be+1 : n_be+n_nk] = nₖ (moles of kinetic species)

• u[n_be+n_nk+1 : n_be+n_nk+M] = ξ (extents of the M kinetic reactions)

• u[end] = T with a SemiAdiabaticCalorimeter, Q with an IsothermalCalorimeter, absent when there is no calorimeter

ODE equations (Leal 2017, Eq. 66):

• dnₖ/dt = νₖᵀ r

• dbₑ/dt = Aₑ νₑᵀ r

• dξ/dt = r

• dT/dt = (q̇ − φ(ΔT)) / Cp_total (semi-adiabatic)

• dQ/dt = q̇ (isothermal)

where nₑ = φ(bₑ) is the equilibrium re-speciation constraint.

ChemistryLab.build_kinetics_params Method

build_kinetics_params(kp::KineticsProblem; ϵ=1e-30) -> NamedTuple

Build the immutable parameter tuple p passed to the ODE function.

Key fields: T, P, ϵ, lna_fn, kin_rxns, species_index, n_initial_full, n_full, cp_fns, rates_buf, index ranges n_be, n_nk, idx_kinetic, idx_equilibrium, νe, νk, Ae.

ChemistryLab.build_u0 Method

build_u0(kp::KineticsProblem) -> Vector{Float64}

Build the initial ODE state vector.

Structure of u:

• Without re-speciation: u = [nₖ₁, …, nₖ_K, ξ₁, …, ξ_M]

• With re-speciation: u = [bₑ₁, …, bₑ_C, nₖ₁, …, nₖ_K, ξ₁, …, ξ_M]

• Semi-adiabatic adds T at the end: u = [..., T₀]

The extents of reaction ξ are carried alongside the kinetic moles, integrating dξ/dt = r. They are redundant with nₖ — the two satisfy nₖ = nₖ(0) + νₖᵀ ξ — but they are what makes the non-kinetic amounts available inside the residual, through n = n(0) + νᵀ ξ. Without them, a rate law gating on a species that is not itself kinetic would read a value frozen at t = 0. They also make reaction_extents and state_at exact rather than quadrature-limited.

The calorimeter's slot stays last and is addressed from the end of the vector, so it is unaffected by the presence of ξ.

ChemistryLab.respeciate! Method

respeciate!(p, u) -> Bool

Solve the equilibrium sub-problem once, and write the result into the running composition p.n_full. Returns true when a solve actually happened.

This is the second half of the operator-splitting step: the ODE advances the kinetic minerals with the speciation held frozen, then this function re-equilibrates the equilibrium partition under the element amounts the ODE has just produced.

The element amounts bₑ carried by the state vector are the constraint of that sub-problem (Leal et al. 2017, Eq. 54). solve conserves A·n, so what has to be handed to it is a composition whose element totals are exactly bₑ — here the previous speciation, projected onto bₑ through the pseudo-inverse of Aₑ. Handing over p.n_full unchanged, as an earlier version did, discards bₑ entirely and leaves the element balance to drift.

ChemistryLab.system_enthalpy Method

system_enthalpy(p, u, T) -> Float64

Σᵢ nᵢ ΔₐH⁰ᵢ(T) [J] over the composition of the ODE state u: the kinetic minerals read from u itself, the equilibrium partition from p.n_full.

Why the heat cannot come from the kinetic reactions here

heat_rate sums rᵢ (−ΔᵣH⁰ᵢ) over the KINETIC reactions, which is right when those reactions produce the hydrates. Under partial equilibrium they do not: they dissolve the anhydrous phases into ions, and the hydrates are precipitated by the Gibbs minimization, whose heat that sum cannot see. Measured on an ordinary Portland cement, counting only the dissolution put the semi-adiabatic temperature rise at 207 K where the test gives some tens of kelvin.

The enthalpy of the whole system has no such blind spot. It is a state function, so the heat released between two states at the same temperature is their difference — reactants, ions and hydrates all counted once, with no reaction stoichiometry to write down. This is Eq. (17)–(21) of Lavergne et al. (2018).

Kinetics solver

ChemistryLab._DEFAULT_KINETICS_SOLVER_FACTORY Constant

_DEFAULT_KINETICS_SOLVER_FACTORY

Ref holding a zero-argument factory () -> KineticsSolver for the default ODE solver used by integrate(kp) (no explicit solver argument).

Priority:

• KineticsOrdinaryDiffEqExt.__init__ registers Rodas5P() on load.

• Only one registration; extensions override it as needed.

Users can override with:

ChemistryLab._DEFAULT_KINETICS_SOLVER_FACTORY[] = () -> KineticsSolver(ode_solver=MyAlg())

ChemistryLab.KineticsSolver Type

struct KineticsSolver{ODE_S, ES}

Bundles the ODE algorithm and optional equilibrium solver used by integrate. Construct once, reuse across multiple KineticsProblem instances.

Fields

• ode_solver: any SciMLBase.AbstractODEAlgorithm (e.g. Rodas5P() from OrdinaryDiffEq), or :auto to let OrdinaryDiffEq choose — its default polyalgorithm detects stiffness at run time and switches. nothing keeps Rodas5P(), which stays the default because on the problems here :auto lands on the same stiff method and returns the same answers: switching would change nothing while removing a reproducible choice. Use nothing before OrdinaryDiffEq is loaded; an error will be raised at solve time.

• equilibrium_solver: optional EquilibriumSolver to re-equilibrate aqueous speciation at each ODE evaluation. When nothing, the kinetic minerals evolve without re-speciation (faster, less accurate).

• kwargs: keyword arguments forwarded to DifferentialEquations.solve (e.g. reltol, abstol, saveat, maxiters).

Examples

using OrdinaryDiffEq          # activates KineticsOrdinaryDiffEqExt
using Optimization, OptimizationIpopt   # needed for equilibrium_solver

es = EquilibriumSolver(cs, HKFActivityModel(), IpoptOptimizer())
ks = KineticsSolver(; ode_solver=Rodas5P(), equilibrium_solver=es,
                     reltol=1e-8, abstol=1e-10)
sol = integrate(kp, ks)

ChemistryLab.KineticsSolver Method

KineticsSolver(; ode_solver=nothing, equilibrium_solver=nothing, kwargs...) -> KineticsSolver

Construct a KineticsSolver.

kwargs are forwarded to DifferentialEquations.solve (e.g. reltol=1e-8, abstol=1e-10, saveat=0:60:3600).

ChemistryLab.integrate Method

integrate(kp::KineticsProblem; kwargs...) -> ODESolution

Shortcut that uses the default solver registered by _DEFAULT_KINETICS_SOLVER_FACTORY.

Requires OrdinaryDiffEq to be loaded (which sets the default to Rodas5P()).

using OrdinaryDiffEq
sol = integrate(kp)              # uses Rodas5P() by default
sol = integrate(kp; reltol=1e-6) # forward kwargs to the ODE solver

Calorimetry

ChemistryLab.AbstractCalorimeter Type

abstract type AbstractCalorimeter end

Base type for calorimeter models that can be coupled to a kinetics simulation.

Concrete subtypes:

• IsothermalCalorimeter: T = constant, tracks Q(t) = ∫q̇dt.

• SemiAdiabaticCalorimeter: variable-T cell (Lavergne et al. 2018).

ChemistryLab.IsothermalCalorimeter Type

struct IsothermalCalorimeter{T} <: AbstractCalorimeter

Isothermal calorimeter: temperature held constant at T [K]; cumulative heat Q(t) = ∫₀ᵗ q̇(τ) dτ [J] integrated as the trailing ODE state, with q̇ the heat of the kinetic reactions — see the caveat on cumulative_heat.

Examples

cal = IsothermalCalorimeter(298.15u"K")
kp = KineticsProblem(cs, reactions, state0, tspan; calorimeter = cal)
sol = integrate(kp, ks)
t, Q    = cumulative_heat(sol, cal)
t, qdot = heat_flow(sol, cal)

ChemistryLab.IsothermalCalorimeter Method

IsothermalCalorimeter(T) -> IsothermalCalorimeter

Plain Real → assumed SI [K]; Quantity → converted to K.

ChemistryLab.SemiAdiabaticCalorimeter Type

struct SemiAdiabaticCalorimeter{C, T, F} <: AbstractCalorimeter

Semi-adiabatic calorimeter following the Lavergne et al. (2018) energy balance:

where:

• q̇(t) [W] is the instantaneous heat-generation rate,

• φ(ΔT) [W] is the heat-loss function (e.g. linear L·ΔT or quadratic a·ΔT + b·ΔT²),

• Cp [J/K] is the fixed calorimeter heat capacity,

• Σᵢ nᵢ Cp°ᵢ(T) is the temperature- and mole-dependent sample heat capacity (computed from p.cp_fns at every ODE step when available).

Fields

• Cp: heat capacity of everything the ODE does not compute for itself — the vessel, the flask, an inert filler — in J/K, stored as a Quantity.

• heat_loss: callable φ(ΔT::Real) -> Real [W].

• T_env: ambient temperature [K] (stored as Quantity).

• T0: initial temperature [K] (stored as Quantity).

Cp must NOT include the sample

The denominator of the balance above is Cp + Σᵢ nᵢ Cp°ᵢ(T), and the second term is recomputed from the database at every step. Adding the sample's own heat capacity to Cp therefore counts it twice and understates the temperature rise — by a factor of about 1.75 on a cement paste at w/c = 0.4, where the paste contributes some 2.5 kJ/K against 0.9 kJ/K for the flask. Three scripts shipped with this package made exactly that mistake, and the field description above used to invite it by saying "calorimeter + sample". Pass the vessel, and let the solver add the paste.

Examples

# Linear heat loss — Newton cooling
cal = SemiAdiabaticCalorimeter(; Cp=4000.0u"J/K", T_env=293.15u"K", L=0.5u"W/K", T0=293.15u"K")

# Quadratic heat loss (Lavergne et al. 2018). `Cp` is the vessel alone;
# the paste's own heat capacity is added by the solver from the database.
cal = SemiAdiabaticCalorimeter(;
    Cp        = 900.0u"J/K",
    T_env     = 293.15u"K",
    heat_loss = ΔT -> 0.3*ΔT + 0.003*ΔT^2,
    T0        = 293.15u"K",
)

kp = KineticsProblem(cs, reactions, state0, tspan; calorimeter = cal)
sol = integrate(kp, ks)
t, T_vec = temperature_profile(sol, cal)
t, qdot  = heat_flow(sol, cal)

References

• Lavergne, F., Ben Fraj, A., Bayane, I. & Barthélémy, J.-F. (2018). Cement and Concrete Research 104, 37–60.

ChemistryLab.SemiAdiabaticCalorimeter Method

SemiAdiabaticCalorimeter(; Cp, T_env, T0, heat_loss=nothing, L=nothing)

Keyword constructor for SemiAdiabaticCalorimeter.

Exactly one of heat_loss or L must be provided:

• heat_loss: callable ΔT -> [W] (e.g. quadratic ΔT -> a*ΔT + b*ΔT^2).

• L: linear Newton cooling coefficient [W/K]. Sets heat_loss = ΔT -> L * ΔT.

All scalar fields accept plain Real (assumed SI) or Quantity:

• Cp → J/K; T_env → K; T0 → K; L → W/K.

ChemistryLab._total_enthalpy Method

_total_enthalpy(n_full, h_fns, T_K) -> Real

Total molar enthalpy H = Σᵢ nᵢ ΔₐH⁰ᵢ(T). Used by the DiscreteCallback in KineticsOrdinaryDiffEqExt.

ChemistryLab.cumulative_heat Method

cumulative_heat(sol, cal::IsothermalCalorimeter) -> (t, Q)

Cumulative heat Q(t) = ∫₀ᵗ q̇(τ) dτ [J], read off the ODE state the isothermal calorimeter adds.

This is the heat of the kinetic reactions only

q̇ is heat_rate, i.e. Σᵢ rᵢ(−ΔᵣH⁰ᵢ) over the kinetic reactions. That is the heat of hydration when those reactions produce the hydrates — the stoichiometric formulation. Under partial equilibrium they only dissolve the anhydrous phases into ions and the hydrates are precipitated by the Gibbs minimization, whose heat this sum cannot see; on an ordinary Portland cement the omission is worth hundreds of joules per gram. Use heat_release, which differences the enthalpy of certified speciations, whenever kp.equilibrium_solver !== nothing.

ChemistryLab.cumulative_heat Method

cumulative_heat(sol, cal::SemiAdiabaticCalorimeter) -> (t, Q)

Integrate the reconstructed heat-flow rate to obtain Q(t) [J].

ChemistryLab.extend_ode! Method

extend_ode!(du, u, p, n_kin, cal::SemiAdiabaticCalorimeter)

Append dT/dt = (q̇ − φ(ΔT)) / Cp_total(T, n) to the ODE right-hand side.

Cp_total = Cp + Σᵢ nᵢ Cp°ᵢ(T) is recomputed at every ODE step from p.cp_fns and p.n_full (Lavergne et al. 2018).

ChemistryLab.heat_flow Method

heat_flow(sol, cal::IsothermalCalorimeter) -> (t, qdot)

Instantaneous heat-generation rate q̇(t) [W], by differencing cumulative_heat — and carrying the same caveat about partial equilibrium.

ChemistryLab.heat_flow Method

heat_flow(sol, cal::SemiAdiabaticCalorimeter) -> (t, qdot)

Reconstruct q̇(t) [W] from the temperature ODE via the energy balance q̇ ≈ Cp × dT/dt + φ(T − T_env).

Note: uses the fixed cal.Cp (not the variable Cp_total) for this post-processing reconstruction.

ChemistryLab.heat_rate Method

heat_rate(kinetic_reactions, rates, T_K) -> Real

Compute the instantaneous heat generation rate [W = J/s]:

q̇ = Σᵢ rᵢ(t) × ΔHᵣ,ᵢ(T)

where rᵢ [mol/s] is the net rate of the i-th kinetic reaction (positive = dissolution/forward) and ΔHᵣ,ᵢ(T) [J/mol] is the enthalpy of reaction.

AD-compatible: the ΔₐH⁰ callables accept ForwardDiff.Dual T inputs.

ChemistryLab.temperature_profile Method

temperature_profile(sol, cal::SemiAdiabaticCalorimeter) -> (t, T)

Extract the temperature profile T(t) [K].

Post-processing a solution

ChemistryLab.degrees_of_hydration Method

degrees_of_hydration(sol, kp::KineticsProblem; times = sol.t) -> OrderedDict{String, Vector{Float64}}

Degree of reaction    of every kinetic species, keyed by species symbol.

Species with a zero initial amount are omitted rather than returned as NaN: a phase that was never present has no degree of reaction.

Examples

α = degrees_of_hydration(sol, kp)
α["C3S"][end]                       # final degree of hydration of alite

See also: mean_degree_of_hydration, state_at.

ChemistryLab.extent_residual Method

extent_residual(sol, kp::KineticsProblem, ξ; times = sol.t) -> Float64

Largest absolute discrepancy, in mol, between the kinetic-species moles carried by the integrator and those implied by the extents ξ, that is   .

The two are redundant by construction — nₖ integrates νₖᵀ r while ξ integrates r — so this measures how far the integrator has let them drift, and should sit at the solver's tolerance. A large value means the solve is too loose, not that the post-processing is inaccurate.

ChemistryLab.heat_release Method

heat_release(sol, kp; times = sol.t, reference = nothing, states = nothing)
    -> (times, Q, q̇)

Cumulative heat Q [J] and heat rate q̇ [W] released along the trajectory, from the certified speciations at times.

Q(t) = H(t₀) − H(t) with H = Σᵢ nᵢ ΔₐH⁰ᵢ(T), which is Eq. (17)–(21) of Lavergne et al. (2018): enthalpy is a state function, so its drop between two states at the same temperature is the heat given off, with reactants, ions and hydrates each counted once and no reaction stoichiometry to write down. Pass reference to measure from a state other than the first. Pass states when the certified replay at those instants is already in hand — it is the expensive part, and a caller that needs the compositions anyway should not pay for it twice.

Why this does not read the running composition

The obvious implementation — accumulate the enthalpy the integrator already has in hand at each accepted step — is wrong here, and not subtly. Under partial equilibrium the assemblage at each step comes from an in-run Gibbs minimization that is warm-started and not certified; the enthalpy is worth hundreds of kilojoules per hydrate, so an assemblage off by one phase moves the curve by more than the entire heat of hydration. Measured on an ordinary Portland cement, the in-run figures ran 12.7, 145, 1174, 936 and 631 J/g at 1 h, 6 h, 12 h, 1 d and 2 d — heat that increases, then decreases, which no calorimeter has ever seen.

speciated_states re-solves each instant and checks it against the KKT conditions, and that is what this reads.

Why the heat rate cannot come from the kinetic reactions either

heat_rate sums rᵢ(−ΔᵣH⁰ᵢ) over the kinetic reactions, which is right when those reactions produce the hydrates. Under partial equilibrium they only dissolve the anhydrous phases into ions; the hydrates are precipitated by the minimization, whose heat that sum cannot see. On the same cement it put a semi-adiabatic temperature rise at 207 K.

ChemistryLab.mean_degree_of_hydration Method

mean_degree_of_hydration(sol, kp::KineticsProblem; times = sol.t, weights = :mass) -> Vector{Float64}

Degree of reaction of the binder as a whole, averaged over the kinetic species.

weights selects the averaging:

• :mass (default) — weighted by the initial mass of each species, the convention used when a single ᾱ is quoted for a cement.

• :mole — weighted by initial moles.

Only species with a non-zero initial amount take part, consistently with degrees_of_hydration.

ChemistryLab.reaction_extents Method

reaction_extents(sol, kp::KineticsProblem; times = sol.t) -> Matrix{Float64}

Extent of reaction    for every kinetic reaction of kp, along the solution sol.

Returns a length(times) × M matrix, M being the number of kinetic reactions, in mol.

The extents are carried by the ODE state — the integrator advances dξ/dt = r alongside the kinetic moles — so this is a read of the solution, not a quadrature. It is exact to the solver's own tolerance at any instant, and costs nothing.

Arguments

• sol: the ODESolution returned by integrate.

• kp: the KineticsProblem that produced it.

• times: instants at which the extents are wanted. Any instants within tspan, in any order; the dense output is interpolated.

See also: state_at, degrees_of_hydration, extent_residual.

ChemistryLab.speciated_states Method

speciated_states(sol, kp::KineticsProblem; times = sol.t) -> Vector{ChemicalState}

The speciated compositions along a solution: the kinetic species read from the ODE state, and the equilibrium partition recovered by re-solving φ(bₑ) at each instant with the element totals the run carried.

This is what state_at deliberately does not do. That function returns the purely kinetic reconstruction n(0) + νᵀξ, because the redistribution performed by the equilibrium solve is not recoverable from the stoichiometry — for a cement that reconstruction is meaningless, putting every dissolved element in solution with not one hydrate. Nor is the composition left in the solver's own buffers usable: that buffer is rewritten at every right-hand-side evaluation, Jacobian differences and rejected steps included, so it is not the accepted composition at any particular time.

Without an equilibrium solver on kp there is nothing to replay, and this returns state_at for each instant.

The instants must be ascending

Each solve is warm-started from the previous one, and the guess is first capped at the element budget and projected back into {Aₑn = bₑ, n ≥ 0}. Both matter, and neither is optional:

• the warm start is the equilibrium of the previous bₑ, so once an element has been spent — the sulfate of an ordinary Portland cement, after the gypsum is gone — it demands more of it than now exists and the interior-point solve begins outside its own feasible set;

• solved instead from a cold guess, φ(bₑ) does not converge at all for a cement: the H⁺ component of bₑ reaches −14 mol, and a solve started from pure water returns no hydrates and a pore solution at pH 6 while the run itself computed 2.2 mol of C-S-H.

The chain is also walked up to the first instant requested, through a few earlier times of the trajectory whose compositions are discarded. That instant is the one with no predecessor of its own, and without the run-up it starts from the cast composition, which carries no active set at all: on a reference OPC its interior-point answer held 56 interior species where the answer has 25 — every candidate hydrate present, four of them at 1e-5 to 1e-6 mol — and neither the certifying solve nor its continuation recovers from that, both inheriting the start. With the run-up that instant is certified, at stationarity 9.1e-13 and element balance 6.1e-15 mol.

Examples

states = speciated_states(sol, kp; times = [1, 7, 28] .* 86400.0)
pH(states[end])                      # pore solution at 28 days
volume_fractions(states[end], groups; reference = state0)

See also: state_at, reaction_extents, volume_fractions.

ChemistryLab.state_at Method

state_at(sol, kp::KineticsProblem, t; ξ = nothing) -> ChemicalState

Reconstruct the full ChemicalState of the system at time t from a kinetics solution.

Every species — kinetic and non-kinetic alike — is obtained from the same stoichiometric replay   , so the returned state satisfies element conservation exactly. Temperature is taken from the solution when a SemiAdiabaticCalorimeter is in use, and from the problem otherwise.

The extents are read from the ODE state, so this is exact to the solver's tolerance. Pass a precomputed extent matrix as ξ — any matrix whose last row holds ξ(t) — to avoid re-interpolating the solution when sweeping many instants.

Re-speciation is not replayed

When the run used an equilibrium solver, the aqueous partition was re-speciated at every accepted step and that redistribution is not recoverable from the stoichiometry alone. The returned state is then the purely kinetic reconstruction; call equilibrate on it to recover the speciated composition.

Examples

ξ = reaction_extents(sol, kp)                      # once
states = [state_at(sol, kp, t; ξ = ξ[i:i, :]) for (i, t) in enumerate(sol.t)]

See also: reaction_extents, degrees_of_hydration, volume_fractions.

The implicit kinetic step

One fully implicit problem per step, with the reaction extents as unknowns of the same Gibbs minimization. See Which route: two ways to advance in time.

ChemistryLab.KineticStepSolver Type

KineticStepSolver(system, model, reactions; opts...)

A solver for one fully implicit kinetic step, following Leal et al. (2017).

reactions are KineticReaction objects. Every species not carried by a declared reaction is at equilibrium; every declared reaction is kinetic and must have a rate law, which KineticReaction already enforces.

Requires OptimaSolver.

Example

kss = KineticStepSolver(cs, DiluteSolutionModel(), [kr_calcite])
st1 = kinetic_step(kss, st0, 60.0u"s")

See also kinetic_step, integrate for the ODEProblem route.

ChemistryLab._certificate_error Method

_certificate_error(c) -> Float64

One number ranking how far a step is from satisfying its own conditions: stationarity, element balance, the worst supersaturation among absent phases, and the residual of the rate equation. Inf when no certificate was produced.

Used to choose between the pinned and free formulations of a step. A ranking is needed rather than the boolean optimal, because two answers may straddle the tolerance and then the choice is decided by rounding — the same code choosing differently on another machine. Ranked, the gap is fifteen orders of magnitude.

ChemistryLab._kinetic_step_eliminated Method

_kinetic_step_eliminated(kss, state, Δt; ...)

The :species step: the pinned species are computed from the extents rather than solved for, and the extents come from a Newton on nr unknowns around a plain equilibrium over the species that remain.

F(Δξ) = Δξ − Δt·r(n(Δξ)), where n(Δξ) puts the imposed amounts where the stoichiometry says and lets the minimization place the rest under the budget those amounts leave. Backward Euler, as kinetic_step is, since r is read at the composition the step produces.

ChemistryLab._normalize_certificate Method

_normalize_certificate(c) -> NamedTuple

One shape for the certificate whichever route produced it.

kinetic_step fills its certificate reference from OptimaSolver's kkt_certificate on the :reactions route and from optimality_certificate on the :species one, and those name the element-balance field differently — feasibility against balance. Returning two shapes from one keyword is a trap for the caller, so both are mapped onto the names optimality_certificate uses.

ChemistryLab._primary_symbols Method

_primary_symbols(system) -> Vector{String}

The symbols of the system's primary components, to rebuild a reduced system over the same basis.

ChemistryLab._reactivity_matrix Method

_reactivity_matrix(reactions, system, kinetic_species) -> Matrix

The matrix K of the reactivity constraint Kᵀn − Δξ = ξ₀: one column per reaction, carrying that reaction's stoichiometry restricted to the kinetic partition — the species whose amounts the reactions alone decide.

The restriction is the modeling decision of the whole scheme. A rate law says how fast a mineral reacts; it says nothing about where matter released into solution ends up, which is the equilibrium problem's business. So an aqueous species must be left OUT of K: it re-speciates, and including it makes Kᵀn a functional of the products, which stops pinning the mineral. Measured on calcite dissolving at 1 µmol/s for 1000 s with the full stoichiometry in K, the amount came out 5.7e-5 mol below n₀ − kΔt, the constraint having absorbed the carbonate's speciation into the extent. The answer solved its own equations exactly and meant something else.

kinetic_species = :auto takes the partition to be every non-aqueous participant of a declared reaction. That covers both styles of model:

• dissolution into solution — Cal → Ca²⁺ + CO₃²⁻. Only Cal is kinetic, so K = [−1], M = KᵀK = 1, Δξ = Δt·r, and the released matter speciates freely under element conservation.

• solid to solid — C₃S + H₂O → C-S-H + CH, the form (Lavergne et al., 2018) uses. Every solid product is kinetic too, so the assemblage is imposed by the stoichiometry rather than found by minimizing G. That is a different model, deliberately, and it is the one to write when the hydrate assemblage is part of what you are prescribing.

Pass an explicit list of species (symbols or indices) to override the partition.

Several reactions may share a mineral — the whole point of a reaction-centric scheme, and what (Lavergne et al., 2018) needs for C₃A reacting through both an ettringite and a monosulfoaluminate pathway. Their columns must be linearly independent, though: two reactions differing only in products the partition excludes would give the same row twice, Kᵀn could not tell their extents apart, and the linear system would be singular. That is checked here, with the offending reactions named.

ChemistryLab._step_rate_closure Method

_step_rate_closure(kss, des, T, P, t, n0) -> (x, params) -> Vector

The reaction rates at composition x, in mol/s, one per declared reaction.

Built once per step so the StateViews and the initial amounts are not rebuilt at every residual evaluation — the residual is called once per Newton iteration and once per column of the difference quotient, so nr + m + … times per step.

ChemistryLab.kinetic_step Method

kinetic_step(kss, state, Δt; t = 0.0, ϵ = 1e-16, parameters = nothing)
    -> ChemicalState

One implicit kinetic step of duration Δt from state.

t is the absolute time handed to rate laws that depend on it. Pass a Ref as parameters to receive the reaction extents Δξ the step found.

warm_start equilibrates the starting GUESS when the system carries solid solutions, leaving the component totals untouched. Keep it on: a mixing phase is admitted by a tangent-plane test, and from a composition where the phase is absent that admission fails — measured, a cold start on C₃S dissolving into a C-S-H solution left one end-member at 2.7e-9 with a stationarity residual of 6.5.

The rate laws are evaluated at the end-of-step composition, which is what makes the step implicit and stable on a stiff system. Δt is the caller's choice; kinetic_step_adaptive chooses it from a local error estimate instead.

One step far beyond the relaxation time can land on the other root

The step is a nonlinear problem, and for a rate law that vanishes at equilibrium it has a second solution: the composition with the mineral wholly dissolved. That one satisfies the element balance and the reactivity row while violating Δξ − Δt·M·r(n) = 0 by the entire extent, so the certificate refuses it — pass certificate and check optimal before trusting a step much larger than 1/k.

Which root the Newton finds is a property of the build, not of the chemistry. Measured on calcite under r = k(1 − Ω) with k = 10⁻⁵ mol/s: steps of 10⁴ s and 10⁶ s converge to the right root on Julia 1.12 and to the other one on 1.13.0-rc4, reported uncertified in both readings. Steps well inside the relaxation time certify on either.

kinetic_step_adaptive is the remedy and is build-independent: it refuses an uncertified step and halves until one certifies, reaching the equilibrium values to eight digits in seven steps on the same case.

ChemistryLab.kinetic_step_adaptive Method

kinetic_step_adaptive(kss, state, Δt; reltol, abstol, ...) -> (state, Δt_used, Δt_next)

One kinetic step with its length chosen from a local error estimate.

Takes one step of Δt and two of Δt/2, compares the extents, and accepts the finer pair when the estimated error is within tolerance — halving and retrying otherwise. Returns the accepted state, the step actually taken, and a suggestion for the next one.

The estimate is Richardson's for a first-order method: ξ_fine − ξ_coarse is the error of the coarse step to leading order. The error controlled is therefore on the reaction extents, which is what a kinetic step advances; the composition follows from them through an exact equilibrium solve, so it carries no separate error of its own.

Δt is the step to attempt, max_halvings bounds the retries, and the tolerance is abstol + reltol·sⱼ per reaction, where sⱼ is the amount the extent acts on — Σᵢ |Kᵢⱼ| nᵢ, the kinetic minerals of reaction j at the start of the step.

Scaling by the extent instead, reltol·|Δξⱼ|, does not work, and the reason is worth stating because it is the obvious thing to write. Δξ ∝ Δt, so that tolerance vanishes with the step while the equilibrium solve's own noise does not: the measured error then behaves as noise / (reltol·Δt) and grows as the step shrinks. Measured, the controller halved to the floor without advancing and the final answer got worse as the tolerance was tightened — 2.5e-5 at 1e-3 against 9.2e-5 at 1e-5. An ODE integrator's reltol multiplies the solution, not the increment, for exactly this reason.

Example

st, dt_used, dt_next = kinetic_step_adaptive(kss, st0, 3600.0u"s"; reltol = 1e-4)

See also kinetic_step for a step of a length you choose, and integrate for the ODEProblem route, which brings SciML's own step-size control and its choice of stiff solvers.