This is the entire script.
*********************************************************************************************************************************************************************************
from __future__ import division
porisity=1.
rho1=1.
rho2=1.
mu1=1.
mu2=1.
K=1
L=1.
nx=100
dx=L/100
valueLeft=0.
valueRight=0.
fluxLeft=0.
timeStepDuration=0.004
from fipy.meshes.grid1D import Grid1D
mesh=Grid1D(dx=dx, nx=nx)
from fipy.variables.cellVariable import CellVariable
phi=CellVariable(name="saturation", mesh=mesh)
p=CellVariable(name="pressure",mesh=mesh)
x=mesh.getCellCenters()[...,0]
phi.setValue(.9, where=x<.5)
phi.setValue(.3, where=x>=.5)
p.setValue(1.)
pcs=0
qin=40*((x>=.1)&(x<=.3))+20*((x>=.7)&(x<=.9))
qout=60*((x>=.4)&(x<=.6))
sin=.8
pin=0
timeStep=1/timeStepDuration
q1s=sin*qin-phi*qout
q2s=(1-sin)*qin-(1-phi)*qout
from fipy.boundaryConditions.fixedValue import FixedValue
from fipy.boundaryConditions.fixedFlux import FixedFlux
BCs1=(FixedValue(faces=mesh.getFacesRight(), value=valueRight),
FixedFlux(faces=mesh.getFacesLeft(), value=fluxLeft))
BCs2=(FixedValue(faces=mesh.getFacesRight(), value=0)),
if __name__=='__main__':
import fipy.viewers
phiviewer=fipy.viewers.make(vars=(phi), limits={'ymin':0.2,
'ymax':1,'datamin':0, 'datamax':1})
pviewer=fipy.viewers.make(vars=(p), limits={'datamin':-5, 'datamax':5})
phiviewer.plot()
pviewer.plot()
#raw_input("initial condition. press <return> to proceed...")
from fipy.terms.transientTerm import TransientTerm
from fipy.terms.powerLawConvectionTerm import PowerLawConvectionTerm
from fipy.terms.explicitSourceTerm import _ExplicitSourceTerm
from fipy.terms.implicitDiffusionTerm import ImplicitDiffusionTerm
eq1=TransientTerm(coeff=+1)==PowerLawConvectionTerm(coeff=1*p.getFaceGrad())+_ExplicitSourceTerm(q1s)
eq2=ImplicitDiffusionTerm(coeff=phi-1)==(phi-phi.getOld())/timeStepDuration+((1-phi)*(1-phi).getGrad()).getDivergence()+_ExplicitSourceTerm(q2s)#PowerLawConvectionTerm(coeff=2*p.getFaceGrad())-p.getFaceGrad().getDivergence()
steps=3
for i in range(steps):
eq1.solve(p, dt=timeStepDuration, boundaryConditions = BCs1)
eq2.solve(phi, dt=timeStepDuration, boundaryConditions=BCs2)
if __name__ =='__main__':
pviewer.plot()
phiviewer.plot()
if __name__ == '__main__':
raw_input("Explicit transient diffusion. Press <return> to proceed...")
*******************************************************************************************************************************************************************************************
2008/9/3 Daniel Wheeler <[EMAIL PROTECTED]>
>
> On Wed, Sep 3, 2008 at 7:07 AM, franck kalala
> <[EMAIL PROTECTED]> wrote:
> > In short
> >
> > I have an equation of this kind in my sets my equations
> >
> >
> > -\nabla\nabla p=nabla(D1(phi)\nabla D2(phi))+q1(phi)+q2(phi)
> >
> > D1(phi)=1-phi
> >
> > D2(phi)=1-phi
> >
> >
> > the unknown is p how do the second term of the above equation?
>
> The right hand side is just a source term. Write it in fipy as you
> would write it mathematically.
>
> Cheers
>
> --
> Daniel Wheeler
>
>
--
***********************************************
***********************************************
Franck Kalala Mutombo
[EMAIL PROTECTED]
African Institut for Mathematical Sciences, Muizenberg Cape Town, South
Africa
[EMAIL PROTECTED]
************************************************
************************************************
