Hello,
I want to solve the Poisson equation by defining its fields on DG spaces. The boundary value problem is
\nabla^2 u = f \, \text{in }\Omega
u = u_0 \text{ on } \partial\Omega
where \Omega is a disk.
I managed to solve this by imposing the boundary condition (BC) with a penalty term, but I would like to enforce it strongly as a DirichletBC, for reasons related to future developments.
Is this possible?
Here is the minimal working example where the BC is imposed with penalty (yes, I am purposedly using a nonlinear solver for future developments):
import numpy as np
from fenics import *
import ufl
# ── Mesh ─────────────────────────────────────────────────────────────────────
r = 1.0 # disk radius
N = 32 # mesh resolution (number of cells along diameter)
mesh = UnitDiscMesh.create(MPI.comm_world, N, 1, 2)
# ── Parameters ───────────────────────────────────────────────────────────────
degree = 2 # polynomial degree of DG space
alpha = 10.0 # penalty parameter (must be large enough)
h = CellDiameter(mesh)
n = FacetNormal(mesh)
# ── Function space ────────────────────────────────────────────────────────────
Q = FunctionSpace(mesh, 'DG', degree)
u = Function(Q)
nu = TestFunction(Q)
J_u = TrialFunction(Q)
# ── Exact solution and RHS ───────────────────────────────────────────────────
class u_exact_expression(UserExpression):
def eval(self, values, x):
values[0] = 1 + x[0]**2 + 2*x[1]**2
def value_shape(self):
return (1,)
class f_expression(UserExpression):
def eval(self, values, x):
values[0] = 6
def value_shape(self):
return (1,)
u_exact = Function(Q)
f = Function(Q)
u_exact.interpolate(u_exact_expression(element=Q.ufl_element()))
f.interpolate(f_expression(element=Q.ufl_element()))
i = ufl.indices(1)[0]
# ── Variational form ─────────────────────────────────────────────────────────
def jump(v, nn): return v('+') * nn('+') + v('-') * nn('-')
def avg(v): return 0.5 * (v('+') + v('-'))
# F_0: bulk term + natural Neumann term
F_0 = (u.dx(i)*nu.dx(i) + f*nu) * dx \
- n[i] * u.dx(i) * nu * ds
# F_I: interior penalty
F_I = (
- avg(u.dx(i)) * jump(nu, n)[i]
+ alpha/h('+') * jump(u, n)[i] * jump(nu, n)[i]
) * dS
F_E = (
alpha/h * (u - u_exact) * nu
) * ds
F = F_0 + F_I + F_E
solve(F == 0, u, [])
error_L2 = errornorm(u_exact, u, 'L2')
print(f"L2 error = {error_L2:.3e}")
and here the one where I try to impose it with DirichletBC
import numpy as np
from fenics import *
import ufl
# ── Mesh ─────────────────────────────────────────────────────────────────────
r = 1.0
N = 32
mesh = UnitDiscMesh.create(MPI.comm_world, N, 1, 2)
# ── Parameters ───────────────────────────────────────────────────────────────
degree = 2
h = CellDiameter(mesh)
n = FacetNormal(mesh)
# ── Function space ────────────────────────────────────────────────────────────
Q = FunctionSpace(mesh, 'DG', degree)
u = Function(Q)
nu = TestFunction(Q)
J_u = TrialFunction(Q)
# ── Exact solution and RHS ───────────────────────────────────────────────────
class u_exact_expression(UserExpression):
def eval(self, values, x):
values[0] = 1 + x[0]**2 + 2*x[1]**2
def value_shape(self):
return (1,)
class f_expression(UserExpression):
def eval(self, values, x):
values[0] = 6
def value_shape(self):
return (1,)
u_exact = Function(Q)
f = Function(Q)
u_exact.interpolate(u_exact_expression(element=Q.ufl_element()))
f.interpolate(f_expression(element=Q.ufl_element()))
alpha = 10.0
i = ufl.indices(1)[0]
# ── Strong Dirichlet BC ───────────────────────────────────────────────────────
class Boundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary
bc = DirichletBC(Q, u_exact, Boundary())
print(f"Constrained DOFs = {len(bc.get_boundary_values())}")
# ── Variational form ─────────────────────────────────────────────────────────
def jump(v, nn): return v('+') * nn('+') + v('-') * nn('-')
def avg(v): return 0.5 * (v('+') + v('-'))
F_0 = (u.dx(i)*nu.dx(i) + f*nu) * dx \
- n[i] * u.dx(i) * nu * ds
F_I = (
- avg(u.dx(i)) * jump(nu, n)[i]
+ alpha/h('+') * jump(u, n)[i] * jump(nu, n)[i]
) * dS
F = F_0 + F_I
solve(F == 0, u, [bc])
error_L2 = errornorm(u_exact, u, 'L2')
print(f"L2 error = {error_L2:.3e}")
The first one works, it gives:
# python3 mwe_penalty.py
No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 2.048e+02 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 1: r (abs) = 1.054e-12 (tol = 1.000e-10) r (rel) = 5.146e-15 (tol = 1.000e-09)
Newton solver finished in 1 iterations and 1 linear solver iterations.
*** Warning: Degree of exact solution may be inadequate for accurate result in errornorm.
L2 error = 1.395e-12
The second one fails (0 constrained DOFs), it gives:
# python3 mwe_dirichlet.py
Constrained DOFs = 0
No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 1.388e-01 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 1: r (abs) = 1.119e-12 (tol = 1.000e-10) r (rel) = 8.058e-12 (tol = 1.000e-09)
Newton solver finished in 1 iterations and 1 linear solver iterations.
*** Warning: Degree of exact solution may be inadequate for accurate result in errornorm.
L2 error = 4.462e+00
How can I make the Dirichlet code work by keeping DG spaces?
Thank you