Hi Jon,
I tried to fix the the value of +v and -v on the right and left by the
following code (for the Robin boundary condition)
(RobinCoeff * g*n).divergence
- fp.ImplicitSourceTerm(coeff=(RobinCoeff).divergence, var=phi)
and I also added
-fp.ImplicitSourceTerm(coeff=(RobinCoeff).divergence, var=dphi)
for the Robin condition ford phi. The code is in
https://bit.ly/3cTWlye
However, the code gave me run time error. Do you know if there is anything
wrong with the code? Your insight into the problem will be greatly appreciated.
Best,
Yuranan
> On Mar 11, 2020, at 3:49 AM, Guyer, Jonathan E. Dr. (Fed) via fipy
> <fipy@nist.gov> wrote:
>
> 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 ]
_______________________________________________
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
