Incomplete steady state solution, not reproducible in commercial software.

This issue was originally posed as a question at here

I'm encountering an issue with FEniCS that I am unable to reproduce in commercial software (Comsol Multiphysics).

I'm solving a diffusion equation, with the flux resulting from a gradient of a form, mu. mu is defined as a function of the unknown variable, c, and another variable p.

The system reaches steady state (implying that the gradient of mu is constant and zero for these boundary conditions), but when I examine mu I see that it is not flat. Meanwhile, if I project mu onto another unknown (mu2 in the second example), then the gradient is eliminated.

I haven't been able to solve this problem by adjusting solver parameters.

MWE of lingering gradient

from dolfin import *

t, dt, t_final = [0, Constant(.01), 1.]
mesh = IntervalMesh(100,0,1)

V1 = FunctionSpace(mesh, "CG", 1)
C = Function(V1)
dU = TrialFunction(V1)
test_C = TestFunction(V1)
C_old = Function(V1)

x = SpatialCoordinate(mesh)
p = .5*(1.+tanh((-x[0]+.5)/.05))

mu = -( -C+p )

J        = -grad( mu )
F_C      = inner(J, grad(test_C))*dx 
Ft_C  = (C-C_old)*test_C*dx
F_C = Ft_C - dt*F_C

C.interpolate(project(((1-p)*.1+p),V1))

t = t+dt(0)
iter, iter_max = 0, 100

C_old.vector()[:] = C.vector()


problem = NonlinearVariationalProblem(F_C, C, [], derivative(F_C, C, dU))
solver = NonlinearVariationalSolver(problem)

while (iter< iter_max):
    iter +=1
    solver.solve()
    #pp.update_all(solution, t, iter)

    t += dt(0)
    C_old.vector()[:] = C.vector()
plot(mu, interactive =True)

MWE projecting mu onto a separate variable, mu2, resulting in correct behaviour.

from dolfin import *

t, dt, t_final = [0, Constant(.01), 1.]
mesh = IntervalMesh(100,0,1)

V1 = FunctionSpace(mesh, "CG", 1)
V = V1*V1
U = Function(V)
C,mu2 = split(U)
dU = TrialFunction(V)
test_C, test_mu2 = TestFunctions(V)

U_old = Function(V)
C_old, mu2_old = split(U_old)

x = SpatialCoordinate(mesh)
p = .5*(1.+tanh((-x[0]+.5)/.05))

mu = -( -C+p )

F_mu2 = (mu2-mu)*test_mu2*dx

J = -grad(mu2)

F_C = inner(J, grad(test_C))*dx 
Ft_C = (C-C_old)*test_C*dx
F_C = Ft_C - dt*F_C

F = F_C+F_mu2

assign(U.sub(0), project(((1-p)*.1+p),V1) )
assign(U.sub(1), project(Constant(0.0),V1))

t = t+dt(0)
iter, iter_max = 0, 100

U_old.vector()[:] = U.vector()


problem = NonlinearVariationalProblem(F, U, [], derivative(F, U, dU))
solver = NonlinearVariationalSolver(problem)

while (iter< iter_max):
    iter +=1
    solver.solve()
    #pp.update_all(solution, t, iter)

    t += dt(0)
    U_old.vector()[:] = U.vector()
plot(mu, interactive =True)

Comments (2)