Re: reaction-advection-diffusion; problems making the solution converge

[Nils Becker]

hi daniel, thanks for looking into this!


> Hi Nils,
> 
> It seems like the code is converging well as far as I can tell. For
> example, when I modify the sweep loop to be
> 
> 
> print
> print 'step: {0}'.format(tt)
> while res > 1e-9: #Solve.sweeps:
> res = eq.sweep(var=rho, dt=Solve.step,
> solver=Solve.customSolver)
> print res
> 
> 
> every step seems to converge in 5 steps to 1e-10.

ok. it seems that you have not changed the iteration, only added print
functions. so my code seems to be converging, yes?

that's also what i thought. however testing the solution externally, it
seems that the solution is not quite accurate (as time progresses, it
seems to get worse).

of course, this external test could be buggy -- however the fact that it
shows that lhs = rhs quite accurately in the beginning makes me think
that it could be correct.

> The non-linear residual only really tells you how well each step
> individually is converging. It seems like you were only using the very
> first non-linear residual, which uses the old value for the
> calculation, which doesn't tell you anything useful. 

so what sweep returns is basically A(n) * x(n-1) - b(n), yes? as far as
i can tell, i always collect the last one of these values in each step.
they do get small, i also see the convergence to below 1e-9 when i set
the tolerance that low.

> 
> One thing that I'm confused about is that it requires 5 sweeps to
> converge for a fully linear equation. Is it actually linear? I
> couldn't tell from the code.

yes. this also surprised me. the potential phi and the growth rate
lambda_ are of course non-linear functions of space. but all
coefficients are independent of the density, and the density rho should
be really appearing only linearly in the equation. so unless i
constructed the equation or the loop wrong, the system should really be
linear. no idea why it would require sweeps. do linear equations with an
implicit source term require sweeping? am i doing the sweep wrong?

n.







0xF0A5C638.asc
Description: application/pgp-keys
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Re: reaction-advection-diffusion; problems making the solution converge

[Daniel Wheeler]

Hi Nils,

It seems like the code is converging well as far as I can tell. For
example, when I modify the sweep loop to be


print
print 'step: {0}'.format(tt)
while res > 1e-9: #Solve.sweeps:
res = eq.sweep(var=rho, dt=Solve.step,
solver=Solve.customSolver)
print res


every step seems to converge in 5 steps to 1e-10.

step: 20.0
0.00192315625797
4.5626520206e-05
1.08935792171e-06
2.60975675747e-08
6.26480480874e-10

step: 21.0
0.00195977695607
4.651749367e-05
1.11103436368e-06
2.66255569394e-08
6.3934852818e-10

step: 22.0
0.00199768415442
4.74384979099e-05
1.13342888153e-06
2.71708361386e-08
6.52634229319e-10

The non-linear residual only really tells you how well each step
individually is converging. It seems like you were only using the very
first non-linear residual, which uses the old value for the
calculation, which doesn't tell you anything useful. Basically "A_new
* x_old -b_new" is not a measure of anything useful. It is only useful
as a quantity for normalization. FiPy uses the old value or the
previous sweep value to calculate the non-linear residual. So the
residual is

   A(n) * x(n-1) - b(n)

where n is the sweep number, and x_(-1) is x_old if the first sweep is
n=0. If only one sweep is executed then the residual isn't a useful
quantity. The linear residual is not returned from the "sweep" method,
which is "A(n) * x(n) - b(n)" and is assumed to be as accurate as the
tolerance of the linear solve. Basically, collecting the first
non-linear residual at each time step is not a useful exercise to
generate a metric for convergence across time steps. In fact, I'm not
even sure how to do that or what it means. I only know how to
understand convergence for a given time step.

One thing that I'm confused about is that it requires 5 sweeps to
converge for a fully linear equation. Is it actually linear? I
couldn't tell from the code.

Cheers,

Daniel




On Fri, Aug 5, 2016 at 8:30 AM, Nils Becker
 wrote:
>>> I'm not sure if using ||lhs|| is a good way to normalize the residual.
>>> Surely that could be very small. Is ||lhs - rhs|| getting smaller
>>> without normalization?
>
> if you don't see any obvious errors in building the solution in fipy, i
> could continue poking around if i knew how to evaluate the various terms
> of the solution and their errors at different times. i see methods
> eq_diffprolif.justErrorvector, but i don't quite understand what exactly
> this returns. is it lhs-rhs? what is justResidualVector? can i get those
> for the individual terms on the rhs as well?
>
> cheers, n.
>
> ___
> fipy mailing list
> fipy@nist.gov
> http://www.ctcms.nist.gov/fipy
>   [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>



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

Re: reaction-advection-diffusion; problems making the solution converge

[Nils Becker]

>> I'm not sure if using ||lhs|| is a good way to normalize the residual.
>> Surely that could be very small. Is ||lhs - rhs|| getting smaller
>> without normalization?

if you don't see any obvious errors in building the solution in fipy, i
could continue poking around if i knew how to evaluate the various terms
of the solution and their errors at different times. i see methods
eq_diffprolif.justErrorvector, but i don't quite understand what exactly
this returns. is it lhs-rhs? what is justResidualVector? can i get those
for the individual terms on the rhs as well?

cheers, n.


0xF0A5C638.asc
Description: application/pgp-keys
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Re: reaction-advection-diffusion; problems making the solution converge

[Nils Becker]

> I'm not sure if using ||lhs|| is a good way to normalize the residual.
> Surely that could be very small. Is ||lhs - rhs|| getting smaller
> without normalization?

ok, i was not very precise in my statement. actually i was doing this:

1. make a plot of lhs pointwise.
2. make a plot of lhs - rhs pointwise, on the same scale.

then compare the size of the deviations lhs-rhs. what i see is that the
absolute deviations lhs-rhs do increase over time and that they reach a
level over 10% of the lhs, at those points in the phase space where the
lhs is largest in magnitude

sorry for the confusion.

n.


0xF0A5C638.asc
Description: application/pgp-keys
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Re: reaction-advection-diffusion; problems making the solution converge

[Daniel Wheeler]

On Wed, Aug 3, 2016 at 9:27 AM, Nils Becker
 wrote:
>> Could you possible post your code so I can play around with it? I
>> can't think offhand what could be the issue or exactly what the error
>> is that you're calculating, it seems like the normalized residual or
>> something.
>
> hi daniel, thanks for your reply. the way i was calculating the residual
> is to export the solution on a grid, import that into mathematica,
> interpolate in space and time and have mathematica calculate
> ||lhs-rhs|| / ||lhs|| where ||...|| is the L2 norm. what i saw was that
> the relative error was low in the beginning but increased over time,
> especially in places where the solution starts to grow exponentially.

I'm not sure if using ||lhs|| is a good way to normalize the residual.
Surely that could be very small. Is ||lhs - rhs|| getting smaller
without normalization?

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

Re: reaction-advection-diffusion; problems making the solution converge

[Nils Becker]

> Could you possible post your code so I can play around with it? I
> can't think offhand what could be the issue or exactly what the error
> is that you're calculating, it seems like the normalized residual or
> something. 

hi daniel, thanks for your reply. the way i was calculating the residual
is to export the solution on a grid, import that into mathematica,
interpolate in space and time and have mathematica calculate
||lhs-rhs|| / ||lhs|| where ||...|| is the L2 norm. what i saw was that
the relative error was low in the beginning but increased over time,
especially in places where the solution starts to grow exponentially.


i am sure this relative error can be gotten from fipy using
residualVector or some other methods of the terms but i have not been
able to understand the docs well enough to replicate that.

so i am posting pretty much the (streamlined but otherwise unchanged)
example script i am having trouble with below. h5py is needed only for
the export, you can comment that out if you dont export. Lambda3 is the
source term coefficient -- by increasing it's values the solution
becomes more rapidly growing and the errors also increase.

i hope the indentation survives emailing

cheers, n.



from __future__ import division
import fipy as fp
from fipy import numerix as nx
# the module for the bare numpy functions
# this is essentially numpy!
NX = nx.NUMERIX
from time import sleep


# this is where we are going to write results
import os
datadir = os.path.join(
os.path.dirname(os.path.realpath(__file__)),
'pde_test_data/'
)
datadir = os.path.realpath(datadir)
try: os.mkdir(datadir)
except OSError: pass

# utility
def nsq(v):
return nx.dot(v, v)

class Domain(object):
upper = 300.
lower = 0
size = upper - lower

class Mesh(Domain):
def __init__(self):
delta = 1.5
Lx, Ly = (self.size,) * 2
dx, dy = (delta,) * 2
nx, ny = (self.size // delta,) * 2
self.mesh = fp.Grid2D(dx=dx, dy=dx, nx=nx, ny=ny, Lx=Lx, Ly=Ly)

# from Mma, the last iteration of the potential was this:

class Phi(object):
scale = 100.
x1, y1 = [100., 100]
x2, y2 = [180., 180]
s1 = 4.
s2 = 12.
dp = 3.
strength = 30.
sharpness_a = 120.
# the following should act on the position field x
# this should be a cell varable in fipy.
@classmethod
def itp(cls, z):
t = nx.tanh(z/cls.sharpness_a**2) + 1
return 0.5 * t
@classmethod
def scalephi(cls, x, y):
[d1, d2] = [(x-cls.x1)**2 + (y-cls.y1)**2,
(x-cls.x2)**2 + (y-cls.y2)**2]
[f2, f1] = [cls.itp(d1 - d2), cls.itp(d2 - d1)]
p1 = 0.5 * cls.s1 * d1 / cls.scale**2
p2 = 0.5 * cls.s2 * d2 / cls.scale**2 - cls.dp
return f1 * p1 + f2 * p2
@classmethod
def __call__(cls, x, y):
return cls.strength * cls.scalephi(x, y)

class Lambda3(object):
def __init__(self, curvature=None, direction=None, start=None):
if curvature is None: curvature = 0.2 / Domain.size**2
if direction is None: direction = nx.array([1, 1.])
direction /= nx.sqrt(nsq(nx.array(direction)))
if start is None: start = nx.array([0.,0])
max_pos = 0.5 * nx.sqrt(2) * Domain.size
max_ht = 1./2 * max_pos**2 * curvature # for 0 at the origin
max_ht = 0.1 * max_ht # tame the beast
(self.curvature, self.direction, self.max_pos, self.max_ht,
self.start) = (
curvature, direction, max_pos, max_ht, start)
def __call__(self, x, y):
distance2 = (NX.dot(self.direction,
(nx.array([x, y]).T - self.start).T) - self.max_pos)**2
lam = self.max_ht - 1./2 * self.curvature * distance2
return lam

class InitialRho(object):
'''Gaussian-like initial condition'''
c = [[100], [100]]
covm = NX.eye(2) * 30**2  # bare numpy
stm = NX.linalg.inv(covm)
Z = (2. * NX.linalg.det(covm))**(1./2)
@classmethod
def __call__(cls, x, y):
p = NX.array([x,y])
return  1. / cls.Z * NX.exp(-1./2 * NX.diag(
cls.stm.dot(p - cls.c).T.dot(p - cls.c)))

mesh = Mesh().mesh
x, y = mesh.cellCenters
x1d, y1d = x[:mesh.nx], y[::mesh.ny] # yes!


# diffusion
Dtensor = nx.eye(2) * 10.
# constant diff coeff fields
Dtensor = fp.CellVariable(mesh, 'dtensor', Dtensor)

# make the input fields
phi = fp.CellVariable(name=r'$\phi$', mesh=mesh, value=Phi()(x, y))
v = - phi.grad  # on the cells, not the faces.

lambda_ = fp.CellVariable(name=r'$\lambda$', mesh=mesh,
value=Lambda3(curvature=1./Domain.size**2, direction=[1.,-1],
start=[0.,Domain.size])(x, y))


# terms
def diff(DD):
# Correction actually needed? in any case it seems to work.
# list of diffusion coefficients of different orders
return fp.DiffusionTermCorrection(coeff=[DD,])
def drift(v):
return - fp.PowerLawConvectionTerm(coeff=v)
def prolif(lambda_):
# implicit = multiplies the var
return 

Re: reaction-advection-diffusion; problems making the solution converge

[Daniel Wheeler]

Hi Nils,

Could you possible post your code so I can play around with it? I
can't think offhand what could be the issue or exactly what the error
is that you're calculating, it seems like the normalized residual or
something. This is a linear equation so I'm not entirely sure why the
non-linear residual wouldn't be very close to zero if the linear
solver is converging. Probably best just to see the code rather than
speculate.

Cheers,

Daniel

On Mon, Aug 1, 2016 at 7:50 AM, Nils Becker
 wrote:
> hi,
>
> i scanned the list and the documentation but have not been able to make
> this work properly yet:
>
> i have a 2D pde of the form
>
> \[
>
> \partial_t \rho =  \nabla \cdot (\rho \nabla \phi) + D \nabla^2 \rho +
> \lambda \rho
>
> \]
>
> where phi is a potential field and lambda is a source field, both are
> constant in time and smoothly varying in space.
>
> i have been solving this with fipy, using
>
> DiffusionTermCorrection(coeff=[D,]) for the diffusion
> PowerLawConvectionTerm(coeff=v ) where v = -phi.grad for the advection
> ImplicitSourceTerm(coeff=lambda) for the source
>
> the boundary conditions are either no-flux or dirichlet. i think
> advection is dominant, although i can't say i have actually tried to
> quantify that.
>
> i am testing the solution by subtracting the rhs and lhs of the equation
> and comparing that to \partial_t \rho in magnitude. i am having trouble
> making this error smaller than around 10%: at points where the density
> grows strongly, the solution deviates.
>
> at the moment i am using a loop like this:
>
> for (st, stnext) in zip(Solve.sample_t, Solve.sample_t[1:]):
> for tt in nx.arange(st, stnext, Solve.step):
> rho.updateOld()
> res = 1.
> while res > Solve.tolerance:
> res = eq.sweep(var=rho, dt=Solve.step)
>
> i've played around with making the mesh smaller -- did not seem to make
> any difference; reducing the tolerance -- some effect, and reducing the
> temporal step size, which did not help much.
>
> should i be doing something differently?
>
> thanks for any help!
>
> n.
> ___
> fipy mailing list
> fipy@nist.gov
> http://www.ctcms.nist.gov/fipy
>   [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]



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

reaction-advection-diffusion; problems making the solution converge

[Nils Becker]

hi,

i scanned the list and the documentation but have not been able to make
this work properly yet:

i have a 2D pde of the form

\[

\partial_t \rho =  \nabla \cdot (\rho \nabla \phi) + D \nabla^2 \rho +
\lambda \rho

\]

where phi is a potential field and lambda is a source field, both are
constant in time and smoothly varying in space.

i have been solving this with fipy, using

DiffusionTermCorrection(coeff=[D,]) for the diffusion
PowerLawConvectionTerm(coeff=v ) where v = -phi.grad for the advection
ImplicitSourceTerm(coeff=lambda) for the source

the boundary conditions are either no-flux or dirichlet. i think
advection is dominant, although i can't say i have actually tried to
quantify that.

i am testing the solution by subtracting the rhs and lhs of the equation
and comparing that to \partial_t \rho in magnitude. i am having trouble
making this error smaller than around 10%: at points where the density
grows strongly, the solution deviates.

at the moment i am using a loop like this:

for (st, stnext) in zip(Solve.sample_t, Solve.sample_t[1:]):
for tt in nx.arange(st, stnext, Solve.step):
rho.updateOld()
res = 1.
while res > Solve.tolerance:
res = eq.sweep(var=rho, dt=Solve.step)

i've played around with making the mesh smaller -- did not seem to make
any difference; reducing the tolerance -- some effect, and reducing the
temporal step size, which did not help much.

should i be doing something differently?

thanks for any help!

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