Kinetics API

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=8.31446)
k(; T = 298.15)   # → 1e-5  (rate constant at reference temperature)

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 initialisation, 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

PK_PARAMS_C2S :: NamedTuple

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

Original paper values (K₁=0.95, K₂=0.0005, K₃=0.0024 d⁻¹). Activation energy from Schindler & Folliard (2005). Reference temperature: 293.15 K (20 °C).

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.

Original paper values (K₁=0.082, K₂=0.00024, K₃=0.0024 d⁻¹). Activation energy from Schindler & Folliard (2005). Reference temperature: 293.15 K (20 °C).

PK_PARAMS_C3S :: NamedTuple

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

Original paper values (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))

PK_PARAMS_C4AF :: NamedTuple

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

Original paper values (K₁=0.165, K₂=0.0015, K₃=0.0024 d⁻¹). Activation energy from Schindler & Folliard (2005). Reference temperature: 293.15 K (20 °C).

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(PK_PARAMS_C3S, "C3S");

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

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)])

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

Construct a RateMechanism with no catalyst contributions.

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)

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

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)",
))

arrhenius_rate_constant(k₀, Ea; T_ref=298.15, R_gas=8.31446261815324) -> 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 8.31446261815324).

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)

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

Build the Parrot & Killoh (1984) cement clinker hydration rate as a KineticFunc.

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):

The rate [mol/s] is n_initial × Aₜ × min(max(r_NG, r_I), r_D) where ξ = α / α_max is the normalised 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

julia> pk = parrot_killoh(PK_PARAMS_C3S, "C3S");

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

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

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

julia> 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.

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

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

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

where stoich[i] is the stoichiometric coefficient (positive for products, negative for reactants), lna[i] is the log-activity of species i, and ΔₐG⁰overT[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⁰overT: 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).

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).

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

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.

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²

FixedSurfaceArea(A) -> FixedSurfaceArea

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

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)

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)

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.

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

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

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)

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).
• ϵ: regularisation 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)

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.

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.

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

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

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.
• ϵ: regularisation 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.

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, [T]]
• With re-speciation: u = [bₑ₁, …, bₑ_C, nₖ₁, …, nₖ_K, [T]]

where bₑ are the element amounts in the equilibrium partition, nₖ are the moles of kinetic species, and T is the temperature (semi-adiabatic only).

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.

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))

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[end] = T (semi-adiabatic only)

ODE equations (Leal 2017, Eq. 66):

• dnₖ/dt = νₖᵀ r
• dbₑ/dt = Aₑ νₑᵀ r
• dT/dt = (q̇ − φ(ΔT)) / Cp_total (semi-adiabatic)

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

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.

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

Build the initial ODE state vector.

Structure of u:

• Without re-speciation: u = [nₖ₁, …, nₖ_K]
• With re-speciation: u = [bₑ₁, …, bₑ_C, nₖ₁, …, nₖ_K]
• Semi-adiabatic adds T at the end: u = [..., T₀]

_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())

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). 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)

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).

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

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).

struct IsothermalCalorimeter{T} <: AbstractCalorimeter

Isothermal calorimeter: temperature held constant at T [K]; cumulative heat Q(t) = ∫₀ᵗ q̇(τ) dτ [J] tracked as an extra ODE state.

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)

IsothermalCalorimeter(T) -> IsothermalCalorimeter

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

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

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

\[\frac{dT}{dt} = \frac{\dot{q}(t) - \varphi(T - T_{\rm env})}{C_p + \sum_i n_i C^\circ_{p,i}(T)}\]

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 calorimeter + sample [J/K] (stored as Quantity).
• heat_loss: callable φ(ΔT::Real) -> Real [W].
• T_env: ambient temperature [K] (stored as Quantity).
• T0: initial temperature [K] (stored as Quantity).

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)
cal = SemiAdiabaticCalorimeter(;
    Cp        = 3449.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.

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.

_total_enthalpy(n_full, h_fns, T_K) -> Real

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

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

Extract cumulative heat Q(t) = ∫₀ᵗ q̇(τ) dτ [J].

When total-enthalpy tracking data are available in sol.prob.p.saved_H, returns Q(t) = H(0) − H(t) (captures both kinetic and equilibrium heat).

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

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

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).

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

Extract instantaneous heat-generation rate q̇(t) [W] from an ODE solution.

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.

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.

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

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