Boundary Condition Issues in FEniCS with Discontinuous Galerkin Method

General
1

I am solving for ψ in a channel using the Discontinuous Galerkin (DG) method in FEniCS. While the solution satisfies boundary conditions for smaller channel heights H (around 10−6), it fails to do so at larger values (e.g., 50×27×10−6). I would appreciate any suggestions or insights into potential errors or incorrect implementations in this code. Below is a code example.

from fenics import *
from mshr import *
import numpy as np
import matplotlib.pyplot as plt

constants

H = 27*50e-6
epsilon_n = 7.0922e-10
zeta0 = 20e-3
Farad = 96485.33289
e = 1.602176634e-19
k = 1.380649e-23
T_flu = 298
c_i0 = 9.44e-5
R = 8.3145
degree = 2

parameter values

p1 = (H**2 * Farad * c_i0 * e) / (epsilon_n * k * T_flu)
a = 1
d = Constant((e * zeta0)/(k * T_flu))

mesh generation

channel_length, channel_height = 1, 1 / 27
patch_size = channel_height
channel = Rectangle(Point(0, 0), Point(channel_length, channel_height))
mesh = generate_mesh(channel, 500)

P = FunctionSpace(mesh, “DG”, degree)

Trial and test functions

psi_trial = TrialFunction(P)
psi = Function(P)
psi_test = TestFunction(P)

Penalty parameters

alpha = Constant(50.0)
alpha_2 = Constant(1000.0)

h = Constant(mesh.hmax())
n = FacetNormal(mesh)

Boundary Conditions

class InletBoundary(SubDomain):
def inside(self, x, on_boundary):
return abs(x[0] - 0) < DOLFIN_EPS and on_boundary

class OutletBoundary(SubDomain):
def inside(self, x, on_boundary):
return abs(x[0] - channel_length) < DOLFIN_EPS and on_boundary

inlet = InletBoundary()
outlet = OutletBoundary()

class DirichletBoundary_upper_left(SubDomain):
def inside(self, x, on_boundary):
return (
between(x[0], (- 3*patch_size/2 + channel_length/2,
-patch_size/2 + channel_length/2))
and abs(x[1] - channel_height) < DOLFIN_EPS
and on_boundary
)
class DirichletBoundary_upper_right(SubDomain):
def inside(self, x, on_boundary):
return (
between(x[0], (patch_size/2 + channel_length/2, 3 * patch_size/2 + channel_length/2))
and abs(x[1] - channel_height) < DOLFIN_EPS
and on_boundary
)

class DirichletBoundary_lower_left(SubDomain):
def inside(self, x, on_boundary):
return (
between(x[0], (- 3*patch_size/2 + channel_length/2,
-patch_size/2 + channel_length/2))
and abs(x[1] - 0) < DOLFIN_EPS
and on_boundary
)

class DirichletBoundary_lower_right(SubDomain):
def inside(self, x, on_boundary):
return (
between(x[0], (patch_size/2 + channel_length/2, 3 * patch_size/2 + channel_length/2))
and abs(x[1] - 0) < DOLFIN_EPS
and on_boundary
)

upper_left = DirichletBoundary_upper_left()
lower_left = DirichletBoundary_lower_left()
upper_right = DirichletBoundary_upper_right()
lower_right = DirichletBoundary_lower_right()

class DirichletBoundary_other(SubDomain):
def inside(self, x, on_boundary):
return (
on_boundary
and not (
upper_left.inside(x, on_boundary) or
lower_left.inside(x, on_boundary) or
upper_right.inside(x, on_boundary) or
lower_right.inside(x, on_boundary) or
inlet.inside(x, on_boundary) or
outlet.inside(x, on_boundary)
)
)

other = DirichletBoundary_other()

Notations for DG SIP

dx = Measure(‘dx’, domain=mesh)
ds = Measure(‘ds’, domain=mesh) # Exterior facets
dS = Measure(‘dS’, domain=mesh) # Interior facets

boundary_markers = MeshFunction(“size_t”, mesh, mesh.topology().dim() - 1)
boundary_markers.set_all(0)
upper_left.mark(boundary_markers, 1)
upper_right.mark(boundary_markers, 1)
lower_left.mark(boundary_markers, 1)
lower_right.mark(boundary_markers, 1)
other.mark(boundary_markers, 2)

ds = Measure(“ds”, domain=mesh, subdomain_data=boundary_markers)

u_patch = Constant(d)

Solve

a_volume = inner(grad(psi_trial), grad(psi_test)) * dx + p1 * a * psi_trial * psi_test * dx

a_interior = (
- inner(avg(grad(psi_trial)), jump(psi_test, n)) * dS
+ inner(avg(grad(psi_test)), jump(psi_trial, n)) * dS
+ alpha/h * inner(jump(psi_trial, n), jump(psi_test, n)) * dS
)

a_boundary = (
- inner(grad(psi_trial), psi_test * n) * ds(1) - inner((psi_trial) * n, grad(psi_test)) * ds(1) + alpha_2/h * (psi_trial) * psi_test * ds(1)
- inner(grad(psi_trial), psi_test * n) * ds(2)
- inner((psi_trial) * n, grad(psi_test)) * ds(2) + alpha_2/h * (psi_trial) * psi_test * ds(2)
)

billinear = a_volume + a_boundary + a_interior

L = p1 * psi_test *dx - u_patch * inner(n, grad(psi_test)) * ds(1) + alpha_2/h * (u_patch) * psi_test * ds(1)

solve(billinear == L, psi)

plt.figure(figsize=(20, 10))
c = plot(psi, cmap = “hot”)
plt.colorbar(c)
plt.title(‘Solution for ψ’)
plt.xlabel(‘x’)
plt.ylabel(‘y’)
plt.show()

2

Please format your code with markdown code syntax, ie

```python
#add code here
```

and ensure that the code is runnable by others.

Additionally, as your mesh is a square mesh, I would advice you to start by using a built in mesh, rather than mshr, as mshr has been deprecated for many years and is not a software many people now have installed on their system.

3

Thank you very much for your reply. I have used a built in mesh in this code. Please find the code below :

from fenics import *
import numpy as np
import matplotlib.pyplot as plt

# constants
H = 27*50e-6
epsilon_n = 7.0922e-10
zeta0 = 20e-3
Farad = 96485.33289
e = 1.602176634e-19
k = 1.380649e-23
T_flu = 298
c_i0 = 9.44e-5
R = 8.3145
degree = 2

# parameter values
p1 = (H**2 * Farad * c_i0 * e) / (epsilon_n * k * T_flu)
a = 1
d = Constant((e * zeta0)/(k * T_flu))

# mesh generation
channel_length, channel_height = 1, 1 / 27
patch_size = channel_height
mesh = RectangleMesh(Point(0, 0), Point(channel_length, channel_height), 100, 100) 

P = FunctionSpace(mesh, "DG", degree)

# Trial and test functions
psi_trial = TrialFunction(P)
psi = Function(P)
psi_test = TestFunction(P)

# Penalty parameters
alpha = Constant(50.0)
alpha_2 = Constant(1000.0)

h = Constant(mesh.hmax())
n = FacetNormal(mesh)

# Boundary Conditions 
class InletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[0] - 0) < DOLFIN_EPS and on_boundary

class OutletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[0] - channel_length) < DOLFIN_EPS and on_boundary

inlet = InletBoundary()
outlet = OutletBoundary()

class DirichletBoundary_upper_left(SubDomain):
    def inside(self, x, on_boundary):
        return (
            between(x[0], (- 3*patch_size/2 + channel_length/2,
                           -patch_size/2 + channel_length/2))
            and abs(x[1] - channel_height) < DOLFIN_EPS
            and on_boundary
        )
class DirichletBoundary_upper_right(SubDomain):
    def inside(self, x, on_boundary):
        return (
            between(x[0], (patch_size/2 + channel_length/2, 3 * patch_size/2 + channel_length/2))
            and abs(x[1] - channel_height) < DOLFIN_EPS
            and on_boundary
        )

class DirichletBoundary_lower_left(SubDomain):
    def inside(self, x, on_boundary):
        return (
            between(x[0], (- 3*patch_size/2 + channel_length/2,
                           -patch_size/2 + channel_length/2))
            and abs(x[1] - 0) < DOLFIN_EPS
            and on_boundary
        )

class DirichletBoundary_lower_right(SubDomain):
    def inside(self, x, on_boundary):
        return (
            between(x[0], (patch_size/2 + channel_length/2, 3 * patch_size/2 + channel_length/2))
            and abs(x[1] - 0) < DOLFIN_EPS
            and on_boundary
        )

upper_left = DirichletBoundary_upper_left()
lower_left = DirichletBoundary_lower_left()
upper_right = DirichletBoundary_upper_right()
lower_right = DirichletBoundary_lower_right()

class DirichletBoundary_other(SubDomain):
    def inside(self, x, on_boundary):
        return (
            on_boundary
            and not (
                upper_left.inside(x, on_boundary) or
                lower_left.inside(x, on_boundary) or
                upper_right.inside(x, on_boundary) or
                lower_right.inside(x, on_boundary) or
                inlet.inside(x, on_boundary) or
                outlet.inside(x, on_boundary)
            )
        )

other = DirichletBoundary_other()

# Notations for DG SIP 
dx = Measure('dx', domain=mesh)
ds = Measure('ds', domain=mesh)                                # Exterior facets
dS = Measure('dS', domain=mesh)                                # Interior facets

boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_markers.set_all(0)
upper_left.mark(boundary_markers, 1)
upper_right.mark(boundary_markers, 1)
lower_left.mark(boundary_markers, 1)
lower_right.mark(boundary_markers, 1)
other.mark(boundary_markers, 2)

ds = Measure("ds", domain=mesh, subdomain_data=boundary_markers)

u_patch = Constant(d)

# Solve 
a_volume = inner(grad(psi_trial), grad(psi_test)) * dx + p1 * a * psi_trial * psi_test * dx

a_interior = (
    - inner(avg(grad(psi_trial)), jump(psi_test, n)) * dS
    + inner(avg(grad(psi_test)), jump(psi_trial, n)) * dS
    + alpha/h * inner(jump(psi_trial, n), jump(psi_test, n)) * dS
)

a_boundary = (
    - inner(grad(psi_trial), psi_test * n) * ds(1) - inner((psi_trial) * n, grad(psi_test)) * ds(1) + alpha_2/h * (psi_trial) * psi_test * ds(1) \
    - inner(grad(psi_trial), psi_test * n) * ds(2) \
    - inner((psi_trial) * n, grad(psi_test)) * ds(2) + alpha_2/h * (psi_trial) * psi_test * ds(2)
)

billinear = a_volume + a_boundary + a_interior

L = p1 * psi_test *dx - u_patch * inner(n, grad(psi_test)) * ds(1) + alpha_2/h * (u_patch) * psi_test * ds(1)

solve(billinear == L, psi)

plt.figure(figsize=(20, 10))
c = plot(psi, cmap = "hot")
plt.colorbar(c)
plt.title('Solution for ψ')
plt.xlabel('x')
plt.ylabel('y')
plt.show()