Skip to content

Commit 9c8e827

Browse files
12Chao12Chao
committed
Add unit tests for thermo coverage dependent models in surface solver tests
1 parent af53b93 commit 9c8e827

File tree

1 file changed

+394
-0
lines changed

1 file changed

+394
-0
lines changed

‎test/rmgpy/solver/solverSurfaceTest.py‎

Lines changed: 394 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -762,3 +762,397 @@ def test_solve_ch3_coverage_dependence(self):
762762
# Check that coverages are different
763763
assert not np.allclose(y, y_off)
764764
assert not np.allclose(species_rates, species_rates_off)
765+
766+
def test_solve_h2_thermo_coverage_dependence(self):
767+
"""
768+
Test the surface batch reactor can properly apply thermo coverage dependent parameters
769+
with the dissociative adsorption of H2.
770+
771+
Here we choose a kinetic model consisting of the dissociative adsorption reaction
772+
H2 + 2X <=> 2 HX
773+
We use a SurfaceArrhenius for the rate expression.
774+
"""
775+
h2 = Species(
776+
molecule=[Molecule().from_smiles("[H][H]")],
777+
thermo=ThermoData(
778+
Tdata=([300, 400, 500, 600, 800, 1000, 1500], "K"),
779+
Cpdata=(
780+
[6.955, 6.955, 6.956, 6.961, 7.003, 7.103, 7.502],
781+
"cal/(mol*K)",
782+
),
783+
H298=(0, "kcal/mol"),
784+
S298=(31.129, "cal/(mol*K)"),
785+
),
786+
)
787+
788+
x = Species(
789+
molecule=[Molecule().from_adjacency_list("1 X u0 p0")],
790+
thermo=ThermoData(
791+
Tdata=([300, 400, 500, 600, 800, 1000, 1500], "K"),
792+
Cpdata=([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], "cal/(mol*K)"),
793+
H298=(0.0, "kcal/mol"),
794+
S298=(0.0, "cal/(mol*K)"),
795+
),
796+
)
797+
798+
hx = Species(
799+
molecule=[Molecule().from_adjacency_list("1 H u0 p0 {2,S} \n 2 X u0 p0 {1,S}")],
800+
thermo=ThermoData(
801+
Tdata=([300, 400, 500, 600, 800, 1000, 1500], "K"),
802+
Cpdata=([1.50, 2.58, 3.40, 4.00, 4.73, 5.13, 5.57], "cal/(mol*K)"),
803+
H298=(-11.26, "kcal/mol"),
804+
S298=(0.44, "cal/(mol*K)"),
805+
),
806+
)
807+
hx.thermo.thermo_coverage_dependence = {hx:{'model':'polynomial', 'enthalpy-coefficients':[1,2,3], "entropy-coefficients":[1,5,3]},}
808+
809+
rxn1 = Reaction(
810+
reactants=[h2, x, x],
811+
products=[hx, hx],
812+
kinetics=SurfaceArrhenius(
813+
A=(9.05e18, "cm^5/(mol^2*s)"),
814+
n=0.5,
815+
Ea=(5.0, "kJ/mol"),
816+
T0=(1.0, "K"),
817+
),
818+
)
819+
820+
rxn2 = Reaction(
821+
reactants=[h2, x, x],
822+
products=[hx, hx],
823+
kinetics=SurfaceArrhenius(
824+
A=(9.05e-18, "cm^5/(mol^2*s)"), # 1e36 times slower
825+
n=0.5,
826+
Ea=(5.0, "kJ/mol"),
827+
T0=(1.0, "K"),
828+
),
829+
)
830+
831+
core_species = [h2, x, hx]
832+
edge_species = []
833+
core_reactions = [rxn1]
834+
edge_reactions = []
835+
836+
T = 600
837+
P_initial = 1.0e5
838+
rxn_system = SurfaceReactor(
839+
T,
840+
P_initial,
841+
n_sims=1,
842+
initial_gas_mole_fractions={h2: 1.0},
843+
initial_surface_coverages={x: 1.0},
844+
surface_volume_ratio=(1e1, "m^-1"),
845+
surface_site_density=(2.72e-9, "mol/cm^2"),
846+
coverage_dependence=True,
847+
termination=[],
848+
)
849+
850+
rxn_system.initialize_model(core_species, core_reactions, edge_species, edge_reactions)
851+
852+
tlist = np.logspace(-13, -5, 81, dtype=float)
853+
854+
assert isinstance(hx.thermo.thermo_coverage_dependence, dict) # check to make sure coverage_dependence is still the correct type
855+
for species, parameters in hx.thermo.thermo_coverage_dependence.items():
856+
assert isinstance(species, Species) # species should be a Species
857+
assert isinstance(parameters, dict)
858+
assert parameters["model"] is not None
859+
assert parameters["enthalpy-coefficients"] is not None
860+
assert parameters["entropy-coefficients"] is not None
861+
# Integrate to get the solution at each time point
862+
t = []
863+
y = []
864+
reaction_rates = []
865+
species_rates = []
866+
start_time = time.time()
867+
for t1 in tlist:
868+
rxn_system.advance(t1)
869+
t.append(rxn_system.t)
870+
# You must make a copy of y because it is overwritten by DASSL at
871+
# each call to advance()
872+
y.append(rxn_system.y.copy())
873+
reaction_rates.append(rxn_system.core_reaction_rates.copy())
874+
species_rates.append(rxn_system.core_species_rates.copy())
875+
run_time = time.time() - start_time
876+
print(f"Simulation took {run_time:.3e} seconds")
877+
878+
# Convert the solution vectors to np arrays
879+
t = np.array(t, float)
880+
y = np.array(y, float)
881+
reaction_rates = np.array(reaction_rates, float)
882+
species_rates = np.array(species_rates, float)
883+
total_sites = y[0, 1]
884+
885+
# Check that we're computing the species fluxes correctly
886+
for i in range(t.shape[0]):
887+
assert abs(reaction_rates[i, 0] - -1.0 * species_rates[i, 0]) < abs(
888+
1e-6 * reaction_rates[i, 0]
889+
)
890+
assert abs(reaction_rates[i, 0] - -0.5 * species_rates[i, 1]) < abs(
891+
1e-6 * reaction_rates[i, 0]
892+
)
893+
assert abs(reaction_rates[i, 0] - 0.5 * species_rates[i, 2]) < abs(
894+
1e-6 * reaction_rates[i, 0]
895+
)
896+
897+
# Check that we've reached equilibrium
898+
assert abs(reaction_rates[-1, 0] - 0.0) < 1e-2
899+
900+
def test_solve_ch3_thermo_coverage_dependence(self):
901+
"""
902+
Test the surface batch reactor can properly apply coverage dependent parameters
903+
with the nondissociative adsorption of CH3
904+
905+
Here we choose a kinetic model consisting of the adsorption reaction
906+
CH3 + X <=> CH3X
907+
We use a sticking coefficient for the rate expression.
908+
"""
909+
910+
ch3 = Species(
911+
molecule=[Molecule().from_smiles("[CH3]")],
912+
thermo=NASA(
913+
polynomials=[
914+
NASAPolynomial(
915+
coeffs=[
916+
3.91547,
917+
0.00184155,
918+
3.48741e-06,
919+
-3.32746e-09,
920+
8.49953e-13,
921+
16285.6,
922+
0.351743,
923+
],
924+
Tmin=(100, "K"),
925+
Tmax=(1337.63, "K"),
926+
),
927+
NASAPolynomial(
928+
coeffs=[
929+
3.54146,
930+
0.00476786,
931+
-1.82148e-06,
932+
3.28876e-10,
933+
-2.22545e-14,
934+
16224,
935+
1.66032,
936+
],
937+
Tmin=(1337.63, "K"),
938+
Tmax=(5000, "K"),
939+
),
940+
],
941+
Tmin=(100, "K"),
942+
Tmax=(5000, "K"),
943+
E0=(135.382, "kJ/mol"),
944+
comment="""Thermo library: primaryThermoLibrary + radical(CH3)""",
945+
),
946+
molecular_weight=(15.0345, "amu"),
947+
)
948+
949+
x = Species(
950+
molecule=[Molecule().from_adjacency_list("1 X u0 p0")],
951+
thermo=NASA(
952+
polynomials=[
953+
NASAPolynomial(coeffs=[0, 0, 0, 0, 0, 0, 0], Tmin=(298, "K"), Tmax=(1000, "K")),
954+
NASAPolynomial(coeffs=[0, 0, 0, 0, 0, 0, 0], Tmin=(1000, "K"), Tmax=(2000, "K")),
955+
],
956+
Tmin=(298, "K"),
957+
Tmax=(2000, "K"),
958+
E0=(-6.19426, "kJ/mol"),
959+
comment="""Thermo library: surfaceThermo""",
960+
),
961+
)
962+
963+
ch3x = Species(
964+
molecule=[
965+
Molecule().from_adjacency_list(
966+
"""1 C u0 p0 c0 {2,S} {3,S} {4,S} {5,S}
967+
2 H u0 p0 c0 {1,S}
968+
3 H u0 p0 c0 {1,S}
969+
4 H u0 p0 c0 {1,S}
970+
5 X u0 p0 c0 {1,S}"""
971+
)
972+
],
973+
thermo=NASA(
974+
polynomials=[
975+
NASAPolynomial(
976+
coeffs=[
977+
-0.552219,
978+
0.026442,
979+
-3.55617e-05,
980+
2.60044e-08,
981+
-7.52707e-12,
982+
-4433.47,
983+
0.692144,
984+
],
985+
Tmin=(298, "K"),
986+
Tmax=(1000, "K"),
987+
),
988+
NASAPolynomial(
989+
coeffs=[
990+
3.62557,
991+
0.00739512,
992+
-2.43797e-06,
993+
1.86159e-10,
994+
3.6485e-14,
995+
-5187.22,
996+
-18.9668,
997+
],
998+
Tmin=(1000, "K"),
999+
Tmax=(2000, "K"),
1000+
),
1001+
],
1002+
Tmin=(298, "K"),
1003+
Tmax=(2000, "K"),
1004+
E0=(-39.1285, "kJ/mol"),
1005+
comment="""Thermo library: surfaceThermoNi111""",
1006+
),
1007+
)
1008+
ch3x.thermo.thermo_coverage_dependence = {ch3x:{'model':'polynomial', 'enthalpy-coefficients':[1,2,3], "entropy-coefficients":[1,5,3]},}
1009+
1010+
rxn1 = Reaction(
1011+
reactants=[ch3, x],
1012+
products=[ch3x],
1013+
kinetics=StickingCoefficient(
1014+
A=0.1,
1015+
n=0,
1016+
Ea=(0, "kcal/mol"),
1017+
T0=(1, "K"),
1018+
Tmin=(200, "K"),
1019+
Tmax=(3000, "K"),
1020+
coverage_dependence={x: {"a": 0.0, "m": -1.0, "E": (0.0, "J/mol")}},
1021+
comment="""Exact match found for rate rule (Adsorbate;VacantSite)""",
1022+
),
1023+
)
1024+
core_species = [ch3, x, ch3x]
1025+
edge_species = []
1026+
core_reactions = [rxn1]
1027+
edge_reactions = []
1028+
1029+
T = 800.0
1030+
P_initial = 1.0e5
1031+
rxn_system = SurfaceReactor(
1032+
T,
1033+
P_initial,
1034+
n_sims=1,
1035+
initial_gas_mole_fractions={ch3: 1.0},
1036+
initial_surface_coverages={x: 1.0},
1037+
surface_volume_ratio=(1.0, "m^-1"),
1038+
surface_site_density=(2.72e-9, "mol/cm^2"),
1039+
coverage_dependence=True,
1040+
termination=[],
1041+
)
1042+
# in chemkin, the sites are mostly occupied in about 1e-8 seconds.
1043+
1044+
rxn_system.initialize_model(core_species, core_reactions, edge_species, edge_reactions)
1045+
1046+
tlist = np.logspace(-13, -5, 81, dtype=float)
1047+
1048+
print("Surface site density:", rxn_system.surface_site_density.value_si)
1049+
1050+
print(
1051+
"rxn1 rate coefficient",
1052+
rxn1.get_surface_rate_coefficient(rxn_system.T.value_si, rxn_system.surface_site_density.value_si),
1053+
)
1054+
1055+
assert isinstance(ch3.thermo.thermo_coverage_dependence, dict) # check to make sure coverage_dependence is still the correct type
1056+
for species, parameters in ch3.thermo.thermo_coverage_dependence.items():
1057+
assert isinstance(species, Species) # species should be a Species
1058+
assert isinstance(parameters, dict)
1059+
assert parameters["model"] is not None
1060+
assert parameters["enthalpy-coefficients"] is not None
1061+
assert parameters["entropy-coefficients"] is not None
1062+
1063+
# Integrate to get the solution at each time point
1064+
t = []
1065+
y = []
1066+
reaction_rates = []
1067+
species_rates = []
1068+
t.append(rxn_system.t)
1069+
# You must make a copy of y because it is overwritten by DASSL at
1070+
# each call to advance()
1071+
y.append(rxn_system.y.copy())
1072+
reaction_rates.append(rxn_system.core_reaction_rates.copy())
1073+
species_rates.append(rxn_system.core_species_rates.copy())
1074+
print("time: ", t)
1075+
print("moles:", y)
1076+
print("reaction rates:", reaction_rates)
1077+
print("species rates:", species_rates)
1078+
start_time = time.time()
1079+
for t1 in tlist:
1080+
rxn_system.advance(t1)
1081+
t.append(rxn_system.t)
1082+
# You must make a copy of y because it is overwritten by DASSL at
1083+
# each call to advance()
1084+
y.append(rxn_system.y.copy())
1085+
reaction_rates.append(rxn_system.core_reaction_rates.copy())
1086+
species_rates.append(rxn_system.core_species_rates.copy())
1087+
run_time = time.time() - start_time
1088+
print(f"Simulation took {run_time:.3e} seconds")
1089+
1090+
# Convert the solution vectors to np arrays
1091+
t = np.array(t, float)
1092+
y = np.array(y, float)
1093+
reaction_rates = np.array(reaction_rates, float)
1094+
species_rates = np.array(species_rates, float)
1095+
V = constants.R * rxn_system.T.value_si * np.sum(y) / rxn_system.P_initial.value_si
1096+
1097+
# Check that we're computing the species fluxes correctly
1098+
for i in range(t.shape[0]):
1099+
assert abs(reaction_rates[i, 0] - -species_rates[i, 0]) < abs(
1100+
1e-6 * reaction_rates[i, 0]
1101+
)
1102+
assert abs(reaction_rates[i, 0] - -species_rates[i, 1]) < abs(
1103+
1e-6 * reaction_rates[i, 0]
1104+
)
1105+
assert abs(reaction_rates[i, 0] - species_rates[i, 2]) < abs(
1106+
1e-6 * reaction_rates[i, 0]
1107+
)
1108+
1109+
# Check that we've reached equilibrium by the end
1110+
assert abs(reaction_rates[-1, 0] - 0.0) < 1e-2
1111+
1112+
# Run model with Covdep off so we can test that it is actually being implemented
1113+
rxn_system = SurfaceReactor(
1114+
T,
1115+
P_initial,
1116+
n_sims=1,
1117+
initial_gas_mole_fractions={ch3: 1.0},
1118+
initial_surface_coverages={x: 1.0},
1119+
surface_volume_ratio=(1.0, "m^-1"),
1120+
surface_site_density=(2.72e-9, "mol/cm^2"),
1121+
termination=[],
1122+
)
1123+
1124+
rxn_system.initialize_model(core_species, core_reactions, edge_species, edge_reactions)
1125+
1126+
tlist = np.logspace(-13, -5, 81, dtype=float)
1127+
1128+
# Integrate to get the solution at each time point
1129+
t = []
1130+
y_off = []
1131+
species_rates_off = []
1132+
t.append(rxn_system.t)
1133+
1134+
# You must make a copy of y because it is overwritten by DASSL at
1135+
# each call to advance()
1136+
y_off.append(rxn_system.y.copy())
1137+
species_rates_off.append(rxn_system.core_species_rates.copy())
1138+
for t1 in tlist:
1139+
rxn_system.advance(t1)
1140+
t.append(rxn_system.t)
1141+
# You must make a copy of y because it is overwritten by DASSL at
1142+
# each call to advance()
1143+
y_off.append(rxn_system.y.copy())
1144+
species_rates_off.append(rxn_system.core_species_rates.copy())
1145+
run_time = time.time() - start_time
1146+
print(f"Simulation took {run_time:.3e} seconds")
1147+
1148+
# Convert the solution vectors to np arrays
1149+
t = np.array(t, float)
1150+
y_off = np.array(y_off, float)
1151+
species_rates_off = np.array(species_rates_off, float)
1152+
1153+
# Check that we've reached equilibrium
1154+
assert abs(species_rates_off[-1, 0] - 0.0) < 1e-2
1155+
1156+
# Check that coverages are different
1157+
assert not np.allclose(y, y_off)
1158+
assert not np.allclose(species_rates, species_rates_off)

0 commit comments

Comments
 (0)