Yuranan - Question: As you can see, shouldn't the value of g be set to -v on the left and +v on the right? The present code gives the solution of phi that is not symmetric. I expect phi to have some symmetric behavior near the left and right boundary.
A: I think you are probably right about the signs. As I recall, I was working through with a right-ward facing normal, so that the signs came out the same, but it should probably be an outward facing normal Questions: How do we treat the boundary condition on the varational terms? Especially, the Robin boundary conditions on this term. I think given a Robin condition: $$ \hat{n}\cdot\left(\vec{a}\phi + b\nabla \phi\right) = g $$ that the boundary condition on the variation would be another Robin condition: $$ \hat{n}\cdot\left(\vec{a}\delts\phi + b\nabla \delta\phi\right) = 0 $$ - Jon > On Mar 7, 2020, at 12:08 AM, Yuranan Hanlumyuang <yurana...@ku.th> wrote: > > 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 > > 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 >> <fipy@nist.gov> 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). >> 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 ] > > _______________________________________________ > 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 ]