Manufactured solution method for multi-material

1

Hi all,

I’m trying to verify a code using the method of manufactured solutions on a multi-material problem.

The equation I’m solving is the poisson equation \frac{du}{dt}=\nabla(D_i \nabla u) + \Gamma_i where D_i is the material property and \Gamma_i is the source term exact solution is u_D = 1 + \sin(2\pi x).

On a 1D problem divided in two subdomains (left and right) the source term should be :
\Gamma_i = - \frac{d}{dx}(D_i \frac{du_D}{dx})

I tried and implemented this in this MWE:

from fenics import *
import sympy as sp
import matplotlib.pyplot as plt

class Left(SubDomain):
    def inside(self, x, on_boundary):
        return between(x[0], (0, 0.5))


class Right(SubDomain):
    def inside(self, x, on_boundary):
        return between(x[0], (0.5, 1))


class Left_bound(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0)


class Right_bound(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 1)


left = Left()
right = Right()
left_bound = Left_bound()
right_bound = Right_bound()

mesh = UnitIntervalMesh(100)

volume_markers = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
volume_markers.set_all(0)
left.mark(volume_markers, 1)
right.mark(volume_markers, 2)

surface_markers = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
surface_markers.set_all(0)
left_bound.mark(surface_markers, 1)
right_bound.mark(surface_markers, 2)


V = FunctionSpace(mesh, 'P', 1)
u = Function(V)
v = TestFunction(V)

# Variational formulation
x = sp.Symbol('x[0]')

u_D = 1 + sp.sin(2*pi*x)  # exact solution

D_01, D_02 = 2, 5  # diffusion coeffs
# source term left
f = - sp.diff(D_01*sp.diff(u_D, x), x)
f = Expression(sp.printing.ccode(f), degree=2)

# source term right
g = - sp.diff(D_02*sp.diff(u_D, x), x)
g = Expression(sp.printing.ccode(g), degree=2)


dx = Measure('dx', domain=mesh, subdomain_data=volume_markers)

F = dot(D_01*grad(u), grad(v))*dx(1)
F += dot(D_02*grad(u), grad(v))*dx(2)

F += f*v*dx(1)
F += g*v*dx(2)

# BCs
bcs = [
    DirichletBC(V, 1, surface_markers, 1),
    DirichletBC(V, 1, surface_markers, 2)]

# Solving
solve(F == 0, u, bcs)
plot(u)
plt.savefig('out.png')

Unfortunately, this doesn’t produce the expected sin function:
image

Surprisingly enough, if D_01 and D_02 are equal, the correct function is obtained
image

Do you guys maybe see what I missed ? Or maybe my understanding of MMS is not great…

Cheers
Rem

2

In your code, you are enforcing:
D_1\frac{\partial u}{\partial x} = D_2\frac{\partial u}{\partial x} on the interface between the two materials.
When D_1\neq D_2, this will cause a discontinuity in the gradient across the interface.
To counteract this behavior, you need to do your integration by parts thoroughly, and will end up with a term like this:

F += -2*pi*cos(2*pi*0.5)*(D_01*v("-")-D_02*v("+"))*dS(3)

where dS(3) is the interface between the two domains.

I’ve attached my adaptation of your code (where I avoid using sympy and only use dolfin imports:

from fenics import *

import matplotlib.pyplot as plt

class Left(SubDomain):
    def inside(self, x, on_boundary):
        return between(x[0], (0, 0.5))


class Right(SubDomain):
    def inside(self, x, on_boundary):
        return between(x[0], (0.5, 1))


class Left_bound(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0)


class Right_bound(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 1)

class Interface(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.5)

left = Left()
right = Right()
left_bound = Left_bound()
right_bound = Right_bound()
interface = Interface()
mesh = UnitIntervalMesh(250)

volume_markers = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
volume_markers.set_all(0)
left.mark(volume_markers, 1)
right.mark(volume_markers, 2)

surface_markers = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
surface_markers.set_all(0)
left_bound.mark(surface_markers, 1)
right_bound.mark(surface_markers, 2)
interface.mark(surface_markers, 3)

V = FunctionSpace(mesh, 'P', 1)
u = TrialFunction(V)
v = TestFunction(V)

# Variational formulation

x= SpatialCoordinate(mesh)
u_D = 1 + sin(2*pi*x[0])  # exact solution

D_01, D_02 = Constant(2), Constant(5)  # diffusion coeffs
# source term left
f = grad(D_01*grad(u_D))[0,0]

# source term right
g = grad(D_02*grad(u_D))[0,0]


dx = Measure('dx', domain=mesh, subdomain_data=volume_markers)
dS = Measure("dS", domain=mesh, subdomain_data=surface_markers)

F = dot(D_01*grad(u), grad(v))*dx(1)
F += dot(D_02*grad(u), grad(v))*dx(2)

F += f*v*dx(1)
F += g*v*dx(2)
F += -2*pi*cos(2*pi*0.5)*(D_01*v("-")-D_02*v("+"))*dS(3)
# BCs
File("sf.pvd") << surface_markers
bcs = [
    DirichletBC(V, 1, surface_markers, 1),
    DirichletBC(V, 1, surface_markers, 2)]
# Solving
uh = Function(V)
solve(lhs(F)==rhs(F), uh, bcs=bcs)
plot(uh)
plt.savefig('out.png')
print(assemble(inner(uh-u_D,uh-u_D)*dx))
3

Hi dokken, thanks for your answer.

Is this needed every time a multi-material simulation is made ? Cause I haven’t seen anything alike in this demo.

https://fenicsproject.org/pub/tutorial/html/._ftut1013.html

Cheers

4

It depends on what interface condition you would like to to prescribe to your problem. By default, you are prescribing

when solving a heat equation, due to the integration by parts.

5

Say we want to prescribe this condition on the interface, why isn’t this additional term in the demo then ?

6

As it is a natural boundary condition, like a homogeneous Neumann condition. By excluding it from the variational form, you are enforcing it variationally.

7

Ok I think I get it. Thanks for the quick recap!