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.

$\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$ Commented 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
    Commented 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
    Commented May 27, 2013 at 17:26
  • $\begingroup$ Read UFL manual. $\endgroup$ Commented May 28, 2013 at 17:47

4 Answers 4

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$.)

$\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$ Commented 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$ Commented 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)
$\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.

$\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
    Commented 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()
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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