I have, as suggested, attacked my 1d compressible flow using Stokes Cavity example in FiPy (I really want to use FiPy for that model). But it still doesn't work. During sweeping, I get the error "The coefficient must be rank 0 for a rank 0 solution variable.". I am very instered in know how to solve that, but also in understanging why this errror appears.
Code used:
from fipy import CellVariable, FaceVariable, Grid1D, DiffusionTerm, Viewer
from fipy.tools import numerix
import matplotlib.pyplot as plt
L = 1.0
N = 50
dL = L / N
viscosity = 1
#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 = 300
else:
sweeps = 5
mesh = Grid1D(nx=N, dx=dL)
pressure = CellVariable(mesh=mesh, name='pressure')
pressureCorrection = CellVariable(mesh=mesh)
xVelocity = CellVariable(mesh=mesh, name='X velocity')
velocity = FaceVariable(mesh=mesh, rank = 1)
xVelocityEq = DiffusionTerm(coeff=viscosity) - pressure.grad.dot([1.])
ap = CellVariable(mesh=mesh, value=1.)
coeff = 1./ ap.arithmeticFaceValue*mesh._faceAreas * mesh._cellDistances
pressureCorrectionEq = DiffusionTerm(coeff=coeff) - velocity.divergence
from fipy.variables.faceGradVariable import _FaceGradVariable
volume = CellVariable(mesh=mesh, value=mesh.cellVolumes, name='Volume')
contrvolume=volume.arithmeticFaceValue
xVelocity.constrain(10., mesh.facesRight | mesh.facesLeft )
X = mesh.faceCenters
pressureCorrection.constrain(0., mesh.facesLeft)
from builtins import range
for sweep in range(sweeps):
## solve the Stokes equations to get starred values
xVelocityEq.cacheMatrix()
xres = xVelocityEq.sweep(var=xVelocity,
underRelaxation=velocityRelaxation)
xmat = xVelocityEq.matrix
## 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 = _FaceGradVariable(pressure)
velocity[0] = xVelocity.arithmeticFaceValue \
+ contrvolume / ap.arithmeticFaceValue * \
(presgrad[0].arithmeticFaceValue-facepresgrad[0])
velocity[..., mesh.exteriorFaces.value] = 0.
## 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)
if __name__ == '__main__':
if sweep%10 == 0:
print('sweep:', sweep, ', x residual:', xres, \
', p residual:', pres, \
', continuity:', max(abs(rhs)))
viewer.plot()
