Convergence problem of linear wave equation

[renol-kth]

equation

I am trying to solve the wave equation \partial_{tt} u - \Delta u = 0 in fenicsx in 2D (but keeping everything constant in the y-axis). I reformulate it into a first order system by introducing a second variable:

\partial_t u+\partial_x p = 0,\quad\partial_t p + \partial_x u = 0

I use Crank-Nicolson to discretize the \partial_t and the weak form is:

\begin{aligned} \begin{pmatrix} M & E \\ F & N \end{pmatrix} \begin{pmatrix} \vec u ^{n+1}\\ \vec p^{n+1} \end{pmatrix} = \begin{pmatrix} M \vec u^n + E \vec p^n \\ F\vec v^n +N\vec p^n \end{pmatrix} \end{aligned}

where v,q are the basis functions corresponding to u,p and I have neglected boundary integrals and constants with \Delta t for readability. The matrices are defined here:

\begin{aligned} M_{ij} &= \left \langle v_j, v_i \right \rangle_{L^2} \\ N_{ij} &= \left \langle q_j,q_i \right \rangle_{L^2} \\ E_{ij} &= \left \langle q_j b, \nabla v_i \right \rangle_{L^2} \\ F_{ij} &= \left \langle v_j b,\nabla q_i\right \rangle_{L^2} \\ \end{aligned}

the choice b:=(1,0) makes E,F represent the partial derivative \partial_x. I use the following analytic solution

u = \sin(\pi x) \cos(\pi t),\quad p = \cos(\pi x) \sin(\pi t),

which determines the initial condition for u,p and the boundary condition for u (p does not need a boudnary condition, it is uniquely determined by u)

numerical tests

I would now expect an L^2 convergence rate of 2 from both first order FEM polynomials and from Crank-Nicolson, which is also what I get.

When I choose Lagrange Finite Elements of degree 2, and I pick a fixed, small \Delta t (i tried \approx 10^{-3}) and then refine only the mesh-size, I expect an L^2 convergence rate of 3. This does not work.

Any suggestions what could be the issue?

---

Here is my code:

In the section with # parameters, one can switch the polynomial degree pdeg and how small the time-step should be.

I just realized, that this code uses V for u and W for v.

from mpi4py import MPI
import numpy as np
import ufl
from basix.ufl import element
from dolfinx.fem import (Constant, Function, dirichletbc,
                         form, functionspace, locate_dofs_topological)
from dolfinx.fem.petsc import assemble_matrix_block,assemble_vector_block
from petsc4py import PETSc
from dolfinx.mesh import *
from ufl import dx, ds, grad, inner, FacetNormal, dot
from dolfinx.fem import Function, FunctionSpace, assemble_scalar

# error
def error_L2(uh, u_ex, degree_raise=3):
    
    degree = uh.function_space.ufl_element().degree()
    family = uh.function_space.ufl_element().family()
    mesh = uh.function_space.mesh
    W = FunctionSpace(mesh, (family, degree + degree_raise)) # high-dim space

    # interpolate
    u_W = Function(W)
    u_W.interpolate(uh)
    u_ex_W = Function(W)
    u_ex_W.interpolate(u_ex)

    # Compute error
    e_W = Function(W)
    e_W.x.array[:] = u_W.x.array - u_ex_W.x.array
    error = form(ufl.inner(e_W, e_W) * ufl.dx)
    error_local = assemble_scalar(error)
    error_global = mesh.comm.allreduce(error_local, op=MPI.SUM)
    return np.sqrt(error_global)

# exact solution
class V_exact:
    def __init__(self):
        self.t = 0.0
        
    def eval(self, x):
        return np.sin(np.pi*x[0])*np.cos(np.pi*self.t)


# Nx = 40
errors = np.array([])
hlist = np.array([])

scalings = [1,2,4,8,16,32,64,128,256,512,1042,2048]
for i in range(len(scalings)-3):

    # parameters
    Nx = 1*scalings[i]*1
    Ny = 1*scalings[i]*1
    Nt = 1*scalings[i]*1
    # Nt = 1024 # UNCOMMENT THIS
    T_end = 2*np.pi+0.5
    h = 1/Nx
    dt = T_end/Nt
    hlist = np.append(hlist,h)
    t = 0
    pdeg = 1  # CHANGE TO 2

    # mesh, spaces, functions
    msh = create_unit_square(MPI.COMM_WORLD, Nx, Nx)
    x = ufl.SpatialCoordinate(msh)
    n = FacetNormal(msh)
    P1 = element("Lagrange", "triangle", pdeg) 
    XV = functionspace(msh, P1) 
    Xp = functionspace(msh, P1) 
    V = ufl.TrialFunction(XV)
    p = ufl.TrialFunction(Xp)
    W = ufl.TestFunction(XV)
    q = ufl.TestFunction(Xp)
    V_old = Function(XV)
    V_new = Function(XV)
    p_old = Function(Xp)
    p_new = Function(Xp)
    b = Constant(msh, (1.,0.))

    # init cond
    def initcond_V(x): return np.sin(np.pi * x[0])
    def initcond_p(x): return 0+0.0*x[0]
    V_old.interpolate(initcond_V)
    p_old.interpolate(initcond_p)
    V_new.interpolate(initcond_V)
    p_new.interpolate(initcond_p)

    # BC
    facetsV = locate_entities_boundary(msh, dim=1,
                                        marker=lambda x: np.logical_or.reduce((
                                            np.isclose(x[0], 0.0),
                                            np.isclose(x[0], 1.0))))
    dofsV = locate_dofs_topological(XV, entity_dim=1, entities=facetsV)
    BCs = [dirichletbc(0.0,dofsV, XV)]

    # weak form
    M = inner(V,W)*dx 
    E = inner(grad(p),b*W)*dx
    F = inner(grad(V),b*q)*dx
    N = inner(p,q)*dx
    a = form([[M, -dt/2*E],
             [-dt/2*F, N]])
    A = assemble_matrix_block(a, bcs=BCs)
    A.assemble()
    Vp_vec = A.createVecLeft()

    # solver
    solver = PETSc.KSP().create(msh.comm)
    solver.setOperators(A)
    solver.setType(PETSc.KSP.Type.PREONLY)
    solver.getPC().setType(PETSc.PC.Type.LU)

    for i in range(Nt):
        t += dt

        # Update the right hand side
        l = form([inner(V_old,W)*dx + dt/2*inner(grad(p_old),b*W)*dx - dt*inner(p_old*dot(b,n),W)*ds,
                    inner(p_old,q)*dx + dt/2*inner(grad(V_old),b*q)*dx - dt*inner(V_old*dot(b,n),q)*ds])        
        L = assemble_vector_block(l,a,bcs=BCs)
        
        # solve
        solver.solve(L, Vp_vec)

        # extract solution
        offset = XV.dofmap.index_map.size_local * XV.dofmap.index_map_bs
        V_new.x.array[:offset] = Vp_vec.array_r[:offset]
        p_new.x.array[:(len(Vp_vec.array_r) - offset)] = Vp_vec.array_r[offset:]
                
        # Update solution at previous time step
        V_old.x.array[:] = V_new.x.array
        p_old.x.array[:] = p_new.x.array

    # error in high-dim polynomial space
    V_ex_expr = V_exact()
    V_ex_expr.t = T_end
    error = error_L2(V_new, V_ex_expr.eval, degree_raise=5) 
    errors = np.append(errors,error)
    print(error)

# print errors and convergence rate
# print(errors)
rates = np.log(errors[1:] / errors[:-1]) / np.log(hlist[1:] / hlist[:-1])
print(rates)

[dokken]

The first thing i did with your code is to check that the initial condition is correct, i.e.

from mpi4py import MPI
import numpy as np
import ufl
from basix.ufl import element
from dolfinx.fem import (Constant, Function, dirichletbc,
                         form, functionspace, locate_dofs_topological)
from dolfinx.fem.petsc import assemble_matrix_block,assemble_vector_block
from petsc4py import PETSc
from dolfinx.mesh import *
from ufl import dx, ds, grad, inner, FacetNormal, dot
from dolfinx.fem import Function, functionspace, assemble_scalar

# error
def error_L2(uh, u_ex, degree_raise=3):
    
    degree = uh.function_space.ufl_element().degree

    family = uh.function_space.ufl_element().element_family.name
    mesh = uh.function_space.mesh
    W = functionspace(mesh, (family, degree + degree_raise)) # high-dim space

    # interpolate
    u_W = Function(W)
    u_W.interpolate(uh)
    u_ex_W = Function(W)
    u_ex_W.interpolate(u_ex)

    # Compute error
    e_W = Function(W)
    e_W.x.array[:] = u_W.x.array - u_ex_W.x.array
    error = form(ufl.inner(e_W, e_W) * ufl.dx)
    error_local = assemble_scalar(error)
    error_global = mesh.comm.allreduce(error_local, op=MPI.SUM)
    return np.sqrt(error_global)

# exact solution
class V_exact:
    def __init__(self):
        self.t = 0.0
        
    def eval(self, x):
        return np.sin(np.pi*x[0])*np.cos(np.pi*self.t)


# Nx = 40
errors = np.array([])
hlist = np.array([])

scalings = [1,2,4,8,16,32,64,128,256,512,1042,2048]
for i in range(len(scalings)-3):

    # parameters
    Nx = 1*scalings[i]*1
    Ny = 1*scalings[i]*1
    Nt = 1*scalings[i]*1
    # Nt = 1024 # UNCOMMENT THIS
    T_end = 2*np.pi+0.5
    h = 1/Nx
    dt = T_end/Nt
    hlist = np.append(hlist,h)
    t = 0
    pdeg = 1  # CHANGE TO 2

    # mesh, spaces, functions
    msh = create_unit_square(MPI.COMM_WORLD, Nx, Nx)
    x = ufl.SpatialCoordinate(msh)
    n = FacetNormal(msh)
    P1 = element("Lagrange", "triangle", pdeg) 
    XV = functionspace(msh, P1) 
    Xp = functionspace(msh, P1) 
    V = ufl.TrialFunction(XV)
    p = ufl.TrialFunction(Xp)
    W = ufl.TestFunction(XV)
    q = ufl.TestFunction(Xp)
    V_old = Function(XV)
    V_new = Function(XV)
    p_old = Function(Xp)
    p_new = Function(Xp)
    b = Constant(msh, (1.,0.))

    # init cond
    def initcond_V(x): return np.sin(np.pi * x[0])
    def initcond_p(x): return 0+0.0*x[0]
    V_old.interpolate(initcond_V)
    p_old.interpolate(initcond_p)
    V_new.interpolate(initcond_V)
    p_new.interpolate(initcond_p)

    # BC
    facetsV = locate_entities_boundary(msh, dim=1,
                                        marker=lambda x: np.logical_or.reduce((
                                            np.isclose(x[0], 0.0),
                                            np.isclose(x[0], 1.0))))
    dofsV = locate_dofs_topological(XV, entity_dim=1, entities=facetsV)
    BCs = [dirichletbc(0.0,dofsV, XV)]

    # weak form
    M = inner(V,W)*dx 
    E = inner(grad(p),b*W)*dx
    F = inner(grad(V),b*q)*dx
    N = inner(p,q)*dx
    a = form([[M, -dt/2*E],
             [-dt/2*F, N]])
    A = assemble_matrix_block(a, bcs=BCs)
    A.assemble()
    Vp_vec = A.createVecLeft()

    # error in high-dim polynomial space
    V_ex_expr = V_exact()
    V_ex_expr.t = t
    error = error_L2(V_new, V_ex_expr.eval, degree_raise=5) 
    errors = np.append(errors,error)
    print(error)

# print errors and convergence rate
# print(errors)
rates = np.log(errors[1:] / errors[:-1]) / np.log(hlist[1:] / hlist[:-1])
print(rates)

which gives the expected rates. So far so good.

Secondly, I would move the form definition out of the loop:

    # Update the right hand side
    l = form([inner(V_old,W)*dx + dt/2*inner(grad(p_old),b*W)*dx - dt*inner(p_old*dot(b,n),W)*ds,
                inner(p_old,q)*dx + dt/2*inner(grad(V_old),b*q)*dx - dt*inner(V_old*dot(b,n),q)*ds])        
    L = create_vector_block(l)
    for i in range(Nt):
        t += dt
        L.zeroEntries()
        L = assemble_vector_block(l,a,bcs=BCs)
        
        # solve
        solver.solve(L, Vp_vec)

        # extract solution
        offset = XV.dofmap.index_map.size_local * XV.dofmap.index_map_bs
        V_new.x.array[:offset] = Vp_vec.array_r[:offset]
        p_new.x.array[:(len(Vp_vec.array_r) - offset)] = Vp_vec.array_r[offset:]
                
        # Update solution at previous time step
        V_old.x.array[:] = V_new.x.array
        p_old.x.array[:] = p_new.x.array

Secondly, have you tried to implement this with a non-blocked method? It is always a good reference to ensure that all variables are updated correctly.

[renol-kth]

Thank you for the quick response and your advice - both on the form-definition and the idea with the matrix assembly.

1. I assume with “non-blocked” you mean the same assembly method that is used in the last (4th) part of the stokes tutorial, i.e. using a “monolithic matrix”, and a mixed functionspace, right?
2. As my implementation yields the correct convergence rate for the order-1 elements but fails for the order-2 ones, what is your intuition why a different matrix assembly would change/fix this?

[dokken]

@renol-kth

assume with “non-blocked” you mean the same assembly method that is used in the last (4th) part of the stokes tutorial, i.e. using a “monolithic matrix”, and a mixed functionspace, right?

Yes.

@renol-kth

As my implementation yields the correct convergence rate for the order-1 elements but fails for the order-2 ones, what is your intuition why a different matrix assembly would change/fix this?

I would make dt even smaller, and do only a few time steps, as the interpolation of the exact solution has the correct rates.

[renol-kth]

I ran the code for 2nd order polynomials and observe this:

\Delta t	N_t	order 3 conv
0.1	1	No
0.01	1	Yes
0.01	10	No
0.001	10	Yes
0.001	100	No
0.0001	100	Yes

Of course, the error should get larger, when simulating longer, but it seems like the convergence order itself is messed up, when simulating longer. I did not expect this… Do you have an idea what could be the problem?

[renol-kth]

For future reference:

1. the convergence rate with p=2 polynomials can be shown to be 3, if one refines the time step faster than the mesh size, rather than choosing Nt fixed like I did in the code above.
2. the boundary integrals are wrong, they should be removed. We have dirichlet BC and therefore they vanish (on both equations), but this was not the reason for the convergence problem.