Hello everyone, I am a beginner in coding, and I don’t know why all my modes are 1. I really appreciate your help （fenicsx0.10x）
Here is the code：
import numpy as np
from mpi4py import MPI
import dolfinx
from dolfinx import mesh, fem
import dolfinx.fem.petsc
import ufl
from petsc4py import PETSc
from scipy.linalg import eigh
# ----------------------------------------------------------------------
# 1. Mesh and function space
# ----------------------------------------------------------------------
comm = MPI.COMM_WORLD
rank = comm.Get_rank()
domain = mesh.create_unit_square(comm, 16, 16) # coarse mesh for speed
V = fem.functionspace(domain, ("Lagrange", 1)) # scalar linear elements
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
# ----------------------------------------------------------------------
# 2. Variational forms: stiffness (Laplacian) and mass (L2)
# ----------------------------------------------------------------------
a_K = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
a_M = ufl.inner(u, v) * ufl.dx
# ----------------------------------------------------------------------
# 3. Dirichlet boundary condition on left edge (x = 0)
# ----------------------------------------------------------------------
def left_boundary(x):
return np.isclose(x[0], 0.0)
fdim = domain.topology.dim - 1
left_facets = mesh.locate_entities_boundary(domain, fdim, left_boundary)
dofs_left = fem.locate_dofs_topological(V, fdim, left_facets)
bc = fem.dirichletbc(PETSc.ScalarType(0.0), dofs_left, V)
# ----------------------------------------------------------------------
# 4. Assemble matrices K and M with consistent boundary conditions
# ----------------------------------------------------------------------
K = dolfinx.fem.petsc.assemble_matrix(fem.form(a_K), bcs=[bc])
K.assemble()
M = dolfinx.fem.petsc.assemble_matrix(fem.form(a_M), bcs=[])
M.assemble()
M.zeroRowsColumns(dofs_left, diag=1.0) # apply Dirichlet to mass matrix
M.assemble()
# Quick check of matrix magnitudes (only rank 0)
if rank == 0:
k_diag = K.getDiagonal().array
m_diag = M.getDiagonal().array
print(f"K diag: min={k_diag.min():.2e}, max={k_diag.max():.2e}, mean={k_diag.mean():.2e}")
print(f"M diag: min={m_diag.min():.2e}, max={m_diag.max():.2e}, mean={m_diag.mean():.2e}")
# ----------------------------------------------------------------------
# 5. Linear solver for K y = b (direct LU with MUMPS)
# ----------------------------------------------------------------------
ksp = PETSc.KSP().create(comm)
ksp.setOperators(K)
ksp.setType(PETSc.KSP.Type.PREONLY)
pc = ksp.getPC()
pc.setType(PETSc.PC.Type.LU)
pc.setFactorSolverType("mumps")
ksp.setFromOptions()
# ----------------------------------------------------------------------
# 6. Subspace iteration parameters
# ----------------------------------------------------------------------
p = 6 # number of eigenpairs sought
q = p # subspace dimension (can be increased, e.g. q=2p)
max_iter = 30
tol = 1e-8
# Initial random subspace X (list of PETSc Vecs)
X_vecs = []
for i in range(q):
vec = K.createVecRight()
while True:
rctx = PETSc.Random().create(comm=comm)
rctx.setType(PETSc.Random.Type.RAND)
rctx.setSeed(i + rank * 100 + 12345)
vec.setRandom(rctx)
if vec.dot(M * vec) > 1e-12:
break
X_vecs.append(vec)
# ----------------------------------------------------------------------
# Helper functions: M‑normalization and M‑orthonormalization
# ----------------------------------------------------------------------
def m_normalize(vec, M):
norm2 = vec.dot(M * vec)
if norm2 < 1e-12:
vec.axpy(1e-8, vec) # tiny perturbation
norm2 = vec.dot(M * vec)
norm = np.sqrt(norm2)
vec.scale(1.0 / norm)
return norm
def m_orthonormalize(vecs, M):
"""Modified Gram–Schmidt with respect to the M inner product."""
n = len(vecs)
for i in range(n):
for j in range(i):
prod = vecs[j].dot(M * vecs[i])
vecs[i].axpy(-prod, vecs[j])
m_normalize(vecs[i], M)
# Initial orthonormalization
m_orthonormalize(X_vecs, M)
# ----------------------------------------------------------------------
# 7. Subspace iteration loop
# ----------------------------------------------------------------------
old_ritz = np.full(p, np.inf)
for it in range(max_iter):
# --- (1) Simultaneous inverse iteration: K Y = M X ---
Y_vecs = []
for i in range(q):
y = K.createVecRight()
Mx = M * X_vecs[i]
ksp.solve(Mx, y)
Y_vecs.append(y)
# --- (2) Rayleigh–Ritz projection ---
# Compute K_r = Yᵀ K Y and M_r = Yᵀ M Y (dense q×q matrices)
K_r_local = np.zeros((q, q))
M_r_local = np.zeros((q, q))
KY_vecs = [K * y for y in Y_vecs] # precompute K*Y
for i in range(q):
for j in range(q):
K_r_local[i, j] = Y_vecs[i].dot(KY_vecs[j])
M_r_local[i, j] = Y_vecs[i].dot(M * Y_vecs[j])
# Global reduction (MPI all‑reduce)
K_r_global = np.zeros_like(K_r_local)
M_r_global = np.zeros_like(M_r_local)
comm.Allreduce(K_r_local, K_r_global, op=MPI.SUM)
comm.Allreduce(M_r_local, M_r_global, op=MPI.SUM)
# Solve the reduced eigenproblem on rank 0
if rank == 0:
K_r_sym = 0.5 * (K_r_global + K_r_global.T)
M_r_sym = 0.5 * (M_r_global + M_r_global.T)
ritz_vals, Q = eigh(K_r_sym, M_r_sym) # ascending eigenvalues
else:
Q = None
ritz_vals = None
ritz_vals = comm.bcast(ritz_vals, root=0)
Q = comm.bcast(Q, root=0)
# --- (3) Update subspace: X_new = Y * Q ---
X_new_vecs = []
for i in range(q):
new_vec = K.createVecRight()
new_vec.set(0.0)
for j in range(q):
coeff = Q[j, i]
if abs(coeff) > 1e-14:
new_vec.axpy(coeff, Y_vecs[j])
X_new_vecs.append(new_vec)
m_orthonormalize(X_new_vecs, M) # numerical stabilization
# --- (4) Convergence test ---
current_ritz = ritz_vals[:p]
if it > 0:
err = np.linalg.norm(current_ritz - old_ritz) / (np.linalg.norm(current_ritz) + 1e-12)
if rank == 0:
print(f"Iter {it+1:2d}: error = {err:.2e}")
if err < tol:
if rank ==`Preformatted text` 0:
print("Subspace iteration converged!")
break
old_ritz = current_ritz.copy()
# Replace old subspace with new one
for vec in X_vecs:
vec.destroy()
X_vecs = X_new_vecs
# ----------------------------------------------------------------------
# 8. Output eigenvalues and frequencies
# ----------------------------------------------------------------------
if rank == 0:
print("\nFinal eigenvalues (λ):")
for i, val in enumerate(current_ritz):
print(f"Mode {i+1}: λ = {val:.8e}")
print("\nFrequencies (ω = sqrt(λ)):")
for i, w in enumerate(np.sqrt(current_ritz)):
print(f"Mode {i+1}: ω = {w:.6f}")
# Cleanup
ksp.destroy()
K.destroy()
M.destroy()
for vec in X_vecs:
vec.destroy()
type or paste code here