View Issue Details

IDProjectCategoryView StatusLast Update
0002028OpenFOAMBugpublic2016-03-20 20:53
ReporterLieven Assigned Tohenry  
Status closedResolutionno change required 
PlatformGNU/LinuxOSUbuntuOS Version14.04
Fixed in Versiondev 
Summary0002028: residualControl ignores U-residuals when using coupled PBiCCCG solver
DescriptionWhen attempting to stop a simulation automatically by setting a residualControl in fvSolution, the U-residuals are ignored when using the coupled PBiCCCG solver.
Steps To ReproduceUsing the incompressible/simpleFoam/pitzDaily tutorial, set fvSolution to:

        solver GAMG;
        tolerance 1e-06;
        relTol 0.1;
        smoother GaussSeidel;
        nPreSweeps 0;
        nPostSweeps 2;
        cacheAgglomeration on;
        agglomerator faceAreaPair;
        nCellsInCoarsestLevel 10;
        mergeLevels 1;

        type coupled;
        solver PBiCCCG;
        preconditioner DILU;
        tolerance (1e-8 1e-8 1e-8);
        relTol (0.001 0.001 0.001);
        solver smoothSolver;
        smoother symGaussSeidel;
        tolerance 1e-05;
        relTol 0.1;
        solver smoothSolver;
        smoother symGaussSeidel;
        tolerance 1e-05;
        relTol 0.1;

    nNonOrthogonalCorrectors 0;
    consistent yes;

        p 1e-1;
        U 1e-2;
        "(k|epsilon|omega|f|v2)" 1e-1;

Note that the residualControl-residuals are set to very tolerant values for illustration purpose.

Additional InformationWhen using the coupled-solver for U, the simulation "converges" in 7 iterations although the final iteration shows:


Time = 7

DILUPBiCCCG: Solving for Ux, Initial residual = 0.0458328, Final residual = 6.89602e-06, No Iterations 9
DILUPBiCCCG: Solving for Uy, Initial residual = 0.0604103, Final residual = 1.84103e-05, No Iterations 9
DILUPBiCCCG: Solving for Uz, Initial residual = 0.398129, Final residual = 0.000116036, No Iterations 9
GAMG: Solving for p, Initial residual = 0.0791176, Final residual = 0.00740836, No Iterations 4
time step continuity errors : sum local = 0.750902, global = 0.10451, cumulative = 0.979509
smoothSolver: Solving for epsilon, Initial residual = 0.0323437, Final residual = 0.00199505, No Iterations 3
smoothSolver: Solving for k, Initial residual = 0.0584042, Final residual = 0.00429584, No Iterations 4
ExecutionTime = 1.02 s ClockTime = 1 s

SIMPLE solution converged in 7 iterations

When using the smoothSolver settings, the solution converges in 223 iterations, taking also the U-residual into account:

Time = 223

smoothSolver: Solving for Ux, Initial residual = 0.00128387, Final residual = 0.000127913, No Iterations 5
smoothSolver: Solving for Uy, Initial residual = 0.0099375, Final residual = 0.000729481, No Iterations 6
GAMG: Solving for p, Initial residual = 0.00680682, Final residual = 0.000619236, No Iterations 6
time step continuity errors : sum local = 0.0299118, global = 0.00248546, cumulative = 1.20842
smoothSolver: Solving for epsilon, Initial residual = 0.000812368, Final residual = 5.19083e-05, No Iterations 3
smoothSolver: Solving for k, Initial residual = 0.00190921, Final residual = 0.00011103, No Iterations 4
ExecutionTime = 19.16 s ClockTime = 23 s

SIMPLE solution converged in 223 iterations

TagsNo tags attached.



2016-03-16 16:12

manager   ~0006039

This should already be fixed in OpenFOAM-dev, could you test it and report back?


2016-03-19 16:48

updater   ~0006042

I've tested with the latest OpenFOAM-dev and the problem is that the Z direction is being solved when it shouldn't :(
Namely, the initial residual for Uz direction doesn't go below 0.28, hence it never falls below 0.01.

This is very strange, I didn't look more closely to coupled matrix solvers a few months ago, as I had assumed this was already taken into account in those situations... give me a few minutes and I'll see if I can figure this one out.


2016-03-19 20:11

updater   ~0006043

After investigating the code and doing some testing, there are two possible solutions that I can think of and none of them are pleasant to implement:

  1- "normFactor" can be neutralized by setting to a high-enough value for the respective invalid directions/components. The problem here is that each "invalid direction" is still solved with just numerical noise (all values near 0.0, but rarely zero) and if something weird occurs in the matrix solver (maybe somehow a division by zero or near infinite!?) while solving this direction, there will be a crash. Because this and other disabled directions will still be "solved" for the same number of iterations as the other directions, which is what this coupled matrix-solving feature was conceived for.

  2- Having an automatic rescaling of number of components in the matrices to be solved. Although this is theoretically possible, it's somewhat of a nightmare situation to develop. This means having 2D vectors, 5 or more crippled tensors (without 1 or more components) and I don't even know what else.

Therefore, the most practical solution that I can think of is this: do not use single-equation-coupled matrix-solvers for 2D cases.

Although this reminds me of Henry's following question on issue #1824: - namely: «There are still some choices to be made, for example should the segregated approach be maintained?»
As far as I can figure out, this is a scenario where the segregated solver might be a more practical implementation for at least 2D simulations over a 3D volume, when solving the momentum equation.
... Although having at least a "vector2D" class would still be practical enough (it's just one more type)... it's the tensors that are a bigger concern... at least as far as I can deduce right now.


2016-03-19 22:23

reporter   ~0006044

Wyldckat, have you also tested this with a 3D case? I'm asking because, although I reported it using a 2D tutorial case, I initially encountered this problem when running a 3D case.


2016-03-19 23:25

updater   ~0006045

@Lieven: Sorry, I should have been clearer:

 1- I used OpenFOAM-dev: - which is the repository with the latest developments for OpenFOAM and which will be part of the next major release. It's the version Henry asked about if you could test with it. I was afraid you wouldn't test with it, so I went ahead and tried it myself.

 2- Using the following settings:

        p 1e-1;
        U 0.30;
        "(k|epsilon|omega|f|v2)" 1e-1;

  would result in simpleFoam stopping at Time = 9. But again, this was with OpenFOAM-dev.

And as I implied, on issue #1824 - - is a similar report to your report and that one was fixed in OpenFOAM-dev, namely on this commit:

This fix cannot be implemented in OpenFOAM 3.0.x, because it would break compatibility with 3.0.1 and 3.0.0, due to the extent of the changes that were needed to fix this issue.


2016-03-20 10:38

manager   ~0006047

This is already fixed in OpenFOAM-dev in which the convergence criterion may be specified for each component.


2016-03-20 20:53

updater   ~0006050

Sorry to reopen+close again, but for the sake of completion, here's the updated example for the residual control per U component:

        p 1e-1;
        Ux 1e-2;
        Uy 1e-2;
        "(k|epsilon|omega|f|v2)" 1e-1;

I had totally forgot about this feature :(

Issue History

Date Modified Username Field Change
2016-03-16 15:58 Lieven New Issue
2016-03-16 16:12 henry Note Added: 0006039
2016-03-19 16:48 wyldckat Note Added: 0006042
2016-03-19 20:11 wyldckat Note Added: 0006043
2016-03-19 22:23 Lieven Note Added: 0006044
2016-03-19 23:25 wyldckat Note Added: 0006045
2016-03-20 10:38 henry Note Added: 0006047
2016-03-20 10:38 henry Status new => closed
2016-03-20 10:38 henry Assigned To => henry
2016-03-20 10:38 henry Resolution open => no change required
2016-03-20 10:38 henry Fixed in Version => dev
2016-03-20 20:53 wyldckat Note Added: 0006050
2016-03-20 20:53 wyldckat Status closed => feedback
2016-03-20 20:53 wyldckat Resolution no change required => reopened
2016-03-20 20:53 wyldckat Status feedback => closed
2016-03-20 20:53 wyldckat Resolution reopened => no change required