Error norm no change for the Poisson equations with two subdomains

1

Hi,

I’m recently working on studying the subdomain in Fenics version 2017.02. I want to test the convergence rate for the Poisson interface problem list below:

image
image1629×1000 132 KB

I want to use the traditional FEM with interface fitted mesh which I suppose to have 2nd order convergence rate with P1 element. When I use the exact solution on domain Ω:[-1,-1]x[1,1], and the interface is a circle centered at the origin with radius r0=1/2, \alpha=2, ß has different values across the interface, where ß=1000 inside of the circle and 1 otherwise.

I will have f=-4 and g=0 with the corresponding exact solution for the problem. So I generated the mesh with marked subdomains and I got the numerical solution that looks similar to the exact solution but the error norm didn’t change too much with the mesh refine N=40, N=80, N=160. I can’t figure it out for a long time. Any help would be really appreciated!

Please find below for the code:

from fenics import *
from mshr import *
# define variables
ax = -1.0
ay = -1.0
bx = 1.0
by = 1.0
Nx=Ny=40
r0=1/2

# Define geometry for background
domain = Rectangle(Point(ax,ay),Point(bx,by))

# Define geometry for subdomains
out = Rectangle(Point(ax,ay),Point(bx,by))-Circle(Point(0, 0), r0)
circle = Circle(Point(0, 0), r0)
# Set subdomains
domain.set_subdomain(1, out)
domain.set_subdomain(2, circle)

mesh = generate_mesh(domain, Nx)
meshfile = File("mesh.pvd") # save mesh data 
meshfile << mesh

hmax = mesh.hmax()
print('max h', hmax, 'for mesh N', Nx)

V = FunctionSpace(mesh, 'P', 1)

alfa = 2
beta1 = 1
beta2 = 1000
class Beta(Expression):
    def __init__(self,markers, **kwargs):
        self.markers = markers
    def eval_cell(self, values, x, cell):
        if self.markers[cell.index] == 1:  
            values[0] = beta1
        else:
            values[0] = beta2
beta = Beta(markers, degree = 1)

# Define boundary condition
u_D = Expression('pow(sqrt(pow(x[0],2) + pow(x[1],2)),alfa)/beta1+ (1/beta1 - 1/beta2)*pow(r0,alfa)',alfa=alfa, beta1=beta1, beta2=beta2,r0=r0, degree=2)

def boundary(x):
    tol = 1e-15
    return abs(x[0]-ax)< tol or abs(x[0]-bx)< tol or abs(x[1]-ay)< tol or abs(x[1]-by)< tol

bc = DirichletBC(V, u_D, boundary)

dx = Measure('dx', domain=mesh, subdomain_data=markers)
# ds = Measure('ds', domain=mesh, subdomain_data=interface_markers)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(-4.0)
g = Constant(0.0)
a = beta*dot(grad(u), grad(v))*dx
L = f*v*dx(1) + f*v*dx(2) #+ g*v*ds

# Compute solution
u = Function(V)
solve(a == L, u, bc)

# Save solution to file in VTK format
vtkfile = File('u.pvd')
vtkfile << u

class u_exact(Expression):
    def __init__(self, subdom, **kwargs):
        self.subdom = subdom
    def eval_cell(self, values, x, cell):
        if self.subdom[cell.index] == 1:
            values[0] = pow(sqrt(pow(x[0],2) + pow(x[1],2)),alfa)/beta1 + (1/beta1 - 1/beta2)*pow(r0,alfa)
        else:
            values[0] = pow(sqrt(pow(x[0],2) + pow(x[1],2)),alfa)/beta2

u_e = u_exact(markers, degree = 2)

u_e_ = interpolate(u_e, V)
file = File("u_e.pvd")
file << u_e_
# Compute error in L2 norm
error_L2 = errornorm(u_e, u, 'L2')

# Compute maximum error at vertices
vertex_values_u_D = u_e.compute_vertex_values(mesh)
vertex_values_u = u.compute_vertex_values(mesh)
import numpy as np
error_max = np.max(np.abs(vertex_values_u_D - vertex_values_u))

# Print errors
# print(u.vector().get_local())
# print(u_e_.vector().get_local())

print('error_L2', error_L2, 'for N =', Nx)
print('error_max', error_max, 'for N =', Nx)
2

See for instance:

for an example with multi material problems with Poisson. Try to first implement a similar approach in your code, and then extend it to 2D when you Get it working.

3

Hi dokken, thanks a lot for your reply!

I read the instance you attached. But I think I already considered the jump conditions in my Variation Form by <g,v> term. Under the exact solution that I want to approach, the g value should be zero so I just omit that term in my code. Did I think something wrong?

Thanks

4

Note that the code you are supplying does not run, as markers are not defined.

I am also wondering why you are using the outdated 2017.2.0 version of FEniCS?

I linked you to that method such that you can use it as a reference on how to implement a multi material poisson problem. As you see in that example, by excluding the internal jump term, you get a wanted discontinuity in the gradient.

5

Hi Dokken,

Sorry that I accidentally deleted one line of the code when I copied here. I should have markers = MeshFunction('size_t', mesh, 2, mesh.domains()) after I defined the FunctionSpace.

Please see the new attached code for your convenience

from fenics import *
from mshr import *
# define variables
ax = -1.0
ay = -1.0
bx = 1.0
by = 1.0
Nx=Ny=40
r0=1/2

# Define geometry for background
domain = Rectangle(Point(ax,ay),Point(bx,by))

# Define geometry for subdomains
out = Rectangle(Point(ax,ay),Point(bx,by))-Circle(Point(0, 0), r0)
circle = Circle(Point(0, 0), r0)
# Set subdomains
domain.set_subdomain(1, out)
domain.set_subdomain(2, circle)

mesh = generate_mesh(domain, Nx)
meshfile = File("mesh.pvd") # save mesh data 
meshfile << mesh

hmax = mesh.hmax()
print('max h', hmax, 'for mesh N', Nx)

V = FunctionSpace(mesh, 'P', 1)

markers = MeshFunction('size_t', mesh, 2, mesh.domains())


alfa = 2
beta1 = 1
beta2 = 1000
class Beta(Expression):
    def __init__(self,markers, **kwargs):
        self.markers = markers
    def eval_cell(self, values, x, cell):
        if self.markers[cell.index] == 1:  
            values[0] = beta1
        else:
            values[0] = beta2
beta = Beta(markers, degree = 1)

# Define boundary condition
u_D = Expression('pow(sqrt(pow(x[0],2) + pow(x[1],2)),alfa)/beta1+ (1/beta1 - 1/beta2)*pow(r0,alfa)',alfa=alfa, beta1=beta1, beta2=beta2,r0=r0, degree=2)

def boundary(x):
    tol = 1e-15
    return abs(x[0]-ax)< tol or abs(x[0]-bx)< tol or abs(x[1]-ay)< tol or abs(x[1]-by)< tol

bc = DirichletBC(V, u_D, boundary)

dx = Measure('dx', domain=mesh, subdomain_data=markers)
# ds = Measure('ds', domain=mesh, subdomain_data=interface_markers)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(-4.0)
g = Constant(0.0)
a = beta*dot(grad(u), grad(v))*dx
L = f*v*dx(1) + f*v*dx(2) #+ g*v*ds

# Compute solution
u = Function(V)
solve(a == L, u, bc)

# Save solution to file in VTK format
vtkfile = File('u.pvd')
vtkfile << u

class u_exact(Expression):
    def __init__(self, subdom, **kwargs):
        self.subdom = subdom
    def eval_cell(self, values, x, cell):
        if self.subdom[cell.index] == 1:
            values[0] = pow(sqrt(pow(x[0],2) + pow(x[1],2)),alfa)/beta1 + (1/beta1 - 1/beta2)*pow(r0,alfa)
        else:
            values[0] = pow(sqrt(pow(x[0],2) + pow(x[1],2)),alfa)/beta2

u_e = u_exact(markers, degree = 2)

u_e_ = interpolate(u_e, V)
file = File("u_e.pvd")
file << u_e_
# Compute error in L2 norm
error_L2 = errornorm(u_e, u, 'L2')

# Compute maximum error at vertices
vertex_values_u_D = u_e.compute_vertex_values(mesh)
vertex_values_u = u.compute_vertex_values(mesh)
import numpy as np
error_max = np.max(np.abs(vertex_values_u_D - vertex_values_u))

# Print errors
# print(u.vector().get_local())
# print(u_e_.vector().get_local())

print('error_L2', error_L2, 'for N =', Nx)
print('error_max', error_max, 'for N =', Nx)

Sorry that the version I used is outdated since I learned FEniCS with it. It’s kind of confusing to use a new version without understanding the version I used thoroughly.

Yeah, I think I understand the example you provided. I enforced my numerical solution to be continuous on the interface under the exact solution(which is continuous on the interface) I have. I have the right variational form, right? So I’m really struggled to find out why the error norm is wrong. I’m wondering whether my subdomains work correctly.

6

This is because you have consistently written the wrong expression in your code:

It should be 1/beta2-1/beta1.
You have done this consistently throughout your code.

As a side note, I would not use eval_cell in an Expression when the function is not discontinuous.

class u_exact(Expression):
    def __init__(self, subdom, **kwargs):
        self.subdom = subdom
    def eval(self, values, x):
        if sqrt(pow(x[0],2) + pow(x[1],2)) > r0:
            values[0] = pow(sqrt(pow(x[0],2) + pow(x[1],2)),alfa)/beta1 + (1/beta2 - 1/beta1)*pow(r0,alfa)
        else:
            values[0] = pow(sqrt(pow(x[0],2) + pow(x[1],2)),alfa)/beta2

If you look at the exact solution produced by your original interpolation (the “u_e.pvd” file) you see that is does not look quite right.

7

Oh, my gosh…Thank you so much to point out that!!! I suffered a lot for this problem.

To use eval instead of eval_cell, is that because of saving calculation time? I see the meaning of eval_cell from the demo https://fenicsproject.org/pub/tutorial/html/._ftut1013.html said ‘’‘This version of the evaluation function has an additional cell argument which we can use to check on which cell we are currently evaluating the function.’’’