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:
where i, D are constants,
with BC
and an IC
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 + \
dt*coeff1*grad(u)*v*dx - \
dt*(1/(2*D))*dot(grad(u), grad(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!