Order of convergence largely oscillates

1

I am using GMSH to produce a sequence of meshes of an annulus (centered at 0, inner radius 1, outer radius 2), with greater resolution on the inner circle to get better accuracy:
Schermata da 2022-07-09 18-11-34

I am solving a poisson equation with exact solution u_ex = 0 if x[0]<=0, x[0]**2 otherwise, Dirichlet BCs on the inner ring, Neumann on the outer ring.

Computing the achieved order of convergence (L2 norm) yields: [ 3.04670631 -0.08825267 2.56832882 2.30941509 1.64144749 0.59804678 2.70940579], which oscillates a lot (the solution is H2 regular). A similar behaviour is obtained with Dirichlet BCs everywhere.

What could be the explanation for this? Is it a problem of the imposition of the boundary conditions?

Thank you in advance…

Let me post the relevant code, mesh generation comes below:



# %% Some data (e.g. exact solution and what not)

# Mesh path
mesh_path = "/home/annulus/"

# Resolutions
resolutions = [0.2, .1, .05, .025, .0125, 0.00625, 0.003125] # mesh width array


# Exact solution and corresponding source
class solution(UserExpression):
    def __init__(self, **kwargs):
        """ Construct the source function """
        super().__init__(self, **kwargs)

    def eval(self, values, x):
        """ Evaluate the source function """
        if x[0] >= 0:
            values[0] = -x[0] ** 2 / 2
        else:
            values[0] = 0


class source(UserExpression):
    def __init__(self, **kwargs):
        """ Construct the source function """
        super().__init__(self, **kwargs)

    def eval(self, values, x):
        """ Evaluate the source function """
        if x[0] >= 0:
            values[0] = 1
        else:
            values[0] = 0


class normal_derivative(UserExpression):
    def __init__(self, **kwargs):
        super().__init__(self, **kwargs)

    def eval(self, values, x):
        if x[0] >= 0:
            if x[0] ** 2 + x[1] ** 2 < 1.5 ** 2:
                values[0] = + x[0]**2
            else:
                values[0] = - x[0]**2 / 2
        else:
            values[0] = 0


u_ex = solution()
f = source()
f_D = solution()  # Dirichlet BC
f_N = normal_derivative() # Neumann BC


# %% Read that mesh

def read_mesh(resolution, mesh_path=mesh_path):
    # The volumetric mesh
    mesh = Mesh()
    with XDMFFile(mesh_path + "mesh_" + str(resolution) + ".xdmf") as infile:
        infile.read(mesh)

    # The boundary conditions
    mvc = MeshValueCollection("size_t", mesh, 1) 
    with XDMFFile(mesh_path + "facet_mesh_" + str(resolution) + ".xdmf") as infile:
        infile.read(mvc, "name_to_read")
    mf = MeshFunctionSizet(mesh, mvc)  # remember, tag 3 is inner ring, tag 2 outer ring

    L1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
    V = FunctionSpace(mesh, L1)

    return mesh, mf, V


# %% The problem

def solve_single_pb(resolution):

    (mesh, mf, V) = read_mesh(resolution)

    plot(mesh)
    plt.show()

    U = TrialFunction(V)
    v = TestFunction(V)

    # Functional formulation
    dOuterRing = Measure("ds", subdomain_data=mf, subdomain_id=2)

    F = inner(grad(U), grad(v)) * dX - f_N * v * dOuterRing - f * v * dX
    a, L = lhs(F), rhs(F)

    # Boundary conditions
    bcs = [DirichletBC(V, f_D, mf, 3)]  # inner ring, outer ring has already Neumann BCs

    u = Function(V)
    solve(a == L, u, bcs)

    err = errornorm(u_ex, u) 

    return err

# %% Error checking

err_vec = []

for res in resolutions:
    err_vec.append(solve_single_pb(res))

err_vec = np.array(err_vec)
ooc = -np.log2( err_vec[1:]/err_vec[:-1])

print(ooc)

And here, the mesh generation code:


# %%
resolution = .0125/8

# An empty geometry
geometry = pygmsh.geo.Geometry()
# Create a model to add data to
model = geometry.__enter__()

# A circle centered at the origin and radius 1
circle = model.add_circle([0.0, 0.0, 0.0], 2.0, mesh_size=5*resolution)  # meshes are always 3D, I will suppress the third component in case

# A hole
hole = model.add_circle([0.0, 0.0, 0.0], 1.0, mesh_size=1*resolution)

# My surface
plane_surface = model.add_plane_surface(circle.curve_loop,[hole.curve_loop])

# Sinchronize, before adding physical entities
model.synchronize()

#%%

# Tagging/marking boundaries and volume. 
# Boundaries with the same tag should be added simultaneously
model.add_physical([plane_surface], "volume")
model.add_physical(circle.curve_loop.curves, "outer_ring")
model.add_physical(hole.curve_loop.curves, "inner_ring")

# Generate the mesh
geometry.generate_mesh(dim=2)
geometry.save_geometry(path+"annulus_"+str(resolution)+".geo_unrolled")  
gmsh.write(path+"mesh_"+str(resolution)+".msh")
gmsh.clear()
geometry.__exit__()

#%%##########################
# Saving to a better format #
#############################

mesh_from_file = meshio.read(path+"mesh_"+str(resolution)+".msh")

def create_mesh(mesh, cell_type, prune_z = False):
    cells = mesh.get_cells_type(cell_type)  # get the cells of some type: it will change!
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    points = mesh.points[:,:2] if prune_z else mesh.points
    out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    return out_mesh

# Using the above function, create line and "plane" mesh
line_mesh = create_mesh(mesh_from_file, "line", prune_z=True)
meshio.write(path+"facet_mesh_"+str(resolution)+".xdmf", line_mesh)
triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True)
meshio.write(path+"mesh_"+str(resolution)+".xdmf", triangle_mesh)
2

Please format your code with markdown syntax, i.e.

```python
def f(x):
    return x
```

Also, without having the complete code, it is very hard to give you any guidance to what is wrong.
There are several things you should consider:

  • Is the source term and/or boundary conditions satisfying the exact equation?
  • what is the magnitude of the actual error?
  • what happens if you use the built in refine function as opposed to reading in new meshes per case?
3

Thank you for the answer: this is the complete code, mesh generation aside, which I have added in my post.

Let me answer to the bullet points:

  • the source and BCs are built from the exact solution
  • here are some values of the L2 error 8.46565131e-02, 1.02449629e-02, 1.08912340e-02, 1.83625535e-03, 3.70450223e-04, 1.18742156e-04, 7.84467417e-05
  • I am not doing this for a new reasons: the rings are not refined in the way I would like by the built-in refine function (i.e. they don’t get smoother), plus, I don’t know how to “refine the facet indicator”
4

Have you considered this?

5

No, since I will be using order 1 finite elements, I should observe order 2 convergence even with a linear approximation of geometry.

After some more tests:

  • the same setup, without hole, has order 2 convergence without oscillations
  • interestingly enough, with respect to my original code, the situation changes to [2.72917419 1.56388078 2.46337273 2.09434059 2.07271966 1.49519038 2.19333781], which looks more promising, if I perform the following modification to the dirichlet BCs:
    def boundary(x):
        return x[0]**2+x[1]**2<1.5**2
    # bcs = [DirichletBC(V, f_D, self.mf, 3)]  # inner ring, outer ring has already Neumann BCs
    bcs = [DirichletBC(V, f_D, boundary)]  # inner ring, outer ring has already Neumann BCs

Do you have any idea why this change is making the situation better? Unfortunately I cannot keep it as is, because I will be doing shape optimization on the inner ring, so, checking the boundary in this way, is no longer doable.

6

As far as I can tell, you have an issue that you are trying to enforce a discontinuous derivative at the line x=0, which is not aligning with your discrete domain.

If you make sure that the line x=0 is align with the facets of the mesh, I would expect that the convergence of the problem will improve.

I would also suggest to use gmsh to mark the physical groups, as it would avoid adding discretization accuracy to affect the selection of degrees of freedom for the boundary conditions.

2 Likes
7

Many thanks for your help.

Indeed, removing the jump from the analytical solution, the OOCs are constantly (above) 2.

By the way, I would be curious of your reasons in arguing that a discontinuous second derivative, not aligning with the mesh, would produce this behaviour.

Useless part
However, I wonder why I get two different behaviours when enforcing dirichlet BCs in the two ways:

``` python
def boundary(x):

  • return x[0]*2+x[1]2<1.52
    # bcs = [DirichletBC(V, f_D, self.mf, 3)] # inner ring, outer ring has already Neumann BCs
    bcs = [DirichletBC(V, f_D, boundary)] # inner ring, outer ring has already Neumann BCs
    ```
    *and *

python* *# def boundary(x):* *# return x[0]**2+x[1]**2<1.5**2* *bcs = [DirichletBC(V, f_D, self.mf, 3)] # inner ring, outer ring has already Neumann BCs* *# bcs = [DirichletBC(V, f_D, boundary)] # inner ring, outer ring has already Neumann BCs* *

The two behaviours are:

** first way: 2.72917419 1.56388078 2.46337273 2.09434059 2.07271966 1.49519038 2.19333781*
** original way (with meshfunction): 3.04670631 -0.08825267 2.56832882 2.30941509 1.64144749 0.59804678 2.70940579*

*Is it possible that the mesh function is wrong? *
I have included the code in the original post, for the mesh generation, but I have tripled checked and everything seems okay.

I am also attaching a plot of the following command, for the two ways of imposing Dirichlet boundary conditions:

python* *plot(mesh)* *u = Function(V)* *bcs[0].apply(u.vector())* *pl = plot(u)* *plt.colorbar(pl)* *plt.show()* *

Results:

Schermata da 2022-07-11 15-30-01
Schermata da 2022-07-11 15-30-12

Edit: my mistake here, DirichletBC will modify even interior nodes if instructed to.