Choosing solvers#

Several simulations available in SimPEG need to numerically solve a partial differential equations system (PDE), such as Simulation3DNodal (DC) Simulation3DMagneticFluxDensity (TDEM) and Simulation3DMagneticFluxDensity (FDEM). A numerical solver is needed to solve the PDEs. SimPEG can make use of the solvers available in pymatsolver, like pymatsolver.Pardiso, pymatsolver.Mumps or pymatsolver.SolverLU. The choice of an appropriate solver can affect the computation time required to solve the PDE. Generally we recommend using direct solvers over iterative solvers for SimPEG, but be aware that direct solvers have much larger memory requirements.

  • The Pardiso solver wraps the oneMKL PARDISO solver available for x86_64 CPUs.

  • The Mumps solver wraps MUMPS, a fast solver available for all CPU brands, including Apple silicon architecture.

  • The SolverLU wraps SciPy’s scipy.sparse.linalg.splu(). The performance of this solver is not up to the level of Mumps and Pardiso. Usage of the SolveLU is recommended only when it’s not possible to use other faster solvers.

The default solver#

We can use simpeg.utils.get_default_solver() to obtain a reasonable default solver available for our system:

import simpeg
import simpeg.electromagnetics.static.resistivity as dc

# Get default solver
solver = simpeg.utils.get_default_solver()
print(f"Solver: {solver}")

which would print out on an x86_64 cpu:

Solver: <class 'pymatsolver.direct.pardiso.Pardiso'>

We can then use this solver in a simulation:

# Define a simple mesh
h = [(1.0, 10)]
mesh = discretize.TensorMesh([h, h, h], origin="CCC")

# And a DC survey
receiver = dc.receivers.Dipole(locations_m=(-1, 0, 0), locations_n=(1, 0, 0))
source = dc.sources.Dipole(
    receiver_list=[receiver], location_a=(-2, 0, 0), location_b=(2, 0, 0)
)
survey = dc.Survey([source])

# Use the default solver in the simulation
simulation = dc.Simulation3DNodal(mesh=mesh, survey=survey, solver=solver)

Note

The priority list used to choose a default solver is:

  1. pymatsolver.Pardiso

  2. pymatsolver.Mumps

  3. pymatsolver.SolverLU

Setting solvers manually#

Alternatively, we can manually set a solver. For example, if we want to use MUMPS in our DC resistivity simulation, we can import pymatsolver.Mumps and pass it to our simulation:

import simpeg.electromagnetics.static.resistivity as dc
from pymatsolver import Mumps

# Manually set Mumps as our solver
simulation = dc.Simulation3DNodal(mesh=mesh, survey=survey, solver=Mumps)

Note

When sharing this code with a colleague, keep in mind that it might not work if MUMPS is not available in their system.

For such scenarios, we recommend using the simpeg.utils.get_default_solver() function, that will always return a suitable solver for the current system.

Ultimately, choosing the best solver is a mixture of the problem you are solving and your current system. Experiment with different solvers yourself to choose the best.