projection method 真的比piso快吗？

mohui

这两天看了一篇文献关于使用四阶龙格库塔时间推进代替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:

mohui

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

东岳

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

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

p残差很高，精度不高，

高到多少，贴log看看，

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

mohui

@李东岳，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%)

mohui

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

wwzhao

@mohui 感谢分享！

试了一下cavity算例，可以正常跑完。
将system/fvSolution里面的压力求解器的relTol设为0就可以使残差降下来，不过相应的矩阵求解迭代步要增加。

mohui

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

wwzhao

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

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

mohui

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

Foamer24

Zephyr

我没看明白啊，他这不还是解了poisson方程吗，那不可能快的起来啊

东岳

@mohui 在 projection method 真的比piso快吗？中说：

这两天看了一篇文献关于使用四阶龙格库塔时间推进代替PISO算法，他的结果显示要比piso快50%,而且精度还高。

只是一篇SCI，后续没准还有一篇SCI出现相反的断论呢，都不好说

mohui

@东岳东岳老师挖了一手好坟啊

Foamer24

这个已经有开源的的openfoam代码了，但是是openfoam2.1.1的，不知道有没有比较会的能做版本适应

mohui

@Foamer24 这个你所说的开源是哪篇文献的？

Foamer24

@mohui 看到rhoEnergyFoam中有三阶四步的龙格库塔迭代。但是OpenFOAM大网格量并行速度相比fluent，不是很快

mohui

@Foamer24 有文献吗？另外这个求解器是官方的吗？我很久没关注新的版本的of了

Foamer24

@mohui 不是，OpenFOAM大网格量计算效率真的比不过fluent，好奇怪

Foamer24

@mohui 终于帮你找到啦！A LOW-DISSIPATIVE SOLVER FOR TURBULENT COMPRESSIBLE FLOWS ON UNSTRUCTURED MESHES, WITH OPENFOAM IMPLEMENTATION