CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Running, Solving & CFD

XiFoam Residuals flattening

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree1Likes
  • 1 Post By wyldckat

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   October 25, 2018, 10:28
Default XiFoam Residuals flattening
  #1
New Member
 
Joaquín Aranciaga
Join Date: Oct 2018
Posts: 15
Rep Power: 3
joaran is on a distinguished road
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.
Attached Images
File Type: jpg Residuals1.jpg (109.9 KB, 16 views)
Attached Files
File Type: txt fvSolution.txt (2.3 KB, 4 views)
joaran is offline   Reply With Quote

Old   March 31, 2019, 17:00
Default
  #2
Super Moderator
 
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 10,955
Blog Entries: 43
Rep Power: 121
wyldckat is a name known to allwyldckat is a name known to allwyldckat is a name known to allwyldckat is a name known to allwyldckat is a name known to allwyldckat is a name known to all
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:
Code:
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 likes this.
__________________

Last edited by wyldckat; March 31, 2019 at 17:01. Reason: There was nothing quick about this... removed "Quick answer" and went with my usual greetings...
wyldckat is offline   Reply With Quote

Old   April 1, 2019, 08:58
Default
  #3
New Member
 
Joaquín Aranciaga
Join Date: Oct 2018
Posts: 15
Rep Power: 3
joaran is on a distinguished road
Hi Bruno, thanks for taking the time to look at this.

Quote:
Originally Posted by wyldckat View Post

So the very first question is: Are you using an SP or DP build? If you run this command:
Code:
echo $WM_PRECISION_OPTION
it should tell you which one you are using. You should be using DP.
I'm using DP.

Quote:
Originally Posted by wyldckat View Post
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.

Quote:
Originally Posted by wyldckat View Post
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
joaran is offline   Reply With Quote

Old   April 2, 2019, 21:44
Default
  #4
Super Moderator
 
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 10,955
Blog Entries: 43
Rep Power: 121
wyldckat is a name known to allwyldckat is a name known to allwyldckat is a name known to allwyldckat is a name known to allwyldckat is a name known to allwyldckat is a name known to all
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:
      Code:
      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.
__________________
wyldckat is offline   Reply With Quote

Old   April 8, 2019, 12:57
Default
  #5
New Member
 
Joaquín Aranciaga
Join Date: Oct 2018
Posts: 15
Rep Power: 3
joaran is on a distinguished road
Hi Bruno,


The mesh seems to be fine:


Code:
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
joaran is offline   Reply With Quote

Old   April 21, 2019, 09:35
Default
  #6
Super Moderator
 
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 10,955
Blog Entries: 43
Rep Power: 121
wyldckat is a name known to allwyldckat is a name known to allwyldckat is a name known to allwyldckat is a name known to allwyldckat is a name known to allwyldckat is a name known to all
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:
Code:
    lduMatrix           1;
to this:
Code:
    lduMatrix           2;
you will get a few more details about each equation being solved, e.g.:
Code:
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...
__________________
wyldckat is offline   Reply With Quote

Reply

Tags
combustion, convergence, openfoam, xifoam

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
plot residuals pigna OpenFOAM Running, Solving & CFD 2 December 19, 2014 03:32
motorBike Residuals for SST k-omega... and mine JR22 OpenFOAM Running, Solving & CFD 6 August 1, 2013 10:08
what to monitor besides residuals? franzdrs FLUENT 5 March 21, 2013 04:59
judging convergence through residuals MachZero Main CFD Forum 7 December 25, 2012 13:18
Convergence - scaled vs unscaled residuals HS FLUENT 1 November 7, 2005 06:45


All times are GMT -4. The time now is 14:50.