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())