### Wrong answer unless I use an 'auxiliary' variable?

157

views

0

I'm solving a simple diffusive steady state equation for c:

$\nabla\cdot\nabla\mu=0$∇·∇

where

$\mu=c-p\left(x\right)$

and p(x) is some other function (for appears irrelevant as long as its order is > 2...)

The solution should be $\mu=constant$

$\nabla\cdot\nabla\mu=0$∇·∇

`μ`=0where

$\mu=c-p\left(x\right)$

`μ`=`c`−`p`(`x`)and p(x) is some other function (for appears irrelevant as long as its order is > 2...)

The solution should be $\mu=constant$

`μ`=`c``o``n``s``t``a``n``t`, 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
#info(prm,True)
solver.solve()
file_c << U.split()[0]
file_mu << U.split()[1]
```

Thanks

Community: FEniCS Project

### 1 Answer

0

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.

\[ \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.

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.

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

\[ \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)
solver.solve()
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)
solver.solve()
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.

$p=.5\cdot\left(1+tanh\left(\frac{\left(.5-x\left[0\right]\right)}{.05}\right)\right)$

p=.5·(1+tanh((.5−x[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.