Enhance the PlasmaPhase to use the cross section data with an example of calculating the elastic collision energy loss rate #1731

BangShiuh · 2024-07-08T13:33:05Z

Changes proposed in this pull request
This PR change the initiation of the PlasmaPhase. The PlasmaPhase after this change takes the ElectronCollisionPlasmaRate cross section data to calculate the plasma thermal properties (specifically the elastic collision energy loss rate). The elastic collision energy loss is the energy transfer from electron to molecule when a fast electron (higher electron temperature) collides to a slower molecule (lower gas temperature). This is useful for modelling a non-equilibrium plasma which electron temperature is much higher than the gas temperature. The elastic collision reaction does not produce a new species, but can heat up the gas: e + O2 => e + O2 + heat. We check all the ElectronCollisionPlasmaReaction to find every reaction which the reactants is equal to the products in the initTherm method in PlasmaPhase. We then calculate the elastic collision energy loss using the cross section data from those reactions we found. The result is compared to BOLOS https://github.com/aluque/bolos. The grid definition and finite difference is different. The result is close when the number of grid points is high. The definition of the electron energy distribution function is different in Cantera and BOLOS. To compare the result, I made a modification to BOLOS, https://github.com/BangShiuh/bolos/tree/mod. You can reproduce the result using the scripts and data in
scripts.zip.

Future work to utilize "effective" cross section:
BOLOS calculates the elastic cross section subtracting the inelastic cross section from the effective cross section. This can be included in the parsing script, lxcat2yaml.py to obtain the elastic cross section. I think it is much simpler to use an elastic cross section in a yaml input file as it does not need further transformation like an effective cross section.

Design philosophy: the original cross section data has the collision type such as attachment, ionization, elastic,... etc. However, the collision type can be determined from the collision equation. For example, elastic collisions have the same reactants and products. Ionization collisions have an extra electron in the products. Cantera manages the reaction equations well and we should use it to our advantage to simplify the input data format.

If applicable, fill in the issue number this pull request is fixing

Closes #

If applicable, provide an example illustrating new features this pull request is introducing

Checklist

The pull request includes a clear description of this code change
Commit messages have short titles and reference relevant issues
Build passes (scons build & scons test) and unit tests address code coverage
Style & formatting of contributed code follows contributing guidelines
The pull request is ready for review

…a phase and elasticElectronEnergyLossRate

codecov · 2024-07-08T13:58:37Z

Codecov Report

Attention: Patch coverage is 74.40000% with 32 lines in your changes missing coverage. Please review.

Project coverage is 75.71%. Comparing base (0c89553) to head (7d866c9).
Report is 397 commits behind head on main.

Files with missing lines	Patch %	Lines
src/thermo/PlasmaPhase.cpp	72.60%	13 Missing and 7 partials ⚠️
src/kinetics/KineticsFactory.cpp	78.57%	3 Missing and 3 partials ⚠️
src/kinetics/ElectronCollisionPlasmaRate.cpp	0.00%	3 Missing and 1 partial ⚠️
src/kinetics/Kinetics.cpp	0.00%	1 Missing and 1 partial ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #1731      +/-   ##
==========================================
+ Coverage   72.76%   75.71%   +2.94%     
==========================================
  Files         378      445      +67     
  Lines       56986    61622    +4636     
  Branches    20691     9675   -11016     
==========================================
+ Hits        41468    46658    +5190     
+ Misses      12463    11894     -569     
- Partials     3055     3070      +15

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

speth

Thanks for working on this, @BangShiuh. I had some suggestions that I think will lead to a simpler implementation. Please let me know if you have any questions.

speth · 2024-07-21T20:52:33Z

interfaces/cython/cantera/thermo.pyx

@@ -1905,6 +1905,24 @@ cdef class ThermoPhase(_SolutionBase):
                raise ThermoModelMethodError(self.thermo_model)
            return pystr(self.plasma.electronSpeciesName())

+    property elastic_electron_energy_loss_rate:
+        """ Elastic electron energy loss rate """


What are the units of this quantity?

speth · 2024-07-21T20:57:31Z

interfaces/cython/cantera/thermo.pxd

+        double elasticElectronEnergyLossRate()
+        double normalizedElasticElectronEnergyLossRate()


These methods should probably have the except +translate_exception specifier -- generally recommended except for trivial getters that can't throw.

speth · 2024-07-21T20:59:27Z

interfaces/cython/cantera/thermo.pyx

+        Normalized elastic electron energy loss rate
+        The elastic electron energy loss rate is normalized
+        by dividing the concentration of electron.


Careful with the formatting here -- the second line will be seen as a continuation of the same sentence by Sphinx.

speth · 2024-07-21T21:30:00Z

include/cantera/thermo/PlasmaPhase.h

+    //! Get elastic electron energy loss rate (eV/s)
+    double elasticElectronEnergyLossRate() {
+        return concentration(m_electronSpeciesIndex) *
+               normalizedElasticElectronEnergyLossRate();
+    }
+
+    //! Get normalized elastic electron energy loss rate (eV-m3/kmol/s)
+    double normalizedElasticElectronEnergyLossRate();
+


Do we need both of these methods, given the relatively simple conversion? Which one is more likely to be used? My guess would be the non-normalized version, but I could easily be convinced otherwise. Also, are you sure about the units for the first method? I think this should be an intrinsic property, so the dimensions should be power per unit mass (or possibly power per unit volume).

Also, can you write out the equation being evaluated, and say something about how that's being done? (since as you note in the PR, the finite difference scheme chosen matters)

speth · 2024-07-21T21:38:01Z

test/python/test_thermo.py

+        self.phase.TPX = 1000, ct.one_atm, "O2:1, E:1e-10"
+        r1 = self.phase.elastic_electron_energy_loss_rate
+        c_electron = self.phase.concentrations[self.phase.species_index("E")]
+        self.assertNear(self.phase.elastic_electron_energy_loss_rate, 4.233683742587e-05)


We've been shifting away from the unittest style assert methods to use plain assert statements that are handled by pytest. Here, you can write:

Suggested change

self.assertNear(self.phase.elastic_electron_energy_loss_rate, 4.233683742587e-05)

assert self.phase.elastic_electron_energy_loss_rate == approx(4.233683742587e-05)

speth · 2024-07-21T22:39:05Z

src/thermo/PlasmaPhase.cpp

+    vector<shared_ptr<Reaction>> reactions;
+    for (AnyMap R : reactionsAnyMapList(*m_kinetics, m_input, m_root)) {
+        reactions.push_back(newReaction(R, *m_kinetics));
+    }
+
+    // add reactions to kinetics object
+    addReactions(*m_kinetics, reactions);
+
+    // init m_collisions
+    m_collisions.resize(0);
+    size_t i = 0;
+    for (shared_ptr<Reaction> reaction : reactions) {
+        if (reaction->type() == "electron-collision-plasma") {
+            m_collisions.push_back(reaction);
+            if (reaction->reactants == reaction->products) {


Instead of building a standalone list of Reactions, you could just call the existing method addReactions(*m_kinetics, m_input, m_root), and then use Kinetics.reaction(i) to loop over the reaction objects and build up m_collisions. Then, you wouldn't need the two new functions reactionsAnyMapList and addReactions(Kinetics&, vector<shared_ptr<Reaction>>).

speth · 2024-07-21T22:50:01Z

include/cantera/thermo/ThermoPhase.h

+    //! The kinetics object associates with ThermoPhase
+    //! Some phase requires Kinetics to perform calculation
+    //! such as PlasmaPhase
+    shared_ptr<Kinetics> m_kinetics;


Because the Kinetics object also holds a shared_ptr<ThermoPhase>, this will create a reference cycle and introduce a memory leak. The solution to this problem is to only hold a weak_ptr in one direction. I'd recommend doing this by having the ThermoPhase object hold a weak_ptr to the Solution (say, m_soln), which will also eliminate the need to inherit from enable_shared_from_this, as you can access the shared_ptr<ThermoPhase> from the Solution.

speth · 2024-07-21T22:51:16Z

include/cantera/thermo/ThermoPhase.h

+    shared_ptr<Kinetics> kinetics() {
+        return m_kinetics;
+    }


Method needs a docstring.

speth · 2024-07-21T22:57:08Z

src/base/Solution.cpp

+    if (!sol->thermo()->kinetics()) {
+        // only add phases that doesn't have kinetics
+        // The phase that has kinetics already
+        phases.push_back(sol->thermo());


This comment could be worded more clearly. Also, I'm not sure I follow the behavior here. Shouldn't the behavior be to skip calling newKinetics if the phase already has a Kinetics object? The logic about adjacent phases only applies to surfaces/edges in any case.

speth · 2024-07-21T23:09:05Z

src/thermo/PlasmaPhase.cpp

+    m_kinetics = newKinetics("bulk");
+    m_kinetics->addThermo(shared_from_this());


I think a lot of this logic might be simpler if you leave initialization of the Kinetics object to the newSolution factory function as it originally stood. The way to connect the Kinetics object to the PlasmaPhase would then be to add here in initThermo a callback to the Solution object (see Solution.registerChangedCallback) that would then update m_collisions when it sees a change in the Kinetics object. You'll also need to call removeChangedCallback from the PlasmaPhase destructor.

BangShiuh added 8 commits July 8, 2024 08:40

[phase, kinetics] make Kinetics belong to ThermoPhase

ddfc210

[kinetics] create reactionsAnyMapList

dbb8b39

[Kinetics] add ready()

8c53d14

[Kinetics] add addReaction(kin rxnList)

abcac91

[phase] implement m_collisions a vector of reaction objects for plasm…

27daa73

…a phase and elasticElectronEnergyLossRate

[test] add test for elasticElectronEnergyLossRate

365e887

add flamespeed.yaml to gitignore

35d0cb2

[kinetics] add assertion

6a8156c

BangShiuh changed the title ~~Collision~~ Enhance the PlasmaPhase to use the cross section data with an example of calculating the elastic collision energy loss rate Jul 8, 2024

BangShiuh changed the title ~~Enhance the PlasmaPhase to use the cross section data with an example of calculating the elastic collision energy loss rate~~ Enhance the PlasmaPhase to use the cross section data with an example of calculating the elastic collision energy loss rate (WIP) Jul 8, 2024

BangShiuh changed the title ~~Enhance the PlasmaPhase to use the cross section data with an example of calculating the elastic collision energy loss rate (WIP)~~ Enhance the PlasmaPhase to use the cross section data with an example of calculating the elastic collision energy loss rate Jul 8, 2024

BangShiuh added 2 commits July 14, 2024 23:10

[phase] use second order central difference scheme

af25b4b

update test for using non-uniform central difference

7d866c9

speth requested changes Jul 22, 2024

View reviewed changes

BangShiuh mentioned this pull request Aug 12, 2024

Non-equilibrium plasma including EEDF computation (two-term approximation) Cantera/enhancements#194

Open

ischoegl added the Thermo label Sep 1, 2024

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Enhance the PlasmaPhase to use the cross section data with an example of calculating the elastic collision energy loss rate #1731

Enhance the PlasmaPhase to use the cross section data with an example of calculating the elastic collision energy loss rate #1731

BangShiuh commented Jul 8, 2024 •

edited

Loading

codecov bot commented Jul 8, 2024 •

edited

Loading

speth left a comment

speth Jul 21, 2024

speth Jul 21, 2024

speth Jul 21, 2024

speth Jul 21, 2024

speth Jul 21, 2024

speth Jul 21, 2024

speth Jul 21, 2024

speth Jul 21, 2024

speth Jul 21, 2024

speth Jul 21, 2024

		double elasticElectronEnergyLossRate()
		double normalizedElasticElectronEnergyLossRate()

	self.assertNear(self.phase.elastic_electron_energy_loss_rate, 4.233683742587e-05)
	assert self.phase.elastic_electron_energy_loss_rate == approx(4.233683742587e-05)

		m_kinetics = newKinetics("bulk");
		m_kinetics->addThermo(shared_from_this());

Enhance the PlasmaPhase to use the cross section data with an example of calculating the elastic collision energy loss rate #1731

Are you sure you want to change the base?

Enhance the PlasmaPhase to use the cross section data with an example of calculating the elastic collision energy loss rate #1731

Conversation

BangShiuh commented Jul 8, 2024 • edited Loading

codecov bot commented Jul 8, 2024 • 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

BangShiuh commented Jul 8, 2024 •

edited

Loading

codecov bot commented Jul 8, 2024 •

edited

Loading