projection method 真的比piso快吗?



  • 这两天看了一篇文献关于使用四阶龙格库塔时间推进代替PISO算法,他的结果显示要比piso快50%,而且精度还高。因为他是有公布代码,所以修改起来很容易(我不太确定改的对不对)。这是我的修改的代码:

    #include "CreatePoissonMatrix.H"
        Info<< "\nStarting time loop\n" << endl;
    
        while (runTime.run())
        {
            #include "readTimeControls.H"
            #include "CourantNo.H"
            #include "setDeltaT.H"
            runTime++;
    
            Info<< "Time = " << runTime.timeName() << nl << endl;
            Uold =U; Uc=U; 
    
            phi =(fvc::interpolate(U) & mesh.Sf());
            dU= runTime.deltaT()*(fvc::laplacian(turbulence->nuEff(),U)-fvc::div(phi,U));
            Uc =Uc +a1*dU; U=Uold+0.5*dU;
            #include "PressureCorrection.H"
    
            phi =(fvc::interpolate(U) & mesh.Sf());
            dU= runTime.deltaT()*(fvc::laplacian(turbulence->nuEff(),U)-fvc::div(phi,U));
            Uc =Uc +a2*dU; U=Uold+0.5*dU;
            #include "PressureCorrection.H"
    
            phi =(fvc::interpolate(U) & mesh.Sf());
            dU= runTime.deltaT()*(fvc::laplacian(turbulence->nuEff(),U)-fvc::div(phi,U));
            Uc =Uc +a3*dU; U=Uold+dU;
            #include "PressureCorrection.H"
    
             phi =(fvc::interpolate(U) & mesh.Sf());
            dU= runTime.deltaT()*(fvc::laplacian(turbulence->nuEff(),U)-fvc::div(phi,U));
            Uc =Uc +a4*dU; U=Uc;
            #include "PressureCorrection.H"
            
            phi =(fvc::interpolate(U) & mesh.Sf());
            laminarTransport.correct();
            turbulence->correct();
             #include "continuityErrs.H"
    
            runTime.write();
    
            Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
                << "  ClockTime = " << runTime.elapsedClockTime() << " s"
                << nl << endl;
        }
    

    其中:PressureCorrection..H

    U.correctBoundaryConditions();
    solve( pEqn == fvc::div(U)/runTime.deltaT());
    U = U -fvc::grad(p)*runTime.deltaT();
    U.correctBoundaryConditions();
    

    CreatePoissonMatrix.H为:

    fvScalarMatrix pEqn
    (
        fvm::laplacian(p) 
    );
    pEqn.setReference(pRefCell,pRefValue);
    

    目前有两个问题,一:我算pityDaily这个算例,对比过速度确实快一些,(前提是piso也是循环四次,但是piso一般两次就够了),但是这里有个很大的问题,p残差很高,精度不高,所以不知道是我改的求解器有问题还是这个算法本身不是很好?
    二:我本来想算cavity这个算例,结果一算就崩了,所以不是很确定我改的对不对。

    在此问下大神们有什么意见或者见解。:big_mouth:



  • @李东岳 ,@赵一铭 ,呼叫大神,:happy: 。


  • 网格教授 OpenFOAM教授 管理员

    他的结果显示要比piso快50%,而且精度还高。

    我建议对比3伟复杂算例,计算若干小时。自带的pitzDaily计算时间太短了。你试试看?找个3D算例。

    p残差很高,精度不高,

    高到多少,贴log看看,

    或者把你做的求解器压缩一下上传我运行运行。



  • @李东岳 ,log文件太大了,我直接贴出最后运行的部分,这个是piso的:

    Courant Number mean: 0.0727094 max: 0.318223
    smoothSolver:  Solving for Ux, Initial residual = 0.00233765, Final residual = 9.36774e-07, No Iterations 2
    smoothSolver:  Solving for Uy, Initial residual = 0.00221123, Final residual = 1.10021e-06, No Iterations 2
    GAMG:  Solving for p, Initial residual = 0.02495, Final residual = 0.00191592, No Iterations 16
    time step continuity errors : sum local = 1.98408e-07, global = -2.5675e-08, cumulative = -5.27605e-07
    GAMG:  Solving for p, Initial residual = 0.0114244, Final residual = 0.000956831, No Iterations 6
    time step continuity errors : sum local = 8.06665e-08, global = -1.4887e-08, cumulative = -5.42492e-07
    GAMG:  Solving for p, Initial residual = 0.000984314, Final residual = 9.17457e-05, No Iterations 12
    time step continuity errors : sum local = 7.76255e-09, global = -1.02027e-09, cumulative = -5.43512e-07
    GAMG:  Solving for p, Initial residual = 0.000126896, Final residual = 8.14763e-07, No Iterations 9
    time step continuity errors : sum local = 6.93864e-11, global = 1.58971e-11, cumulative = -5.43496e-07
    smoothSolver:  Solving for k, Initial residual = 0.00183417, Final residual = 6.75215e-07, No Iterations 2
    bounding k, min: 0 max: 0.61915 average: 0.0500166
    ExecutionTime = 880.24 s  ClockTime = 907 s
    

    这个是rk4projectionfoam:

    Courant Number mean: 0.0691815 max: 0.254882
    Time = 0.1
    
    GAMG:  Solving for p, Initial residual = 0.197725, Final residual = 0.014632, No Iterations 4
    GAMG:  Solving for p, Initial residual = 0.0293635, Final residual = 0.00242183, No Iterations 8
    GAMG:  Solving for p, Initial residual = 0.215681, Final residual = 0.0184253, No Iterations 3
    GAMG:  Solving for p, Initial residual = 0.0182562, Final residual = 0.00127419, No Iterations 8
    smoothSolver:  Solving for k, Initial residual = 0.00192867, Final residual = 7.26218e-07, No Iterations 2
    bounding k, min: 0 max: 1.08318 average: 0.0514211
    time step continuity errors : sum local = 1.74224e-05, global = -1.07657e-07, cumulative = -0.0014665
    ExecutionTime = 621.75 s  ClockTime = 640 s
    

    [0_1506409559511_rk4projectionfoam.rar](正在上传 100%)



  • @李东岳 ,我没有权限上传所以就给百度云链接 http://pan.baidu.com/s/1c8kFlO


  • OpenFOAM副教授

    @mohui 感谢分享!

    1. 试了一下cavity算例,可以正常跑完。

    2. 将system/fvSolution里面的压力求解器的relTol设为0就可以使残差降下来,不过相应的矩阵求解迭代步要增加。



  • @wwzhao ,恩,我试试看。另外,请问你对这个有什么见解呢?噱头大于实际效果?


  • OpenFOAM副教授

    @mohui 正如你所说的,piso一般循环两次即可收敛,也就是要求解两次压力泊松方程,而rk4则需要求解4次。实际上对于不可压缩流体,求解压力泊松方程所需的时间占总计算时间的70%-80%左右,所以我觉得这方法对于一般问题没太大意义。

    但那篇文章作者声称是这方法是low-dissipative的,我想对于一些特定问题应该还是有用的:cheeky:



  • @wwzhao, 恩,确实是。看来这个问题可以不用去深究,我本来打算移植到多相流求解器看下是否节省时间,从目前结果来看,貌似也没省多少时间。