Wave equation in 2D solved in the frequency domain


81
views
0
6 weeks ago by

Hello,

I write here a new question related to my previous question https://www.allanswered.com/post/vnebp/pde-complex-numbers/.

I am interested in solving the wave equation in 2D, not in the time domain, but rather in the frequency domain. The equation we're interested in is the following:  $\nabla^2u\left(x,t\right)-\frac{1}{c^2}\frac{\partial^2}{\partial t^2}u\left(x,t\right)=f\left(x,t\right)$2u(x,t)1c2 2t2 u(x,t)=ƒ (x,t) with zero Dirichlet initial conditions and boundary conditions. The source term  $f\left(x,t\right)=h\left(x\right)g\left(t\right)$ƒ (x,t)=h(x)g(t), where  $g\left(t\right)=sin\left(2\pi f_1t\right)+3sin\left(2\pi f_2t\right)$g(t)=sin(2πƒ 1t)+3sin(2πƒ 2t). By computing the fft of g(t), we can solve the problem in the frequency domain only for the most relevant modes (only 2 for this simple case).

The code below is well commented and speaks for itself.

The problem is 
1) When I try to reconstruct the solution    $u\left(x,t\right)=\sum_ke^{-i\frac{2\pi}{T}k}u_k$u(x,t)=kei2πT kuk   I don't retrieve the same solution I obtain solving the problem in time.
2) What I find strange is that the real part of solution using the fft(source)_k is equal to the imaginary part of the solution using the conjugate of fft(source)_k and viceversa. I don't think it should be like this. If any, I'd expect the two solutions to be complex conjugates.

Has anyone encountered similar problems? Or, as anyone solve the wave equation in the frequency domain in a different way?

Thank you very much.


from dolfin import *
from mshr import *
import numpy  as np



# Mesh details 
xlim, ylim = 2., 4.
size = 60
domain = Rectangle(Point(0., 0.), Point(xlim, ylim))


# Sensor location - where we'll save the solution
xs = 1.
ys = 1.


# Order of the Finite Elements 
pdim = 1

# Create mesh
mesh = generate_mesh(domain,size)
hmin = mesh.hmin()


# Velocity - CFL number - Time constants
c_ = 5.
CFL = 0.5

dt = CFL*hmin/(sqrt(2.)*c_)
print "dt = "+str(dt)
T = 3#100*dt
nT = int(T/dt)+1
tvec = np.linspace(0,3,nT)



####################################################################
#                         FE Space + BC                            #
#################################################################### 
# Define Function space, trial and test functions
V2 = VectorFunctionSpace(mesh, 'CG', pdim, dim=2)
(u_r, u_i) = TrialFunctions(V2)
(phi_r, phi_i) = TestFunctions(V2)


# 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) 

bcs = DirichletBC(V2, Constant((0.,0.)), boundaries, 1)

 
####################################################################
#                         Source definition                        #
#################################################################### 
# fft(source_time_vec)_k * h(x), where h(x) is a Gaussian centered in (xSource,ySource)
freq_gauss = 'fr_kk * (1.0/(2.*pi*sigmax*sigmay) * exp(- pow((x[0]-mux),2)/(2.*pow(sigmax,2)) - pow((x[1]-muy),2)/(2.*pow(sigmay,2))))' 

sigma_space = 0.05 # variance of the space gaussian
xSource = 1.
ySource = 2.

fr_kk = 0.

source_termR = Expression(freq_gauss, fr_kk = fr_kk, mux=xSource, muy=ySource, sigmax=sigma_space, sigmay=sigma_space, degree=pdim)
source_termI = Expression(freq_gauss, fr_kk = fr_kk, mux=xSource, muy=ySource, sigmax=sigma_space, sigmay=sigma_space, degree=pdim)


a1 = 1.
a2 = 3.
fc1 = 4.
fc2 = 5.

def source_time_function(t):
    eval_source = a1*sin(2*pi*fc1*t) + a2*sin(2*pi*fc2*t) # start with a simpler soruce term 
    return eval_source

source_term_time = np.zeros((nT))
for tt in range(0,nT):
    source_term_time[tt] = source_time_function(tvec[tt])
 
source_term_freq = np.fft.fft(source_term_time)


# Select only the most significant modes
def indices(a, func):
    return [i for (i, val) in enumerate(a) if func(val)]

treshold = 0.1
inds = indices(np.abs(source_term_freq)/nT, lambda x: x>treshold) 
bigK = len(inds)
# Question 1  - why should I divide by nT? Values are very big otherwise. It really depends on the source function 



####################################################################
#                         Initializations                          #
#################################################################### 

# Save initial displacement (actually pressure here) and acceleration  
vtk_file_u_r = File("result_source/p"+str(pdim)+"_"+str(i)+"_real.pvd")
vtk_file_u_i = File("result_source/p"+str(pdim)+"_"+str(i)+"_imag.pvd")
vtk_file_source = File("result_source/p"+str(pdim)+"_"+str(i)+"_source.pvd")


# initialize the vectors where we will save the solution (real and imag part) for all the modes 
sensor_xR = np.zeros((bigK)) 
sensor_xI = np.zeros((bigK))

u = Function(V2)
u_real = Function(V2)
u_imag = Function(V2)


####################################################################
#          Define LHS and RHS & loop over the modes                #
#################################################################### 
def A_block(value_trial,value_test,omegaa):
    return (pow(omegaa/c_,2) * inner(value_trial, value_test) - inner(grad(value_trial), grad(value_test))) * dx

def b_block(ssource,value_test):
    return inner(ssource, value_test) * dx



for kk in range(0,int(bigK/2)): # we can loop over half of the frequencies (the others will simply need a minus!)
    print inds[kk]

    #################################################################################
    #                Frist run + SAVE SOLUTION AT SENSOR POINT                      #
    ################################################################################# 
    source_kk = source_term_freq[inds[kk]]/nT
    # Question 2 - why do I need to divide by nT? Or don't I?
    print source_kk

    source_termR.fr_kk = source_kk.real
    source_termI.fr_kk = source_kk.imag

    omega_ = 2*np.pi/T*inds[kk]
    
    A_tot = A_block(u_r,phi_r,omega_) - A_block(u_i,phi_i,omega_) + A_block(u_r,phi_i,omega_) + A_block(u_i,phi_r,omega_)
    b_tot = b_block(source_termR,phi_r) + b_block(source_termI,phi_i)

    A = assemble(A_tot)
    b = assemble(b_tot)
    bcs.apply(A,b)
    solve(A, u.vector(), b)

    (u_real, u_imag) = u.split(True)
    vtk_file_u_r << u_real
    vtk_file_u_i << u_imag

    # Save the value of u 
    exact_point = Point(np.array((xs,ys))) 
    print u_real(exact_point)
    print u_imag(exact_point)
    
    sensor_xR[kk] = u_real(exact_point)
    sensor_xI[kk] = u_imag(exact_point)
 
   
    #################################################################################
    #                Second run + SAVE SOLUTION AT SENSOR POINTS                     #
    ################################################################################# 
    # This time the rhs is the complex conjugate, omega is the same because is squared
    
    source_termI.fr_kk = -source_kk.imag # equivalent to : source_term_freq[inds[bigK-kk-1]]


    b = assemble(b_tot)
    bcs.apply(A,b)
    solve(A, u.vector(), b)

    (u_real, u_imag) = u.split(True)
    vtk_file_u_r << u_real
    vtk_file_u_i << u_imag

    # Save the value of u 
    exact_point = Point(np.array((xs,ys))) 
    print u_real(exact_point)
    print u_imag(exact_point)
    
    sensor_xR[int(bigK/2)+kk] = u_real(exact_point)
    sensor_xI[int(bigK/2)+kk] = u_imag(exact_point)

    # Question 3 - why can't I use u_r and u_i (trial functions) and I need to define a new function on V2? (I get error using la_solve instead...)
    # Question 4 - why the real part of solution using the fft(source)_k is equal to the imaginary part of the solution using the conjugate of fft(source)_k and viceversa?
Community: FEniCS Project
Please login to add an answer/comment or follow this question.

Similar posts:
Search »