
[Sponsors] 
October 25, 2018, 09:28 
XiFoam Residuals flattening

#1 
New Member
Joaquín Aranciaga
Join Date: Oct 2018
Posts: 19
Rep Power: 3 
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.holzmanncfd.de. 

March 31, 2019, 16:00 

#2 
Super Moderator
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 10,956
Blog Entries: 43
Rep Power: 122 
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 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 1e2 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:
Best regards, Bruno
__________________
Last edited by wyldckat; March 31, 2019 at 16:01. Reason: There was nothing quick about this... removed "Quick answer" and went with my usual greetings... 

April 1, 2019, 07:58 

#3  
New Member
Joaquín Aranciaga
Join Date: Oct 2018
Posts: 19
Rep Power: 3 
Hi Bruno, thanks for taking the time to look at this.
Quote:
Quote:
Quote:
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 

April 2, 2019, 20:44 

#4 
Super Moderator
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 10,956
Blog Entries: 43
Rep Power: 122 
Quick answer: Then my last suspects are:
__________________


April 8, 2019, 11:57 

#5 
New Member
Joaquín Aranciaga
Join Date: Oct 2018
Posts: 19
Rep Power: 3 
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 zipup check OK. Faceface connectivity OK. <<Writing 4 cells with two nonboundary 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 (nonclosed singly connected) (0 0 0) (0 0.035 0.001) right 35 72 ok (nonclosed singly connected) (0.07 0 0) (0.07 0.035 0.001) top 70 142 ok (nonclosed singly connected) (0 0.035 0) (0.07 0.035 0.001) bottom 70 142 ok (nonclosed singly connected) (0 0 0) (0.07 0 0.001) frontAndBack 4900 5112 ok (nonclosed 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 (nonempty/wedge) directions (1 1 0) Mesh has 2 solution (nonempty) directions (1 1 0) All edges aligned with or perpendicular to nonempty directions. Boundary openness (2.98368e18 5.63584e18 2.38694e16) OK. Max cell openness = 1.05879e16 OK. Max aspect ratio = 1 OK. Minimum face area = 1e06. Maximum face area = 1e06. Face area magnitudes OK. Min volume = 1e09. Max volume = 1e09. Total volume = 2.45e06. Cell volumes OK. Mesh nonorthogonality Max: 0 average: 0 Nonorthogonality check OK. Face pyramids OK. Max skewness = 8.32667e14 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 

April 21, 2019, 08:35 

#6 
Super Moderator
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 10,956
Blog Entries: 43
Rep Power: 122 
Quick answer: I've left this on my todo 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 "OpenFOAM2.2.x/etc/controlDict" and change this: Code:
lduMatrix 1; Code:
lduMatrix 2; 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.98411e08, No Iterations 1 Normalisation factor = 0.328259 DILUPBiCG: Solving for Uy, Initial residual = 0.00505369, Final residual = 4.55551e08, No Iterations 1 Max StCourant 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.67585e08, No Iterations 1 min(b) = 0.000161812 Normalisation factor = 0.116263 DILUPBiCG: Solving for Xi, Initial residual = 0.000866927, Final residual = 1.42646e08, 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.40677e06, 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.5201e06, 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.97874e06, global = 1.00026e06, cumulative = 2.29312e05 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.5932e05, Final residual = 3.18415e10, No Iterations 1 Normalisation factor = 0.32927 DILUPBiCG: Solving for Uy, Initial residual = 3.37783e05, Final residual = 4.70605e10, No Iterations 1 Max StCourant 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.03922e05, Final residual = 1.16209e07, No Iterations 1 min(b) = 0.000161812 Normalisation factor = 0.11664 DILUPBiCG: Solving for Xi, Initial residual = 1.99508e05, Final residual = 1.76309e09, 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.04963e07, 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.03208e07, 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.74103e05 DICPCG: Solving for p, Initial residual = 0.0442298, Final residual = 2.8683e07, No Iterations 12 diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0 time step continuity errors : sum local = 3.6708e11, global = 1.58646e11, cumulative = 2.29312e05 Normalisation factor = 1.39986 DILUPBiCG: Solving for epsilon, Initial residual = 0.0381131, Final residual = 1.69447e07, No Iterations 1 Normalisation factor = 0.0180495 DILUPBiCG: Solving for k, Initial residual = 0.0129366, Final residual = 9.14256e07, No Iterations 1 ExecutionTime = 0.68 s ClockTime = 0 s 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...
__________________


Tags 
combustion, convergence, openfoam, xifoam 
Thread Tools  Search this Thread 
Display Modes  


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