Re: Charge dynamics and Neumann boundary conditions

[Yuranan Hanlumyuang]

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 
>  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  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 
>>>  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 
  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  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.   
> 

Re: Charge dynamics and Neumann boundary conditions

[Guyer, Jonathan E. Dr. (Fed) via fipy]

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  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 
>>  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 
>>>  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  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.faceVa

Re: Charge dynamics and Neumann boundary conditions

[Yuranan Hanlumyuang]

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 
>  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 
>> mailto: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  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

Re: Charge dynamics and Neumann boundary conditions

[Yuranan Hanlumyuang]

Hi Jon,

Since I am new to FiPy, it may take me awhile to follow through your solution. 
I am sure will have further questions. In the mean time, thank you so much for 
your prompt response.

Best regard,
Yuranan




___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Re: Charge dynamics and Neumann boundary conditions

[Guyer, Jonathan E. Dr. (Fed) via fipy]

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 
>  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  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 ]

Re: Charge dynamics and Neumann boundary conditions

[Guyer, Jonathan E. Dr. (Fed) via fipy]

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  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://bit.ly/32fIA8a
> 
> 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 ]

Charge dynamics and Neumann boundary conditions

[Yuranan Hanlumyuang]

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://bit.ly/32fIA8a 

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 ]