Setting interior boundary condition for elliptic problem with DG

[plugged]

Hello everyone,

I have a question how to incorporate Dirichlet-type boundary conditions on interior surfaces for an elliptic problem using a discontinuous Galerkin (DG) discretization. Basically, I have to solve some problem on some holdall domain consisting of several sub-domains. I have an elliptic equation on some of the subdomains, and nothing on the others. On the exterior boundaries, I prescribe Dirichlet boundary conditions (and I know how to implement these), and on the interior boundaries I also want to prescribe Dirichlet boundary conditions. For the usual continuous Galerkin method, I know how to achieve this, with the help of the strong form DirichletBC. For the discontinuous Galerkin method, I am using a symmetric interior penalty approach which allows me to incorporate the exterior boundary conditions, but I am not sure how to treat the interior boundary conditions.

I have create the following MWE to illustrate this. The mesh is created using Gmsh and looks like this

The green subdomain is marked with 1 and the orange interior subdomain is marked with 2. Moreover, all exterior boundaries are marked with 1 and all interior boundaries with 2. I will attach the code for the Gmsh .geo file at the end of the post for reproducability.

Now, I can use the continuous Galerkin method as follows to achieve what I want

from fenics import *
import cashocs

mesh, subdomains, boundaries, dx, ds, dS = cashocs.import_mesh("mesh.xdmf")
u_out = Constant(1.0)
u_int = Constant(2.0)

V = FunctionSpace(mesh, "CG", 1)
u = TrialFunction(V)
v = TestFunction(V)

a = dot(grad(u), grad(v)) * dx(1)
L = Constant(0.0) * v * dx(1)

bc_out = DirichletBC(V, u_out, boundaries, 1)
bc_int = DirichletBC(V, u_int, boundaries, 2)
bcs = [bc_out, bc_int]

u_cg = Function(V)
A, b = assemble_system(a, L, bcs, keep_diagonal=True)
A.ident_zeros()
solve(A, u_cg.vector(), b)

and the result looks like this

Now, I don’t really care about the solution in the inner circle, in this case it is extended with 0, but I would also not mind if it were extended with u_int or any other (finite) value, it does not change the solution in subdomain 1.

Note, that I have used my package cashocs to import the mesh, but that is just for convenience as this creates the MeshFunctions and Measures which have the physical tags from the gmsh mesh, which I have converted to XDMF. Any other reader which can handle Gmsh meshes for FEniCS and incorporates the physical tags will also work.

My starting point for solving this problem with the discontinuous Galerkin method is the following, which only incorporates the boundary conditions on the exterior faces

from fenics import *
import cashocs

mesh, subdomains, boundaries, dx, ds, dS = cashocs.import_mesh("mesh.xdmf")
u_out = Constant(1.0)
u_int = Constant(2.0)

V = FunctionSpace(mesh, "DG", 1)
u = TrialFunction(V)
v = TestFunction(V)

sigma = 1e1
n = FacetNormal(mesh)
h = CellDiameter(mesh)
h_avg = (h("+") + h("-")) / 2.0

F = (
    dot(grad(u), grad(v)) * dx(1)
    - dot(grad(u), n) * v * ds
    - dot(grad(v), n) * (u - u_out) * ds
    + Constant(sigma) / h * (u - u_out) * v * ds
    - dot(avg(grad(u)), jump(v, n)) * dS
    - dot(avg(grad(v)), jump(u, n)) * dS
    + Constant(sigma) / h_avg * dot(jump(u, n), jump(v, n)) * dS
)
a, L = system(F)

u_dg = Function(V)
A, b = assemble_system(a, L, [], keep_diagonal=True)
A.ident_zeros()
solve(A, u_dg.vector(), b)

This code actually applies the exterior boundary conditions correctly. Of course, it does not do anything with the interior boundary at the moment, so that the result looks something like this

I have tried several ideas to also add the interior boundary conditions with this formulation, but had no success. Can anyone help me with this?

Thanks a lot in advance,
Sebastian

PS: For the sake of completeness, here is the code for the Gmsh .geo file

SetFactory("OpenCASCADE");
Rectangle(1) = {0, 0, 0, 1, 1, 0};
Disk(2) = {0.5, 0.5, 0, 0.25, 0.25};
BooleanDifference{ Surface{1}; Delete;}{ Surface{2};}

Mesh.MeshSizeMax = 5e-2;

Physical Surface(1) = {1};
Physical Surface(2) = {2};

Physical Curve(1) = {6,7,8,9};
Physical Curve(2) = {5};

[nate]

This is entirely possible, but a bit tricky even with UFL. You need to carefully manage the ("+") and ("-") sides and impose the exterior Nitsche terms accordingly. I have an implementation for the DG formulation of - \nabla \cdot \left(G(u, \nabla u) \nabla u \right) here: dolfin_dg/dolfin_dg/dg_form.py at master · nate-sime/dolfin_dg · GitHub. Unfortunately this work has become extremely messy as I’ve moved on to other things and not had time to tidy up the interface, but hopefully it gives you the gist of what you need to do.

[plugged]

Thanks a lot for your reply. I have just tested this and it seems to work. I will post the code containing the solution on monday for everyone who is interested in this.

Could you also, perhaps after I post my code, elaborate briefly why this approach works? In my understanding, you just treat the terms the same way they would if they arise on an actual outer boundary (your code uses the same code for this), but then you just sum up the contributions from the ("+") and ("-") sides - why does this work? Sorry for asking, but I am still new to DG methods - and in textbooks (at least the ones I am familiar with) the notation used is different to UFL.

[plugged]

As promised, here is my solution, based on the link you have sent, @nate.

The code is as follows:

from fenics import *
import cashocs

mesh, subdomains, boundaries, dx, ds, dS = cashocs.import_mesh("mesh.xdmf")
u_out = Constant(1.0)
u_int = Constant(2.0)

V = FunctionSpace(mesh, "DG", 1)
u = TrialFunction(V)
v = TestFunction(V)

sigma = 1e1
n = FacetNormal(mesh)
h = CellDiameter(mesh)
h_avg = (h("+") + h("-")) / 2.0

F = (
    dot(grad(u), grad(v)) * dx(1)
    - dot(grad(u), n) * v * ds
    - dot(grad(v), n) * (u - u_out) * ds
    + Constant(sigma) / h * (u - u_out) * v * ds
    - dot(avg(grad(u)), jump(v, n)) * dS
    - dot(avg(grad(v)), jump(u, n)) * dS
    + Constant(sigma) / h_avg * dot(jump(u, n), jump(v, n)) * dS
    + (
        -dot(grad(u), n) * v
        - dot(grad(v), n) * (u - u_int)
        + Constant(sigma) / h * (u - u_int) * v
    )("+")
    * dS(2)
    + (
        -dot(grad(u), n) * v
        - dot(grad(v), n) * (u - u_int)
        + Constant(sigma) / h * (u - u_int) * v
    )("-")
    * dS(2)
    + Constant(0.0) * v * dx(1)
)
a, L = system(F)

u_dg = Function(V)
A, b = assemble_system(a, L, [], keep_diagonal=True)
A.ident_zeros()
solve(A, u_dg.vector(), b)

and produces the following output

I have also thought about why this formulation is correct, but now it makes sense to me: Originally, I thought that I only had to prescribe the value on the interior side of the inner boundary, as this is what happens physically. However, doing so is not so simple as the + and - numbering of the cells is (somewhat) arbitrary as I understand it. So what is way easier is to prescribe the boundary condition from both sides, so that both the + and the - integrals are there, in the same fashion as they would occur for an outer boundary, but now they, of course, are present on both the interior as well as the exterior side of the inner boundary - so that makes sense.

Thanks a lot for you help @nate.

[nate]

@plugged

However, doing so is not so simple as the + and - numbering of the cells is (somewhat) arbitrary as I understand it.

I’m glad you have the result you desire. However I have a comment: This formulation is not arbitrary. By setting a Dirichlet boundary on the interior you are essentially splitting your problem into two domains \Omega_1 and \Omega_2 with shared interface \Gamma_{1,2}. Due to the discontinuity of the basis at that interface \Gamma_{1,2} one needs to impose appropriate boundary data in order for the problem to be well posed. You can also think of \Gamma_{1,2} being a component of the exterior boundary of both \Omega_1 and \Omega_2, independently. Therefore to impose Dirichlet data for the problems in both domains we impose the boundary data weakly on both sides (("+") and ("-")) of \Gamma_{1,2}.

Perhaps this could be elucidated further by considering that one could impose different boundary data for \Gamma_{1,2} for each side. This is permitted since your basis is discontinuous. It’s just a bit tricky to keep track of which side ("+") and ("-") correspond to in the physical problem.

[plugged]

Thanks for your answer - yes that is exactly how I understood this. And as I am (at the moment) not interested in using different boundary data on both sides of the interface, this is totally fine.