Hi FiPy team, I have a few subsequent questions about the (1) the Robin boundary condition, and (2) the Newton-Method implementation. I listed the questions in
https://bit.ly/2TI8jlW <https://bit.ly/2TI8jlW> along with some relevant equations. The questions concern (1) the values of ‘g’ in the Robin boundary conditions, and (2) the boundary condition for the variational of parameters in the Newton method. I would be greatly appreciated any guidance you may offer regarding other points as well. Best, Yuranan > On Feb 25, 2020, at 10:32 PM, Guyer, Jonathan E. Dr. (Fed) via fipy > <[email protected]> wrote: > > 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 > > <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 >> <[email protected] <mailto:[email protected]>> 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 <[email protected]> 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 >>> >>> <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 >>> [email protected] >>> http://www.ctcms.nist.gov/fipy >>> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] >> >> >> _______________________________________________ >> fipy mailing list >> [email protected] >> http://www.ctcms.nist.gov/fipy >> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] > > > _______________________________________________ > fipy mailing list > [email protected] <mailto:[email protected]> > http://www.ctcms.nist.gov/fipy <http://www.ctcms.nist.gov/fipy> > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy > <https://email.nist.gov/mailman/listinfo/fipy> ]
_______________________________________________ fipy mailing list [email protected] http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
