Integration of problem with Robin and periodic boundary conditions gives unpredictable results

Question

I perform a sequence of runs of the following program either from the command line or from a shell (using Ubuntu 14.10, Fenics version 1:1.5.0.1~pp).

I randomly get either the correct (symmetric) solution, or wrong (asymmetric) solutions.

During the sequence of runs JIT is called randomly. After a JIT call I most often get the wrong solution.

Turning the resolution to res=10 helps. Higher resolutions (e.g. res=20) are more likely to produce the wrong solution.

#!/usr/bin/env python
#
#  Use finite elements to solve Poisson equation 
#  with Robin and periodic boundary conditions 
#
#


from dolfin import *           
from mshr import *        


#  Define domain

side     =  2.
radius   =  0.35
res      =  20                   #  number of nodes along boundary side     

radius2  =  radius*radius
halfside =  0.5 * side
tol      =  side/res/100.        #  tolerance for coordinate comparisons

domain   =  Rectangle(dolfin.Point(-halfside, -halfside),  
                      dolfin.Point(halfside,halfside))     -   \
            Circle(dolfin.Point(0.0, 0.0), .35)


#  Define source and fluxes

q = Constant(  0)                #  for the moment, no production on the border
p = Constant( 10)                #  degradation on the border
f = Constant( .1)                # production in the bulk


#  Define boundaries

class InnerCircle(SubDomain):
   def inside(self, x, on_boundary):
      return on_boundary and (x[0]*x[0]+x[1]*x[1] - radius2 < tol)


rbc = InnerCircle()


class SquareBorder(SubDomain):

   def inside(self, x, on_boundary):
      # return True if on left or bottom boundary 
      # AND NOT on one of the two corners (-L/2, L/2) and (L/2, -L/2)
      return (on_boundary and 
              ( (abs(x[0]+halfside)<tol) or  (abs(x[1]+halfside)<tol) ) and not
            ( ( (abs(x[0]+halfside)<tol) and (abs(x[1]-halfside)<tol) ) or 
              ( (abs(x[0]-halfside)<tol) and (abs(x[1]+halfside)<tol) )) ) 

   def map(self, x, xn):
      # identify points on the opposite sides
      if ( (abs(x[0]-halfside)<tol) and (abs(x[1]-halfside)<tol) ):
         xn[0] = x[0] - side
         xn[1] = x[1] - side
      elif ( abs(x[0]-halfside)<tol ):
         xn[0] = x[0] - side
         xn[1] = x[1]
      else:   # abs(x[1]-halfside)<tol )
         xn[0] = x[0]
         xn[1] = x[1] - side


pbc = SquareBorder()


#  Create mesh and define function space

mesh=generate_mesh(domain,res)
V = FunctionSpace(mesh, "Lagrange", 1, constrained_domain=pbc)


# Define boundary segments for Robin and periodic boundary conditions

boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim()-1)

rbc.mark(boundary_parts, 0)
pbc.mark(boundary_parts, 1)


# Define variational problem

u = TrialFunction(V)
v = TestFunction(V)
dss=ds(domain=mesh, subdomain_data=boundary_parts)
a = inner(nabla_grad(u), nabla_grad(v))*dx + p*u*v*dss(0)
l = f*v*dx + p*q*v*dss(0)           


# Assemble conditions

A = assemble(a)
L = assemble(l)


# Solve A U = b 
#
#  lu    :  LU factorization (Gaussian elimination)
#  ilu   :  incomplete LU factorization (preconditioner)
#  cg    :  conjugate gradient (Krylov solver)

set_log_level(PROGRESS)
u = Function(V)
U = u.vector()
#solve(A,U,L)    
solve(A,U,L,'lu')    
#solve(A,U,L,'cg','ilu') 


#  Compute nodal values

u_ar = u.vector().array()
coor = mesh.coordinates()


# Plot solution and mesh

offset=-u_ar.min()+0.5*(u_ar.max()-u_ar.min())
plot(u+offset)
interactive()

Comments

Hi, could you post an image of the wrong solution?

---

Sorry for not having answered immediately.
Here are typical examples of the "wrong" and "right" solution.
They keep coming out completely at random, I am puzzled.

Answer 1

Hi, the following snippet illustrates that the mesh generation with mshr is not deterministic, which explains why you randomly get 'correct/incorrect' solution

from dolfin import *
from mshr import *
import numpy as np

res = 10
mesh_type = 'mshr'

get_mesh = lambda domain, res:\
           generate_mesh(domain, res) if mesh_type == 'mshr' else\
          UnitSquareMesh(20, 20)

domain = Rectangle(Point(0, 0), Point(1, 1)) - Circle(Point(0.5, 0.5), 0.125)

x_ = get_mesh(domain,res).coordinates().reshape((-1, 2))

n_mesh_changes = 0
for i in range(100):
    x = get_mesh(domain,res).coordinates().reshape((-1, 2))
    if np.linalg.norm(x - x_) > 1E-14:
        n_mesh_changes += 1
    np.copyto(x_, x)

print 'Mesh changed', n_mesh_changes 

That said, as far as I know, periodic boundary conditions work best if the mesh in the periodic directions is such that the target boundary exactly mirrors the periodic boundary. In your case this is not true, as some of the edges lack a periodic counterpart. This is demonstrated here

class PeriodicBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] < DOLFIN_EPS

    def map(self, x, y):
        y[0] = x[0] - 1.0
        y[1] = x[1]

pbc = PeriodicBoundaryComputation()
mesh = get_mesh(domain,res)
mesh.init(1)

paired_edges = pbc.compute_periodic_pairs(mesh, PeriodicBoundary(), 1)
edge_f = FacetFunction('size_t', mesh, 0)

for i, (edge0, edge1) in enumerate(paired_edges.iteritems(), 1):
    edge_f[edge0] = i
    edge_f[edge1[1]] = i

plot(edge_f) 
interactive()

Comments

Thank you Mirok,

what you say seems quite likely! I couldn't think of it.

But then:

• is there a canonical way to generate a mesh that is compatible with periodic boundary conditions?

• is there a way to instruct mshr to create every time the same mesh?