Thanks a lot Daniel. I updated the script with more realistic values for
the Buckley Leverett equation. In fact it is not strange that fipy could
not sense the left boundary. We fixed saturation (sw) on the left side
of the domain, i.e., a face value, While the velocity value in the
convection term is calculated using the cell value, which is equal to
zero every where except at the left boundary. I fixed the value of a
cell next to the left boundary and it works now, even without the
diffusion term. We might as well have used the sw.faceValue as the
coefficient of the the convection term, which works but I don't like its
result.
Here's the updated script for Buckley Leverett equation:
from fipy import *
u = 1.e-3
L = 100.
nx = 200
dt = 100.
dx = L/nx
mesh = Grid1D(dx = L/nx, nx = nx)
x = mesh.cellCenters
sw = CellVariable(mesh=mesh, name="saturation", hasOld=True, value = 0.)
sw.setValue(1,where = x<=dx)
sw.constrain(1.,mesh.facesLeft)
#sw.constrain(0., mesh.facesRight)
sw.faceGrad.constrain([0], mesh.facesRight)
eq = TransientTerm(coeff=1) + UpwindConvectionTerm(coeff = u
*(sw**3.)/(sw**4.+(1-sw)**2.*(1-sw**2.)) * [[1]]) == 0
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-5:
swres = eq.sweep(dt = dt, var = sw)
print(swres)
viewer.plot()
On Fri, 2012-09-14 at 11:12 -0400, Daniel Wheeler wrote:
> 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 ]
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
