PETSc solver crashes with MPI: `_solver_petsc` passes full global matrix without partitioning

[ricardo0115]

Description

When using set_solver('petsc', ...) with MPI (e.g., mpiexec -n 4 python script.py), the PETSc solver crashes because _solver_petsc in fedoo/core/base.py passes the full global CSR matrix to PETSc.Mat().createAIJWithArrays() on every rank:

# fedoo/core/base.py, line ~675
A_petsc = PETSc.Mat().createAIJWithArrays(A.shape, (A.indptr, A.indices, A.data))

When PETSc is initialized with MPI, createAIJWithArrays expects each rank to pass only its local rows. Since fedoo passes the full A.indptr (length n_rows + 1) on every rank, but PETSc distributes the rows across ranks, each rank expects local_rows + 1 entries in indptr and gets n_rows + 1 instead.

Error

File "petsc4py/PETSc/Mat.pyx", line 1077, in petsc4py.PETSc.Mat.createAIJWithArrays
ValueError: A matrix with 9701 rows requires a row pointer of length 9702 (given: 38805)

(38804 = total DOFs, 9701 = 38804 / 4 ranks)

Reproducer

import fedoo as fd
import numpy as np

mesh = fd.mesh.hole_plate_mesh(nr=50, nt=50, length=100, height=100, radius=20, elm_type="quad4", sym=False)
props = np.array([1e5, 0.3, 1e-5, 300, 1000, 0.3])
material = fd.constitutivelaw.Simcoon("EPICP", props)
material.use_elastic_lt = True

fd.Assembly.delete_memory()
fd.ModelingSpace("2Dstress")
center = mesh.nearest_node(mesh.bounding_box.center)

wf = fd.weakform.StressEquilibrium(material)
wf.fbar = True
assemb = fd.Assembly.create(wf, mesh, mesh.elm_type)

pb = fd.problem.NonLinear(assemb)
pb.set_solver("petsc", solver_type="preonly", pc_type="lu", pc_factor_mat_solver_type="mumps")

bc_periodic = fd.constraint.PeriodicBC(periodicity_type="small_strain")
pb.bc.add(bc_periodic)
pb.bc.add("Dirichlet", center, "Disp", 0, name="center")

pb.bc.add("Dirichlet", "E_xx", 0.005, name="_Strain")
pb.bc.add("Dirichlet", "E_yy", -0.003, name="_Strain")
pb.bc.add("Dirichlet", "E_xy", 0.001, name="_Strain")

pb.nlsolve(dt=0.1, t0=0, tmax=1, update_dt=True, print_info=1, interval_output=0.1)

# Single rank works fine:
python reproducer.py

# Multiple ranks crashes:
mpiexec -n 4 python reproducer.py
# ValueError: A matrix with 9701 rows requires a row pointer of length 9702 (given: 38805)

Root cause

The set_solver docstring mentions MPI support via PCMPI:

To use petsc with MPI parallelization (PCMPI), the script mpi4py needs to be installed, and the script should be launched using mpiexec.

However, _solver_petsc does not account for MPI at all: it passes the full global scipy sparse matrix to PETSc without any row partitioning. For MPI to work, the function would need to either:

1. Partition the CSR arrays per rank before calling createAIJWithArrays, or
2. Use PETSc.Mat().createAIJ() with comm=PETSc.COMM_SELF to force sequential operation on each rank (which would bypass the benefit of MPI), or
3. Use PETSc's PCMPI preconditioner, which runs a sequential KSP on rank 0 and distributes only the preconditioner — but this still requires the matrix to be created correctly on the communicator.

Since the assembly itself is not distributed (no mesh partitioning, no distributed DOF numbering), option (1) would require significant architectural changes beyond just the solver function.

Context

This was found while benchmarking parallelization strategies for nonlinear FEM simulations (elastoplastic periodic homogenization, ~20k nodes) on an HPC cluster with PBS Pro. The single-rank PETSc path works correctly (both MUMPS direct and CG iterative), and the OpenMP parallelization via pypardiso/MKL works well (4x speedup at 4–8 threads). The MPI path is the only configuration that fails.

Environment

• fedoo 0.7.0
• petsc 3.24.5, petsc4py, mpi4py, OpenMPI
• Python 3.12, Ubuntu 24.04 (Apptainer container)

Suggestion

At minimum, the set_solver docstring should clarify that MPI parallelization of the linear solve is not currently supported (the existing PCMPI note is misleading). Ideally, _solver_petsc could detect MPI.COMM_WORLD.size > 1 and either raise a clear error or fall back to COMM_SELF.