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


96
views
0
3 months ago by
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.
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.

Similar posts:
Search »