Hi James,
Thanks for this example. I think I understand what is happening. Let's
start with the code for solving the problem,
~~~~
import fipy
import numpy
mesh = fipy.Grid1D(nx=10)
var = fipy.CellVariable(mesh=mesh)
var[:] = 2 * numpy.random.random(mesh.numberOfCells)
var.constrain(1., where=mesh.facesRight)
var.faceGrad.constrain(0, where=mesh.facesLeft)
eqn = fipy.ConvectionTerm([-1]) == fipy.DiffusionTerm(1.0)
for _ in range(100):
eqn.sweep(var)
print var
~~~
Note that the "var.faceGrad.constrain" is required with convection
terms to make sure that there is an outlet flow. This is a weird
decision in many ways, but was made for backwards compatibility with
previous FiPy versions, but this is not the main issue here.
Notice that when the diffusion coefficient is small or large then this
system solves in one sweep or close to one sweep (set it to 1000.0 or
0.0001 for example). The number of sweeps required to solve the system
is at its worst when the diffiusion and convection coefficients are
close in magnitude. So what is going on? Basically, the system is not
fully implicit.
The left hand boundary condition is dependent on the left hand cell
value but gets added to the b vector in the case of non upwind
convection terms. In the case of fully upwind convection terms this
contribution is zero so as the diffusion coefficient goes to zero the
CovectionTerm becomes fully upwind (or if a UpwindConvectionTerm is
selected). As the diffusion coefficient gets large the contribution of
the convection term decreases and thus the contribution from the
explicit boundary condition is reduced.
The discretization for the left most cell for a central difference
scheme (we're using power law above, but it has some central
difference contribution at D=1 and u=-1) is
D * phi0 - (u / 2 + D) * phi1 = -u * phiB / 2
where phi0 is the value at the left most cell and phi1 is the cell to
the right. In FiPy, phiB (the boundary value of phiB) is not
calculated implicitly, but is explicitly set for the convection term
boundary condition. Basically phiB is the previous value of phi0 in
the FiPy code. That may not be the best way to deal with this. It may
be that it needs to be explicit for stability reasons.
A possible way to deal with this issue is to use terms for the outlet
condition that are upwind at the outflow, but power law everywhere
else.
Anyway, I hope that helps with understanding what's going on even if
the implementation is sub-optimal.
Cheers,
Daniel
On Thu, 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
>
>
> _______________________________________________
> fipy mailing list
> fipy@nist.gov
> http://www.ctcms.nist.gov/fipy
> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>
--
Daniel Wheeler
_______________________________________________
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
