update in runProtonicMembrane doc

Documentation/publishedExamples/runProtonicMembrane.rst

@@ -11,7 +11,7 @@ Protonic Membrane model
 
 Load and parse input from given json files
 ==========================================
-The source of the json files can be seen in :battmofile:`protonicMembrane<ProtonicMembrane/jsonfiles/protonicMembrane.json>` and :battmofile:`1d-PM-geometry.json<ProtonicMembrane/jsonfiles/1d-PM-geometry.json>`
+The source of the json files can be seen in :battmofile:`protonicMembrane.json<ProtonicMembrane/jsonfiles/protonicMembrane.json>` and :battmofile:`1d-PM-geometry.json<ProtonicMembrane/jsonfiles/1d-PM-geometry.json>`
 
 .. code-block:: matlab
 
@@ -31,8 +31,11 @@ We setup the input parameter structure which will we be used to instantiate the
 .. code-block:: matlab
 
   inputparams = ProtonicMembraneCellInputParams(jsonstruct);
-  
-  % We setup the grid, which is done by calling the function :battmo:`setupProtonicMembraneCellGrid`
+
+We setup the grid, which is done by calling the function :battmo:`setupProtonicMembraneCellGrid`
+
+.. code-block:: matlab
+
   [inputparams, gen] = setupProtonicMembraneCellGrid(inputparams, jsonstruct);
 
 
@@ -62,19 +65,13 @@ We setup the initial state using a default setup included in the model
 
 Schedule
 ========
-We setup the schedule, which means the timesteps and also the control we want to use. In this case we use current control and the current equal to zero (see here :battmofile:`here<ProtonicMembrane/protonicMembrane.json#86>`).
-We compute the steady-state solution so that the time stepping here is more an artifact to reach the steady-state solution. In particular, it governs the pace at which we increase the non-linearity (not detailed here).
+We setup the schedule, which means the timesteps and also the control we want to use. In this case we use current control and the current equal to zero (see :battmofile:`here<ProtonicMembrane/jsonfiles/protonicMembrane.json#118>`).
+We compute the steady-state solution and the time stepping here does not correspond to time values but should be seen as step-wise increase of the effect of the non-linearity (in particular in the expression of the conductivity which includes highly nonlineaer effect with the exponential terms. We do not detail here the method).
 
 .. code-block:: matlab
 
   schedule = model.Control.setupSchedule(inputparams.jsonstruct);
 
-We change the default tolerance
-
-.. code-block:: matlab
-
-  model.nonlinearTolerance = 1e-8;
-
 
 Simulation
 ==========
@@ -87,7 +84,7 @@ We run the simulation
 
 Plotting
 ========
-We setup som shortcuts for convenience
+We setup som shortcuts for convenience and introduce plotting options
 
 .. code-block:: matlab
 
@@ -105,7 +102,7 @@ We recover the position of the mesh cell of the discretization grid. This is use
 
   xc = model.(elyte).grid.cells.centroids(:, 1);
 
-We consider the solution obtained at the last time step, which corresponds to the solution at steady-state.
+We consider the solution obtained at the last time step, which corresponds to the solution at steady-state. The second line adds to the state variable all the variables that are derived from our primary unknowns.
 
 .. code-block:: matlab
 

Documentation/publishedExamples/runProtonicMembranePreamble.rst

@@ -6,21 +6,23 @@ We consider the model of a mixed proton and electron conducting membrane, as des
 Governing equations
 -------------------
 
-We have three components given by the proton (:math:`H^+`), the :math:`p` and :math:`n` type. We want to find expressions for the fluxes for each of the components
+In the membrane, we have three components or *species* given by the proton (:math:`H^+`), and the :math:`p` and :math:`n`
+type charge carriers. We need expressions for the fluxes for each of those. The governing equations will then be given
+by charge and mass conservation equations.
 
 We denote by :math:`\phi` the electrostatic potential. For each of the components :math:`\alpha=\{H^+, p, n\}`, we
 introduce the electrochemical potential denoted :math:`\bar\mu_\alpha` and the chemical potential denoted
 :math:`\mu_\alpha`. The potentials are related with each other through the relation
 
 .. math::
 
-   \bar\mu_\alpha = \mu_\alpha + z_\alpha F \phi
+   \bar\mu_\alpha = \mu_\alpha + z_\alpha F \phi.
 
 The fluxes are governed by the gradient of the electrochemical potential. We have
 
 .. math::
    
-   j_{\alpha} = -k_\alpha\nabla\bar\mu_\alpha
+   j_{\alpha} = -k_\alpha\nabla\bar\mu_\alpha,
 
 for some coefficient :math:`k_\alpha` which is not necessarily a constant.
 
@@ -43,18 +45,17 @@ For the :math:`p` and :math:`n` type conductivities,  we use the empirical relat
 
    \sigma_p(E) = \sigma_p^0\exp\left(\frac{F(E - E_{\text{ref},p})}{RT}\right)\quad\text{ and }\quad\sigma_n(E) = \sigma_n^0\exp\left(\frac{-F(E - E_{\text{ref},n})}{RT}\right).
 
-Here :math:`E_{\text{ref},p}` and :math:`E_{\text{ref},n}` are two reference potentials (see below).   
+Here :math:`E_{\text{ref},p}` and :math:`E_{\text{ref},n}` are two reference potentials. 
 
 We consider the steady state. We could introduce later charge and mass capacitors. The unknowns are the functions :math:`\phi(x)` and :math:`E(x)` in the electrolyte. The governing equations are given by the mass conservation for the proton and the charge conservation.
 
-The mass conservation for $H^+$ is given by
+The **mass conservation equation** for :math:`H^+` is given by
 
 .. math::
 
    \nabla\cdot j_{H^+} = 0.
 
-The current density is given by :math:`i = F (j_{H^+} + j_p - j_n)` and the \textbf{charge
-conservation equation} is
+The total current density is given by :math:`i = F (j_{H^+} + j_p - j_n)` and the **charge conservation equation** is
 
 .. math::
 
@@ -72,7 +73,8 @@ We can rewrite :math:`i_{\text{el}}` as
    
    i_{\text{el}} = - (\sigma_p(E) + \sigma_n(E))\nabla ( E + \phi ).
 
-The governing equations for :math:`\phi(x)` and :math:`\pi(x)` are therefore the differential equations
+We finally obtain the governing equations for :math:`\phi(x)` and :math:`\pi(x)` as the following system of differential
+equations
 
 .. math::
 
@@ -91,17 +93,18 @@ We define the over-potential :math:`\eta` as
    \eta_\text{elde} = \pi_\text{elde} - \phi_\text{elde} - \text{OCP}_\text{elde}
 
 where :math:`\text{OCP}_\text{elde}` is the open-circuit potential for the given electrode. The value of the
-$\text{OCP}$ at each electrode depend on the composition at the electrode (see :cite:`V_llestad_2019` for the expressions).
+OCP at each electrode depend on the composition at the electrode (see :cite:`V_llestad_2019` for the expressions).
 
 At the anode, we imposte that the proton current  is given through the following Buttler-Volmer type expression
 
 .. math::
    i_{H^+, \text{an}} = -i_0\frac{e^{-\beta\frac{ z F \eta_{\text{an}}}{RT}} - e^{( 1- \beta)\frac{ z F \eta_{\text{an}}}{RT}}}{ 1+ \frac{i_0}{i_{l,c}}e^{-\beta\frac{ z F \eta_{\text{an}}}{RT}} -  \frac{i_0}{i_{l,a}}e^{( 1- \beta)\frac{ z F \eta_{\text{an}}}{RT}}}
 
 
-Here, :math:`i_{l,c}` and :math:`i_{l,a}` are given constants. The value of the reference current density $i_0$ is also constant 
+Here, :math:`i_{l,c}` and :math:`i_{l,a}` are given constants. The value of the reference current density :math:`i_0` is
+also constant.
 
-The total current is given by :math:`i_{\text{an}} = I` for some constant current $I$.
+The total current is given by :math:`i_{\text{an}} = I` for some constant current :math:`I`.
 
 At the cathode, we impose that the electrostatic potential is equal to zero and 
 a relation between the :math:`H^+` current and the over-potential that takes a linear form,

Documentation/publishedExamples/runProtonicMembrane_source.rst

@@ -11,7 +11,7 @@ Source code for runProtonicMembrane
   %% Protonic Membrane model
   
   %% Load and parse input from given json files
-  % The source of the json files can be seen in :battmofile:`protonicMembrane<ProtonicMembrane/jsonfiles/protonicMembrane.json>` and
+  % The source of the json files can be seen in :battmofile:`protonicMembrane.json<ProtonicMembrane/jsonfiles/protonicMembrane.json>` and
   % :battmofile:`1d-PM-geometry.json<ProtonicMembrane/jsonfiles/1d-PM-geometry.json>`
   
   filename = fullfile(battmoDir(), 'ProtonicMembrane', 'jsonfiles', 'protonicMembrane.json');
@@ -27,6 +27,7 @@ Source code for runProtonicMembrane
   
   inputparams = ProtonicMembraneCellInputParams(jsonstruct);
   
+  %%
   % We setup the grid, which is done by calling the function :battmo:`setupProtonicMembraneCellGrid`
   [inputparams, gen] = setupProtonicMembraneCellGrid(inputparams, jsonstruct);
   
@@ -44,17 +45,14 @@ Source code for runProtonicMembrane
   
   %% Schedule
   % We setup the schedule, which means the timesteps and also the control we want to use. In this case we use current
-  % control and the current equal to zero (see here :battmofile:`here<ProtonicMembrane/protonicMembrane.json#86>`).
+  % control and the current equal to zero (see :battmofile:`here<ProtonicMembrane/jsonfiles/protonicMembrane.json#118>`).
   %
-  % We compute the steady-state solution so that the time stepping here is more an artifact to reach the steady-state
-  % solution. In particular, it governs the pace at which we increase the non-linearity (not detailed here).
+  % We compute the steady-state solution and the time stepping here does not correspond to time values but should be seen
+  % as step-wise increase of the effect of the non-linearity (in particular in the expression of the conductivity which
+  % includes highly nonlineaer effect with the exponential terms. We do not detail here the method).
   
   schedule = model.Control.setupSchedule(inputparams.jsonstruct);
   
-  %%
-  % We change the default tolerance
-  model.nonlinearTolerance = 1e-8;
-  
   %% Simulation
   % We run the simulation
   
@@ -64,7 +62,7 @@ Source code for runProtonicMembrane
   %
   
   %%
-  % We setup som shortcuts for convenience
+  % We setup som shortcuts for convenience and introduce plotting options
   an    = 'Anode';
   ct    = 'Cathode';
   elyte = 'Electrolyte';
@@ -79,7 +77,8 @@ Source code for runProtonicMembrane
   xc = model.(elyte).grid.cells.centroids(:, 1);
   
   %%
-  % We consider the solution obtained at the last time step, which corresponds to the solution at steady-state.
+  % We consider the solution obtained at the last time step, which corresponds to the solution at steady-state. The second
+  % line adds to the state variable all the variables that are derived from our primary unknowns.
   state = states{end};
   state = model.addVariables(state, schedule.control);
   

ProtonicMembrane/examples/runProtonicMembrane.m

@@ -1,7 +1,7 @@
 %% Protonic Membrane model
 
 %% Load and parse input from given json files
-% The source of the json files can be seen in :battmofile:`protonicMembrane<ProtonicMembrane/jsonfiles/protonicMembrane.json>` and
+% The source of the json files can be seen in :battmofile:`protonicMembrane.json<ProtonicMembrane/jsonfiles/protonicMembrane.json>` and
 % :battmofile:`1d-PM-geometry.json<ProtonicMembrane/jsonfiles/1d-PM-geometry.json>`
 
 filename = fullfile(battmoDir(), 'ProtonicMembrane', 'jsonfiles', 'protonicMembrane.json');
@@ -35,17 +35,14 @@
 
 %% Schedule
 % We setup the schedule, which means the timesteps and also the control we want to use. In this case we use current
-% control and the current equal to zero (see here :battmofile:`here<ProtonicMembrane/protonicMembrane.json#86>`).
+% control and the current equal to zero (see :battmofile:`here<ProtonicMembrane/jsonfiles/protonicMembrane.json#118>`).
 %
-% We compute the steady-state solution so that the time stepping here is more an artifact to reach the steady-state
-% solution. In particular, it governs the pace at which we increase the non-linearity (not detailed here).
+% We compute the steady-state solution and the time stepping here does not correspond to time values but should be seen
+% as step-wise increase of the effect of the non-linearity (in particular in the expression of the conductivity which
+% includes highly nonlineaer effect with the exponential terms. We do not detail here the method).
 
 schedule = model.Control.setupSchedule(inputparams.jsonstruct);
 
-%%
-% We change the default tolerance
-model.nonlinearTolerance = 1e-8;
-
 %% Simulation
 % We run the simulation
 
@@ -55,7 +52,7 @@
 %
 
 %%
-% We setup som shortcuts for convenience
+% We setup som shortcuts for convenience and introduce plotting options
 an    = 'Anode';
 ct    = 'Cathode';
 elyte = 'Electrolyte';
@@ -70,7 +67,8 @@
 xc = model.(elyte).grid.cells.centroids(:, 1);
 
 %%
-% We consider the solution obtained at the last time step, which corresponds to the solution at steady-state.
+% We consider the solution obtained at the last time step, which corresponds to the solution at steady-state. The second
+% line adds to the state variable all the variables that are derived from our primary unknowns.
 state = states{end};
 state = model.addVariables(state, schedule.control);