- removed milestone
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)
-
-
- changed status to invalid
This is too complex mathematical issue and can hardly help identifying an eventual bug. We can't base bug reports on comparisons with commercial software.
- Log in to comment
Removing milestone: 1.7 (automated comment)