Merge pull request #267 from ludwig-cf/prerelease-0.20.0

Prerelease 0.20.0

CHANGES.md

@@ -3,11 +3,17 @@
 
 version 0.20.0
 
+- IMPORTANT: The input file can no longer be specified as a command
+  line argument. It must be called "input" in the current directory.
+- The electrokinetics sector has been updated and information is
+  available at
+  https://ludwig.epcc.ed.ac.uk/tutorials/electrokinetics/electrokinetics.html
+- A D3Q27 model is available
 - Added coverage via https://about.codecov.io/
 - The free energy is now reported at t = 0 for the initial state.
 - An extra "first touch" option has been added. For details, see
   https://ludwig.epcc.ed.ac.uk/inputs/parallel.html
-
+- Various minor changes and code quality improvements.
 
 version 0.19.1
 - Fix bug in io_subfile to prevent failure at iogrid > 1.

README.md

@@ -44,7 +44,8 @@ $ make test
 ```
 
 
-Full details of the build process are available at
+Full details of the build process, and tutorials on how to
+use the code are available at
 <a href = "https://ludwg.epcc.ed.ac.uk/">https://ludwig.epcc.ed.ac.uk/</a>.
 
 #### Background

docs/electrokinetic.tex

@@ -273,180 +273,6 @@ \subsubsection{Gouy-Chapman}
 Debye length $l_D=3.514$, the surface potential $\Psi_D=2.267\e{-4}$
 and the centre potential $\Psi_c=-3.395\e{-5}$. 
 
-\subsubsection{Liquid-junction potential}
-
-The liquid junction potential 
-is a charge separation process that 
-occurs when electrolytes with slightly different concentrations
-whose species have different diffusivities are brought into contact.
-Charges from the regions of higher concentration diffuse   
-into the parts with lower concentration. Due to the difference 
-in diffusivity they migrate at different speeds, leaving parts of
-the system charged. This leads to a build-up of a potential
-which balances the diffusive flux.
-
-After the initial build-up phase the potential decreases slowly 
-again until the charge concentration has become homogeneous throughout 
-the system. Both timescales of emergence and decay of the potential
-can be separated by chosing a sufficiently large system size.
-
-This problem allowed us to verify the correct temporal 
-behaviour of the Nernst-Planck equation solver by resolving the transient 
-dynamics without having to account for advective terms.
-
-\begin{figure}[h!t]
-\includegraphics[width=0.495\textwidth]{./pics/test_lj_zoom1.pdf}
-\includegraphics[width=0.495\textwidth]{./pics/test_lj_zoom3.pdf}
-\caption{Time evolution of the liquid junction potential for $L_x=128$ (blue), $L_x=256$ (green) and $L_x=512$ (blue). The dashed black curves represent the approximate solution in the limit $l_D/L_x<<1$. Deviations can be seen which are presumably due
-to the approximate nature of the analytical solution and the fact that we
-have a different asymptotic flux close to the boundary - the theoretical 
-analysis assumes no flux at a boundary which is infinitely far away from
-the diffusive region.} 
-\label{fig3} 
-\end{figure}
-
-For simplicity we considered systems of size 
-$L_x\times L_y\times Lz=128\times4\times4$ and 
-$L_x\times L_y\times Lz=256\times4\times4$ with
-periodic boundary conditions at either end.
-The two halfs were electroneutral and had ionic concentrations 
-$\rho_{L,\pm}=\rho_{0,\pm} + \delta\rho$ and 
-$\rho_{R,\pm}=\rho_{0,\pm} - \delta\rho$ 
-with $\rho_{0,\pm}=1\e{-2}$ and $\delta\rho = 0.01$.
-
-The potential difference between both sides during the build-up 
-obeys approximately
-
-\begin{equation}
-\Delta\Psi(t)\simeq\Delta\Psi_P \left\{1-\exp\left(-\frac{t}{\tau_e}\right)\right\}\\
-\end{equation}
-
-with 
-
-\begin{equation}
-\Delta\Psi_P=\frac{(D_+ D_-)}{\beta e (D_+ + D_-)} \frac{\delta\rho}{\rho_0}\\ 
-\end{equation}
-
-as saturation value of the potential difference.
-The saturation time scale is given by
-
-\begin{equation}
-\tau_e=\frac{\varepsilon}{\beta \, e^2 (D_+ + D_-) \rho_0}.
-\end{equation}
-
-\begin{figure}[htpb]
-\includegraphics[width=0.9\textwidth]{./pics/test_lj_decay1.pdf}
-\caption{Rescaled plot of the decay: Times have been Time evolution of the liquid junction potential for $L_x=128$ (blue), $L_x=256$ (green) and $L_x=512$ (blue). The dashed black curves represent the approximate solution in the limit $l_D/L_x<<1$. Deviations can be seen which are due
-to the approximate nature of the analytical solution and the diffusive fluxes close to the boundary.} 
-\label{fig4} 
-\end{figure}
-
-A more exact solution can be derived in the limit $l_D/L_x<<1, N\to\infty$: 
-
-\begin{eqnarray}
-\Delta\Psi(t)&=&\Delta\Psi_P \left\{1-\exp\left(-\frac{t}{\tau_e}\right)\right\}\frac{4}{\pi}\left\{\sum_{m=1}^N \frac{\sin^3(m\pi/2)}{m} \exp\left(-\frac{ m^2\, t}{\tau_d}\right)\right\}\\
-\tau_d&=&\frac{L^2}{2\pi^2 (D_+ + D_-)}.\label{taud}
-\end{eqnarray}
-
-It contains as well the dependence on the decay timescale $\tau_d$.
-Only odd indices $m$ contribute to the sum:
-
-\begin{eqnarray}
-\sum_{m=1}^N \frac{\sin^3(m\pi/2)}{m} \exp\left(-\frac{m^2\, t}{\tau_d}\right)&=&\nonumber\\
-&&\hspace*{-4cm} \exp\left(-\frac{t}{\tau_d}\right)-\frac{1}{3} \exp\left(-\frac{9 t}{\tau_d}\right)+\frac{1}{5}\exp\left(-\frac{25t}{\tau_d}\right)-\frac{1}{9}\exp\left(-\frac{81 t}{\tau_d}\right)+\ldots
-\end{eqnarray}
-
-A complete discussion of the solution can be found in \cite{Mafe}. 
-There, the upper limit of significant modes has been also estimated as $N_{max} = L/\pi l_D$.
-Note the factor 2 difference between Eq. \ref{taud} and the corresponding expression in \cite{Mafe}.
-
-The following parameters were used:
-dielectric permittivity $\epsilon=3.3\e{3}$, temperature $\beta^{-1}=3.333\e{-5}$, unit charge $e=1$, valency $z_\pm=\pm1$, diffusivities $D_+=0.0125$ and $D_-=0.0075$.
-We obtained
-$\Delta\Psi_P=1.6667\e{-7}$, $\tau_e=550$, $\tau_d=41501.2\, (L_x=128)$, $\tau_d=166004.6\, (L_x=256)$ and $\tau_d=664018.5 (L_x=512)$.
-The results for the potential difference over time are shown in Fig. \ref{fig3}.
-
-Fig. \ref{fig4} shows results with times rescaled to the decay 
-time scale $\tau_d$ (cf. Eq. \ref{taud}). Obviously the 
-deviations we observe are not due to the limited system size 
-and have a more systematic origin. 
-
-The curves coincide if 
-the theoretic limit for $\tau_d$ is rescaled by a factor $1.067$,
-suggesting the effective system length for this sort of setup is
-actually about 3\% larger than the numerical value.
-
-A reason for this might be the approximate nature of the analytical solution 
-and the fact that it was gained
-for an infinitely large system with constant charge concentrations,
-vanishing currents at both ends and finite diffusive zone of size $L_x$.
-In our situation the entire system is within the diffusive zone.
-This may lead to smaller effective diffusivities or larger effective
-system sizes.
-Interestingly, all runs with solid walls at both ends resulted in 
-oscillatory behaviour and an effective system size of $2L_x$.
-
-
-\subsubsection{Electroosmotic Flow}
-
-To test the implementation with all couplings to external and 
-internal forces we consider a forced charged fluid in a slit
-of size $L_z$. An electrostatic field $E_{||}$ is allied
-parallel to the walls. The entire system is electroneutral with 
-each wall having the surface charge density $\sigma$ 
-and compensationg counterions with total charge $2 \sigma A_{wall}$
-in the fluid.
-
-In equilibrium the charge density at a distance $x$ from the wall obeys
-
-\begin{equation}
-\rho(x)=\frac{\rho_0}{\cos^2(K\,x)}
-\end{equation}
-
-with 
-
-\begin{equation}
-\rho_0=K^2/2\pi l_B
-\end{equation}
-
-and 
-
-\begin{eqnarray}
-K \,L_x \tan\left(\frac{K\, L_x}{2}\right)&=&\pi\, l_B\, L_x\, 2\sigma\label{kex} \\
-K \,L_x&\simeq&\sqrt{4\pi \,l_B\,L_x\,2\sigma}\label{klin}.
-\end{eqnarray}
-
-
-The liniarised version Eq. \ref{klin} has only a limited range of applicabilty.
-We solved Eq. \ref{kex} numerically and found solutions 
-$K=0.01959\; (\sigma=0.003125)$ and $K=0.03311\; (\sigma=0.00125)$, 
-which is reasonably far away from the 
-theoretical limit $K_{max}=\pi/L_x$ set by the tangent. 
-
-Note the factor 2 difference on the lhs of Eq. \ref{kex} with respect 
-to \cite{Capuani, Rotenberg}. There is also a factor $L_x$ missing on 
-the rhs of Eq. \ref{klin}.
-
-The steady state velocity of the fluid can be derived from the 
-force balance of the gradient of the stresses and the electrostatic
-forces:
-
-\begin{eqnarray}
-v_y(x)&=&\hat{v} \ln\left(\frac{\cos(K\,x)}{\cos(K\,L_x/2)}\right)\label{vy}\\
-\hat{v}&=&\frac{e \,E_{||}\rho_o}{\eta\, K^2}=\frac{e \,E_{||}}{2\pi\eta l_B}\label{vhat}
-\end{eqnarray}
-
-The result for two different charge densities is shown in Fig. \ref{fig6}.
-The accuracy is acceptable with deviations for high surface 
-charged potentially being caused by the chosen discretisation or 
-by the numerical solution of Eq. \ref{kex} approaching the limit 
-of $\pi/L_z\simeq0.049$.  
-
-\begin{figure}[htpb]
-\includegraphics[width=0.9\textwidth]{./pics/test_eo.pdf}
-\caption{Steady state flow profile across a slit of width $L_x=63$ for an applied field of magnitude $e \beta E L_x=1.89$ for two different charge densities $\sigma=0.003125$ (red) and $\sigma=0.0125$ (blue), Bjerrum length $l_B=0.7234$, viscosity $\eta=0.1$ and unit charge $e=1$. The dashed black lines correspond to theoretical prediction according to Eqs. \ref{vy} and \ref{vhat}.} 
-\label{fig6} 
-\end{figure}
 
 \subsubsection{Debye-H\"uckel Theory}