Add missing reference element method on entities

[SoilRos]

The DUNE grid interface requires that entities are able to give a reference element directly from the entity and without necessarily going through the geometry.

[aritorto]

jenkins build this please

[aritorto]

@SoilRos thanks for noticing the lack of this method in Entity interface :) It looks great, except for a warning that brought up the only question/comment I have

[SoilRos]

jenkins build this please

[SoilRos]

@aritorto Thanks for the feedback :) This should be ready from my side now.

[aritorto]

Thanks for the update!

[bska]

I'm a little on the fence about this. The CpGrid does not really have reference elements in the traditional sense. It is true that each cell nominally is a cube, but pinch-outs and faults mean that the number of (distinct) vertices and intersections cannot be predicted/classified ahead of time.

[SoilRos]

@bska I agree with you about geometry types, but note that this PR does not add cubes as reference elements of the CpGrid entities. This was already part of the implementation, just only accessible from the geometry. What this PR is doing is to conform to the DUNE grid interface for the existing functionality.

[bska]

What this PR is doing is to conform to the DUNE grid interface for the existing functionality.

Right, but we can't really conform to that interface because the underlying entities don't really support the requisite notions.

[SoilRos]

One step at a time ;)

You are right in the sense that a cube is not a proper reference element of the entity geometry object. And while its topology is wrong, the mapping between reference element and geometry object is correct¹. This is indeed what is used in OPM today so it makes sense to be provided to the user until a proper reference element is implemented: In this way any solver (maybe beyond OPM) that is independent of the reference element topology (e.g. FV solvers) will benefit from this grid as a DUNE grid. As said, I am not introducing anything new. Or are you proposing of removing the reference element of the Geometry type as well until this is completely solved?

Footnotes

1. EDIT: the current bilinear geometry is a proper geometry, but is not correct in the sense that it does not form a proper grid. See discussion below. ↩

[bska]

Or are you proposing of removing the reference element of the Geometry type as well until this is completely solved?

I'm not entirely sure what I'm proposing, but it's more of the latter. The corner-point grid geometry itself does not support a reference element and it never will. The local coordinate implementation you mention is an approximation. Any method that depends on a reference element will never work on a corner-point geometry.

[SoilRos]

I'm not entirely sure what I'm proposing, but it's more of the latter.

Ok, then this needs to discussed with @blattms. I got the task of testing the grid against the DUNE grid checks (See #903). Without this, I do not know how we could achieve that.

The corner-point grid geometry itself does not support a reference element and it never will.

What prevents the implementation of the 5 missing reference elements (out of the 9) that arise from collapsing vertices of a cube?

Any method that depends on a reference element will never work on a corner-point geometry.

Out of interest, what makes this (general) statement true? Is it the lack of affine mappings, or the lack of elements that can be posed on the topology of the reference element?

[bska]

The corner-point grid geometry itself does not support a reference element and it never will.

What prevents the implementation of the 5 missing reference elements (out of the 9) that arise from collapsing vertices of a cube?

For one, there's no guarantee that the cell's barycentre is inside the element and, as I mentioned earlier, there's no compile-time bound on the number sub-entities/intersections of the cell's codim-1 entities.

If you're not already familiar with these stratigraphic grids, I recommend reading the overview in the 3rd chapter of the first MRST book (open access).

[SoilRos]

@bska Thanks for that reference, that was definitely a nice read!

To be fair, I am not sure I get your point. It seems to me that we are talking about different things here. I am talking about the implementation of reference elements (in the sense of the DUNE paper), which is by definition a convex polytope. I do not see how the barycenter of an entity, or the number of runtime sub-entities comes to play here.

Indeed, my impression is that the definition of a geometric realization (in the sense of the DUNE paper, again) of the CPGrid is kind of loosely defined, on purpose and for very good reasons, but this has nothing to do with the reference elements of the deformed entities.

[bska]

My point is that even if you pretend that there is a reference element in the Dune sense, the real geometry of corner-point grids does not support a reference element. Pretty much any operation for which you'd use a reference element in a numerical method will fail in the context of corner-point geometries. That's why it really doesn't make any sense to put a lot of effort into defining those operations.

[SoilRos]

[...] the real geometry of corner-point grids does not support a reference element.

I guess this is precisely the part that I do not understand. Could you please extend why? The MRST book that you provided says the following:

Since each face of a grid cell is specified by four (arbitrary) points, the cell interfaces in the
grid will generally be bilinear and possibly strongly curved surfaces. Geometrically, this
can lead to several complications. Cell faces on different sides of a fault may intersect each
other so that cells overlap volumetrically, or they may leave void spaces in between them.
There may be tiny overlap areas between cell faces on different sides a fault, and so on.
All these factors contribute to make fault geometries hard to interpret in a geometrically
consistent way: a subdivision into triangles is, for instance, not unique. Likewise, top and
bottom surfaces may intersect for highly curved cells with high aspect ratios, cell centroids
may be outside the cell volume, etc.

Which basically admits that the proposed interpolation for the faces (bilinear between corners of the entity, or triangular subdivision) does not provide a valid geometric realization, hence, not a proper grid. But this does not mean that there is not a suitable geometric realization of the grid. From what I can see for example, making a bi- or tri-linear interpolation based on the end points of the vertical pillars¹ and not on the corners of the entity will yield a mapping that does not have such gaps and overlaps between geometries. Perhaps am I wrong here? Or, is there perhaps something beyond these gaps that I do not know?

Footnotes

1. Assuming that the pillars are straight. ↩

[atgeirr]

From what I can see for example, making a bi- or tri-linear interpolation based on the end points of the vertical pillars

This is true, and it is possible to find such a representation also if the pillars are non-vertical. We observed this many years ago, and used this in our work on streamline tracing (where obviously you need a valid embedding of the grid into space...) but did not publish on that specifically. A researcher from NR came across the same problem in a different setting and published this paper: https://link.springer.com/article/10.1007/s10596-015-9500-0, I believe that would be the go-to source for an implementation taking into account these problems.

However: For compatibility with E100, we cannot really adopt this for OPM Flow. For compatibility we need to agree on the number of active cells in grids, and with features like MINPV that lets us deactivate cells with less than a certain pore volume, it becomes important to calculate the volume in identical fashion. The current volume calculation is therefore using the trilinear mapping, to get identical pore volumes as in E100.

[SoilRos]

Yes, I also meant non-vertical pillars. I was trying to not be very general because I have no idea how this would look like in curved pillars. But yes, very interesting paper, that's exactly what I meant!

I understand that compatibility with E100 is very very important for this project. However, wrong is wrong no matter the angle you see it, here and in E100. In this form, OPM (and E100) will never have a proper corner point grid, and that will bite us (them) sooner or later for more advanced features like local grid refinement, DG discretizations, streamline calculations, etc. To not drop compatibility: would it be perhaps possible to think about a (non-)permissive flag or tag that enables/disables such geometry? That is for example how Microsoft managed to catch up with the ISO C++ standard after having a pretty broken non-conforming compiler, which was the status quo for Windows programmers. Could it be that aside of the more advanced features this ends up be more beneficial in other measures too? (e.g., less stiff problem, fewer hard to catch errors due to overlaps/voids, etc).

[atgeirr]

I understand that compatibility with E100 is very very important for this project. However, wrong is wrong no matter the angle you see it, here and in E100. In this form, OPM (and E100) will never have a proper corner point grid, and that will bite us (them) sooner or later for more advanced features like local grid refinement, DG discretizations, streamline calculations, etc. To not drop compatibility: would it be perhaps possible to think about a (non-)permissive flag or tag that enables/disables such geometry?

I have considered this, but not implemented it. Note that for some things that one would usually get from the grid, such as cell depths and cell volumes, OPM Flow does not get this from the grid, to preserve E100 compatibility while letting the CpGrid class have more reasonable geometry. So in part what you are looking for is already there. We have done DG in the past with CpGrid for example, quite successfully.

Another example: explicit NNCs can be added to connect arbitrary grid cells. This is obviously not possible to realize as a proper embedding of the grid in space. We have chosen to put these connections in the grid, as faces/Intersections, to avoid having separate processing of NNCs in many places. So a grid with those is not geometrically sensible.

In the end, if you use CpGrid (without NNCs or volume/depth shenanigans) through the Dune interface, you should get something reasonable enough to do interesting numerics on, if the grid otherwise conforms to the requirements of the discretization method.

I do not think this leads to less stiff problems though: the MINPV feature for example is there specifically to make the problems less stiff, a true discretization of the complete space without MINPV would have worse condition numbers.

[SoilRos]

@atgeirr thanks for the detailed explanation.

We have done DG in the past with CpGrid for example, quite successfully.

Amazing! I was under the impression that one needed a proper geometric realization for this. Is there a paper, or a source where I can learn more about this?

I do not think this leads to less stiff problems though: the MINPV feature for example is there specifically to make the problems less stiff, a true discretization of the complete space without MINPV would have worse condition numbers.

I was meaning this in the sense of having the same active cells for simulation (e.g., same MINPV criterion) but different geometry so that the discretization only differs by the cell geometries being used. But I get what you mean too, enabling/disabling MINPV is probably more important in the end.

In the end, if you use CpGrid (without NNCs or volume/depth shenanigans) through the Dune interface, you should get something reasonable enough to do interesting numerics on, if the grid otherwise conforms to the requirements of the discretization method.

So, coming back to @bska's question. I think this settles that proper geometries and reference elements exist for CpGrids, they are just not entirely implemented yet. But his question is still valid: Is it worth having the referenceElement (and for the same effect obj.type() -> Dune::GeometryType) interface on Entity and Geometry even if we know that they are technically wrong?

[atgeirr]

Amazing! I was under the impression that one needed a proper geometric realization for this. Is there a paper, or a source where I can learn more about this?

https://andreas.folk.ntnu.no/papers/diag-ecmor-xiv.pdf