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 used “eqns.solve()” as well, but still it didn’t work. I would
appreciate any guidance that you can offer.
Best regards,
Yuranan
PS. 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)
Yuranan Hanlumyuang, Ph. D
Assist. Professor
Department of Materials Engineering
Kasetsart University
50 Ngamwongwan Rd, Ladyao
Chatuchak, Bangkok 10900
Email: [email protected]
Tel: (66)2-7970999 Ext. 2119
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
