Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

interface_current method and accompanying test added #1398

Merged
merged 8 commits into from
Oct 10, 2022

Conversation

c-randall
Copy link
Contributor

@c-randall c-randall commented Sep 27, 2022

Changes proposed in this pull request

  • Here, a new method is added for the interface kinetics class. The method interface_current(phase) is useful when charge transfer reactions occur at an interface. The method returns the net positive charge moving into phase phase via these reactions (Units: A/m^2).
  • To test the new method, a new test_interface_current test has been added. Prior to the pull request, the package has been built and tested locally. All tests were passed successfully.
  • At present, the new method is implemented in the cython/python folders. If there is interest in implementing this functionality via the full C++ code, please reach out to let me know.

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

A short python script is included below to demonstrate how to use this new method. It references the existing lithium_ion_battery.yaml file already available in the Cantera package.

"""
Interface current: net positve current into a specified phase (A/m^2)
"""

import cantera as ct

file = 'lithium_ion_battery.yaml'

anode = ct.Solution(file,"anode")
elect = ct.Solution(file,"electron")
elyte = ct.Solution(file,"electrolyte")
anode_interface = ct.Interface(file,"edge_anode_electrolyte",[anode,elect,elyte])

anode.X = [0.9, 0.1]
elyte.X = [0.4, 0.3, 0.15, 0.15]

anode.electric_potential = 0.
elyte.electric_potential = 3.

int_current = anode_interface.interface_current(elyte)

print(int_current)

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

@codecov
Copy link

codecov bot commented Sep 27, 2022

Codecov Report

Merging #1398 (76bfe1f) into main (332bb9b) will increase coverage by 0.02%.
The diff coverage is 100.00%.

@@            Coverage Diff             @@
##             main    #1398      +/-   ##
==========================================
+ Coverage   71.03%   71.06%   +0.02%     
==========================================
  Files         363      363              
  Lines       51900    51911      +11     
  Branches    17369    17374       +5     
==========================================
+ Hits        36867    36890      +23     
+ Misses      12686    12674      -12     
  Partials     2347     2347              
Impacted Files Coverage Δ
include/cantera/kinetics/InterfaceKinetics.h 85.71% <ø> (ø)
interfaces/cython/cantera/kinetics.pyx 88.77% <100.00%> (+6.30%) ⬆️
src/kinetics/InterfaceKinetics.cpp 74.51% <100.00%> (+0.76%) ⬆️

📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more

@speth
Copy link
Member

speth commented Sep 27, 2022

Thanks for this contribution, @c-randall.

As a procedural note, please don't introduce merge commits as part of your PR. We will rebase this on top of the current main branch when accepting the PR.

@decaluwe
Copy link
Member

@speth, @bryanwweber @ischoegl etal - in addition to his comment on the C++ code base, @c-randall and I were wondering if there is any desire to also implement in the Matlab interface, at this point, or just wait for the toolbox redesign to be complete?

@ischoegl
Copy link
Member

ischoegl commented Sep 28, 2022

@speth, @bryanwweber @ischoegl etal - in addition to his comment on the C++ code base, @c-randall and I were wondering if there is any desire to also implement in the Matlab interface, at this point, or just wait for the toolbox redesign to be complete?

@decaluwe / @c-randall ... If you'd like to add it to MATLAB, I'd suggest to create the functionality in C++ (as suggested at the top), so it can be called from any API (either directly from C++ or from clib). I.e. the API would just add a wrapper for core functionality.

@decaluwe
Copy link
Member

decaluwe commented Sep 28, 2022

Thanks @ischoegl --yes, that's what we would do, but the question is more whether we are holding off on Matlab updates for the time being until the toolbox redesign is complete.

(also, if/when we do put it in C++, would the correct method be something like void InterfaceKinetics::getInterfaceCurrent(double* i_int))? TIA.

@c-randall c-randall marked this pull request as ready for review September 28, 2022 17:21
@ischoegl
Copy link
Member

ischoegl commented Sep 28, 2022

Thanks @ischoegl --yes, that's what we would do, but the question is more whether we are holding off on Matlab updates for the time being until the toolbox redesign is complete.

(also, if/when we do put it in C++, would the correct method be something like void InterfaceKinetics::getInterfaceCurrent(double* i_int))? TIA.

Since there are plans to include this in other API's, I believe it should be implemented in C++ regardless of whether it's added to the current MATLAB toolbox or not. The signature should be

double InterfaceKinetics::interfaceCurrent() const;

where the get prefix can be omitted (while it is not accessing a member variable, it's still a simple calculation). If you add it to clib as well, adding the method to either of the two MATLAB toolboxes (or Fortran, or the .NET interface) will be very simple, i.e. it's ready to go.

Copy link
Member

@decaluwe decaluwe left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi @c-randall, excellent work moving this into the C++. I have left a few comments below, most of them are nitpicky and related to naming conventions and nomenclature, where I think consistency with other places in the code or with the literature at large might make it (slightly) easier for a reader to follow.

The only one which might matter is for the function declaration in the header file. TBQH, I've never really understood what the size_t data type is or what it's for, but have tried to mimic how it gets used, elsewhere in the code. Given that, you might change the data type of your input phase index.

I also don't know if the declaration needs an const; appended to the end. Nearby functions that return variables do have it, but my C++ is not such that I know why or when this is done. @bryanwweber @speth @ischoegl , etc., can one of you weigh in?

Again, very nice job, and thanks for including a test of these capabilities.

@@ -287,6 +287,8 @@ class InterfaceKinetics : public Kinetics
*/
int phaseStability(const size_t iphase) const;

double InterfaceCurrent(int iPhase);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Given the other declarations in this file, I wonder if this might be preferred:

Suggested change
double InterfaceCurrent(int iPhase);
double InterfaceCurrent(const size_t iphase);

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I read about the const qualifier and size_t type and agree with this. I have updated the code here and in other areas to account for this.

@@ -69,6 +69,7 @@ cdef extern from "cantera/kinetics/InterfaceKinetics.h":
cdef cppclass CxxInterfaceKinetics "Cantera::InterfaceKinetics":
void advanceCoverages(double, double, double, double, size_t, size_t) except +translate_exception
void solvePseudoSteadyStateProblem() except +translate_exception
double InterfaceCurrent(int) except +translate_exception
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If the declaration in InterfaceKinetics.h is changed, then change accordingly, here.

an interface. It is defined here as the net positive charge entering the
phase ``phase`` (Units: A/m^2 for a surface, A/m for an edge reaction).
"""
iPhase = self.phase_index(phase)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i_phase would be more python-ic.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Agreed, this change has been made.

@@ -595,4 +595,30 @@ void InterfaceKinetics::setPhaseStability(const size_t iphase, const int isStabl
}
}

double InterfaceKinetics::InterfaceCurrent(int iPhase)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Again, change accordingly if the declaration changes.

@@ -595,4 +595,30 @@ void InterfaceKinetics::setPhaseStability(const size_t iphase, const int isStabl
}
}

double InterfaceKinetics::InterfaceCurrent(int iPhase)
{
int sp = thermo(iPhase).nSpecies();
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is nitpicky, but I'd recommend something like nSp or n_sp. sp in this class typically refers to an actual species object.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Agreed, I changed to nSp.

{
int sp = thermo(iPhase).nSpecies();
double charge_k[sp];
double net_k[sp];
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Again, nitpicky, but something like sdot_k would be more consistent with the nomenclature in the literature, which makes it a little easier for an interested reader to comprehend.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Agreed, this has been changed.

phases = [anode_int, anode, elect, elyte]

for p in phases:
productions = anode_int.get_net_production_rates(p)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Again, a name like sdot_k or even just net_prod_rates would probably be more instructive.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Agreed, variable names are now updated. I went with net_prod_rates here.

Copy link
Member

@ischoegl ischoegl left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi @c-randall ... thanks for the updates!

I have some minor comments that should ensure that the code complies with our C++ coding style (see CONTRIBUTING.md) and avoids unnecessary loops.

include/cantera/kinetics/InterfaceKinetics.h Outdated Show resolved Hide resolved
src/kinetics/InterfaceKinetics.cpp Outdated Show resolved Hide resolved
src/kinetics/InterfaceKinetics.cpp Outdated Show resolved Hide resolved
src/kinetics/InterfaceKinetics.cpp Outdated Show resolved Hide resolved
src/kinetics/InterfaceKinetics.cpp Outdated Show resolved Hide resolved
Copy link
Member

@ischoegl ischoegl left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@c-randall ... thanks for the quick turn-around. I located some more missing white spaces, but otherwise this looks great!

src/kinetics/InterfaceKinetics.cpp Outdated Show resolved Hide resolved
test/python/test_kinetics.py Outdated Show resolved Hide resolved
test/python/test_kinetics.py Outdated Show resolved Hide resolved
Copy link
Member

@ischoegl ischoegl left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

LGTM - thanks!

@ischoegl ischoegl merged commit f1d86d1 into Cantera:main Oct 10, 2022
@decaluwe
Copy link
Member

@ischoegl - many thanks for looking this over, and for the improvements.

One question on one of your comments, for my own knowledge. I nearly commented on this, but since Interface::netProdRates returns a doublereeal, wasn't sure: does the function calling it not need to accept it's return type as such? Thanks again.

@ischoegl
Copy link
Member

Sure. doublereal and double can be used interchangeably as they really are the same type - there are some historical Fortran-related definitions in config.h, see

// define types doublereal, integer, and ftnlen to match the
// corresponding Fortran data types on your system. The defaults
// are OK for most systems

typedef double doublereal;   // Fortran double precision
typedef int integer;      // Fortran integer
typedef int ftnlen;       // Fortran hidden string length type

As far as I am aware, there isn't a reason to keep doublereal around any longer. At the same time, replacing them all - while easy - would create unnecessary churn that should be avoided.

@c-randall c-randall deleted the faradaic_current branch December 7, 2022 17:59
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants