Add flux-limiter FEM-TVD discretisation to PorousFlow

[WilkAndy]

Tagging @cpgr @hsheldon

Rationale

The full upwinding used in PorousFlow is too diffusive for problems involving the tracking of fronts, so a new type of upwinding needs to be implemented to simulate these problems.

The Richards module has both SUPG and full upwinding. SUPG is less diffusive than full upwinding, but has proven to be useless in many porous-flow simulations because it doesn't prevent fluid from being withdrawn from a fluid-less node. Therefore only full upwinding has been implemented in the PorousFlow module.

Unfortunately, full upwinding is also very diffusive, so sharp fronts get smoothed. @hsheldon and I discussed implementing SUPG for these types of problems but eventually decided it wasn't worth it: (1) it may not really solve the problem, and indeed we found off-hand comments in papers that suggested it wouldn't; (2) it's very complicated, being dependent on element geometry (things like dxi_dxyz) and derivatives of nonlinear functions (so second-derivatives enter the Jacobian, uggh).

So we settled on trying to implement a flux-limiter TVD discretisation. This may or may not solve the front-tracking problem, but at least will offer users the possibility of another discretisation that could improve convergence, especially in close-to-steadystate situations.

Description

I'm hoping to implement the approach described in

D Kuzmin and S Turek "High-resolution FEM-TVD schemes based on a fully multidimensional flux limiter" Journal of Computational Physics, Volume 198 Issue 1, 20 July 2004, Pages 131 - 158

This may end up being a large and complicated issue to implement, involving Kernels, BCs and Materials.

Impact

Ability to use another type of upwinding in PorousFlow

[WilkAndy]

One of the first questions is what is the KuzminTurek "u" in the PorousFlow setting? They start with a nonlinear continuity equation (Eqn(1)) but the paper is written for the form Eqn(15) and Eqn(31). Eqn(31)'s "k" can depend on u, but still the PorousFlow equations have to be massaged into that form.

The PorousFlow equations involve time derivatives of the mass of a species. That is a sum over all phases of density saturation mass_fraction * porosity. This would be the natural definition of "u". However, those terms don't appear in the Darcy flux. Therefore, i'm going to use as my "u" variable

u = density * mass_fraction

There is a different "u" for each phase of the species. A sum (over phases) of all these "u" variables, multiplied by saturation and porosity, appears in the PorousFlow time derivative, but i don't think this added complexity invalidates the KuzminTurek approach. For if the advection is correctly discretised for each phase then the summed version should display the same mass conservation and convergence properties. Nevertheless, this is simply my choice and it could lead to bad results!

In @hsheldon 's case where she has 1 phase and almost constant density and unity saturation, then choosing u ~ mass_fraction is probably a good thing.

[WilkAndy]

Another question is how exactly to form the K matrix in Eqn(31). The standard MOOSE approach differs from KuzminTurek Eqn(17). MOOSE expands the NonlinearVariables, eg porepressure, as

P = P_j * phi_j

where there is a sum over "j", and P_j is the value of porepressure at node j, and phi_j is the j^th shape function. Using phi_j evaluated at quadpoints, functions such as saturation, relative permeability, etc, are evaluated at the quadpoints.

The problem with using KuzminTurek's Eqn(17) is that we have to evaluate the Darcy velocity at the nodes, and currently in PorousFlow we only evaluate things like permeability and grad(P) at the quadpoints. Eqn(17) is quite advantageous because it leads to Eqn(20) which means expensive quad-point integration only has to be performed once. Nevertheless, i don't really want to screw around with MOOSE this much, so i will not use Eqn(17) exactly.

Instead, the c matrix in Eqn(20) will be evaluated using sum-over-quadpoints of (darcyQp() * phi_j), and this will be performed every Residual and Jacobian evaluation. The nodal term (called "v") in Eqn(17) and Eqn(20) then involves mobility() / u. To evaluate this expression i'll have to write another function in PorousFlowDarcyBase and all derived classes.

[WilkAndy]

For the heat equation, i suppose u = enthalpy * density

[cpgr]

Some light holiday reading!

[WilkAndy]

You lucky thing, i'm not on holidays yet

[hsheldon]

As far as my understanding goes (which is not very far!), I think your choices of u for mass and heat seem reasonable. Let's hope it works!

[WilkAndy]

So i've found out the the RDG module might be helpful to accurately model advection in PorousFlow. Currently RDG is more of a set of tools rather than a thing we can use immediately and directly. But @yidongxiainl has got an App that uses it to do advection of species without overshoots and diffusion. He, and/or a student of @permcody , and/or other MOOSE team people are planning to work on this.

I'm hoping we can start something before the middle of April. I know this is not altogether straightforward because it involves screwing around with MOOSE at a fundamental level.

Probably we can start small - just setting up a PorousFlow simulation with advection in 1D that uses the RDG idea, and see whether we get good results. But maybe @yidongxiainl or @permcody or @rpodgorney or Rich Martineau have other ideas.

The accurate advection of (chemical) species without diffusion is also needed by the geothermal group at ETH, and I believe it will eventually act as a fundamental stumbling block to everyone who uses PorousFlow with nontrivial fluid components.

At the end of the day, i think what we need is the advection equation (probably usually with extra Kernels for diffusion, sources/sinks corresponding to chemistry, etc) where the velocity is given as a Material property that depends on all sorts of other things in the problem.

Tagging @cpgr once again in case he doesn't get auto notification.

Tagging @rpodgorney so that he's aware of this development.

[hsheldon]

That sounds promising! Do you need anything from me at this stage?

[hsheldon]

Is there any documentation for RDG that I can look at, or perhaps some examples?

[yidongxiainl]

Hi Andy, Heather and all,

Sorry for the delayed response, as I got no email access for the past few days.

The currently optimal RDG implementation is the multi-dimensional Euler equations for computational fluid dynamics, which is currently in an in-house app called Bighorn. The source code was recently moved out of the moose physics module collection back into house, as it was needed for some export-control related new app development.

As I told Andy last week, another RDG implementation is an in-house porous flow & heat transport app called Badger. This is the very first RDG based app developed by me and Hai Huang. In order to move the code outside to GitHub, I would to consult with Cody about necessary paperwork of open-source licensing.

For the time being, I recommend this article < https://link.springer.com/article/10.1007/s00603-016-0951-y> for reference.

Yidong Xia, Ph.D. | Computational Scientist Materials Science & Engineering Department Idaho National Laboratory 1955 N. Fremont Ave., PO Box 2211, Idaho Falls, ID 83415-2025 Tel: (208) 526-7490

On Mon, Mar 19, 2018 at 4:56 PM, Heather Sheldon notifications@github.com wrote:

Is there any documentation for RDG that I can look at, or perhaps some examples?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_idaholab_moose_issues_10426-23issuecomment-2D374413195&d=DwMCaQ&c=54IZrppPQZKX9mLzcGdPfFD1hxrcB__aEkJFOKJFd00&r=cyGw8U5koLkfOR3LoksE-oYORShC_AQreo7-FFvKKzw&m=5cMLqctMaKgn5hqXhkav5T91GH-ejGQycDCEDC5bg1E&s=XjZ5Wb36qX6iihgvpnEs6ktAv641JBH7ichxFSFeRcQ&e=, or mute the thread https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_notifications_unsubscribe-2Dauth_AJ7EvPLVjrklIdVUQJk-2DZeSuKdBLet3-2Dks5tgDebgaJpZM4RGSuG&d=DwMCaQ&c=54IZrppPQZKX9mLzcGdPfFD1hxrcB__aEkJFOKJFd00&r=cyGw8U5koLkfOR3LoksE-oYORShC_AQreo7-FFvKKzw&m=5cMLqctMaKgn5hqXhkav5T91GH-ejGQycDCEDC5bg1E&s=xhvsYDt-E19-Qb4F9Tt43Zne0LQAhoLqbeQdJqo0Mgc&e= .

[WilkAndy]

Thanks for the article, @yidongxiainl . I can't download it from the hotel, but the article looks exciting. I discussed with this with @permcody and Rich M, explaining that this is a critical (or THE critical) aspect missing from PorousFlow, and their opinion was that since RDG itself is not export controlled, it should be fine to move it from Badger (or Bighorn? i'm not sure) into modules. Of course i'm not interested in any secret stuff - i just want to solve \dot{u}+grad(v u) = 0, for my nasty v.

[WilkAndy]

Yes, @hsheldon , look at http://www.mooseframework.org/moose/documentation/modules/rdg/index.html

[yidongxiainl]

Andy, yes, the example in the RDG module was so far the "optimal" implementation based on the moose framework. Together with David Andrs, we removed all the overheads we could identify.

The full-text of the referred article can be downloaded at my ResearchGate page < https://www.researchgate.net/publication/298912690_Assessment_of_a_Hybrid_ContinuousDiscontinuous_Galerkin_Finite_Element_Code_for_Geothermal_Reservoir_Simulations

Yidong Xia, Ph.D. | Computational Scientist Materials Science & Engineering Department Idaho National Laboratory 1955 N. Fremont Ave., PO Box 2211, Idaho Falls, ID 83415-2025 Tel: (208) 526-7490

On Mon, Mar 19, 2018 at 11:33 PM, Andy Wilkins notifications@github.com wrote:

Yes, @hsheldon https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_hsheldon&d=DwMCaQ&c=54IZrppPQZKX9mLzcGdPfFD1hxrcB__aEkJFOKJFd00&r=cyGw8U5koLkfOR3LoksE-oYORShC_AQreo7-FFvKKzw&m=WU1X7y788NkRVa_qMzjjffVsl-pIsWhNOOB5JiXesQ0&s=WWG1W4Xk1fhZnhvKFd3VWOVLJuNr5crAEKb7K1JFNaY&e= , look at http://www.mooseframework.org/moose/documentation/modules/ rdg/index.html https://urldefense.proofpoint.com/v2/url?u=http-3A__www.mooseframework.org_moose_documentation_modules_rdg_index.html&d=DwMCaQ&c=54IZrppPQZKX9mLzcGdPfFD1hxrcB__aEkJFOKJFd00&r=cyGw8U5koLkfOR3LoksE-oYORShC_AQreo7-FFvKKzw&m=WU1X7y788NkRVa_qMzjjffVsl-pIsWhNOOB5JiXesQ0&s=fgFFNDv2L3EXJTvBkAbvBuyFZT3i6whP9tK7qUaoqVk&e=

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_idaholab_moose_issues_10426-23issuecomment-2D374481387&d=DwMCaQ&c=54IZrppPQZKX9mLzcGdPfFD1hxrcB__aEkJFOKJFd00&r=cyGw8U5koLkfOR3LoksE-oYORShC_AQreo7-FFvKKzw&m=WU1X7y788NkRVa_qMzjjffVsl-pIsWhNOOB5JiXesQ0&s=-cB4LGL3OgXl-xnVlGnlghnxvUfRaSKPSy3niI7Cw-Y&e=, or mute the thread https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_notifications_unsubscribe-2Dauth_AJ7EvDMQwT8s7H5lVsDfOgPqenIaBFVsks5tgJSQgaJpZM4RGSuG&d=DwMCaQ&c=54IZrppPQZKX9mLzcGdPfFD1hxrcB__aEkJFOKJFd00&r=cyGw8U5koLkfOR3LoksE-oYORShC_AQreo7-FFvKKzw&m=WU1X7y788NkRVa_qMzjjffVsl-pIsWhNOOB5JiXesQ0&s=kKoBbtC4-A2N-qvXUwYVIotagmyj3kgSGp_frZXukF8&e= .

[WilkAndy]

Great, thank you for the paper @yidongxiainl . I'll probably read it on my way back home