mkin4py

mkin(microkinetics) 4 py(thon)

Copyright © 2015 - Gabriel Sabença Gusmão

A general package for linearly defining and solving microkinetic catalytic systems.

Description

A microkinetic package translated from Linearized Microkinetic Catalytic System Solver

By the author:

Gabriel S. Gusmão gusmaogabriels@gmail.com under Dr. Phillip Christopher christopher@engr.ucr.edu advisement.

For detailed information, refer to code comments or associated publication. Gusmão, G. S. & Christopher, P., A general and robust approach for defining and solving microkinetic catalytic systems. AIChE J. 00, (2014).; http://dx.doi.org/10.1002/aic.14627

The 17-Step Ethylene Epoxidation by Stegelmann et al. has been used as example. Stegelmann, C., Schiødt, N. C., Campbell, C. T. & Stoltze, P. Microkinetic modeling of ethylene oxidation over silver. J. Catal. 221, 630–649 (2004).

1. Set-up the environment conditions (temperature, pressure, gas constant)
2. Create a MK (microkinetic model object)
   • Define its dimensions: number of reactants (rows) and elementary reactions (columns) involved in the stoichiometry matrix, and parse the rows that refer to free-species (non-adsorbed)
   • Parse the stoichiometry matrix (must be of size number of reactants × number of elementary reactions)
   • Set the kinetic parameters: Activation Energies and Pre-exponential factors (must be of the size of the involved elementary reactions)
   • Set the fixed concentration of free-species (molar fraction in non-adsorbed phase)
   • Parse the string-labels of involved species (array of size of number of species)
3. Solve the nonlinear steady-state system using successive linearized corrections.
   • The original Newton-type method integrates its correction flow with RK4; it does not integrate physical time.
   • The inner linear system uses QMR (quasi-minimal residual) by default, with a dense fallback. A direct dense solver is also selectable.

Defaults are in mkin4py.solver.params.convergence_params. Override them per call, for example mkin4py.solver.solve.rk4(linear_solver='dense', criteria=1e-8).

Features

Linearization

The project makes use of explicit routines for the calculation of the MK model derivatives

- *Jacobian*: Available as standard.
- *Hessian*: Used in the convex two-step method (`rk4(param=2)`; details in the aforementioned reference)

JAX execution

Model arrays and returned coverages/rates are JAX arrays. Numerical solver attempts are JIT-compiled; model configuration and the restart driver run in Python. The original analytic derivatives remain available, and explicit array kernels support JAX transformations; see JAX compatibility.

On the way

1. Additional linear-system solvers beyond the selectable QMR and dense methods.
2. Evolutionary methods for the definition of best convergence parameters for stiff problems (when TOF`s are close to the machine precision)

Instructions

See JAX execution, local CLI, and compatibility details.

• Installation

  Requires Python 3.11 or later and JAX 0.4.38 or later. Pip installs JAX as a dependency. Install or upgrade to the JAX release from PyPI:

  python -m pip install --upgrade mkin4py
  

  Or install the current source using the original installation route:

  python -m pip install --upgrade https://github.com/gusmaogabriels/mkin4py/zipball/master
  

  Version 2.0.0 includes JAX execution, steady-state sensitivities and the local CLI. Version 1.0 remains available for legacy environments.

• Local CLI

  mkin4py methods --json
  mkin4py example --json > model.json
  mkin4py solve model.json --json
  mkin4py solve model.json --linear-solver dense --json
  

  The CLI reads and solves the model locally and enables 64-bit precision by default. The Python example below enables it explicitly for stiff kinetics.

• Example: Stoltze's 17-Step Ethylene Epoxidation MK system

  import mkin4py
  import jax
  import jax.numpy as np
  jax.config.update('jax_enable_x64', True)
  
  # Environment Conditions
  T = 500; #K
  P = 1; #bar
  gas_constant = 8.31456e-3 # Gas Constant - kJ/(mol×K)
  
  # Set the environment conditions
  mkin4py.environment.set_temperature(T)
  mkin4py.environment.set_gas_constant(gas_constant)
  mkin4py.environment.set_pressure(P)
  
  # Stoichiometric Matrix
  ms = [
  [-1, 1, 0, 0,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\
  [ 0, 0, 0, 0, 0, 0,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\
  [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\
  [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\
  [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,-1, 0, 0],\
  [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,-1],\
  [-1, 1,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,-1,-1, 1, 1,-1, 0, 0, 1,-1, 0, 0, 1,-1, 1,-1],\
  [ 1,-1,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\
  [ 0, 0, 2,-2,-2, 2,-1, 1, 1,-1,-1, 1, 0, 0, 0, 0, 1,-1,-6, 6, 0, 0,-1, 1,-1, 1,-5, 5, 1,-1, 0, 0, 0, 0],\
  [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4,-4, 0, 0, 0, 0, 1,-1, 3,-3,-2, 2, 0, 0, 0, 0],\
  [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,-1, 0, 0,-1, 1],\
  [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0,-1, 1, 0, 0],\
  [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\
  [ 0, 0, 0, 0, 2,-2, 0, 0,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\
  [ 0, 0, 0, 0, 0, 0, 1,-1,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\
  [ 0, 0, 0, 0, 0, 0, 0, 0, 1,-1, 0, 0,-1, 1,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\
  [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,-1, 1,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\
  [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,-1,-1, 1,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\
  [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,-1,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0],\
  [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,-1,-1, 1, 0, 0, 0, 0, 0, 0],\
  ];
  
  nreac = [0, 1] # Reactant Rows in mS
  nprod = [2, 3, 4, 5] # Product Rows in mS
  stoichs = np.array(nreac + nprod) # Reactants and Products are not under PSSA
  mkin4py.mkmodel.create(len(ms),len(ms[0]),stoichs) # Initialize the model
  mkin4py.mkmodel.set_ms(ms) # Set the stoichiometry matrix
  
  # Species labels 
  splabels = ['O2','C2H4','C2H4O','CH3CHO','CO2','H2O','*','O2*','O*','OH*',\
  'H2O*','CO2*','C2H4*','O∙O*','C2H4∙O*','CH2CH2O∙O*','C2H4O∙O*','CH3CHO∙O*',\
  'CH2CHOH∙O*','CH2CHO∙O*']
  
  mkin4py.mkmodel.set_splabels(splabels) # Set species labels
  
  # Pre-exponential Factors of Eelementary Reactions (1/s)
  va =[2.71e5, 1.1e12, 4.0e12, 8.0e14, 2.0e7, 1.3e15, 7.2e7,\
  2.2e11, 9.0e14, 5.3e14, 1.95e8, 4.8e12, 1.13e13, 2.11e12,\
  9.0e12, 4.5e10, 2.9e13, 2.6e9, 2.0e20, 5.3e13, 7.2e7, 2.2e11,\
  4.0e11, 3.1e14, 2.6e13, 1.3e9, 1.0e20, 5.5e13, 1.4e10, 1.0e11,\
  3.6e14, 1.0e8, 5.9e14, 1.4e9]
  
  # Activation Barriers for Elementary ReactionS (kJ/mol)
  vea = [5.7000, 47.3000, 75.0000, 157.5000, 20.0000, 96.9000, 0, 37.1000, 112.0000,\
  183.3000, 0, 39.1000, 95.0000, 93.5000, 95.0000, 204.3000, 41.9000, 4.4000,\
  11.0000, 791.6000, 0, 30.1000, 32.0000, 42.8000, 86.0000, 106.1000, 0, 906.6000,\
  65.6000, 50.0000, 38.9000, 0, 46.6000, 0]
  
  # Set the kinetic parameters
  mkin4py.mkmodel.set_kinetic_params(np.array(va,ndmin=2).T,np.array(vea,ndmin=2).T)
  
  y = [0.5, 0.5, 0, 0, 0, 0] # Reactants and Products Initial Fraction
  mkin4py.mkmodel.set_concentrations(y) # Set the *free*-species concentrations
  

• Evaluation:

  sol = mkin4py.solver.solve.rk4() # RK4 Newton correction flow; QMR inner linear solve
  # Outputs
  print ('...')
  print (sol['msg'], 'time: ', sol['time'])
  print ('Coverage')
  print (sol['coverage'])
  print ('Rates')
  print (sol['rates'])
  

• Output:

  The original output below is retained as a numerical reference. Elapsed time depends on the machine and includes any JAX compilation incurred by the call; the historical time shown here is not a JAX benchmark. Convergence requires the maximum absolute surface-species rate to be at most criteria (default 1e-8); tiny residuals may differ from the reference. Check sol['success'] and sol['status'] when handling results programmatically.

  ...
  Convergence achieved time:  2.25999999046
  Coverage
  [[  5.00000000e-01]
  [  5.00000000e-01]
  [  0.00000000e+00]
  [  0.00000000e+00]
  [  0.00000000e+00]
  [  0.00000000e+00]
  [  4.39342950e-01]
  [  1.19743307e-03]
  [  1.07992516e-01]
  [  1.10447591e-01]
  [  2.99730332e-09]
  [  7.70711567e-10]
  [  1.00256049e-01]
  [  9.78269843e-02]
  [  1.32727193e-01]
  [  9.64368428e-03]
  [  3.28419897e-08]
  [  4.59118359e-13]
  [  5.65562897e-04]
  [  1.12150752e-15]]
  Rates
  [[ -4.24252190e+01]
  [ -2.49532633e+01]
  [  1.29732696e+01]
  [  5.58723011e-04]
  [  2.39588699e+01]
  [  2.39588699e+01]
  [  0.00000000e+00]
  [  3.65929509e-13]
  [  0.00000000e+00]
  [  4.32857086e-11]
  [  2.76796815e-16]
  [  2.87485591e-11]
  [  0.00000000e+00]
  [ -2.27373675e-13]
  [ -4.65661287e-10]
  [  9.86479981e-16]
  [  0.00000000e+00]
  [ -1.60491195e-13]
  [  4.65661287e-10]
  [ -1.42115222e-11]]