Add get/set equivalenceRatio/mixtureFraction functions #851

g3bk47 · 2020-05-09T15:34:19Z

I accidentally deleted the original pull request. The exact same changes are provided again in this PR. The original discussion can be found here: #811

codecov · 2020-05-09T15:53:31Z

Codecov Report

Merging #851 into main will increase coverage by 0.16%.
The diff coverage is 94.01%.

@@            Coverage Diff             @@
##             main     #851      +/-   ##
==========================================
+ Coverage   70.95%   71.11%   +0.16%     
==========================================
  Files         376      376              
  Lines       45889    46201     +312     
==========================================
+ Hits        32560    32858     +298     
- Misses      13329    13343      +14

Impacted Files	Coverage Δ
include/cantera/thermo/Phase.h	`86.04% <ø> (ø)`
include/cantera/thermo/ThermoPhase.h	`28.28% <ø> (ø)`
src/thermo/Phase.cpp	`82.63% <66.66%> (-1.27%)`	⬇️
src/thermo/ThermoPhase.cpp	`76.43% <94.05%> (+6.08%)`	⬆️
samples/cxx/flamespeed/flamespeed.cpp	`94.59% <100.00%> (-0.24%)`	⬇️
test/thermo/ThermoPhase_Test.cpp	`100.00% <100.00%> (ø)`

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 047e74b...bebb6e2. Read the comment docs.

speth

Thanks for submitting this pull request. In general, I think this will be a useful addition.

One concern I have is the sheer number of new functions introduced here, especially in the C++ interface. What I think would help would be to use an argument to switch between molar and mass basis, rather than having separate function names for each. This could just be an optional basis argument, which would mirror what you've done in the Python interface.

speth · 2020-06-05T20:42:58Z

AUTHORS

@@ -18,7 +18,7 @@ Steven DeCaluwe (@decaluwe), Colorado School of Mines
 Vishesh Devgan (@vdevgan)
 Thomas Fiala (@thomasfiala), Technische Universität München
 David Fronczek
-@g3bk47
+Thorsten Zirwes, Karlsruhe Institute of Technology


We generally include the GitHub handle for contributors who have one, i.e.

Suggested change

Thorsten Zirwes, Karlsruhe Institute of Technology

Thorsten Zirwes (@g3bk47), Karlsruhe Institute of Technology

but it's up to you if you don't want to include it.

speth · 2020-06-05T20:43:27Z

src/thermo/ThermoPhase.cpp

+    auto iC = elementIndex("C");
+    auto iS = elementIndex("S");
+    auto iO = elementIndex("O");
+    auto iH = elementIndex("H");


For very simple types like size_t and double, I would suggest using the type names explicitly rather than auto.

speth · 2020-06-05T20:44:29Z

src/thermo/ThermoPhase.cpp

+    for (size_t i=0; i!=m_kk; ++i)
+    {


A couple of style notes applicable here and elsewhere:

Suggested change

for (size_t i=0; i!=m_kk; ++i)

{

for (size_t i = 0; i != m_kk; ++i) {

The opening brace should be on the same line for if statements as well.

speth · 2020-06-05T20:45:31Z

include/cantera/thermo/ThermoPhase.h

+    //! Set the mixture composition according to the equivalence ratio.
+    //! Fuel and oxidizer compositions are specified by their mass fractions (do
+    //! not need to be normalized). Pressure is kept constant.
+    //! Elements C, S, H and O are considered for the oxidation. 
+    // 
+    //
+    /*!
+     * @param phi        equivalence ratio
+     * @param fuelComp   mass fractions of species in the fuel
+     * @param oxComp     mass fractions of species in the oxidizer 
+     */


To get Doxygen to do the right thing, the formatting for the longer comment blocks like this should be

//! One or two line //! summary. /*! *! (extended description, including parameters) */

speth · 2020-06-05T20:46:29Z

include/cantera/thermo/ThermoPhase.h

+    void setMixtureFraction_X(double mixFrac, const double* fuelComp, const double* oxComp);
+
+    //! @copydoc ThermoPhase::setMixtureFraction_X(double mixFrac, const double* fuelComp, const double* oxComp)
+    void setMixtureFraction_X(double mixFrac, const vector_fp& fuelComp, const vector_fp& oxComp);


Cantera doesn't generally provide functions which take vector_fp& as arguments -- users are expected to use the functions which take simple double*. So I think this overload and the other similar ones should be removed. The string and compsitionMap ones are of course fine.

speth · 2020-06-05T21:39:35Z

src/thermo/ThermoPhase.cpp

+        throw CanteraError("ThermoPhase::getMixtureFraction_Y",
+                "fuel and oxidizer have the same composition");
+
+    auto mf = ( 2.0*(Z_C -Z_Co) + 0.5*(Z_H -Z_Ho) + 2.0*(Z_S -Z_So) - (Z_O -Z_Oo) ) 


Suggested change

auto mf = ( 2.0*(Z_C -Z_Co) + 0.5*(Z_H -Z_Ho) + 2.0*(Z_S -Z_So) - (Z_O -Z_Oo) )

auto mf = (2.0*(Z_C - Z_Co) + 0.5*(Z_H - Z_Ho) + 2.0*(Z_S - Z_So) - (Z_O - Z_Oo))

speth · 2020-06-05T21:41:17Z

src/thermo/ThermoPhase.cpp

+    // mass fractions of fuel and oxidizer
+    vector_fp massFractions(m_kk * 2);


Using a single vector to store both the fuel and oxidizer compositions is an unnecessary complication.

speth · 2020-06-05T21:48:10Z

src/thermo/ThermoPhase.cpp

+    // mean molecular weights of fuel and oxidizer
+    auto rmmw_f = std::inner_product(fuelComp, fuelComp+m_kk, molecularWeights().data(),
+            0.0, std::plus<double>(), std::divides<double>());
+    auto rmmw_o = std::inner_product(oxComp, oxComp+m_kk, molecularWeights().data(), 0.0,
+            std::plus<double>(), std::divides<double>());


I think I would prefer the readability of the straightforward loop for doing this calculation.

speth · 2020-06-05T21:48:23Z

src/thermo/ThermoPhase.cpp

+    auto YtoX = [=](const double& Xk, const double& Mk){return Xk/(rmmw_f*Mk);};
+
+    // compute mole fractions for fuel and oxidizer
+    std::transform(fuelComp, fuelComp+m_kk, molecularWeights().data(), X.data(), YtoX);
+    std::transform(oxComp, oxComp+m_kk, molecularWeights().data(), X.data()+m_kk, YtoX);


I think I would prefer the readability of the straightforward loop for doing this calculation.

speth · 2020-06-05T21:54:50Z

include/cantera/thermo/ThermoPhase.h

+     * @param returns    equivalence ratio 
+     */
+
+    double getEquivalenceRatio() const;


So this is essentially an implementation of the previous definition of getEquivalenceRatio, is that correct? In that case, is it necessary or useful?

This function is the special case of assuming that fuel and oxidizer are pure (what is sometimes called "local equivalence ratio", because it only considers the locally available oxygen and required oxygen and is therefore independent of fuel or oxidizer compositions). I introduced this function for convenience, because you can just write gas.getEquivalenceRatio() instead of specifying a concrete oxidizer and fuel composition like gas.getEquivalenceRatio("H2:1","O2:1"), where the actual values do not matter as long as the compositions only contain fuel and oxidizer elements, respectively.
In the python interface, I solved this by defaulting the oxidizer and fuel compositions to None. However, defaulting the compositionMap& oxComp to compositionMap{} or double* oxComp to nullptr feels a but clunky in the C++ interface. I could rename the function to getLocalEquivalenceRatio, but this might be more confusing. What do you think?

I think what you have here makes sense. Maybe it would help to add a suggestion to use the versions that take fuel and oxidizer compositions if either of these conditions for the fuel or oxidizer isn't met.

g3bk47 · 2020-06-06T13:57:59Z

Thanks for the comments. I will make the requested changes. Just a few questions:

About the behavior='old': This would mean that the default behavior of the python interface is different from default behavior in the C++ interface (I think it doesn't make sense to bring the incorrect "old" python implementation to C++). And should the use of the old behavior be marked as deprecated?
Regarding the bases="mole/mass": Do you want to specify this as a std::string argument or should I introduce an enum to the C++ side similar to

cantera/interfaces/cython/cantera/thermo.pyx

Lines 8 to 10 in 01a4e90

cdef enum ThermoBasis:

mass_basis = 0

molar_basis = 1

I saw that there is already a mixture_fraction function for couterflow flames:

cantera/interfaces/cython/cantera/onedim.py

Lines 1358 to 1378 in 113bd4b

    
               def mixture_fraction(self, m): 
        
                   r""" 
        
                   Compute the mixture fraction based on element *m* 
        
                   The mixture fraction is computed from the elemental mass fraction of 
        
                   element *m*, normalized by its values on the fuel and oxidizer 
        
                   inlets: 
        
                   .. math:: Z = \frac{Z_{\mathrm{mass},m}(z) - 
        
                                       Z_{\mathrm{mass},m}(z_\mathrm{oxidizer})} 
        
                                      {Z_{\mathrm{mass},m}(z_\mathrm{fuel}) - 
        
                                       Z_{\mathrm{mass},m}(z_\mathrm{oxidizer})} 
        
                   :param m: 
        
                       The element based on which the mixture fraction is computed, 
        
                       may be specified by name or by index 
        
                   >>> f.mixture_fraction('H') 
        
                   """ 
        
                   emf = self.elemental_mass_fraction(m) 
        
                   return (emf - emf[-1]) / (emf[0] - emf[-1])

This is based on a single element, whereas the new functions are technically the Bilger mixture fraction based on all elements. I could rename the new functions to BilgerMixtureFraction. Or maybe there could be another argument like element="O", element="C" or element="Bilger" (where "Bilger" is the default), so the user can switch between the single element mode (which might be useful for tracking differential diffusion effects) and the Bilger mode, which considers all elements.

speth · 2020-06-08T15:22:26Z

Thanks for the comments. I will make the requested changes. Just a few questions:

About the behavior='old': This would mean that the default behavior of the python interface is different from default behavior in the C++ interface (I think it doesn't make sense to bring the incorrect "old" python implementation to C++). And should the use of the old behavior be marked as deprecated?

Yes, I think the "old" definition can be deprecated. I agree that bringing the "old" definition into C++ is not worthwhile, so we will just have to deal with the discrepancy for this version.

Regarding the bases="mole/mass": Do you want to specify this as a std::string argument or should I introduce an enum to the C++ side similar to

cantera/interfaces/cython/cantera/thermo.pyx

Lines 8 to 10 in 01a4e90

cdef enum ThermoBasis:

mass_basis = 0

molar_basis = 1

For the C++ interface, probably an enum class would be the most idiomatic, something like

enum class ThermoBasis
{
    mass,
    molar
};

I saw that there is already a mixture_fraction function for couterflow flames:

cantera/interfaces/cython/cantera/onedim.py

Lines 1358 to 1378 in 113bd4b

def mixture_fraction(self, m):

r"""

Compute the mixture fraction based on element *m*

The mixture fraction is computed from the elemental mass fraction of

element *m*, normalized by its values on the fuel and oxidizer

inlets:

.. math:: Z = \frac{Z_{\mathrm{mass},m}(z) -

Z_{\mathrm{mass},m}(z_\mathrm{oxidizer})}

{Z_{\mathrm{mass},m}(z_\mathrm{fuel}) -

Z_{\mathrm{mass},m}(z_\mathrm{oxidizer})}

:param m:

The element based on which the mixture fraction is computed,

may be specified by name or by index

>>> f.mixture_fraction('H')

"""

emf = self.elemental_mass_fraction(m)

return (emf - emf[-1]) / (emf[0] - emf[-1])

This is based on a single element, whereas the new functions are technically the Bilger mixture fraction based on all elements. I could rename the new functions to BilgerMixtureFraction. Or maybe there could be another argument like element="O", element="C" or element="Bilger" (where "Bilger" is the default), so the user can switch between the single element mode (which might be useful for tracking differential diffusion effects) and the Bilger mode, which considers all elements.

That sounds like an interesting generalization, and having the option of getting the Bilger mixture fraction for the counterflow flame would be nice. I like the idea of doing this with an extra, optional argument as you've suggested.

g3bk47 · 2020-06-15T07:56:51Z

I think I addressed all requested changes in the most recent commits. Just one comment:
I changed the name from ThermoBasis to ThermoBasisType for the enum class in the C++ interface because the name clashed with the enum in the python interface.
https://github.com/Cantera/cantera/pull/851/files#diff-6e2333a075896e676562c6cacb247be3R142-R146
and

cantera/interfaces/cython/cantera/thermo.pyx

Lines 8 to 10 in d02f2a7

    
           cdef enum ThermoBasis: 
        
               mass_basis = 0 
        
               molar_basis = 1

One is an enum and the other an enum class, so there is no implicit conversion between them. I guess the ThermoBasis could be changed to use the ThermoBasisType definition so that everything is named ThermoBasis. But I didn't want to change the existing code in this PR.

speth

Thanks for providing the updates, I think they have helped a lot. I have some more comments, but I think most of these are relatively minor issues.

speth · 2020-06-16T20:31:29Z

include/cantera/thermo/ThermoPhase.h

+/*!
+ * @name CONSTANTS - Specification of the input mixture composition 
+ */
+//@{


Despite the error in the existing code, there should only be one block (delimited by @{ and @}) with the name "constants", and this enum should go inside that block.

include/cantera/thermo/ThermoPhase.h

interfaces/cython/cantera/thermo.pyx

src/thermo/ThermoPhase.cpp

test/thermo/ThermoPhase_Test.cpp

g3bk47 · 2020-06-17T00:19:47Z

There are still some inconsistencies, good catch. I will take care of making all requested changes ASAP. To address your questions:

regarding the output of the stoichiometric air fuel ratio
gas.get_stoich_air_fuel_ratio('CH4', 'O2', 'mass') -> 2.0
This is technically a minor inconsistency, which I briefly mentioned before (Add get/set equivalenceRatio/mixtureFraction functions #811 (comment))
The ThermoBasis in the python interface affects either the output of some functions (like density) or the input of the setters (like HP). In the functions I introduced, the ThermoBasis only affects how the inputs (fuel and oxidizer concentrations) are interpreted, never the output. There are basically three functions I introduced:

get equivalence ratio: doesn't the matter if user specifies mass or mole because the concept of equivalence ratio pretty independent of it
get mixture fraction: returns always kg fuel / kg mixture, regardless of ThermoBase. I virtually do not know anyone who uses a mixture fraction in terms of mol fuel / mol mixture so I think a user might expect this behavior, therefore my choice (which makes it technically inconsistent to the python interface).
get stoich air to fuel ratio: always returns mol oxidizer / mol air, regardless of ThermoBase, just because it was convenient for my implementation of the other functions. Therefore, the function above treats the input as 1 kg CH4/1 kg fuel and 1 kg O2 /1 kg oxidizer, but returns the stoichiometric air to fuel ratio as mol oxidizer / mol fuel. If you want, I can change the output when the ThermoBasis is fuel, but then what about the mixture fraction function? Or maybe the name should be changed to something like "molarStoichiometricAirToFuelRatio"? Or alternatively, only the mixture fraction function could be renamed to massMixtureFraction, but I think this might be more confusing.

About the changed output {'CO': 0.14784006249290751, 'NH3': 0.3595664554540104, 'O2': 0.49259348205308207} in the gas.mole_fraction_dict() example: I will double check it.

speth · 2020-06-17T03:27:05Z

Ah yes, I see your point about the stoichiometric air fuel ratio, and I realize now that the documentation does indeed say that it provides the AFR on a molar basis. Like the mixture fraction, though, I think that the mass-based definition is the one more often used.

If you switched this to use the mass-based definition, you would just need to swap mole fractions for mass fractions in the equivalence ratio methods, and they will work the same way.

g3bk47 · 2020-06-20T03:52:21Z

In addition to the requested changes, I introduced two helper functions on the C++ side (o2Required and o2Present). Together with the change, that stoichAirToFuelRatio returns kg oxidizer/kg fuel, this led to significant simplification of the code.
As a nice side effect, all functions now work with mass fractions internally, just like the thermo object itself. And an additional error is detected in some cases, where fuel and oxygen compositions are mixed up.

Regarding the doxygen comment for ThermoBasis: I removed the open block //@{, but left everything else in place.

g3bk47 · 2020-06-20T05:06:42Z

One python test in the CI has failed ("test_ignition_delay_sensitivity").
However, I am unable to reproduce this on my system. When I check the output, I get:

assertNear: -1.26076767538e-09 vs. -1.26487098268e-09 (OK)

When I look at the result of

gas.set_equivalence_ratio(0.4, 'H2', 'O2:1.0, AR:4.0')

I get e.g. for the mass fraction of H2:

y(H2) = 0.00833872772358

When I calculate this by hand, the correct mass fraction should be

y(H2) = 4032/483527 = 0.0083387277235811...

which is pretty much the same. So maybe someone else should double check why this test fails, before merging the PR.

speth · 2020-06-30T22:41:12Z

I just rebased this onto the current main branch, and didn't get any test failures locally. Hopefully, it will pass the checks on the CI servers as well.

bryanwweber · 2020-07-01T01:06:45Z

The failures seem to be on macOS only... I can take a look tomorrow now that the conda packages seem sorted out

g3bk47 · 2020-07-01T10:53:09Z

The test failure persists on macOS. At least, the "wrong" result seems so to be always the same numerical value (across all different python versions and also from the last CI run), which might make debugging easier. Is there a way to get an image of the macOS VM (if it is a VM), so I can help debugging?

bryanwweber · 2020-07-01T14:22:36Z

I can't reproduce this locally on my mac. I'm going to push some debugging statements to this branch to see what the CI is doing.

speth · 2020-07-02T23:56:22Z

I pushed a commit to tighten the tolerances on the sensitivity variables, which at least on my computer gives better agreement between the finite difference and sensitivity coefficient-based calculations in this test (worst relative agreement is around 1e-3, which is much less than the test tolerance of 5e-2). Hopefully this will be enough to pass the CI tests.

The problem seems to stem from the fact that the default tolerances on the sensitivity variables is fairly loose (1e-4), but tightening them has a tendency to lead to integration failures. I suspect this is another reminder that we need to take another look at the details of how we solve the equations for the sensitivities, since it's not as efficient or robust as I think should be possible. But that's obviously well outside the scope of this PR.

bryanwweber · 2020-07-03T01:13:41Z

Thanks @speth! Feel free to drop my commit with the debug statements (either you or @g3bk47).

I suspect this is another reminder that we need to take another look at the details of how we solve the equations for the sensitivities, since it's not as efficient or robust as I think should be possible.

Can you file an issue about this so we don't forget?

The test_ignition_delay_sensitivity test needs tighter intergrator tolerances on the sensitivity coefficients in order to give results that reliably match the finite difference approach on all platforms.

The old version of the sample was using a strange definition for the composition of air which corresponded to neither of the typical approximations, for example "O2: 1.0, N2: 3.76" or "O2: 0.21, N2: 0.79".

speth · 2020-07-05T15:09:54Z

@bryanwweber I documented the issues I've observed with sensitivity analysis in Cantera/enhancements#55.

speth requested changes Jun 5, 2020

View reviewed changes

g3bk47 requested a review from speth June 16, 2020 15:41

speth requested changes Jun 16, 2020

View reviewed changes

g3bk47 requested a review from speth June 20, 2020 05:06

speth force-pushed the GetSetEquilPR branch from a18147c to 3eb5aaf Compare June 30, 2020 22:40

speth changed the base branch from master to main June 30, 2020 23:14

bryanwweber force-pushed the GetSetEquilPR branch from 5de2dc3 to f41fdaf Compare July 1, 2020 15:20

np4075@kit.edu and others added 5 commits July 4, 2020 12:27

Add get/set equivalenceRatio/mixtureFraction functions

d453581

Add python example for get/set equivalence ratio

20fd418

MixtureFraction/EquilvalenceRatio: add basis paramter and element option

db5fbd0

Simplify computation of equivalence ratio and mixture fraction

d20ce60

[Test] Tighten integration tolerances to get repeatable test results

e8a6fac

The test_ignition_delay_sensitivity test needs tighter intergrator tolerances on the sensitivity coefficients in order to give results that reliably match the finite difference approach on all platforms.

speth force-pushed the GetSetEquilPR branch from 742b010 to e8a6fac Compare July 4, 2020 16:28

[Test] Update output files for C++ flamespeed sample

bebb6e2

The old version of the sample was using a strange definition for the composition of air which corresponded to neither of the typical approximations, for example "O2: 1.0, N2: 3.76" or "O2: 0.21, N2: 0.79".

speth approved these changes Jul 4, 2020

View reviewed changes

speth merged commit 07dac7b into Cantera:main Jul 5, 2020

g3bk47 mentioned this pull request Jul 6, 2020

Get/Set Equivalence Ratio/Mixture Fraction Functions for C++ Cantera/enhancements#28

Closed

speth mentioned this pull request Jan 26, 2021

test_reactor.TestReactorSensitivities fails on Gentoo Linux #937

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add get/set equivalenceRatio/mixtureFraction functions #851

Add get/set equivalenceRatio/mixtureFraction functions #851

g3bk47 commented May 9, 2020

codecov bot commented May 9, 2020 •

edited

Loading

speth left a comment

speth Jun 5, 2020

speth Jun 5, 2020

speth Jun 5, 2020

speth Jun 5, 2020

speth Jun 5, 2020

speth Jun 5, 2020

speth Jun 5, 2020

speth Jun 5, 2020

speth Jun 5, 2020

speth Jun 5, 2020

g3bk47 Jun 6, 2020

speth Jun 16, 2020

g3bk47 commented Jun 6, 2020

speth commented Jun 8, 2020

g3bk47 commented Jun 15, 2020

speth left a comment

speth Jun 16, 2020

g3bk47 commented Jun 17, 2020

speth commented Jun 17, 2020

g3bk47 commented Jun 20, 2020 •

edited

Loading

g3bk47 commented Jun 20, 2020 •

edited

Loading

speth commented Jun 30, 2020

bryanwweber commented Jul 1, 2020

g3bk47 commented Jul 1, 2020

bryanwweber commented Jul 1, 2020

speth commented Jul 2, 2020

bryanwweber commented Jul 3, 2020

speth commented Jul 5, 2020

	Thorsten Zirwes, Karlsruhe Institute of Technology
	Thorsten Zirwes (@g3bk47), Karlsruhe Institute of Technology

	for (size_t i=0; i!=m_kk; ++i)
	{
	for (size_t i = 0; i != m_kk; ++i) {

	auto mf = ( 2.0(Z_C -Z_Co) + 0.5(Z_H -Z_Ho) + 2.0*(Z_S -Z_So) - (Z_O -Z_Oo) )
	auto mf = (2.0(Z_C - Z_Co) + 0.5(Z_H - Z_Ho) + 2.0*(Z_S - Z_So) - (Z_O - Z_Oo))

		// mass fractions of fuel and oxidizer
		vector_fp massFractions(m_kk * 2);

Add get/set equivalenceRatio/mixtureFraction functions #851

Add get/set equivalenceRatio/mixtureFraction functions #851

Conversation

g3bk47 commented May 9, 2020

codecov bot commented May 9, 2020 • edited Loading

Codecov Report

speth left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

g3bk47 commented Jun 6, 2020

speth commented Jun 8, 2020

g3bk47 commented Jun 15, 2020

speth left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

g3bk47 commented Jun 17, 2020

speth commented Jun 17, 2020

g3bk47 commented Jun 20, 2020 • edited Loading

g3bk47 commented Jun 20, 2020 • edited Loading

speth commented Jun 30, 2020

bryanwweber commented Jul 1, 2020

g3bk47 commented Jul 1, 2020

bryanwweber commented Jul 1, 2020

speth commented Jul 2, 2020

bryanwweber commented Jul 3, 2020

speth commented Jul 5, 2020

codecov bot commented May 9, 2020 •

edited

Loading

g3bk47 commented Jun 20, 2020 •

edited

Loading

g3bk47 commented Jun 20, 2020 •

edited

Loading