2
$\begingroup$

This presentation by the Imperial College in London has a nice example in it, on page 8, Burgers Equations. The first part of their code reads like this:

from dolfin import*

n = 30
mesh = UnitInterval(n)
V = FunctionSpace(mesh, "CG",2)
ic = project(Expression("sin(2*pi*x[0])"), V)
u = Function(ic)
u_next = Function(V)
v= TestFunction(V)
nu = Constant(0.0001)
timestep = Constant(1.0/n)

F1 = ((u_next - u)/timestep*v)*dx
F2 = (u_next*grad(u_next)*v)*dx
F3 = (nu*grad(v)*grad(u_next))*dx

On my system, this code fails, giving shape mismatch errors for the sums F2 and F3. I don't understand why.

Improve this question
$\endgroup$
4
  • $\begingroup$ Are you sure that is the correct equation? u grad(u) = 0? That's a vector-valued equation for a scalar variable. $\endgroup$
    – Anders Logg
    May 27, 2013 at 17:09
  • $\begingroup$ sorry, I'm being really sloppy here. what I mean is $\vec{u} ( \nabla \cdot \vec{u}) = 0$. I'll adjust it. $\endgroup$
    – seb
    May 27, 2013 at 17:23
  • $\begingroup$ @AndersLogg ok, I realize that this doesn't make too much sense. I need to put more thought into it. $\endgroup$
    – seb
    May 27, 2013 at 17:26
  • $\begingroup$ Read UFL manual. $\endgroup$
    – Jan Blechta
    May 28, 2013 at 17:47

4 Answers 4

Reset to default
4
$\begingroup$

So the problem is that FEniCS doesn't recognize that in 1D, u and grad(u) have the same shape. (You could file an issue about this on the Dolfin issue tracker.)

If you replace grad by div, your code works. (The general form of the inviscid Burgers equation is $\mathop{\mathrm{div}}(u^2) = f$.)

Improve this answer
$\endgroup$
2
  • $\begingroup$ I guess this is not issue. u is rank 0, shape () but grad(u) is rank 1, shape (1,). This is consistent and is consistetly solved by using div(u) instead of grad(u). $\endgroup$
    – Jan Blechta
    May 28, 2013 at 17:42
  • $\begingroup$ On the other hand F3 = (nu*dot(grad(v), grad(u_next)))*dx is valid in arbitrary dimension. $\endgroup$
    – Jan Blechta
    May 28, 2013 at 17:47
3
$\begingroup$

If your domain is one dimensional you should use:

mesh = UnitIntervalMesh(10)

For partial derivative in x direction use:

u.dx(0)
Improve this answer
$\endgroup$
2
$\begingroup$

The error is in the line

a = inner(inner(nabla_grad(u),u),v)*dx

Here you take the inner product between the vector nabla_grad(u) and the scalar u. That does not make sense.

Improve this answer
$\endgroup$
1
  • $\begingroup$ ok, so how would I tell fenics to take the derivative with respect to x instead of the gradient?, i.e. a = (d/dx u)*uvdx ? $\endgroup$
    – seb
    May 27, 2013 at 17:20
0
$\begingroup$

Thanks for all the tips and hints. So after some dist-upgrading, debugging, manual reading, I finally got this script running. I'm including my solution here, just for the sake of completeness. Note that this code models the inviscid Burgers equations in one dimension.

from dolfin import*

n = 100
mesh = UnitIntervalMesh(n)
#mesh = IntervalMesh(500,-1.,10.)
V = FunctionSpace(mesh,'CG',1)
u = Function(V)
u0 = Expression("1./(sqrt(2*pi)*0.1) * exp(-0.5*pow((x[0]-0.5),2)/0.01)")
u.interpolate(u0)

u_next = Function(V)
v = TestFunction(V)

nu = Constant(0.01)
timestep = Constant(0.01)

F = ((u_next-u)/timestep*v
     + inner(u_next.dx(0)*u_next, v))*dx
bc = DirichletBC(V, 0.0, "on_boundary")

t = 0.0
end =10
while(t <= end):
    solve(F==0, u_next, bc)
    u.assign(u_next)
    t+=float(timestep)
    plot(u)
interactive()
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.