python 3
hi, i'm considering updating my scripts from 2.7 to python 3.4 or 3.6. should i expect any problems? is fipy tested on py 3? thanks! 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 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
>> 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
> 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
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 ]
