Skip to content

Commit d1f9807

Browse files
sevyharrissevyharris
committed
Cache inv omega for faster jacobian execution
the jacobian function was populating an array with 1/V or 1/A for each gas or surface species each time it called the jacobian. Now it caches the values. 1/V and 1/A do not change because it's a constant volume reactor, but the total number of core species could change as RMG grows the model, and this enables the use of cached values that get updated only as required
1 parent bddb491 commit d1f9807

1 file changed

Lines changed: 50 additions & 35 deletions

File tree

‎rmgpy/solver/surface.pyx‎

Lines changed: 50 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -70,6 +70,8 @@ cdef class SurfaceReactor(ReactionSystem):
7070
cdef public np.ndarray species_on_surface # (catalyst surface, not core/edge surface)
7171
cdef public np.ndarray thermo_coeff_matrix
7272
cdef public np.ndarray stoi_matrix
73+
cdef public np.ndarray _jacobian_inv_omega
74+
cdef public bint _jacobian_inv_omega_num_core_species
7375

7476
cdef public bint coverage_dependence
7577
cdef public dict coverage_dependencies
@@ -115,6 +117,8 @@ cdef class SurfaceReactor(ReactionSystem):
115117
self.constant_volume = True
116118
self.sens_conditions = sens_conditions
117119
self.n_sims = n_sims
120+
self._jacobian_inv_omega = None # array of 1/V (gas-phase) or 1/A (surface-phase) for corresponding gas or surface species
121+
self._jacobian_inv_omega_num_core_species = 0
118122

119123
self.coverage_dependencies = {}
120124

@@ -380,11 +384,15 @@ cdef class SurfaceReactor(ReactionSystem):
380384
i = self.get_species_index(spec)
381385
self.y0[i] = total_surface_sites * coverage # moles in reactor
382386

387+
# array of stored 1/V (gas-phase) or 1/A (surface-phase) for corresponding gas or surface species to help compute Jacobian
388+
self._jacobian_inv_omega = np.full(self.num_core_species, 1.0 / V, dtype=np.float64)
383389
for j, isSurfaceSpecies in enumerate(self.species_on_surface): # should only go up to core species
384390
if isSurfaceSpecies:
385391
self.core_species_concentrations[j] = self.y0[j] / V / surface_volume_ratio_si # moles per m2 of surface
392+
self._jacobian_inv_omega[j] = 1.0 / (V * surface_volume_ratio_si)
386393
else:
387394
self.core_species_concentrations[j] = self.y0[j] / V # moles per m3 of gas
395+
self._jacobian_inv_omega_num_core_species = self.num_core_species
388396

389397
def compute_network_variables(self, pdep_networks=None):
390398
# ToDo: this should allow pressure to vary?
@@ -688,16 +696,23 @@ cdef class SurfaceReactor(ReactionSystem):
688696
V = self.V
689697
A = V * self.surface_volume_ratio.value_si # surface area in m^2
690698

691-
# Per-species normalization: Omega[i] = A for surface species, V for gas
692-
Omega = np.full(num_core_species, V, dtype=np.float64)
693-
for i in range(num_core_species):
694-
if self.species_on_surface[i]:
695-
Omega[i] = A
699+
# Cache inverse per-species normalization on the instance since
700+
# species_on_surface is static for SurfaceReactor and V/A are constant
701+
inv_omega = self._jacobian_inv_omega
702+
cached_num_core_species = self._jacobian_inv_omega_num_core_species
703+
704+
if cached_num_core_species != num_core_species:
705+
inv_omega = np.full(num_core_species, 1.0 / V, dtype=np.float64)
706+
for i in range(num_core_species):
707+
if self.species_on_surface[i]:
708+
inv_omega[i] = 1.0 / A
709+
self._jacobian_inv_omega = inv_omega
710+
self._jacobian_inv_omega_num_core_species = num_core_species
696711

697712
# Concentrations — surface species in mol/m^2, gas in mol/m^3
698713
C = np.zeros_like(self.core_species_concentrations)
699714
for j in range(num_core_species):
700-
C[j] = y[j] / Omega[j]
715+
C[j] = y[j] * inv_omega[j]
701716

702717
for j in range(num_core_reactions):
703718

@@ -707,8 +722,8 @@ cdef class SurfaceReactor(ReactionSystem):
707722
k = kf[j]
708723

709724
if ir[j, 1] == -1: # unimolecular forward
710-
# deriv = k * V / Omega[s]; gas: k*1 = k; surface: k*V/A = k/svr
711-
deriv = k * V / Omega[ir[j, 0]]
725+
# deriv = k * V * inv_omega[s]; gas: k*1 = k; surface: k*V/A = k/svr
726+
deriv = k * V * inv_omega[ir[j, 0]]
712727
pd[ir[j, 0], ir[j, 0]] -= deriv
713728

714729
pd[ip[j, 0], ir[j, 0]] += deriv
@@ -719,7 +734,7 @@ cdef class SurfaceReactor(ReactionSystem):
719734

720735
elif ir[j, 2] == -1: # bimolecular forward
721736
if ir[j, 0] == ir[j, 1]: # A + A -> products
722-
deriv = 2 * k * C[ir[j, 0]] * V / Omega[ir[j, 0]]
737+
deriv = 2 * k * C[ir[j, 0]] * V * inv_omega[ir[j, 0]]
723738
pd[ir[j, 0], ir[j, 0]] -= 2 * deriv
724739

725740
pd[ip[j, 0], ir[j, 0]] += deriv
@@ -730,7 +745,7 @@ cdef class SurfaceReactor(ReactionSystem):
730745

731746
else: # A + B -> products
732747
# derivative with respect to n[A]
733-
deriv = k * C[ir[j, 1]] * V / Omega[ir[j, 0]]
748+
deriv = k * C[ir[j, 1]] * V * inv_omega[ir[j, 0]]
734749
pd[ir[j, 0], ir[j, 0]] -= deriv
735750
pd[ir[j, 1], ir[j, 0]] -= deriv
736751

@@ -741,7 +756,7 @@ cdef class SurfaceReactor(ReactionSystem):
741756
pd[ip[j, 2], ir[j, 0]] += deriv
742757

743758
# derivative with respect to n[B]
744-
deriv = k * C[ir[j, 0]] * V / Omega[ir[j, 1]]
759+
deriv = k * C[ir[j, 0]] * V * inv_omega[ir[j, 1]]
745760
pd[ir[j, 0], ir[j, 1]] -= deriv
746761
pd[ir[j, 1], ir[j, 1]] -= deriv
747762

@@ -754,7 +769,7 @@ cdef class SurfaceReactor(ReactionSystem):
754769
else: # trimolecular forward
755770

756771
if (ir[j, 0] == ir[j, 1]) and (ir[j, 0] == ir[j, 2]): # A+A+A
757-
deriv = 3 * k * C[ir[j, 0]] * C[ir[j, 0]] * V / Omega[ir[j, 0]]
772+
deriv = 3 * k * C[ir[j, 0]] * C[ir[j, 0]] * V * inv_omega[ir[j, 0]]
758773
pd[ir[j, 0], ir[j, 0]] -= 3 * deriv
759774

760775
pd[ip[j, 0], ir[j, 0]] += deriv
@@ -765,7 +780,7 @@ cdef class SurfaceReactor(ReactionSystem):
765780

766781
elif ir[j, 0] == ir[j, 1]: # A+A+B
767782
# derivative with respect to n[A]
768-
deriv = 2 * k * C[ir[j, 0]] * C[ir[j, 2]] * V / Omega[ir[j, 0]]
783+
deriv = 2 * k * C[ir[j, 0]] * C[ir[j, 2]] * V * inv_omega[ir[j, 0]]
769784
pd[ir[j, 0], ir[j, 0]] -= 2 * deriv
770785
pd[ir[j, 2], ir[j, 0]] -= deriv
771786

@@ -776,7 +791,7 @@ cdef class SurfaceReactor(ReactionSystem):
776791
pd[ip[j, 2], ir[j, 0]] += deriv
777792

778793
# derivative with respect to n[B]
779-
deriv = k * C[ir[j, 0]] * C[ir[j, 0]] * V / Omega[ir[j, 2]]
794+
deriv = k * C[ir[j, 0]] * C[ir[j, 0]] * V * inv_omega[ir[j, 2]]
780795
pd[ir[j, 0], ir[j, 2]] -= 2 * deriv
781796
pd[ir[j, 2], ir[j, 2]] -= deriv
782797

@@ -788,7 +803,7 @@ cdef class SurfaceReactor(ReactionSystem):
788803

789804
elif ir[j, 1] == ir[j, 2]: # A+B+B
790805
# derivative with respect to n[A]
791-
deriv = k * C[ir[j, 1]] * C[ir[j, 1]] * V / Omega[ir[j, 0]]
806+
deriv = k * C[ir[j, 1]] * C[ir[j, 1]] * V * inv_omega[ir[j, 0]]
792807
pd[ir[j, 0], ir[j, 0]] -= deriv
793808
pd[ir[j, 1], ir[j, 0]] -= 2 * deriv
794809

@@ -799,7 +814,7 @@ cdef class SurfaceReactor(ReactionSystem):
799814
pd[ip[j, 2], ir[j, 0]] += deriv
800815

801816
# derivative with respect to n[B]
802-
deriv = 2 * k * C[ir[j, 0]] * C[ir[j, 1]] * V / Omega[ir[j, 1]]
817+
deriv = 2 * k * C[ir[j, 0]] * C[ir[j, 1]] * V * inv_omega[ir[j, 1]]
803818
pd[ir[j, 0], ir[j, 1]] -= deriv
804819
pd[ir[j, 1], ir[j, 1]] -= 2 * deriv
805820

@@ -811,7 +826,7 @@ cdef class SurfaceReactor(ReactionSystem):
811826

812827
elif ir[j, 0] == ir[j, 2]: # A+B+A
813828
# derivative with respect to n[A]
814-
deriv = 2 * k * C[ir[j, 0]] * C[ir[j, 1]] * V / Omega[ir[j, 0]]
829+
deriv = 2 * k * C[ir[j, 0]] * C[ir[j, 1]] * V * inv_omega[ir[j, 0]]
815830
pd[ir[j, 0], ir[j, 0]] -= 2 * deriv
816831
pd[ir[j, 1], ir[j, 0]] -= deriv
817832

@@ -822,7 +837,7 @@ cdef class SurfaceReactor(ReactionSystem):
822837
pd[ip[j, 2], ir[j, 0]] += deriv
823838

824839
# derivative with respect to n[B]
825-
deriv = k * C[ir[j, 0]] * C[ir[j, 0]] * V / Omega[ir[j, 1]]
840+
deriv = k * C[ir[j, 0]] * C[ir[j, 0]] * V * inv_omega[ir[j, 1]]
826841
pd[ir[j, 0], ir[j, 1]] -= 2 * deriv
827842
pd[ir[j, 1], ir[j, 1]] -= deriv
828843

@@ -834,7 +849,7 @@ cdef class SurfaceReactor(ReactionSystem):
834849

835850
else: # A+B+C, all distinct
836851
# derivative with respect to n[A]
837-
deriv = k * C[ir[j, 1]] * C[ir[j, 2]] * V / Omega[ir[j, 0]]
852+
deriv = k * C[ir[j, 1]] * C[ir[j, 2]] * V * inv_omega[ir[j, 0]]
838853
pd[ir[j, 0], ir[j, 0]] -= deriv
839854
pd[ir[j, 1], ir[j, 0]] -= deriv
840855
pd[ir[j, 2], ir[j, 0]] -= deriv
@@ -846,7 +861,7 @@ cdef class SurfaceReactor(ReactionSystem):
846861
pd[ip[j, 2], ir[j, 0]] += deriv
847862

848863
# derivative with respect to n[B]
849-
deriv = k * C[ir[j, 0]] * C[ir[j, 2]] * V / Omega[ir[j, 1]]
864+
deriv = k * C[ir[j, 0]] * C[ir[j, 2]] * V * inv_omega[ir[j, 1]]
850865
pd[ir[j, 0], ir[j, 1]] -= deriv
851866
pd[ir[j, 1], ir[j, 1]] -= deriv
852867
pd[ir[j, 2], ir[j, 1]] -= deriv
@@ -858,7 +873,7 @@ cdef class SurfaceReactor(ReactionSystem):
858873
pd[ip[j, 2], ir[j, 1]] += deriv
859874

860875
# derivative with respect to n[C]
861-
deriv = k * C[ir[j, 0]] * C[ir[j, 1]] * V / Omega[ir[j, 2]]
876+
deriv = k * C[ir[j, 0]] * C[ir[j, 1]] * V * inv_omega[ir[j, 2]]
862877
pd[ir[j, 0], ir[j, 2]] -= deriv
863878
pd[ir[j, 1], ir[j, 2]] -= deriv
864879
pd[ir[j, 2], ir[j, 2]] -= deriv
@@ -875,7 +890,7 @@ cdef class SurfaceReactor(ReactionSystem):
875890
k = kr[j]
876891

877892
if ip[j, 1] == -1: # unimolecular reverse
878-
deriv = k * V / Omega[ip[j, 0]]
893+
deriv = k * V * inv_omega[ip[j, 0]]
879894
pd[ip[j, 0], ip[j, 0]] -= deriv
880895

881896
pd[ir[j, 0], ip[j, 0]] += deriv
@@ -886,7 +901,7 @@ cdef class SurfaceReactor(ReactionSystem):
886901

887902
elif ip[j, 2] == -1: # bimolecular reverse
888903
if ip[j, 0] == ip[j, 1]: # P+P -> reactants
889-
deriv = 2 * k * C[ip[j, 0]] * V / Omega[ip[j, 0]]
904+
deriv = 2 * k * C[ip[j, 0]] * V * inv_omega[ip[j, 0]]
890905
pd[ip[j, 0], ip[j, 0]] -= 2 * deriv
891906

892907
pd[ir[j, 0], ip[j, 0]] += deriv
@@ -897,7 +912,7 @@ cdef class SurfaceReactor(ReactionSystem):
897912

898913
else: # P1+P2 -> reactants
899914
# derivative with respect to n[P1]
900-
deriv = k * C[ip[j, 1]] * V / Omega[ip[j, 0]]
915+
deriv = k * C[ip[j, 1]] * V * inv_omega[ip[j, 0]]
901916
pd[ip[j, 0], ip[j, 0]] -= deriv
902917
pd[ip[j, 1], ip[j, 0]] -= deriv
903918

@@ -908,7 +923,7 @@ cdef class SurfaceReactor(ReactionSystem):
908923
pd[ir[j, 2], ip[j, 0]] += deriv
909924

910925
# derivative with respect to n[P2]
911-
deriv = k * C[ip[j, 0]] * V / Omega[ip[j, 1]]
926+
deriv = k * C[ip[j, 0]] * V * inv_omega[ip[j, 1]]
912927
pd[ip[j, 0], ip[j, 1]] -= deriv
913928
pd[ip[j, 1], ip[j, 1]] -= deriv
914929

@@ -921,7 +936,7 @@ cdef class SurfaceReactor(ReactionSystem):
921936
else: # trimolecular reverse
922937

923938
if (ip[j, 0] == ip[j, 1]) and (ip[j, 0] == ip[j, 2]): # P+P+P
924-
deriv = 3 * k * C[ip[j, 0]] * C[ip[j, 0]] * V / Omega[ip[j, 0]]
939+
deriv = 3 * k * C[ip[j, 0]] * C[ip[j, 0]] * V * inv_omega[ip[j, 0]]
925940
pd[ip[j, 0], ip[j, 0]] -= 3 * deriv
926941

927942
pd[ir[j, 0], ip[j, 0]] += deriv
@@ -932,7 +947,7 @@ cdef class SurfaceReactor(ReactionSystem):
932947

933948
elif ip[j, 0] == ip[j, 1]: # P1+P1+P2
934949
# derivative with respect to n[P1]
935-
deriv = 2 * k * C[ip[j, 0]] * C[ip[j, 2]] * V / Omega[ip[j, 0]]
950+
deriv = 2 * k * C[ip[j, 0]] * C[ip[j, 2]] * V * inv_omega[ip[j, 0]]
936951
pd[ip[j, 0], ip[j, 0]] -= 2 * deriv
937952
pd[ip[j, 2], ip[j, 0]] -= deriv
938953

@@ -943,7 +958,7 @@ cdef class SurfaceReactor(ReactionSystem):
943958
pd[ir[j, 2], ip[j, 0]] += deriv
944959

945960
# derivative with respect to n[P2]
946-
deriv = k * C[ip[j, 0]] * C[ip[j, 0]] * V / Omega[ip[j, 2]]
961+
deriv = k * C[ip[j, 0]] * C[ip[j, 0]] * V * inv_omega[ip[j, 2]]
947962
pd[ip[j, 0], ip[j, 2]] -= 2 * deriv
948963
pd[ip[j, 2], ip[j, 2]] -= deriv
949964

@@ -955,7 +970,7 @@ cdef class SurfaceReactor(ReactionSystem):
955970

956971
elif ip[j, 1] == ip[j, 2]: # P1+P2+P2
957972
# derivative with respect to n[P1]
958-
deriv = k * C[ip[j, 1]] * C[ip[j, 1]] * V / Omega[ip[j, 0]]
973+
deriv = k * C[ip[j, 1]] * C[ip[j, 1]] * V * inv_omega[ip[j, 0]]
959974
pd[ip[j, 0], ip[j, 0]] -= deriv
960975
pd[ip[j, 1], ip[j, 0]] -= 2 * deriv
961976

@@ -966,7 +981,7 @@ cdef class SurfaceReactor(ReactionSystem):
966981
pd[ir[j, 2], ip[j, 0]] += deriv
967982

968983
# derivative with respect to n[P2]
969-
deriv = 2 * k * C[ip[j, 0]] * C[ip[j, 1]] * V / Omega[ip[j, 1]]
984+
deriv = 2 * k * C[ip[j, 0]] * C[ip[j, 1]] * V * inv_omega[ip[j, 1]]
970985
pd[ip[j, 0], ip[j, 1]] -= deriv
971986
pd[ip[j, 1], ip[j, 1]] -= 2 * deriv
972987

@@ -978,7 +993,7 @@ cdef class SurfaceReactor(ReactionSystem):
978993

979994
elif ip[j, 0] == ip[j, 2]: # P1+P2+P1
980995
# derivative with respect to n[P1]
981-
deriv = 2 * k * C[ip[j, 0]] * C[ip[j, 1]] * V / Omega[ip[j, 0]]
996+
deriv = 2 * k * C[ip[j, 0]] * C[ip[j, 1]] * V * inv_omega[ip[j, 0]]
982997
pd[ip[j, 0], ip[j, 0]] -= 2 * deriv
983998
pd[ip[j, 1], ip[j, 0]] -= deriv
984999

@@ -989,7 +1004,7 @@ cdef class SurfaceReactor(ReactionSystem):
9891004
pd[ir[j, 2], ip[j, 0]] += deriv
9901005

9911006
# derivative with respect to n[P2]
992-
deriv = k * C[ip[j, 0]] * C[ip[j, 0]] * V / Omega[ip[j, 1]]
1007+
deriv = k * C[ip[j, 0]] * C[ip[j, 0]] * V * inv_omega[ip[j, 1]]
9931008
pd[ip[j, 0], ip[j, 1]] -= 2 * deriv
9941009
pd[ip[j, 1], ip[j, 1]] -= deriv
9951010

@@ -1001,7 +1016,7 @@ cdef class SurfaceReactor(ReactionSystem):
10011016

10021017
else: # P1+P2+P3, all distinct
10031018
# derivative with respect to n[P1]
1004-
deriv = k * C[ip[j, 1]] * C[ip[j, 2]] * V / Omega[ip[j, 0]]
1019+
deriv = k * C[ip[j, 1]] * C[ip[j, 2]] * V * inv_omega[ip[j, 0]]
10051020
pd[ip[j, 0], ip[j, 0]] -= deriv
10061021
pd[ip[j, 1], ip[j, 0]] -= deriv
10071022
pd[ip[j, 2], ip[j, 0]] -= deriv
@@ -1013,7 +1028,7 @@ cdef class SurfaceReactor(ReactionSystem):
10131028
pd[ir[j, 2], ip[j, 0]] += deriv
10141029

10151030
# derivative with respect to n[P2]
1016-
deriv = k * C[ip[j, 0]] * C[ip[j, 2]] * V / Omega[ip[j, 1]]
1031+
deriv = k * C[ip[j, 0]] * C[ip[j, 2]] * V * inv_omega[ip[j, 1]]
10171032
pd[ip[j, 0], ip[j, 1]] -= deriv
10181033
pd[ip[j, 1], ip[j, 1]] -= deriv
10191034
pd[ip[j, 2], ip[j, 1]] -= deriv
@@ -1025,7 +1040,7 @@ cdef class SurfaceReactor(ReactionSystem):
10251040
pd[ir[j, 2], ip[j, 1]] += deriv
10261041

10271042
# derivative with respect to n[P3]
1028-
deriv = k * C[ip[j, 0]] * C[ip[j, 1]] * V / Omega[ip[j, 2]]
1043+
deriv = k * C[ip[j, 0]] * C[ip[j, 1]] * V * inv_omega[ip[j, 2]]
10291044
pd[ip[j, 0], ip[j, 2]] -= deriv
10301045
pd[ip[j, 1], ip[j, 2]] -= deriv
10311046
pd[ip[j, 2], ip[j, 2]] -= deriv

0 commit comments

Comments
 (0)