Wrong answer unless I use an 'auxiliary' variable?

3 months ago by
I'm solving a simple diffusive steady state equation for c:
and p(x) is some other function (for appears irrelevant as long as its order is > 2...)
The solution should be  $\mu=constant$μ=constant , but the solution to the MWE below gives variation in \mu. Meanwhile, if I define a new unknown,  $\mu2=\mu$μ2=μ , and use it in the diffusion equation, I get the correct answer (uncomment the line below). Can anyone please kshed some light on this? I've played with quadrature degree, element order, tolerances, etc to no avail.

from dolfin import *

mesh = IntervalMesh(100,0,1)

P1 = FiniteElement("Lagrange",mesh.ufl_cell(), 1)
P = P1*P1
V = FunctionSpace(mesh, P)

U = Function(V)
C,mu2 = split(U)
dU = TrialFunction(V)
test_C, test_mu = TestFunctions(V)

x = SpatialCoordinate(mesh)
p = x[0]**3

mu = C-p

F_C      = inner(grad(mu), grad(test_C))*dx  #correct: mu is flat
#F_C      = inner(grad(mu2), grad(test_C))*dx  #Wrong: mu is not flat
F_mu = (mu2-mu)*test_mu*dx
F = F_C+F_mu

bc = DirichletBC(V.sub(0),0, CompiledSubDomain("(std::abs(x[0]) < DOLFIN_EPS) && on_boundary"))

file_c = File("c.pvd")
file_mu = File("mu.pvd")

problem = NonlinearVariationalProblem(F, U, bc, derivative(F, U, dU))
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters
prm["newton_solver"]["absolute_tolerance"] = 1e-10
prm["newton_solver"]["relative_tolerance"] = 1e-10

file_c << U.split()[0]
file_mu << U.split()[1]​


Community: FEniCS Project

1 Answer

3 months ago by
In the first case, you substitute in the first equation for \(\mu\) from the second one. This gives you
\[ \left( \nabla (c - p(x), \nabla \bar{c}) \right) = 0, \]
which is a correct equation for the unknown \( c \). But you also add the term \( (\mu - c + p(x), \bar{\mu}) \), i.e., you have
\[ \left( \nabla (c - p(x), \nabla \bar{c}) \right) + (\mu - c + p(x), \bar{\mu}) = 0. \]
The first term is satisfied as a equation for \(c\) and from the second term then follows that \(\mu = 0\), which can be checked more accurately by rising the element order.

In the second case, you do not substitute and solve the system as a whole for \((c, \mu) \)
\[ \left( \nabla \mu, \nabla \bar{c}) \right) + (\mu - c + p(x), \bar{\mu}) = 0, \]
giving you the constant solution.
Thanks for the reply but this doesn't really help. I include the second term to make it convenient to switch between using mu2 or not, but it doesn't make any difference to the solution. Either way, the solution without the secondary variable (i.e.: as written above) does not produce \mu = a constant as it should. Higher order elements get better, but still appear to be dependant on the form of p(x). E.g. if
$p=.5\cdot\left(1+tanh\left(\frac{\left(.5-x\left[0\right]\right)}{.05}\right)\right)$p=.5·(1+tanh((.5x[0]).05 ))
then higher order elements get better, but still not perfect. Meanwhile the projection gets it bang on (< 10^-10 variation).
Therefore, I am hazarding a guess that this is related to quadrature / sampling of the function p....

Incidentally, the same model in other FEM software produces \mu = constant without the second variable at all. It is only FEniCS that I have encountered this.
written 3 months ago by Mike Welland  
The first case, i.e.,
\[ \left( \nabla (c - p(x)), \nabla \bar{c} \right) + \left( \mu - c + p(x), \bar{\mu} \right) = 0, \]
is equivalent to solving just the equation for \(c\)
\[ \left( \nabla (c - p(x)), \nabla \bar{c} \right) = 0, \]
and then finding \( \mu \) as the projection of \( c-p(x) \) onto the function space, which agrees to the second term.
from fenics import *

mesh = UnitIntervalMesh(100)

V = FunctionSpace(mesh, 'P', 1)

c = Function(V)
test_c = TestFunction(V)

x = SpatialCoordinate(mesh)
# p = x[0]**3
p = 0.5*(1.0 + tanh((0.5-x[0])/0.05))

F = inner(grad(c-p), grad(test_c))*dx

bc = DirichletBC(V, 0.0, 'near(x[0], 0.0) && on_boundary')

problem = NonlinearVariationalProblem(F, c, bc, derivative(F, c))
solver = NonlinearVariationalSolver(problem)

mu = project(c-p, V)
print mu.vector().get_local()​
As is written in the FEniCS tutorial, "projection results in approximate values at the nodes". That's why you don't get a perfect constant solution.

In the second case, you compute for \( \mu \) (first term) that satisfies \( \mu = c-p(x) \) (second term)
from fenics import *

mesh = UnitIntervalMesh(100)

P = FiniteElement('Lagrange', mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, MixedElement([P, P]))

u = Function(V)
c, mu = split(u)
test_c, test_mu = TestFunctions(V)

x = SpatialCoordinate(mesh)
# p = x[0]**3
p = 0.5*(1.0 + tanh((0.5-x[0])/0.05))

# F_c = inner(grad(c-p), grad(test_c))*dx # CASE 1
F_c  = inner(grad(mu), grad(test_c))*dx # CASE 2
F_mu = (mu-c+p)*test_mu*dx
F = F_c + F_mu

bc = DirichletBC(V.sub(0), 0.0, 'near(x[0], 0.0) && on_boundary')

problem = NonlinearVariationalProblem(F, u, bc, derivative(F, u))
solver = NonlinearVariationalSolver(problem)

mu = u.sub(1, deepcopy=True)
print mu.vector().get_local()​
written 12 weeks ago by Adam Janecka  
Thanks Adam, so are you saying that the solution is actually the same, but it is just my visualization of the solution (projection of mu) that is the difference..?
written 11 weeks ago by Mike Welland  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »