Example using SLEPc for a complex eigenvalue problem¶
We consider a beam-like structure clamped at one end, assuming a linear elastic material. For illustration purpose, we try to solve the problem
with \(B=K\), where \(\mathrm{I}\) is the imaginary unit.
Note
You can download the complete python script
TestModalComplex.py
Imports and Environment Check
import numpy as np
import ufl
from dolfinx import fem, mesh, default_scalar_type
from dolfinx.fem.petsc import assemble_matrix
from mpi4py import MPI
from petsc4py import PETSc
from slepc4py import SLEPc
print("================START===================")
print(PETSc.ScalarType)
print("If not <class 'numpy.complex128'> run:")
print("source dolfinx-complex-mode")
print("in your Docker environment")
print("========================================")
#--------------------------------------------------------------------
Imports essential libraries for mesh creation, variational formulation, and eigenvalue solving.
Ensures complex arithmetic (numpy.complex128) is enmathrm{I}abled, as required for this problem.
Function to Create Complex Stiffness Matrix
# FUNCTION TO CREATE COMPLEX MATRIX A = K+j*B (j imag. unit)
#--------------------------------------------------------------------
def create_complex_stiffness(K,B):
"""Create a complex stiffness matrix A = K + j*B."""
# Ensure K and B are fully assembled
if not K.assembled:
K.assemble()
if not B.assembled:
B.assemble()
# Check if K and B are of same size
if not K.getSize()== B.getSize():
print("Error in create_complex_stiffness: K and B must be of same size")
# Create a new PETSc matrix for the complex stiffness
A = PETSc.Mat().createAIJ(size=K.getSize(), comm=MPI.COMM_WORLD)
A.setUp()
# Get the CSR representation of K
rows, cols = K.getSize()
ai, aj, av = K.getValuesCSR()
bi, bj, bv = K.getValuesCSR()
# Create a complex array for the new values
complex_values = av + 1j * bv
# Assemble A with the complex values
A.setValuesCSR(ai, aj, complex_values)
A.assemble()
return A
Avoids MatCopy Issue: Instead of copying K and modifying it, create_complex_stiffness directly constructs the matrix using setValuesCSR, ensuring compatibility with PETSc’s complex arithmetic.
Checks Matrix Dimensions: Ensures K and B are compatible.
Domain and Mesh Definition
#--------------------------------------------------------------------
# DOMAIN AND MESH
#--------------------------------------------------------------------
L = [5.0, 0.6, 0.4]
N = [25, 3, 2]
beam_mesh = mesh.create_box(MPI.COMM_WORLD, [np.zeros(3), np.array(L)],
N, mesh.CellType.hexahedron)
Defines a 3D hexahedral mesh for a beam-like domain.
Constitutive Law
#--------------------------------------------------------------------
# CONSTITUTIVE LAW (LINEAR ELASTICITY)
#--------------------------------------------------------------------
E, nu, rho = 2e11, 0.3, 7850
# Define complex-valued constants for linear elasticity
mu = fem.Constant(beam_mesh, PETSc.ScalarType(E / (2 * (1 + nu))))
lamda = fem.Constant(beam_mesh, PETSc.ScalarType(E * nu / ((1 + nu) * (1 - 2 * nu))))
def epsilon(u):
return ufl.sym(ufl.grad(u))
def sigma(u):
return lamda * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)
Defines the material properties and constitutive relations for linear elasticity, assuring compatibility with complex arithmetic using PETSc.ScalarType.
Function Space and Boundary Conditions
#--------------------------------------------------------------------
# FUNCTION SPACE, TRIAL AND TEST FUNCTION
#--------------------------------------------------------------------
V = fem.functionspace(beam_mesh, ("CG", 1,(beam_mesh.geometry.dim, )))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
#--------------------------------------------------------------------
# DIRICHLET BOUNDARY CONDITIONS
#--------------------------------------------------------------------
def clamped_boundary(x): return np.isclose(x[0], 0)
u_zero = np.array((0,) * beam_mesh.geometry.dim, dtype=default_scalar_type)
bc = fem.dirichletbc(u_zero, fem.locate_dofs_geometrical(V, clamped_boundary), V)
Defines the function space, trial and test functions, and clamped boundary conditions.
Assemble Stiffness and Mass Matricess and Complex Matrix A
#--------------------------------------------------------------------
# VARIATIONAL FORMS, STIFFNESS AND MASS MATRIX
#--------------------------------------------------------------------
k_form = fem.form(ufl.inner(sigma(u), epsilon(v)) * ufl.dx)
m_form = fem.form(rho * ufl.inner(u, v) * ufl.dx)
# Assemble matrices
K = assemble_matrix(k_form, bcs=[bc], diagonal=62831)
K.assemble()
M = assemble_matrix(m_form, bcs=[bc], diagonal=1./62831)
M.assemble()
#--------------------------------------------------------------------
# CREATE COMPLEX STIFFNESS MATRIX
#--------------------------------------------------------------------
A = create_complex_stiffness(K,K)
Assembles the stiffness and mass matrices using variational forms and creates a complex-valued stiffness matrix using the previously defined function.
Solve the Eigenvalue Problem and Writes the Results
#--------------------------------------------------------------------
# SLEPc EIGENSOLVER
#--------------------------------------------------------------------
#
# Set up SLEPc eigenvalue solver
eigensolver = SLEPc.EPS().create(MPI.COMM_WORLD)
eigensolver.setOperators(A, M)
eigensolver.setProblemType(SLEPc.EPS.ProblemType.GHEP)
eigensolver.setDimensions(16) # Number of eigenpairs to compute
eigensolver.setFromOptions()
# Solve the eigenvalue problem
eigensolver.solve()
#--------------------------------------------------------------------
# OUTPUT RESULTS
#--------------------------------------------------------------------
n_converged = eigensolver.getConverged()
if n_converged > 0:
vr, vi = A.getVecs()
for i in range(min(n_converged, 16)):
eigval = eigensolver.getEigenpair(i, vr, vi)
omega = np.sqrt(eigval.real)
print(f"Eigenvalue {i}: λ = {eigval:.6e}, ω = {omega:.6f} rad/s")
Remarks
Direct matrix operations (e.g., MatCopy) may fail if the matrix is not in a fully compatible state. This was resolved by using setValuesCSR to explicitly construct the complex matrix.
Enabling complex128 mode ensures compatibility for problems involving imaginary units.
Note
In case of an error message try first source dolfinx-complex-mode in your Docker environment.