Closed 321zhangtao closed 6 months ago
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.
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.
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?
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.
Please show the code you are running now.
xVelocityEq
is evidently assigned toNone
, 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)
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, becauseself.xVelocityEq
is notNone
.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, becauseself.xVelocityEq
is notNone
.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."
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
.
Thanks for your kind advice, I'll try it!
When you have a
TransientTerm
, the time stepdt
is a required argument tosweep()
orsolve()
, 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
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.
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)))
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.