Enforce boundary condition on Cahn Hilliard - Navier Stokes couple

General
1

Dear FeniCs community, I hope everyone is doing well.

I have a code that couples the Cahn Hilliard and Navier Stokes equations. The way in which the initial conditions are defined is via a mapping of a periodic subdomain, as shown:

  class PeriodicBoundaryX(SubDomain):
      def inside(self, x, on_boundary):
          return near(x[0], 0.0) and on_boundary
  
      def map(self, x, y):
          y[0] = x[0] - 1.0
          y[1] = x[1]

This code works perfectly for a (1 x 1) mesh, however, whenever I try to implement the code for a different geometry (look at image), the code diverges. I tried to modify the class this way:

mesh
mesh602×565 8.63 KB

  class PeriodicBoundaryX(SubDomain):
      def inside(self, x, on_boundary):
          return DOLFIN_EPS > x[0] > (- DOLFIN_EPS) and on_boundary
  
      def map(self, x, y):
        if (x[1] < 0.8):
          y[0] = x[0] - 1
          y[1] = x[1]
        else:
          y[0] = x[0] - 0.2
          y[1] = x[1]

However, the code diverges. The PB is enforced like this on the code:

  pbc = PeriodicBoundaryX()
  P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
  P2 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
  elementCH = MixedElement([P1, P1]) 
  elementNS = MixedElement([P2, P1])  
  MEch = FunctionSpace(mesh, elementCH, constrained_domain=pbc)
  MEns = FunctionSpace(mesh, elementNS, constrained_domain=pbc)

Can anyone give me any pointers into how to define this class so that my code converges?