使用周期性网格时,python并行和顺序运行在fipy中产生不同的结果

2021-01-08 53点热度 0人点赞 0条评论

我运行了以下代码来模拟二维网格中圆柱体周围的流动:

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)
.




然而,当我以并行模式运行它时(比如使用两个处理器)
档案看起来很奇怪。对于速度剖面,y=0.2和y=0.8之间的曲线与y=0和y=1之间顺序计算得出的曲线大致相似。不过,压力分布有很大不同。为单元格510打印的值为-3.163(压力)、1.209(xvelocity)和-0.044(yvelocity)
.




使用网格2d
而不是周期网格2底端
,我包含了额外的边界条件,我相信在这个例子中,这相当于使用周期网格。新的BCS是:

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)
.




当我使用周期网格时,我的边界条件是否未被指定?(我想这可以解释同一个问题的两种不同的解决方案)或者并行计算和周期网格由于某种原因而不和谐?

查看隐藏内容需要支付:¥1
查看

未经允许不得转载!使用周期性网格时,python并行和顺序运行在fipy中产生不同的结果

本文地址:https://ans.52learn.online/2173

ANS52LEARN

DO BEST