Periodic Boundary Conditions for Annulus

Legacy DOLFIN
1

Hi everyone,

I have a quick question regarding periodic boundary conditions in (legacy) FEniCS. A MWE for the problem is sketched below. Assume that I have some annulus and I want to simulate only a quarter. Then I would have to apply periodic boundary conditions at the orange boundaries.

image
image1628×842 41.1 KB

For a scalar quantity, I think this is straightforward to do. However, I am interested in simulating fluid flow, the simplest example being Stokes flow. For this, one would also need to “rotate” the direction vectors - what exits at the top to the left has to come into the domain from below in an upward direction on the right side. Can this be handled with FEniCS? If so, I am really happy to provide a MWE for the scalar problem to get this working in the vector case as well.

Thanks a lot in advance and best regards,
Sebastian

2

I do not think this is straightforward to do in legacy FEniCS, as the construct for periodicity only considers dof to dof equality (with no option of rotation), as it eliminates dofs for the dofmap.

This kind of constraint should be possible on DOLFINx using DOLFINx_MPC (disclaimer I am the author of DOLFINx MPC), as it allows for u_i = \sum_{j=0, j\neq}^N c_ju_j constraints on the DOF level.

3

You could also transfer to a cylindrical or spherical coordinate system. Then you don’t need to worry about adjusting for rotations.

4

So in order to do this with Legacy FEniCS, I could assemble the system matrix and rhs myself and then would have to modify the rows of the “slave” side of the system to incorporate the rotation, something like

u_x_slave = 0
u_y_slave = - u_x_master

for my example (so that the system matrix changes and the rhs for the corresponding row is 0). I guess this is also what fenics does for the “traditional” periodic boundary conditions.

Is there something to worry about the sparsity pattern for the matrix? Like fenics not allocating space as the mesh nodes are not neighbors?

5

Yeah, for the simple problem this works, however my application has a complex three dimensional geometry which cannot be adequately transformed to cylindrical coordinates.

6

Yes, a whole ton of things that will be problematic in parallel.

One is allocation of sparisty pattern, second is identitication of master/puppet pairs across processes.

In the case of a pure equality, one can simply change the dofmap. However, when you have coefficients, such as negation, some care has to be taken. See for instance slide 3 of:

The further slides might also be useful.

7

Thanks a lot. So it’s doable, but might involve some effort, especially when doing this in parallel.

One further question: In the slides you have linked, you use some prolongation matrix to incorporate these multi-point constraints. Is there any other benefit of using this approach rather than modifying the rows of the matrix directly other than that it can preserve the symmetry of the matrix?

8

Note that I never form a full prolongation matrix.

The method computes the local contribution Ae and then computes K^TAeK (where only rows of K that are relevant is used in the product). This forms a local part of Ae and a nonlocal part, which in turn is inserted into the global matrix

1 Like