Sorry. I never said what the few things were. FiPy has no-flux boundary conditions by default. It's a characteristic of the finite-volume method. So, your Eqs. (3) are satisfied without applying any constraints.
Eqs. (4) and (5) are a [Robin condition](https://www.ctcms.nist.gov/fipy/documentation/USAGE.html#applying-robin-boundary-conditions). I've modified your notebook to apply a Robin condition on \phi. The singular matrix arises because \rho is initialized to zero, which means the coefficient matrix for DiffusionTerm(coeff=Drho, var=phi) is initially zero. The total matrix is still block diagonal, so I don't really understand why this gives a singular matrix, but initializing \rho to something other than zero seems to solve that problem. Unfortunately, as I said before, it still diverges. - Jon > On Feb 25, 2020, at 10:04 AM, Guyer, Jonathan E. Dr. (Fed) via fipy > <fipy@nist.gov> wrote: > > 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://gcc01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fbit.ly%2F32fIA8a&data=02%7C01%7Cjonathan.guyer%40nist.gov%7C18104f76caae4eec145e08d7ba042d06%7C2ab5d82fd8fa4797a93e054655c61dec%7C1%7C1%7C637182399411065840&sdata=h%2FJN%2BFiGcsAGDTNclX%2FqQImXHUgA3%2Fc5LEb5P0F8PTQ%3D&reserved=0 >> >> 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 ] _______________________________________________ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
