Merge branch 'main' into joss

docs/source/conf.py

@@ -3,7 +3,7 @@
 # -- Project information
 
 project = 'EchemFEM'
-copyright = u'2022-2023, Lawrence Livermore National Laboratory'
+copyright = u'2022-2024, Lawrence Livermore National Laboratory'
 author = 'Thomas Roy'
 
 release = '0.1'
@@ -21,6 +21,11 @@
     'sphinx.ext.mathjax',
 ]
 
+mathjax3_config = {
+  'loader': {'load': ['[tex]/mhchem']},
+  'tex': {'packages': {'[+]': ['mhchem']}},
+}
+
 intersphinx_mapping = {
     'python': ('https://docs.python.org/3/', None),
     'numpy': ('https://docs.scipy.org/doc/numpy/', None),

docs/source/examples.rst

@@ -0,0 +1,37 @@
+Examples
+========
+
+Here is a brief overview of the examples in `echemfem/examples <https://github.com/LLNL/echemfem/tree/main/examples>`_.
+
+* **Planar flow reactor** with electroneutral Nernst-Planck. A Butler-Volmer expression is used for the redox reaction, and there are no homogeneous bulk reactions. A custom mesh is used with refinements close to the electrodes.
+
+    * `2D case for binary electrolyte <https://github.com/LLNL/echemfem/tree/main/examples/bortels_twoion.py>`_ (:math:`\ce{CuSO4}`)
+    * `Nondimensional version of binary case <https://github.com/LLNL/echemfem/tree/main/examples/bortels_twoion_nondim.py>`_ 
+    * `2D case for three-ion electrolyte <https://github.com/LLNL/echemfem/tree/main/examples/bortels_threeion.py>`_ (:math:`\ce{CuSO4}` and :math:`\ce{H2SO4}`)
+    * `Nondimensional version of three ion case <https://github.com/LLNL/echemfem/tree/main/examples/bortels_threeion_nondim.py>`_ 
+    * `3D version of the three-ion case using custom preconditioners for a Discontinuous Galerkin (DG) scheme, and nondimensional <https://github.com/LLNL/echemfem/tree/main/examples/bortels_threeion_extruded_3D_nondim.py>`_ 
+
+* Simple 1D reaction-diffusion system for :math:`\ce{CO2}` electrolysis in :math:`\ce{KHCO3}`
+
+    * `Simplified bicarbonate bulk reactions <https://github.com/LLNL/echemfem/tree/main/examples/gupta.py>`_
+    * `Full bicarbonate bulk reactions <https://github.com/LLNL/echemfem/tree/main/examples/carbonate.py>`_
+    * `Different implementation of the full bicarbonate bulk reactions <https://github.com/LLNL/echemfem/tree/main/examples/carbonate_homog_params.py>`_ (using the :doc:`homog_params <user_guide/homog_params>` interface)
+
+* 2D flow-past the electrode with advection-diffusion
+
+    * `Two species toy model with shear flow <https://github.com/LLNL/echemfem/tree/main/examples/tworxn.py>`_
+    * `CO2 electrolysis in bicarbonate with linear charge-transfer kinetics and shear flow <https://github.com/LLNL/echemfem/tree/main/examples/bicarb.py>`_
+    * `Two species toy model with irregular electrode with Navier-Stokes for the flow <https://github.com/LLNL/echemfem/tree/main/examples/tworxn_irregular.py>`_
+
+* 1D model for the :math:`\ce{CO2}` electrolysis in a copper catalyst layer of a gas diffusion electrode (GDE). The model uses electroneutral Nernst-Planck in a porous medium.
+
+    * `Using substitution of the electroneutrality equation to eliminate a species <https://github.com/LLNL/echemfem/tree/main/examples/catalyst_layer_Cu.py>`_
+    * `Solving the electroneutrality equation explicitly <https://github.com/LLNL/echemfem/tree/main/examples/catalyst_layer_Cu_full.py>`_
+
+* `A simple Vanadium flow battery using advection-diffusion-reaction, Poisson for the ionic potential with a predefinied conductivity and Navier-Stokes-Brinkman for the flow <https://github.com/LLNL/echemfem/tree/main/examples/simple_flow_battery.py>`_
+
+* `A symmetric cylindrical pore model for CO2 electrolysis using electroneutral Nernst-Planck and simplified bicarbonate bulk reactions <https://github.com/LLNL/echemfem/tree/main/examples/cylindrical_pore.py>`_
+
+* `A tandem flow-past the electrode system with an Ag catalyst placed in front of a Cu catalyst with electroneutral Nernst-Planck and Tafel kinetics <https://github.com/LLNL/echemfem/tree/main/examples/bicarb_Ag_Cu_tandem_example.py>`_
+
+* `Simple example for the BMCSL model for finite size ion effects <https://github.com/LLNL/echemfem/tree/main/examples/BMCSL.py>`_

docs/source/index.rst

@@ -24,6 +24,7 @@ Contents
 .. toctree::
    :titlesonly:
 
+   examples
    api
 
 Copyright and License

docs/source/user_guide/homog_params.rst

@@ -7,17 +7,20 @@ each reaction. Below is a list of different keys that should appear in each
 dictionary
 
 * :Key: ``"stoichiometry"``
-  :Type: :py:class:`dict`
+  :Type: :py:class:`dict` of :py:class:`int`
   :Description: This entry defines the stoichiometry of the reaction. Each key in this dictionary should be the name of a species (as defined in the ``"name"`` key values of ``conc_params``). The corresponding value is the stoichemetry coefficient for the species in this reaction. It should be an integer: negative for a reactant, positive for a product.
 * :Key: ``"forward rate constant"``
   :Type: :py:class:`float` or Firedrake expression
   :Description: Rate constant for the forward reaction.
-* :Key: ``"backward rate constant"`` or ``"equilibrium constant"``
+* :Key: ``"backward rate constant"`` or ``"equilibrium constant"`` (optional)
   :Type: :py:class:`float` or Firedrake expression
-  :Description: Rate constant for the backward reaction or equilibrium constant of the reaction. Only of one these is required.
+  :Description: Rate constant for the backward reaction or equilibrium constant of the reaction. For a reversible reaction, only of one these is required. For an irreversible reaction, leave unset.
 * :Key: ``"reference concentration"`` (optional)
   :Type: :py:class:`float` or Firedrake expression
   :Description: Value of 1M in the appropriate units. By setting this, the rate constants are assumed to have the same units (typically 1/s), and the equilibrium constant to be nondimensional.
+* :Key: ``"reaction order"`` (optional)
+  :Type: :py:class:`dict` of :py:class:`int`
+  :Description: This entry can be used to enforce reaction order. Each key in this dictionary should be the name of a species (as defined in the ``"name"`` key values of ``conc_params``). The corresponding value (positive integer) is reaction order of the species in this reaction. If unset, the absolute value of the stoichiometry coeffcient is the reaction order.
 
 Assuming an arbitrary homogeneous reaction system where the forward and
 backward rate constants for reaction :math:`i` are :math:`k_i^f` and
@@ -32,6 +35,9 @@ of mass action, the reaction for species :math:`j` is given by
 where :math:`c_n` is the concentration of species :math:`n`. If the equilibrium
 constant :math:`K_i` is provided instead of the backward rate constant, then it
 is automatically recovered via :math:`k_i^b = \dfrac{k_i^f}{K_i}`. 
+The above formula assumes single-step reactions so that the reaction orders are
+given by the stoichiometry coefficients. If that is not the case, the exponents
+can be replaced by the ``"reaction order"`` inputs.
 
 Note that here the rate constants can have different units, depending on the
 number of reactants or products. It is also common to write the reactions in
@@ -54,9 +60,9 @@ Here is an example of a reaction system that can be implement via
 
 .. math::
 
-   \mathrm{CO}_2 + \mathrm{OH}^- \xrightarrow[k_1^b]{k_1^f} \mathrm{HCO}_3^-
+   \mathrm{CO}_2 + \mathrm{OH}^- \xrightleftharpoons[k_1^b]{k_1^f} \mathrm{HCO}_3^-
 
-   \mathrm{HCO}_3^- + \mathrm{OH}^- \xrightarrow[k_2^b]{k_2^f} \mathrm{CO}_3^{2-}
+   \mathrm{HCO}_3^- + \mathrm{OH}^- \xrightleftharpoons[k_2^b]{k_2^f} \mathrm{CO}_3^{2-} + \mathrm{H}_2\mathrm{O}
 
 Here, we assume the dilute solution case where :math:`\mathrm{H}_2\mathrm{O}`
 concentration is not tracked. In ``conc_params``, we will have entries for the
@@ -79,8 +85,31 @@ list:
       "backward rate constant": k2b
      }]
 
+Now here is another bicarbonate approximation with a high pH approximation,
+where the reaction system is simplified into a single irreversible reaction:
+
+.. math::
+
+   \mathrm{CO}_2 + 2\mathrm{OH}^- \xrightarrow[]{k_1^f} \mathrm{CO}_3^{2-} + \mathrm{H}_2\mathrm{O}
+
+Since this is not actually a single-step reaction, we need to specify the
+reaction order for :math:`\mathrm{OH}^-`. We can use the following
+``homog_params`` list:
+
+.. code::
+
+    [{"stoichiometry": {"CO2": -1,
+                       "OH": -2,
+                       "CO3": 1},
+     "reaction order": {"OH": 1},
+     "forward rate constant": k1f
+     }]
+
+
 As an alternative to the ``homog_params`` interface, the reactions can be
 written directly and passed as a function in ``physical_params["bulk
-reaction"]`` (:doc:`physical_params`). See ``examples/carbonate.py`` and
+reaction"]`` (:doc:`physical_params`). See ``examples/cylindrical_pore.py`` for
+the direct implementation of the high pH approximation of the bicarbonate bulk
+reactions. See ``examples/carbonate.py`` and
 ``examples/carbonate_homog_params.py`` for two equivalent implementations of
-the bicarbonate bulk reactions.
+the full bicarbonate bulk reactions.

echemfem/solver.py

@@ -368,7 +368,7 @@ def _internal_dS(subdomain_id=None, domain=None):
                 us[self.num_liquid:self.num_liquid + self.num_gas], self.pg, gas_params)
 
         # setup applied voltage
-        if self.flow["poisson"] or self.flow["electroneutrality"] or self.flow["electroneutrality full"]:
+        if physical_params.get("U_app"):
             U_app = physical_params["U_app"]
             if isinstance(U_app, (Constant, Function)):
                 self.U_app = U_app
@@ -636,7 +636,8 @@ def setup_solver(self, initial_guess=True, initial_solve=True):
 
             elif (self.flow["electroneutrality"] or self.flow["poisson"] or self.flow["electroneutrality full"]):
                 U0 = Function(self.Vu)
-                U0.assign(self.U_app / 2)
+                if hasattr(self, "U_app"):
+                    U0.assign(self.U_app / 2)
                 if initial_solve:
                     v0 = TestFunction(self.Vu)
                     u0 = [us[i] for i in range(self.num_mass)]
@@ -1223,7 +1224,7 @@ def bulk_reaction(u):
                 reaction_f[idx] = reaction["forward rate constant"]
                 if reaction.get("backward rate constant"):
                     reaction_b[idx] = reaction["backward rate constant"]
-                else:
+                elif reaction.get("equilibrium constant"):
                     reaction_b[idx] = reaction["forward rate constant"] / \
                                         reaction["equilibrium constant"]
                 if reaction.get("reference concentration"):
@@ -1232,11 +1233,15 @@ def bulk_reaction(u):
                     c_ref = 1.0
                 for name in reaction["stoichiometry"]:
                     s = reaction["stoichiometry"][name]
+                    order = abs(s) # default
+                    if reaction.get("reaction order"):
+                        if reaction["reaction order"].get(name):
+                            order = reaction["reaction order"][name]
                     a = u[self.i_c[name]] / c_ref
                     if s < 0:
-                        reaction_f[idx] *= a ** (-s)
+                        reaction_f[idx] *= a ** order
                     if s > 0:
-                        reaction_b[idx] *= a ** s
+                        reaction_b[idx] *= a ** order
 
             for i in self.idx_c:
                 name = self.conc_params[i]["name"]

examples/BMCSL.py

@@ -3,7 +3,12 @@
 
 @author: Nitish Govindarajan
 
+Simple example for the BMCSL model for finite size ion effects
+
 Code to reproduce Figure 1 from doi:10.1016/j.jcis.2007.08.006
+Biesheuvel, P.M. and Van Soestbergen, M., 2007. Counterion volume effects in
+mixed electrical double layers. Journal of Colloid and Interface Science,
+316(2), pp.490-499.
 """
 
 # Import required libraries

examples/bicarb.py

@@ -5,7 +5,8 @@
 class CarbonateSolver(EchemSolver):
     def __init__(self):
         """
-        Bicarbonate example reproduced from:
+        Flow past an electrode with bicarbonate bulk reactions and linear CO2 electrolysis
+        Example reproduced from:
         Lin, T.Y., Baker, S.E., Duoss, E.B. and Beck, V.A., 2021. Analysis of
         the Reactive CO2 Surface Flux in Electrocatalytic Aqueous Flow
         Reactors. Industrial & Engineering Chemistry Research, 60(31),

examples/bicarb_Ag_Cu_tandem_example.py

@@ -11,19 +11,19 @@
 comm = MPI.COMM_WORLD
 rank = comm.Get_rank()
 
-##Description:
-#This example case runs a tandem electrode system, with an Ag catalyst placed
-#in front of a Cu catalyst, at a Peclet number of 1E3. 
-#The advection/diffusion/migration/electroneutrality equations are solved 
-#for the bulk electrolyte, and a Tafel description is applied for the surface 
-#catalyst reactions.
+"""
+This example case runs a tandem electrode system, with an Ag catalyst placed
+in front of a Cu catalyst, at a Peclet number of 1E3. 
+The advection/diffusion/migration/electroneutrality equations are solved 
+for the bulk electrolyte, and a Tafel description is applied for the surface 
+catalyst reactions.
 
-#Assumed products of Cu catalyst are C2H4, C2H6O, and H2.
-#Assumed products of Ag catalysis are CO and H2.
+Assumed products of Cu catalyst are C2H4, C2H6O, and H2.
+Assumed products of Ag catalysis are CO and H2.
+
+--------------------------------------------------------------------------
+list of references
 
-##--------------------------------------------------------------------------
-#list of references
-"""
 Bicarbonate flow reactor setup given in:
 Lin, T.Y., Baker, S.E., Duoss, E.B. and Beck, V.A., 2021. Analysis of
 the Reactive CO2 Surface Flux in Electrocatalytic Aqueous Flow
@@ -53,10 +53,9 @@
 Microkinetics with Continuum Transport Models to Understand Electrochemical 
 CO2 Reduction in Flow Reactors. PRX Energy, 2(3), pp.033010.
 
+--------------------------------------------------------------------------
+end of list of references
 """
-##--------------------------------------------------------------------------
-#end of list of references
-
 
 
 #---------------------------------------------------------------------------

examples/bortels_threeion.py

@@ -7,6 +7,20 @@
 parser.add_argument("--family", type=str, default='DG')
 args, _ = parser.parse_known_args()
 
+"""
+2D Flow-plate reactor with electroneutral Nernst-Planck. Using a custom gmsh
+mesh with refinement close to the electrodes. GMSH is used to get the mesh file
+from the .geo file as follows:
+
+    gmsh -2 bortels_structuredquad.geo
+
+Three ion case taken from: 
+Bortels, L., Deconinck, J. and Van Den Bossche, B., 1996. The multi-dimensional
+upwinding method as a new simulation tool for the analysis of multi-ion
+electrolytes controlled by diffusion, convection and migration. Part 1. Steady
+state analysis of a parallel plane flow channel. Journal of Electroanalytical
+Chemistry, 404(1), pp.15-26.
+"""
 
 class BortelsSolver(EchemSolver):
     def __init__(self):

examples/bortels_threeion_extruded_3D_nondim.py

@@ -8,6 +8,24 @@
 
 set_log_level(DEBUG)
 
+"""
+3D Flow-plate reactor with electroneutral Nernst-Planck. Using a custom gmsh
+mesh with refinement close to the electrodes. The mesh is then extruded in the
+z-direction. Custom block preconditioners are used. The file submitter.pl can
+be adapted to launch slurm jobs using perl.
+
+Three ion case adapted from the 2D case in: 
+Bortels, L., Deconinck, J. and Van Den Bossche, B., 1996. The multi-dimensional
+upwinding method as a new simulation tool for the analysis of multi-ion
+electrolytes controlled by diffusion, convection and migration. Part 1. Steady
+state analysis of a parallel plane flow channel. Journal of Electroanalytical
+Chemistry, 404(1), pp.15-26.
+
+Nondimensionalization and preconditioning from:
+Roy, T., Andrej, J. and Beck, V.A., 2023. A scalable DG solver for the
+electroneutral Nernst-Planck equations. Journal of Computational Physics, 475,
+p.111859.
+"""
 
 parser = argparse.ArgumentParser()
 parser.add_argument('-ref_levels', type=int, default=0)
@@ -55,6 +73,7 @@ def coarsen_form(self, form, args):
 class BortelsSolver(EchemSolver):
     def __init__(self):
         # Create an initial coarse mesh
+        # Note: will need to use "gmsh -2 mesh_name.geo" command to convert to .msh format
         if args.mesh == 'structured':
             plane_mesh = Mesh('bortels_structuredquad_nondim_coarse1664.msh')
             layers = 8
@@ -197,7 +216,6 @@ def set_velocity(self):
 
 # Custom solver parameters
 
-
 class CoarsenPenaltyPMGPC(P1PC):
     def coarsen_form(self, form, replace_dict):
         k = 1

examples/bortels_threeion_nondim.py

@@ -1,6 +1,25 @@
 from firedrake import *
 from echemflow import EchemSolver
 
+"""
+2D Flow-plate reactor with electroneutral Nernst-Planck. Using a custom gmsh
+mesh with refinement close to the electrodes. GMSH is used to get the mesh file
+from the .geo file as follows:
+
+    gmsh -2 bortels_structuredquad_nondim.geo
+
+Three ion case taken from: 
+Bortels, L., Deconinck, J. and Van Den Bossche, B., 1996. The multi-dimensional
+upwinding method as a new simulation tool for the analysis of multi-ion
+electrolytes controlled by diffusion, convection and migration. Part 1. Steady
+state analysis of a parallel plane flow channel. Journal of Electroanalytical
+Chemistry, 404(1), pp.15-26.
+
+Nondimensionalization from:
+Roy, T., Andrej, J. and Beck, V.A., 2023. A scalable DG solver for the
+electroneutral Nernst-Planck equations. Journal of Computational Physics, 475,
+p.111859.
+"""
 
 class BortelsSolver(EchemSolver):
     def __init__(self):

examples/bortels_twoion.py

@@ -5,6 +5,20 @@
 parser.add_argument("--family", type=str, default='DG')
 args, _ = parser.parse_known_args()
 
+"""
+2D Flow-plate reactor with electroneutral Nernst-Planck. Using a custom GMSH
+mesh with refinement close to the electrodes. GMSH is used to get the mesh file
+from the .geo file as follows:
+
+    gmsh -2 bortels_structuredquad.geo
+
+Two ion case taken from: 
+Bortels, L., Deconinck, J. and Van Den Bossche, B., 1996. The multi-dimensional
+upwinding method as a new simulation tool for the analysis of multi-ion
+electrolytes controlled by diffusion, convection and migration. Part 1. Steady
+state analysis of a parallel plane flow channel. Journal of Electroanalytical
+Chemistry, 404(1), pp.15-26.
+"""
 
 class BortelsSolver(EchemSolver):
     def __init__(self):

examples/bortels_twoion_nondim.py

@@ -1,6 +1,25 @@
 from firedrake import *
 from echemfem import EchemSolver
 
+"""
+2D Flow-plate reactor with electroneutral Nernst-Planck. Using a custom gmsh
+mesh with refinement close to the electrodes. GMSH is used to get the mesh file
+from the .geo file as follows:
+
+    gmsh -2 bortels_structuredquad_nondim.geo
+
+Two ion case taken from: 
+Bortels, L., Deconinck, J. and Van Den Bossche, B., 1996. The multi-dimensional
+upwinding method as a new simulation tool for the analysis of multi-ion
+electrolytes controlled by diffusion, convection and migration. Part 1. Steady
+state analysis of a parallel plane flow channel. Journal of Electroanalytical
+Chemistry, 404(1), pp.15-26.
+
+Nondimensionalization from:
+Roy, T., Andrej, J. and Beck, V.A., 2023. A scalable DG solver for the
+electroneutral Nernst-Planck equations. Journal of Computational Physics, 475,
+p.111859.
+"""
 
 class BortelsSolver(EchemSolver):
     def __init__(self):