Nodes belonging to multiple periodic boundary conditions are improperly handled

[fmonteghetti]

Hi Jørgen,

I came across this issue while investigating regression in some of my models.

Issue: Problems with multiple periodic boundary conditions may not be accurately discretized if there are nodes common to several periodic boundaries.

Reproduction steps: Run the attached demo_periodic_gep_modified.py in dokken92/dolfinx_mpc:v0.4.1 and

• compare the computed eigenvalues to the exact ones
• compare the computed eigenvectors (eigenvectors_with_dolfinx_mpc.vtu) to the ones obtained with fenics (eigenvectors_with_fenics.xdmf).

Background: demo_periodic_gep_modified.py computes the eigenvalues of the Laplacian on the unit square with two periodic boundary conditions: the left edge is mapped to the right edge and the bottom edge is mapped to the top edge. The four corner nodes (0,0), (0,1), (1,0), and (1,1) belong to two different periodic boundary conditions.

Comparing the fenics and dolfinx_mpc results shows that dolfinx_mpc is accurate only for eigenfunctions that are null at all of the corner nodes, while fenics has no such issue. This can be seen on the first eigenfunction u_00, which should be constant, as well as on u_04 and u_08.

Observations and suggestions:

While investigating this I have made two observations:

• Running demo_periodic_gep_modified.py in an older version of dolfinx_mpc [1] leads to the same problem but milder.
• All my models with two periodic boundary conditions regressed between [1] and v0.4.1. All of them have spurious behavior at corner nodes. Most of them have a different sparsity pattern in v0.4.1.

[1] dolfinx 2019.2.9.99 / dolfinx_mpc 4c9b04c4678d0f36c0a819c2892cc73c7a5f85ab.

Here are some uneducated suggestions:

• The logic is correct but there is a loose tolerance somewhere that manifests itself more sensibly in newer version of dolfinx/dolfinx_mpc.
• There is an issue in the logic that can overlook nodes belonging to multiple periodic boundary conditions.
• create_periodic_constraint_topological() is not designed to be called twice, and an alternative strategy should be used to deal with multiple periodic boundary conditions.

Attachment: issue-double-pbc.zip

• demo_periodic_gep_modified.py: modification of the demo demo_periodic_gep.py available in dokken92/dolfinx_mpc:v0.4.1 to use two periodic boundary conditions instead of one.
• eigenvectors_with_dolfinx_mpc.vtu: eigenfunctions computed by running demo_periodic_gep_modified.py in dokken92/dolfinx_mpc:v0.4.1.
• eigenvectors_with_fenics.xdmf: eigenfunctions computed using fenics on the same problem.

[jorgensd]

I am aware of this issue. The MPCs, being as general as they are, does not sanity check the constraints. i.e., if you constrain a degree of freedom twice, you need to remove one of the constraints.

As create_periodic_constraint_* simply appends data to a numpy array: https://github.com/jorgensd/dolfinx_mpc/blob/master/python/dolfinx_mpc/multipointconstraint.py#L162

An additional check could be added here that checks if the input slaves are already in another constraint. However, this is a slippery slope, as the next sane step would be to check that no slave degree of freedom is a master in another constraint, i.e.

c_i = b c_j
c_j = d c_k

which should be simplified to

c_i = c_k
c_j = c_k 

If one wants this kind of nested behaviour, the setup-time for the assemblers would grow with the number of constraints.

To summarize:

1. The current MPCs are limited to one equation per DOF, i.e. x_i= a x_j and x_i = b x_k is not supported
2. Nested constraints: x_i = a x_j, x_j = b x_k is not supported.

A way to work around this is by creating a single periodic constraint for both of the edges in your problem, and create a mapping that depends on the geometrical coordinates. I.e. for a square: (mapping left ->right), (bottom to top)

def mapping(x):
     left_constraint = x[1] > 1e-14
     mapped_values = np.zeros((2, x.shape[1]))
     mapped_values[0] = left_constraint + x[0]
     mapped_values[1] = np.invert(left_constraint) +  x[1]
     return mapped_values

This means that if the coordinates of the dofs on the facet is on the left side (excluding top right corner) they will be mapped from (0, y) to (1, y) and (x, 0) to (x, 1).

[fmonteghetti]

Thanks a lot for your answer, I will try to make this work on the unit square demo and report back. It sounds like it would be worth turning this into a demo, since studies of periodic media commonly encounter this case.

I typically work with models where the slave and master boundaries are identified with tags provided by gmsh, rather than geometrical equations. So, in principle, I know directly the slave and master dofs that are to be mapped by each periodic boundary condition, as well as the dofs that are over-constrained.

Is there a way to use this information to avoid embedding geometrical tolerances in the definition of the periodic boundary condition?

For example, a lower-level way could be to provide dolfinx_mpc with the slave and master dofs directly. It would then be the responsibility of the user to ensure that no dofs appear twice, but this is easy to check. Is this possible right now?

[jorgensd]

Thanks a lot for your answer, I will try to make this work on the unit square demo and report back. It sounds like it would be worth turning this into a demo, since studies of periodic media commonly encounter this case.

I typically work with models where the slave and master boundaries are identified with tags provided by gmsh, rather than geometrical equations. So, in principle, I know directly the slave and master dofs that are to be mapped by each periodic boundary condition, as well as the dofs that are over-constrained.

Is there a way to use this information to avoid embedding geometrical tolerances in the definition of the periodic boundary condition?

For example, a lower-level way could be to provide dolfinx_mpc with the slave and master dofs directly. It would then be the responsibility of the user to ensure that no dofs appear twice, but this is easy to check. Is this possible right now?

The lower level way of adding a constraint is: https://github.com/jorgensd/dolfinx_mpc/blob/8abd427d6d909cc2883d4e45ee04cb080db2fe40/python/dolfinx_mpc/multipointconstraint.py#L39-L60 The complicated bit is to get the global number of the master dof, as this might live on another process that then slave dof. What the wrappers do is to do this in a scalable way, with MPI neighborhood communicators.

[fmonteghetti]

Thanks a lot for you help, I really appreciate it. I need to read your MPI magic! :)

The suggested workaround works: in demo_periodic_gep_modified.py, we can keep two calls to create_periodic_constraint_topological() as long as the slave_to_master_map of the second periodic boundary condition excludes dofs that belong to either the slave or master boundary of the previous periodic boundary condition. I will send a modified demo when I have a clean implementation.

By the way, is there a clean way to test x[0]>1e-14? In the fenics days I relied on fenics.near(x[0],0).

For future readers, here is a description of the workaround.

--- Workaround to enforce N periodic boundary conditions Currently, a user calling create_periodic_constraint_topological() several times can end up with a wrong computation without any warnings or errors. One way to apply N periodic boundary conditions reliably is to define:

1. N functions I_slave[n]: x -> bool that identifies the n-th slave boundary.
2. N functions I_master[n]: x -> bool that identifies the n-th master boundary.
3. N functions slave_to_master[n]: x -> x that maps the n-th slave boundary to the n-th master boundary.

The user must ensure that slave_to_master[n] excludes dofs that are master or slaves of one of the n-1 previously applied periodic boundary conditions (i.e. dofs that satisfy I_slave[k][x]==True or I_master[k]==True for k<n). This can be done by sending excluded dofs outside the mesh (there may be a cleaner way). ---

The workaround is consistent with fenics usage, where one must provide (1,2,3). However, in fenics the user does not need to manually exclude slave or master dofs of the previously applied boundary conditions.

My understanding is that if create_periodic_constraint_topological() is modified to accept a way of identifying the master boundary, then all the information would be there to safely apply N periodic boundary conditions within dolphin_mpc. But changing the API is clearly not ideal.

[fmonteghetti]

demo_periodic_gep_modified.zip

Here is a modified demo_periodic_gep.py that demonstrates the workaround described above. It produces the correct eigenvalues (including multiplicities) in dokken92/dolfinx_mpc:v0.4.1. The two modifications are:

• Slave/master dofs are excluded just before the call to create_periodic_constraint_topological() by re-defining pbc_vert_slave_to_master_map. This syntax generalizes to N periodic boundary conditions easily (put the periodic boundary conditions in a list and loop over the previous ones).
• Dof exclusion is done by setting the first coordinate to numpy.nan.

If you agree with these modifications, then perhaps this could be added as a demo, or I could merge it with demo_periodic_gep.py.

The existing demo_periodic_gep.py tests whether diagval works, while the attached demo_periodic_gep_modified.py demonstrates how to use multiple periodic boundary conditions, which could be useful to people interested in periodic media moving from fenics to fenicsx.

[jorgensd]

Feel free to make a pull request

[jorgensd]

Demo illustrating how to handle this has been added in #22