Read boundary faces from mesh and set boundary condition in a demo

[jeff-cohere]

Instead of developing a separate demo to illustrate the tagging of boundary faces as in #182, @bishtgautam and I decided to modify a few of our existing demos to be able to read HDF5 mesh files and figure out boundary faces from labels within. Then the demos can assign boundary conditions as instructed by the user.

Specifically:

1. If the user supplies an HDF5 mesh file with labels indicating boundary faces, those boundary faces are read into a DMLabel (either automatically by DMPlex or by using HDF5 logic--I'm not sure how much PETSc does for us in this case).
2. The DMLabel (either explicitly or from within DMPlex) is used to identify boundary faces to which boundary condition is assigned.
3. The user provides a boundary condition to be assigned to the boundary faces.

We've been using the terms "Dirichlet" and "Neumann" for boundary conditions, but it seems to me that we might benefit from discussing "pressure", "flow/no-flow", and "seepage" boundary conditions instead. As @nocollier has pointed out, the FE and FV methods actually express these boundary conditions differently in terms of the underlying Dirichlet and Neumann constraints because of differences in the math.

Initially, we intend to support only a single boundary condition for the entire domain boundary. Later, we can enhance this as needed, perhaps allowing different boundary conditions on different "domain boundaries" (the "south" boundary condition vs the "east" boundary condition, say).

I'm going to try this out on the transient_snes_mpaof90 demo first to see how it looks.

[jedbrown]

Boundaries should be labeled in the HDF5 file (and we can expand the schema if it's insufficient) and the PETSc reader can populate the "Face Sets" label. I'm more familiar with details here when using Gmsh and to some extent, Exodus. I have this on my plate to add, and I think it'll be useful for us and (especially) those applications with lots of different boundary conditions on different surfaces.

https://gitlab.com/petsc/petsc/-/issues/915

[jeff-cohere]

Awesome. Thanks, Jed.

[jeff-cohere]

I'm making a survey of our demos, trying to figure out how we enforce boundary conditions. It looks like a few things are going on:

• The two-point-flux method loops over faces with a support size of 1
• The BDM method checks the "marker" label to see whether a given face falls on the boundary
• The MPFA-O method (and whoever else makes use of the TDyMesh logic) interprets faces with a support size != 2 as boundary faces and marks them appropriately.
• It looks like the boundary conditions in the WY code are commented out.

Does this sound right?

We've been talking about setting up or just using a DMLabel that marks the set of boundary faces, which is close to what the BDM logic does. In particular, we're talking about starting with a single label containing all boundary faces in an undifferentiated way--the "domain boundary". If this is right, I can modify the code above to use that label. I assume the label can be created with a call to DMPlexMarkBoundaryFaces. This is nice and simple, as it works both for meshes read from files and generated on the fly.

Anyone have a good name for our label? "boundary"?

@bishtgautam @jedbrown @knepley

[bishtgautam]

@nocollier ⬆️

[knepley]

I like "boundary" and starting using it after "marker. We can start with everything marked as 1, and then if we need differentiation, we can change the value on some points.

[nocollier]

I think that the WY and BDM implementations will end up getting removed so I wouldn't worry terribly about how the boundary conditions are implemented. If we later want the quad/hex BDM/BDDF spaces, they are now available in PETSc and so we can push much of that implementation into PETSc and use whatever they do for boundary conditions. The WY code does not have an analog in PETSc that I know of and at the moment is a 'stunt' with boundary conditions hard coded. I could work to make it more general, but am not highly motivated to do so as it doesn't solve our problem.

[jeff-cohere]

@nocollier, what do you envision for how the WY and BDM implementations are represented in TDycore? I suppose this gets back to our collective vision about whether we want a bunch of "demos" or a single "driver" that offers the functionality of all "modes" and "discretizations." I don't know whether it's too early to start talking about such things again, but I'm interested to hear what people think of this topic at this point.

[nocollier]

In my mind, TDycore is a library from which we can write various types of subsurface flow problem implementations--linear, nonlinear, time-dependent non-linear. I think that this is important to being usable by earth system models. I don't have anything against a single monolithic driver, but if we push too much into the driver, then it becomes harder to be used by E3SM or other earth system models. Ideally all 'modes' and features would work for all discretization techniques. But as the project has progressed, some of the methods we explored we have learned aren't going to work for the kinds of problems/meshes we need to run on. So in my mind they are not something that I think is worthwhile to pursue and develop at this time. I think that it is fine for those methods to sit in the code-base but if you try to do something that isn't implemented for that method, then it just throws a error. For example, if you try to solve a transient problem with WY, it should error, or set boundary conditions using the function calls.