Condition Profiles

For an introduction into Kinetica's arbitrary simulation condition framework, see the tutorial on Arbitrary Simulation Conditions. This page aims to explain the underlying implementation of Kinetica's variable condition profiles, such that users can easily extend Kinetica by adding their own profiles to represent new macroscopic state variable changes.

As noted previously, variable condition profiles can take on one of two types of definition which effect how they are represented and interacted with inside Kinetica's solvers:

• Directly variable condition profiles (subtypes of Kinetica.AbstractDirectProfile) encode a condition $X$ as a continuous function of time $X(t)$. The result of this function is then computed directly whenever it is needed at a given time $t$.
• Gradient-variable condition profiles (subtypes of Kinetica.AbstractGradientProfile) encode a condition $X$ within a gradient with respect to time $\frac{dX}{dt}$. These condition profiles must be integrated with time using DifferentialEquations.jl to yield an interpolable ODESolution which can be queried to produce values of the condition at time $t$.

Both types of variable condition profile are simple to implement within Kinetica, provided a few criteria are fulfilled:

• Both are implemented as structs with supertypes as above that indicate how Kinetica should handle them. These structs must have a number of fields, but are otherwise entirely customisable for whatever data needs to be held. These usually come with outer constructors that create variables that are needed within the profile's condition function, as well as binding a condition function to one of the struct's fields.
• A condition function (usually called _f_[struct name] for direct condition profiles or _grad_[struct_name] for gradient condition profiles) must exist such that it can be bound to a condition profile struct. This function takes two arguments - the time t and the profile that it will ultimately be bound to.
• A method of Kinetica.create_discrete_tstops! must be created for the new profile struct which tells Kinetica where rate constants should be updated in discrete rate constant update simulations.

We'll demonstrate this by implementing both directly variable and gradient-variable versions of a simple sinusoidal condition profile.

Directly Variable Profile Implementation

When a condition profile can be represented exactly as a simple function of time, it is usually easiest to implement as a directly variable condition profile. In the case of our sinusoidal profile, this function of time will be a sinusoid:

\[X(t)=A\sin\left( 2\pi f t \right)\]

where $A$ is the amplitude and $f$ is the number of oscillations per second.

We will begin by creating the struct for our condition profile, which we'll call SinusoidDirectProfile. Aside from the parameters above, we'll also need to implement fields for some basics which Kinetica expects in every directly variable profile:

• f: The function which will be bound in the struct's outer constructor.
• X_start: The initial value of the condition profile.
• t_end: The time at which the condition profile should stop varying. This can either be provided by the user or calculated within the outer constructor.
• tstops: Time points at which ODE solvers should ensure to stop at and recalculate. This array is useful for handling discontinuities in condition profiles, and is usually modified by Kinetica.create_discrete_tstops!.
• sol: Profile solution over the requested timespan, usually created internally by calling Kinetica.solve_variable_condition!(::Kinetica.AbstractDirectProfile, ::ODESimulationParams).

For this profile, we will also implement a field X_end to signify the value of the condition once it stops varying. In some profiles, this is more easily calculated by the profile itself, e.g. as a function of time when t_end is given. Conversely, in other profiles, it can be convenient to only provide X_end and let t_end be calculated instead. As a sinusoidal profile is periodic and may reach the same value multiple times, users will provide t_end directly and we will calculate and save the value of X_end in the outer constructor.

An implementation of our condition profile's struct may therefore look something like:

using Kinetica

mutable struct SinusoidDirectProfile{uType, tType} <: Kinetica.AbstractDirectProfile
    f::Function
    A::uType
    freq::uType
    X_start::uType
    X_end::uType
    t_end::tType
    tstops::Vector{tType}
    sol
end

We've used a parametric type here to ensure that all values that may affect the value of our condition are using consistent numeric types.

With this definition complete, we can now define our profile's condition function and outer constructor. The condition function will take a time t and an instance of the SinusoidDirectProfile. The outer constructor will take a minimal set of parameters needed to populate and calculate the entire condition profile, and create an instance of the profile with the condition function bound to the f field:

function _f_SinusoidDirectProfile(t, profile::SinusoidDirectProfile)
    return typeof(profile.X_start)(
        ((t <= 0.0) * profile.X_start) +
        ((t > 0.0 && t <= profile.t_end) * (profile.X_start + profile.A*sin(2*pi*profile.freq*t))) +
        ((t > profile.t_end) * profile.X_end)
    )
end

function SinusoidDirectProfile(;
    A::uType,
    freq::uType,
    X_start::uType,
    t_end::tType
) where {uType <: AbstractFloat, tType <: AbstractFloat}

    X_end = X_start + A*sin(2*pi*freq*t_end)
    tstops = [t_end]

    return SinusoidDirectProfile(_f_SinusoidDirectProfile, A, freq, X_start, X_end, t_end, tstops, nothing)
end

How we've defined _f_SinusoidDirectProfile(t, profile) here may look a little strange. This is actually a way of writing if-else statements that is compatible with Symbolics.jl, and is equivalent to writing

function _f_SinusoidDirectProfile(t, profile::SinusoidDirectProfile)
    if t <= 0.0
        return profile.X_start
    elseif t > 0.0 && t <= profile.t_end
        return profile.X_start + profile.A*sin(2*pi*profile.freq*t)
    else
        return profile.X_end
    end
end

However, Symbolics isn't compatible with this form of if-else statement so instead we use a bit of boolean algebra to embed the same meaning into our condition function. We've therefore assumed that simulation time will start at t=0.0, before which time the condition will be X_start. After t_end, the condition profile will stop varying and be set to the value of the sinusoid at this time. In between these two times, we take the instantaneous value of the sinusoid instead.

Finally, we need to define a method of Kinetica.create_discrete_tstops! for our new profile. This function creates time stopping points for discrete rate update simulations. Some profiles contain periods of time where their condition is stationary, during which rate constant updates are unnecessary, so their method of this function can avoid creating time points during these periods. However, since the gradient of a sinusoid is only ever momentarily stationary, this optimisation is not necessary. Our implementation of this method is therefore as follows:

function Kinetica.create_discrete_tstops!(profile::SinusoidDirectProfile, ts_update::AbstractFloat)
    if ts_update > profile.t_end
        throw(ArgumentError("Error defining tstops, `ts_update` is too large."))
    end
    profile.tstops = create_savepoints(0.0, profile.t_end, ts_update)
    return
end

This is simply checking that the rate update timestep ts_update is not larger than the total duration of our condition profile, then setting the tstops of our profile to an evenly spaced array of time points between t=0.0 and t=t_end. create_savepoints is just a Kinetica utility function for creating evenly spaced ranges of values while avoiding floating point rounding errors.

With this all implemented, our directly variable condition profile is now ready for use! Let's start by constructing a ConditionSet for a discrete rate update simulation with a sinusoidally-varying temperature:

conditions = ConditionSet(Dict(
    :T => SinusoidDirectProfile(;
        A = 75.0,
        freq = 0.25,
        X_start = 900.0,
        t_end = 20.0
    )),
    ts_update=tconvert(0.1, "ms", "s")
)

Because we've made SinusoidDirectProfile a subtype of Kinetica.AbstractDirectProfile, Kinetica's ConditionSet has done a bunch of work in the background and set this profile up for use within Symbolics.jl by registering the generated SinusoidDirectProfile.f function (which is really just a reference to _f_SinusoidDirectProfile) into Symbolics' computation graph. This means it's ready for use within kinetic simulations, just like that!

# Load in existing CRN...
# pars = ODESimulationParams(...)
# calc = PrecalculatedArrheniusCalculator(...)

solvemethod = VariableODESolve(pars, conditions, calc)
res = solve_network(solvemethod, sd, rd)

[ Info:  - Calculating variable condition profiles.
[ Info:  - Filtering reactions...
[ Info:    - Removed 0 filtered reactions from network
[ Info:  - Performing calculator-specific network setup.
[ Info:  - Removing low-rate reactions
[ Info:    - Low rate cutoff: automatic (cutoff = 5.0e-11)
[ Info:    - Removed 0 low-rate reactions from network.
[ Info:  - Setting up ReactionSystem
[ Info:  - Created ReactionSystem
[ Info:  - Pre-calculating rate constants at discrete time intervals.
[ Info:  - Created callback for discrete rate constant updates.
[ Info:  - Formulating ODEProblem
[ Info:    - Sparse? true
[ Info:    - Analytic Jacobian? true
[ Info:  - Setting up integrator...
[ Info:  - Solving network...
[ Info:  - Solved.

using Plots

p1 = plot(res)
p2 = conditionsplot(res, :T)
plot(p1, p2, layout=(2, 1))

As you can see, the condition profile has been seamlessly integrated into Kinetica. Within solve_network, a new method of Kinetica.solve_variable_condition!(::Kinetica.AbstractDirectProfile, ::ODESimulationParams) was created and run for our condition profile, which generated an interpolable DiffEqArray that is now available in the solved profile's sol field:

profile = get_profile(res.conditions, :T)
profile.sol

t: 1001-element Vector{Float64}:
  0.0
  0.02
  0.04
  0.06
  0.08
  0.1
  0.12
  0.14
  0.16
  0.18
  ⋮
 19.84
 19.86
 19.88
 19.9
 19.92
 19.94
 19.96
 19.98
 20.0
u: 1001-element Vector{Vector{Float64}}:
 [900.0]
 [902.3558069308597]
 [904.7092889646985]
 [907.0581234988886]
 [909.3999925173229]
 [911.7325848780173]
 [914.0535985939293]
 [916.3607431047407]
 [918.6517415373642]
 [920.9243329529422]
 ⋮
 [881.348258462636]
 [883.6392568952592]
 [885.9464014060704]
 [888.2674151219825]
 [890.6000074826774]
 [892.9418765011116]
 [895.2907110353015]
 [897.6441930691404]
 [899.9999999999999]

Gradient-Variable Profile Implementation

Implementing a gradient-variable condition profile is a very similar procedure to the above directly variable profile - we'll still define a new struct that needs some specific fields to be compatible with the rest of Kinetica, we'll still make an outer constructor for this struct that creates a specialised function which will be automatically registered within Symbolics.jl, and we'll need a new method for Kinetica.create_discrete_tstops!.

Unlike the directly variable profile, we'll need to implement a function for our sinusoid's gradient with respect to time. In the case of this profile, a gradient-based definition is not required since we already have an exact function for mapping our sinusiodal condition to simulation time. However, it turns out that the gradient-based representation actually helps us include an additional parameter to describe the phase of the sinusoid $\varphi$:

\[\frac{dX}{dt}=2\pi fA \cos \left( 2\pi f t + φ \right)\]

This would've been inconvenient to include into our directly variable form, as the user's input value of X_start would've needed to be changed within the outer constructor to take the starting phase into account. However, a gradient-based representation lends itself naturally to this addition.

The struct for this condition profile, which we'll call SinusoidGradientProfile, requires the same fields as its directly variable counterpart. The exception to this is the f field, which previously held the actual time-dependent condition function. In gradient-based profiles, this is replaced with a field called grad, which holds the gradient function. We also need to subtype Kinetica.AbstractGradientProfile this time, and we'll add a field for our starting phase $\varphi$. Our new struct might therefore look something like this:

mutable struct SinusoidGradientProfile{uType, tType} <: Kinetica.AbstractGradientProfile
    grad::Function
    A::uType
    freq::uType
    φ::uType
    X_start::uType
    t_end::tType
    tstops::Vector{tType}
    sol
end

The condition function and outer constructor for this profile will also be similar to their directly variable counterparts, again with an extra parameter φ added:

function _grad_SinusoidGradientProfile(t, profile::SinusoidGradientProfile)
    return typeof(profile.X_start)(
        ((t <= 0.0) * 0.0) +
        ((t > 0.0 && t <= profile.t_end) * (2*pi*profile.freq*profile.A*cos(2*pi*profile.freq*t + profile.φ))) +
        ((t > profile.t_end) * 0.0)
    )
end

function SinusoidGradientProfile(;
    A::uType,
    freq::uType,
    φ::uType,
    X_start::uType,
    t_end::tType
) where {uType <: AbstractFloat, tType <: AbstractFloat}

    tstops = [t_end]
    return SinusoidGradientProfile(_grad_SinusoidGradientProfile, A, freq, φ, X_start, t_end, tstops, nothing)
end

The gradient function we've defined here describes a function that is stationary (gradient of zero) before t=0.0 and after t=t_end, but that varies sinusoidally in between. This sinusoidal variation now also takes into account the starting phase $\varphi$.

The implementation of Kinetica.create_discrete_tstops! can be the same as it was for the directly variable profile, as we are describing the same underlying function:

function Kinetica.create_discrete_tstops!(profile::SinusoidGradientProfile, ts_update::AbstractFloat)
    if ts_update > profile.t_end
        throw(ArgumentError("Error defining tstops, `ts_update` is too large."))
    end
    profile.tstops = create_savepoints(0.0, profile.t_end, ts_update)
    return
end

And with that, we're done implementing another condition profile! Let's see it in action:

conditions = ConditionSet(Dict(
    :T => SinusoidGradientProfile(;
        A = 75.0,
        freq = 0.25,
        φ = pi/2,
        X_start = 975.0,
        t_end = 20.0
    )),
    ts_update=tconvert(0.1, "ms", "s")
)

solvemethod = VariableODESolve(pars, conditions, calc)
res = solve_network(solvemethod, sd, rd)

[ Info:  - Calculating variable condition profiles.
[ Info:  - Filtering reactions...
[ Info:    - Removed 0 filtered reactions from network
[ Info:  - Performing calculator-specific network setup.
[ Info:  - Removing low-rate reactions
[ Info:    - Low rate cutoff: automatic (cutoff = 5.0e-11)
[ Info:    - Removed 0 low-rate reactions from network.
[ Info:  - Setting up ReactionSystem
[ Info:  - Created ReactionSystem
[ Info:  - Pre-calculating rate constants at discrete time intervals.
[ Info:  - Created callback for discrete rate constant updates.
[ Info:  - Formulating ODEProblem
[ Info:    - Sparse? true
[ Info:    - Analytic Jacobian? true
[ Info:  - Setting up integrator...
[ Info:  - Solving network...
[ Info:  - Solved.

This time, because we're using a gradient-based profile, when solve_network is called Kinetica goes away and assembles a ModelingToolkit.jl ODESystem and solves it with respect to time, just like it does when performing the kinetic simulation of the CRN immediately afterwards. As such, the sol field of our profile now contains a full ODESolution:

profile = get_profile(res.conditions, :T)
profile.sol

retcode: Success
Interpolation: 1st order linear
t: 201001-element Vector{Float64}:
  0.0
  0.0001
  0.0002
  0.0003
  0.0004
  0.0005
  0.0006
  0.0007
  0.0008
  0.0009
  ⋮
 19.9993
 19.9994
 19.9995
 19.9996
 19.9997
 19.9998
 19.9999
 20.0
 20.0
u: 201001-element Vector{Vector{Float64}}:
 [975.0]
 [974.9999990747245]
 [974.9999962988983]
 [974.9999916725213]
 [974.9999851955938]
 [974.9999768681158]
 [974.9999666900876]
 [974.9999546615094]
 [974.9999407823814]
 [974.999925052704]
 ⋮
 [974.9999546614908]
 [974.9999666900691]
 [974.9999768680973]
 [974.9999851955753]
 [974.9999916725028]
 [974.9999962988798]
 [974.999999074706]
 [974.9999999999815]
 [974.9999999999815]

We can also look at the simulation results from using our gradient-based profile. Notice how both the concentrations and the temperature have changed since we introduced the phase parameter - now temperature starts from the top of the sinusiod rather than the middle:

p1 = plot(res)
p2 = conditionsplot(res, :T)
plot(p1, p2, layout=(2, 1))