pyamg.gallery#

Matrix gallery of model problems.

Functions#

  • poisson() : Poisson problem using Finite Differences

  • linear_elasticity() : Linear Elasticity using Finite Elements

  • stencil_grid() : General stencil generation from 1D, 2D, and 3D

  • diffusion_stencil_2d() : 2D rotated anisotropic FE/FD stencil

pyamg.gallery.advection_2d(grid, theta=0.7853981633974483, l_bdry=1.0, b_bdry=1.0)[source]#

Generate matrix and RHS for upwind FD discretization of 2D advection.

The 2D advection equation

(cos(theta),sin(theta)) dot grad(u) = 0,

with inflow boundaries on the left and bottom of the domain. Assume uniform grid spacing, dx=dy, even for grid[0] != grid[1].

Parameters:
gridtuple

Number of points in y and x, (ny, nx), note the ordering.

thetafloat, optional

Rotation angle theta in radians defines direction of advection (cos(theta),sin(theta)).

l_bdryfloat, array

Left boundary value. If float, then constant in-flow boundary value applied. If array, then length of array must be equal to ny=grid[0], and this array defines non-constant boundary value on the left.

b_bdryfloat, array

Bottom boundary value. If float, then constant in-flow boundary value applied. If array, then length of array must be equal to nx=grid[1], and this array defines non-constant boundary value on the bottom.

Returns:
csr_array

Defines 2D FD upwind discretization, with boundary.

array

Defines right-hand-side with boundary contributions.

See also

poisson

Examples

>>> from numpy import pi
>>> from pyamg.gallery import advection_2d
>>> A, rhs = advection_2d( (4,4), theta=pi/4)
>>> print(A.toarray().round(4))
[[ 1.4142  0.      0.     -0.7071  0.      0.      0.      0.      0.    ]
 [-0.7071  1.4142  0.      0.     -0.7071  0.      0.      0.      0.    ]
 [ 0.     -0.7071  1.4142  0.      0.     -0.7071  0.      0.      0.    ]
 [ 0.      0.      0.      1.4142  0.      0.     -0.7071  0.      0.    ]
 [ 0.      0.      0.     -0.7071  1.4142  0.      0.     -0.7071  0.    ]
 [ 0.      0.      0.      0.     -0.7071  1.4142  0.      0.     -0.7071]
 [ 0.      0.      0.      0.      0.      0.      1.4142  0.      0.    ]
 [ 0.      0.      0.      0.      0.      0.     -0.7071  1.4142  0.    ]
 [ 0.      0.      0.      0.      0.      0.      0.     -0.7071  1.4142]]
pyamg.gallery.demo()[source]#

Outline basic demo.

pyamg.gallery.diffusion_stencil_2d(epsilon=1.0, theta=0.0, type='FE')[source]#

Rotated Anisotropic diffusion in 2d of the form.

-div Q A Q^T grad u = - (C^2 + eps S^2) u_xx - 2(1 - eps) C S u_xy - (eps C^2 + S^2) u_yy

where C=cos(theta), S=sin(theta), and

Q = [cos(theta) -sin(theta)]

[sin(theta) cos(theta)]

A = [1 0 ]

[0 eps ]

Parameters:
epsilonfloat

Anisotropic diffusion coefficient: -div A grad u, where A = [1 0; 0 epsilon]. The default is isotropic, epsilon=1.0.

thetafloat

Rotation angle theta from the positive x-axis in radians. Defines -div Q A Q^T grad, where Q = [cos(theta) -sin(theta); sin(theta) cos(theta)]. The default is theta = 0.0.

type{‘FE’,’FD’}

Specifies the discretization as Q1 finite element (FE) or 2nd order finite difference (FD).

Returns:
array

Diffusion stencil of size (3,3).

See also

stencil_grid, poisson

Notes

Not all combinations are supported.

The stencil is ordered with y varying first; see stencil_grid for more details.

Examples

>>> import numpy as np
>>> import scipy as sp
>>> from pyamg.gallery.diffusion import diffusion_stencil_2d
>>> from pyamg.gallery import stencil_grid
>>> sten = diffusion_stencil_2d(epsilon=0.0001,theta=np.pi/6,type='FD')
>>> print(sten)
[[-0.2164847 -0.750025   0.2164847]
 [-0.250075   2.0002    -0.250075 ]
 [ 0.2164847 -0.750025  -0.2164847]]

Consider a 2 x 4 grid ([x0, x1] x [y0, y1, y2, y3]). The first dimension of the stencil defines x.

>>> nx, ny = (2, 4)
>>> sten = diffusion_stencil_2d(epsilon=0.1, type='FD')
>>> A = stencil_grid(sten, (nx, ny)).toarray()
>>> print(sten)
[[-0.  -1.   0. ]
 [-0.1  2.2 -0.1]
 [ 0.  -1.  -0. ]]
>>> print(A)
[[ 2.2 -0.1  0.   0.  -1.   0.   0.   0. ]
 [-0.1  2.2 -0.1  0.   0.  -1.   0.   0. ]
 [ 0.  -0.1  2.2 -0.1  0.   0.  -1.   0. ]
 [ 0.   0.  -0.1  2.2  0.   0.   0.  -1. ]
 [-1.   0.   0.   0.   2.2 -0.1  0.   0. ]
 [ 0.  -1.   0.   0.  -0.1  2.2 -0.1  0. ]
 [ 0.   0.  -1.   0.   0.  -0.1  2.2 -0.1]
 [ 0.   0.   0.  -1.   0.   0.  -0.1  2.2]]
pyamg.gallery.gauge_laplacian(npts, spacing=1.0, beta=0.1)[source]#

Construct a Gauge Laplacian from Quantum Chromodynamics for regular 2D grids.

Note that this function is not written efficiently, but should be fine for N x N grids where N is in the low hundreds.

Parameters:
nptsint

Number of pts in x and y directions.

spacingfloat

Grid spacing between points.

betafloat

Temperature. If beta=0, then the result is the typical 5pt Laplacian stencil.

Returns:
csr_array

A is Hermitian positive definite for beta > 0.0. A is Symmetric semi-definite for beta = 0.0.

References

[1]

MacLachlan, S. and Oosterlee, C., “Algebraic Multigrid Solvers for Complex-Valued Matrices”, Vol. 30, SIAM J. Sci. Comp, 2008

Examples

>>> from pyamg.gallery import gauge_laplacian
>>> A = gauge_laplacian(10)
pyamg.gallery.linear_elasticity(grid, spacing=None, E=100000.0, nu=0.3, format=None)[source]#

Linear elasticity problem with Q1 finite elements on a regular rectangular grid.

Parameters:
gridtuple

Length 2 tuple of grid sizes, e.g. (10, 10).

spacingtuple

Length 2 tuple of grid spacings, e.g. (1.0, 0.1).

Efloat

Young’s modulus.

nufloat

Poisson’s ratio.

formatstr

Format of the returned sparse matrix (eg. ‘csr’, ‘bsr’, etc.).

Returns:
csr_array

FE Q1 stiffness matrix.

array

Rigid body modes.

See also

linear_elasticity_p1

Notes

Only 2d.

References

[1]

J. Alberty, C. Carstensen, S. A. Funken, and R. KloseDOI “Matlab implementation of the finite element method in elasticity” Computing, Volume 69, Issue 3 (November 2002) Pages: 239 - 263 http://www.math.hu-berlin.de/~cc/

Examples

>>> from pyamg.gallery import linear_elasticity
>>> A, B = linear_elasticity((4, 4))
pyamg.gallery.linear_elasticity_p1(vertices, elements, E=100000.0, nu=0.3, format=None)[source]#

P1 elements in 2 or 3 dimensions.

Parameters:
verticesarray

Array of vertices of a triangle or tets.

elementsarray

Array of vertex indices for tri or tet elements.

Efloat

Young’s modulus.

nufloat

Poisson’s ratio.

formatstr

Sparse array format: ‘csr’, ‘csc’, ‘coo’, ‘bsr’.

Returns:
csr_array

FE Q1 stiffness matrix.

Notes

Both 2d and 3d.

References

[1]

J. Alberty, C. Carstensen, S. A. Funken, and R. KloseDOI “Matlab implementation of the finite element method in elasticity” Computing, Volume 69, Issue 3 (November 2002) Pages: 239 - 263 http://www.math.hu-berlin.de/~cc/

Examples

>>> from pyamg.gallery import linear_elasticity_p1
>>> import numpy as np
>>> E = np.array([[0, 1, 2],[1, 3, 2]])
>>> V = np.array([[0.0, 0.0],[1.0, 0.0],[0.0, 1.0],[1.0, 1.0]])
>>> A, B = linear_elasticity_p1(V, E)
pyamg.gallery.load_example(name)[source]#

Load an example problem by name.

Parameters:
namestr

Name of the example to load.

Returns:
dict

Dictionary with variable names and data.

Notes

Each example is stored in a dictionary with the following keys:

  • A : sparse matrix

  • B : near-nullspace candidates

  • vertices : dense array of nodal coordinates

  • elements : dense array of element indices

Current example names are:%s

  • airfoil

  • bar

  • helmholtz_2D

  • knot

  • local_disc_galerkin_diffusion

  • recirc_flow

  • unit_cube

  • unit_square

Examples

>>> from pyamg.gallery import load_example
>>> ex = load_example('knot')
pyamg.gallery.poisson(grid, dtype=<class 'float'>, format=None, type='FD')[source]#

Return a sparse matrix for the N-dimensional Poisson problem.

The matrix represents a finite difference approximation to the Poisson problem on a regular n-dimensional grid with unit grid spacing and Dirichlet boundary conditions.

The last dimension is iterated over first: z, then y, then x. This should be used with nupy.mgrid() or numpy.ndenumerate().

Parameters:
gridtuple of integers

Grid dimensions e.g. (100,100).

dtypedtype

Target dtype.

formatstr

Sparse array format to return, e.g. “csr”, “coo”, etc.

typestr

Discretization type, either finite difference (FD) or finite element (FE).

Returns:
sparray

Sparse matrix.

Notes

The matrix is symmetric and positive definite (SPD).

Examples

>>> from pyamg.gallery import poisson
>>> # 4 nodes in one dimension
>>> poisson((4,)).toarray()
array([[ 2., -1.,  0.,  0.],
       [-1.,  2., -1.,  0.],
       [ 0., -1.,  2., -1.],
       [ 0.,  0., -1.,  2.]])

>>> # rectangular two dimensional grid
>>> poisson((2,3)).toarray()
array([[ 4., -1.,  0., -1.,  0.,  0.],
       [-1.,  4., -1.,  0., -1.,  0.],
       [ 0., -1.,  4.,  0.,  0., -1.],
       [-1.,  0.,  0.,  4., -1.,  0.],
       [ 0., -1.,  0., -1.,  4., -1.],
       [ 0.,  0., -1.,  0., -1.,  4.]])
pyamg.gallery.regular_triangle_mesh(nx, ny)[source]#

Construct a regular triangular mesh in the unit square.

Parameters:
nxint

Number of nodes in the x-direction.

nyint

Number of nodes in the y-direction.

Returns:
Vertarray

Vertex list of size (nx*ny, 2).

E2Varray

Element list of size (Nex, 3).

Examples

>>> from pyamg.gallery import regular_triangle_mesh
>>> E2V,Vert = regular_triangle_mesh(3, 2)
pyamg.gallery.sprand(m, n, density, format='csr')[source]#

Return a random sparse matrix.

Parameters:
m, nint, int

Shape of the result.

densityfloat

Target a matrix density with nnz(A) = m*n*density and 0<=density<=1.

formatstr

Sparse array format to return, e.g. ‘csr’, ‘coo’, etc.

Returns:
sparray

Sparse matrix of size (m, n).

Examples

>>> from pyamg.gallery import sprand
>>> A = sprand(5,5,3/5.0)
pyamg.gallery.stencil_grid(S, grid, dtype=None, format=None)[source]#

Construct a sparse matrix form a local matrix stencil.

Parameters:
Sndarray

Matrix stencil stored in n-d array.

gridtuple

Tuple containing the grid dimensions.

dtypedtype

Data type of the result.

formatstr

Sparse array format to return, e.g. “csr”, “coo”, etc.

Returns:
sparray

Sparse array which represents the operator given by applying stencil S at each vertex of a regular grid with given dimensions.

Notes

The grid vertices are enumerated as arange(prod(grid)).reshape(grid). This implies that the last grid dimension cycles fastest, while the first dimension cycles slowest. For example, if grid=(2,3) then the grid vertices are ordered as (0,0), (0,1), (0,2), (1,0), (1,1), (1,2).

This coincides with the ordering used by the NumPy functions numpy.ndenumerate() and numpy.mgrid().

Examples

>>> from pyamg.gallery import stencil_grid
>>> stencil = [-1,2,-1]  # 1D Poisson stencil
>>> grid = (5,)          # 1D grid with 5 vertices
>>> A = stencil_grid(stencil, grid, dtype=float, format='csr')
>>> A.toarray()
array([[ 2., -1.,  0.,  0.,  0.],
       [-1.,  2., -1.,  0.,  0.],
       [ 0., -1.,  2., -1.,  0.],
       [ 0.,  0., -1.,  2., -1.],
       [ 0.,  0.,  0., -1.,  2.]])

>>> stencil = [[0,-1,0],[-1,4,-1],[0,-1,0]] # 2D Poisson stencil
>>> grid = (3,3)                            # 2D grid with shape 3x3
>>> A = stencil_grid(stencil, grid, dtype=float, format='csr')
>>> A.toarray()
array([[ 4., -1.,  0., -1.,  0.,  0.,  0.,  0.,  0.],
       [-1.,  4., -1.,  0., -1.,  0.,  0.,  0.,  0.],
       [ 0., -1.,  4.,  0.,  0., -1.,  0.,  0.,  0.],
       [-1.,  0.,  0.,  4., -1.,  0., -1.,  0.,  0.],
       [ 0., -1.,  0., -1.,  4., -1.,  0., -1.,  0.],
       [ 0.,  0., -1.,  0., -1.,  4.,  0.,  0., -1.],
       [ 0.,  0.,  0., -1.,  0.,  0.,  4., -1.,  0.],
       [ 0.,  0.,  0.,  0., -1.,  0., -1.,  4., -1.],
       [ 0.,  0.,  0.,  0.,  0., -1.,  0., -1.,  4.]])