cellLimited<cubic> limiter is not differentiable

The cubic polynomial in cubicGradientLimiter.H is

$$((ar + b)r + 1)r \\ a = \frac{2}{r_t^2} - \frac{2}{r_t^3} \\ b = -\frac32 a r_t$$

but surely it should be

$$((ar + b)r + 1)r \\ a = \frac{1}{r_t^2} - \frac{2}{r_t^3} \\ b = -\frac32 a r_t - \frac{1}{2r_t}$$

to satisfy the constraints from Michalak and Ollivier-Gooch (2008, p.5), specifically that the derivative at the transition point

$r_t$

should be zero so it joins up smoothly with 1. This problem is illustrated in the graph below:

Linking this cfd-online thread here just in case someone replies to it after I create this issue.

assigned to @kuti

added bug label

@scleakey, thank you very much for this Wakabayashi-style catching! :)

One q: have you experienced any practical issues in your sims due to this bug? If so, could you elaborate on them a bit?

I assume you would be OK if we put your name into the bug-fix commit. (please feel free to provide a git patch if you want).

@kuti I've been programming a Godunov method solver in OpenFOAM and realised I could do MUSCL reconstruction using fvc::grad() as long as I used one of cellLimited, cellLimited<Venkatakrishnan>, cellLimited<cubic>, cellMDLimited, faceLimited, or faceMDLimited.

The problem was that my solver uses artificial compressibility which means it has to converge in pseudo-time each real-time step, and it was having trouble converging with any of the limiter options. I read that the limiter needs to be differentiable for convergence to steady-state, so had a delve into the source code... cellLimited is not differentiable, cellLimited<Venkatakrishnan> is not differentiable as it clips the limiter at 1 and, as established above, cellLimited<cubic> is not differentiable either even though it claims to be.

I have implemented a correct version of cellLimited<cubic> myself which does converge a little bit better, but not actually all that much. I think the real game-changer is implementing the other suggestion in the Venkatakrishnan and Michalak-Ollivier-Gooch papers: stopping the limiter from activating in uniform flow regions. I tested this for like 30 mins on Friday and it seems very promising but obviously I need to do a bit more testing.

To be honest, I don't think my git skills stretch far enough to do a git patch, but I would be fine if you did it and put my name into the bug-fix commit.

Thanks @scleakey, I will prepare the commit on behalf of you.

I think the real game-changer is implementing the other suggestion in the Venkatakrishnan and Michalak-Ollivier-Gooch papers: stopping the limiter from activating in uniform flow regions. I tested this for like 30 mins on Friday and it seems very promising but obviously I need to do a bit more testing.

If those separate changes could be promising, could you please lead us to make those changes in the software as well?

@kuti Sure, I will get back to you once I tested those other changes a bit more - I programmed them in a Friday afternoon rush and I am very suspicious that they worked the first time I tried!

Many thanks @scleakey

• Could you please check, and if possible confirm, the content of the patch (with your name - I have skipped your e-mail, let me know it if need be): 0001-BUG-cellLimited-cubic-ensure-the-limiter-is-differen.patch.
  • You can open it with a text editor, or apply it as a git patch onto your local branch (e.g. git am < 0001-*).
• Some optional qs:
  • Could you please elaborate a bit more on how you have derived the equation you stated?
  • Is there any minimal test case or test scenario from a paper that we can observe/test any improvements?

Thanks a lot.

mentioned in commit fe89c07c43f4fb99920d20225436fbf5c7a3ce2f

@kuti, thanks for putting together the patch!

Patch

I think it should be

a_((rt_ - 2)/pow3(rt_)),
b_(-(3*a_*sqr(rt_) + 1)/(2*rt_))

instead of

a_((rt_ - 2)/pow3(rt_)),
b_(-(3*a_*sqr(rt_) + 1)/(2*a_))

(the very last variable is different).

Derivation

Let me know if you spot any errors!

Test case

Off the top of my head, I can't really think of any case to demonstrate this that doesn't require making a new solver. And even then, you would have to hope that the limiter would want to oscillate between the two sides of the non-differentiable point. In the new solver I'm making, sometimes this happens towards the end of the simulation after the whole thing going fine until that point

Excellent explanations, highly appreciated.

I think the patch is good to go. Please don't worry about the test case.

Many thanks!

mentioned in commit 77e98b0e

mentioned in commit e861a19f

mentioned in commit ab49eaf9

@scleakey ,

Many thanks for your help and work. Highly appreciated. Added your commit into the develop branch and your name into the contributors' list.

Please consider to mention in the CFD-Online forum that the issue has been resolved, if possible.

Also please do not hesitate to lead us to make improvements in the code on anything.

Many thanks again.

closed