I am trying to build understanding on how to deal with the reflection of acoustic waves using FEniCS. It is a fairly simple problem, yet I can’t manage to get it working for I am new to FEniCS, so any basic comments on my code will help me a lot!
The domain is a unit square with a point source in the middle of it. I have:
The Wave Equation \partial u^2 / \partial t^2 = c^2 \nabla^2 u + f
The Weak Form a(u,v) = \int_{\Omega} (vu^{n+1} - c^2 \Delta t^2 \nabla v \nabla u^{n+1})dx L_{n+1}(v) = \int_{\Omega}(2u^{n} - u^{n-1} + \Delta t^2 f^{n+1})vdx - \int_{\Gamma^N} (gv)ds
Neumann boudary condition on all the 4 walls \partial u / \partial n = g = 0
The code I have written with some pieces of examples
from fenics import *
import numpy as np
import matplotlib.pyplot as plt
T = 5.0
num_steps = 250
dt = T / num_steps
c = 1
# Create mesh and define function space
nx = 80
ny = 80
mesh = RectangleMesh(Point(-1,-1), Point(1,1),nx,ny)
V = FunctionSpace(mesh, "Lagrange", 1)
u = TrialFunction(V)
v = TestFunction(V)
#Initial pressure
init = "0"
g_ini = Expression(init,t=0, degree=0)
u_0 = interpolate(g_ini,V)
u_1 = interpolate(g_ini,V)
#Null source term
f = Constant(0)
#Neumann boundary condition
g = Constant('0.0')
#Variational Form
a = u*v*dx - c*c*dt*dt*dot(grad(u), grad(v))*dx
L = (2*u_1 - u_0 + dt*dt*f)*v*dx - g*v*ds
# Time-stepping
u = Function(V)
t = 0
file = File("wave/wave.pvd")
for n in range(num_steps):
#Point source
PS = PointSource(V, Point(0, 0), sin(10 * t))
PS.apply(u.vector())
#Solve, no boundaries specified as I have only Neumann
solve(a == L, u)
t += dt
plot(u)
#Update solutions
u_0.assign(u_1)
u_1.assign(u)
file<<u
The output is completely static. No point source and no reflection. What am I missing ?
You should define a VectorFunctionSpace NOT a FunctionSpace as displacement is not scalar but is a vector
In your weak form you should use inner product for multiplication of the trial and test functions (Same thing for the linear part of the variational form)
The terms in the bilinear part of the variational form, should be added NOT subtracted
I did some modifications in your code. The below should work:
from fenics import *
import numpy as np
t = 0
T = 5.0
num_steps = 250
dt = T / num_steps
c = 1
nx = 80
ny = 80
mesh = RectangleMesh(Point(-1, -1), Point(1, 1), nx, ny)
V = VectorFunctionSpace(mesh, 'CG', 1)
u = TrialFunction(V)
v = TestFunction(V)
g_ini = Constant((0.0, 0.0))
u_0 = interpolate(g_ini, V)
u_1 = interpolate(g_ini, V)
a = inner(u, v) * dx + c * c * dt * dt * inner(grad(u), grad(v)) * dx
L = 2*inner(u_1,v)*dx-inner(u_0,v)*dx
bc = DirichletBC(V,g_ini, "on_boundary")
u = Function(V)
file = File("wave/wave.pvd")
while t <= T:
A, b = assemble_system(a, L, bc)
PS = PointSource(V, Point(0, 0), sin(10 * t))
PS.apply(b)
solve(A, u.vector(), b)
u_0.assign(u_1)
u_1.assign(u)
t += dt
file << u
Thanks a lot for the explanations and the code @Leo, that did trick!
May I ask a second question?
For the next step I’d like to have a front wave propagating from the left boundary to the right (no point source). I tried redefining the expression for the left boundary and maintaining the others:
I realized something unexpected. Take a look at the output, it shows the displacement in the X coordinate. The amplitudes are getting lower as they drive away from the source.
Well, based on what I am seeing here, the problem could be the size of the time step. You may want to use smaller time steps (0.001 instead of 0.004). Please look at the results of this quick test:
Hi, sorry to bother again, but I was reading back through this and thinking what if I wanted to declare a pressure wave equation? Pressure is a scalar, so it should be OK to define a FunctionSpace and use regular multiplication in the weak form, right?
I did a quick test and I show a comparison between defining a FunctionSpace (left) and a VectorFunctionSpace (right). They seem very similar to me. Do you have any comment on that?
I’m asking this because my next step is to define absorbing boundary conditions, maybe through perfectly-matched-layers, and the only references I found use pressure to deduce the PMLs.
