diff --git a/rmgpy/solver/base.pyx b/rmgpy/solver/base.pyx index 1aeebc63150..3470c95088d 100644 --- a/rmgpy/solver/base.pyx +++ b/rmgpy/solver/base.pyx @@ -714,7 +714,8 @@ cdef class ReactionSystem(DASx): if sensitivity: time_array = [] norm_sens_array = [[] for spec in self.sensitive_species] - RTP = constants.R * self.T.value_si / self.P.value_si + if not self.constant_volume: + RTP = constants.R * self.T.value_si / self.P.value_si # identify sensitive species indices sens_species_indices = np.array([species_index[spec] for spec in self.sensitive_species], int) # index within core_species list of the sensitive species diff --git a/rmgpy/solver/surface.pyx b/rmgpy/solver/surface.pyx index a9e0cd22b49..22272c4bf07 100644 --- a/rmgpy/solver/surface.pyx +++ b/rmgpy/solver/surface.pyx @@ -26,8 +26,8 @@ ############################################################################### """ -Contains the :class:`SimpleReactor` class, providing a reaction system -consisting of a homogeneous, isothermal, isobaric batch reactor. +Contains the :class:`SurfaceReactor` class, providing a reaction system +consisting of a heterogeneous, isothermal, constant volume batch reactor. """ import itertools @@ -70,6 +70,8 @@ cdef class SurfaceReactor(ReactionSystem): cdef public np.ndarray species_on_surface # (catalyst surface, not core/edge surface) cdef public np.ndarray thermo_coeff_matrix cdef public np.ndarray stoi_matrix + cdef public np.ndarray _jacobian_inv_omega + cdef public bint _jacobian_inv_omega_num_core_species cdef public bint coverage_dependence cdef public dict coverage_dependencies @@ -115,6 +117,8 @@ cdef class SurfaceReactor(ReactionSystem): self.constant_volume = True self.sens_conditions = sens_conditions self.n_sims = n_sims + self._jacobian_inv_omega = None # array of 1/V (gas-phase) or 1/A (surface-phase) for corresponding gas or surface species + self._jacobian_inv_omega_num_core_species = 0 self.coverage_dependencies = {} @@ -380,11 +384,15 @@ cdef class SurfaceReactor(ReactionSystem): i = self.get_species_index(spec) self.y0[i] = total_surface_sites * coverage # moles in reactor + # array of stored 1/V (gas-phase) or 1/A (surface-phase) for corresponding gas or surface species to help compute Jacobian + self._jacobian_inv_omega = np.full(self.num_core_species, 1.0 / V, dtype=np.float64) for j, isSurfaceSpecies in enumerate(self.species_on_surface): # should only go up to core species if isSurfaceSpecies: self.core_species_concentrations[j] = self.y0[j] / V / surface_volume_ratio_si # moles per m2 of surface + self._jacobian_inv_omega[j] = 1.0 / (V * surface_volume_ratio_si) else: self.core_species_concentrations[j] = self.y0[j] / V # moles per m3 of gas + self._jacobian_inv_omega_num_core_species = self.num_core_species def compute_network_variables(self, pdep_networks=None): # ToDo: this should allow pressure to vary? @@ -467,9 +475,9 @@ cdef class SurfaceReactor(ReactionSystem): surface_volume_ratio_si = self.surface_volume_ratio.value_si C = np.zeros_like(self.core_species_concentrations) - V = self.V # constant volume reactor - A = self.V * surface_volume_ratio_si # area - total_sites = self.surface_site_density.value_si * A # todo: double check units + V = self.V # constant volume reactor in m^3 + A = self.V * surface_volume_ratio_si # area in m^2 + total_sites = self.surface_site_density.value_si * A # moles of sites for j in range(num_core_species): if species_on_surface[j]: C[j] = (N[j] / V) / surface_volume_ratio_si @@ -614,7 +622,7 @@ cdef class SurfaceReactor(ReactionSystem): res = core_species_rates * V # mol/s - if self.sensitivity and False: + if self.sensitivity: delta = np.zeros(len(N), float) delta[:num_core_species] = res if self.jacobian_matrix is None: @@ -633,3 +641,417 @@ cdef class SurfaceReactor(ReactionSystem): # Return DELTA, IRES. IRES is set to 1 in order to tell DASPK to evaluate the sensitivity residuals return delta, 1 + + @cython.boundscheck(False) + def jacobian(self, double t, np.ndarray[np.float64_t, ndim=1] y, np.ndarray[np.float64_t, ndim=1] dydt, + double cj, np.ndarray[np.float64_t, ndim=1] senpar = np.zeros(1, float)): + """ + Return the analytical Jacobian for the reaction system. + + Constant volume (isochoric), variable pressure, with mixed gas-phase + and surface species. + + kf[j] and kr[j] are given in volume units (mol / m^3 / s), even for + surface reactions. generate_rate_coefficients premultiplies the raw + surface rate coefficient by surface_volume_ratio (= A/V, units 1/m) so + that the residual can uniformly compute: + + reaction_rate [mol/m^3/s] = k * prod_r C[r] + dn_i/dt [mol/s] = nu_ij * reaction_rate * V + + for every reaction regardless of phase. + + Concentrations + -------------- + Gas-phase species: C[i] = y[i] / V (mol / m^3) + Surface species: C[i] = y[i] / V / svr = y[i] / A (mol / m^2) + + where svr = surface_volume_ratio (1/m) + + Derivative rule + --------------- + deriv w.r.t. n[s] = k * (prod_{r != s} C[r]) * V / Omega[s] + + For gas-phase species: V / Omega[s] = V / V = 1 + For surface species: V / Omega[s] = V / A = 1 / svr + + No corr terms (like in SimpleReactor) because V and A are both constant. + """ + cdef np.ndarray[np.int_t, ndim=2] ir, ip + cdef np.ndarray[np.float64_t, ndim=1] kf, kr, C + cdef np.ndarray[np.float64_t, ndim=2] pd + cdef np.ndarray[np.float64_t, ndim=1] Omega # per-species normalization + cdef int num_core_reactions, num_core_species, i, j + cdef double k, V, A, deriv + + ir = self.reactant_indices + ip = self.product_indices + + kf = self.kf # already in volume units for all reactions + kr = self.kb # already in volume units for all reactions + num_core_reactions = len(self.core_reaction_rates) + num_core_species = len(self.core_species_concentrations) + + pd = np.zeros((num_core_species, num_core_species), dtype=np.float64) + np.fill_diagonal(pd, -cj) + + V = self.V + A = V * self.surface_volume_ratio.value_si # surface area in m^2 + + # Cache inverse per-species normalization on the instance since + # species_on_surface is static for SurfaceReactor and V/A are constant + inv_omega = self._jacobian_inv_omega + cached_num_core_species = self._jacobian_inv_omega_num_core_species + + if cached_num_core_species != num_core_species: + inv_omega = np.full(num_core_species, 1.0 / V, dtype=np.float64) + for i in range(num_core_species): + if self.species_on_surface[i]: + inv_omega[i] = 1.0 / A + self._jacobian_inv_omega = inv_omega + self._jacobian_inv_omega_num_core_species = num_core_species + + # Concentrations — surface species in mol/m^2, gas in mol/m^3 + C = np.zeros_like(self.core_species_concentrations) + for j in range(num_core_species): + C[j] = y[j] * inv_omega[j] + + for j in range(num_core_reactions): + + # ------------------------------------------------------------------ + # Forward reaction + # ------------------------------------------------------------------ + k = kf[j] + + if ir[j, 1] == -1: # unimolecular forward + # deriv = k * V * inv_omega[s]; gas: k*1 = k; surface: k*V/A = k/svr + deriv = k * V * inv_omega[ir[j, 0]] + pd[ir[j, 0], ir[j, 0]] -= deriv + + pd[ip[j, 0], ir[j, 0]] += deriv + if ip[j, 1] != -1: + pd[ip[j, 1], ir[j, 0]] += deriv + if ip[j, 2] != -1: + pd[ip[j, 2], ir[j, 0]] += deriv + + elif ir[j, 2] == -1: # bimolecular forward + if ir[j, 0] == ir[j, 1]: # A + A -> products + deriv = 2 * k * C[ir[j, 0]] * V * inv_omega[ir[j, 0]] + pd[ir[j, 0], ir[j, 0]] -= 2 * deriv + + pd[ip[j, 0], ir[j, 0]] += deriv + if ip[j, 1] != -1: + pd[ip[j, 1], ir[j, 0]] += deriv + if ip[j, 2] != -1: + pd[ip[j, 2], ir[j, 0]] += deriv + + else: # A + B -> products + # derivative with respect to n[A] + deriv = k * C[ir[j, 1]] * V * inv_omega[ir[j, 0]] + pd[ir[j, 0], ir[j, 0]] -= deriv + pd[ir[j, 1], ir[j, 0]] -= deriv + + pd[ip[j, 0], ir[j, 0]] += deriv + if ip[j, 1] != -1: + pd[ip[j, 1], ir[j, 0]] += deriv + if ip[j, 2] != -1: + pd[ip[j, 2], ir[j, 0]] += deriv + + # derivative with respect to n[B] + deriv = k * C[ir[j, 0]] * V * inv_omega[ir[j, 1]] + pd[ir[j, 0], ir[j, 1]] -= deriv + pd[ir[j, 1], ir[j, 1]] -= deriv + + pd[ip[j, 0], ir[j, 1]] += deriv + if ip[j, 1] != -1: + pd[ip[j, 1], ir[j, 1]] += deriv + if ip[j, 2] != -1: + pd[ip[j, 2], ir[j, 1]] += deriv + + else: # trimolecular forward + + if (ir[j, 0] == ir[j, 1]) and (ir[j, 0] == ir[j, 2]): # A+A+A + deriv = 3 * k * C[ir[j, 0]] * C[ir[j, 0]] * V * inv_omega[ir[j, 0]] + pd[ir[j, 0], ir[j, 0]] -= 3 * deriv + + pd[ip[j, 0], ir[j, 0]] += deriv + if ip[j, 1] != -1: + pd[ip[j, 1], ir[j, 0]] += deriv + if ip[j, 2] != -1: + pd[ip[j, 2], ir[j, 0]] += deriv + + elif ir[j, 0] == ir[j, 1]: # A+A+B + # derivative with respect to n[A] + deriv = 2 * k * C[ir[j, 0]] * C[ir[j, 2]] * V * inv_omega[ir[j, 0]] + pd[ir[j, 0], ir[j, 0]] -= 2 * deriv + pd[ir[j, 2], ir[j, 0]] -= deriv + + pd[ip[j, 0], ir[j, 0]] += deriv + if ip[j, 1] != -1: + pd[ip[j, 1], ir[j, 0]] += deriv + if ip[j, 2] != -1: + pd[ip[j, 2], ir[j, 0]] += deriv + + # derivative with respect to n[B] + deriv = k * C[ir[j, 0]] * C[ir[j, 0]] * V * inv_omega[ir[j, 2]] + pd[ir[j, 0], ir[j, 2]] -= 2 * deriv + pd[ir[j, 2], ir[j, 2]] -= deriv + + pd[ip[j, 0], ir[j, 2]] += deriv + if ip[j, 1] != -1: + pd[ip[j, 1], ir[j, 2]] += deriv + if ip[j, 2] != -1: + pd[ip[j, 2], ir[j, 2]] += deriv + + elif ir[j, 1] == ir[j, 2]: # A+B+B + # derivative with respect to n[A] + deriv = k * C[ir[j, 1]] * C[ir[j, 1]] * V * inv_omega[ir[j, 0]] + pd[ir[j, 0], ir[j, 0]] -= deriv + pd[ir[j, 1], ir[j, 0]] -= 2 * deriv + + pd[ip[j, 0], ir[j, 0]] += deriv + if ip[j, 1] != -1: + pd[ip[j, 1], ir[j, 0]] += deriv + if ip[j, 2] != -1: + pd[ip[j, 2], ir[j, 0]] += deriv + + # derivative with respect to n[B] + deriv = 2 * k * C[ir[j, 0]] * C[ir[j, 1]] * V * inv_omega[ir[j, 1]] + pd[ir[j, 0], ir[j, 1]] -= deriv + pd[ir[j, 1], ir[j, 1]] -= 2 * deriv + + pd[ip[j, 0], ir[j, 1]] += deriv + if ip[j, 1] != -1: + pd[ip[j, 1], ir[j, 1]] += deriv + if ip[j, 2] != -1: + pd[ip[j, 2], ir[j, 1]] += deriv + + elif ir[j, 0] == ir[j, 2]: # A+B+A + # derivative with respect to n[A] + deriv = 2 * k * C[ir[j, 0]] * C[ir[j, 1]] * V * inv_omega[ir[j, 0]] + pd[ir[j, 0], ir[j, 0]] -= 2 * deriv + pd[ir[j, 1], ir[j, 0]] -= deriv + + pd[ip[j, 0], ir[j, 0]] += deriv + if ip[j, 1] != -1: + pd[ip[j, 1], ir[j, 0]] += deriv + if ip[j, 2] != -1: + pd[ip[j, 2], ir[j, 0]] += deriv + + # derivative with respect to n[B] + deriv = k * C[ir[j, 0]] * C[ir[j, 0]] * V * inv_omega[ir[j, 1]] + pd[ir[j, 0], ir[j, 1]] -= 2 * deriv + pd[ir[j, 1], ir[j, 1]] -= deriv + + pd[ip[j, 0], ir[j, 1]] += deriv + if ip[j, 1] != -1: + pd[ip[j, 1], ir[j, 1]] += deriv + if ip[j, 2] != -1: + pd[ip[j, 2], ir[j, 1]] += deriv + + else: # A+B+C, all distinct + # derivative with respect to n[A] + deriv = k * C[ir[j, 1]] * C[ir[j, 2]] * V * inv_omega[ir[j, 0]] + pd[ir[j, 0], ir[j, 0]] -= deriv + pd[ir[j, 1], ir[j, 0]] -= deriv + pd[ir[j, 2], ir[j, 0]] -= deriv + + pd[ip[j, 0], ir[j, 0]] += deriv + if ip[j, 1] != -1: + pd[ip[j, 1], ir[j, 0]] += deriv + if ip[j, 2] != -1: + pd[ip[j, 2], ir[j, 0]] += deriv + + # derivative with respect to n[B] + deriv = k * C[ir[j, 0]] * C[ir[j, 2]] * V * inv_omega[ir[j, 1]] + pd[ir[j, 0], ir[j, 1]] -= deriv + pd[ir[j, 1], ir[j, 1]] -= deriv + pd[ir[j, 2], ir[j, 1]] -= deriv + + pd[ip[j, 0], ir[j, 1]] += deriv + if ip[j, 1] != -1: + pd[ip[j, 1], ir[j, 1]] += deriv + if ip[j, 2] != -1: + pd[ip[j, 2], ir[j, 1]] += deriv + + # derivative with respect to n[C] + deriv = k * C[ir[j, 0]] * C[ir[j, 1]] * V * inv_omega[ir[j, 2]] + pd[ir[j, 0], ir[j, 2]] -= deriv + pd[ir[j, 1], ir[j, 2]] -= deriv + pd[ir[j, 2], ir[j, 2]] -= deriv + + pd[ip[j, 0], ir[j, 2]] += deriv + if ip[j, 1] != -1: + pd[ip[j, 1], ir[j, 2]] += deriv + if ip[j, 2] != -1: + pd[ip[j, 2], ir[j, 2]] += deriv + + # ------------------------------------------------------------------ + # Reverse reaction (ip <-> ir relative to forward block above) + # ------------------------------------------------------------------ + k = kr[j] + + if ip[j, 1] == -1: # unimolecular reverse + deriv = k * V * inv_omega[ip[j, 0]] + pd[ip[j, 0], ip[j, 0]] -= deriv + + pd[ir[j, 0], ip[j, 0]] += deriv + if ir[j, 1] != -1: + pd[ir[j, 1], ip[j, 0]] += deriv + if ir[j, 2] != -1: + pd[ir[j, 2], ip[j, 0]] += deriv + + elif ip[j, 2] == -1: # bimolecular reverse + if ip[j, 0] == ip[j, 1]: # P+P -> reactants + deriv = 2 * k * C[ip[j, 0]] * V * inv_omega[ip[j, 0]] + pd[ip[j, 0], ip[j, 0]] -= 2 * deriv + + pd[ir[j, 0], ip[j, 0]] += deriv + if ir[j, 1] != -1: + pd[ir[j, 1], ip[j, 0]] += deriv + if ir[j, 2] != -1: + pd[ir[j, 2], ip[j, 0]] += deriv + + else: # P1+P2 -> reactants + # derivative with respect to n[P1] + deriv = k * C[ip[j, 1]] * V * inv_omega[ip[j, 0]] + pd[ip[j, 0], ip[j, 0]] -= deriv + pd[ip[j, 1], ip[j, 0]] -= deriv + + pd[ir[j, 0], ip[j, 0]] += deriv + if ir[j, 1] != -1: + pd[ir[j, 1], ip[j, 0]] += deriv + if ir[j, 2] != -1: + pd[ir[j, 2], ip[j, 0]] += deriv + + # derivative with respect to n[P2] + deriv = k * C[ip[j, 0]] * V * inv_omega[ip[j, 1]] + pd[ip[j, 0], ip[j, 1]] -= deriv + pd[ip[j, 1], ip[j, 1]] -= deriv + + pd[ir[j, 0], ip[j, 1]] += deriv + if ir[j, 1] != -1: + pd[ir[j, 1], ip[j, 1]] += deriv + if ir[j, 2] != -1: + pd[ir[j, 2], ip[j, 1]] += deriv + + else: # trimolecular reverse + + if (ip[j, 0] == ip[j, 1]) and (ip[j, 0] == ip[j, 2]): # P+P+P + deriv = 3 * k * C[ip[j, 0]] * C[ip[j, 0]] * V * inv_omega[ip[j, 0]] + pd[ip[j, 0], ip[j, 0]] -= 3 * deriv + + pd[ir[j, 0], ip[j, 0]] += deriv + if ir[j, 1] != -1: + pd[ir[j, 1], ip[j, 0]] += deriv + if ir[j, 2] != -1: + pd[ir[j, 2], ip[j, 0]] += deriv + + elif ip[j, 0] == ip[j, 1]: # P1+P1+P2 + # derivative with respect to n[P1] + deriv = 2 * k * C[ip[j, 0]] * C[ip[j, 2]] * V * inv_omega[ip[j, 0]] + pd[ip[j, 0], ip[j, 0]] -= 2 * deriv + pd[ip[j, 2], ip[j, 0]] -= deriv + + pd[ir[j, 0], ip[j, 0]] += deriv + if ir[j, 1] != -1: + pd[ir[j, 1], ip[j, 0]] += deriv + if ir[j, 2] != -1: + pd[ir[j, 2], ip[j, 0]] += deriv + + # derivative with respect to n[P2] + deriv = k * C[ip[j, 0]] * C[ip[j, 0]] * V * inv_omega[ip[j, 2]] + pd[ip[j, 0], ip[j, 2]] -= 2 * deriv + pd[ip[j, 2], ip[j, 2]] -= deriv + + pd[ir[j, 0], ip[j, 2]] += deriv + if ir[j, 1] != -1: + pd[ir[j, 1], ip[j, 2]] += deriv + if ir[j, 2] != -1: + pd[ir[j, 2], ip[j, 2]] += deriv + + elif ip[j, 1] == ip[j, 2]: # P1+P2+P2 + # derivative with respect to n[P1] + deriv = k * C[ip[j, 1]] * C[ip[j, 1]] * V * inv_omega[ip[j, 0]] + pd[ip[j, 0], ip[j, 0]] -= deriv + pd[ip[j, 1], ip[j, 0]] -= 2 * deriv + + pd[ir[j, 0], ip[j, 0]] += deriv + if ir[j, 1] != -1: + pd[ir[j, 1], ip[j, 0]] += deriv + if ir[j, 2] != -1: + pd[ir[j, 2], ip[j, 0]] += deriv + + # derivative with respect to n[P2] + deriv = 2 * k * C[ip[j, 0]] * C[ip[j, 1]] * V * inv_omega[ip[j, 1]] + pd[ip[j, 0], ip[j, 1]] -= deriv + pd[ip[j, 1], ip[j, 1]] -= 2 * deriv + + pd[ir[j, 0], ip[j, 1]] += deriv + if ir[j, 1] != -1: + pd[ir[j, 1], ip[j, 1]] += deriv + if ir[j, 2] != -1: + pd[ir[j, 2], ip[j, 1]] += deriv + + elif ip[j, 0] == ip[j, 2]: # P1+P2+P1 + # derivative with respect to n[P1] + deriv = 2 * k * C[ip[j, 0]] * C[ip[j, 1]] * V * inv_omega[ip[j, 0]] + pd[ip[j, 0], ip[j, 0]] -= 2 * deriv + pd[ip[j, 1], ip[j, 0]] -= deriv + + pd[ir[j, 0], ip[j, 0]] += deriv + if ir[j, 1] != -1: + pd[ir[j, 1], ip[j, 0]] += deriv + if ir[j, 2] != -1: + pd[ir[j, 2], ip[j, 0]] += deriv + + # derivative with respect to n[P2] + deriv = k * C[ip[j, 0]] * C[ip[j, 0]] * V * inv_omega[ip[j, 1]] + pd[ip[j, 0], ip[j, 1]] -= 2 * deriv + pd[ip[j, 1], ip[j, 1]] -= deriv + + pd[ir[j, 0], ip[j, 1]] += deriv + if ir[j, 1] != -1: + pd[ir[j, 1], ip[j, 1]] += deriv + if ir[j, 2] != -1: + pd[ir[j, 2], ip[j, 1]] += deriv + + else: # P1+P2+P3, all distinct + # derivative with respect to n[P1] + deriv = k * C[ip[j, 1]] * C[ip[j, 2]] * V * inv_omega[ip[j, 0]] + pd[ip[j, 0], ip[j, 0]] -= deriv + pd[ip[j, 1], ip[j, 0]] -= deriv + pd[ip[j, 2], ip[j, 0]] -= deriv + + pd[ir[j, 0], ip[j, 0]] += deriv + if ir[j, 1] != -1: + pd[ir[j, 1], ip[j, 0]] += deriv + if ir[j, 2] != -1: + pd[ir[j, 2], ip[j, 0]] += deriv + + # derivative with respect to n[P2] + deriv = k * C[ip[j, 0]] * C[ip[j, 2]] * V * inv_omega[ip[j, 1]] + pd[ip[j, 0], ip[j, 1]] -= deriv + pd[ip[j, 1], ip[j, 1]] -= deriv + pd[ip[j, 2], ip[j, 1]] -= deriv + + pd[ir[j, 0], ip[j, 1]] += deriv + if ir[j, 1] != -1: + pd[ir[j, 1], ip[j, 1]] += deriv + if ir[j, 2] != -1: + pd[ir[j, 2], ip[j, 1]] += deriv + + # derivative with respect to n[P3] + deriv = k * C[ip[j, 0]] * C[ip[j, 1]] * V * inv_omega[ip[j, 2]] + pd[ip[j, 0], ip[j, 2]] -= deriv + pd[ip[j, 1], ip[j, 2]] -= deriv + pd[ip[j, 2], ip[j, 2]] -= deriv + + pd[ir[j, 0], ip[j, 2]] += deriv + if ir[j, 1] != -1: + pd[ir[j, 1], ip[j, 2]] += deriv + if ir[j, 2] != -1: + pd[ir[j, 2], ip[j, 2]] += deriv + + self.jacobian_matrix = pd + cj * np.identity(num_core_species, float) + return pd diff --git a/test/rmgpy/solver/solverSurfaceTest.py b/test/rmgpy/solver/solverSurfaceTest.py index d4af80264f4..1e484f916f1 100644 --- a/test/rmgpy/solver/solverSurfaceTest.py +++ b/test/rmgpy/solver/solverSurfaceTest.py @@ -39,9 +39,151 @@ from rmgpy.solver.surface import SurfaceReactor from rmgpy.species import Species from rmgpy.thermo import ThermoData, NASA, NASAPolynomial +from rmgpy.solver.termination import TerminationTime +from rmgpy.rmg.settings import ModelSettings, SimulatorSettings +def get_i_thing(thing, thing_list): + for i in range(len(thing_list)): + if thing.is_isomorphic(thing_list[i]): + return i + return -1 + class SurfaceReactorTest: + def setup_class(self): + # Define a simple surface mechanism that we can use for testing + # This is faster than loading in chemkin + self.species_list = [ + Species( + molecule=[Molecule().from_smiles("[Ar]")], + 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"""), + ), + Species( + molecule=[Molecule().from_smiles("[Ne]")], + 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"""), + ), + Species( + molecule=[Molecule().from_smiles("N#N")], + 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"""), + ), + Species( + molecule=[Molecule().from_smiles("C")], + 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"""), + ), + Species( + molecule=[Molecule().from_smiles("[O][O]")], + 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"""), + ), + Species( + molecule=[Molecule().from_smiles("O")], + 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"""), + ), + Species( + molecule=[Molecule().from_smiles("[H][H]")], + 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"""), + ), + Species( + molecule=[Molecule().from_smiles("*")], + 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"""), + ), + Species( + molecule=[Molecule().from_smiles("[H]*")], + 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"""), + ), + Species( + molecule=[Molecule().from_smiles("[O]=*")], + 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"""), + ), + Species( + molecule=[Molecule().from_smiles("[CH3]*")], + 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"""), + ), + Species( + molecule=[Molecule().from_smiles("[H][H].*")], + 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"""), + ), + Species( + molecule=[Molecule().from_smiles("[OH]*")], + 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"""), + ), + Species( + molecule=[Molecule().from_smiles("O.*")], + 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"""), + ), + Species( + molecule=[Molecule().from_smiles("[CH2]=*")], + 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"""), + ), + Species( + molecule=[Molecule().from_smiles("[CH]#*")], + 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"""), + ), + Species( + molecule=[Molecule().from_smiles("[C]$*")], + 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"""), + ) + ] + X = self.species_list[get_i_thing(Species(molecule=[Molecule().from_smiles("*")]), self.species_list)] + O2 = self.species_list[get_i_thing(Species(molecule=[Molecule().from_smiles("[O][O]")]), self.species_list)] + OX = self.species_list[get_i_thing(Species(molecule=[Molecule().from_smiles("O=*")]), self.species_list)] + H2 = self.species_list[get_i_thing(Species(molecule=[Molecule().from_smiles("[H][H]")]), self.species_list)] + HX = self.species_list[get_i_thing(Species(molecule=[Molecule().from_smiles("[H]*")]), self.species_list)] + CH4 = self.species_list[get_i_thing(Species(molecule=[Molecule().from_smiles("C")]), self.species_list)] + CH3X = self.species_list[get_i_thing(Species(molecule=[Molecule().from_smiles("C*")]), self.species_list)] + OHX = self.species_list[get_i_thing(Species(molecule=[Molecule().from_smiles("[OH]*")]), self.species_list)] + H2O = self.species_list[get_i_thing(Species(molecule=[Molecule().from_smiles("O")]), self.species_list)] + H2OX = self.species_list[get_i_thing(Species(molecule=[Molecule().from_smiles("O.*")]), self.species_list)] + CHX = self.species_list[get_i_thing(Species(molecule=[Molecule().from_smiles("[CH]#*")]), self.species_list)] + CH2X = self.species_list[get_i_thing(Species(molecule=[Molecule().from_smiles("[CH2]=*")]), self.species_list)] + + self.reaction_list = [ + Reaction( + reactants=[X, X, O2], + products=[OX, OX], + kinetics=StickingCoefficient(A=0.07, n=0, Ea=(0,'kcal/mol'), T0=(1,'K')) + ), + Reaction( + reactants=[X, X, H2], + products=[HX, HX], + kinetics=StickingCoefficient(A=0.046, n=0, Ea=(0,'kcal/mol'), T0=(1,'K')) + ), + Reaction( + reactants=[HX, CH3X], + products=[X, X, CH4], + kinetics=SurfaceArrhenius(A=(3.3e+21,'cm^2/(mol*s)'), n=0, Ea=(11.95,'kcal/mol'), T0=(1,'K')) + ), + Reaction( + reactants=[X, OHX], + products=[OX, HX], + kinetics=SurfaceArrhenius(A=(7.39e+19,'cm^2/(mol*s)'), n=0, Ea=(18.475,'kcal/mol'), T0=(1,'K')) + ), + Reaction( + reactants=[X, H2O], + products=[H2OX], + kinetics=StickingCoefficient(A=0.75, n=0, Ea=(0,'kcal/mol'), T0=(1,'K')) + ), + Reaction( + reactants=[X, H2OX], + products=[HX, OHX], + kinetics=SurfaceArrhenius(A=(1.15e+19,'cm^2/(mol*s)'), n=0, Ea=(24.235,'kcal/mol'), T0=(1,'K')) + ), + Reaction( + reactants=[OHX, CHX], + products=[OX, CH2X], + kinetics=SurfaceArrhenius(A=(4.4e+22,'cm^2/(mol*s)'), n=0.101, Ea=(10.143,'kcal/mol'), T0=(1,'K')) + ), + Reaction( + reactants=[OHX, CH2X], + products=[OX, CH3X], + kinetics=SurfaceArrhenius(A=(1.39e+21,'cm^2/(mol*s)'), n=0.101, Ea=(4.541,'kcal/mol'), T0=(1,'K')) + ), + Reaction( + reactants=[X, OX, CH4], + products=[OHX, CH3X], + kinetics=SurfaceArrhenius(A=(5e+18,'cm^4/(mol^2*s)'), n=0.7, Ea=(10.038,'kcal/mol'), T0=(1,'K')) + ) + ] + def test_solve_h2(self): """ Test the surface batch reactor with a dissociative adsorption of H2 @@ -1175,3 +1317,82 @@ def test_solve_ch3_thermo_coverage_dependence(self): # Check that coverages are different assert not np.allclose(y, y_off) assert not np.allclose(species_rates, species_rates_off) + + def test_jacobian_with_finite_differences(self): + + # make a simple mechanism + x_O2 = 3e-5 + x_CH4 = 1e-5 + T = (900, 'K') + P = (1, 'atm') + + surface_volume_ratio = (1.0, "m^-1") # TODO, try higher surface volume ratios + surface_site_density = (2.483e-8, "kmol/m^2") # read from Cantera yaml just to be sure these match + termination = TerminationTime((1.0, 's')) + + + Ar = self.species_list[get_i_thing(Species(smiles='[Ar]'), self.species_list)] + CH4 = self.species_list[get_i_thing(Species(smiles='C'), self.species_list)] + CO2 = self.species_list[get_i_thing(Species(smiles='O=C=O'), self.species_list)] + O2 = self.species_list[get_i_thing(Species(smiles='[O][O]'), self.species_list)] + X = self.species_list[get_i_thing(Species(smiles='*'), self.species_list)] + + initial_gas_mole_fractions = {O2: x_O2, CH4: x_CH4, Ar: 1.0 - x_O2 - x_CH4} + initial_surface_coverages = {X: 1.0} + sensitive_species = [CO2] + + reaction_system = SurfaceReactor( + T, + P, + n_sims=1, + initial_gas_mole_fractions=initial_gas_mole_fractions, + initial_surface_coverages=initial_surface_coverages, + surface_volume_ratio=surface_volume_ratio, + surface_site_density=surface_site_density, + termination=[termination], + sensitive_species=sensitive_species, + ) + + # Need to list some sens_worksheets so this doesn't crash when sensitivity is turned on + reaction_system.simulate( + core_species=self.species_list, + core_reactions=self.reaction_list, + edge_species=[], + edge_reactions=[], + surface_species=[], + surface_reactions=[], + model_settings=ModelSettings(tol_move_to_core=1e5), # tol_move_to_core isn't set by default which causes an error + simulator_settings=SimulatorSettings(), # defaults + sensitivity=True, + sens_worksheet=['temp_sensitivity.csv'], + ) + + # now comapare the jacobian to a finite difference approximation of the jacobian for the sensitive species + t = reaction_system.t + y = reaction_system.y + dydt = reaction_system.dydt + + J_analytical = reaction_system.jacobian(t, y, dydt, cj=0.0) + + # turn off sensitivity so that it only handles the first n values of y and doesn't get bogged + # down in all the partial derivatives that get stored in y for sensitivity calculations + reaction_system.sensitivity = False + + n = reaction_system.num_core_species + eps_machine = np.sqrt(np.finfo(float).eps) + n_scale = 1e-10 + + # Get baseline residual f(n) + f0, _ = reaction_system.residual(t, y[:n], dydt[:n]) + + # Build finite difference Jacobian by column + J_fd = np.zeros((n, n)) + for s in range(n): + y_perturbed = y.copy() + delta_n_s = eps_machine * max(abs(y[s]), n_scale) # figure out how much to perturb the value + y_perturbed[s] += delta_n_s # add \Delta n_s + + f_perturbed, _ = reaction_system.residual(t, y_perturbed[:n], dydt[:n]) + J_fd[:, s] = (f_perturbed - f0) / delta_n_s + + assert np.allclose(J_analytical, J_fd, rtol=1e-2, atol=1e-6)