Central Upwind Flux (rhoCentralFoam) in Implicit LUSGS Solver
 

I'm trying to use the central upwind flux featured in rhoCentralFoam in an implicit LUSGS based time integration scheme. The code already uses flux splitting schemes like your classic roe, HLLC etc. and as far as I understand it I should be able to use the central upwind scheme in rhoCenralFoam to calculate the inviscid fluxes and plug them into my LUSGS sweeps. Upon implementation this would appear to not be the case. I get appropriate fluxes (identical to roe) on initialization but when the code iterates a bit the fluxes become incorrect (in magnitude and direction) and the solution diverges.

I've also tried to do it the other way, put an LUSGS time integration sweep in rhoCentralFoam but again I run into issues. The wedge test case for example converges but I end up with an incorrect shock angle. Note that this is an inviscid test case so we're only looking at the inviscid fluxes, nothing to do with the predictor/corrector method rhoCentralFoam uses for viscous fluxes.

I appear to be missing a subtly in the way the central upwind scheme works. Has anyone else has success using it in an implicit solve?

Hello Chris,

I am also facing convergence issues with LUSGS methods. Initially, I tried to use the in-built TVD methods like you but I eventually implemented the Venkatakrishnan limiter as a cellLimtedGrad scheme and used the limited gradient to interpolate the pos & neg values. I tried running the wedge case you mentioned and I managed to achieve convergence (res < 1e-10) within 500 iterations. You may refer to this post where I have posted part of my code if you're interested. I tried running the wedge case using vanLeer and vanAlbada limited methods but the solution did not converge until 1000 iterations. I didn't run beyond that.

On the other hand, I am having still convergence issues with other cases (e.g. low-speed flow past cylinder). I have tried using several flux schemes (HLL, Rusanov, AUSM+, SLAU, ...) but the problem still persists. I am not sure why. It would help me if you could share your thoughts, insights or experience regarding implementing implicit methods since I could be missing something...

USV

Are you getting the correct solution with the wedge case and the default rhoCentral flux scheme? I was able to obtain convergence but the shock angle was completely wrong. If you've got an LUSGS sweep that does obtain a good solution I'd be very interested in taking a look at it. I'm specifically interested in the rhoCentral scheme because i'd like to tests the merits of the central-upwind scheme in a number of cases where pure upwind may not be the best (but these cases require an implicit solver for the sake of time).

Regarding convergence, I've never had a whole lot of luck getting LUSGS methods to converge with more complicated cases. Simple cases are typically fine but once you throw in some complexity residuals can be difficult to drive low. Switching the limiters is likely influencing convergence simply because they're bouncing back and forth between cells setting them first/2nd order and not arriving at a steady value. Particularly vanLeer which is quite sharp. VK's limiter is smoother and thus would provide better convergence.

What we've typically done with LUSGS is if residuals end up stalling, but the solution is stable, you can judge steady state on your own and judge if the solution is converged or not. Not ideal but it gets you a solution.

Your A-symmetry in the cylinder case is another beast altogether, I don't have any ideas on that one.

Hello Chris,

I am able to get the correct solution using LUSGS sweeps. The theoretical shock angle is marked with the white line. In fact, the shock angle is more accurate than rhoCentralFoam at the set endTime=0.2. The LUSGS results are obtained at the 500th iteration. It was simulated using the HLL scheme and Venkatakrishnan limiter.

I am under the impression that with the current implementation of rhoCentralFoam, the Kurganov and Tadmor flux schemes correspond to HLL and Rusanov schemes, respectively. Please correct me if I am wrong.

I agree with you. In fact, the wedge case is the only one for which I am able to achieve convergence using 2nd order methods. The solver converges nicely for 1st order methods for all other cases. So, I am quite sure that the solver itself is implemented correctly. I wonder how others manage to drive their residuals low when using higher order schemes...

Could anybody help us on this point?

USV

Generally speaking, if you want to drive your residuals lower you're going to need a better implicit scheme. LUSGS simply bottoms out or takes a long time to do so. The link below is an excellent paper that goes over how LUSGS converges relative to a fully implicit matrix solver as well as describes a block-coupled solver. (not sure the matrix libraries required are available in your version of openfoam though, extend should, in theory, be able to implement the block-coupled approach.) How you define your residuals also plays a role. Ultimately residuals are an estimate of how good your solution is, not the be-all end-all of determining if you have a good solution, at least in a practical sense.

https://www.researchgate.net/profile...8e348389f6.pdf

Regarding the flux schemes I do believe they are different. The HLL scheme is a Riemann solver and one of the main proposed advantages of the central-upwind scheme described is that it is Reimann free. Also, the paper by Kurganov and Tadmor certainly describes a novel scheme. (I could also be mistaken, it has been a while since I did a proper review)

This really is more of an academic exercise for myself though, wanting to look deeper into the advantages and disadvantages of the central-upwind scheme. I would be curious, if you have the time and feel like indulging me, as to what results you get on the wedge when you use the fluxes obtained following the exact original rhoCentralFoam approach. You should be able to apply the VK limiter as well if needed. You should (I think) be able to basically just use the phi, phiUp, phiEp for the inviscid mass, momentum and energy fluxes. This really goes back to my original question as to what I did wrong with my LUSGS attempts and I could very well be making an incorrect interpretation of the code here.

Hello Chris,

Thank you for the reply.

Could you give an example of a different way to define residuals? Are you referring to other quantities like CL and CD which I've seen being used to check for convergence in some papers? These quantities may not be applicable to all problems though.

I agree with you on the point that, conceptually, HLL and Kurganov-Tadmor schemes are different but their implementations are identical when using 2nd order TVD schemes. This is what I realized when I started to read up on rhoCentralFoam a couple of years ago. Perhaps someone could back me up on this...

I tried doing that replacing all the TVD interpolations with my VK limiter. You were right, I couldn't get the correct shock angle either. Please refer to the first attached image(white line is exact shock angle). However, once I explicitly applied simple reflective boundary conditions across the 'bottom' and 'obstacle' patches, I was able to get the correct shock angle (see second attached image)! So, you may try doing the same and it might solve your problem.

USV

Thanks for taking a look at that. Your initial "incorrect" shock results were considerably closer than the ones I got. My shock ended up impinging on the top boundary it was so far off. Looks like something is incorrect or missing in my implementation (likely proper limiters)! Good to know, thank you.

As far as residual definition goes, yes you could look at a variable of interest. I was mostly referring to how you define your initial residual (assuming you're normalizing), you could take your very first value or average over the first few iterations etc. Most days it usually just comes down to converging until the residuals flatline or the solution stops changing in a meaningful way. It's non-ideal in that it requires user intervention and judgement but even commercial codes like Star-CCM (which uses a fully implicit matrix method) don't always drive residuals to zero.

As a side note, you could also try checking the Jacobian approximations (spectral radius in the original LUSGS), that's what's going to influence your convergence rate the most. My knowledge here isn't the strongest but an error in how you define those or even finding some literature on alternate definitions may help out (if they exist).

That's interesting regarding HLL vs. Kurganov-Tadmor, I'll look into it deeper.

Is there any implicit LUSGS Solver available in OpenFOAM ?
 

Dear Usv100, Dear Chris,

Is it possible for you to share your Implicit LUSGS Solver ? I am running my cases with rhoCentralFoam and as I explained here :
https://www.cfd-online.com/Forums/openfoam-solving/216125-help-convergence-time-rhocentralfoam-help.html

I have real troubles with the convergence time of rhoCentralFoam. Usv100 has suggested me to use an implicit LUSGS Solver. I contacted the first author of the first paper being mentioned in the comments ( C. Shen, X.-L. Xia, Y.-Z. Wang, F. Yu, and Z.-W. Jiao, Implementation of density-based implicit LU-SGS solver in the framework of OpenFOAM, Adv. Eng. Softw. 91, (2016) ). However, I did not get any answer. I am wondering if you have submitted your solver to OpenFoam in order to integrate it within their package. I think that would help lots of people like me.

Best regards,
Mary