The governing equations are listed here of my notes on page 4. It's a reproduction of other's paper which solves the equations with COMSOL. The problems arise when I want to solve for the consistent initial condition of the algebraic components, $\phi_1$ and $\phi_2$. They consist of a kinetic reaction term for five reactants
\begin{equation} a \sum_j i_j \end{equation} \begin{equation} i_j = i_{o,j}\biggl[ \prod_i \biggl(\frac{C_i}{C_{i,o}}\biggr)^{p_{i,j}}\exp\biggl(\frac{\alpha_{aj} F}{RT}(\eta-U_j)\biggr) - \prod_i \biggl(\frac{C_i}{C_{i,o}}\biggr)^{q_{i,j}}\exp\biggl(-\frac{\alpha_{cj} F}{RT}(\eta-U_j)\biggr) \biggr] \end{equation} where $\eta=\phi_1-\phi_2$.
First try for these two algebraic equations on MINPACK gives $l^2$ norm of the residual around $1.0E-2$ for one dimension case $=500$. Then I twisted the equations so they become one equation of variable $\eta$ (not two coupled equations), the equations listed here on page 4. Remarkably the $l^2$ norm drops to $1.0E-5$. Nevertheless, when I integrate over time they failed on the very first step. It seems that when the equations are all coupled together with this kinetic term, the convergence of the Newton's method failed (I've only tried on the solver DASKR (DASSL's variant) and RADAU5). I'm using method of line and the space discretization is second order. And the Neumann boundary conditions are all coupled back to the governing equations with the ghost points. The Jacobian matrices are generated by Tapenade. The problem is very close to this one.
Furthermore I think (with my little numerical knowledge) the equation of this kind \begin{equation} \sigma \nabla^2 \phi_1 = a\sum_j i_j \end{equation} resembles a singular perturbed equation for the coefficient $\sigma/a$ can be of order $1.0E-10$ on the sulfur cathode. $a$ is the surface area of the pore walls per unit volume of the total electrode. From the output data at the boundary there is a one point drop from a relatively smooth inner domain. I'm thinking pseudo-spectral method might solve the problem?
Another question is that is it perfectly fine to impose boundary conditions (specifically Neumann type) on PDE as a discretized algebraic equations without putting back into the governing equation with ghost point? I feel like somehow the boundary point won't satisfy the governing equation this way. And when implementing in the IDE or DAE solver I need to explicitly specify they are algebraic term though.
One last question: do the initial conditions need to satisfy the boundary conditions? I think most system don't. Or else I have to let the applied current to be zero first and then at $t>0$ crank up to some value, making a discontinuity at the first step. I've found one paper describe this discontinuous initial condition, DAEs not ODEs page 371, but I still don't know the technique in solving them.
I've been working this whole summer vacation from scratch and this is my very first time solving a real (non-toy) numerical problems. Any suggestion would be greatly appreciated.
EDIT
Quad precision without precondition just solves my first problem regarding ill-conditioned matrix. The third problem (inconsistent initial and boundary condition) is also solved by making initial guessing on the algebraic variables sufficiently close to the values when the source term jumps. But the problem is quad precision implementation is more than 10x slower than 8 bytes double precision. Will it help if I reorder the Jacobian matrix so that it is banded (currently the linear solver is LAPACK with dense matrix)? Or I just need to learn PETSc and make use of the parallel feature and sparse solver? Also, how fined should it be so that the solution can considered fine-grid solution (and hence make error estimation from that)? The second problem above in the main context still puzzles me though.
EDIT.II
Reaction current $\nabla \cdot i_2$ figure from a review paper on porous electrode. The author define a dimensionless current density $\delta = \frac{\alpha_a FIL}{RT}(\frac{1}{\kappa}+\frac{1}{\sigma})$ which monitor how non-uniform the current distribution can be. When $\delta=100$, the current at the boundary is very steep already. The system I studied has even higher of this parameter though.
EDIT.III
I solve the discontinuity problem regarding the Poisson equation. From the dimensionless parameter above, there is a term $L$. I put $1$ in my code which is non-physical for it means 1m electrode length for consistent unit. After lowering it down to around $1E-5$, the source term becomes smooth. Thanks for Geoff's a priori thought that made me look closer to the physical quantity.
