DOLFIN

Spatial dimension in assembling matrices

Asked by Derek Thomas

I'm trying to solve the wave equation in FEniCS. I've adapted the Poisson demos and the advection demo to try to integrate it. When I run this script I get an error "Need to know the spatial dimension to compute the shape of derivatives."

The error is associated with the assembly of the matrix A. I'm not sure how to go about debugging this error. Thanks,

Derek

from dolfin import *
from numpy import pi
# Create mesh and define function space
mesh = UnitSquare(100,100)
V = FunctionSpace(mesh, "CG", 1)
W = V*V
c0 = 1.0

# Parameters
T = 5.0
dt = 0.1
t = dt

# Define Dirichlet boundary (x = 0 or x = 1
def boundary(x):
    return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS or x[1] < DOLFIN_EPS or x[1] > 1.0 - DOLFIN_EPS

# Define boundary condition
# (p0, g0) = Function(W)
p0 = Constant(0.0)
g0 = Constant(0.0)
bc = DirichletBC(W, (p0,g0), boundary)

# Define variational problem
(p, g) = TrialFunction(W)
(v, w) = TestFunction(W)

p_mid = 0.5*(p+p0)
g_mid = 0.5*(g+g0)

f = Constant(0.0)

a = v*(g-g0)*dx -w*(p-p0)*dx + dt * ( dot(grad(p_mid),grad(v))*dx + w*g_mid*dx)

L = f*dx

p1 = Point(0.1,0.1)

A = assemble(a)
bc.apply(A)

solver = LUSolver(A)
solver.parameters["reuse_factorization"] = True
out_file = File("test.pvd")

while t<T:
    b = assemble(L)
    bc.apply(b)
    delta = PointSource(V, p1, sin(2*pi*t))
    delta.apply(b)

    solver.solve(p.vector(), b)
    p0 = p
    plot(p)
    out_file<<(p,t)

    t += dt

interactive()

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Praveen C
Solved:
Last query:
Last reply:
Revision history for this message
Praveen C (cpraveen) said : 		#1

I think there are too many mistakes in your code.

You must define p0, g0 as of type Function, since these are solution values at previous time step.

u = Function(W)
uold = Function(W)
p0 = uold[0]
g0 = uold[1]

Then you also need the TrialFunction and TestFunction.

Write down your weak formulation like B=0

 B = ....

then do

a = lhs(B)
L = rhs(B)

When you solve, you are passing only "p" but your solution is (p,g) in your notation.

Your L = f*dx is missing the test function.

praveen

Revision history for this message
Derek Thomas (derekcthomas) said : 		#2

Thanks for your comments Praveen, I apologize for posting code with so many mistakes. Following your comments, I know have a script that runs up to the solver.solve point. Now I get an error that the argument to the LUSolver_solve method is of the wrong type. I originally had solver.solve(p.vector(), b), I understand that I was neglecting to pass in the g portion of the solution. How would I go about doing that? If I use solver.solve((p.vector(), g.vector(),b), I get an error that an 'Indexed' object has no attribute 'vector'. I've tried solver.solve((p,g), b), but get a type error. What is the correct method for converting my p and g functions to a vector that can be passed to the solver? Thanks,

Derek

Revision history for this message
Best Praveen C (cpraveen) said : 		#3

Maybe something like should work.

v = TrialFunction(W)
w = TestFunction(W)

u = Function(W)
uold = Function(W)

# set initial condition
...
uold.assign(u)

p0 = uold[0]
g0 = uold[1]

pmid = 0.5*(v[0] + p0)
gmid = 0.5*(v[1] + g0)

B = (1/dt)*(v[0] - p0)*w[0] + (1/dt)*(v[1] - g0)*w[1] + ....

etc.

solver.solve(u.vector(), b)

Revision history for this message
Derek Thomas (derekcthomas) said : 		#4

Thanks Praveen C, that solved my question.