Allowing separated InterfaceKernels using periodic boundaries

[ttruster]

Reason

Meshes written by --mesh-only that contain split blocks from BreakMeshByBlock cannot be utilized for physical analyses with InterfaceKernels in a subsequent input file which reads the mesh, due to the missing element-neighbor connectivity information. This issue impacts the ability of cohesive zone modeling for distributed meshes. Also, there is not a mechanism for inputting element-neighbor pairs other than the automatically detected pairs due to adjacency (shared nodes).

Design

Generalize Periodic Boundary Condition (PBC) input to allow weak and strong enforcement. Strong enforcement uses DOF constraint equations (current approach); weak enforcement would use surface integrals as part of the weak form of the boundary value problem.
Provide the names of the two topologically conforming (potentially separated) boundaries for weak enforcement by using the AddPeriodicBCAction, so that the pair is registered in the DOF map but does not get a constraint matrix.
Add method onPeriodicBoundary within thread loops of assembly with analogous behavior to onInterface with handles boundary ID of an element that lack an adjacent neighbor but do have a separated neighbor.
Utilize existing topological_neighbor and inverse_map methods to handle reinitialization and existing InterfaceKernels for physics.

Impact

Allows PBC class to help with solving interface problems like cohesive zone modeling and to treat periodic conforming meshes using interface kernels instead of nodal constraints, allowing Nitsche/DG/mortar/penalty treatments.

References

Originated from Discussion #22948 (reply in thread)
Requires enhancements of Libmesh from PR libMesh/libmesh#3525

[ttruster]

Reposting some text to start the discussion.

From working on interface methods for a while, my opinion is that there are three somewhat well-delineated categories:

1. Continuous/conforming meshes and no interface; using multi-point constraints to tie blocks for separate parts or periodic conditions
2. Conforming or adaptively refined meshes with hanging nodes that have Penalty/Nitsche/Lagrange-multiplier enforcement of weak continuity along the interface. Here, the corner nodes on the coarse mesh are conforming everywhere, even on external boundaries.
3. Nonconforming interfaces, mesh tying, and contact/friction modeling as well as embedded mesh methods using Penalty/Nitsche/Mortar enforcement of the interface conditions.

These 3 categories increase in generality of problem class from top to bottom while also increasing in implementation/architecture complexity and projection errors. I've implemented all 3 categories before in my matlab FEM codes and in FEAP. Within MOOSE, I see these 3 mapping onto Kernels, InterfaceKernels, and MortarConstraints respectively. So, the middle category is a compromise in difficulty and generality. The difficulty is moderate since the quadrature point location on each side of the interface is identical once an offset vector (potentially zero) is introduced; I assume that idea works even for adaptively refined meshes and the InterfaceKernels for those. This means that no extra projection or copying of stateful materials from quadrature point locations on opposite sides is needed. So, my suggested restriction is to say that non-adjacent element-neighbor concept would only be permitted for two parent element faces that geometrically conform after a constant offset vector. This would allow usage of existing InterfaceKernel classes to solve problems with "gaps" between the two surfaces, for example finite-thickness cohesive zone (CZM) interfaces in solid mechanics.

I finally got some time to work on this within the past week. The source code is here https://github.com/Truster-Research-Group/moose/tree/topological-neighbor-PBC and https://github.com/ttruster/libmesh/tree/weak-enforcement-PBC. Perhaps we can open an issue to continue the discussion of this? I had also heard that the discussion of stateful mortar constraints was continuing, so I noticed a new issue was started here for that #24030.

Some of the discussions to occur in the issue include

• Error checks to verify which variables are periodic and that they are attached to the interface kernel
• A test of whether any PBC exist in the model, to test at the top of onPeriodicBoundary
• Devise a caching scheme in MOOSE (seems better than libmesh) to list the order that the PBC pairs will be integrated, so that the results from topological_neighbor can be stored instead of recomputed at each nonlinear iteration. Perhaps a sorted table based on element and side could be used, so that a binary search could be used?
• What happens if the two sides of an "interface" live on two processors? I suppose the algebraic ghosting is not too difficult to handle with the relationship managers, but how about other data structures such as the element-neighbor pairs and their offsets?
• How would you differentiate a standard interface kernel from the separated interface kernels? The former requires only a single boundary argument whereas the latter conceptually requires two boundaries right?
• Do we want to strictly enforce that all meshes resulting from BreakMeshByBlock be coupled with a PBC object? And if so, then should we improve the automation of inputting PBC pairs (or rather, can BreakMeshByBlock provide a list of all the boundary ID/names that it generated?)
• any other concerns with my crude implementation

[lindsayad]

@ttruster this would resolve #21161 and #21154, wouldn't it? And then consequently make the hacks we have to do for BreakMeshByBlockGenerator in #24018 unnecessary?

[ttruster]

@lindsayad by glancing over these, I believe this approach would (have a good chance to) fix all three of these.

From #21161 I agree that traditional theory of CZM is that the two material points that "talk to each other" (have force interaction) remain fixed until complete fracture/separation; this is the difference from contact mechanics with friction where the interacting points changes during sliding. The approach I'm proposing here, using the terminology of #21161, is to use the existing periodic boundaries system with a "weak-enforcement" flag to serve as a data-container or at least a means to look up the "fake neighbor" (after breaking the mesh or just from two pairs of boundaries that were never broken).

I haven't dealt with distributed meshes yet, but I would hope that by removing the "hack" and using "supported" functions from libmesh such as topological_neighbor, that the Ghosting system could then be used properly.

Whether the ex-neighbor information is stored when the mesh is broken in BreakMeshByBlock (which that MeshGenerator should probably be finalized now to call nullify neighbors or whatever the correct API is to decouple them) or if it is stored after the first pass through the InterfaceKernel assembly loop (actually, the material property initialization phase would be the most logical place probably; its the first call of the threaded loops, I assume), it would save on run-time to not have to re-compute the neighbor at each iteration since it will never change. But, we'd need a good idea for where to store that (since the neighbor pointers on Elem * are out of the question because of definition of being adjacent); instead of an unordered list of element-neighbor pairs, could it be sorted by the element-face ID to allow a binary search, or could it be stored in the same order that thread loops will pass over the element-faces?

Because that's the drawback of this: while it makes the handling of CZM more robust, it also makes it "slower". Though I have no idea what fraction of time this is compared to solving Kd=F or computing K and F; maybe it doesn't "matter"...

[ttruster]

I also forgot, I'm including in here the ability of the two sides to not be physically adjacent; the two boundaries just need to be offset by a constant vector. That's why to me this fits "logically" with periodic boundary conditions and why it motivates doing something "different", i.e. not just forcing the CZM stuff to fit within the existing InterfaceKernel (IK) system.

In short, the difference is:

1. Classical IK are for coupling two different variables along an internal boundary/interface (single boundary), such that the total solution field is discontinuously interpolated along the interface but has weak continuity or other condition imposed
2. Proposed IK is for coupling the same variable (or potentially two different ones) along two boundaries that topologically conform, such that the total solution field is discontinuously interpolated along the interface but has weak continuity or other condition imposed

Because the "result" (second phrase) is the same for both case, it's logical to me to use an InterfaceKernel object to perform the physical integration in both cases. The weak form terms appear identical for both cases; the issue is the discrete representation of the two cases is different.

[ttruster]

One more idea, from looking over those 90+ comments at #12033, and based on my list of bullets: these are a "little" more complicated to describe in the input file because the user has to "know" the list of paired boundaries to mention in the Periodic BC (PBC) object. And then they need to add an InterfaceKernel (IK) object to the "element" side boundary.

Maybe we can write a composite action that writes both the PBC and the IK. Then the boundaries are only mentioned once. The downside is, how to pass all the IK parameter values since that action needs to potentially handle all possible types of IK. So we are probably stuck with just needing to add a list of the PBC using the existing AddPeriodicBCAction.

Anyways, for my physical problems, I'll be generating my .i files using Matlab scripts from my code https://bitbucket.org/trusterresearchgroup/deiprogram which identifies the element-neighbor pairs both for adjacent and for separated interfaces in an originally contiguous periodic model.

[lindsayad]

I think an action makes sense. We could leverage the existing cohesive zone action

[ttruster]

I can work on modifying the existing master action to create both the Interface Kernel and the periodic boundary.

I will also try out Distributed Meshes on my own.

Are we going to enforce that, all meshes generated with the BreakMeshByBlock will need a periodic boundary in order to assign an Interface Kernel? That will mean that all current CZM input files need to be revised to include two boundary names. This means that backward compatibility is lost; but since the old way wasn't fully robust, then is that an acceptable change?

[lindsayad]

is that an acceptable change?

Yes, absolutely. Allowing the application of BreakMeshByBlockGenerator to CZM is one of my biggest regrets in MOOSE. I would do anything to have it be more robust.

[lindsayad]

You will be a hero to me if I can remove the hacks in #24018

[ttruster]

@lindsayad and @roystgnr: Well I have some interesting information to report about the current version of the source code with these weak periodic boundaries. I'll have to be very specific, because I will need your help to know what in the rest of the code will need to change. If you want to see the files I can push them into my branch.

1. I verified that my current executables pass the tests in https://github.com/idaholab/moose/blob/next/test/tests/mesh/splitting/tests.
2. I made a version of bilinear_mixed.i from the TM/CZM test cases that has an 8x4 mesh and uses named interfaces (block0-block1) and split that mesh into 2 processors. The reason for this: with the automatic mesh splitting feature, the split line goes VERTICAL, putting half the model in each .cpr file. And, since the CZM interface is HORIZONTAL, that means that the element-neighbor pairs in both distributed meshes are in the file; I verified that by plotting the mesh and looking at the node numbering (no case where neighbor and element from a pair are in different meshes; for any mesh where the pairs are separated, then I get errors when running the analysis).
3. The model from BreakMeshByBlockGenerator runs on distributed mesh both from t=0 and from recover; however, the computed physics is wrong because the opposing side is not found so the Interface Kernel doesn't compute anything (i.e. the disp_x variable is zero for all time)
4. The model using this proposed periodic boundary object fails with a libmesh error message "... src/mesh/mesh_base.C, line 1458 ..." using tensor-mechanics-debug and starting at t=0
5. The model using this proposed periodic boundary object SUCCEEDS using tensor-mechanics-opt from t=0, and the computed physics is identical to the serial mesh result (from either BreakMeshByBlockGenerator or the periodic case)
6. Also, BOTH the debug and opt versions can recover from step t=2 and go to step t=4 or whatever future time and compute the right physics.

So in summary:

• the opt version doesn't run some error check that is run in the debug version, which permits the periodic boundary case to run and do its job. So I got lucky and was able to get a sample working tonight!
• for topological_neighbor to work, the split meshes must contain the element and neighbor pair. My guess is, for the cases with more splits and the discontinuous meshes from BreakMeshByBlockGenerator, then the pairs are not ghosted properly so they don't get listed in the split files although they are needed
• My hypothesis: The reason the ghosts are not written is because the AddPeriodicBCAction is called in MOOSE not in the mesh phase but in the analysis phase, so the splitting algorithm doesn't notice the required ghosting that the periodic boundary would be requesting.

Thoughts on how to get those ghosts written out?

[ttruster]

Before I go to bed, to un-informed ideas to get the ghosts written

1. have BreakMeshByBlockGenerator do something that makes a ghost listing for whatever it breaks up
2. store the periodic boundary data somewhere in the MOOSE mesh object instead, so that the ghost dependencies needed for them would be registered

The second case makes sense to me, as right now I'm never checking that the periodic boundary condition is actually assigned to the variable that the Interface Kernel is acting on, I'm only using it to look up boundary pairs.
Also, I assume mesh splitting uses Metis or a similar library. Which might want to know which element-neighbors need ghosting at best to create a better partitioning or at worst to make sure all the necessary pairs are written to the separated files.

[roystgnr]

I've forgotten - what is MOOSE currently doing with PeriodicBoundary objects on a DistributedMesh? You're exactly right about the problem - DistributedMesh only knows how to delete remote elements and repartition/redistribute already-properly-ghosted elements; it doesn't yet know how to get improperly-deleted elements back if the required ghosting level is increased after the fact by something like the post facto addition of a new periodic boundary). But IIRC we hit this issue years ago, and although MOOSE was unable to make the obvious "add all periodic BCs before prepare_for_use()" fix, IIRC there was some workaround for it ... what was that workaround and why would it not be working here too?

[lindsayad]

When there are "late" geometric ghosting functors to be added like in the case of periodic boundaries, we tell the mesh to not delete any remote elements until the periodic boundaries have been added