Extensible jacobian and network jacobian interfaces #1634
Conversation
anthony-walker
force-pushed
the
ext_prec
branch
from
eb36c56
to
205567a
Compare
Codecov ReportAttention: Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #1634 +/- ##
==========================================
+ Coverage 72.80% 72.86% +0.06%
==========================================
Files 381 382 +1
Lines 54020 54204 +184
Branches 9210 9238 +28
==========================================
+ Hits 39328 39495 +167
- Misses 11707 11714 +7
- Partials 2985 2995 +10 ☔ View full report in Codecov by Sentry. |
anthony-walker
force-pushed
the
ext_prec
branch
3 times, most recently
from
1c005f9
to
ab83aa3
Compare
There was a problem hiding this comment.
While I don't have enough experience in extensible things for a complete review, I noticed that a lot of test are failing with:
*******************************************************************************
NotImplementedError thrown by MassFlowController::buildReactorJacobian:
Not implemented.
*******************************************************************************
which doesn't appear to be .NET related.
anthony-walker
force-pushed
the
ext_prec
branch
2 times, most recently
from
982eb65
to
96bd8af
Compare
|
Before you get any further, I just wanted to suggest undoing the merge commit and rebasing onto the current |
anthony-walker
force-pushed
the
ext_prec
branch
from
a78547c
to
6f61796
Compare
anthony-walker
force-pushed
the
ext_prec
branch
2 times, most recently
from
4857e53
to
84ff012
Compare
anthony-walker
force-pushed
the
ext_prec
branch
from
decbbac
to
feb862c
Compare
anthony-walker
force-pushed
the
ext_prec
branch
from
feb862c
to
a525724
Compare
anthony-walker
force-pushed
the
ext_prec
branch
from
a525724
to
9a2d1b2
Compare
There was a problem hiding this comment.
Thanks for putting this together, @anthony-walker. I think it helps in laying out the necessary pieces for handling the cross-reactor terms in the Jacobian -- in particular, the mechanics surrounding Wall::buildNetworkJacobian seem like a good basis to keep building on.
I think there is some room for simplification of this structure, which I've left specific comments on. However, I think there is a significant issue with sorting out exactly which derivatives need to be computed and what the right formulas are for the different cases. The challenge I see us needing to resolve to generalize this across the different reactor types is that what derivative is needed for a single interaction (for example, wall heat transfer) depends on both reactor types and what variables they use to represent each term, and we need a way of selecting the correct formula for each Jacobian element based on that pair of reactor types.
| // derivative of temperature transformed by ideal gas law | ||
| vector<double> moles(m_nsp); | ||
| getMoles(moles.data()); | ||
| double dTdni = pressure() * m_vol / GasConstant / std::accumulate(moles.begin(), moles.end(), 0.0); |
There was a problem hiding this comment.
Like in IdealGasMoleReactor,
test/python/test_reactor.py
Outdated
| # check for values | ||
| for i in range(gas1.n_species): | ||
| assert np.isclose(jac[0, i + 1 + r1.n_vars], - U * A * gas2.T) | ||
| assert np.isclose(jac[r1.n_vars, i + 1], U * A * gas1.T) | ||
| if (i + 1 != 2): | ||
| assert np.isclose(jac[0, i + 1], U * A * gas1.T, rtol=1e-2) | ||
| assert np.isclose(jac[r1.n_vars, i + 1 + r1.n_vars], - U * A * gas2.T, rtol=1e-2) |
There was a problem hiding this comment.
I think you should check against the finite difference Jacobian. This test is just checking against the same incorrect formulas in Python as what's in the C++ implementation.
test/python/test_reactor.py
Outdated
| gas2 = ct.Solution("h2o2.yaml", "ohmech") | ||
| gas2.TPX = 300, ct.one_atm, "O2:1.0" | ||
| r2 = self.reactorClass(gas2) | ||
| r2.chemistry_enabled = False |
There was a problem hiding this comment.
The analytical Jacobian doesn't currently respect the chemistry_enabled flag, so there are a lot of non-zero elements here that should be zero based on the setting of this flag.
|
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: |
|
@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. |
It can sometimes be useful to add on to the default evaluation in such a way that the standard before or after delegations cannot handle. A default evaluation function allows for this.
in regards to equations, handling volumes, etc.
|
@speth pushed up some more changes, I will have to put a bit of thought into some of the other changes. |
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
ExtensibleReactorsystem 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
scons build&scons test) and unit tests address code coverage