Re: FiPy sweep convergence bottoms out

[Daniel Wheeler]

Hi Krishna,

According to this, fixed point iteration only gets you linear
convergence in general,

http://www.math-cs.gordon.edu/courses/ma342/handouts/rate.pdf

I think you need Newton for quadratic convergence.

Another issue to be aware of is that the error and residual are two
separate things. Determining the convergence rate seems to be done
with the error as opposed to the residual. They may well be
interchangeable for these purposes, but you may want to check that.

A further issue is talking about linear versus quadratic convergence
and first, second etc convergence. It seems that one refers to linear
or quadratic in reference to non-linear or linear iterative solutions,
but first, second order convergence when referring to improvements in
accuracy by refining the grid. This is a little confusing.

Cheers,

Daniel





On Tue, Sep 27, 2016 at 7:53 PM, Gopalakrishnan, Krishnakumar
 wrote:
> Hi Ray,
>
>
>
> This looks like linear convergence to me.  Anyway, the bottom-line is that,
> this is too slow. We need something better – is there any acceleration
> routine available?  I tried Jonathan’s gist notebook showing Newton
> implementation using the ResidualTerm, but couldn’t get past a bunch of
> errors from the Python interpreter.
>
>
>
> Or could there be a more fundamental problem in our code
> formulation/structure itself ? The solutions loo correct, when compared to a
> commercial PDE package.
>
>
>
> Krishna
>
>
>
> From: fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] On Behalf Of
> Raymond Smith
> Sent: Wednesday, September 28, 2016 12:40 AM
>
>
> To: fipy@nist.gov
> Subject: Re: FiPy sweep convergence bottoms out
>
>
>
> I am confused as well now. Comparing with the plots on this Wiki page which
> are also semi-log, it looks to me like Krishna is seeing linear convergence.
>
>
>
> On Tue, Sep 27, 2016 at 3:57 PM, Guyer, Jonathan E. Dr. (Fed)
>  wrote:
>
> I don't understand what you mean by "supra-linear trend in the semiology
> plot". You show clear 2nd order convergence, which is what I would expect.
>
>> On Sep 27, 2016, at 4:37 PM, Krishna  wrote:
>>
>> As you can see,  we need supra-linear trend in the semiology plot, such as
>> to continue with the  linear drops achieved in the first few sweeps.
>>
>> I.e. the solver is effectively bottoming out. Under-relaxation factors, or
>> solver-changes don't seem to work.
>>
>> In fact, for certain under-relaxation factors (including 1.0 for 2 of the
>> variables), it breaks the simulation, and produces NaNs right from the first
>> sweep.
>>
>> Krishna
>>
>>
>>
>>  Original Message 
>> From: "Gopalakrishnan, Krishnakumar" 
>> Sent: Tuesday, September 27, 2016 09:02 PM
>> To: fipy@nist.gov
>> Subject: RE: FiPy sweep convergence bottoms out
>>
>> Thank you Ray, Thanks for pointing that out.
>>
>>
>>
>> Here’s the link to Semilog plot. It takes nearly 22 sweeps to achieve a
>> tolerance of 10^-4 for \phi_e and \phi_s_neg.
>>
>>
>>
>> Furthermore, the time spent in sweeping (within each time-step) increases
>> as time progresses.
>>
>>
>>
>> https://imperialcollegelondon.box.com/s/4ix6pozs1h9syt1r3fbkw2pi05ooicmy
>>
>>
>>
>> Krishna
>>
>>
>>
>> From: fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] On Behalf Of
>> Raymond Smith
>> Sent: Tuesday, September 27, 2016 8:51 PM
>> To: fipy@nist.gov
>> Subject: Re: FiPy sweep convergence bottoms out
>>
>>
>>
>> Hi, Krishna.
>>
>> It would be more clear to plot the residuals on a semi-log plot (or
>> equivalently plot the log of residual vs sweep number) to more clearly show
>> the value of the small residuals, as the plots in that link make it look to
>> me like the residuals all go to zero.
>>
>> Ray
>>
>>
>>
>> On Tue, Sep 27, 2016 at 12:42 PM, Gopalakrishnan, Krishnakumar
>>  wrote:
>>
>>
>>
>
>> Hi,
>>
>>
>>
>> We are solving a system of 5 coupled non-linear PDEs.
>>
>>
>>
>> As shown in this plot of residuals vs. sweep count
>> https://imperialcollegelondon.box.com/s/9davbq2gq5eani98xuuj2cw9tmz4mbu3 ,
>> our residuals die down very slowly, i.e. the solver bottoms out. The drop in
>> all the residuals is linear at first, and then asymptotically bottoms out to
>> a value.
>>
>>
>>
>> How do we get our residuals to drop faster, i.e. with lesser sweeps and
>> faster convergence ? I tried changing solvers and tolerances, but curiously
>> enough the results remain identical.
>>
>>
>>
>> Any pointers on this will be much appreciated.
>>
>>
>>
>>
>>
>> Krishna
>>
>>
>>
>>
>>
>>
>>
>>
>> ___
>> 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: 

Re: problem with my initial conditions

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

Good suggestion, but `u.value = ___` and `u.setValue(___)` *should* do the same 
thing in this instance.

> On Oct 5, 2016, at 11:16 AM, Raymond Smith  wrote:
> 
> Hi, Damien.
> 
> Try using u.setValue(0.5 + 0.5*sin(pi*x)*sin(pi*y)) and similarly for v. 
> These lines could naturally be placed right after the definition of u and v 
> (certainly before the time stepping begins).
> 
> See
> http://www.ctcms.nist.gov/fipy/examples/diffusion/generated/examples.diffusion.mesh1D.html
> 
> Ray
> 
> On Wed, Oct 5, 2016 at 5:32 AM, Damien CHUDEAU  
> wrote:
> Hello everybody,
> 
> I am a French student and I am trying to model the emergence of patterns
> on the fur of animals using the Turing equations (involving the
> concentrations u and v) and to code it thanks to the FiPy library. A
> pattern appearing only if the system is unstable, it is necessary that
> the initial conditions on u and v can be changed. Here is my program
> where everything is initialized to 0.5 ... which does not suit me so. I
> have tried to find help on google and on French forums but no one could
> help me. Can anyone explain me how to solve my problem ?
> 
> from fipy import *
> from math import *
> 
> a = 0.25
> b = 0.75
> d = 20
> nx = 20
> ny = nx
> dx = 1.
> dy = dx
> L = dx * nx
> 
> mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)
> 
> u = CellVariable(name = "concentration u", mesh = mesh, hasOld=True,
> value = 0.5)
> v = CellVariable(name = "concentration v", mesh = mesh, hasOld=True,
> value = 0.5)
> 
> D = 1.
> 
> eqn0 = TransientTerm(var=u) == a - u + u*u*v + DiffusionTerm(1, var=u)
> eqn1 = TransientTerm(var=v) == b - u*u*v + DiffusionTerm(D, var=v)
> 
> eqn = eqn0 & eqn1
> 
> #pas d'échange avec l'extérieur -> conditions de Neumann
> 
> u.faceGrad.constrain(((0,),(0,)), where=mesh.facesTop)
> u.faceGrad.constrain(((0,),(0,)), where=mesh.facesBottom)
> u.faceGrad.constrain(((0,),(0,)), where=mesh.facesRight)
> u.faceGrad.constrain(((0,),(0,)), where=mesh.facesLeft)
> v.faceGrad.constrain(((0,),(0,)), where=mesh.facesTop)
> v.faceGrad.constrain(((0,),(0,)), where=mesh.facesBottom)
> v.faceGrad.constrain(((0,),(0,)), where=mesh.facesRight)
> v.faceGrad.constrain(((0,),(0,)), where=mesh.facesLeft)
> 
> vi = Viewer((u, v))
> 
> for t in range(1):
> u.updateOld()
> v.updateOld()
> eqn.solve(dt=1.e-3)
> vi.plot()
> 
> I tried adding the following lines of code ... but it does not work.
> 
> x, y = mesh.cellCenters
> u.value = 0.5 + 0.5*sin(pi*x)*sin(pi*y)
> v.value = 1 + sin(pi*x)*sin(pi*y)
> 
> Thank you in advance
> 
> Damien
> ___
> 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 ]

from math import *

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

I see the standard python math library imported a lot in examples provided by 
fipy users. 

Please don't do this. 

The math library supplies nothing useful for fipy (anything useful in math is 
available from numpy, which we provide as numerix). As seen in the example of 
M. Chudeau, math can override things that fipy depends upon from numpy and 
spicy.

If you insist, please keep clean namespaces with `import math` instead of `from 
math import *`. This is good advice for all your imports. We've hopefully 
removed all instances of `from ___ import *` from our own examples.
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Re: problem with my initial conditions

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

I suspect the problem is that by doing `from math import *`, you are importing 
a version of the `sin` function that only operates on single numbers, not on 
fields. As a result, u and v get assigned uniform values.

`from numpy import *` would be better.
`from fipy.tools.numerix import *` would be better, still.
`from fipy import numerix as nmx` would be best of all.

> On Oct 5, 2016, at 8:32 AM, Damien CHUDEAU  wrote:
> 
> Hello everybody,
> 
> I am a French student and I am trying to model the emergence of patterns 
> on the fur of animals using the Turing equations (involving the 
> concentrations u and v) and to code it thanks to the FiPy library. A 
> pattern appearing only if the system is unstable, it is necessary that 
> the initial conditions on u and v can be changed. Here is my program 
> where everything is initialized to 0.5 ... which does not suit me so. I 
> have tried to find help on google and on French forums but no one could 
> help me. Can anyone explain me how to solve my problem ?
> 
> from fipy import *
> from math import *
> 
> a = 0.25
> b = 0.75
> d = 20
> nx = 20
> ny = nx
> dx = 1.
> dy = dx
> L = dx * nx
> 
> mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)
> 
> u = CellVariable(name = "concentration u", mesh = mesh, hasOld=True, 
> value = 0.5)
> v = CellVariable(name = "concentration v", mesh = mesh, hasOld=True, 
> value = 0.5)
> 
> D = 1.
> 
> eqn0 = TransientTerm(var=u) == a - u + u*u*v + DiffusionTerm(1, var=u)
> eqn1 = TransientTerm(var=v) == b - u*u*v + DiffusionTerm(D, var=v)
> 
> eqn = eqn0 & eqn1
> 
> #pas d'échange avec l'extérieur -> conditions de Neumann
> 
> u.faceGrad.constrain(((0,),(0,)), where=mesh.facesTop)
> u.faceGrad.constrain(((0,),(0,)), where=mesh.facesBottom)
> u.faceGrad.constrain(((0,),(0,)), where=mesh.facesRight)
> u.faceGrad.constrain(((0,),(0,)), where=mesh.facesLeft)
> v.faceGrad.constrain(((0,),(0,)), where=mesh.facesTop)
> v.faceGrad.constrain(((0,),(0,)), where=mesh.facesBottom)
> v.faceGrad.constrain(((0,),(0,)), where=mesh.facesRight)
> v.faceGrad.constrain(((0,),(0,)), where=mesh.facesLeft)
> 
> vi = Viewer((u, v))
> 
> for t in range(1):
>u.updateOld()
>v.updateOld()
>eqn.solve(dt=1.e-3)
>vi.plot()
> 
> I tried adding the following lines of code ... but it does not work.
> 
> x, y = mesh.cellCenters
> u.value = 0.5 + 0.5*sin(pi*x)*sin(pi*y)
> v.value = 1 + sin(pi*x)*sin(pi*y)
> 
> Thank you in advance
> 
> Damien
> ___
> 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: problem with my initial conditions

[Raymond Smith]

Hi, Damien.

Try using u.setValue(0.5 + 0.5*sin(pi*x)*sin(pi*y)) and similarly for v.
These lines could naturally be placed right after the definition of u and v
(certainly before the time stepping begins).

See
http://www.ctcms.nist.gov/fipy/examples/diffusion/generated/examples.diffusion.mesh1D.html

Ray

On Wed, Oct 5, 2016 at 5:32 AM, Damien CHUDEAU 
wrote:

> Hello everybody,
>
> I am a French student and I am trying to model the emergence of patterns
> on the fur of animals using the Turing equations (involving the
> concentrations u and v) and to code it thanks to the FiPy library. A
> pattern appearing only if the system is unstable, it is necessary that
> the initial conditions on u and v can be changed. Here is my program
> where everything is initialized to 0.5 ... which does not suit me so. I
> have tried to find help on google and on French forums but no one could
> help me. Can anyone explain me how to solve my problem ?
>
> from fipy import *
> from math import *
>
> a = 0.25
> b = 0.75
> d = 20
> nx = 20
> ny = nx
> dx = 1.
> dy = dx
> L = dx * nx
>
> mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)
>
> u = CellVariable(name = "concentration u", mesh = mesh, hasOld=True,
> value = 0.5)
> v = CellVariable(name = "concentration v", mesh = mesh, hasOld=True,
> value = 0.5)
>
> D = 1.
>
> eqn0 = TransientTerm(var=u) == a - u + u*u*v + DiffusionTerm(1, var=u)
> eqn1 = TransientTerm(var=v) == b - u*u*v + DiffusionTerm(D, var=v)
>
> eqn = eqn0 & eqn1
>
> #pas d'échange avec l'extérieur -> conditions de Neumann
>
> u.faceGrad.constrain(((0,),(0,)), where=mesh.facesTop)
> u.faceGrad.constrain(((0,),(0,)), where=mesh.facesBottom)
> u.faceGrad.constrain(((0,),(0,)), where=mesh.facesRight)
> u.faceGrad.constrain(((0,),(0,)), where=mesh.facesLeft)
> v.faceGrad.constrain(((0,),(0,)), where=mesh.facesTop)
> v.faceGrad.constrain(((0,),(0,)), where=mesh.facesBottom)
> v.faceGrad.constrain(((0,),(0,)), where=mesh.facesRight)
> v.faceGrad.constrain(((0,),(0,)), where=mesh.facesLeft)
>
> vi = Viewer((u, v))
>
> for t in range(1):
> u.updateOld()
> v.updateOld()
> eqn.solve(dt=1.e-3)
> vi.plot()
>
> I tried adding the following lines of code ... but it does not work.
>
> x, y = mesh.cellCenters
> u.value = 0.5 + 0.5*sin(pi*x)*sin(pi*y)
> v.value = 1 + sin(pi*x)*sin(pi*y)
>
> Thank you in advance
>
> Damien
> ___
> 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 ]

problem with my initial conditions

[Damien CHUDEAU]

Hello everybody,

I am a French student and I am trying to model the emergence of patterns 
on the fur of animals using the Turing equations (involving the 
concentrations u and v) and to code it thanks to the FiPy library. A 
pattern appearing only if the system is unstable, it is necessary that 
the initial conditions on u and v can be changed. Here is my program 
where everything is initialized to 0.5 ... which does not suit me so. I 
have tried to find help on google and on French forums but no one could 
help me. Can anyone explain me how to solve my problem ?

from fipy import *
from math import *

a = 0.25
b = 0.75
d = 20
nx = 20
ny = nx
dx = 1.
dy = dx
L = dx * nx

mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)

u = CellVariable(name = "concentration u", mesh = mesh, hasOld=True, 
value = 0.5)
v = CellVariable(name = "concentration v", mesh = mesh, hasOld=True, 
value = 0.5)

D = 1.

eqn0 = TransientTerm(var=u) == a - u + u*u*v + DiffusionTerm(1, var=u)
eqn1 = TransientTerm(var=v) == b - u*u*v + DiffusionTerm(D, var=v)

eqn = eqn0 & eqn1

#pas d'échange avec l'extérieur -> conditions de Neumann

u.faceGrad.constrain(((0,),(0,)), where=mesh.facesTop)
u.faceGrad.constrain(((0,),(0,)), where=mesh.facesBottom)
u.faceGrad.constrain(((0,),(0,)), where=mesh.facesRight)
u.faceGrad.constrain(((0,),(0,)), where=mesh.facesLeft)
v.faceGrad.constrain(((0,),(0,)), where=mesh.facesTop)
v.faceGrad.constrain(((0,),(0,)), where=mesh.facesBottom)
v.faceGrad.constrain(((0,),(0,)), where=mesh.facesRight)
v.faceGrad.constrain(((0,),(0,)), where=mesh.facesLeft)

vi = Viewer((u, v))

for t in range(1):
u.updateOld()
v.updateOld()
eqn.solve(dt=1.e-3)
vi.plot()

I tried adding the following lines of code ... but it does not work.

x, y = mesh.cellCenters
u.value = 0.5 + 0.5*sin(pi*x)*sin(pi*y)
v.value = 1 + sin(pi*x)*sin(pi*y)

Thank you in advance

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