Initial condition for wave equations

[Miraboreasu]

If I have two initial conditions, one of them is related to the space. They are
u(x<1, t=0)=x,u(x>=1, t=0)=-x,
\frac{\partial u(x, t=0)}{\partial t}=0
I cannot find how to handle the initial conditions like these from the book of Solving PDEs in Python, anyone can please give me some material to learn?
Thank you.

[kamensky]

The most concise way to do this would be to interpolate an Expression, where you can use the C++ ternary operator to implement piecewise functions:

from dolfin import *
N = 100
mesh = IntervalMesh(N,0,2)
deg = 1
V = FunctionSpace(mesh,"CG",deg)
ic = interpolate(Expression("x[0] < 1? x[0] : -x[0]",degree=deg),V)

from matplotlib import pyplot as plt
plot(ic)
plt.show()

(The string passed as an argument to the Expression constructor is a snippet of C++ code.) Alternatively, if you have a more complicated initial condition that is difficult to write in a single line of C++, you could express it in UFL, using a SpatialCoordinate, the project it onto the desired function space. However, projection of discontinuous functions with the default L^2 projection implemented by the project function leads to oscillations (and requires solving a linear system), so you’d likely want to use a lumped-mass projection, e.g.,

from dolfin import *
N = 100
mesh = IntervalMesh(N,0,2)
deg = 1
V = FunctionSpace(mesh,"CG",deg)
x = SpatialCoordinate(mesh)
ic_ufl = conditional(lt(x[0],1),x[0],-x[0])

# Not good for projecting discontinuous data; leads to oscillation:
#ic = project(ic_ufl,V)

def lumpedProject(f):
    v = TestFunction(V)
    lhs = assemble(inner(Constant(1.0),v)*dx)
    rhs = assemble(inner(f,v)*dx)
    u = Function(V)
    as_backend_type(u.vector())\
        .vec().pointwiseDivide(as_backend_type(rhs).vec(),
                               as_backend_type(lhs).vec())
    return u

ic = lumpedProject(ic_ufl)

from matplotlib import pyplot as plt
plot(ic)
plt.show()

You can see some more discussion of implementing mass lumping in @bleyerj’s tutorial here, including an extension to P2 elements.

[Miraboreasu]

Thank you for your nice introduction, before I read the lumped projection, I don’t understand very well. My PDE is
\frac{\partial^2u}{\partial t^2}=c^2\frac{\partial^2u}{\partial x^2}
I apply a central discretization to time,
(\frac{\partial^2u}{\partial t^2})=\frac{u-{2u1}+u0}{\Delta t^2}

so I have u for the current solution, and u1, u0 for the previous, and the one before the previous. How to connect u1, u0 to the

ic = interpolate(Expression("x[0] < 1? x[0] : -x[0]",degree=deg),V)

Left and right boundary conditions are zero. My current code is

from dolfin import *
import numpy as np

Length=50
c = 1


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

mesh = IntervalMesh(100,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)



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

bc1= DirichletBC(V, 0.0, boundaries, 1)
bc2= DirichletBC(V, 0.0, boundaries, 2)
bc=[bc1, bc2]

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)
    #print(t)
    t += dt
    vtkfile1 << (u,t)
    vtkfile2 << (velocity, t)

[kamensky]

You can use

u1.assign(ic)

if ic is a Function in V (as the output of interpolate would be). The argument may also be a constant-coefficient linear combination of Functions in V. Setting u0 depends on having another initial condition for the velocity, or knowing u(t_0-\Delta t) directly.

You might also consider re-arranging the time integrator to use a displacement and velocity from the previous time step, as implemented in my answer here:

Second order time discretization problem

You can convert to a first order system but then formally eliminate the extra velocity unknown to avoid solving in a mixed function space. The implicit midpoint rule gives you \frac{u^{n+1} - u^n}{\Delta t} = \frac{1}{2}(\dot{u}^{n+1} + \dot{u}^{n})~~~\Rightarrow~~~\dot{u}^{n+1} = 2(\frac{u^{n+1} - u^n}{\Delta t}) - \dot{u}^n which can be substituted into the left-hand side of \frac{\dot{u}^{n+1} - \dot{u}^n}{\Delta t} = F((u^{n+1} + u^{n})/2) to get a problem for u^{n+1} alone, given bot…

Another note: If you initialize the wave equation with discontinuous data, you will likely get spurious oscillations in the solution without some form of additional stabilization or shock capturing. (That was one of the demonstrative examples considered in this paper. FEniCS code for the wave equation test is here, but it uses an additional library I wrote to do isogeometric analysis with FEniCS.)