Re: reaction-advection-diffusion; problems making the solution converge
On Fri, Aug 5, 2016 at 12:15 PM, Nils Becker wrote: > > ok. it seems that you have not changed the iteration, only added print > functions. so my code seems to be converging, yes? Absolutely. Appears to be converging well. > 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). It could converge well to the wrong solution. Convergence and accuracy are different. > 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. I'm not sure that comparing the residual across time steps is indicative of anything useful. For example, the residual for a large and small time step could be entirely different, but converge equally well and also be just as accurate in terms of the actually variable values. It's more the amount of convergence at a given time step that's useful (i.e. the residual after normalization by the first sweep). > so what sweep returns is basically A(n) * x(n-1) - b(n), yes? Yup. > 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. I see. The version of the code that I got from you only ever did one sweep. > 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? No they shouldn't require any sweeps. It's probably that the LU solver has a slightly more accurate starting condition and thus the normalization of the error is different. The matrix and b vector aren't changing. See, https://github.com/usnistgov/fipy/blob/develop/fipy/solvers/pysparse/linearLUSolver.py#L80 -- 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
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
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
>> 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
> 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
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
> 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 fp.Im
Re: reaction-advection-diffusion; problems making the solution converge
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
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 ]
