Constraint to enforce matching of gradients across boundaries

[dschwen]

I need a construct to enforce the "periodicity" of the gradients of a non-linear variable without enforcing the periodicity of the values. Applications are for physics where the value does not appear in the model, but only the gradient does. Like for example the electric field (as gradient of a potential) or the strain (as gradient of the displacements).

I was told, the periodic boundary system should be helpful (to provide the matching nodes/sides) but that such a constraint was not really a boundary condition, but rather would have to be implemented using a lagrange multiplier (using a boundary constrained face variable).

Is the periodic BC module in libMesh flexible enough to be (ab)used to extract the matching nodes? Is libMesh/MOOSE even ready to implement this sort of constraint?

@roystgnr , @jwpeterson , @benkirk, @dknez , @permcody , @andrsd

P.S.: we had a somewhat related request (periodicity with offsets) here: https://github.com/idaholab/moose/issues/1345

[roystgnr]

I assume you're using C1 elements? Otherwise we don't even enforce continuity of gradients between elements when the elements aren't on opposite sides of the mesh.

The PeriodicBoundary is not currently capable of doing what you want.

It might not be too hard a feature to add, though. Let me look.

See the tests for C_ZERO in compute_periodic_constraints() in fe_base.C? Instead of testing the continuity of the variable type, we could keep track of the continuity type to enforce as a PeriodicBoundary member variable, which would default to the continuity of the variable type but which could be overridden by the user in the C1 case.

The obvious trouble is that, if you tried to enforce only-value or only-gradient periodicity on a C1 element, the cholesky_solve() will fail due to the singular matrix. We'd need a solver that could not only give us a solution to the underconstrained projection problem but also let us know which dofs were "free" and which would need constraint equations. In the periodic-gradients case, we'd simultaneously have to cope with the diagonal problem (where neither gradient dof is individually constrained, a linear combination is) and the corner problem (on corner vertices, you may not get to pick which dof gets assigned your constraint equation because there may be conflicts with existing equations). Getting that all working correctly (including on ParallelMesh) was nasty enough even in the simple all-side-dofs-are-constrained cases.

[dschwen]

I assume you're using C1 elements?

Well, color me naive! I assumed that if I can obtain a gradient (for example on first order lagrange) I can just whip up a constraint such as |grad(u_l-u_r)| and force it to go to zero using a lagrange multiplier (or whatever works).

[roystgnr]

You can, actually. You just don't need nearly as much libMesh support to do so. ;-) As a very first crack you can add "(|grad_u(x)-grad_u(opposite(x))|,v)/epsilon" to residuals on both boundaries. It won't give you the correct answer (not even in the weak sense that a lagrange multiplier would; you'll converge to an O(epsilon) error), but it'll be as good as the "penalty boundary condition" hack we used to use for all our Dirichlet BCs in the years before DirichletBoundary.

I'm still a little "iffy" on what that really means if you don't have a problem with solutions in H2/C1, though. It sounds overconstrained in that case to me.

In any case, now I think I understand your question better. You want to enforce the periodic gradient continuity weakly, but you need libMesh support for finding opposite(x), making sure the element containing it is semilocal on a ParallelMesh, enforcing level-one rules during AMR, etc? Creating a PeriodicBoundary with no periodic variables (which will require some library changes; currently the internals assume "if the set of periodic vars is empty, all vars are periodic" to enable users to not micromanage that common case) would do all that, I think - all the code checking for a topological_neighbor() across a PBC would then find the right neighbor, and your physics code could query get_corresponding_pos() from the PeriodicBoundary, but there wouldn't be any constraints generated.

[dschwen]

Hmm, I guess if my elements are not C1 then the gradient is not defined at nodes and edges anyways, but for those elements I could get the gradient in the interiors of the neighboring elements across the boundary. That does sound like a clean definition for me. But then again what do you do for a nodal BC?! Ok, I guess I need to think about this a bit harder...

[jwpeterson]

Sounds like more of a math issue than a code issue.

[dschwen]

The code issue is setting up the pbc without any variables.

On Mon, Oct 12, 2015, 3:22 PM John W. Peterson notifications@github.com wrote:

Closed #545 https://github.com/libMesh/libmesh/issues/545.

— Reply to this email directly or view it on GitHub https://github.com/libMesh/libmesh/issues/545#event-433381475.

[dschwen]

I reopened this because it is a code issue. As Roy said the PBC code needs to be adapted to allow the part of the PBC setup that finds the opposite nodes and ghosts the DOFs without actually applying the PBC for a variable.

[dschwen]

I think there may be a simpler way to get what I want in therms of representative volume element (RVE) periodicity. It is described in this paper:

http://www.mate.tue.nl/mate/pdfs/167.pdf

The gist is that we have a simulation cell with a bunch of nodes

1---2---3---4
|   |   |   |
*---*---*---*
|   |   |   |
*---*---*---*
|   |   |   |
5---6---7---8

Now rather than having periodicity as in u1=u5, u2=u6, ..., u4=u8 we move each side into a frame of reference relative to its first node! So the BC is

u2-u1=u6-u5, u3-u1=u7-u5, u4-u1=u8-u5 (the first term drops out as 0=0, which results in the extra unconstrained degrees of freedom).

This formulation retains the shape of the opposing sides, but allows them to move relative to each other.

What makes this attractive is that this BC is structurally similar to the regular Periodic BC (no gradients).

[roystgnr]

Oh, I was totally misunderstanding what you wanted, then. I thought you wanted all components of the gradient to be periodic, not just the tangential component. I'm still not sure under what conditions you'll have a properly constrained BVP, but at least now the question is "won't the operator have a 1-D kernel?" rather than "are you even in the right Hilbert space?"

For the moment I think the fastest way to get up and running is going to be weak enforcement. In the long term, I've wanted to add the ability to do a user-defined "transformation" in variable space, not just in physical space. My goal was to enable periodicity in polar coordinates of cartesian velocity vectors, but the same capability would allow you to add a single SCALAR variable and use it to represent the jump from one side to the other.

[dschwen]

I want all components to be periodic. But my last post outlines an alternative approach that does not involve gradients.

On Thu, Apr 21, 2016, 7:46 PM roystgnr notifications@github.com wrote:

Oh, I was totally misunderstanding what you wanted, then. I thought you wanted all components of the gradient to be periodic, not just the tangential component. I'm still not sure under what conditions you'll have a properly constrained BVP, but at least now the question is "won't the operator have a 1-D kernel?" rather than "are you even in the right Hilbert space?"

For the moment I think the fastest way to get up and running is going to be weak enforcement. In the long term, I've wanted to add the ability to do a user-defined "transformation" in variable space, not just in physical space. My goal was to enable periodicity in polar coordinates of cartesian velocity vectors, but the same capability would allow you to add a single SCALAR variable and use it to represent the jump from one side to the other.

— You are receiving this because you modified the open/close state. Reply to this email directly or view it on GitHub https://github.com/libMesh/libmesh/issues/545#issuecomment-213207128

[dschwen]

Oh, I think I understand your comment now. The normal component of the gradient is not readily available in the face elements, is it?

[friedmud]

Sure it is... Just dot with the outward facing normal.

Also I don't understand your "RVE" approach. It really just looks like poorly constraining the gradient. On Fri, Apr 22, 2016 at 10:01 AM Daniel Schwen notifications@github.com wrote:

Oh, I think I understand your comment now. The normal component of the gradient is not readily available in the face elements, is it?

— You are receiving this because you are subscribed to this thread.

Reply to this email directly or view it on GitHub https://github.com/libMesh/libmesh/issues/545#issuecomment-213440334

[dschwen]

Ah ok, good to know (makes sense). The RVE approach should perfectly enforce the gradient periodicity only for first order shape functions. On the plus side this should be doable by constraining dofs directly. Or am I missing something?

[friedmud]

Don't forget that it will also only really be correct for structured grids... and I can't see how it would work with adaptivity.

Basically: I don't really think it's a good idea... it's too "discrete" for my tastes. But: go nuts :-) On Fri, Apr 22, 2016 at 10:41 AM Daniel Schwen notifications@github.com wrote:

Ah ok, good to know (makes sense). The RVE approach should perfectly enforce the gradient periodicity only for first order shape functions. On the plus side this should be doable by constraining dofs directly. Or am I missing something?

— You are receiving this because you commented. Reply to this email directly or view it on GitHub https://github.com/libMesh/libmesh/issues/545#issuecomment-213453400