Comparing Simulation Modes
This tutorial explains the two simulation modes in MKMCXX 3 from the numerical point of view.
Microkinetic models are written as a system of ordinary differential equations
where y contains the gas-phase amounts and surface coverages.
There are two natural ways to solve this system:
transient_solve: integrate the ODE system forward in timesteady_state_solve: solve the algebraic steady-state system
These are the two simulation modes exposed by MKMCXX.
The Two Numerical Viewpoints
transient_solve
This mode performs time integration.
Use it when you want:
- the physical time evolution of the system
- a robust route toward a stationary state
- a diagnostic view of whether the system is still changing
The run stops when either:
t_endis reached, or- the maximum absolute derivative drops below
steady_state_threshold
steady_state_solve
This mode does not integrate in time. It solves the nonlinear algebraic system
dy/dt = 0 directly with KINSOL.
Use it when you want:
- the stationary solution itself
- no transient trajectory
- a faster direct solve when the nonlinear problem is well conditioned
Its main limitation is that a direct nonlinear solve can be more sensitive to the initial guess than time integration.
Example 1: A Realistic Surface Mechanism
We first use the CO-methanation-inspired mechanism from
examples/simulation-modes/input.mkx.
This is a useful comparison case because it is realistic enough to show the different numerical behavior of the two modes.
Part 1: Time Integration
The runnable example file is:
examples/simulation-modes/input.mkx
It runs a transient sweep at:
475 K500 K525 K550 K
Run it with:
<path-to-executable>/mkmcxx -i examples/simulation-modes/input.mkx -f
The transient run completes successfully for all four temperatures. In the
convergence report, the Stop column is final_time, not steady_state:
Job ID SimID End t Stop Path Threshold Max|dy/dt|
------------------------------------------------------------------------------------------------
1 T=475 2.000e+07 final_time - 1.000e-08 1.683e-05
2 T=500 1.000e+07 final_time - 1.000e-08 3.234e-05
3 T=525 8.000e+06 final_time - 1.000e-08 5.468e-05
4 T=550 6.000e+06 final_time - 1.000e-08 8.323e-05
This means:
- the integration is numerically stable
- the system is evolving in a sensible way
- the requested time horizon was reached
- the chosen steady-state threshold was not yet met
So the transient mode works well here, but it is still reporting a state that is only approximately stationary.
Part 2: Direct Steady-State Solve
Now keep the same components and reactions blocks, but replace the
simulation block with:
simulation {
steady_ref @ steady_state_solve(
model=thermal,
solver={method=kinsol_newton, fnormtol=1e-10, scsteptol=1e-10, max_iter=400, strategy=none},
conditions=[
{T=475},
{T=500},
{T=525},
{T=550}
]
)
}
When this is run, all four temperatures converge to steady state. However, the
important detail is the Path column:
Job ID SimID End t Stop Path Threshold Max|dy/dt|
------------------------------------------------------------------------------------------------
1 T=475 n/a steady_state warm_start 1.000e-10 2.564e-11
2 T=500 n/a steady_state warm_start 1.000e-10 5.642e-11
3 T=525 n/a steady_state warm_start 1.000e-10 3.017e-11
4 T=550 n/a steady_state warm_start 1.000e-10 4.268e-11
The key observation is that Path = warm_start for all four temperatures.
This tells us that the algebraic steady-state solve did not converge cleanly from the raw default initial state using the direct KINSOL route alone. MKMCXX therefore triggered its internal fallback mechanism: it first preconditioned the state with a transient warm-up and then retried the steady-state solve.
Pedagogically, this is useful because it shows a realistic limitation of the algebraic formulation:
- the transient method is robust on this mechanism
- the steady-state method still reaches the same stationary chemistry
- but the nonlinear solve needs help to get there
The Path entry is exactly where you see that.
Comparing The Results
Even though the direct steady-state route needed the fallback path, the final stationary results are very close to the transient results.
For the methane-forming step CH3* + H* => CH4 + 2*, the rates are:
| Temperature | Transient rate | Steady-state rate | Relative difference |
|---|---|---|---|
475 K |
5.609174e-06 |
5.609536e-06 |
6.45e-05 |
500 K |
1.078062e-05 |
1.078295e-05 |
2.16e-04 |
525 K |
1.822803e-05 |
1.823359e-05 |
3.05e-04 |
550 K |
2.774318e-05 |
2.774168e-05 |
5.41e-05 |
The maximum absolute concentration difference is also small:
| Temperature | Max absolute concentration difference |
|---|---|
475 K |
1.62e-05 |
500 K |
5.09e-05 |
525 K |
6.60e-05 |
550 K |
1.26e-05 |
So both modes are describing the same surface state within numerical
approximation. The remaining difference is expected, because the transient runs
stopped at final_time rather than at the requested steady-state threshold.
Example 2: A Simple System With No Fallback
The previous mechanism is intentionally realistic. Now let us look at a much simpler example:
components {
A* {phase=surface, init=0.2},
B* {phase=surface, init=0.8}
}
reactions {
{A*}=>{B*} @ ArrheniusDefault(Vf=1.0, Vb=1.0, Eaf=0.0, Eab=0.0)
}
This is just a reversible interconversion between two surface species.
Transient Solve
simulation {
transient_ref @ transient_solve(
model=thermal,
solver={method=cvode_bdf, atol=1e-12, rtol=1e-12, t_end=100.0, steady_state_threshold=1e-10},
conditions=[{T=600}]
)
}
The transient run reaches steady state directly:
Job ID SimID End t Stop Path Threshold Max|dy/dt|
------------------------------------------------------------------------------------------------
1 T=600 1.259e+01 steady_state - 1.000e-10 1.084e-11
Steady-State Solve
simulation {
steady_ref @ steady_state_solve(
model=thermal,
solver={method=kinsol_newton, fnormtol=1e-12, scsteptol=1e-12, max_iter=50, strategy=line_search},
conditions=[{T=600}]
)
}
Here the steady-state solve also succeeds immediately, and now the Path
column reads direct:
Job ID SimID End t Stop Path Threshold Max|dy/dt|
------------------------------------------------------------------------------------------------
1 T=600 n/a steady_state direct 1.000e-12 0.000e+00
This is the simplest case of the algebraic approach working exactly as one would hope: no fallback, no transient preconditioning, just a direct nonlinear solve.
Comparing The Surface Coverages
Both runs report the same final coverages:
| Mode | A* |
B* |
|---|---|---|
transient_solve |
0.5 |
0.5 |
steady_state_solve |
0.5 |
0.5 |
This is the cleanest possible demonstration that the two modes are solving the same underlying model, just from different numerical viewpoints.
What To Take Away
The main lesson is:
transient_solveis the time-integration view of the ODE systemsteady_state_solveis the algebraic view of the same system
On simple problems, both modes converge directly and give the same answer. On more nonlinear surface mechanisms, the transient mode is often the more robust route, while the steady-state mode may require an internal fallback path before it reaches the same stationary chemistry.
So when reading MKMCXX output:
- use
Stopto see whether the requested end condition was met - use
Pathto see how the steady-state solver actually converged - compare the final concentrations or rates when you want to verify that the two modes agree