Solution breaks down for finer mesh/higher degree of FE

mesh
1

Hello, this is a continuation of efficient-way-to-sum-up-fem-functions, where I was solving the homogenous linear wave equation.

While doing so, I have noticed that the solution breaks down for Newmark’s time discretization scheme (vs central scheme) in the case of a finer mesh and a higher degree of finite element.

Code:

from mpi4py import MPI
from dolfinx import fem, mesh, io

import dolfinx.fem.petsc
import ufl
import numpy as np

n_cells = 100
degree = 2
TIME_SCHEME = "Newmark" # central | Newmark

msh = mesh.create_unit_square(MPI.COMM_WORLD, n_cells, n_cells, cell_type=mesh.CellType.quadrilateral)

t = fem.Constant(msh, 0.)
dt = 5e-3
T = 1 
c = 1

V = fem.FunctionSpace(msh, ("Lagrange", degree))

uh = fem.Function(V)
v = ufl.TestFunction(V)
u = ufl.TrialFunction(V)

def u_init(x, sigma=0.05, mu=0.5):
    return np.exp(-0.5*((x[0]-mu)/sigma)**2)*np.exp(-0.5*((x[1]-mu)/sigma)**2) 

if TIME_SCHEME == "central":
    u_n = fem.Function(V)
    u_nn = fem.Function(V)
        
    dudtdt = (u - 2 * u_n + u_nn) / dt**2
    du = (u + u_nn) / 2

    u_n.interpolate(u_init)
    u_nn.interpolate(u_init)

elif TIME_SCHEME == "Newmark":
    u_n = fem.Function(V)
    dudt_n = fem.Function(V)
    dudtdt_n = fem.Function(V)

    u_temp = fem.Function(V)
    dudt_temp = fem.Function(V)

    beta = 0.25
    gamma = 0.5

    dudtdt = u
    du = u_n + dudt_n*dt + ((1-beta)*dudtdt_n + 2*beta*dudtdt) * dt**2/2

    u_n.interpolate(u_init)
    dudt_n.interpolate(u_init)
    dudtdt_n.interpolate(u_init)
    
F = 1 / c**2 * ufl.inner(dudtdt, v) * ufl.dx \
    + ufl.inner(ufl.grad(du), ufl.grad(v)) * ufl.dx \

a, L = ufl.system(F)
petsc_options = {"ksp_type": "preonly",
                "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"}
problem = dolfinx.fem.petsc.LinearProblem(a, L, u=uh, bcs=[], petsc_options=petsc_options)

from pathlib import Path
folder = Path("FENICS - results/mine_out_wave_equation_issue")
vtx_u = io.VTXWriter(msh.comm, folder / f"results_{TIME_SCHEME}_cells{n_cells}_degree{degree}.bp", u_n, engine="BP4") 
vtx_u.write(0)

while t.value < T:
    t.value += dt
    problem.solve()
                            
    if TIME_SCHEME == "central":
        u_nn.x.array[:] = u_n.x.array
        u_n.x.array[:] = uh.x.array  
        
    elif TIME_SCHEME == "Newmark":
        u_temp.x.array[:] = u_n.x.array + dudt_n.x.array*dt + ((1-beta)*dudtdt_n.x.array + 2*beta*uh.x.array) * dt**2/2
        dudt_temp.x.array[:] = dudt_n.x.array + ((1-gamma)*dudtdt_n.x.array + gamma*uh.x.array) * dt

        u_n.x.array[:] = u_temp.x.array
        dudt_n.x.array[:] = dudt_temp.x.array
        dudtdt_n.x.array[:] = uh.x.array

    vtx_u.write(t.value)

Results for central scheme:

image.jpg566×698 94.8 KB
image.jpg569×692 134 KB

image.jpg567×694 92.2 KB
image.jpg569×700 136 KB

Results for Newmark scheme:

image.jpg568×693 96.7 KB
image.jpg556×693 131 KB

image.jpg560×696 95 KB
image.jpg567×691 130 KB

2

Do those time stepping methods need to satisfy a CFL-like condition, i.e. if you decrease h you also need to decrease appropriately dt?

3

Newmark’s method is implicit and should be unconditionally stable, whereas the central scheme is explicit. So from the results with the central scheme on the finest mesh, I assume that the time step is okay.

EDIT: So sorry for the misconception, but lowering time step actually helped to solve the problem. However, I am confused by that.