### (Deleted) Wave equation 2D - solution with wrong phase?

93
views
0
9 months ago by
Hello,

I reformulate my question to which I got no reply.

I want to solve the wave equation   $u_{tt}\left(x,t\right)-u_{xx}\left(x,t\right)=f\left(x,t\right)$utt(x,t)uxx(x,t)=ƒ (x,t)  with zero initial conditions and zero boundary conditions.

Given  $f\left(x,t\right)=sin\left(2\pi t\right)h\left(x\right)$ƒ (x,t)=sin(2πt)h(x) , where h(x) is a Gaussian centered in the middle of the domain, I expect the solution to be 'essentially' equal to  $sin\left(2\pi t\right)$sin(2πt) in every point of the domain. This is what we can obtain by solving the wave equation in the frequencies domain.

Now, the code here below doesn't give such solution, but rather a sinus with a longer period. Does anyone have any experience with waves in the FEniCS community? Can anyone post a relevant solution? Any numerical time integration method is welcome.

I have attached first the code where I use simple central difference for the time integration. The following one is using the Newmark method. I'd attach plots of the result but it's not possible on this platform (I guess). However, both solutions have the same period, but the one computed with the central difference is dissipative and the amplitude of the wave decreases with time.

from dolfin import *
from mshr import *
import numpy  as np
import scipy.io as sio

# Velocity
c_ = 1.

# Create mesh
xlim, ylim = 1., 1.
size = 60
domain = Rectangle(Point(0., 0.), Point(xlim, ylim))
mesh = generate_mesh(domain,size)
hmin = mesh.hmin()
hmax = mesh.hmax()

# Order of the Finite Elements
pdim = 1

# Time Constants
t = 0.
CFL = 0.5
dt = CFL*hmax/(sqrt(2.)*c_)
print dt
T = 10.
nT = int(T/dt)+1
tvec = np.linspace(0,T,nT)

# Boundaries
class Bound(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and ((abs(x[1] - 0.) < DOLFIN_EPS) or (abs(x[1] - ylim) < DOLFIN_EPS) or (abs(x[0] - 0.) < DOLFIN_EPS) or (abs(x[0] - xlim) < DOLFIN_EPS))

boundaries = FacetFunction("size_t", mesh, 0)

Bound().mark(boundaries, 1)

# Define Function space, trial and test functions
V = FunctionSpace(mesh, 'P', pdim)
u, phi = TrialFunction(V), TestFunction(V)

# Define Dirichlet BC
bcs = DirichletBC(V, Constant(0.), boundaries, 1)

# Fields from previous time step (displacement, velocity)
v0IC = Constant(0.)
u0 = interpolate(v0IC, V)
u1 = interpolate(v0IC, V)

# Define sensor locations - where I save the solution:
xs = xlim/4
ys = ylim/2

xs2 = xlim*0.85
ys2 = ylim*0.67

XS = [xs, xs2]
YS = [ys, ys2]

sensor_x = np.zeros((1,nT))
sensor_x2 = np.zeros((1,nT))

# Define the source location and expression
gauss = '(1.0/(2.*ppi*sigmax*sigmay) * exp(- pow((x[0]-mux),2)/(2.*pow(sigmax,2)) - pow((x[1]-muy),2)/(2.*pow(sigmay,2))))'

sigma_space = min(xlim,ylim)/200 # variance of the space gaussian
xSource = xlim/2
ySource = ylim/2

toneburst = 'sin(2*ppi*fc1*t)'
t = 0.
fc1 = 1.
source_term = Expression(toneburst+'*'+gauss, fc1=fc1, ppi=np.pi, mux=xSource, muy=ySource, sigmax=sigma_space, sigmay=sigma_space, t=t, degree=pdim)

# Define the system to solve
a_lhs = u * phi * dx + pow(dt*c_,2) * inner(grad(u),grad(phi)) * dx
A = assemble(a_lhs)

a_rhs = source_term * pow(dt,2) * phi * dx  + (2. * u1 - u0) * phi * dx

u = Function(V)

tstep = 0
# Time loop
while t < T:

# Update time
t += dt
print "t="+str(t)

# infinite source (for the whole simulation)
source_term.t = t

# source only for a little time and then 0
# if t<=2.:
#     source_term.t = t # inject source all the time
# else:
#     source_term.t = 0.

# Assemble the RHS & apply BCs
b = assemble(a_rhs)
bcs.apply(A,b)
# Solve system
solve(A, u.vector(), b)

# Update solution, velocity and acceleration:
u0.assign(u1)
u1.assign(u)

# save the solution at two relevant positions:
exact_point = Point(np.array((xs,ys)))
u1value = u(exact_point)
sensor_x[0,tstep] = u1value

exact_point2 = Point(np.array((xs2,ys2)))
u1value = u(exact_point2)
sensor_x2[0,tstep] = u1value

tstep+=1

# SAVE ON FILE
l = []
l.append("solutions/frtime/centraldiff_Lsourcesin_size"+str(hmin)+"_T"+str(T)+"_c"+str(c_)+"_pdim"+str(pdim))
l.append('.mat')

filenamex = ''.join(l)

sio.savemat(filenamex, {'solution_x':sensor_x,'solution_x2':sensor_x2,'solution_x_coord':XS,'solution_y_coord':YS,'dt':dt,'T':T,\
'mesh_size':size,'hmin':hmin,'xlim':xlim,'ylim':ylim,'source_pos_x':xSource,'source_pos_y':ySource,'velocity':c_,})​

from dolfin import *
from mshr import *
import numpy  as np
import scipy.io as sio

# Velocity
c_ = 1.

# Create mesh
xlim, ylim = 1., 1.
size = 60
domain = Rectangle(Point(0., 0.), Point(xlim, ylim))
mesh = generate_mesh(domain,size)
hmin = mesh.hmin()
hmax = mesh.hmax()

# Order of the Finite Elements
pdim = 1

# Time Constants
t = 0.
CFL = 0.5
dt = CFL*hmax/(sqrt(2.)*c_)
print dt
T = 10.
nT = int(T/dt)+1
tvec = np.linspace(0,T,nT)

# Newmark Constants
beta_, gamma_ = 0.25, 0.5

# Boundaries
class Bound(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and ((abs(x[1] - 0.) < DOLFIN_EPS) or (abs(x[1] - ylim) < DOLFIN_EPS) or (abs(x[0] - 0.) < DOLFIN_EPS) or (abs(x[0] - xlim) < DOLFIN_EPS))

boundaries = FacetFunction("size_t", mesh, 0)

Bound().mark(boundaries, 1)

# Define Function space, trial and test functions
V = FunctionSpace(mesh, 'P', pdim)
a, phi = TrialFunction(V), TestFunction(V)

# Define Dirichlet BC
bcs = DirichletBC(V, Constant(0.), boundaries, 1)

# Fields from previous time step (displacement, velocity)
v0IC = Constant(0.)
u0 = interpolate(v0IC, V)
v0 = interpolate(v0IC, V)

# Define sensor locations - where I save the solution:
xs = xlim/4
ys = ylim/2

xs2 = xlim*0.85
ys2 = ylim*0.67

XS = [xs, xs2]
YS = [ys, ys2]

sensor_x = np.zeros((1,nT))
sensor_x2 = np.zeros((1,nT))

# Define the source location and expression
gauss = '(1.0/(2.*ppi*sigmax*sigmay) * exp(- pow((x[0]-mux),2)/(2.*pow(sigmax,2)) - pow((x[1]-muy),2)/(2.*pow(sigmay,2))))'

sigma_space = min(xlim,ylim)/200 # variance of the space gaussian
xSource = xlim/2
ySource = ylim/2

toneburst = 'sin(2*ppi*fc1*t)'
t = 0.
fc1 = 1.
source_term = Expression(toneburst+'*'+gauss, fc1=fc1, ppi=np.pi, mux=xSource, muy=ySource, sigmax=sigma_space, sigmay=sigma_space, t=t, degree=pdim)

# Compute a0 from u0
a0_lhs = -1/pow(c_,2) * a * phi * dx

A0 = assemble(a0_lhs)
b0 = assemble(a0_rhs)

bcs.apply(A0,b0)

a0 = Function(V)
solve (A0,a0.vector(),b0)

# Define LHS and RHS
def m_block(value):
return 1/pow(c_,2) * value * phi * dx
def k_block(value):

M = assemble(m_block(a))
K = assemble(k_block(a))
A = M + (pow(dt,2)*beta_) * K

a_rhs = - (k_block(u0) + dt*k_block(v0) + (pow(dt,2)*(0.5-beta_))*k_block(a0)) + source_term*phi*dx

a1 = Function(V)
u1_proj = Function(V)

tstep = 0
# Time loop
while t < T:

# Update time
t += dt
print "t="+str(t)
source_term.t = t
# if t<=2.:
#     source_term.t = t # inject source all the time
# else:
#     source_term.t = 0.

# Assemble the RHS & apply BCs
b = assemble(a_rhs)
bcs.apply(A,b)
# Solve system
solve(A, a1.vector(), b)

# Newmark method to compute displacement and velocity
u1 = u0 + dt*v0 + pow(dt,2)*((0.5-beta_)*a0 +beta_*a1)
v1 = v0 + dt*((1-gamma_)*a0 + gamma_*a1)

# Save solution to VTK format

u1_proj.vector()[:] = project(u1, V).vector()
# vtk_file_u << u1_proj

# Update solution, velocity and acceleration:
u0.assign(u1)
v0.assign(v1)
a0.assign(a1)

# save the solution at two relevant positions:
exact_point = Point(np.array((xs,ys)))
u1value = u1_proj(exact_point)
sensor_x[0,tstep] = u1value

exact_point2 = Point(np.array((xs2,ys2)))
u1value = u1_proj(exact_point2)
sensor_x2[0,tstep] = u1value

tstep+=1

# SAVE ON FILE
l = []
l.append("solutions/frtime/newmark_Lsourcesin_size"+str(hmin)+"_T"+str(T)+"_c"+str(c_)+"_pdim"+str(pdim))
l.append('.mat')

filenamex = ''.join(l)

sio.savemat(filenamex, {'solution_x':sensor_x,'solution_x2':sensor_x2,'solution_x_coord':XS,'solution_y_coord':YS,'dt':dt,'T':T,\
'mesh_size':size,'hmin':hmin,'xlim':xlim,'ylim':ylim,'source_pos_x':xSource,'source_pos_y':ySource,'velocity':c_,})​
Community: FEniCS Project