Selecting a subset of species from a `SolutionArray` object does not work as expected when extracting `net_production_rates`

[dylanrubini]

Problem description

Selecting a subset of species from a SolutionArray object does not work as expected when extracting the net_production_rates

Steps to reproduce

Consider the following Minimum Working Example

import cantera as ct

gas = ct.Solution("./Kinetic-Model/chem_annotated-propane.yaml")
gas.TPY = 900, 1.0e5, {"propane(1)": 1.0}

state = ct.SolutionArray(gas, shape=10)

print(f"Y ---> Shape = {state('propane(1)').Y.shape}")
print(f"w_dot ----> Shape = {state('propane(1)').net_production_rates.shape}")

Behavior

The first print statement prints:

Y species Shape = (10, 1)

as expected because we intended to extract 1 species, in this case propane(1)

However, the second print statement gives:

w_dot Shape = (10, 65)

This is the shape of the entire net_production_rates array including all the species (which totals to 65 in this case).

I expected that extracting a subset of species state('propane(1)') would apply for all properties that have dimension equal to the number of species, not just for the species mass fractions Y and mole fractions X.

Is this the expected behaviour?

System information

• Cantera version: 2.6.0
• OS: CentOS7
• Python/MATLAB/other software versions: Python

[dylanrubini]

It is also worth noting that the same behaviour occurs in Cantera 3.0.0

[ischoegl]

Cantera 3.0 reimplements the 2.6 interface, so it is a pre-existing condition. I'm not sure there ever was a rationale - at the back-end, it is a slicing operation that appears to not be implemented.

[dylanrubini]

Thank you, I appreciate the clarification

[ischoegl]

Fwiw, the implementation should follow the rationale of the underlying Solution object:

In [1]: import cantera as ct
   ...: gas = ct.Solution("h2o2.yaml")
   ...: gas.TPY = 900, 1.e05, {"H2": 1.0, "O2": 1.0}

In [2]: gas["H2"].Y
Out[2]: array([0.5])

In [3]: gas["H2"].net_production_rates
Out[3]: array([-2.58564267e-08])

[speth]

I think the underlying issue here comes from an incomplete attempt to implement SolutionArray for surface phases. For surfaces, the number of species that the Interface object computes kinetics-related properties for (n_total_species) is more than just the number of surface phase species (n_species). Besides the slicing issue noted here, there are other problems with using SolutionArray on Interface objects that I think are tricky to fix (most importantly, reaction rates depend on composition information from multiple phases, and the state of those Solution objects is not synchronized when looping over the SolutionArray's states).

For the time being, we might be best off disabling the use of SolutionArray with Interface.

[ischoegl]

For the time being, we might be best off disabling the use of SolutionArray with Interface

Unfortunately, this would break serialization of 1D simulation results (example: catalytic_combustion.py)

As an alternative, I disabled problematic methods in an update to #1726. Fwiw, I believe the lack of synchronization for Interface should be tracked/fixed as a completely separate issue.