### How to capture a jump in the solution along an internal boundary using discontinuous galerkin methods

96

views

0

I am trying to solve a poisson equation in a unit square domain $\Omega$Ω with homogeneous Dirichlet boundary conditions on $\partial\Omega$∂Ω, along with an internal boundary at x=0.5, $\Gamma$Γ, where I enforce a jump in the solution (as in https://fenicsproject.org/qa/10927/neumann-bc-depending-on-solution/). The problem is therefore written as

\begin{align}

-\nabla^2 u &= f, \,\text{in }\Omega,\\

u &= 0, \,\text{on }\partial\Omega,\\

\nabla u \cdot \mathbf{n} &= k[\![u]\!],\,\text{on }\Gamma.

\end{align}

I want to solve this problem using discontinuous galerkin methods. Using the comments in https://fenicsproject.org/qa/10927/neumann-bc-depending-on-solution/, along with the tutorial given in https://github.com/FEniCS/dolfin/blob/master/demo/undocumented/dg-poisson/python/demo_dg-poisson.py, I have the following code. I add in a term in the variational form that I believe should be the necessary addition for the facets along the internal boundary.

\begin{align}

-\nabla^2 u &= f, \,\text{in }\Omega,\\

u &= 0, \,\text{on }\partial\Omega,\\

\nabla u \cdot \mathbf{n} &= k[\![u]\!],\,\text{on }\Gamma.

\end{align}

I want to solve this problem using discontinuous galerkin methods. Using the comments in https://fenicsproject.org/qa/10927/neumann-bc-depending-on-solution/, along with the tutorial given in https://github.com/FEniCS/dolfin/blob/master/demo/undocumented/dg-poisson/python/demo_dg-poisson.py, I have the following code. I add in a term in the variational form that I believe should be the necessary addition for the facets along the internal boundary.

```
from dolfin import *
class DirichletBoundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary
class RobinInterface(SubDomain):
def inside(self, x, on_boundary):
return near(0.5 - x[0], 0)
N = 32
mesh = UnitSquareMesh(N,N)
V = FunctionSpace(mesh, "DG", 1)
u = TrialFunction(V)
v = TestFunction(V)
n = FacetNormal(mesh)
h = CellDiameter(mesh)
h_avg = (h('+') + h('-'))/2
x = SpatialCoordinate(mesh)
I0 = conditional(x[0] < 0.5, 1, 0)
I1 = conditional(x[0] > 0.5, 1, 0)
k = Constant(8.0)
U0 = 4*x[0]**2
U1 = (1 + 4 * (1/k)) * 4 * x[0] * (1 - x[0]) \
+ 4 * (2*x[0] - 1) * (1 - x[0])
U = (I0 * U0 + I1 * U1) * sin(pi * x[1])
f = - div(grad(U))
u0 = Constant(0.0)
g = Constant(0.0)
boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1, 0)
boundaries.set_all(0)
DirichletBoundary().mark(boundaries, 1)
RobinInterface().mark(boundaries, 4)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
dS = Measure('dS', domain=mesh, subdomain_data=boundaries)
alpha = 4.0
gamma = 8.0
a = dot(grad(v), grad(u))*dx \
- dot(avg(grad(u)), jump(v, n))*dS(0) \
- dot(avg(grad(v)), jump(u, n))*dS(0) \
+ alpha/h_avg * dot(jump(u, n), jump(v, n))*dS(0) \
- dot(grad(u), v*n)*ds(1) \
- dot(grad(v), u*n)*ds(1) \
+ gamma/h *u*v*ds(1) \
- dot(k*jump(u, n), jump(v, n))*dS(4) \
L = v*f*dx - u0*dot(grad(v), n)*ds(1) + (gamma/h)*u0*v*ds(1)
u = Function(V)
solve(a == L, u)
u_exact = project(U, V)
diff = u-u_exact
error = assemble((inner(diff,diff)*dx))**0.5/assemble(inner(u,u)*dx)**0.5
```

However, I do not get error convergence as I refine the mesh, as I do not think that this is correctly capturing the jump at the internal interface $\Gamma$Γ. I would be so grateful if someone could see where I've gone wrong here. Thank you!
Community: FEniCS Project

Please login to add an answer/comment or follow this question.