Hi Altan, See the attached script. It probably does something closer to what you're after. I only solve it in the top left quadrant as the source term changes sign in the upper left and lower right quadrants which will cause issues (FiPy should handle this, but doesn't seem to be). Also I changed from the central difference to the upwind convection term. I needed to add a constraint to the left hand faces due to FiPy's weirdness with inlet and outlet boundary conditions.
The problem doesn't seem like a very good boundary value problem as it is just like solving a whole bunch of 1D problems. The numerical solution will eventually just go to zero everywhere as that is actually a solution to the steady state problem. Cheers On Fri, Dec 7, 2012 at 6:21 PM, Allawala, Altan <[email protected]>wrote: > Hi all, > > Adam, thanks for your efforts. Yes, adding the "+ [[-xmax],[-xmax]]" > after defining the mesh shifts the mesh by -xmax in the x and y direction. > Or at least this is my understanding. I was hoping that I was wrong so that > I could make a correction to my code but at this point I don't see why the > solution is diverging. > > Regarding the indentation, good spot. It was happening because the > indentation was a mixture of spaces and a tab. But as you said, this > unfortunately is not the problem. > > I find it interesting that the solution does not diverge when solved in > the 1st/3rd quadrant. However, the solution blows up in quite an odd way > when solved in the 2nd/4th quadrant. And it blows up in exactly the same > way as when the solution is solved in the 1st/3rd quadrant but running time > in reverse. I can see why this would happen from the signs of the terms in > the equation. But I'm not sure how to get around this. > > Regards, > Altan > > > On Fri, Dec 7, 2012 at 2:28 PM, Adam Stone <[email protected]> wrote: > >> 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 [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 ] >> >> > > _______________________________________________ > fipy mailing list > [email protected] > http://www.ctcms.nist.gov/fipy > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] > > -- Daniel Wheeler
<<image/png>>
test3.py
Description: Binary data
_______________________________________________ fipy mailing list [email protected] http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
