Я запустил следующий код для имитации обтекания цилиндра в 2D-сетке:
from fipy import CellVariable, FaceVariable, Grid2D, DiffusionTerm, ImplicitSourceTerm, PeriodicGrid2DTopBottom, DistanceVariable, Viewer
from fipy.tools import numerix
L = 1.0
N = 50
dL = L / N
viscosity = 1
U = 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 = 500
else:
sweeps = 5
mesh = PeriodicGrid2DTopBottom(nx=N, ny=N, dx=dL, dy=dL)
pressure = CellVariable(mesh=mesh, name='pressure')
pressureCorrection = CellVariable(mesh=mesh)
xVelocity = CellVariable(mesh=mesh, name='X velocity')
yVelocity = CellVariable(mesh=mesh, name='Y velocity')
velocity = FaceVariable(mesh=mesh, rank=1)
pfi=3000.
lfi=0.01
x, y = mesh.cellCenters
var1 = DistanceVariable(name='distance to center', mesh=mesh, value=numerix.sqrt((x-N*dL/2.)**2+(y-N*dL/2.)**2))
rad=0.1
pi_fi= CellVariable(mesh=mesh, value=0.,name='Fluid-interface energy map')
pi_fi.setValue(pfi*numerix.exp(-1.*(var1-rad)/lfi), where=(var1 > rad) )
pi_fi.setValue(pfi, where=(var1 <= rad))
xVelocityEq = DiffusionTerm(coeff=viscosity) - pressure.grad.dot([1., 0.]) - ImplicitSourceTerm(pi_fi)
yVelocityEq = DiffusionTerm(coeff=viscosity) - pressure.grad.dot([0., 1.]) - ImplicitSourceTerm(pi_fi)
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(U, mesh.facesLeft | mesh.facesRight)
yVelocity.constrain(0., mesh.facesLeft | mesh.facesRight)
X, Y = mesh.faceCenters
pressureCorrection.constrain(0., mesh.facesLeft & (Y < dL))
if __name__ == '__main__':
viewer = Viewer(vars=(pressure, xVelocity, yVelocity, velocity),
xmin=0., xmax=1., ymin=0., ymax=1., colorbar=True)
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
yres = yVelocityEq.sweep(var=yVelocity,
underRelaxation=velocityRelaxation)
## 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[1] = yVelocity.arithmeticFaceValue \
+ contrvolume / ap.arithmeticFaceValue * \
(presgrad[1].arithmeticFaceValue-facepresgrad[1])
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)
if __name__ == '__main__':
if sweep%10 == 0:
print('sweep:', sweep, ', x residual:', xres, \
', y residual', yres, \
', p residual:', pres, \
', continuity:', max(abs(rhs)))
viewer.plot()
print(pressure.globalValue[..., 510])
print(xVelocity.globalValue[..., 510])
print(yVelocity.globalValue[..., 510])
Это должно решить уравнения Навье-Стокса с периодическими условиями сверху / снизу и цилиндром в центресетки (вот почему в моих уравнениях есть неявный исходный термин). Скорость равна единице слева и справа и параллельна оси X. Этот пример примерно соответствует потоку через множество цилиндрических препятствий, расположенных на одинаковом расстоянии. Когда я запускаю его без параллельных вычислений , профили давления и скорости выглядят нормально. Значения, напечатанные для ячейки 510: 2.788 (давление), 1.104 (xvelocity) и -0.289 (yvelocity) .
Однако при запуске в параллельном режиме (скажем, с использованием 2 процессоров), профили выглядят странно. Для профиля скорости график между y = 0.2 и y = 0.8 более или менее похож на график из последовательных вычислений между y = 0 и y = 1. Профиль давления совсем другой, хотя. Значения, которые печатаются для ячейки 510: -3,163 (давление), 1,209 (xvelocity) и -0,044 (yvelocity) .
Чтобы использовать Grid2D вместо PeriodicGrid2DTopBottom, я включил дополнительные граничные условия, которые, как я считаю, будут эквивалентны использованию периодической сетки в этом примере. Тогда новые BC будут:
xVelocity.constrain(U, mesh.facesLeft | mesh.facesRight)
xVelocity.faceGrad.constrain(mesh.faceNormals * 0., where=mesh.facesBottom | mesh.facesTop)
yVelocity.constrain(0., mesh.exteriorFaces)
Таким образом, я получаю одинаковые выходные данные, работая последовательно или параллельно с двумя процессорами: 2.784 (давление), 1.104 (xvelocity) и -0.290(yvelocity) .
Были ли мои граничные условия занижены, когда я использовал периодическую сетку? (Полагаю, это объясняет два разных решения одной и той же проблемы) Или параллельные вычисления и периодические сетки по какой-то причине не ладят?