DG FEM Flux approximation using TPFA scheme

variational formulation
1

Hi,
I am trying to implement DG-Fem for the Poisson equation. Here flux term is calculated using the Two-Point Flux Approximation scheme.(Reference Paper)

Find u_h \in V_h such that
\int_{\Gamma_{int}} [v]\frac{[u]}{|| h^+ - h^- ||} dS = \int_{\Gamma_N} v_h ds + \int_{\Omega} fv_h dx for all v_h \in V

Here h^+ and h^- is the center of the current cell and neighboring cell respectively.
Now I have implemented the same:

V = FunctionSpace(mesh,"DG",0)
bc = DirichletBC(V, Constant(0), "near(x[0],0.0) || near(x[0],1.0)")
f = Expression("10*exp(-(pow(x[0]-0.5,2) + pow(x[1]-0.5,2))/0.02)", degree=1)
g = Expression("sin(5*x[0])", degree=1)

x_ = interpolate(Expression("x[0]", degree=1),V)
y_ = interpolate(Expression("x[1]",degree=1),V)
Delta_h = sqrt(jump(x_)**2 + jump(y_)**2)

u = TrialFunction(V)
v = TestFunction(V)

a = jump(u)/Delta_h*jump(v)*dS
L = f*v*dx + g*v*ds
u = Function(V)
solve(a==L, u, bc)

plt_u = plot(u)
cbar = plt.colorbar(plt_u)

Output:
output2

The same is implemented using CFEM as follows:

from dolfin import *
import matplotlib.pyplot as plt
plt.rcParams['image.cmap'] = 'jet'

mesh = UnitSquareMesh(16,16,"left")
V = FunctionSpace(mesh, "CG", 1)
bc = DirichletBC(V, Constant(0), "near(x[0],0) || near(x[0],1)")

f = Expression("10*exp(-(pow(x[0]-0.5,2) + pow(x[1]-0.5,2))/0.02)", degree=1)
g = Expression("sin(5*x[0])", degree=1)

u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u), grad(v))*dx
L = f*v*dx + g*v*ds

u = Function(V)
solve(a==L, u, bc)
plt_u = plot(u)
cbar = plt.colorbar(plt_u)

Output:
output1

In the case of DG-Fem, output is wrong. Could you please tell me what wrong I am doing here?

Thank you,
Ranjeet

2

After having a quick look at your reference paper, I made some modifications to your code, and it should work:

from ufl import *
from dolfin import *
import matplotlib.pyplot as plt
plt.rcParams['image.cmap'] = 'jet'

mesh_param = 16

mesh = UnitSquareMesh(mesh_param,mesh_param,"crossed")
x = SpatialCoordinate(mesh)

class Dirichlet_left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0],0)
    
class Dirichlet_right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0],1)
    
class Neumann(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1],0) or near(x[1],1)
    
dirichlet_left = Dirichlet_left()
dirichlet_right = Dirichlet_right()
neumann = Neumann()
    
boundaries = MeshFunction("size_t", mesh, 1, 0)
dirichlet_left.mark(boundaries, 1)
dirichlet_right.mark(boundaries, 2)
neumann.mark(boundaries, 3)
ds = Measure("ds", subdomain_data = boundaries)

V = FunctionSpace(mesh,"DG",0)
f = 10 * exp( -( ( x[0]-0.5 )**2 + (x[1]-0.5)**2 ) / 0.02 )
g = sin(5*x[0])

x_ = interpolate(Expression("x[0]", degree=1),V)
y_ = interpolate(Expression("x[1]",degree=1),V)
Delta_h = sqrt(jump(x_)**2 + jump(y_)**2)

u = TrialFunction(V)
v = TestFunction(V)

a = jump(u)*jump(v)/Delta_h*dS + (u*v/x_) * ds(1) + (u*v/(1-x_)) * ds(2)
L = f*v*dx + g*v*ds(3) + (f*v/x_)*ds(1) + (f*v/(1-x_))*ds(2)
u = Function(V)
solve(a==L, u)

plt_u = plot(u)
cbar = plt.colorbar(plt_u)
plt.show()

DirichletBCs can not be applied to DG spaces, since the dofs are not on the boundary.

1 Like
3

Thank you for your reply. Now I got the correct solution.