There are many iterative methods for this kind of system. Finite element methods, geometric multigrid, algebraic multigrid. Stochastic methods…

## Software

### taichi

“Sparse simulator” Tai Chi is presumably also able to solve PDEs? 🤷🏼♂️ If so that would be nifty because it is also differentiable

### TenFEM

TenFEM offers a small library of differentiable FEM solvers fpr Tensorflow.

### JuliaFEM

Julaifem is an umbrella organisation supporting julia-backed FEM solvers.

### Firedrake

Firedrake seems popular for Python projects that do not need GPU?

Firedrake is an automated system for the solution of partial differential equations using the finite element method (FEM). Firedrake uses sophisticated code generation to provide mathematicians, scientists, and engineers with a very high productivity way to create sophisticated high performance simulations.

- Expressive specification of any PDE using the Unified Form Language from the FEniCS Project.
- Sophisticated, programmable solvers through seamless coupling with PETSc.…
- Sophisticated automatic optimisation, including sum factorisation for high order elements, and vectorisation.

### FEniCS

Also seems to be a friendly PDE solver. Also lacking in GPU support.

FEniCS is a popular open-source (LGPLv3) computing platform for solving partial differential equations (PDEs). FEniCS enables users to quickly translate scientific models into efficient finite element code. With the high-level Python and C++ interfaces to FEniCS, it is easy to get started, but FEniCS offers also powerful capabilities for more experienced programmers. FEniCS runs on a multitude of platforms ranging from laptops to high-performance clusters.

As an illustration of how to program a simple PDE model with FEniCS, consider the Stokes equations in variational form:

\[\int_{\Omega} \mathrm{grad} \, u : \mathrm{grad} \, v \,\mathrm{d}x \, – \int_{\Omega} p \, \mathrm{div} \, v \,\mathrm{d}x + \int_{\Omega} \mathrm{div} \, u \, q \,\mathrm{d}x = \int_{\Omega} f \cdot v \,\mathrm{d}x.\]

The variational problem is easily transcribed into Python using mathematical operators in FEniCS:

```
# Define function space
P2 = VectorElement('P', tetrahedron, 2)
P1 = FiniteElement('P', tetrahedron, 1)
TH = P2 * P1
W = FunctionSpace(mesh, TH)
# Define variational problem
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
a = inner(grad(u), grad(v))*dx - p*div(v)*dx + div(u)*q*dx
L = dot(f, v)*dx
# Compute solution
w = Function(W)
solve(a == L, w, [bc1, bc0])
```

The above code snippet also shows how to define a suitable finite element function space, using continuous piecewise quadratic vector-valued functions for the velocity and continuous piecewise linear functions for the pressure (Taylor-Hood). The computational domain and mesh are also easily created with FEniCS, here defined by three spheres immersed in a 3D channel:

```
# Define domain
h = 0.25
r = 0.3*h
box = Box(Point(0, 0, 0), Point(1, h, h))
s0 = Sphere(Point(0.3, 0.50*h, 0.50*h), r)
s1 = Sphere(Point(0.5, 0.65*h, 0.65*h), r)
s2 = Sphere(Point(0.7, 0.35*h, 0.35*h), r)
domain = box - s0 - s1 - s2
# Generate mesh
mesh = generate_mesh(domain, 32)
```

### PolyFEM

PolyFEM is a simple C++ and Python finite element library. We provide a wide set of common PDEs including:

- Laplace
- Helmholtz
- Linear Elasticity
- Saint-Venant Elasticity
- Neo-Hookean Elasticity
- Stokes
- Navier Stokes [beta]
PolyFEM simplicity lies on the interface: just pick a problem, select some boundary condition, and solve. No need to construct complicated function spaces, or learn a new scripting language: everything is set-up trough a JSON interface or through the Setting class in python.

### SfePy

SfePy is a software for solving systems of coupled partial differential equations (PDEs) by the finite element method in 1D, 2D and 3D. It can be viewed both as black-box PDE solver, and as a Python package which can be used for building custom applications. The word “simple” means that complex FEM problems can be coded very easily and rapidly.

### fipy

fipy is original school:

FiPy is an object oriented, partial differential equation (PDE) solver, written in Python, based on a standard finite volume (FV) approach. … The FiPy framework includes terms for transient diffusion, convection and standard sources, enabling the solution of arbitrary combinations of coupled elliptic, hyperbolic and parabolic PDEs. Currently implemented models include phase field treatments of polycrystalline, dendritic, and electrochemical phase transformations as well as a level set treatment of the electrodeposition process.

## PyAMG

pyamg is an algebraic multigrid solver, which is … a kind of multi-resolution method?

AMG is a multilevel technique for solving large-scale linear systems with optimal or near-optimal efficiency. Unlike geometric multigrid, AMG requires little or no geometric information about the underlying problem and develops a sequence of coarser grids directly from the input matrix. This feature is especially important for problems discretised on unstructured meshes and irregular grids:

```
A = poisson((500,500), format='csr') # 2D Poisson problem on 500x500 grid
ml = ruge_stuben_solver(A) # construct the multigrid hierarchy
print ml # print hierarchy information
b = rand(A.shape[0]) # pick a random right hand side
x = ml.solve(b, tol=1e-10) # solve Ax=b to a tolerance of 1e-8
print("residual norm is", norm(b - A*x)) # compute norm of residual vector
```

See Yousef Saad’s textbooks on algebraic multigrid methods. (free online)