Re: 2D PDE

Adam Stone Fri, 07 Dec 2012 09:21:41 -0800

Altan,

There seem to be some problems with the final loop's indentation:

for i in range(10):
     # move forward in time by one time step
     psi.updateOld()

# sweep until you get the desired residue
     res = 100.
     while res > 0.001:
         res = eqn.sweep(dt=dt, solver=GeneralSolver(iterations=20, tolerance=1e-5))
    print res
         vi1.plot()

It seems like the plot() call has an extra indent that it shouldn't which for some reason is being treated as if it's part of the sweep loop rather than just causing an error. It also shouldn't be printing more than 10 res values unless 'print res' is part of the sweep loop. If I do this:

    while res > 0.001:
        res = eqn.sweep(dt=dt, solver=GeneralSolver(iterations=20, tolerance=1e-5))
        vi1.plot()
    print i, res

the plot shows the same divergence happening at the corners, but only this is printed:

0 0.000944978555
1 0.000885568500159

So your original code is actually behaving like this:

    while res > 0.001:
        res = eqn.sweep(dt=dt, solver=GeneralSolver(iterations=20, tolerance=1e-5))
        print i, res
        vi1.plot()

I don't quite understand why, but in any case it seems like what you want is:

    while res > 0.001:
        res = eqn.sweep(dt=dt, solver=GeneralSolver(iterations=20, tolerance=1e-5))
    print i, res
    vi1.plot()

This no longer shows the divergent behavior in the plot, but doesn't address the bigger problem that it's diverging on the third timestep. I believe this is related to your mesh. I don't understand what exactly the + [[-xmax],[-xmax]] term is doing, because if I do 'print mesh' I get 'UniformGrid2D(dx=0.2, nx=45, dy=0.2, ny=45)' whether or not this last term is commented out. I infer that this syntax is meant to mirror the +XY quadrant into the other 3 quadrants without having to solve redundant cells in a 4x larger grid. But if I run the code with this term commented out, the 10 steps are able to complete without divergence, so somehow this term is causing the solution to diverge.

The image in the plot doesn't appear to change after 10 steps, but if I go to 100 steps, I get this:

So something is indeed happening now, and without divergence, but it's not the expected behavior. Maybe someone else can see where to go next with this.

Adam

On 12/06/2012 09:42 PM, Allawala, Altan wrote:

Hello,

I'm trying to numerically solve for the steady-state solution of the following PDE:

\frac{\partial \psi}{\partial t}=-y\frac{\partial \psi}{\partial x}-2xy(y^2+1)\psi

It can be seen that the analytical steady state solution is:

\psi=exp(-(x^2+1)(y^2+1))

Attached is the program written in FiPy to do this. I've set the initial function to be a Gaussian and evolved it forward in time. For the first 10 seconds of running the code, very little happens. Then the top left and lower right portions of the grid develop non-zero values and the solution diverges.

I'm quite sure that I'm expressing this problem in the right way (ie: using a convection term) and that the residues and tolerances etc have been set sufficiently small. Any idea why the solution is diverging instead of converging to the analytic value?

Thanks for you help.

Regards,
Altan
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

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

Re: 2D PDE

Reply via email to