Skip to content

Instantly share code, notes, and snippets.

@speth
Last active March 14, 2023 18:39
Show Gist options
  • Save speth/6f9f47dba815ae45597e42169a88b6d6 to your computer and use it in GitHub Desktop.
Save speth/6f9f47dba815ae45597e42169a88b6d6 to your computer and use it in GitHub Desktop.
Test problem for reactor Jacobian with surface species
import cantera as ct
import numpy as np
np.set_printoptions(precision=2, linewidth=144)
surf = ct.Interface('dummysurf2.yml', 'surf')
surf2 = ct.Interface('dummysurf2.yml', 'surf2')
gas = surf.adjacent['gas']
gas.TPX = 1000, 2e5, 'A:0.6, B:0.3, C:0.2, D:0.1'
surf.coverages = 'A(S):0.1, B(S):0.2, C(S):0.3, D(S):0.2, (S):0.2'
surf2.coverages = 'A(S):0.1, D(S):0.2, (S):0.2'
surf.derivative_settings = {'skip-coverage-dependence': True, 'skip-electrochemistry': True}
surf2.derivative_settings = {'skip-coverage-dependence': True, 'skip-electrochemistry': True}
r = ct.IdealGasMoleReactor(gas)
r.volume = 3
rsurf2 = ct.ReactorSurface(surf2, r, A=9e-4)
rsurf1 = ct.ReactorSurface(surf, r, A=5e-4)
net = ct.ReactorNet([r])
net.step()
print('Analytical Jacobian:')
print(r.jacobian)
print()
print('Finite difference Jacobian')
print(r.finite_difference_jacobian)
units: {length: cm, quantity: mol, activation-energy: J/mol}
phases:
- name: gas
thermo: ideal-gas
species: [A, B, C, D]
kinetics: gas
reactions: [gas-reactions]
state:
T: 300.0
P: 1.01325e+05
- name: surf
thermo: ideal-surface
adjacent-phases: [gas]
species: [A(S), B(S), C(S), D(S), (S)]
kinetics: surface
reactions: [surf-reactions]
state:
T: 900.0
coverages: {(S): 1.0}
site-density: 2.7063e-09
species:
- name: A
composition: {C: 1}
thermo:
model: NASA7
temperature-ranges: [300.0, 1000.0, 3000.0]
data:
- [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]
- name: B
composition: {C: 1}
thermo:
model: NASA7
temperature-ranges: [300.0, 1000.0, 3000.0]
data:
- [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]
- name: C
composition: {C: 1}
thermo:
model: NASA7
temperature-ranges: [300.0, 1000.0, 3000.0]
data:
- [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]
- name: D
composition: {C: 1}
thermo:
model: NASA7
temperature-ranges: [300.0, 1000.0, 3000.0]
data:
- [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]
- name: A(S)
composition: {C: 1}
thermo:
model: NASA7
temperature-ranges: [300.0, 1000.0, 3000.0]
data:
- [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]
- name: B(S)
composition: {C: 1}
thermo:
model: NASA7
temperature-ranges: [300.0, 1000.0, 3000.0]
data:
- [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]
- name: C(S)
composition: {C: 1}
thermo:
model: NASA7
temperature-ranges: [300.0, 1000.0, 3000.0]
data:
- [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]
- name: D(S)
composition: {C: 1}
thermo:
model: NASA7
temperature-ranges: [300.0, 1000.0, 3000.0]
data:
- [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]
- name: (S)
composition: {}
thermo:
model: NASA7
temperature-ranges: [300.0, 1000.0, 3000.0]
data:
- [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]
gas-reactions:
- equation: A + B = C + D
rate-constant: {A: 1200, b: 0, Ea: 0}
- equation: A + C = B + D
rate-constant: {A: 800, b: 0, Ea: 0}
surf-reactions:
- equation: A + (S) => A(S)
rate-constant: {A: 5e3, b: 0, Ea: 0}
- equation: B + (S) => B(S)
rate-constant: {A: 2e4, b: 0, Ea: 0}
- equation: A(S) => A + (S)
rate-constant: {A: 500, b: 0, Ea: 0}
- equation: B(S) => B + (S)
rate-constant: {A: 800, b: 0, Ea: 0}
- equation: C(S) => C + (S)
rate-constant: {A: 500, b: 0, Ea: 0}
- equation: D(S) => D + (S)
rate-constant: {A: 400, b: 0, Ea: 0}
- equation: A(S) + B(S) => C(S) + D(S)
rate-constant: {A: 1200, b: 0, Ea: 0}
Analytical Jacobian:
[[ 1.14e+05 0.00e+00 4.88e+08 -8.34e+08 5.06e+08 -3.51e+09 -8.91e+14 -2.48e+14 1.99e+14 -1.50e+15 -1.32e+14 -3.03e+14 -5.54e+12 2.71e+14]
[ 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00]
[-1.72e+00 0.00e+00 -1.04e+03 1.83e+04 -1.61e+02 6.07e+04 3.00e+09 0.00e+00 -1.62e+09 5.00e+09 0.00e+00 0.00e+00 0.00e+00 -6.01e+08]
[ 1.70e+00 0.00e+00 -4.01e+02 -2.11e+04 1.76e+03 -5.75e+04 0.00e+00 0.00e+00 0.00e+00 0.00e+00 8.00e+08 0.00e+00 0.00e+00 -1.20e+09]
[-1.70e+00 0.00e+00 4.01e+02 2.11e+04 -1.76e+03 5.75e+04 0.00e+00 0.00e+00 0.00e+00 3.24e+09 0.00e+00 1.12e+09 0.00e+00 0.00e+00]
[ 1.72e+00 0.00e+00 1.04e+03 -1.83e+04 1.61e+02 -6.07e+04 7.54e+09 4.36e+09 0.00e+00 0.00e+00 0.00e+00 0.00e+00 4.00e+08 0.00e+00]
[ 0.00e+00 0.00e+00 1.31e-01 0.00e+00 0.00e+00 0.00e+00 -1.12e+10 -3.76e+09 1.62e+09 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00]
[ 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 -6.90e+09 -4.36e+09 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00]
[ 0.00e+00 0.00e+00 -1.31e-01 0.00e+00 0.00e+00 0.00e+00 1.81e+10 8.11e+09 -1.62e+09 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00]
[ 0.00e+00 0.00e+00 4.54e-02 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 -8.32e+09 -4.03e+07 -1.07e+09 0.00e+00 6.01e+08]
[ 0.00e+00 0.00e+00 0.00e+00 1.81e-01 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 -8.12e+07 -8.40e+08 0.00e+00 0.00e+00 1.20e+09]
[ 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 -3.16e+09 4.03e+07 -1.12e+09 0.00e+00 0.00e+00]
[ 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 8.12e+07 4.03e+07 0.00e+00 -4.00e+08 0.00e+00]
[ 0.00e+00 0.00e+00 -4.54e-02 -1.81e-01 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 1.15e+10 8.00e+08 2.20e+09 4.00e+08 -1.80e+09]]
Finite difference Jacobian
[[ 1.14e+05 7.55e+06 3.47e+08 -9.75e+08 3.65e+08 -3.65e+09 -4.47e+14 -1.33e+14 1.17e+14 -7.57e+14 -7.77e+13 -1.37e+14 -1.22e+13 1.60e+14]
[ 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00]
[-1.72e+00 -1.09e+02 -1.04e+03 1.83e+04 -1.61e+02 6.07e+04 3.00e+09 0.00e+00 -1.62e+09 5.00e+09 0.00e+00 0.00e+00 0.00e+00 -6.01e+08]
[ 1.70e+00 1.20e+02 -4.01e+02 -2.11e+04 1.76e+03 -5.75e+04 0.00e+00 0.00e+00 0.00e+00 0.00e+00 7.99e+08 0.00e+00 0.00e+00 -1.20e+09]
[-1.70e+00 -1.20e+02 4.01e+02 2.11e+04 -1.76e+03 5.75e+04 0.00e+00 0.00e+00 0.00e+00 3.24e+09 0.00e+00 1.12e+09 0.00e+00 0.00e+00]
[ 1.72e+00 1.09e+02 1.04e+03 -1.83e+04 1.61e+02 -6.07e+04 7.54e+09 4.36e+09 0.00e+00 0.00e+00 0.00e+00 0.00e+00 4.00e+08 0.00e+00]
[ 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 -1.12e+10 -3.76e+09 1.62e+09 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00]
[ 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 -6.90e+09 -4.36e+09 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00]
[ 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 1.81e+10 8.11e+09 -1.62e+09 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00]
[ 0.00e+00 -5.46e-04 4.54e-02 -6.45e-09 0.00e+00 -1.94e-08 0.00e+00 0.00e+00 0.00e+00 -8.32e+09 -4.03e+07 -1.07e+09 0.00e+00 6.01e+08]
[ 0.00e+00 -1.09e-03 -8.07e-10 1.81e-01 0.00e+00 -4.84e-09 0.00e+00 0.00e+00 0.00e+00 -8.12e+07 -8.40e+08 1.44e+01 2.15e+01 1.20e+09]
[ 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 -3.16e+09 4.03e+07 -1.12e+09 0.00e+00 0.00e+00]
[ 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 8.12e+07 4.03e+07 -7.18e+00 -4.00e+08 0.00e+00]
[ 0.00e+00 1.64e-03 -4.54e-02 -1.81e-01 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 1.15e+10 8.00e+08 2.20e+09 4.00e+08 -1.80e+09]]
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment