scipy 3-D steady conduction: spsolve stalls on 100k+ voxel grids, Jacobi-preconditioned CG solves in seconds
Context: steady heat conduction in a heated copper plate cooled by drilled coolant holes, a 7-point finite-volume grid of about 150k to 300k cells at 0.1 mm (1.2M to 2.5M at 0.05 mm).
Problem: scipy.sparse.linalg.spsolve (SuperLU) had not finished a batch of these 3-D cases after 10 minutes. Fill-in makes direct factorisation of 3-D Laplacians very slow.
Fix: keep the matrix symmetric positive definite and use conjugate gradient with a Jacobi preconditioner:
Mc = M.tocsr(); dinv = 1.0 / Mc.diagonal() pre = spla.LinearOperator(Mc.shape, matvec=lambda v: dinv * v) T, info = spla.cg(Mc, b, x0=np.full(N, guess), rtol=1e-9, maxiter=40000, M=pre) assert info == 0
Two solves at about 156k cells took about 1.5 s in total, and the half-spacing grids (up to 2.5M cells) still ran in reasonable time. In SciPy 1.12 and later the keyword is rtol; tol was removed.
To keep the matrix symmetric, do not put fixed-temperature cells (coolant voxels, Dirichlet walls) in the metal rows. For a metal cell next to a coolant cell, add the series conductance g = A / (dx/(2k) + 1/h) to the metal cell's diagonal and g*T_coolant to its right-hand side. Give the coolant cells identity rows with b = T_coolant. A heat flux boundary goes only into b.
Check: halving the grid spacing changed peak temperatures by about 2 K, so 0.1 mm was converged for holes of about 0.4 to 0.9 mm.
Side lesson: a quick 1-D fin estimate that holds all the heat in a thin surface layer overstated the mid-gap temperature between widely spaced cooling holes by several times. In 3-D the heat spreads deeper before it turns toward the holes. Run the 3-D model before rejecting a wider hole spacing on thermal grounds.
Reproduced on a comparable model (SciPy 1.18.1): 7-point FV, copper k=400, two drilled holes with your series-conductance coupling, identity coolant rows, flux on top; max|M - M.T| = 0. At N=157,464: spsolve (SuperLU, COLAMD) 118 s, splu(permc_spec='MMD_AT_PLUS_A', SymmetricMode=True) 41 s, cg 0.37 s. At N=512,000, splu was killed for memory and cg took 2.0 s. So the main point holds. Three corrections. (1) On a single-material grid the Jacobi preconditioner does nothing: the diagonal is almost constant, giving 315 vs 324 iterations at 157k and 454 vs 459 at 512k, and Jacobi-CG was slightly slower in wall time. It only helps when k varies strongly between materials. If iteration count matters, pyamg.smoothed_aggregation_solver(M).aspreconditioner() takes 8-9 iterations, with similar wall time at these sizes because of the setup cost. (2) tol was deprecated in 1.12 and removed in 1.14: 1.12.0 and 1.13.1 raise DeprecationWarning, 1.14.1 raises TypeError. (3) Timing caveat: on a shared 8-core box at load ~4-5, default OpenBLAS threading made the same cg solve 13x slower (157k: 4.9 s vs 0.37 s with OPENBLAS_NUM_THREADS=1). The vector dot products and norms inside cg oversubscribe: 1000 dots of a 27k-element vector took 3.0 s with threads unset or 8, and 0.004 s with 1-4 threads. Set the BLAS thread count before you benchmark iterative solvers. Accuracy is fine: rtol=1e-9 gave max|T - T_direct| = 5e-5 K.