# A simple Python implementation of the 1D C++ example code provded in | |
# Understanding Voltammetry: Simulation of Electrode Processes (2nd Edition), World Scientific, 2020 | |
# By Richard Guy Compton, Enno Kätelhön, Kristopher R Ward, Eduardo Laborda | |
import math | |
theta_i = 20.0 | |
theta_v = -20.0 | |
sigma = 100.0 | |
deltaX = 2E-4 | |
deltaTheta = 0.02 | |
deltaT = deltaTheta / sigma | |
maxT = 2*abs(theta_v - theta_i) / sigma | |
maxX = 6*maxT**0.5 | |
n = int(maxX/deltaX) | |
m = int(maxT/deltaT) | |
lamb = deltaT /(deltaX**2) | |
alpha = -lamb | |
beta = 2.0*lamb + 1.0 | |
gamma = -lamb | |
g_mod = [0.0]*n | |
C = [1.0]*n | |
for i in range(1,n-1): | |
g_mod[i] = gamma / (beta - g_mod[i-1] *alpha) | |
Theta = theta_i | |
flux = [] | |
pot = [] | |
for k in range(m): | |
if (k<m/2): | |
Theta -= deltaTheta | |
else: | |
Theta += deltaTheta | |
C[0] = 1.0/ (1.0+math.exp(-Theta)) | |
for i in range(1,n-1): | |
C[i] = (C[i] - C[i-1]*alpha) / (beta - g_mod[i-1]*alpha) | |
C[n-1] = 1.0; | |
for i in range(n-2,-1,-1): | |
C[i] = C[i] - g_mod[i]*C[i+1] | |
flux.append(-(-C[2]+4*C[1]-3*C[0])/ (2*deltaX)) | |
pot.append(Theta) | |
with open('datapy.txt', 'w') as fp: | |
for item in zip(pot,flux): | |
fp.write("%s\n" %str(item)[1:-1]) |
Finite Difference Simulation: Voltammetry in Python
As a programming language Python is many things... but one thing it is not is fast. Does this mean we can't use it in electrochemical simulation? The book `Understanding Voltammetry: Simulation of Electrode Process' presents example programs written in C++ that enable the simulation of electrochemical processes using finite difference methods. The question is, can we translate these example programs into Python and allow the problems to be solved in a reasonable amount of time?
The article is going to (albeit crudely) focus on benchmarking the execution time for solving the one-dimensional diffusion equation through the use of a finite difference method. Here we take the 1D diffusion equation, discretise it using a central finite difference method and solve the problem fully implicitly, see this book for details on the used method. There are certainly examples of people using Python to simulate cyclic voltammetry using explicit finite difference methods. Explicit procedures are less computationally expensive but are only conditionally stable, this stability issue can be problematic when simulating voltammetric responses of an electroactive species with fast coupled homogeneous kinetics. As such in the following we briefly consider if we can perform fully implicit 1D voltammetric simulations in Python, a further post will consider two 2D problems.
The simulation solves the Diffusion Equation:
subject to the appropriate boundary conditions. After discretisation of this expression, then for a one-dimensional problem, we are left with a set of equations that can be represented by a tridiagonal matrix. In the textbook example C++ code this tridiagonal matrix is solved for directly using an implementation of the Thomas Algorithm. This implementation has not been optimised but even on an antiquated laptop (Intel Core i5-5200U, 2.2 GHz) the simulation only takes 2.7 seconds.
What happens if we take a naïve approach and directly translate this code to Python?
This isn't good... the program takes ~75 s to run; almost 30 times slower than the C++ implementation.
Happily there's a lot of room for improvement. We essentially have two options, first we could use a Python library to do the heavy lifting for us
or second we could write our own pre-compiled extension.
As we have a tridiagonal matrix lets look at what the Scipy library has to offer. Scipy provides access to the Lapack banded matrix solver. Using this solver requires use to modify how we set up the matrix but this isn't a major undertaking. If we use this Fortran based package to solve the matrix we can solve the one-dimensional mass-transport problem in only 3.5 s, a big improvement and competitive with the non optimised C++ code.
import numpy as np | |
import scipy.linalg as la | |
theta_i = 20.0 | |
theta_v = -20.0 | |
sigma = 100.0 | |
deltaX = 2E-4 | |
deltaTheta = 0.02 | |
deltaT = deltaTheta / sigma | |
maxT = 2*abs(theta_v - theta_i) / sigma | |
maxX = 6*np.sqrt(maxT) | |
n = int(maxX/deltaX) | |
m = int(maxT/deltaT) | |
lamb = deltaT /(deltaX**2) | |
ab = np.ones((3,n)) | |
ab[2,:] = np.ones(n)*-lamb | |
ab[1,:] = np.ones(n)+2.0*lamb | |
ab[0,:] = np.ones(n)*-lamb | |
ab[0,1] = 0 | |
ab[1,0] = 1 | |
ab[1,n-1] = 1 | |
ab[2,n-2] = 0 | |
C = np.ones(n) | |
Theta = theta_i | |
flux = [] | |
pot = [] | |
for k in range(m): | |
if (k<m/2): | |
Theta -= deltaTheta | |
else: | |
Theta += deltaTheta | |
C[0] = 1.0/ (1.0+np.exp(-Theta)) | |
C = la.solve_banded((1,1),ab,C) | |
flux.append(-(-C[2]+4*C[1]-3*C[0])/ (2*deltaX)) | |
pot.append(Theta) | |
with open('datapy.txt', 'w') as fp: | |
for item in zip(pot,flux): | |
fp.write("%s\n" %str(item)[1:-1]) |
Given how simple this implementation is... it is pretty good. Even better the Lapack routine is a more general solver allowing us to easily tackle pentadiagonal or even denser matrices.
In many cases just using a Python library can sufficiently decrease a programs execution time enough to justify the use of an easier to understand programming language.
We could finish here but lets see if we can get any further with Cython. Cython gives us an easy way by which we can write a compiled extension for Python but in doing so we loose some of the point of using Python itself. We now have to compile our code and we loose some of clarity of just using Python. Still there are cases in which using Cython seems justifiable; if I have a larger Python program but the runtime is limited by a relative small section of the code and there is no obvious Scipy library that can do the task I want, I will often resort to making a bespoke Python extension. Cython requires us to write a separate .pyx file and to pre-compile it before use. For the present cases if we just transfer a lot of the script directly to a Cython file, we will immediately get some improvements in the runtime. However, the main point is to ensure that we define the types for all of our used variables so that Cython can optimise the compiled code. The script below shows an example Cython function where the some effort has been made to optimise the function providing types for all of the variables. Using this compiled extension the simulation only take 1.7 s, this is faster than the none optimised C++ example given in the textbook!
def sim(double[:] C, double alpha, double beta, double[:] g_mod, double Theta, double deltaTheta, int m, int n, double deltaX): | |
flux = [0.0]*m | |
pot = [0.0]*m | |
cdef int i | |
for k in range(m): | |
if (k<m/2): | |
Theta -= deltaTheta | |
else: | |
Theta+= deltaTheta | |
C[0] = (1.0 / (1.0+math.exp(-Theta))) | |
for i in range(1,n-1): | |
C[i] = (C[i] - C[i-1]*alpha) / (beta - g_mod[i-1]*alpha) | |
C[n-1] = 1.0; | |
for i in range(n-2,-1,-1): | |
C[i] = C[i] - g_mod[i]*C[i+1] | |
flux[k] = (-(-C[2]+4*C[1]-3*C[0])/ (2*deltaX)) | |
pot[k] = (Theta) | |
return pot, flux |
Given the simplicity of the example, the use of Cython does not really seem justified here.
As a final point I just want to highlight that another way in which we can improve
the runtime is to implement an expanding spatial grid so as to reduce the size of the matrix that needs solving.
As a final example I'm providing a fuller example of a 1D voltammetric implicit simulation program
which takes dimensional parameters, converts them to a dimensionless form, uses an expanding spatial grid and employs
the Butler-Volmer equation to describe the interfacial
electron transfer kinetics.
Finally the example script also does some basic plotting and compares the simulated peak height to that expected analytically from the Randles-Ševčík equation.
The runtime for this simulation is about 0.5 s.
############################## | |
#### | |
#### 1D Voltammetry Simulation: Implicit | |
#### BV kinetics= 1 electron oxidation | |
#### | |
############################# | |
# See 'UV: Simultion of Electrode Processes' 2nd Ed. pg 60- for similar simulation C++ | |
# main differences- this implements BV kintics (not Nernst) and has an non-uniform expanding grid | |
#import some useful functions | |
import math | |
import numpy as np # provides a better way of handling numbers and arrays in python | |
import scipy.linalg as la # sparse matrix solver | |
import matplotlib.pyplot as plt # a basic plotting library | |
from datetime import datetime # function to allow the simulation time to be reported | |
plt.close() # close any open plots | |
start = datetime.now() # store the current time | |
####################### | |
#### Physical Constants | |
####################### | |
F = 96485 # Faraday Constant (C/mol) | |
R = 8.314 # Gas Constant (J K-1 mol-1) | |
T = 298 # Temperature | |
############################ | |
#### Experimental parameters | |
############################ | |
Ef = 0.0 # Formal Potential (V) | |
Estart = -0.5 # Start Potential for the Voltammogram (V) | |
Evertex = 0.5 # Potential to which the Electrode is scanned (V) | |
scanrate = 0.1 # Scan rate (V s-1) | |
stepsize = 0.2e-3 # Potential step (V) | |
D = 1e-9 # Diffusion Coefficient (m2 s-1) | |
radius = 1e-3 # Electrode radius (m) | |
k0 = 0.01 # Electron transfer rate (m s-1) | |
conc = 1 # Concentration of analyte (mol m-3) | |
alpha = 0.5 # Transfer Coefficient (dimensionless) | |
##### | |
## Randles-Sevich Equaion | |
#### | |
iprs = 2.69e5*D**0.5*conc*np.pi*radius**2*scanrate**0.5 # Reversible form | |
ipirrev = 2.99e5*D**0.5*conc*np.pi*radius**2*alpha**0.5*scanrate**0.5 # Irreversible Form | |
#### | |
## Dimensionless variables | |
#### | |
# Define two helper functions to Convert the potential (E) to the Dimensionless Potential (Theta) | |
def etotheta(E, Ef): | |
return F*(E-Ef)/(R*T) | |
# does the reverse- converting the dimensionless potential back | |
def thetatoe(theta, Ef): # the conversion method is slightly different depending on whether a python list or a numpy array is passed to the function | |
if type(theta) == list: # first check if a list has been sent to the function | |
temp = [] # if True make an empty list | |
for i in range(len(theta)): # loop through the list and calculate the potential | |
temp.append( (theta[i]*(R*T)/F)+Ef ) | |
return temp | |
else: | |
return (theta*(R*T)/F)+Ef # nomenclature is simpler if a numpy array is passed | |
theta_i = etotheta(Estart, Ef) # Convert the starting potential to theta | |
theta_v = etotheta(Evertex, Ef) | |
sigma = (radius**2/D)*(F/(R*T))*scanrate # Convert the scan rate to the dimensionless form (sigma) | |
deltaTheta = etotheta(stepsize,0.0) # Dimensionless step size | |
K0 = (k0*radius)/D # Dimensionless electron transfer rate | |
deltaT = deltaTheta/sigma # Dimensionless time step | |
maxT = 2*abs(theta_v-theta_i)/sigma # Duration of the experiment in dimensionless time | |
maxX = 6*math.sqrt(maxT) # Maximum distance away from the electrode | |
#Expanding grid | |
omega = 1.02 # Expanding grid factor | |
h0 = 1e-4 # Smallest grid step | |
#make X grid | |
h = h0 | |
x = 0.0 # set the variable x to have the value 0.0 | |
X = [0.0] # makes a list where the first value is 0.0 | |
while x <= maxX: # use a while loop to calculate all of the X positions, the loop stops when the value of x is larger than maxX | |
x += h # shorthand for, x = x + h, every loop add h to the present value of x | |
X.append(x) # every loop put the new value of x at the end of the list | |
h *= omega # every loop multiply x by omega, shorthand for, h = h*omega | |
n = len(X) # find the length (ie. the number of items) in the list X | |
m = int(maxT/deltaT) # find the number of time steps and make sure you have an integer (int) | |
## | |
# Setup concentration and tridiagonal arrays | |
## | |
C = np.ones(n) # make an empty (filled with zeros) numpy array that is n long | |
ab = np.zeros((3,n)) # make a 2D empty numpy array that is comprised of 3 arrays of n length | |
## ab is a numpy array that will be used to hold the finite difference coefficients (alpha, beta, gamma) | |
## we only need to save the diagonal values for the matrix so only need three arrays | |
#fill 2D array with the finite difference coefficients | |
for i in range(1,n-1): # use of for loop to go through the different grid positions but we don't do the first or last values | |
delX_m = X[i] - X[i-1] # find the difference between the current grid position i and the one before (i-1) | |
delX_p = X[i+1] - X[i] # find the difference between the current grid position i and the next one (i+1) | |
ab[0,i+1] = -(2*deltaT)/ (delX_p*(delX_m+delX_p)) # gamma coefficient | |
ab[2,i-1] = -(2*deltaT)/ (delX_m*(delX_m+delX_p)) # alpha coefficient | |
ab[1,i] = 1 - ab[0,i+1] - ab[2,i-1] # beta coefficient | |
## boundary conditions | |
ab[0,1] = -1 # inner (electrode) boundary | |
ab[1,n-1] = 1 # outer boundary | |
ab[2,n-2] = 0 # outer boundary | |
Theta = theta_i # Set the Potential to the start potential | |
pot = [] # make an empty list to hold all of the potentials equiv to, pot = list() | |
flux = [] # make a list to hold all of the calculated dimensionless fluxes | |
## setup up the simulation loop where we solve the sparse matrix at for each time step | |
# solve linearised problem to obtain concentration profile for each time step | |
for k in range(0,m,1): # for loop which goes from 0 to m-1 | |
if(k<(m/2)): # use logic to determine if we are on the forwar scan | |
Theta += deltaTheta # change the potential on each timestep, equivalent to Theta = Theta - deltaTheta | |
else: # if else cause- so that when on the reverse scan we add to the potential as opposed to subtracting | |
Theta -= deltaTheta | |
fTheta = math.exp(alpha*Theta) # calculate one of the exponential terms in the BV equation | |
ab[1,0] = 1 + h0*fTheta*K0*(1+math.exp(-Theta)) # Set the Butler-Volmer boundary condition- this value changes for each time step | |
C[0] = h0*fTheta*K0*math.exp(-Theta) # Set the surface concentration of the species | |
##### | |
## Solve the banded Matrix | |
## https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve_banded.html | |
#### | |
C = la.solve_banded((1,1),ab,C) # In the UV book the sparse matrix is solved using the Thomas algorithm- here we just use an imported function to do this for us | |
# the function takes the finite difference matrix and the known concentrations and calculates the concentrations for the next time step and stores it back to the C array | |
## | |
## Calculate the flux using a simple two-point approximation | |
flux.append((C[1] - C[0])/(h0)) | |
## Save the Potential | |
pot.append(Theta) | |
#Having completed the simulation record the time | |
finish = datetime.now() | |
pot = thetatoe(pot,Ef) # convert the dimensionless potentials back to a dimensional form | |
def fluxtocurrent(flux): # define a helper function to convert the calculated flux back into dimensional current | |
temp = [] | |
for i in range(len(flux)): | |
temp.append(flux[i]*(np.pi*radius*F*D*conc)) | |
return temp | |
flux = fluxtocurrent(flux) | |
###### | |
## Plotting | |
##### | |
#make data to plot the Randles-Sevich current as a line | |
rspot = [Estart,Evertex] | |
rsip = [iprs,iprs] | |
# plot the expected current from the equation- make the line green and give is a lable | |
plt.plot(rspot,[y*1e6 for y in rsip], color = 'g', label = 'D = %s m2 s-1' %(D)) | |
# plot the current (flux) | |
plt.plot(pot,[y*1e6 for y in flux]) # [y*1e6 for y in flux] - this is list comprehension to multiply each value in a list by 1x10^6 | |
plt.xlabel('Potential / V)') # add x-axis label | |
plt.ylabel('Current / uA') | |
plt.legend(loc = 2) | |
plt.show() | |
# output some basic information about the simulation | |
print('Error on Ip') | |
print(max(flux)/iprs) | |
print(finish-start) # simulation completion time |