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, andpyfvm will create the matrix/right-hand side or the nonlinear system for it. Thispackage 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 withDirichlet boundary conditions, simply do
import meshpleximport meshzooimport numpy as npfrom scipy.sparse import linalgimport pyfvmfrom pyfvm.form_language import Boundary, dS, dV, integrate, n_dot_gradclass 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 meshzoovertices, 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 simplemesh, but anything else that provides vertices and cells works as well. For example,reading from a wide variety of mesh files is supported (viameshio):
mesh = meshplex.read("pacman.e")Likewise, PyAMG is a much faster solverfor this problem
import pyamgml = 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 andobviously the solver call is different. For Bratu's problem:
import pyfvmfrom pyfvm.form_language import *import meshzooimport numpyfrom sympy import expimport meshplexclass 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 = numpy.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 residualcomputation f.eval, e.g.,
import scipy.optimizeu = scipy.optimize.newton_krylov(f.eval, u0)Installation
pyfvm is available from the Python PackageIndex, so simply type
pip install pyfvmto install.
Testing
To run the tests, check out this repository and type
pytestLicense
This software is published under the GPLv3 license.
