Divergence of Raviart-Thomas and de Rham complex

Question

Dear all,

I am setting up a solver and at a certain point I need to have a divergence free vector field. I can get a divergence free vector field (up to machine precision) if I adequately select my function spaces. I start with a stream function (psi) discretized with CG (continuous Lagrange interpolants) and compute my vector field (u) as:

u = curl psi

If u is discretized by Raviart-Thomas (RT) elements, the curl of elements of CG are mapped to elements of RT. This should be exact (up to machine precision) and therefore the divergence of my velocity field should be zero (up to machine precision).

I have tried this but I am unable to reproduce it.

I have set up a small code in which I compute the velocity from a stream function and then its divergence. I do this for different number of elements and polynomial degree of the basis functions. As you can see, the error goes up quite quickly.

Is there a problem in the implementation of RT elements?

import dolfin
import numpy
import pylab

# define an exact stream function
psi_exact_str = 'x[1]<=pi ? epsilon*cos(x[0])-(1.0/(cosh((x[1]-0.5*pi)/delta)*cosh((x[1]-0.5*pi)/delta)))/delta : epsilon*cos(x[0]) + (1.0/(cosh((1.5*pi-x[1])/delta)*cosh((1.5*pi-x[1])/delta)))/delta'   
epsilon = 0.05
delta = numpy.pi/15.0

psi_exact = dolfin.Expression(psi_exact_str,epsilon=epsilon,delta=delta)

# define the number of elements and polynomial degree of basis functions
n_elements = [5,10,20,40,80]

pp = numpy.arange(1,5)

# define the mesh
xBounds = numpy.array([0.,2.0*numpy.pi])
yBounds = numpy.array([0.,2.0*numpy.pi])

div_max = numpy.zeros([len(pp),len(n_elements)])

for kp,p in enumerate(pp):
    for kn,n in enumerate(n_elements):
        print 'p = ' + str(p) + '   k = ' + str(n)
        mesh = dolfin.RectangleMesh(dolfin.Point(xBounds[0],yBounds[0],0.0),dolfin.Point(xBounds[1],yBounds[1],0.0),n,n,'crossed')

        # define the function spaces
        U = dolfin.FunctionSpace(mesh,'RT',p)      # velocity
        PSI = dolfin.FunctionSpace(mesh,'CG',p)    # stream function
        P = dolfin.FunctionSpace(mesh,'DG',p-1)    # divergence of velocity


        # compute the finite element approximation of the analytical stream function
        psi_h = dolfin.project(psi_exact,PSI,solver_type='lu')

        # initialize the velocity and its divergence
        u_h = dolfin.Function(U)
        p_h = dolfin.Function(P)


        # compute the velocity 
        u_h.assign(dolfin.project(dolfin.curl(psi_h),U,solver_type='lu'))

        # compute the divergence 
        p_h.assign(dolfin.project(dolfin.div(u_h),P,solver_type='lu'))


        # get the maximum value of the divergence
        div_max[kp,kn] = numpy.abs(p_h.vector().array()).max()



pylab.semilogy(pp,div_max)
pylab.legend(n_elements,loc=4)
pylab.xlabel('p')
pylab.show()

Answer 1

Maybe it's too late, but your function psi does not vanish at the boundary. Then $ div rot \psi \neq \theta$. Another problem is that you are assuming that
div (project curl(u))=div(curl(project(u)))
and I am not quite sure about that.

Comments

Thank you for your comment.

I do not agree with you on both points.

1.An example that contradicts your statement is: \psi(x,y) = 1. Another example is \psi(x,y) = x. div rot \psi = 0 is an identity, regardless of psi vanishing at the boundary or not. It is valid on simply connected domains.

2.div (project curl(u))=div(curl(project(u)))

I agree with what you state if u is a general function but not in the case I present.

In the case I present u \in CG_p. Therefore (and this is the reason why RT elements have been constructed) curl u \in RT_{p-1} (may be p, depending on how RT is defined). This means:

a) project(u) = u
b) project(curl u) = curl u (or at least is should be)

Another property is that if q \in RT_{p-1} then div q \in DG_{p-1}.