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 <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 ]
