Re: 2D PDE

[Adam Stone]

I just realized what the extra mesh term is doing, just shifting all the values by -xmax, so by removing the term only the first quadrant is solved, and this eliminates the instabilities because they don't occur in this quadrant. Not sure why that wasn't obvious this morning, but I see that it doesn't help resolve why there are instabilities in two of the other quadrants. I guess someone more skilled with fipy will have to answer that. Adam On 12/07/2012 12:21 PM, Adam Stone wrote: 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 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 ]