[joaran]

Dear CFD Community,



I'm using OpenFoam 2.2.0, trying to solve the XiFoam tutorial (moriyoshiHomogeneous). The only file I modified is the fvSolution, where I changed the nOuterCorrectors to 200 and nCorrectors to 2, and added the ResidualControl and relaxationFactors subdicts. The file is attached. The problem is that in many time steps, all the fields' residuals start oscillating around a certain value, i.e. their mean values "flatten", and the PIMPLE algorithm does not converge to the tolerances I set for them (being the pressure field the most affected one, with the highest initial residuals). It doesn't seem to matter whether I refine the mesh, or coarsen it, or even change the time step or BCs. See attached a residual plot for the first time steps of the run. I also tried adjusting the relaxation factors for all the fields, each at a time, and altogether, and perhaps it converges well where it didn't before, but at some point it will show this same behaviour again.
I've visited many threads where this topic is discussed for different solvers, and everywhere the cause of the problem is pointed towards the quality of the mesh and/or appropiateness of BCs.

I don't consider it to be applicable for this case I'm bringing, because the mesh is just composed of squares, and there are no inlets or outlets!

I also noticed that other solvers using PIMPLE carry this "issue" too. Say, e.g. in Tobi's book*, 4th ed., figure 11.19. There, the second time step shows the same I'm experiencing. I also tried with reactingFoam, and the tutorial case of this solver also shows such issue.


Perhaps this behaviour is unavoidable, and it's something we must learn to live with. I'm not sure I'm comfortable though, ending the simulation with many time steps having initial residuals on the order of 10^-4 or 10^-3, or sometimes even 10^-2 just because I don't have the skills or knowledge to lessen them.



Any word of advice or help will be very much highly appreciated, folks!




Best,
Joaquín





*Tobias Holzmann. Mathematics, Numerics, Derivations and OpenFOAM(R), Holz-
mann CFD, Leoben, fourth edition, February 2017. URL www.holzmann-cfd.de.

[wyldckat]

Greetings Joaquín,

I've only now managed to get to the PM you sent me on the 11th of February about this thread...

I'm not familiar with this tutorial case and the solver, but chemistry is a pretty complicated thing, specially if you're trying to use single precision.

So the very first question is: Are you using an SP or DP build? If you run this command:

echo $WM_PRECISION_OPTION

it should tell you which one you are using. You should be using DP.

Now, using "nOuterCorrectors" set to 200 is akin to trying to beat numbers into a very fine pulp and expect it to still make sense... mmm... how else can I explain this... you are discretizing time and space, as well as using discrete models. But if there are any possible inconsistencies in your boundary conditions or in the models, it is possible that instability arises from using the models in an unexpected way.

So, for example, by using "nOuterCorrectors" set to 200, you are telling the solver that "this mesh is perfect and I don't care what you think it is, it's 100% correct"... but the solver goes and does that and the result is that the time step solved right now, does not correlate well with the next time step.
For example, imagine that the domain was initialized as a perfect vacuum. Then you open an inlet at 1 bar... how small of a time step do you believe would be necessary for it to even reach 0.1s without crashing? If you use 1e-2 and 200 nOuterCorrectors, the changes between time step are so massive, that each time step is completely disjointed from the previous time step, that enforcing 100% accuracy in each time step would be nearly irrelevant, given the sudden change.

As for the OpenFOAM version you are using, I strongly suggest that you try with a more recent version. Things have changed a lot since then and one particular improvement comes to mind: PBiCGStab, a more stable implementation of PBiCG.


So my guess here is that several things might be at play here:

1. The original tutorial case might be incorrectly set-up or that it was only a proof of concept of some sort... namely that it worked just fine as it were, but any changes like this one will result in problems.
2. If SP is being used, accuracy goes out the door and never comes back... so 200 or 20 nOuterCorrectors would almost give you the same result...
3. Unfit for purpose, such as the old PBiCG or any other modelling details may be unfit for what you want to do.
4. It is possible that for so many nOuterCorrectors, the solver needs to be rewritten, so that it can properly handle solving each group of equations in stages, instead of all at once...

Either way, using 200 nOuterCorrectors usually doesn't make much sense, unless you are doing some special coupling mechanisms that require refining the result for the current time step. Because making a very large time step to take advantage of it, will require a fairly simple geometry.

Best regards,
Bruno

[joaran]

Hi Bruno, thanks for taking the time to look at this.

Originally Posted by wyldckat

So the very first question is: Are you using an SP or DP build? If you run this command:

echo $WM_PRECISION_OPTION

it should tell you which one you are using. You should be using DP.

I'm using DP.

Originally Posted by wyldckat

Now, using "nOuterCorrectors" set to 200 is akin to trying to beat numbers into a very fine pulp and expect it to still make sense... mmm... how else can I explain this... you are discretizing time and space, as well as using discrete models. But if there are any possible inconsistencies in your boundary conditions or in the models, it is possible that instability arises from using the models in an unexpected way.

I just set to 200 the "nOuterCorrectors" to tell the solver that I wanted it to end the time step only when the indicated initial residuals were below certain levels. I actually don't expect it to take so longer, and as you say 200 is a big number, when the time step converges it usually doesn't take more than 20 or 30 iterations, depending on URFs.

Originally Posted by wyldckat

As for the OpenFOAM version you are using, I strongly suggest that you try with a more recent version. Things have changed a lot since then and one particular improvement comes to mind: PBiCGStab, a more stable implementation of PBiCG.

I've installed the latest OF version (6) but this particular issue remains more or less the same.

Something to be noticed is that when I turn off the thermo.correct() method from the solver, it doens't have any trouble to converge (but the solution can't be trusted anymore because it stops updating the thermophysical properties with the temperature, and in an IC application those changes are substantial) so I think the problem may be linked somehow to an oscillation due to thermo properties. Anyway, I'll keep trying to tackle this problem and I'll post it here if I arrive to a solution.

Best,
Joaquín

[wyldckat]

Quick answer: Then my last suspects are:

1. Mesh: any flaws in the mesh can have drastic consequences. Just because simpleFoam seems to work with a mesh, doesn't meant that the fluid flowed properly in it.
   • Run a full check with:
     

     checkMesh -allTopology -allGeometry

2. Boundary conditions: Main suspect would be the pressure field.
   • Sometimes the case is initialized as if we suddenly turn on gravity or turn on velocity.
   • Try first converging the flow field. Then turn on thermodynamics and chemistry after the flow is converged.

[joaran]

Hi Bruno,


The mesh seems to be fine:

Checking topology...
    Boundary definition OK.
    Cell to face addressing OK.
    Point usage OK.
    Upper triangular ordering OK.
    Face vertices OK.
    Topological cell zip-up check OK.
    Face-face connectivity OK.
  <<Writing 4 cells with two non-boundary faces to set twoInternalFacesCells
    Number of regions: 1 (OK).

Checking patch topology for multiply connected surfaces...
    Patch               Faces    Points   Surface topology                   Bounding box
    left                35       72       ok (non-closed singly connected)   (0 0 0) (0 0.035 0.001)
    right               35       72       ok (non-closed singly connected)   (0.07 0 0) (0.07 0.035 0.001)
    top                 70       142      ok (non-closed singly connected)   (0 0.035 0) (0.07 0.035 0.001)
    bottom              70       142      ok (non-closed singly connected)   (0 0 0) (0.07 0 0.001)
    frontAndBack        4900     5112     ok (non-closed singly connected)   (0 0 0) (0.07 0.035 0.001)

Checking geometry...
    Overall domain bounding box (0 0 0) (0.07 0.035 0.001)
    Mesh has 2 geometric (non-empty/wedge) directions (1 1 0)
    Mesh has 2 solution (non-empty) directions (1 1 0)
    All edges aligned with or perpendicular to non-empty directions.
    Boundary openness (-2.98368e-18 5.63584e-18 2.38694e-16) OK.
    Max cell openness = 1.05879e-16 OK.
    Max aspect ratio = 1 OK.
    Minimum face area = 1e-06. Maximum face area = 1e-06.  Face area magnitudes OK.
    Min volume = 1e-09. Max volume = 1e-09.  Total volume = 2.45e-06.  Cell volumes OK.
    Mesh non-orthogonality Max: 0 average: 0
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 8.32667e-14 OK.
    Coupled point location match (average 0) OK.
    Face tets OK.
    Min/max edge length = 0.001 0.001 OK.
    All angles in faces OK.
    Face flatness (1 = flat, 0 = butterfly) : min = 1  average = 1
    All face flatness OK.
    Cell determinant (wellposedness) : minimum: 1 average: 3.8302
    Cell determinant check OK.
    Concave cell check OK.
    Face interpolation weight : minimum: 0.5 average: 0.5
    Face interpolation weight check OK.
    Face volume ratio : minimum: 1 average: 1
    Face volume ratio check OK.

Mesh OK.

BC's are symmetryPlane for every field, and one thing to notice is that it starts converging well, so it doesn't seem to be an initialization problem. I also have tried converging the flow and applying thermodynamics afterwards, but it doesn't converge anyway.


I appreciate your interest!


Best,
Joaquín

[wyldckat]

Quick answer: I've left this on my to-do list and finally managed to take another look into this... although I don't have a solution.

So the first problem is that this case is not entirely well defined... if everything is set to symmetry, it can get borderline undefined, in the sense that one small error in modelling can make the simulation loose it's north... for example, any accidental distortions to the pressure field could make the case evolve in an entirely different direction...

The residuals we see don't always tell the whole story, given that they are normalized results. If you edit the file "OpenFOAM-2.2.x/etc/controlDict" and change this:

    lduMatrix           1;

to this:

    lduMatrix           2;

you will get a few more details about each equation being solved, e.g.:

Courant Number mean: 0.00395639 max: 0.0164715
Time = 0.0051

diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
PIMPLE: iteration 1
--> FOAM Warning : 
    From function lduMatrix::operator-=(const lduMatrix& A)
    in file matrices/lduMatrix/lduMatrix/lduMatrixOperations.C at line 286
    Unknown matrix type combination
    this : diagonal:0 symmetric:0 asymmetric:1
    A    : diagonal:0 symmetric:0 asymmetric:0
   Normalisation factor = 0.409449
DILUPBiCG:  Solving for Ux, Initial residual = 0.00433255, Final residual = 3.98411e-08, No Iterations 1
   Normalisation factor = 0.328259
DILUPBiCG:  Solving for Uy, Initial residual = 0.00505369, Final residual = 4.55551e-08, No Iterations 1
Max St-Courant Number = 0.0401747
--> FOAM Warning : 
    From function lduMatrix::operator-=(const lduMatrix& A)
    in file matrices/lduMatrix/lduMatrix/lduMatrixOperations.C at line 286
    Unknown matrix type combination
    this : diagonal:0 symmetric:0 asymmetric:1
    A    : diagonal:0 symmetric:0 asymmetric:0
   Normalisation factor = 0.0664362
DILUPBiCG:  Solving for b, Initial residual = 0.000856083, Final residual = 2.67585e-08, No Iterations 1
min(b) = 0.000161812
   Normalisation factor = 0.116263
DILUPBiCG:  Solving for Xi, Initial residual = 0.000866927, Final residual = 1.42646e-08, No Iterations 1
max(Xi) = 2.94208
max(XiEq) = 5.07463
Combustion progress = 3.9669%
--> FOAM Warning : 
    From function lduMatrix::operator-=(const lduMatrix& A)
    in file matrices/lduMatrix/lduMatrix/lduMatrixOperations.C at line 286
    Unknown matrix type combination
    this : diagonal:0 symmetric:0 asymmetric:1
    A    : diagonal:0 symmetric:0 asymmetric:0
   Normalisation factor = 48.3544
DILUPBiCG:  Solving for hau, Initial residual = 0.152665, Final residual = 2.40677e-06, No Iterations 1
--> FOAM Warning : 
    From function lduMatrix::operator+=(const lduMatrix& A)
    in file matrices/lduMatrix/lduMatrix/lduMatrixOperations.C at line 207
    Unknown matrix type combination
    this : diagonal:0 symmetric:0 asymmetric:1
    A    : diagonal:0 symmetric:0 asymmetric:0
   Normalisation factor = 297.45
DILUPBiCG:  Solving for ha, Initial residual = 0.0257083, Final residual = 2.5201e-06, No Iterations 1
--> FOAM Warning : 
    From function lduMatrix::operator-=(const lduMatrix& A)
    in file matrices/lduMatrix/lduMatrix/lduMatrixOperations.C at line 286
    Unknown matrix type combination
    this : diagonal:0 symmetric:1 asymmetric:0
    A    : diagonal:0 symmetric:0 asymmetric:0
   Normalisation factor = 0.000134327
DICPCG:  Solving for p, Initial residual = 0.694001, Final residual = 0.0202307, No Iterations 2
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 3.97874e-06, global = 1.00026e-06, cumulative = 2.29312e-05
PIMPLE: iteration 2
--> FOAM Warning : 
    From function lduMatrix::operator-=(const lduMatrix& A)
    in file matrices/lduMatrix/lduMatrix/lduMatrixOperations.C at line 286
    Unknown matrix type combination
    this : diagonal:0 symmetric:0 asymmetric:1
    A    : diagonal:0 symmetric:0 asymmetric:0
   Normalisation factor = 0.410598
DILUPBiCG:  Solving for Ux, Initial residual = 2.5932e-05, Final residual = 3.18415e-10, No Iterations 1
   Normalisation factor = 0.32927
DILUPBiCG:  Solving for Uy, Initial residual = 3.37783e-05, Final residual = 4.70605e-10, No Iterations 1
Max St-Courant Number = 0.0400679
--> FOAM Warning : 
    From function lduMatrix::operator-=(const lduMatrix& A)
    in file matrices/lduMatrix/lduMatrix/lduMatrixOperations.C at line 286
    Unknown matrix type combination
    this : diagonal:0 symmetric:0 asymmetric:1
    A    : diagonal:0 symmetric:0 asymmetric:0
   Normalisation factor = 0.0666275
DILUPBiCG:  Solving for b, Initial residual = 2.03922e-05, Final residual = 1.16209e-07, No Iterations 1
min(b) = 0.000161812
   Normalisation factor = 0.11664
DILUPBiCG:  Solving for Xi, Initial residual = 1.99508e-05, Final residual = 1.76309e-09, No Iterations 1
max(Xi) = 2.94208
max(XiEq) = 5.07445
Combustion progress = 3.96692%
--> FOAM Warning : 
    From function lduMatrix::operator-=(const lduMatrix& A)
    in file matrices/lduMatrix/lduMatrix/lduMatrixOperations.C at line 286
    Unknown matrix type combination
    this : diagonal:0 symmetric:0 asymmetric:1
    A    : diagonal:0 symmetric:0 asymmetric:0
   Normalisation factor = 42.1606
DILUPBiCG:  Solving for hau, Initial residual = 0.00717122, Final residual = 3.04963e-07, No Iterations 1
--> FOAM Warning : 
    From function lduMatrix::operator+=(const lduMatrix& A)
    in file matrices/lduMatrix/lduMatrix/lduMatrixOperations.C at line 207
    Unknown matrix type combination
    this : diagonal:0 symmetric:0 asymmetric:1
    A    : diagonal:0 symmetric:0 asymmetric:0
   Normalisation factor = 304.134
DILUPBiCG:  Solving for ha, Initial residual = 0.00202257, Final residual = 1.03208e-07, No Iterations 1
--> FOAM Warning : 
    From function lduMatrix::operator-=(const lduMatrix& A)
    in file matrices/lduMatrix/lduMatrix/lduMatrixOperations.C at line 286
    Unknown matrix type combination
    this : diagonal:0 symmetric:1 asymmetric:0
    A    : diagonal:0 symmetric:0 asymmetric:0
   Normalisation factor = 8.74103e-05
DICPCG:  Solving for p, Initial residual = 0.0442298, Final residual = 2.8683e-07, No Iterations 12
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 3.6708e-11, global = -1.58646e-11, cumulative = 2.29312e-05
   Normalisation factor = 1.39986
DILUPBiCG:  Solving for epsilon, Initial residual = 0.0381131, Final residual = 1.69447e-07, No Iterations 1
   Normalisation factor = 0.0180495
DILUPBiCG:  Solving for k, Initial residual = 0.0129366, Final residual = 9.14256e-07, No Iterations 1
ExecutionTime = 0.68 s  ClockTime = 0 s

If you plot the "Normalisation factor" for each equation, it might give you a sense of the actual level of error that is going on here in each solution... namely that if both values are fairly small, then the normalization isn't as accurate...

The way the "fvSchemes" file is configured does let me wondering if this case is just an approximate model. What I mean is that when the residuals stay so high, it usually is because the modelling is not accurate enough and is missing a few relevant details... and the extensive use of "limited" in this file leads me to suspect that the real case should have a much finer mesh and that it should not be defined with symmetry boundaries... or at least the boundaries opposite to the center of combustion should not be symmetry and should instead be farther away...


Lastly, I have no idea who worked on this solver originally and perhaps the OpenFOAM Foundation could provide some more insight into it if you contact them. For example, if this was just a proof of concept, or if there was a thesis associated to it or even if this solver is even finished... for example, fireFoam is still in ongoing development or pending development, depending on who you ask about it...