Low accurate rate in dynamic mechanics problem

\rho u_{tt}=\frac{\partial \sigma}{\partial x}
There is a constant applied to the left end of a 1D bar. Right end is fixed, the system is initial with zero displacements.
After I tried some different time step size, the accurate rate about L1 error is only 1, the rate is calculated by
r=\frac{log(E_1/E_2)}{log(h_1/h2)}

the code is

from dolfin import *
import numpy as np

E=3300000000.0
rho=1300.0
c = sqrt(E/rho)

Po=30000000.0
t1=0.00015
t2=0.001
b0=1.0
Length=10.0

# Time variables
dt = 0.00005; t = 0.0; T = 0.02

mesh = IntervalMesh(400,0,Length)
V=FunctionSpace(mesh, "Lagrange", 1)

class LeftBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near (x[0], 0.0)

class RightBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near (x[0], Length)

RB=RightBoundary()
LB=LeftBoundary()

boundaries = MeshFunction("size_t", mesh, 0)
boundaries.set_all(0)
RB.mark(boundaries, 1)
LB.mark(boundaries, 2)

F=30000000

ds = Measure("ds", subdomain_data=boundaries)

# Previous and current solution
u1= interpolate(Constant(0.0), V)
u0= interpolate(Constant(0.0), V)

# Variational problem at each time
u = TrialFunction(V)
v = TestFunction(V)

a = u*v*dx + dt*dt*c*c*inner(grad(u), grad(v))*dx
L = 2*u1*v*dx - u0*v*dx + dt*dt*c*c*F*v*ds(2)

bc = DirichletBC(V, 0.0, boundaries, 1)

u=Function(V)
velocity = Function(V)

vtkfile1=File('wave/u.pvd')
vtkfile2=File('wave/v.pvd')

#while t <= T:
while t < T + DOLFIN_EPS:
    A, b = assemble_system(a, L, bc)
    solve(A, u.vector(), b)
    velocity.vector()[:] = (u.vector() - u1.vector()) / dt
    u0.assign(u1)
    u1.assign(u)
    t += dt
    vtkfile1 << (u,t)
    vtkfile2 << (velocity, t)
1 Like

The rate depends on both your temporal and spatial discretization scheme, and you would have to supply an minimal working code for anyone to be of any help.

1 Like

Please take a look.

from dolfin import *
import numpy as np

E=3300000000.0
rho=1300.0
c = sqrt(E/rho)

Po=30000000.0
t1=0.00015
t2=0.001
b0=1.0
Length=10.0

# Time variables
dt = 0.00005; t = 0.0; T = 0.02

mesh = IntervalMesh(400,0,Length)
V=FunctionSpace(mesh, "Lagrange", 1)

class LeftBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near (x[0], 0.0)

class RightBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near (x[0], Length)

RB=RightBoundary()
LB=LeftBoundary()

boundaries = MeshFunction("size_t", mesh, 0)
boundaries.set_all(0)
RB.mark(boundaries, 1)
LB.mark(boundaries, 2)

F=30000000

ds = Measure("ds", subdomain_data=boundaries)

# Previous and current solution
u1= interpolate(Constant(0.0), V)
u0= interpolate(Constant(0.0), V)

# Variational problem at each time
u = TrialFunction(V)
v = TestFunction(V)

a = u*v*dx + dt*dt*c*c*inner(grad(u), grad(v))*dx
L = 2*u1*v*dx - u0*v*dx + dt*dt*c*c*F*v*ds(2)

bc = DirichletBC(V, 0.0, boundaries, 1)

u=Function(V)
velocity = Function(V)

vtkfile1=File('wave/u.pvd')
vtkfile2=File('wave/v.pvd')

#while t <= T:
while t < T + DOLFIN_EPS:
    A, b = assemble_system(a, L, bc)
    solve(A, u.vector(), b)
    velocity.vector()[:] = (u.vector() - u1.vector()) / dt
    u0.assign(u1)
    u1.assign(u)
    t += dt
    vtkfile1 << (u,t)
    vtkfile2 << (velocity, t)

Please attach your convergence rate estimation code as well, and write the code such that it produces the rates you mention in the first post

The variational form you’ve written looks like it’s a first order \theta-scheme. So achieving first order convergence seems pretty good. Consider the elastodyanmics demo for an example of a (at best 2nd order) \theta scheme for hyperbolic problems.

1 Like

Hi nate, I tried general-alpha method in another platform, which has a high order discretization in space, however, if I perform some refinement study including time, the order of accuracy of velocity is below 1, so I want to know is there a high order time discretization for wave equation with Naumann and Dirichlet boundary conditions?