Yuranan - There are a few things going on here.
I've modified your notebook and posted it at https://github.com/guyer/sharing-github/blob/master/ion-transport-diagnosis.ipynb. The singular matrix has gone away, but the solution diverges. In general, I find the drift-diffusion equations pretty challenging to solve. Newton-Raphson iterations are almost always called for and under-relaxation or other damping can be needed, too. - Jon > On Feb 24, 2020, at 7:11 PM, Yuranan Hanlumyuang <yurana...@ku.th> wrote: > > Dear the FiPy Team, > > We are trying to solve the coupled equations in the paper titled > "Diffuse-charge dynamics in electrochemical systems" by Bazant, Thornton and > Ajdari using FiPy. I am posting a copy of the python code and the equations > here https://bit.ly/32fIA8a > > The example of Binary_Electrolytes posted by Prof. Guyer has been followed > closely. However, we keep encountering the errors of self.factorizeFnc(*args, > **kwargs) with the RuntimeError: Factor is exactly singular. We were unable > to diagnose if there is anything wrong with the way we set (1) the boundary > conditions, or (2) the coupled equations. > > We tried to use “eqns.solve()” as well, but still it didn’t work. I would > appreciate any guidance that you can offer. > > The code in the text format is listed below: > from fipy import * > > m = Grid1D(nx=100,Lx=2.0) > > v = 1 > epsilon= 0.05 > delta = 0.1 > x = m.cellCenters[0] > > > ##Initial conditions > c = CellVariable(name="c", mesh=m , value=1.0, hasOld=True) > > rho = CellVariable(name="rho", mesh=m , value=0.0, hasOld=True) > > phi = CellVariable(name="phi",mesh=m, hasOld=True) > phi.setValue(v*(x-1), where=m.cellCenters) > > #vi =Viewer((phi, c, rho),datamin=-1.1, datamax=1.1) > #vi.plot() > > > ##Boundary conditions > c.faceGrad.constrain(-rho.faceValue*phi.faceGrad, where=m.facesLeft) > c.faceGrad.constrain(-rho.faceValue*phi.faceGrad, where=m.facesRight) > > rho.faceGrad.constrain(-c.faceValue*phi.faceGrad, where=m.facesLeft) > rho.faceGrad.constrain(-c.faceValue*phi.faceGrad, where=m.facesRight) > > > phi.faceGrad.constrain((phi.faceValue+v)/(epsilon*delta), where=m.facesLeft) > phi.faceGrad.constrain((phi.faceValue-v)/(-epsilon*delta), where=m.facesRight) > > > ##Equations > > Drho = epsilon*rho > Dc = epsilon*c > > eq1 = TransientTerm(var=c) == DiffusionTerm(coeff=epsilon, var=c) > +DiffusionTerm(coeff=Drho, var=phi) > eq2 = TransientTerm(var=rho) == DiffusionTerm(coeff=epsilon, var=rho) > +DiffusionTerm(coeff=Dc, var=phi) > eq3 = DiffusionTerm(coeff=epsilon*epsilon, > var=phi)+ImplicitSourceTerm(coeff=1.0, var=rho) == 0 > eqns = eq1 & eq2 & eq3 > > > ##Solve > #vi =Viewer((phi, c, rho),datamin=-1.1, datamax=1.1) > #from builtins import range > > res = 1e+10 > restol= 1e-13 > > for t in range(100): > c.updateOld() > rho.updateOld() > phi.updateOld() > while res > restol: > res = eqns.sweep(dt=1e-15) > print(res) > > Best, > Yuranan > > > _______________________________________________ > fipy mailing list > fipy@nist.gov > http://www.ctcms.nist.gov/fipy > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] _______________________________________________ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
