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,
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 - pi/2 > 0 ? exp(-D * pow(sin(x), 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)) * cos(x) + pow(sin(x), 2)', degree=2, i=i) #coefficient for zero order derivative term of u coeff1 = Expression('( sin(x)*(i-cos(x)) - cos(x)/(sin(x)*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!