# Solution breaks down for finer mesh/higher degree of FE

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 \

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 Newmark scheme:

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?

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.