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 

- 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
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Reply via email to