# Variation formulation for a second order PDE with first order derivative term

Hi, I am a new user to FEniCS project and intending to solve a second order 2d non-linear PDE. Before that, I am trying on a reduced case, which is a 1d case as following:

\frac{\partial u(x, t)}{\partial t}=-\frac{1}{sin(x)}\frac{\partial}{\partial x}[sin^2(x)(i-cos(x))u-Dsin(x)\frac{\partial u}{\partial x}]

where i, D are constants,
with BC

\frac{\partial u}{\partial x}|_{x=0, \pi} =0

and an IC

u(t=0) = u_0(x)

After reading the several examples from the tutorial, I try to convert this to variational formulation by expanding the right-hand-side and using the backward method in time steps. I came up with the following codes

from fenics import *
import scipy.integrate as scipyintegrate

T = 4
num_steps = 10
dt = T / num_steps

i = Constant(1.5)
D = Constant(60)
denorm, _ = scipyintegrate.quad(lambda x: sin(x) * exp(-D*sin(x)*sin(x)), 0, pi/2) #use for initial condition

mesh = IntervalMesh(200, 0, pi)
V = FunctionSpace(mesh, 'CG', 1)

# initial condition
u_0 = Expression('x[0] - pi/2 > 0 ? exp(-D * pow(sin(x[0]), 2)) / denorm : 0', \
degree=2, D=D, denorm=denorm)
u_n = interpolate(u_0, V)

u = TrialFunction(V)
v = TestFunction(V)

coeff0 = Expression('2 * (i-cos(x[0])) * cos(x[0]) + pow(sin(x[0]), 2)', degree=2, i=i)  #coefficient for zero order derivative term of u

coeff1 = Expression('( sin(x[0])*(i-cos(x[0])) - cos(x[0])/(sin(x[0])*2*D) )', degree=2, i=i, D=D) # #coefficient for first order derivative term of u (du / dt)

# coeff2 = Expression('1 / (2 * D)', degree=0, D=D)   #coefficient for second order derivative term of u (d2u/dt2)

F = u*v*dx + \
dt*coeff0*u*v*dx + \
u_n*v*dx
a, L = lhs(F), rhs(F)

vtkfile = File('solution/solution.pvd')

u = Function(V)
t = 0
for n in range(num_steps):
t += dt
solve(a == L, u)
vtkfile << (u, t)

u_n.assign(u)


I think I have some wrong with the first order derivative term and it has the error:

ufl.log.UFLException: Can only integrate scalar expressions. The integrand is a tensor expression with value shape (1,) and free indices with labels ().


I wonder how I can fix this problem and make it work, any suggestion is appreciated. Thanks in advance!

The gradient produces a vector of shape (1, ). Therefore, you should take the first component of it:

  dt*coeff1*grad(u)[0]*v*dx


You could also use u.dx(0) instead of grad(u)[0].

1 Like