Skip to content

Commit 4008277

Browse files
sevyharrissevyharris
committed
add SurfaceReactor jacobian test with finite differences
1 parent 3b1f125 commit 4008277

1 file changed

Lines changed: 221 additions & 0 deletions

File tree

‎test/rmgpy/solver/solverSurfaceTest.py‎

Lines changed: 221 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -39,9 +39,151 @@
3939
from rmgpy.solver.surface import SurfaceReactor
4040
from rmgpy.species import Species
4141
from rmgpy.thermo import ThermoData, NASA, NASAPolynomial
42+
from rmgpy.solver.termination import TerminationTime
43+
from rmgpy.rmg.settings import ModelSettings, SimulatorSettings
4244

4345

46+
def get_i_thing(thing, thing_list):
47+
for i in range(len(thing_list)):
48+
if thing.is_isomorphic(thing_list[i]):
49+
return i
50+
return -1
51+
4452
class SurfaceReactorTest:
53+
def setup_class(self):
54+
# Define a simple surface mechanism that we can use for testing
55+
# This is faster than loading in chemkin
56+
self.species_list = [
57+
Species(
58+
molecule=[Molecule().from_smiles("[Ar]")],
59+
thermo=NASA(polynomials=[NASAPolynomial(coeffs=[2.5,0,0,0,0,-745.375,4.37967], Tmin=(200,'K'), Tmax=(1000,'K')), NASAPolynomial(coeffs=[2.5,0,0,0,0,-745.375,4.37967], Tmin=(1000,'K'), Tmax=(6000,'K'))], Tmin=(200,'K'), Tmax=(6000,'K'), comment="""Thermo library: primaryThermoLibrary"""),
60+
),
61+
Species(
62+
molecule=[Molecule().from_smiles("[Ne]")],
63+
thermo=NASA(polynomials=[NASAPolynomial(coeffs=[2.5,0,0,0,0,-745.375,3.35532], Tmin=(200,'K'), Tmax=(1000,'K')), NASAPolynomial(coeffs=[2.5,0,0,0,0,-745.375,3.35532], Tmin=(1000,'K'), Tmax=(6000,'K'))], Tmin=(200,'K'), Tmax=(6000,'K'), comment="""Thermo library: primaryThermoLibrary"""),
64+
),
65+
Species(
66+
molecule=[Molecule().from_smiles("N#N")],
67+
thermo=NASA(polynomials=[NASAPolynomial(coeffs=[3.53101,-0.000123661,-5.02999e-07,2.43531e-09,-1.40881e-12,-1046.98,2.96747], Tmin=(200,'K'), Tmax=(1000,'K')), NASAPolynomial(coeffs=[2.95258,0.0013969,-4.92632e-07,7.8601e-11,-4.60755e-15,-923.949,5.87189], Tmin=(1000,'K'), Tmax=(6000,'K'))], Tmin=(200,'K'), Tmax=(6000,'K'), comment="""Thermo library: primaryThermoLibrary"""),
68+
),
69+
Species(
70+
molecule=[Molecule().from_smiles("C")],
71+
thermo=NASA(polynomials=[NASAPolynomial(coeffs=[4.20542,-0.00535559,2.51124e-05,-2.13763e-08,5.97526e-12,-10161.9,-0.921283], Tmin=(100,'K'), Tmax=(1084.12,'K')), NASAPolynomial(coeffs=[0.908259,0.0114541,-4.57174e-06,8.29193e-10,-5.66316e-14,-9719.97,13.9931], Tmin=(1084.12,'K'), Tmax=(5000,'K'))], Tmin=(100,'K'), Tmax=(5000,'K'), comment="""Thermo library: primaryThermoLibrary"""),
72+
),
73+
Species(
74+
molecule=[Molecule().from_smiles("[O][O]")],
75+
thermo=NASA(polynomials=[NASAPolynomial(coeffs=[3.53732,-0.00121572,5.3162e-06,-4.89446e-09,1.45846e-12,-1038.59,4.68368], Tmin=(100,'K'), Tmax=(1074.55,'K')), NASAPolynomial(coeffs=[3.15382,0.00167804,-7.69974e-07,1.51275e-10,-1.08782e-14,-1040.82,6.16756], Tmin=(1074.55,'K'), Tmax=(5000,'K'))], Tmin=(100,'K'), Tmax=(5000,'K'), comment="""Thermo library: primaryThermoLibrary"""),
76+
),
77+
Species(
78+
molecule=[Molecule().from_smiles("O")],
79+
thermo=NASA(polynomials=[NASAPolynomial(coeffs=[4.05764,-0.000787929,2.90875e-06,-1.47516e-09,2.12833e-13,-30281.6,-0.311362], Tmin=(100,'K'), Tmax=(1130.23,'K')), NASAPolynomial(coeffs=[2.84325,0.00275108,-7.81028e-07,1.07243e-10,-5.79385e-15,-29958.6,5.9104], Tmin=(1130.23,'K'), Tmax=(5000,'K'))], Tmin=(100,'K'), Tmax=(5000,'K'), comment="""Thermo library: primaryThermoLibrary"""),
80+
),
81+
Species(
82+
molecule=[Molecule().from_smiles("[H][H]")],
83+
thermo=NASA(polynomials=[NASAPolynomial(coeffs=[3.43536,0.000212712,-2.78629e-07,3.4027e-10,-7.76039e-14,-1031.36,-3.90842], Tmin=(100,'K'), Tmax=(1959.07,'K')), NASAPolynomial(coeffs=[2.78819,0.000587616,1.59022e-07,-5.52763e-11,4.34328e-15,-596.156,0.112618], Tmin=(1959.07,'K'), Tmax=(5000,'K'))], Tmin=(100,'K'), Tmax=(5000,'K'), comment="""Thermo library: primaryThermoLibrary"""),
84+
),
85+
Species(
86+
molecule=[Molecule().from_smiles("*")],
87+
thermo=NASA(polynomials=[NASAPolynomial(coeffs=[0,0,0,0,0,0,0], Tmin=(298,'K'), Tmax=(1000,'K')), NASAPolynomial(coeffs=[0,0,0,0,0,0,0], Tmin=(1000,'K'), Tmax=(3000,'K'))], Tmin=(298,'K'), Tmax=(3000,'K'), comment="""Thermo library: surfaceThermoPt111"""),
88+
),
89+
Species(
90+
molecule=[Molecule().from_smiles("[H]*")],
91+
thermo=NASA(polynomials=[NASAPolynomial(coeffs=[-2.0757,0.0173581,-2.60921e-05,1.89282e-08,-5.38836e-12,-3166.19,8.15362], Tmin=(298,'K'), Tmax=(1000,'K')), NASAPolynomial(coeffs=[2.72248,-0.00106817,1.98654e-06,-1.12048e-09,2.09812e-13,-4218.24,-15.3207], Tmin=(1000,'K'), Tmax=(2000,'K'))], Tmin=(298,'K'), Tmax=(2000,'K'), comment="""Thermo library: surfaceThermoPt111"""),
92+
),
93+
Species(
94+
molecule=[Molecule().from_smiles("[O]=*")],
95+
thermo=NASA(polynomials=[NASAPolynomial(coeffs=[-0.294476,0.0144163,-2.61323e-05,2.19006e-08,-6.98019e-12,-16461.9,-0.199446], Tmin=(298,'K'), Tmax=(1000,'K')), NASAPolynomial(coeffs=[2.90245,-0.000338584,6.43373e-07,-3.66327e-10,6.90094e-14,-17049.7,-15.256], Tmin=(1000,'K'), Tmax=(2000,'K'))], Tmin=(298,'K'), Tmax=(2000,'K'), comment="""Thermo library: surfaceThermoPt111"""),
96+
),
97+
Species(
98+
molecule=[Molecule().from_smiles("[CH3]*")],
99+
thermo=NASA(polynomials=[NASAPolynomial(coeffs=[-0.0444549,0.0194368,-1.91029e-05,1.11269e-08,-2.73736e-12,-6388.04,-0.173376], Tmin=(298,'K'), Tmax=(1000,'K')), NASAPolynomial(coeffs=[8.65705,-0.00790308,1.401e-05,-7.40016e-09,1.31517e-12,-8635.99,-44.3353], Tmin=(1000,'K'), Tmax=(2000,'K'))], Tmin=(298,'K'), Tmax=(2000,'K'), comment="""Thermo library: surfaceThermoPt111"""),
100+
),
101+
Species(
102+
molecule=[Molecule().from_smiles("[H][H].*")],
103+
thermo=NASA(polynomials=[NASAPolynomial(coeffs=[3.86406,0.000753456,-1.65571e-06,1.55223e-09,-4.46782e-13,-1689.28,-8.85807], Tmin=(298,'K'), Tmax=(1000,'K')), NASAPolynomial(coeffs=[4.0688,-0.000495807,6.59234e-07,-1.72598e-10,7.62965e-15,-1700.7,-9.71918], Tmin=(1000,'K'), Tmax=(2000,'K'))], Tmin=(298,'K'), Tmax=(2000,'K'), comment="""Thermo library: surfaceThermoPt111"""),
104+
),
105+
Species(
106+
molecule=[Molecule().from_smiles("[OH]*")],
107+
thermo=NASA(polynomials=[NASAPolynomial(coeffs=[1.4236,0.015783,-2.91659e-05,2.50433e-08,-8.04088e-12,-18999.3,-3.1523], Tmin=(298,'K'), Tmax=(1000,'K')), NASAPolynomial(coeffs=[5.03574,-0.00134422,2.25916e-06,-1.08548e-09,1.77876e-13,-19634.4,-20.0345], Tmin=(1000,'K'), Tmax=(2000,'K'))], Tmin=(298,'K'), Tmax=(2000,'K'), comment="""Thermo library: surfaceThermoPt111"""),
108+
),
109+
Species(
110+
molecule=[Molecule().from_smiles("O.*")],
111+
thermo=NASA(polynomials=[NASAPolynomial(coeffs=[2.72971,0.00871052,-1.29132e-05,1.07295e-08,-3.39434e-12,-32612.7,-6.0448], Tmin=(298,'K'), Tmax=(1000,'K')), NASAPolynomial(coeffs=[5.85496,-0.00328847,5.56991e-06,-2.73008e-09,4.55898e-13,-33304.6,-21.3518], Tmin=(1000,'K'), Tmax=(2000,'K'))], Tmin=(298,'K'), Tmax=(2000,'K'), comment="""Thermo library: surfaceThermoPt111"""),
112+
),
113+
Species(
114+
molecule=[Molecule().from_smiles("[CH2]=*")],
115+
thermo=NASA(polynomials=[NASAPolynomial(coeffs=[-2.23007,0.0292223,-4.33155e-05,3.31428e-08,-9.96471e-12,-222.256,8.30173], Tmin=(298,'K'), Tmax=(1000,'K')), NASAPolynomial(coeffs=[6.8346,-0.00514926,9.15491e-06,-4.84917e-09,8.63767e-13,-2258.98,-36.2215], Tmin=(1000,'K'), Tmax=(2000,'K'))], Tmin=(298,'K'), Tmax=(2000,'K'), comment="""Thermo library: surfaceThermoPt111"""),
116+
),
117+
Species(
118+
molecule=[Molecule().from_smiles("[CH]#*")],
119+
thermo=NASA(polynomials=[NASAPolynomial(coeffs=[-2.66805,0.0290693,-4.82654e-05,3.87589e-08,-1.19749e-11,-2918.16,9.72941], Tmin=(298,'K'), Tmax=(1000,'K')), NASAPolynomial(coeffs=[4.9043,-0.00263865,4.71729e-06,-2.51267e-09,4.49659e-13,-4464.41,-26.7108], Tmin=(1000,'K'), Tmax=(2000,'K'))], Tmin=(298,'K'), Tmax=(2000,'K'), comment="""Thermo library: surfaceThermoPt111"""),
120+
),
121+
Species(
122+
molecule=[Molecule().from_smiles("[C]$*")],
123+
thermo=NASA(polynomials=[NASAPolynomial(coeffs=[-1.94351,0.0197767,-3.36337e-05,2.69027e-08,-8.27959e-12,7000.57,7.1747], Tmin=(298,'K'), Tmax=(1000,'K')), NASAPolynomial(coeffs=[2.81347,-0.000693952,1.30308e-06,-7.387e-10,1.38796e-13,6060.03,-15.5738], Tmin=(1000,'K'), Tmax=(2000,'K'))], Tmin=(298,'K'), Tmax=(2000,'K'), comment="""Thermo library: surfaceThermoPt111"""),
124+
)
125+
]
126+
X = self.species_list[get_i_thing(Species(molecule=[Molecule().from_smiles("*")]), self.species_list)]
127+
O2 = self.species_list[get_i_thing(Species(molecule=[Molecule().from_smiles("[O][O]")]), self.species_list)]
128+
OX = self.species_list[get_i_thing(Species(molecule=[Molecule().from_smiles("O=*")]), self.species_list)]
129+
H2 = self.species_list[get_i_thing(Species(molecule=[Molecule().from_smiles("[H][H]")]), self.species_list)]
130+
HX = self.species_list[get_i_thing(Species(molecule=[Molecule().from_smiles("[H]*")]), self.species_list)]
131+
CH4 = self.species_list[get_i_thing(Species(molecule=[Molecule().from_smiles("C")]), self.species_list)]
132+
CH3X = self.species_list[get_i_thing(Species(molecule=[Molecule().from_smiles("C*")]), self.species_list)]
133+
OHX = self.species_list[get_i_thing(Species(molecule=[Molecule().from_smiles("[OH]*")]), self.species_list)]
134+
H2O = self.species_list[get_i_thing(Species(molecule=[Molecule().from_smiles("O")]), self.species_list)]
135+
H2OX = self.species_list[get_i_thing(Species(molecule=[Molecule().from_smiles("O.*")]), self.species_list)]
136+
CHX = self.species_list[get_i_thing(Species(molecule=[Molecule().from_smiles("[CH]#*")]), self.species_list)]
137+
CH2X = self.species_list[get_i_thing(Species(molecule=[Molecule().from_smiles("[CH2]=*")]), self.species_list)]
138+
139+
self.reaction_list = [
140+
Reaction(
141+
reactants=[X, X, O2],
142+
products=[OX, OX],
143+
kinetics=StickingCoefficient(A=0.07, n=0, Ea=(0,'kcal/mol'), T0=(1,'K'))
144+
),
145+
Reaction(
146+
reactants=[X, X, H2],
147+
products=[HX, HX],
148+
kinetics=StickingCoefficient(A=0.046, n=0, Ea=(0,'kcal/mol'), T0=(1,'K'))
149+
),
150+
Reaction(
151+
reactants=[HX, CH3X],
152+
products=[X, X, CH4],
153+
kinetics=SurfaceArrhenius(A=(3.3e+21,'cm^2/(mol*s)'), n=0, Ea=(11.95,'kcal/mol'), T0=(1,'K'))
154+
),
155+
Reaction(
156+
reactants=[X, OHX],
157+
products=[OX, HX],
158+
kinetics=SurfaceArrhenius(A=(7.39e+19,'cm^2/(mol*s)'), n=0, Ea=(18.475,'kcal/mol'), T0=(1,'K'))
159+
),
160+
Reaction(
161+
reactants=[X, H2O],
162+
products=[H2OX],
163+
kinetics=StickingCoefficient(A=0.75, n=0, Ea=(0,'kcal/mol'), T0=(1,'K'))
164+
),
165+
Reaction(
166+
reactants=[X, H2OX],
167+
products=[HX, OHX],
168+
kinetics=SurfaceArrhenius(A=(1.15e+19,'cm^2/(mol*s)'), n=0, Ea=(24.235,'kcal/mol'), T0=(1,'K'))
169+
),
170+
Reaction(
171+
reactants=[OHX, CHX],
172+
products=[OX, CH2X],
173+
kinetics=SurfaceArrhenius(A=(4.4e+22,'cm^2/(mol*s)'), n=0.101, Ea=(10.143,'kcal/mol'), T0=(1,'K'))
174+
),
175+
Reaction(
176+
reactants=[OHX, CH2X],
177+
products=[OX, CH3X],
178+
kinetics=SurfaceArrhenius(A=(1.39e+21,'cm^2/(mol*s)'), n=0.101, Ea=(4.541,'kcal/mol'), T0=(1,'K'))
179+
),
180+
Reaction(
181+
reactants=[X, OX, CH4],
182+
products=[OHX, CH3X],
183+
kinetics=SurfaceArrhenius(A=(5e+18,'cm^4/(mol^2*s)'), n=0.7, Ea=(10.038,'kcal/mol'), T0=(1,'K'))
184+
)
185+
]
186+
45187
def test_solve_h2(self):
46188
"""
47189
Test the surface batch reactor with a dissociative adsorption of H2
@@ -1175,3 +1317,82 @@ def test_solve_ch3_thermo_coverage_dependence(self):
11751317
# Check that coverages are different
11761318
assert not np.allclose(y, y_off)
11771319
assert not np.allclose(species_rates, species_rates_off)
1320+
1321+
def test_jacobian_with_finite_differences(self):
1322+
1323+
# make a simple mechanism
1324+
x_O2 = 3e-5
1325+
x_CH4 = 1e-5
1326+
T = (900, 'K')
1327+
P = (1, 'atm')
1328+
1329+
surface_volume_ratio = (1.0, "m^-1") # TODO, try higher surface volume ratios
1330+
surface_site_density = (2.483e-8, "kmol/m^2") # read from Cantera yaml just to be sure these match
1331+
termination = TerminationTime((1.0, 's'))
1332+
1333+
1334+
Ar = self.species_list[get_i_thing(Species(smiles='[Ar]'), self.species_list)]
1335+
CH4 = self.species_list[get_i_thing(Species(smiles='C'), self.species_list)]
1336+
CO2 = self.species_list[get_i_thing(Species(smiles='O=C=O'), self.species_list)]
1337+
O2 = self.species_list[get_i_thing(Species(smiles='[O][O]'), self.species_list)]
1338+
X = self.species_list[get_i_thing(Species(smiles='*'), self.species_list)]
1339+
1340+
initial_gas_mole_fractions = {O2: x_O2, CH4: x_CH4, Ar: 1.0 - x_O2 - x_CH4}
1341+
initial_surface_coverages = {X: 1.0}
1342+
sensitive_species = [CO2]
1343+
1344+
reaction_system = SurfaceReactor(
1345+
T,
1346+
P,
1347+
n_sims=1,
1348+
initial_gas_mole_fractions=initial_gas_mole_fractions,
1349+
initial_surface_coverages=initial_surface_coverages,
1350+
surface_volume_ratio=surface_volume_ratio,
1351+
surface_site_density=surface_site_density,
1352+
termination=[termination],
1353+
sensitive_species=sensitive_species,
1354+
)
1355+
1356+
# Need to list some sens_worksheets so this doesn't crash when sensitivity is turned on
1357+
reaction_system.simulate(
1358+
core_species=self.species_list,
1359+
core_reactions=self.reaction_list,
1360+
edge_species=[],
1361+
edge_reactions=[],
1362+
surface_species=[],
1363+
surface_reactions=[],
1364+
model_settings=ModelSettings(tol_move_to_core=1e5), # tol_move_to_core isn't set by default which causes an error
1365+
simulator_settings=SimulatorSettings(), # defaults
1366+
sensitivity=True,
1367+
sens_worksheet=['temp_sensitivity.csv'],
1368+
)
1369+
1370+
# now comapare the jacobian to a finite difference approximation of the jacobian for the sensitive species
1371+
t = reaction_system.t
1372+
y = reaction_system.y
1373+
dydt = reaction_system.dydt
1374+
1375+
J_analytical = reaction_system.jacobian(t, y, dydt, cj=0.0)
1376+
1377+
# turn off sensitivity so that it only handles the first n values of y and doesn't get bogged
1378+
# down in all the partial derivatives that get stored in y for sensitivity calculations
1379+
reaction_system.sensitivity = False
1380+
1381+
n = reaction_system.num_core_species
1382+
eps_machine = np.sqrt(np.finfo(float).eps)
1383+
n_scale = 1e-10
1384+
1385+
# Get baseline residual f(n)
1386+
f0, _ = reaction_system.residual(t, y[:n], dydt[:n])
1387+
1388+
# Build finite difference Jacobian by column
1389+
J_fd = np.zeros((n, n))
1390+
for s in range(n):
1391+
y_perturbed = y.copy()
1392+
delta_n_s = eps_machine * max(abs(y[s]), n_scale) # figure out how much to perturb the value
1393+
y_perturbed[s] += delta_n_s # add \Delta n_s
1394+
1395+
f_perturbed, _ = reaction_system.residual(t, y_perturbed[:n], dydt[:n])
1396+
J_fd[:, s] = (f_perturbed - f0) / delta_n_s
1397+
1398+
assert np.allclose(J_analytical, J_fd, rtol=1e-2, atol=1e-6)

0 commit comments

Comments
 (0)