Steady state reached while gradients remain.

Question

Hello, I'm solving a diffusion equations, 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 project and examine mu I see that it is not flat.

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

Does anyone have an explanation for this please?

Thanks, MWE below.

from dolfin import *
from cbcpost 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)

pp = PostProcessor(dict(casedir="Results", clean_casedir=True))
pp.add_field(SolutionField("C", dict(save=True, save_as=["xdmf"])))
pp.add_field(SolutionField("mu", dict(save=True, save_as=["xdmf"])))

solution = {"C": lambda: project(C,V1), "mu": lambda: project(mu,V1)  } 

pp.update_all(solution, 0, 0)

while (iter< iter_max):
    iter +=1
    solver.solve()
    pp.update_all(solution, t, iter)
    t += dt(0)
    C_old.vector()[:] = C.vector()

pp.finalize_all()

Comments

What is the magnitude of the solution gradient? Does the effect vanish when you neglect p?

---

Yes the effect vanishes when I neglect p.
The maximal magnitude of grad(C) is -10.

I noticed that the shape of mu is the second derivative of the function p. Therefore I wonder if there is a subtlety with the differentiation of the form or something along those lines?

(Forgive my ascii art). Mu looks like:

             /\
---------   /   ---------
          \/

where the flat line is .05, and the extrema are +- 0.0013

Thanks

---

UPDATE:

I implemented this model in the COMSOL Multiphysics (commercial) and it does not have this artefact. Does anyone have any possible explanations?

Thanks