pyfvm

Creating finite volume equation systems with ease.

pyfvm provides everything that is needed for setting up finite volume equation
systems. The user needs to specify the finite volume formulation in a
configuration file, and pyfvm will create the matrix/right-hand side or the
nonlinear system for it. This package is for everyone who wants to quickly
construct FVM systems.

Examples

Linear equation systems

pyfvm works by specifying the residuals, so for solving Poisson’s equation with
Dirichlet boundary conditions, simply do

import meshplex
import meshzoo
import numpy as np
from scipy.sparse import linalg
import pyfvm
from pyfvm.form_language import Boundary, dS, dV, integrate, n_dot_grad
class Poisson:
    def apply(self, u):
        return integrate(lambda x: -n_dot_grad(u(x)), dS) - integrate(lambda x: 1.0, dV)
    def dirichlet(self, u):
        return [(lambda x: u(x) - 0.0, Boundary())]
# Create mesh using meshzoo
vertices, cells = meshzoo.rectangle_tri(
    np.linspace(0.0, 2.0, 401), np.linspace(0.0, 1.0, 201)
)
mesh = meshplex.Mesh(vertices, cells)
matrix, rhs = pyfvm.discretize_linear(Poisson(), mesh)
u = linalg.spsolve(matrix, rhs)
mesh.write("out.vtk", point_data={"u": u})

This example uses meshzoo for creating a
simple mesh, but anything else that provides vertices and cells works as well.
For example, reading from a wide variety of mesh files is supported (via
meshio):

mesh = meshplex.read("pacman.e")

Likewise, PyAMG is a much faster solver for
this problem

import pyamg
ml = pyamg.smoothed_aggregation_solver(matrix)
u = ml.solve(rhs, tol=1e-10)

More examples are contained in the examples directory.

Nonlinear equation systems

Nonlinear systems are treated almost equally; only the discretization and
obviously the solver call is different. For Bratu’s problem:

import pyfvm
from pyfvm.form_language import *
import meshzoo
import numpy as np
from sympy import exp
import meshplex
class Bratu:
    def apply(self, u):
        return integrate(lambda x: -n_dot_grad(u(x)), dS) - integrate(
            lambda x: 2.0 * exp(u(x)), dV
        )
    def dirichlet(self, u):
        return [(u, Boundary())]
vertices, cells = meshzoo.rectangle_tri(
    np.linspace(0.0, 2.0, 101), np.linspace(0.0, 1.0, 51)
)
mesh = meshplex.Mesh(vertices, cells)
f, jacobian = pyfvm.discretize(Bratu(), mesh)
def jacobian_solver(u0, rhs):
    from scipy.sparse import linalg
    jac = jacobian.get_linear_operator(u0)
    return linalg.spsolve(jac, rhs)
u0 = np.zeros(len(vertices))
u = pyfvm.newton(f.eval, jacobian_solver, u0)
mesh.write("out.vtk", point_data={"u": u})

Note that the Jacobian is computed symbolically from the Bratu class.

Instead of pyfvm.newton, you can use any solver that accepts the residual
computation f.eval, e.g.,

import scipy.optimize
u = scipy.optimize.newton_krylov(f.eval, u0)