Update expression and time-dependent assembly

Question

Hello, I'm trying to solve the time-dependent elastic-acousitic equation

rho (lambda+2mu) grad(div(u)) - mu curl(curl(u)) = f

where u is 2D vector and represents the displacement.

f is the time-dependent source term defined as follows:
f = h(x,y)*g(t), where h(x,y) is a 2D Gaussian centered in (0.5,0.2) and g(t) is a Morlet wavelet (oscillatory behaviour between -1 and 1). Moreover, the two components of the f vector are the same, i.e. f_x = f_y

I'm solving this problem for zero initial conditions using the Newmark method for time integration.

I've used the incompressible Navier-Stokes as tutorial https://fenicsproject.org/olddocs/dolfin/1.5.0/python/demo/documented/navier-stokes/python/documentation.html, but the solution is monotonically increasing instead of oscillatory (as the source term).

Can anyone spot why is so? It'd be very helpful! Also, is it possible to save the source term in a vtk file? I've tried to update the time, but the way I'm doing it now doesn't seem to work (the saved file is fixed to the first value of the source)

Here is my code:

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



# Physics Constants for low alloy steel 
rho_ = Constant(7.8)            # Mg*m^(-3)
E = Constant(205)       # GPa
nu = Constant(0.3)      # constant
lambda_ = Constant((E*nu)/((1.0+nu)*(1.0-2.0*nu))) #~109.6 Gpa
mu_ = Constant(E/(2.0*(1.0+nu))) #~73.1 GPa

# Newmark Constants
beta_ = 0.25
gamma_ = 0.5

# Time Constants
dt = 0.0025
T = 1.


# Create mesh & Define Function space
domain = Rectangle(Point(0., 0.), Point(1., 1.))
mesh = generate_mesh(domain, 16)

V = VectorFunctionSpace(mesh, 'CG', 1, dim=2)
u, phi = TrialFunction(V), TestFunction(V)


# Define natural boundary conditions: <u,n>=0 on each boundary (also internal)

def horizontal_boundary(x, on_boundary):
    return on_boundary and  (abs(x[1] - 1.) < DOLFIN_EPS or  abs(x[1] - 0.) < DOLFIN_EPS) 

def vertical_boundary(x, on_boundary):
    return on_boundary and  (abs(x[0] - 1.) < DOLFIN_EPS or  abs(x[0] - 0.) < DOLFIN_EPS)

ux_0 = Constant(0.0)
uy_0 = Constant(0.0) 

bcH = DirichletBC(V.sub(1), uy_0, horizontal_boundary)
bcV = DirichletBC(V.sub(0), ux_0, vertical_boundary)

bcs = [bcH, bcV]



# Initial conditions - all zero

zeroIC = Constant((0.0,0.0))
u0 = interpolate(zeroIC, V)
v0 = interpolate(zeroIC, V)
a0 = interpolate(zeroIC, V)

u1 = Function(V)

# Define time-dependent source
fM = 5. # frequency
s = 6./(2*pi*fM)
tD = 0.6 # shift in time to start with an approximately zero source
wavelet =  '(exp(-pow(t-tD,2)/(2.*pow(s,2))) * cos(2.*pi*fM*(t-tD)))'
gauss = '(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))))'
source_expr=Expression((wavelet+'*'+gauss, wavelet+'*'+gauss), mux=0.5, muy=0.2, sigmax=0.3, sigmay=0.3, fM = fM, s=s, tD =tD, t=0.0, degree=1)

source = interpolate(source_expr,V)


# Matrices
s = lambda_*div(u)*div(phi)*dx + mu_*inner( grad(u)+transpose(grad(u)), grad(phi) )*dx 

m = rho_ * inner(u,phi) * dx


# Assemble mass matrix and right hand side
A = assemble(m + beta_*pow(dt,2)*s)
L = rho_ * inner( (u0 + dt*v0 + (0.5-beta_)*pow(dt,2)*a0) ,phi) * dx + beta_*pow(dt,2)*inner(source,phi)*dx



# Time loop
t = 0.0
vtk_file = File("displacement.pvd")


vtk_file2 = File("source.pvd")

while t <= T:

    # Update the source that is time depedent and save it
    source_expr.t = t
    if t <= 0.3 : 
        vtk_file2 << source # THIS SEEMS NOT TO WORK?

    # Assemble the RHS
    b = assemble(L)

    # Axpply BCs
    [bc.apply(A) for bc in bcs] 

    # Solve system for displacement
    solve(A, u1.vector(), b)

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


    # Save solution to VTK format
    vtk_file << u1

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

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

Answer 1

You have an unnecessary line

source = interpolate(source_expr,V)

This will make a "static" function from source_expr as it is right then, but the source function will not be updated when you update source_expr (unless you run this line again inside the time loop).

Just use source_expr instead of source and it will probably work, I see you are updating t in source_expression. There is no real need to interpolate source_expr to a function before putting it inside the form

Comments

Thank a lot! It works :)

I thought I needed the interpolation over V to do the inner product of the source term with the test function.