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 <[email protected]> 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
>> <[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).
>> 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]> 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
>>>>
>>>> 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]
>> 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]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
