Skip to content

Commit

Permalink
[unittest] Add Kinetics derivative tests
Browse files Browse the repository at this point in the history
Unit tests match example in Kinetics derivative documentation.
  • Loading branch information
ischoegl authored and speth committed Aug 4, 2023
1 parent 6ba497c commit ceec659
Showing 1 changed file with 51 additions and 0 deletions.
51 changes: 51 additions & 0 deletions test/python/test_jacobian.py
Expand Up @@ -368,6 +368,48 @@ def test_net_rop_ddP(self):
drop_num = self.rop_ddP(mode="net")
self.assertNear(drop[self.rxn_idx], drop_num, self.rtol)

def rate_ddT(self, mode=None, const_p=False, rtol=1e-6):
# numerical derivative for production rates with respect to temperature
def calc():
if mode == "creation":
return self.gas.creation_rates
if mode == "destruction":
return self.gas.destruction_rates
if mode == "net":
return self.gas.net_production_rates

dt = self.tpx[0] * rtol
dp = 0 if const_p else self.tpx[1] * rtol
self.gas.TP = self.tpx[0] + dt, self.tpx[1] + dp
rate1 = calc()
self.gas.TP = self.tpx[:2]
rate0 = calc()
return (rate1 - rate0) / dt

def test_net_rate_ddT(self):
# check equivalence of numerical and analytical derivatives of net creation
# rates with respect to temperature

# constant pressure - need to account for density change
# numeric: d(omegadot)/dT =
# analytic: d(omegadot)/dT + dC/dT d(omegadot)/dC
dcdt = - self.gas.density_mole / self.gas.T
drate = self.gas.net_production_rates_ddT
drate += self.gas.net_production_rates_ddC * dcdt
drate_num = self.rate_ddT(mode="net", const_p=True)
for spc_ix in self.rix + self.pix:
assert drate[spc_ix] == pytest.approx(drate_num[spc_ix], self.rtol)

# constant density (volume) - need to account for pressure change
# numeric: d(omegadot)/dT =
# analytic: d(omegadot)/dT + dP/dT d(omegadot)/dP
dpdt = self.gas.P / self.gas.T
drate = self.gas.net_production_rates_ddT
drate += self.gas.net_production_rates_ddP * dpdt
drate_num = self.rate_ddT(mode="creation") - self.rate_ddT(mode="destruction")
for spc_ix in self.rix + self.pix:
assert drate[spc_ix] == pytest.approx(drate_num[spc_ix], self.rtol)

def rate_ddX(self, spc_ix, mode=None, const_t=True, rtol_deltac=1e-6,
atol_deltac=1e-20, ddX=True):
# numerical derivative for production rates with respect to mole fractions
Expand Down Expand Up @@ -603,6 +645,10 @@ def test_reverse_rop_ddT(self):
def test_net_rop_ddT(self):
super().test_net_rop_ddT()

@pytest.mark.usefixtures("has_temperature_derivative_warnings")
def test_net_rate_ddT(self):
super().test_net_rate_ddT()


class TestPlog(FromScratchCases, utilities.CanteraTest):
# Plog reaction
Expand Down Expand Up @@ -639,6 +685,11 @@ def test_reverse_rop_ddT(self):
def test_net_rop_ddT(self):
super().test_net_rop_ddT()

@pytest.mark.xfail(reason="Change of reaction enthalpy is not considered")
@pytest.mark.filterwarnings("ignore:.*does not consider.*(electron|enthalpy).*:UserWarning")
def test_net_rate_ddT(self):
super().test_net_rate_ddT()


class FullTests:
# Generic test class to check derivatives evaluated for an entire reaction mechanisms
Expand Down

0 comments on commit ceec659

Please sign in to comment.