how to establish the pressure correction equation of the stokes

[321zhangtao]

I want to solve the following set of stokes equations, and I have written an expression for the momentum equation. However, when using the simple algorithm, the required pressure correction equation does not seem to have found a reference,even in www.ctcms.nist.gov/fipy/documentation/. I hope you can help me.

    self.mu = 1e-3  # dynamic viscosity
    self.rho = 1e3  #fluid density
    self.dt = 1
    self.pressureRelaxation = 0.8
    self.velocityRelaxation = 0.5
    self.mesh = fp.Grid3D(nx=10, ny=20, nz=10,dx=0.01, dy=0.01, dz=0.01)
    self.ap = CellVariable(mesh=mesh, value=1.)
    self.coeff = 1./ ap.arithmeticFaceValue*mesh._faceAreas * mesh._cellDistances
    self.pressure = fp.CellVariable(mesh=self.mesh,name='pressure', value=0.0)
    self.pressureCorrection = fp.CellVariable(mesh=self.mesh)
    self.gravity = fp.CellVariable(mesh=self.mesh,name='gravity', elementshape=(3,),value=[0,9.81*10**(-6),0])
    self.porosity = fp.CellVariable(mesh=self.mesh,name='porosity', value=0.0, hasOld=True)       
    self.xVelocity = fp.CellVariable(mesh=mesh, name='X velocity', value=0.0, hasOld=True)
    self.yVelocity = fp.CellVariable(mesh=mesh, name='Y velocity', value=0.0, hasOld=True)
    self.zVelocity = fp.CellVariable(mesh=mesh, name='Z velocity', value=0.0, hasOld=True)
    self.face_velocity = fp.FaceVariable(mesh=mesh, rank=1)
    self.xVelocityEq = fP.TransientTerm(coeff=self.porosity*self.rho)\
                       -fP.DiffusionTerm(coeff=self.porosity*self.mu)\
                       +fP.ConvectionTerm(coeff=self.porosity*self.xVelocity.old)\
                       -self.pressure.grad.dot([1., 0., 0.])\
                       +self.rho*self.porosity*self.gravity.dot([1., 0., 0.]) == 0
    self.yVelocityEq = fP.TransientTerm(coeff=self.porosity*self.rho)\
                       -fP.DiffusionTerm(coeff=self.porosity*self.mu)\
                       +fP.ConvectionTerm(coeff=self.porosity*self.yVelocity.old)\
                       -self.pressure.grad.dot([0., 1., 0.])\
                       +self.rho*self.porosity*self.gravity.dot([0., 1., 0.]) == 0
    self.zVelocityEq = fP.TransientTerm(coeff=self.porosity*self.rho)\
                       -fP.DiffusionTerm(coeff=self.porosity*self.mu)\
                       +fP.ConvectionTerm(coeff=self.porosity*self.zVelocity.old)\
                       -self.pressure.grad.dot([0., 0., 1.])\
                       +self.rho*self.porosity*self.gravity.dot([0., 0., 1.]) == 0
    self.pressureCorrectionEq = DiffusionTerm(coeff=coeff) - self.face_velocity.divergence\
                                -fP.TransientTerm(self.porosity,var=self.face_velocity)\ # it haven‘t been complished’

[wd15]

Try and simplify the problem. Start with only 2D and neglect the convection terms for now. Use the Stokes flow example as a basis and adapt it to the problem that you'd like to solve.

[321zhangtao]

Thank you. I know the example and have learned it, and also the Linear Equations . However, the convection is necessary in this example, and it is what confused me, and cannot be handled in pressure correction equation by myself.

[321zhangtao]

Try and simplify the problem. Start with only 2D and neglect the convection terms for now. Use the Stokes flow example as a basis and adapt it to the problem that you'd like to solve.

I followed your advice and tried it like the stokes flow example, but when program run to:
xVelocityEq.cacheMatrix()
self.xres = self.xVelocityEq.sweep(var=xVelocity,underRelaxation=velocityRelaxatio)
it reports liek this:
AttributeError: 'NoneType' object has no attribute 'cacheMatrix'
AttributeError: 'NoneType' object has no attribute 'sweep'
What happened?

[guyer]

Please show the code you are running now. xVelocityEq is evidently assigned to None, but that does not happen in either the Stokes flow example or in your code above.

[321zhangtao]

Please show the code you are running now. xVelocityEq is evidently assigned to None, but that does not happen in either the Stokes flow example or in your code above.

        self.mu = 1e-3  # dynamic viscosity
        self.rho = 1e3  #fluid density
        self.dt = 1
        self.pressureRelaxation = 0.8
        self.velocityRelaxation = 0.5
        self.mesh = fp.Grid3D(nx=10, ny=20, nz=10,dx=0.01, dy=0.01, dz=0.01)
        self.ap = fp.CellVariable(mesh=self.mesh, value=1.)
        self.coeff = 1./ self.ap.arithmeticFaceValue*self.mesh._faceAreas * self.mesh._cellDistances
        self.pressure = fp.CellVariable(mesh=self.mesh,name='pressure')
        self.pressureCorrection = fp.CellVariable(mesh=self.mesh)
        self.volume = fp.CellVariable(mesh=self.mesh, value=self.mesh.cellVolumes, name='Volume')
        self.contrvolume=self.volume.arithmeticFaceValue
        self.force= fp.CellVariable(mesh=self.mesh,name='force', elementshape=(3,))
        self.gravity = fp.CellVariable(mesh=self.mesh,name='gravity', elementshape=(3,),value=[0,-6.938,-6.938])
        self.porosity = fp.CellVariable(mesh=self.mesh,name='porosity', value=0.0, hasOld=True)       
        self.xVelocity = fp.CellVariable(mesh=self.mesh, name='X velocity', value=0.0, hasOld=True)
        self.yVelocity = fp.CellVariable(mesh=self.mesh, name='Y velocity', value=0.0, hasOld=True)
        self.zVelocity = fp.CellVariable(mesh=self.mesh, name='Z velocity', value=0.0, hasOld=True)
        self.face_velocity = fp.FaceVariable(mesh=self.mesh,name='face_velocity', rank=1)
        self.xVelocityEq = fp.TransientTerm(coeff=self.porosity*self.rho)\
                           -fp.DiffusionTerm(coeff=self.porosity*self.mu)\
                           +fp.ConvectionTerm(coeff=(self.porosity*self.xVelocity.old).faceGrad)\
                           -self.porosity*self.pressure.grad.dot([1., 0., 0.])\
                           +self.rho*self.porosity*self.gravity.dot([1., 0., 0.]) +self.force.dot([1., 0., 0.])
        self.yVelocityEq = fp.TransientTerm(coeff=self.porosity*self.rho)\
                           -fp.DiffusionTerm(coeff=self.porosity*self.mu)\
                           +fp.ConvectionTerm(coeff=(self.porosity*self.yVelocity.old).faceGrad)\
                           -self.porosity*self.pressure.grad.dot([0., 1., 0.])\
                           +self.rho*self.porosity*self.gravity.dot([0., 1., 0.])+self.force.dot([0., 1., 0.]) 
        self.zVelocityEq = fp.TransientTerm(coeff=self.porosity*self.rho)\
                           -fp.DiffusionTerm(coeff=self.porosity*self.mu)\
                           +fp.ConvectionTerm(coeff=(self.porosity*self.yVelocity.old).faceGrad)\
                           -self.porosity*self.pressure.grad.dot([0., 0., 1.])\
                           +self.rho*self.porosity*self.gravity.dot([0., 0., 1.]) +self.force.dot([0., 0., 1.])== 0
        self.pressureCorrectionEq = fp.DiffusionTerm(coeff=self.coeff)\
                                    - self.face_velocity.divergence
        X, Y , Z= self.mesh.faceCenters
        self.pressureCorrection.constrain((0.,), self.mesh.facesTop & (Y < 0.01))
        sweeps=100
        for sweep in range(sweeps):
            self.xVelocity.updateOld()
            self.yVelocity.updateOld()
            self.zVelocity.updateOld()
         ## solve the Stokes equations to get starred values
            self.xVelocityEq.cacheMatrix()
            self.xres = self.xVelocityEq.sweep(var=xVelocity,underRelaxation=velocityRelaxatio)
            self.xmat = self.xVelocityEq.matrix
            self.yres = self.yVelocityEq.sweep(var=yVelocity,underRelaxation=velocityRelaxation)
            self.zres = self.zVelocityEq.sweep(var=zVelocity,underRelaxation=velocityRelaxation)
            ## update the ap coefficient from the matrix diagonal
            self.ap[:] = -numerix.asarray(self.xmat.takeDiagonal())   ## update the face velocities based on starred values with the Rhie-Chow correction.
            ## cell pressure gradient
            presgrad = self.pressure.grad
            ## face pressure gradient
            facepresgrad = _FaceGradVariable(self.pressure)
            self.face_velocity[0] = self.xVelocity.arithmeticFaceValue + self.contrvolume / self.ap.arithmeticFaceValue * \
               (self.presgrad[0].arithmeticFaceValue-facepresgrad[0])
            self.face_velocity[1] = self.yVelocity.arithmeticFaceValue + self.contrvolume / self.ap.arithmeticFaceValue * \
               (self.presgrad[1].arithmeticFaceValue-facepresgrad[1])
            self.face_velocity[2] = self.zVelocity.arithmeticFaceValue + self.contrvolume / self.ap.arithmeticFaceValue * \
               (self.presgrad[2].arithmeticFaceValue-facepresgrad[2])
           ## solve the pressure correction equation
            self.pressureCorrectionEq.cacheRHSvector()
           ## left bottom point must remain at pressure 0, so no correction
            self.pres = self.pressureCorrectionEq.sweep(var=pressureCorrection,cacheRHSvector=True)
            self.rhs = self.pressureCorrectionEq.RHSvector
            ## update the pressure using the corrected value
            self.pressure.setValue(self.pressure + self.pressureRelaxation * self.pressureCorrection )
            ## update the velocity using the corrected pressure
            self.xVelocity.setValue(self.xVelocity - self.pressureCorrection.grad[0] / self.ap * self.mesh.cellVolumes)
            self.yVelocity.setValue(self.yVelocity - self.pressureCorrection.grad[1] / self.ap * self.mesh.cellVolumes)
            self.zVelocity.setValue(self.zVelocity - self.pressureCorrection.grad[2] / self.ap * self.mesh.cellVolumes)

[guyer]

This code is incomplete. If I add a class definition at the top:

class something():
    def __init__(self):
        pass
    
    def do_something(self):

and then run

something().do_something()

I get

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [3], line 1
----> 1 something().do_something()

Cell In [2], line 51, in something.do_something(self)
     49 ## solve the Stokes equations to get starred values
     50    self.xVelocityEq.cacheMatrix()
---> 51    self.xres = self.xVelocityEq.sweep(var=xVelocity,underRelaxation=velocityRelaxatio)
     52    self.xmat = self.xVelocityEq.matrix
     53    self.yres = self.yVelocityEq.sweep(var=yVelocity,underRelaxation=velocityRelaxation)

NameError: name 'xVelocity' is not defined

self.xVelocityEq.cacheMatrix() does not raise and error, because self.xVelocityEq is not None.

Please show us what you're actually doing.

[321zhangtao]

This code is incomplete. If I add a class definition at the top:

class something():
    def __init__(self):
        pass
    
    def do_something(self):

and then run

something().do_something()

I get

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [3], line 1
----> 1 something().do_something()

Cell In [2], line 51, in something.do_something(self)
     49 ## solve the Stokes equations to get starred values
     50    self.xVelocityEq.cacheMatrix()
---> 51    self.xres = self.xVelocityEq.sweep(var=xVelocity,underRelaxation=velocityRelaxatio)
     52    self.xmat = self.xVelocityEq.matrix
     53    self.yres = self.yVelocityEq.sweep(var=yVelocity,underRelaxation=velocityRelaxation)

NameError: name 'xVelocity' is not defined

self.xVelocityEq.cacheMatrix() does not raise and error, because self.xVelocityEq is not None.

Please show us what you're actually doing.

This code is incomplete. If I add a class definition at the top:

class something():
    def __init__(self):
        pass
    
    def do_something(self):

and then run

something().do_something()

I get

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [3], line 1
----> 1 something().do_something()

Cell In [2], line 51, in something.do_something(self)
     49 ## solve the Stokes equations to get starred values
     50    self.xVelocityEq.cacheMatrix()
---> 51    self.xres = self.xVelocityEq.sweep(var=xVelocity,underRelaxation=velocityRelaxatio)
     52    self.xmat = self.xVelocityEq.matrix
     53    self.yres = self.yVelocityEq.sweep(var=yVelocity,underRelaxation=velocityRelaxation)

NameError: name 'xVelocity' is not defined

self.xVelocityEq.cacheMatrix() does not raise and error, because self.xVelocityEq is not None.

Please show us what you're actually doing.

I found your advice was right, the self.xVelocityEq was written wrong. The new question was dt in TransientTerm, how can i set in my codes?
The report was listed:
File "C:\Program Files\Itasca\PFC500\exe64\python27\lib\site-packages\fipy\terms\transientTerm.py", line 156, in _checkDt
raise TypeError, "dt must be specified."

[guyer]

When you have a TransientTerm, the time step dt is a required argument to sweep() or solve(), e.g.,

self.xres = self.xVelocityEq.sweep(dt=<...>, var=xVelocity,underRelaxation=velocityRelaxatio)

Please work through some of the the introductory FiPy examples, especially examples.diffusion.mesh1D.

[321zhangtao]

Thanks for your kind advice, I'll try it!

[321zhangtao]

When you have a TransientTerm, the time step dt is a required argument to sweep() or solve(), e.g.,

self.xres = self.xVelocityEq.sweep(dt=<...>, var=xVelocity,underRelaxation=velocityRelaxatio)

Please work through some of the the introductory FiPy examples, especially examples.diffusion.mesh1D.

Hi guyer
I tryied it in Python before combined with Praticle Flow Code，the code is as followed:

import numpy as np
import fipy as fp
from fipy.tools import numerix
mu = 1e-3  # dynamic viscosity
rho = 1e3  #fluid density
dt = 0.1
total_t=10
pressureRelaxation = 0.8
velocityRelaxation = 0.5
mesh = fp.Grid3D(nx=10, ny=20, nz=10,dx=0.01, dy=0.01, dz=0.01)
ap = fp.CellVariable(mesh=mesh, value=1.)
time = fp.Variable()
coeff = 1./ ap.arithmeticFaceValue*mesh._faceAreas * mesh._cellDistances
pressure = fp.CellVariable(mesh=mesh,name='pressure')
pressureCorrection = fp.CellVariable(mesh=mesh)
volume = fp.CellVariable(mesh=mesh, value=mesh.cellVolumes, name='Volume')
contrvolume=volume.arithmeticFaceValue
force= fp.CellVariable(mesh=mesh,name='force', elementshape=(3,))
gravity = fp.CellVariable(mesh=mesh,name='gravity', elementshape=(3,),value=[0.,-9.81,0.])
porosity = fp.CellVariable(mesh=mesh,name='porosity', value=0.,  hasOld=True)       
xVelocity = fp.CellVariable(mesh=mesh, name='X velocity', value=0.,  hasOld=True)
yVelocity = fp.CellVariable(mesh=mesh, name='Y velocity',  value=0., hasOld=True)
zVelocity = fp.CellVariable(mesh=mesh, name='Z velocity', value=0., hasOld=True)
face_velocity = fp.FaceVariable(mesh=mesh,name='face_velocity', rank=2)
Uin1 = 2e-4
xVelocity.constrain((0.,),mesh.facesTop)
xVelocity.constrain((0.,),mesh.facesBottom)
yVelocity.constrain((-Uin1,),mesh.facesTop)
yVelocity.constrain((-Uin1,),mesh.facesBottom)
pressure.faceGrad.constrain([0.,0.,0.],mesh.facesBack)
steps=100
X, Y , Z= mesh.faceCenters
pressureCorrection.constrain((0.,), mesh.facesTop & (Y < 0.2))
t=0 
while t < total_t:
    t +=dt 
    xVelocity.updateOld()
    yVelocity.updateOld()
    zVelocity.updateOld()
    porosity.updateOld()
    xVelocityEq = fp.TransientTerm(coeff=rho*porosity)\
                      -fp.DiffusionTerm(coeff=porosity*mu)\
                      +fp.ConvectionTerm(coeff=(porosity*xVelocity.old).faceGrad)\
                      -porosity*pressure.grad.dot([1.,0.,0.])\
                      +rho*porosity*gravity.dot([1.,0.,0.]) 
    yVelocityEq = fp.TransientTerm(coeff=porosity*rho)\
                       -fp.DiffusionTerm(coeff=porosity*mu)\
                       +fp.ConvectionTerm(coeff=(porosity*yVelocity.old).faceGrad)\
                       -porosity*pressure.grad.dot([0.,1.,0.])\
                       +rho*porosity*gravity.dot([0.,1.,0.])
    zVelocityEq = fp.TransientTerm(coeff=porosity*rho)\
                     -fp.DiffusionTerm(coeff=porosity*mu)\
                     +fp.ConvectionTerm(coeff=(porosity*yVelocity.old).faceGrad)\
                     -porosity*pressure.grad.dot([0.,0.,1.])\
                      +rho*porosity*gravity.dot([0.,0.,1.]) 
    pressureCorrectionEq = fp.DiffusionTerm(coeff=coeff)\
                                - face_velocity.divergence
    for step in range(steps):
           xVelocityEq.cacheMatrix()
           xres = xVelocityEq.sweep(var=xVelocity,underRelaxation=velocityRelaxation,dt=dt)
           xmat = xVelocityEq.matrix
           yres = yVelocityEq.sweep(var=yVelocity,underRelaxation=velocityRelaxation,dt=dt)
           zres = zVelocityEq.sweep(var=zVelocity,underRelaxation=velocityRelaxation,dt=dt)
    ## update the ap coefficient from the matrix diagonal
           ap[:] = -numerix.asarray(xmat.takeDiagonal())   ## update the face velocities based on starred values with the Rhie-Chow correction.
    ## cell pressure gradient
           presgrad = pressure.grad
    ## face pressure gradient
           facepresgrad = fp._FaceGradVariable(pressure)
           face_velocity[0] = xVelocity.arithmeticFaceValue + contrvolume / ap.arithmeticFaceValue * \
                            (presgrad[0].arithmeticFaceValue-facepresgrad[0])
           face_velocity[1] = yVelocity.arithmeticFaceValue + contrvolume / ap.arithmeticFaceValue * \
                            (presgrad[1].arithmeticFaceValue-facepresgrad[1])
           face_velocity[2] = zVelocity.arithmeticFaceValue + contrvolume / ap.arithmeticFaceValue * \
                             (presgrad[2].arithmeticFaceValue-facepresgrad[2])
           pressureCorrectionEq.cacheRHSvector()
           pres = pressureCorrectionEq.sweep(var=pressureCorrection)
           rhs = pressureCorrectionEq.RHSvector
           pressure.setValue(pressure + pressureRelaxation * pressureCorrection )
           xVelocity.setValue(xVelocity - pressureCorrection.grad[0] / ap * mesh.cellVolumes)
           yVelocity.setValue(yVelocity - pressureCorrection.grad[1] / ap * mesh.cellVolumes)
           zVelocity.setValue(zVelocity - pressureCorrection.grad[2] / ap * mesh.cellVolumes)

there was an error, and I exmanied it for whole day and cannot figure out, can you help me ?

Traceback (most recent call last):

  File "H:\conda\lib\site-packages\spyder_kernels\py3compat.py", line 356, in compat_exec
    exec(code, globals, locals)

  File "e:\pfc\cfd\pychaxun.py", line 62, in <module>
    xres = xVelocityEq.sweep(var=xVelocity,underRelaxation=velocityRelaxation,dt=dt)

  File "H:\conda\lib\site-packages\fipy\terms\term.py", line 237, in sweep
    solver._solve()

  File "H:\conda\lib\site-packages\fipy\solvers\scipy\scipySolver.py", line 26, in _solve
    self.var[:] = numerix.reshape(self._solve_(self.matrix, self.var.ravel(), numerix.array(self.RHSvector)), self.var.shape)

  File "H:\conda\lib\site-packages\fipy\solvers\scipy\linearLUSolver.py", line 31, in _solve_
    LU = splu(L.matrix.asformat("csc"), diag_pivot_thresh=1.,

  File "H:\conda\lib\site-packages\scipy\sparse\linalg\_dsolve\linsolve.py", line 366, in splu
    return _superlu.gstrf(N, A.nnz, A.data, A.indices, A.indptr,

RuntimeError: Factor is exactly singular

[guyer]

This means your system of equations is not well posed. There could be a lot of reasons for that.

As @wd15 said:

Start with only 2D and neglect the convection terms for now

You've added transient behavior, convection, and a 3rd dimensions to our Stokes example. Add one at a time and get each working before you add the next.

[321zhangtao]

Hi， guyer
Your advice is valueble，and i have tried it these days. I firstly applied it in 3D, and then added some source term, the next is Convectionterm, it was all performed normally. But when i added xVelocity.old in Convectionterm, the residual turn out to be inf. Can you give me some tips? Is the Convectionterm described as ▽(εVV) handled properly? my code is followed:

from fipy import CellVariable, FaceVariable, Grid3D, DiffusionTerm, ConvectionTerm,Variable
from fipy.tools import numerix
from fipy.variables.faceGradVariable import _FaceGradVariable
L = 1.0
N = 10
dL = L / N
U = 1.
time = Variable()
dt=1
viscosity = 1  # dynamic viscosity
density = 1e3#fluid density
#0.8 for pressure and 0.5 for velocity are typical relaxation values for SIMPLE
pressureRelaxation = 0.8
velocityRelaxation = 0.5
if __name__ == '__main__':
    sweeps = 100
else:
    sweeps = 5
mesh = Grid3D(nx=N, ny=N, nz=N, dx=dL, dy=dL, dz=dL)
pressure = CellVariable(mesh=mesh, name='pressure')
pressureCorrection = CellVariable(mesh=mesh)
xVelocity = CellVariable(mesh=mesh, name='X velocity', value=0.,  hasOld=True)
yVelocity = CellVariable(mesh=mesh, name='Y velocity', value=0.,  hasOld=True)
zVelocity = CellVariable(mesh=mesh, name='z velocity', value=0.,  hasOld=True)
porosity = CellVariable(mesh=mesh,name='porosity', value=0.1,  hasOld=True)
force= CellVariable(mesh=mesh,name='force', elementshape=(3,),value=([0.1,0.1,0.1]),  hasOld=True)
gravity = CellVariable(mesh=mesh,name='gravity', elementshape=(3,),value=[0.,-9.81,0.])
velocity = FaceVariable(mesh=mesh,elementshape=(3,))
ap = CellVariable(mesh=mesh, value=1.)
coeff = 1./ ap.arithmeticFaceValue*mesh._faceAreas * mesh._cellDistances
volume = CellVariable(mesh=mesh, value=mesh.cellVolumes, name='Volume')
contrvolume=volume.arithmeticFaceValue
xVelocity.constrain(0., mesh.facesTop | mesh.facesBottom | mesh.facesFront | mesh.facesBack)
xVelocity.constrain(U, mesh.facesLeft)
xVelocity.constrain(U, mesh.facesRight)
yVelocity.constrain(0., mesh.exteriorFaces)
zVelocity.constrain(0., mesh.exteriorFaces)
X, Y, Z= mesh.faceCenters
pressureCorrection.constrain(0., mesh.facesBottom & (Y < dL))
#pressureCorrection.constrain(0., mesh.facesFront & (Z < dL))
while time() < 10:
      time.setValue(time() + dt)
          #porosity.setValue(0.5+time*0.02)
          #force.setValue(np.tile([0.1+time*0.1,0.32,0.2+time**2*0.01],(1000,1)).T)
      xVelocity.updateOld()
      yVelocity.updateOld()
      zVelocity.updateOld()
      porosity.updateOld() 
      xVelocityEq =-DiffusionTerm(coeff=viscosity*porosity)\
                 +ConvectionTerm(coeff=(density*porosity*xVelocity.old).faceGrad)\
                 -porosity*pressure.grad[0]\
                 +density*porosity*gravity[0]+force[0]
      yVelocityEq =-DiffusionTerm(coeff=porosity*viscosity)\
                 +ConvectionTerm(coeff=(density*porosity*yVelocity.old).faceGrad)\
                 -porosity*pressure.grad[1]\
                 +density*porosity*gravity[1]+force[1] 
      zVelocityEq =-DiffusionTerm(coeff=porosity*viscosity)\
                 +ConvectionTerm(coeff=(density*porosity*zVelocity.old).faceGrad)\
                 -porosity*pressure.grad[2]\
                 +density*porosity*gravity[2] +force[2]
      pressureCorrectionEq = DiffusionTerm(coeff=coeff)\
                                        - velocity.divergence
      for sweep in range(sweeps):
          xVelocityEq.cacheMatrix()
          xres = xVelocityEq.sweep(var=xVelocity,underRelaxation=velocityRelaxation)
          xmat = xVelocityEq.matrix
          yres = yVelocityEq.sweep(var=yVelocity,underRelaxation=velocityRelaxation)
          zres = zVelocityEq.sweep(var=zVelocity,underRelaxation=velocityRelaxation)   
    ## update the ap coefficient from the matrix diagonal
          ap[:] = -numerix.asarray(xmat.takeDiagonal())
          presgrad = pressure.grad
          facepresgrad = _FaceGradVariable(pressure)
          velocity[0] = xVelocity.arithmeticFaceValue + contrvolume / ap.arithmeticFaceValue * \
           (presgrad[0].arithmeticFaceValue-facepresgrad[0])
          velocity[1] = yVelocity.arithmeticFaceValue + contrvolume / ap.arithmeticFaceValue * \
           (presgrad[1].arithmeticFaceValue-facepresgrad[1])
          velocity[2] = zVelocity.arithmeticFaceValue + contrvolume / ap.arithmeticFaceValue * \
           (presgrad[2].arithmeticFaceValue-facepresgrad[2])
          velocity[..., mesh.exteriorFaces.value] = 0.
          velocity[0, mesh.facesLeft.value] = U
          velocity[0, mesh.facesRight.value] = U
    ## solve the pressure correction equation
          pressureCorrectionEq.cacheRHSvector()
    ## left bottom point must remain at pressure 0, so no correction
          pres = pressureCorrectionEq.sweep(var=pressureCorrection)
          rhs = pressureCorrectionEq.RHSvector
    ## update the pressure using the corrected value
          pressure.setValue(pressure + pressureRelaxation * pressureCorrection )
    ## update the velocity using the corrected pressure
          xVelocity.setValue(xVelocity - pressureCorrection.grad[0] / \
                                               ap * mesh.cellVolumes)
          yVelocity.setValue(yVelocity - pressureCorrection.grad[1] / \
                                               ap * mesh.cellVolumes)
          zVelocity.setValue(zVelocity - pressureCorrection.grad[2] / \
                                               ap * mesh.cellVolumes)
      if __name__ == '__main__':
          if time%2 == 0:
             print('sweep:', sweep, ', x residual:', xres, \
                                    ', y residual', yres, \
                                    ', z residual', zres, \
                                   ', p residual:', pres, \
                                   ', continuity:', max(abs(rhs)))