This seems to work okay:
from fipy import *
u = 1.e-3
L = 100.
nx = 200
dt = 1000.e-1
mesh = Grid1D(dx = L/nx, nx = nx)
sw = CellVariable(mesh=mesh, name="saturation", hasOld=True, value = 0.)
sw.constrain(1.,mesh.facesLeft)
#sw.constrain(0., mesh.facesRight)
sw.faceGrad.constrain([0], mesh.facesRight)
eq = TransientTerm(coeff=1) + UpwindConvectionTerm(coeff = 0.5 *
sw * [[1]]) == DiffusionTerm(1e-10)
sw.constrain(1.,mesh.facesLeft)
#sw.constrain(0., mesh.facesRight)
sw.faceGrad.constrain([0], mesh.facesRight)
steps = 500
viewer = Viewer(vars = sw, datamax=1.1, datamin=-0.1)
for step in range(steps):
sw.updateOld()
swres = 1.0e6
while swres > 1e-2:
swres = eq.sweep(dt = dt, var = sw)
print(swres)
viewer.plot()
I had to add the diffustion term so the equation could "feel" the
boundary condition. It is sort of crazy, but an upwind convection term
with a zero velocity will not interact with the boundary (I think, I
didn't think it through carefully). It sort of isn't a boundary value
problem at that point. It's confusing.
Also I used $\frac{\partial S^2}{\partial x}$ in place of
$S\frac{\partial S}{\partial x}$.
Hope that helps a bit.
On Thu, Sep 13, 2012 at 4:08 PM, A. A. Eftekhari <[email protected]> wrote:
> Hello all,
>
> I'm using fipy to solve flow in porous media problems, and I must say
> with the new version (3.0) there is no problem in solving single phase
> multi-component convection diffusion equation coupled with Darcy's law.
> My difficulty is with solving two-phase flow and specially
> Buckley-Leverett problem. For the sake of simplicity, let's assume that
> we can write it as
>
> \frac{\partial S}{\partial t}+S\frac{\partial S}{\partial x}=0
>
> Fipy does not accept multiplication of anything other than integers and
> floats to a ConvectionTerm and it issues an error when I use the
> following code. Is there a way to solve this problem?
>
> Thanks a lot in advance,
>
> Ehsan
>
> from fipy import *
>
> u = 1.e-3
> L = 100.
> nx = 200
> dt = 1000.
>
> mesh = Grid1D(dx = L/nx, nx = nx)
>
> sw = CellVariable(mesh=mesh, name="saturation", hasOld=True, value = 0.)
> sw.constrain(1.,mesh.facesLeft)
> #sw.constrain(0., mesh.facesRight)
> sw.faceGrad.constrain([0], mesh.facesRight)
>
> eq = TransientTerm(coeff=1) + sw*UpwindConvectionTerm(coeff = (u,)) ==
> 0.
>
> sw.constrain(1.,mesh.facesLeft)
> #sw.constrain(0., mesh.facesRight)
> sw.faceGrad.constrain([0], mesh.facesRight)
>
> steps = 500
> viewer = Viewer(vars = sw)
> for step in range(steps):
> sw.updateOld()
> swres = 1.0e6
> while swres > 1e-6:
> swres = eq.sweep(dt = dt, var = sw)
> print(swres)
> viewer.plot()
>
> _______________________________________________
> fipy mailing list
> [email protected]
> http://www.ctcms.nist.gov/fipy
> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
--
Daniel Wheeler
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
