Diffusion on symmetric problem becomes asymmetric.

Question

Hi.

Consider a simple heat diffusion problem in a domain Ω=[0,1] consisting of three subdomains
Ω_0 = [0 , 0.4)
Ω_1 = [0.4 , 0.6]
Ω_2 = (0.6 , 1]

Diffusion coefficient k is piecewise constant
k = 1 in Ω_0
k = 0.1 in Ω_1
k = 1 in Ω_2

Initially the temperature u is set to
u = 1 in Ω_0
u = 0 in Ω_1
u = 1 in Ω_2

On this link you will find plots of the initial values and a plot of the FEniCS solution at first time step.
Images and python file
The problem is that the initial values are symmetric but the diffusion is asymmetric.

We had some problems with setting the initial values too. We need to use 0<= x <= 0.4 to define Ω_0, such that our initial values is u = 1 in Ω_0 = [0 , 0.4) and thereby symmetric. If we set 0<= x < 0.4 the initial condition becomes asymmetric!

We have tried a lot of different ways to define the initial values, but cannot get this problem to be symmetric from start to end.

""" Heat diffusion into an element with different diffusion coefficient"""

from dolfin import *
import time

# set parameters
t = 0.0;            #Start time
dt = 0.005;           #Step size
T = 0.05;            #Total simulation time


# define mesh
mesh = UnitIntervalMesh(30)

# define function space
V = FunctionSpace(mesh, "CG", 1)
V0 = FunctionSpace(mesh, 'DG', 0)

# define functions
c1 = Function(V)
c_trial = TrialFunction(V)
phi = TestFunction(V)


# left part
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return ( (x[0] <= 0.4) )
left = Left()

# right part
class Right(SubDomain):
    def inside(self, x, on_boundary):
        return ( (x[0] > 0.6) )
right = Right()


# Initialize mesh function for interior domains
domains = CellFunction("size_t", mesh)
domains.set_all(1)
left.mark(domains, 0)
right.mark(domains, 2)



# Define new measures associated with the interior domains
dxx = Measure("dx")[domains]

# set initial condition
c0_trial = TrialFunction(V0)
phi0 = TestFunction(V0)
c0 = Function(V0)


a = c0_trial*phi0*dx
l = phi0*dxx(0)+phi0*dxx(2)
solve(a == l, c0)
c0 = interpolate(c0,V)


# define form
F = (1/dt)*(c_trial-c0)*phi*dx \
    + dot(grad(c_trial), grad(phi))*dxx(0) \
    + 0.2*dot(grad(c_trial), grad(phi))*dxx(1) \
    + dot(grad(c_trial), grad(phi))*dxx(2)


a = lhs(F)
l = rhs(F)

# preassembly
A = assemble(a)
b = None

plot1 = plot(c0); time.sleep(5)
plot1.write_png("initial.png")

while t < T:

    b = assemble(l, tensor=b)
    solve(A, c1.vector(), b)

    c0.assign(c1)

    print " time = %g" %t

    t += dt
    plot2 = plot(c1 , title = "c1"); time.sleep(5)
    if t == dt:
        plot2.write_png("assymetric.png")

Anyone have an idea to solve this?
Using OSX 10.8.5 and FEniCS V1.2
Thank you in advance.
A & C

Update: 23/1-14 at 12.50
We now found out that if we initialize the problem with
Ω_0 = [0 , 0.4]
Ω_1 = (0.4 , 0.6]
Ω_2 = (0.6 , 1]

And solve the problem with
Ω_0 = [0 , 0.4)
Ω_1 = [0.4 , 0.6]
Ω_2 = (0.6 , 1]

The solution becomes symmetric. But this cannot be the right way to go? There must be an error in FEniCS otherwise I don't understand this. Anyone who knows more about this?

mikael-mortensen · Answer 1 · Jan 23, 2014

Hi,

It's a lot of code to look through, but I think this could be because you're using interpolate instead of project when moving c0 from DG0 to CG1. Interpolate will use only one of the values from element left or right of the node, depending on which element is visited last. Change

c0 = interpolate(c0, V)

to

c0 = project(c0, V)

and I think you should be good.

Ok, I tested and it seems like project is symmetric, but leads to some initial oscillations. To get rid of those you could do as suggested by miroK. To get it truly symmetric for long runs, though, you need also use an odd number of elements.

MiroK · Answer 2 · Jan 23, 2014

Hi, apparently asymmetry is introduced into the system at some point. If you inspect your c0

print c0.vector().array()

you will see that it is in fact not symmetric; the left domain has more ones then the right one. Instead of solving a problem to get c0 I defined it as follows

c0 = interpolate(Expression("(x[0] >= 0.4 - DOLFIN_EPS) & (x[0] <= 0.6 + DOLFIN_EPS) ? 0.0 : 1"), V)

The rest of the code remained the same and the solution is symmetric all the way to steady state at T=1 or so.

Diffusion on symmetric problem becomes asymmetric.

2 Answers