I want to implement code that uses the subspace iteration method to obtain modes and then perform topology optimization using MMA. But it is too slow and too complicated, and I hope someone can give me some advice.
Here is the code
import shutil
import os
cache_dir = os.path.expanduser(“~/.cache/fenics”)
if os.path.exists(cache_dir):
shutil.rmtree(cache_dir)
print(“Cleared FEniCS cache”)
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
import time
comm = MPI.COMM_WORLD
rank = comm.Get_rank()
lx, ly = 2.0, 1.0
nx, ny = 80, 40
domain = mesh.create_rectangle(comm, [[0.0, 0.0], [lx, ly]], [nx, ny],
cell_type=mesh.CellType.quadrilateral)
V = fem.functionspace(domain, (“Lagrange”, 1))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a_K0 = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx # unit stiffness matrix
a_M0 = ufl.inner(u, v) * ufl.dx # unit mass matrix
Boundary condition (fixed left end)
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)
Penalty parameter for boundary DOFs
BIG_PENALTY = 1e12
================================================================
2. Design variable space (DG-0, constant per cell)
================================================================
V_design = fem.functionspace(domain, (“DG”, 0))
rho = fem.Function(V_design, name=“Density”)
vol_frac = 0.3
rho.interpolate(lambda x: np.full(x[0].shape, vol_frac))
SIMP interpolation parameters
p_simp = 3.0 # stiffness interpolation exponent
q_simp = 3.0 # mass interpolation exponent (use >p_simp to suppress localized modes)
E0 = 1.0
Emin = 1e-9
def E_simp(rho_e):
return Emin + rho_e**p_simp * (E0 - Emin)
def M_simp(rho_e):
return rho_e**q_simp # density interpolation factor
================================================================
3. Assemble stiffness matrix K and mass matrix M (with penalty BC)
================================================================
def assemble_matrices(rho):
# Interpolated Young’s modulus (cell-wise, evaluated at quadrature points)
E_eff = Emin + rhop_simp * (E0 - Emin)
# Reconstruct variational forms
a_K_rho = E_eff * ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
a_M_rho = rhoq_simp * ufl.inner(u, v) * ufl.dx
K_mat = dolfinx.fem.petsc.assemble_matrix(fem.form(a_K_rho), bcs=[])
K_mat.assemble()
M_mat = dolfinx.fem.petsc.assemble_matrix(fem.form(a_M_rho), bcs=[])
M_mat.assemble()
# Apply penalty BC (push boundary DOF eigenvalues to high frequency)
K_mat.zeroRowsColumns(dofs_left, diag=BIG_PENALTY)
M_mat.zeroRowsColumns(dofs_left, diag=1.0)
K_mat.assemble()
M_mat.assemble()
return K_mat, M_mat
================================================================
4. Subspace iteration solver (fully encapsulated from user code)
================================================================
def subspace_iteration(K_mat, M_mat, p, q, max_iter=50, tol=1e-8):
mpi_comm = MPI.COMM_WORLD
rank = mpi_comm.Get_rank()
# ----- 4.1 Create KSP solver (for K y = M x) -----
# Create KSP solver using PETSc communicator
petsc_comm = K_mat.comm
ksp = PETSc.KSP().create(petsc_comm)
ksp.setOperators(K_mat)
ksp.setType(PETSc.KSP.Type.PREONLY)
pc = ksp.getPC()
pc.setType(PETSc.PC.Type.LU)
pc.setFactorSolverType("mumps")
ksp.setFromOptions()
# ----- 4.2 Initialize subspace basis (random vectors, M-orthonormalized) -----
X_vecs = []
for i in range(q):
vec = K_mat.createVecRight()
while True:
rctx = PETSc.Random().create(comm=comm)
rctx.setType(PETSc.Random.Type.RAND)
rctx.setSeed(i + rank * 100 + 12345)
vec.setRandom(rctx)
norm2 = vec.dot(M_mat * vec)
if norm2 > 1e-12:
break
X_vecs.append(vec)
# M-orthonormalization functions
def m_normalize(vec, M):
norm2 = vec.dot(M * vec)
if norm2 < 1e-12:
vec.axpy(1e-8, vec)
norm2 = vec.dot(M * vec)
norm = np.sqrt(norm2)
vec.scale(1.0 / norm)
return norm
def m_orthonormalize(vecs, M):
qlen = len(vecs)
for i in range(qlen):
for j in range(i):
prod = vecs[j].dot(M * vecs[i])
vecs[i].axpy(-prod, vecs[j])
m_normalize(vecs[i], M)
m_orthonormalize(X_vecs, M_mat)
old_ritz = np.full(p, np.inf)
for it in range(max_iter):
# Simultaneous inverse iteration: solve K Y = M X
Y_vecs = []
for i in range(q):
y = K_mat.createVecRight()
Mx = M_mat * X_vecs[i]
ksp.solve(Mx, y)
Y_vecs.append(y)
# Rayleigh-Ritz projection
K_r_local = np.zeros((q, q))
M_r_local = np.zeros((q, q))
KY_vecs = [K_mat * y for y in Y_vecs]
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_mat * Y_vecs[j])
# MPI global sum
K_r_global = np.zeros_like(K_r_local)
M_r_global = np.zeros_like(M_r_local)
mpi_comm.Allreduce(K_r_local, K_r_global, op=MPI.SUM)
mpi_comm.Allreduce(M_r_local, M_r_global, op=MPI.SUM)
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)
eps = 1e-10
success = False
for _ in range(10):
try:
M_reg = M_r_sym + eps * np.eye(q)
ritz_vals, Q = eigh(K_r_sym, M_reg)
# Check positivity
if np.min(ritz_vals) > -1e-12:
success = True
break
except np.linalg.LinAlgError:
eps *= 10
if not success:
print("Warning: generalized eigenproblem ill-conditioned, fallback to standard eig.")
ritz_vals, Q = eigh(K_r_sym)
# M-orthonormalize eigenvectors
for i in range(q):
norm_m = np.sqrt(Q[:, i] @ M_r_sym @ Q[:, i])
if norm_m > 1e-12:
Q[:, i] /= norm_m
else:
ritz_vals = None
Q = None
# Broadcast (must be outside if-else)
ritz_vals = mpi_comm.bcast(ritz_vals, root=0)
Q = mpi_comm.bcast(Q, root=0)
# Update subspace basis
X_new_vecs = []
for i in range(q):
new_vec = K_mat.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_mat)
# Convergence check
current_ritz = ritz_vals[:p]
if it > 0:
err = np.linalg.norm(current_ritz - old_ritz) / (np.linalg.norm(current_ritz) + 1e-12)
if err < tol:
if rank == 0:
print(f" Subspace iteration converged at step {it+1}")
break
old_ritz = current_ritz.copy()
# Clean old vectors
for vec in X_vecs:
vec.destroy()
X_vecs = X_new_vecs
# ----- 4.4 Extract first p eigenvectors -----
eigenvalues = current_ritz[:p]
eigenvectors = X_vecs[:p]
# Clean KSP and remaining vectors
ksp.destroy()
for vec in X_vecs[p:]:
vec.destroy()
return eigenvalues, eigenvectors
def compute_sensitivity(rho, phi, eigenvalue):
phi_func = fem.Function(V)
phi_func.x.array[:] = phi.array
# Create DG-0 space (constant per cell)
V0 = fem.functionspace(domain, ("DG", 0))
v0 = ufl.TestFunction(V0)
# Element stiffness energy ∫_e ∇φ·∇φ dx, multiplied by v0 to obtain cell-wise integrals
energy_K0_form = ufl.inner(ufl.grad(phi_func), ufl.grad(phi_func)) * v0 * ufl.dx
cell_energy_K0_vec = fem.assemble_vector(fem.form(energy_K0_form))
cell_energy_K0 = cell_energy_K0_vec.array.copy()
# Element mass energy ∫_e φ·φ dx
energy_M0_form = ufl.inner(phi_func, phi_func) * v0 * ufl.dx
cell_energy_M0_vec = fem.assemble_vector(fem.form(energy_M0_form))
cell_energy_M0 = cell_energy_M0_vec.array.copy()
rho_e = rho.x.array
dK_drho = p_simp * (rho_e**(p_simp-1)) * (E0 - Emin) * cell_energy_K0
dM_drho = q_simp * (rho_e**(q_simp-1)) * cell_energy_M0
sens = dK_drho - eigenvalue * dM_drho
return sens
def sensitivity_filter(sens, rho, rmin):
mesh = domain
num_cells = len(sens)
# Obtain cell centers
geometry = mesh.geometry
cells = mesh.topology.connectivity(mesh.topology.dim, 0)
cell_centers = np.zeros((num_cells, 2))
for i in range(num_cells):
verts = cells.links(i)
center = np.mean(geometry.x[verts], axis=0) # shape (3,)
cell_centers[i] = center[:2]
filtered = np.zeros_like(sens)
for i in range(num_cells):
w_sum = 0.0
val = 0.0
xi = cell_centers[i]
for j in range(num_cells):
xj = cell_centers[j]
dist = np.linalg.norm(xi - xj)
if dist < rmin:
w = rmin - dist
w_sum += w
val += w * sens[j]
filtered[i] = val / w_sum if w_sum > 0 else sens[i]
return filtered
Hand-written MMA solver (single constraint, dual method)
class MMA_SingleConstraint:
def init(self, x0, move=0.2):
self.n = len(x0)
self.x = x0.copy()
self.move = move
# Initial asymptotes
self.L = np.maximum(x0 - move, 0.0)
self.U = np.minimum(x0 + move, 1.0)
self.alpha = np.zeros(self.n)
self.beta = np.ones(self.n)
self.x_prev = x0.copy()
self.x_prev2 = x0.copy()
def update_asymptotes(self, x_curr, iter_num):
for j in range(self.n):
if iter_num <= 2:
Lj = x_curr[j] - self.move
Uj = x_curr[j] + self.move
else:
s = (x_curr[j] - self.x_prev[j]) * (self.x_prev[j] - self.x_prev2[j])
if s < 0:
Lj = x_curr[j] - 0.7 * (self.U[j] - self.L[j])
Uj = x_curr[j] + 0.7 * (self.U[j] - self.L[j])
else:
Lj = x_curr[j] - 1.2 * (self.U[j] - self.L[j])
Uj = x_curr[j] + 1.2 * (self.U[j] - self.L[j])
Lj = max(Lj, 0.0)
Uj = min(Uj, 1.0)
if Uj - Lj < 1e-6:
Lj = x_curr[j] - 0.5
Uj = x_curr[j] + 0.5
Lj = max(Lj, 0.0)
Uj = min(Uj, 1.0)
self.L[j] = Lj
self.U[j] = Uj
self.x_prev2 = self.x_prev.copy()
self.x_prev = x_curr.copy()
def compute_pq(self, grad_f0, grad_f1):
p0 = np.zeros(self.n)
q0 = np.zeros(self.n)
p1 = np.zeros(self.n)
q1 = np.zeros(self.n)
for j in range(self.n):
xj = self.x[j]
Lj = self.L[j]
Uj = self.U[j]
df0 = grad_f0[j]
p0[j] = max(0.0, (Uj - xj)**2 * df0) + 1e-12
q0[j] = max(0.0, -(xj - Lj)**2 * df0) + 1e-12
df1 = grad_f1[j]
p1[j] = max(0.0, (Uj - xj)**2 * df1) + 1e-12
q1[j] = max(0.0, -(xj - Lj)**2 * df1) + 1e-12
return p0, q0, p1, q1
def solve_subproblem(self, p0, q0, p1, q1):
# Dual variable μ Newton iteration
mu = 0.0
for _ in range(20):
x_new, f1_approx = self._compute_x_and_f1(mu, p0, q0, p1, q1)
if f1_approx <= 0:
break
eps = 1e-6
_, f1_plus = self._compute_x_and_f1(mu + eps, p0, q0, p1, q1)
df1 = (f1_plus - f1_approx) / eps
if df1 < 1e-12:
break
mu = mu - f1_approx / df1
mu = max(mu, 0.0)
x_new, _ = self._compute_x_and_f1(mu, p0, q0, p1, q1)
return x_new
def _compute_x_and_f1(self, mu, p0, q0, p1, q1):
x_new = np.zeros(self.n)
f1_val = 0.0
for j in range(self.n):
Pj = p0[j] + mu * p1[j]
Qj = q0[j] + mu * q1[j]
Lj = self.L[j]
Uj = self.U[j]
sqrtP = np.sqrt(Pj)
sqrtQ = np.sqrt(Qj)
x_star = (sqrtP * Lj + sqrtQ * Uj) / (sqrtP + sqrtQ)
xj = np.clip(x_star, self.alpha[j], self.beta[j])
x_new[j] = xj
f1_val += p1[j] / (Uj - xj) + q1[j] / (xj - Lj)
return x_new, f1_val
================================================================
7. Main optimization loop
================================================================
max_iters = 80
tol_change = 1e-4
filter_radius = 0.04
move_limit = 0.2
Volume constraint
total_volume = fem.assemble_scalar(fem.form(1.0 * ufl.Measure(“dx”, domain=domain)))
Vmax = vol_frac * total_volume
Initialize design variables
rho_vec = np.full(domain.topology.index_map(domain.topology.dim).size_local, vol_frac)
rho.x.array[:] = rho_vec
MMA object
mma = MMA_SingleConstraint(rho_vec, move=move_limit)
Record history
history_lambda =
for iter_opt in range(max_iters):
# ----- 7.1 Assemble K, M for current density -----
K_mat, M_mat = assemble_matrices(rho)
# ----- 7.2 Subspace iteration for fundamental frequency and mode shape -----
# Note: only p=1 eigenvalue needed, subspace dimension q=6~8
evals, evecs = subspace_iteration(K_mat, M_mat, p=1, q=6, max_iter=20, tol=1e-8)
lambda1 = evals[0]
phi1 = evecs[0] # PETSc Vec
obj = -lambda1 # MMA minimizes -λ₁
# ----- 7.3 Volume constraint value -----
vol_curr = fem.assemble_scalar(fem.form(rho * ufl.dx))
g_val = vol_curr / Vmax - 1.0
# ----- 7.4 Sensitivity analysis -----
sens = compute_sensitivity(rho, phi1, lambda1)
# Volume constraint gradient (linear)
grad_g = np.full(len(rho_vec), 1.0 / Vmax)
# Sensitivity filtering (optional, stabilizes optimization)
sens = sensitivity_filter(sens, rho, filter_radius)
# Objective gradient (since f0 = -λ, grad_f0 = -sens)
grad_f0 = -sens
if rank == 0:
print(f"Iter {iter_opt+1:3d}: λ₁ = {lambda1:.6e}, obj = {obj:.6e}, Vol = {vol_curr/Vmax:.4f}")
# ----- 7.5 Update MMA asymptotes -----
mma.update_asymptotes(rho_vec, iter_opt+1)
# ----- 7.6 Compute MMA subproblem coefficients -----
p0, q0, p1, q1 = mma.compute_pq(grad_f0, grad_g)
# ----- 7.7 Solve subproblem to obtain new design variables -----
rho_new = mma.solve_subproblem(p0, q0, p1, q1)
# ----- 7.8 Convergence check -----
change = np.max(np.abs(rho_new - rho_vec))
if rank == 0:
print(f" Max change = {change:.3e}")
if change < tol_change and iter_opt > 10:
if rank == 0:
print(f"Optimization converged at iteration {iter_opt+1}")
break
# Update design variables
rho_vec = rho_new
rho.x.array[:] = rho_vec
history_lambda.append(lambda1)
# Clean matrices and vectors
K_mat.destroy()
M_mat.destroy()
phi1.destroy()
================================================================
8. Save final results
================================================================
if rank == 0:
print(“\n========== Optimization completed ==========”)
# Finally compute the fundamental frequency one more time
K_final, M_final = assemble_matrices(rho)
evals_final, _ = subspace_iteration(K_final, M_final, p=1, q=6, max_iter=50, tol=1e-8)
print(f"Final fundamental frequency (circular ω) = {np.sqrt(evals_final[0]):.6f} rad/s")
# Save density field
with dolfinx.io.XDMFFile(comm, “topopt_frequency_max.xdmf”, “w”) as f:
f.write_mesh(domain)
f.write_function(rho)
print(“Optimized density field saved to topopt_frequency_max.xdmf”)
K_final.destroy()
M_final.destroy()