Mixed boundary conditions silently violated

dolfinx
1

I am using DOLFINx v0.9.0 to solve the basic one-dimensional wave equation by space-time FEM (formulated as first-order system). My script runs without warnings or error messages, the results look qualitatively good, but quantitatively they are scaled by a factor that depends on the discretization. Suprisingly, the only non-zero BC (initial displacement) is scaled to by the same factor too. The exact solution for this model problem is just a product of sine and cosine with values between -1…1, but the numerical results are scaled by

  • 1e53, 1.7e4, 7e4 for nx=20 and nt=20, 50, 100, respectively;
  • 5e4, 2e5, 8e5 for nx=50 and nt=50, 100, 200, respectively.

What can be the reason for the BC not enforced correctly?

from mpi4py import MPI
import ufl
from dolfinx import fem, mesh, plot, default_scalar_type
from dolfinx.fem import dirichletbc, locate_dofs_geometrical
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import locate_entities_boundary
from basix.ufl import element, mixed_element
import numpy as np
#import matplotlib.pyplot as plt
#import pyvista as pv

# --- 1. Problem Parameters ---
c = 1.0  # Wave speed
nx = 50  # Number of elements in space
# dx = 1/50 = 0.02. dt should be smaller. T/nt < dx/c -> 2/nt < 0.02 -> nt > 100 (CFL condition)
nt = 200 # Number of elements in time

# --- 2. Create the Space-Time Mesh ---
# The domain is a 2D rectangle [0, L] x (0, T)
domain = mesh.create_rectangle(MPI.COMM_WORLD,
                               [np.array([0, 0]), np.array([1.0, 1.0])],
                               [nx, nt],
                               mesh.CellType.quadrilateral)

# x[0] is space, x[1] is time
x = ufl.SpatialCoordinate(domain)

# --- 3. Define the Function Space ---
P1 = element("Lagrange", domain.basix_cell(), 1)
W_element = mixed_element([P1, P1])
W = fem.functionspace(domain, W_element)

(u, v) = ufl.TrialFunctions(W)
(Du, Dv) = ufl.TestFunctions(W)

# Equation 1: du/dt - v = 0
a1 = (u.dx(1) * Dv - v * Dv) * ufl.dx
L1 = fem.Constant(domain, default_scalar_type(0.0)) * Dv * ufl.dx

# Equation 2: dv/dt - c^2 * d^2u/dx^2 = f
# Integrate by parts in space: -c^2 * d^2u/dx^2 -> c^2 * du/dx * dw_v/dx
f = fem.Constant(domain, default_scalar_type(0.0)) # Source term
a2 = (v.dx(1) * Du + c**2 * u.dx(0) * Du.dx(0)) * ufl.dx
L2 = f * Du * ufl.dx

# Combine into a single system
a = a1 + a2 
L_form = L1 + L2

# For x=0 (left boundary)
facets_left = locate_entities_boundary(domain, 1, lambda x: np.isclose(x[0], 0.0))
dofs_u_left = fem.locate_dofs_topological(W.sub(0), 1, facets_left)

# For x=1 (right boundary)
facets_right = locate_entities_boundary(domain, 1, lambda x: np.isclose(x[0], 1.0))
dofs_u_right = fem.locate_dofs_topological(W.sub(0), 1, facets_right)

# For t=0 (initial time)
facets_t0 = locate_entities_boundary(domain, 1, lambda x: np.isclose(x[1], 0.0))
dofs_u_initial = fem.locate_dofs_topological(W.sub(0), 1, facets_t0)
dofs_v_initial = fem.locate_dofs_topological(W.sub(1), 1, facets_t0)

# For t=1 (final time)
facets_t1 = locate_entities_boundary(domain, 1, lambda x: np.isclose(x[1], 1.0))
dofs_v_final = fem.locate_dofs_topological(W.sub(1), 1, facets_t1)

u_init = fem.Function(W.sub(0).collapse()[0])
u_init.interpolate(lambda x: np.sin(np.pi * x[0]))    # x[0]*(1-x[0])

# Create Dirichlet BCs for u and v at both boundaries
bc_u_left = dirichletbc(default_scalar_type(0), dofs_u_left, W.sub(0))
bc_u_right = dirichletbc(default_scalar_type(0), dofs_u_right, W.sub(0))
bc_u_initial = dirichletbc(u_init, dofs_u_initial)
bc_v_initial = dirichletbc(default_scalar_type(0), dofs_v_initial, W.sub(1))

# Collect all BCs in a list
bcs = [bc_u_left, bc_u_right, bc_u_initial, bc_v_initial]

# 1. Set up the linear problem
problem = LinearProblem(a, L_form, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})

# 2. Solve for the solution (mixed function: (u, v))
w_sol = problem.solve()

A = fem.petsc.assemble_matrix(fem.form(a), bcs)
A.assemble()
print("Matrix norm:", A.norm())

# Extract u and v components
u_sol = w_sol.sub(0).collapse()
v_sol = w_sol.sub(1).collapse()

# Print statistics for u
print("u min:", np.nanmin(u_sol.x.array))
print("u max:", np.nanmax(u_sol.x.array))
print("u mean:", np.nanmean(u_sol.x.array))
print("u contains NaN:", np.isnan(u_sol.x.array).any())

# Print statistics for v
print("v min:", np.nanmin(v_sol.x.array))
print("v max:", np.nanmax(v_sol.x.array))
print("v mean:", np.nanmean(v_sol.x.array))
print("v contains NaN:", np.isnan(v_sol.x.array).any())
2

The BCs were formally valid, but incorrectly applied. How to apply them correctly, is explained in detail in the tutorial Mixed finite element problems of the FEniCS Workshop by @dokken. For sake of completeness, I provide the corrected script:

from mpi4py import MPI
import dolfinx
from dolfinx.fem.petsc import LinearProblem
import basix.ufl
import ufl
import numpy as np
import matplotlib.pyplot as plt
import pyvista as pv

domain = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 4, 8, dolfinx.mesh.CellType.quadrilateral)

el_u = basix.ufl.element("Lagrange", domain.basix_cell(), 1)
el_v = basix.ufl.element("Lagrange", domain.basix_cell(), 1)
el_mixed = basix.ufl.mixed_element([el_u, el_v])

W = dolfinx.fem.functionspace(domain, el_mixed)
u, v = ufl.TrialFunctions(W)
Du, Dv = ufl.TestFunctions(W)

# Equation 1: du/dt - v = 0
a1 = (u.dx(1) * Dv - v * Dv) * ufl.dx
L1 = dolfinx.fem.Constant(domain, dolfinx.default_scalar_type(0.0)) * Dv * ufl.dx

# Equation 2: dv/dt - c^2 * d^2u/dx^2 = f
# Integrate by parts in space:   d^2u/dx^2   -->  du/dx * dw_v/dx
a2 = (v.dx(1) * Du + u.dx(0) * Du.dx(0)) * ufl.dx
f = dolfinx.fem.Constant(domain, dolfinx.default_scalar_type(0.0)) # Source term
L2 = f * Du * ufl.dx

# Combine into a single system
a = a1 + a2 
L_form = L1 + L2

def left_right_marker(x):
    return np.logical_or(np.isclose(x[0], 0.0), np.isclose(x[0], 1.0))
left_right_facets = dolfinx.mesh.locate_entities_boundary(domain, domain.topology.dim - 1, left_right_marker)

def initial_time_marker(x):
    return np.isclose(x[1], 0.0)
initial_time_facets = dolfinx.mesh.locate_entities_boundary(domain, domain.topology.dim - 1, initial_time_marker)

W0 = W.sub(0)
U, U_to_W0 = W0.collapse()

W1 = W.sub(1)
V, V_to_W1 = W1.collapse()

leftright_u_dofs   = dolfinx.fem.locate_dofs_topological((W0,U), domain.topology.dim - 1, left_right_facets)
initial_u_dofs = dolfinx.fem.locate_dofs_topological((W0,U), domain.topology.dim - 1, initial_time_facets)
initial_v_dofs = dolfinx.fem.locate_dofs_topological((W1,V), domain.topology.dim - 1, initial_time_facets)

leftright_u = dolfinx.fem.Function(U)
def leftright_f(x):
    values = np.zeros((1, x.shape[1]))
    values[0, :] = 0.0*x[0]
    return values
leftright_u.interpolate(leftright_f)
leftright_u_bc = dolfinx.fem.dirichletbc(leftright_u, leftright_u_dofs, W0)

initial_u = dolfinx.fem.Function(U)
def initial_u_f(x):
    values = np.zeros((1, x.shape[1]))
    values[0, :] = np.sin(np.pi * x[0])
    return values
initial_u.interpolate(initial_u_f)
initial_u_bc = dolfinx.fem.dirichletbc(initial_u, initial_u_dofs, W0)

initial_v = dolfinx.fem.Function(V)
def initial_v_f(x):
    values = np.zeros((1, x.shape[1]))
    values[0, :] = 0.0*x[0]
    return values
initial_v.interpolate(initial_v_f)
initial_v_bc = dolfinx.fem.dirichletbc(initial_v, initial_v_dofs, W1)

bcs = [leftright_u_bc, initial_u_bc, initial_v_bc]
problem = LinearProblem(a, L_form, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
w_sol = problem.solve()

u_plot = w_sol.sub(0).collapse()
#v_plot = w_sol.sub(1).collapse()

u_grid = pv.UnstructuredGrid(*dolfinx.plot.vtk_mesh(u_plot.function_space))
u_grid.point_data["u"] = u_plot.x.array
plotter_u = pv.Plotter()
plotter_u.add_mesh(u_grid, show_edges=True)
plotter_u.view_xy()
plotter_u.show()