This repo contains discretization methods for the elliptic PDE
,
in two spatial dimensions that handle general quadrilateral grids. In particular, the MPFA-O and MPFA-L methods are implemented, they yield locally mass conservative discretizations that are consistent for rough grids. The MPFA-L method has particulary good monotonicity properties, i.e., the maximum principle is respected for a wide range of grids. It's therefore considered as the state of the art method for porous media flow problems on quadrilatereal grids.
The code requires numpy for assembling discretization matrices, for running the convergence tests one also needs scipy and sympy.
The way to define a domain and discretize it into control volumes, is to discretize the unit square with rectangles, then perturb the grid points. This approach is very flexible and alows for complicated domains and control volume discretizations.
from discretization.mesh import Mesh
import numpy as np
nx = ny = 6 #Number of grid points in x and y direction on the unit square
perturbation = lambda p:np.array([p[0],0.5*p[0]+p[1]]) #perturbation on every grid point p
mesh = Mesh(nx,ny,perturbation,ghostboundary=True) #Creates a mesh object with a strip of ghostcells for boundary handling
mesh.plot()This would result in the parallelogram discretization, note that we have 8 grid points (in orange) in each direction, 2 more than 6, as we have a strip of ghost cells.
For solving the Poisson equation on this, one would define the problem data with numpy arrays and python functions, then pass it to the compute_matrix and compute_vector functions, together with the mesh object. In this example we have a homogeneous domain (permeability is a matrix of ones), and a isotropic medium (tensor is diagonal).
from discretization.FVML import compute_matrix,compute_vector
import math
source = lambda x , y : math.sin(y)*math.cos(x)
boundary_condition = lambda x , y :0
tensor = np.eye(2)
permeability = np.ones(( mesh.num_unknowns,mesh.num_unknowns))
A = np.zeros((mesh.num_unknowns,mesh.num_unknowns))# stiffness matrix
f = np.zeros(mesh . num_unknowns)# load vector
compute_matrix(mesh,A,tensor,permeability)
compute_vector(mesh,f,source,boundary_condition)
u = np.linalg.solve(A,f)
mesh.plot_vector(u)This would result in the solution
For more interesting equations, such as the time dependent, non-linear Richards' equation, see this file.