James - I think Daniel will have more insight into why this is happening and if there is anything to be done about it besides artificial relaxation, but I just want to say how much I appreciate your putting this together. This is a very lucid illustration.
- Jon > On Sep 15, 2016, at 5:13 PM, James Pringle <jprin...@unh.edu> wrote: > > Dear FiPy users -- > > This is a simple example of how and why fipy may fail to solve a > steady advection diffusion problem, and how solving the transient > problem can reduce the error. I also found something that was a > surprise -- the "initial" condition of a steady problem can affect > the solution for some solvers. > > At the end are two interesting questions for those who want to > understand what FiPy is actually doing.... I admit to being a bit > lost > > The equation I am solving is > > \Del\dot (D\del psi + u*psi) =0 > > Where the diffusion D is 1, and u is a vector (1,0) -- so there is > only a flow of speed -1 in the x direction. I solve this equation > on a 10 by 10 grid. The boundary conditions are no normal gradient > on the y=0 and y=10 boundary: > > psi_y =0 at y=0 and y=10 > > For the x boundary, I impose a value of x=1 on the inflow boundary at x=10 > (this is a little tricky -- the way the equation is written, u is > the negative of velocity). > > psi=1 at x=10 > > and a no-normal-gradient condition at x=0. > > psi_x=0 at x=0 > > since all of the domain and boundary is symmetrical with respect to > y, we can assume psi_y=0 is zero everywhere. This reduces the equation to > > psi_xx + psi_x =0 > > The general solution to this equation is > > psi=C1+C2*exp(-x) > > Where C1 and C2 are constants. For these boundary conditions, C1=1 > and C2=0, so psi=1 everywhere. > > Now run the code SquareGrid_HomemadeDelaunay and look at figure(3) > -- this is the plot of psi versus x, and you can see that it does > not match the true solution of psi=1 everywhere! Instead, it > appears to be decaying exponential. The blue line is actually just > (1+exp(-x)). What is going on? It appears that small errors in the > boundary condition at x=0 are allowing C2 to not be exactly 0, and > this error is this exponential mode. The error is the artificial > exiting of a correct mode of the interior equation, albeit one that > should not be excited by these BC's. > > Interestingly, the size of this error depends on the value the psi > is initially set to. If the line > > psi=fipy.CellVariable(name='psi',mesh=mesh,value=0.0) > > is changed so psi is initially 1, the error goes away entirely; if > it is set to some other value, you get different errors. I do not > know if this is a bug, or just numerical error in a well programmed > solver. > > Now if you run SquareGrid_HomemadeDelaunay_transient which implements > > psi_t = \Del\dot (D\del psi + u*psi) > > you can see that the error in the numerical solution is advected > out of the x=0 boundary, and the solution approaches the true > solution of psi=1 rapidly. > > The interesting question is if the formulation of the boundary > condition at x=0 could be altered to less excite the spurious mode? > > Also, why does the "initial" condition have any effect on the > steady equation? > > Cheers, > Jamie > > <JMP_Make_Grids.py><SquareGrid_HomemadeDelaunay_transient.py><SquareGrid_HomemadeDelaunay.py>_______________________________________________ > 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 ]
