Implement least-squares slope reconstruction in RDG

[WilkAndy]

Rationale

There is no slope reconstruction in RDG at the moment, apart from what is implicit in AEFVSlopeLimitingOneD. Implement least-squares slope reconstruction for multi-dimensions. Focus only on reconstructing linear solutions from constant, ie the P0P1 RDG method.

Discussion point I think that slope reconstruction is the first step in enhancing the RDG model to be useful to other modules like PorousFlow. Is this a good place to start, or should i spend my energies on something else: either another sort of slope reconstruction, or a DGKernel, or something else?

Tagging @permcody , @yidongxiainl , @andrsd , @joshuahansel , @mrodb , @rpodgorney , @lingzou , @cpgr , @YaqiWang , @philippschaedle , @hsheldon

Description

This slope reconstruction will implement Eqns(28)-(30) in the bighorn doco https://github.com/idaholab/moose/files/2440636/xia2016bighorn.pdf .

1. I think this implementation should be a UserObject inheriting from SlopeReconstructionMultiD.
2. I think i'll use a Householder transform, as suggested in https://www.researchgate.net/publication/2770667_High-Order_ENO_Schemes_for_Unstructured_Meshes_Based_on_Least-Squares_Reconstruction to solve the least-squares system.
3. I think i should implement this for one Variable only, and if the user has multiple constant-monomials then they will have multiple UserObjects. Or should I implement this to reconstruct one Material property? I think PorousFlow might like the latter.
4. To test this beast, i suspect i'll need another piece of code that extracts the results to an AuxVariable or something that i can CSVDiff.
5. I don't know what to do about Jacobians: the reconstructed slope will depend on the Variables (or Material properties) at this element and neighbouring elements.

Discussion points Please comment on the above statements!

Impact

Ability to do reconstruct linear Variables or Material properties from piecewise-constant versions, as a start of doing P0P1 RDG.

[joshuahansel]

I think it'd be great to make some base classes for slope reconstruction and slope limiting that are general to different systems of equations. You mentioned the implicit slope reconstruction in AEFVSlopeLimitingOneD - if I recall correctly, the 1-D slope limiters operated on 3 different reconstructed slopes per element - a left slope, a central slope, and a right slope. I think the multi-D limiters operate on slopes from the current element and its neighbors, but only 1 per element. So I think that the implementations are inherently different, and perhaps this is why the reconstruction+limiting were both implemented in AEFVSlopeLimitingOneD. @yidongxiainl knows this better.

Comments on your numbered items:

1. Yes, I think so.
2. No idea - I personally have no experience with ENO/least-squares reconstruction.
3. Not sure about this. As far as a know, in general there is a transformation to a set of primitive variables, where the slopes are reconstructed, limited, and then a reverse transformation is made, so it may be less efficient to take a single-variable approach (I don't know how significant). And I don't understand your comment about a material property - I thought you were only talking about the user objects - the material just uses the slope limiter UO interface. I think there was a reason why it couldn't all (reconstruction+limiting) be done in a material (maybe because of the need of access to neighbor cell center values). This might be a good discussion for @andrsd.
4. Perhaps the user object(s) plus a material using them and output the material to CSV to compare to hand-calculated gold values? AuxVariable seems wrong to me because of the reconstruction nature.
5. I believe @yidongxiainl told me that it's acceptable to approximate the Jacobians by ignoring the slopes (just using cell center values). I'm not sure myself - I'd like to test this more rigorously when I have time.

[yidongxiainl]

I think it'd be great to make some base classes for slope reconstruction and slope limiting that are general to different systems of equations. You mentioned the implicit slope reconstruction in AEFVSlopeLimitingOneD - if I recall correctly, the 1-D slope limiters operated on 3 different reconstructed slopes per element - a left slope, a central slope, and a right slope. I think the multi-D limiters operate on slopes from the current element and its neighbors, but only 1 per element. So I think that the implementations are inherently different, and perhaps this is why the reconstruction+limiting were both implemented in AEFVSlopeLimitingOneD. @yidongxiainl knows this better.

Comments on your numbered items:

1. Yes, I think so.
2. No idea - I personally have no experience with ENO/least-squares reconstruction.
3. Not sure about this. As far as a know, in general there is a transformation to a set of primitive variables, where the slopes are reconstructed, limited, and then a reverse transformation is made, so it may be less efficient to take a single-variable approach (I don't know how significant). And I don't understand your comment about a material property - I thought you were only talking about the user objects - the material just uses the slope limiter UO interface. I think there was a reason why it couldn't all (reconstruction+limiting) be done in a material (maybe because of the need of access to neighbor cell center values). This might be a good discussion for @andrsd.
4. Perhaps the user object(s) plus a material using them and output the material to CSV to compare to hand-calculated gold values? AuxVariable seems wrong to me because of the reconstruction nature.
5. I believe @yidongxiainl told me that it's acceptable to approximate the Jacobians by ignoring the slopes (just using cell center values). I'm not sure myself - I'd like to test this more rigorously when I have time.

In response to @joshuahansel The slope reconstruction & limiting is based on what system of equations we are solving. It can be a single variable scalar transport equation (i.e. the one in the moose module) or the multi-variable Navier-Stokes / Euler equations. So in the relevant user object(s), we are doing loop over the elements, and in each element, we treat the variable(s) all together, but not a sub-loop over the variable(s). That was designed (a hard coded way) just to get our word done back then. It could be designed for more general purpose, allowing for arbitrary number of variables to go through the slope reconstruction & limiting. This is a code design problem, but not physics problem. @WilkAndy may be in a better position to make it such general purpose.

Again, yes, 1D is methodologically different than 2D and 3D for doing limiting. 1D (the TVD slope limiter case) requires only a one-step treatment to obtain the reconstructed & limiting slope. In the manual, it explains the related formula.

[yidongxiainl]

So in other words, I endorse the idea of making the slope reconstruction + limiting a more general purpose functionality in moose module. It is basically a design problem. It can be simple in some applications, but can also be tricky in some other. Let me make an example with the Navier-Stokes equations.

o-- For compressible Navier-Stokes equations, the variables we solve are the conserved variables, i.e. (rho, rhou, rhov, rhow, rhoe). But they are NOT recommended to be directly used in slope reconstruction & limiting. The reason is that the effect is not good. Instead, it is much better to do the slope reconstruction & limiting using the primitive variables, e.g. (p, u, v, w, T), or (rho, u, v, w, p). The process is, 1) calculate the primitive variables from the conserved variables for each element; 2) do the slope reconstruction & limiting based on the primitive variables; 3) now with cell-averages and their limiting slopes, the primitive variables at the cell faces can be calculated used to compute all needed variables for calculating fluxes.

This is an example to explain, slope reconstruction & limiting is sometimes more than a consideration of programming. Actual effect has to be also considered if a system of equations are dealt with.

For porous flow & heat transport, my past experience is just do pressure and temperature, which I found OK.

[permcody]

Note to all here. The RelationshipManagers are in and working well (both Geometric and Algebraic). It should be very easy now to have your UserObject register (tell MOOSE) that it needs neighbor information. Then you'll always have the necessary information as you run your algorithms.

[joshuahansel]

Note to all here. The RelationshipManagers are in and working well (both Geometric and Algebraic). It should be very easy now to have your UserObject register (tell MOOSE) that it needs neighbor information. Then you'll always have the necessary information as you run your algorithms.

Currently, we've only done only constant monomial (P0) solutions, and we're getting neighbor information via getElementalValue(neighbor), where we get neighbor via elem->neighbor_ptr(i_side). Is there an advantage to using a RelationshipManager if this is the only thing we need to do?

[permcody]

Is there an advantage to using a RelationshipManager if this is the only thing we need to do?

Absolutely! As we've been discussing here in the cube over the last few days. We've discovered that normal FE assembly doesn't require algebraic ghosting at all! When you start doing adaptivity, DG, or several other things, then you start needing it. There's an opportunity to optimize MOOSE for a lot of common cases where we can just turn off ghosting. Therefore objects that do intend to use neighboring information because it's alway been there shouldn't count on that in the future. Please register the appropriate relationship manager for your use.

[WilkAndy]

@permcody - i'm going to need help with the RelationshipManager stuff. Could you point me to a good example (as @joshuahansel said - i think i'm only going to need neighbour information).

@yidongxiainl - i was not planning on using P0P1 for P or T, so this needs further discussion, please. Currently P and T are solved well using the schemes already in PorousFlow. What is not solved well is the advection of fluid species, which are mass concentrations (kg of chemical / m^3 of fluid). That is: Darcy's law gives an advective velocity for the fluid species, and i just want to advect the mass-concentrations using this velocity. One complication is that the Darcy velocity can depend on the mass concentrations, because things like the fluid density and viscosity can depend on the concentrations.

This is why i mentioned Materials above: i thought i could calculate Darcy velocity as a Material, element-average it (P0) and reconstruct the linear form (P0P1), then use it to advect the mass concentrations (which are P0). I thought if RDG assumed the user provided the advective velocity as a Material Property (calculated as however they desired) the RDG module could do the remainder of the work.

Please tell me what you think, everyone. Yes, this is not just a code design question, but a physics question too, and i'm a real newbie in this field.

(Aside: i disagree about the Jacoabians as i can almost 100% guarantee that we'll need the full Jacobian for convergence.)

[yidongxiainl]

@permcody - i'm going to need help with the RelationshipManager stuff. Could you point me to a good example (as @joshuahansel said - i think i'm only going to need neighbour information).

@yidongxiainl - i was not planning on using P0P1 for P or T, so this needs further discussion, please. Currently P and T are solved well using the schemes already in PorousFlow. What is not solved well is the advection of fluid species, which are mass concentrations (kg of chemical / m^3 of fluid). That is: Darcy's law gives an advective velocity for the fluid species, and i just want to advect the mass-concentrations using this velocity. One complication is that the Darcy velocity can depend on the mass concentrations, because things like the fluid density and viscosity can depend on the concentrations.

This is why i mentioned Materials above: i thought i could calculate Darcy velocity as a Material, element-average it (P0) and reconstruct the linear form (P0P1), then use it to advect the mass concentrations (which are P0). I thought if RDG assumed the user provided the advective velocity as a Material Property (calculated as however they desired) the RDG module could do the remainder of the work.

Please tell me what you think, everyone. Yes, this is not just a code design question, but a physics question too, and i'm a real newbie in this field.

(Aside: i disagree about the Jacoabians as i can almost 100% guarantee that we'll need the full Jacobian for convergence.)

@WilkAndy Andy, surely you decide the form of governing equations, and we'll see how it goes from there. I use P and T based formulation as an example to explain in the case of porous flow we do not need to consider whether we should reconstruct & limit the derived variables instead of the primary variables. Once we look at your governing equations, then it is clear what to do next.

Regarding the discussion of Jacobians, I agree with you 100% if it is under the context of continuous FEM. But with regard to P0P1, we need to test and see the effect. I am holding an open mind in this.

[permcody]

If you just need side neighbor information. Just add this line to your validParams: https://github.com/idaholab/moose/blob/devel/test/src/userobjects/ElemSideNeighborLayersTester.C#L30

If you need some variable number of side neighbor layers, then you can either derive from the existing class and set it in code somewhere, or you can require that the user supply a parameter to indicate the number of layers like I've done in this particular object: https://github.com/idaholab/moose/blob/devel/test/src/userobjects/ElemSideNeighborLayersTester.C#L31-L32

That's it - pretty painless. MOOSE will take care of the rest.

On Mon, Oct 15, 2018 at 2:15 PM Andy Wilkins notifications@github.com wrote:

@permcody https://github.com/permcody - i'm going to need help with the RelationshipManager stuff. Could you point me to a good example (as @joshuahansel https://github.com/joshuahansel said - i think i'm only going to need neighbour information).

@yidongxiainl https://github.com/yidongxiainl - i was not planning on using P0P1 for P or T, so this needs further discussion, please. Currently P and T are solved well using the schemes already in PorousFlow. What is not solved well is the advection of fluid species, which are mass concentrations (kg of chemical / m^3 of fluid). That is: Darcy's law gives an advective velocity for the fluid species, and i just want to advect the mass-concentrations using this velocity. One complication is that the Darcy velocity can depend on the mass concentrations, because things like the fluid density and viscosity can depend on the concentrations.

This is why i mentioned Materials above: i thought i could calculate Darcy velocity as a Material, element-average it (P0) and reconstruct the linear form (P0P1), then use it to advect the mass concentrations (which are P0). I thought if RDG assumed the user provided the advective velocity as a Material Property (calculated as however they desired) the RDG module could do the remainder of the work.

Please tell me what you think, everyone. Yes, this is not just a code design question, but a physics question too, and i'm a real newbie in this field.

(Aside: i disagree about the Jacoabians as i can almost 100% guarantee that we'll need the full Jacobian for convergence.)

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/idaholab/moose/issues/12310#issuecomment-429997232, or mute the thread https://github.com/notifications/unsubscribe-auth/AC5XIH28S4-msXyH7rc1TwTd0rkC2iR9ks5ulOzegaJpZM4Xblfs .

[hsheldon]

@wilkandy – you say you are not planning to use P0P1 for P or T as they are already solved well – I may be missing something here (as in I haven’t followed this in detail so don’t actually now what P0P1 means!), but I would like to point out that T is not always solved well. I see oscillations in T without upwinding, but over-diffusive behaviour with upwinding, which is an issue for the problems I am trying to solve as it masks the behaviour of interest.

Heather

Subject: Re: [idaholab/moose] Implement least-squares slope reconstruction in RDG (#12310)

@permcodyhttps://github.com/permcody - i'm going to need help with the RelationshipManager stuff. Could you point me to a good example (as @joshuahanselhttps://github.com/joshuahansel said - i think i'm only going to need neighbour information).

@yidongxiainlhttps://github.com/yidongxiainl - i was not planning on using P0P1 for P or T, so this needs further discussion, please. Currently P and T are solved well using the schemes already in PorousFlow. What is not solved well is the advection of fluid species, which are mass concentrations (kg of chemical / m^3 of fluid). That is: Darcy's law gives an advective velocity for the fluid species, and i just want to advect the mass-concentrations using this velocity. One complication is that the Darcy velocity can depend on the mass concentrations, because things like the fluid density and viscosity can depend on the concentrations.

This is why i mentioned Materials above: i thought i could calculate Darcy velocity as a Material, element-average it (P0) and reconstruct the linear form (P0P1), then use it to advect the mass concentrations (which are P0). I thought if RDG assumed the user provided the advective velocity as a Material Property (calculated as however they desired) the RDG module could do the remainder of the work.

Please tell me what you think, everyone. Yes, this is not just a code design question, but a physics question too, and i'm a real newbie in this field.

(Aside: i disagree about the Jacoabians as i can almost 100% guarantee that we'll need the full Jacobian for convergence.)

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHubhttps://github.com/idaholab/moose/issues/12310#issuecomment-429997232, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AKJ_8naKhCTvmIVb_Mf7p2rKDfVBxw45ks5ulOzigaJpZM4Xblfs.

[WilkAndy]

@permcody - thanks for the Relationship manager stuff, i'll let you know if i need any more help

@hsheldon - thanks for the reminder about T advection. This is because T is being advected by the fluid, right? Re P0P1: it means T (and fluid species) will be constant monomial (that's the P0 part, the "P" stands for piecewise and the "zero" for constant: T is piecewise constant over the mesh) but that to compute fluxes, we "reconstruct" a piecewise linear function for T by looking at the elemental value of T and its nearest neighbours' elemental values of T (that's the P1 part - the "one" meaning "linear"). Using this reconstructed P0P1 function, we can "limit" it (by, eg, placing hard bounds on it to avoid crazy things like T<0) and also compute the fluxes by looking at the function at the element boundaries, but those things are not part of this github issue. Someone can correct me if i'm wrong.

[hsheldon]

@wilkandy – yes, because T is being advected by the fluid.

Subject: Re: [idaholab/moose] Implement least-squares slope reconstruction in RDG (#12310)

@permcodyhttps://github.com/permcody - thanks for the Relationship manager stuff, i'll let you know if i need any more help

@hsheldonhttps://github.com/hsheldon - thanks for the reminder about T advection. This is because T is being advected by the fluid, right? Re P0P1: it means T (and fluid species) will be constant monomial (that's the P0 part, the "P" stands for piecewise and the "zero" for constant: T is piecewise constant over the mesh) but that to compute fluxes, we "reconstruct" a piecewise linear function for T by looking at the elemental value of T and its nearest neighbours' elemental values of T (that's the P1 part - the "one" meaning "linear"). Using this reconstructed P0P1 function, we can "limit" it (by, eg, placing hard bounds on it to avoid crazy things like T<0) and also compute the fluxes by looking at the function at the element boundaries, but those things are not part of this github issue. Someone can correct me if i'm wrong.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHubhttps://github.com/idaholab/moose/issues/12310#issuecomment-430029410, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AKJ_8nD4Gh9jEzjmCuW_wKyoBcKUzO9Lks5ulQVdgaJpZM4Xblfs.

[WilkAndy]

Closing this issue: I have built Kuzmin-Turek flux-limiting TVD advection into PorousFlow, and this is almost identical to RDG(P0P1) except it operates using nodal information rather than cell-centered information, so is more suitable in PorousFlow. So i'm no longer interested in RDG.