Cantera / cantera

Chemical kinetics, thermodynamics, and transport tool suite
https://cantera.org
Other
581 stars 342 forks source link

Extensible jacobian and network jacobian interfaces #1634

Open anthony-walker opened 8 months ago

anthony-walker commented 8 months ago

Changes proposed in this pull request

In this pull request I have setup basic structure for adding contributions to the jacobian from flow devices and walls. I have also added an extensible jacobian interface to the ExtensibleReactor system and modified the existing jacobian system to pass a vector of triplets as is needed in lieu of creating a large number of separate sparse matrices.

@speth @ischoegl some additional tests are needed for the wall and flow contributions but I wanted to get the content officially into a PR.

Checklist

codecov[bot] commented 8 months ago

Codecov Report

Attention: Patch coverage is 82.24638% with 49 lines in your changes are missing coverage. Please review.

Project coverage is 75.74%. Comparing base (06b3500) to head (3193fa3).

Files Patch % Lines
src/zeroD/ReactorNet.cpp 81.25% 7 Missing and 5 partials :warning:
include/cantera/zeroD/ReactorBase.h 25.00% 6 Missing :warning:
include/cantera/zeroD/Wall.h 33.33% 4 Missing :warning:
src/zeroD/Reactor.cpp 66.66% 2 Missing and 2 partials :warning:
include/cantera/base/Delegator.h 62.50% 2 Missing and 1 partial :warning:
include/cantera/zeroD/ReactorNet.h 75.00% 2 Missing and 1 partial :warning:
interfaces/cython/cantera/delegator.pyx 75.00% 3 Missing :warning:
src/zeroD/Wall.cpp 91.89% 0 Missing and 3 partials :warning:
test/zeroD/test_zeroD.cpp 95.08% 0 Missing and 3 partials :warning:
include/cantera/zeroD/Reactor.h 71.42% 2 Missing :warning:
... and 4 more
Additional details and impacted files ```diff @@ Coverage Diff @@ ## main #1634 +/- ## ========================================== + Coverage 75.69% 75.74% +0.05% ========================================== Files 443 444 +1 Lines 60970 61226 +256 Branches 9551 9589 +38 ========================================== + Hits 46153 46378 +225 - Misses 11786 11803 +17 - Partials 3031 3045 +14 ```

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

anthony-walker commented 8 months ago

@speth @ischoegl I am not entirely sure why the .NET tests are failing but it doesn't seem to be with this PR, is this a known issue?

Otherwise I believe this code is ready for review.

anthony-walker commented 8 months ago

@ischoegl these failures are from my latest bug fixes there is a separate commit which addresses the .NET failure issue, 0a34092.

speth commented 6 months ago

Before you get any further, I just wanted to suggest undoing the merge commit and rebasing onto the current main branch instead.

speth commented 3 months ago

For an example of where the current derivatives are clearly not correct, consider the following example, adapted from one of your new tests:

import cantera as ct
import numpy as np

# create first reactor
gas1 = ct.Solution("h2o2.yaml", "ohmech")
gas1.TPX = 600, ct.one_atm, "O2:1.0"
r1 = ct.IdealGasMoleReactor(gas1)
r1.name = 'R1'
gas1.set_multiplier(0.0)

# create second reactor
gas2 = ct.Solution("h2o2.yaml", "ohmech")
gas2.TPX = 300, ct.one_atm, "O2:1.0"
gas2.set_multiplier(0.0)
r2 = ct.IdealGasMoleReactor(gas2)
r2.name = 'R2'

# create wall
U = 500.0
A = 3.0
w = ct.Wall(r1, r2, U=U, A=A)
net = ct.ReactorNet([r1, r2])

def print_jac(J):
    for i,j in zip(*np.nonzero(J)):
        name_i = net.component_name(i).replace('temperature', 'T')
        name_j = net.component_name(j).replace('temperature', 'T')
        net.component_name(i).replace('temperature', 'T')
        #print(f"J[{i: 3d}, {j: 3d}] = {J[i,j]:.3e}")
        print(f"J[{name_i:8s}, {name_j:8s}] = {J[i,j]:.3e}")

print('Semi-analytical Jacobian:')
print_jac(net.jacobian)
print('\nFinite difference Jacobian:')
print_jac(net.finite_difference_jacobian)

which outputs:

Semi-analytical Jacobian:
J[R1: T   , R1: T   ] = -2.727e+00
J[R1: T   , R1: H2  ] = 9.000e+05
J[R1: T   , R1: H   ] = 9.000e+05
J[R1: T   , R1: O   ] = 9.000e+05
J[R1: T   , R1: O2  ] = 9.000e+05
J[R1: T   , R1: OH  ] = 9.000e+05
J[R1: T   , R1: H2O ] = 9.000e+05
J[R1: T   , R1: HO2 ] = 9.000e+05
J[R1: T   , R1: H2O2] = 9.000e+05
J[R1: T   , R1: AR  ] = 9.000e+05
J[R1: T   , R1: N2  ] = 9.000e+05
J[R1: T   , R2: H2  ] = -4.500e+05
J[R1: T   , R2: H   ] = -4.500e+05
J[R1: T   , R2: O   ] = -4.500e+05
J[R1: T   , R2: O2  ] = -4.500e+05
J[R1: T   , R2: OH  ] = -4.500e+05
J[R1: T   , R2: H2O ] = -4.500e+05
J[R1: T   , R2: HO2 ] = -4.500e+05
J[R1: T   , R2: H2O2] = -4.500e+05
J[R1: T   , R2: AR  ] = -4.500e+05
J[R1: T   , R2: N2  ] = -4.500e+05
J[R2: T   , R1: H2  ] = 9.000e+05
J[R2: T   , R1: H   ] = 9.000e+05
J[R2: T   , R1: O   ] = 9.000e+05
J[R2: T   , R1: O2  ] = 9.000e+05
J[R2: T   , R1: OH  ] = 9.000e+05
J[R2: T   , R1: H2O ] = 9.000e+05
J[R2: T   , R1: HO2 ] = 9.000e+05
J[R2: T   , R1: H2O2] = 9.000e+05
J[R2: T   , R1: AR  ] = 9.000e+05
J[R2: T   , R1: N2  ] = 9.000e+05
J[R2: T   , R2: T   ] = -1.887e+00
J[R2: T   , R2: H2  ] = -4.500e+05
J[R2: T   , R2: H   ] = -4.500e+05
J[R2: T   , R2: O   ] = -4.500e+05
J[R2: T   , R2: O2  ] = -4.500e+05
J[R2: T   , R2: OH  ] = -4.500e+05
J[R2: T   , R2: H2O ] = -4.500e+05
J[R2: T   , R2: HO2 ] = -4.500e+05
J[R2: T   , R2: H2O2] = -4.500e+05
J[R2: T   , R2: AR  ] = -4.500e+05
J[R2: T   , R2: N2  ] = -4.500e+05

Finite difference Jacobian:
J[R1: T   , R1: T   ] = -2.727e+00
J[R1: T   , R1: O2  ] = 4.589e+04
J[R1: T   , R2: T   ] = 3.107e+00
J[R2: T   , R1: T   ] = 1.752e+00
J[R2: T   , R2: T   ] = -1.887e+00
J[R2: T   , R2: O2  ] = -1.294e+04
anthony-walker commented 3 weeks ago

@speth apologies for the delay, I have made some updates but still need to do some of the larger changes. It is still on my radar.