proteus.LinearAlgebraTools module

Tools for n-dimensional linear algebra

Vectors are just numpy arrays, as are dense matrices. Sparse matrices are CSR matrices. Parallel vector and matrix are built on top of those representations using PETSc.

Inheritance diagram of proteus.LinearAlgebraTools

proteus.LinearAlgebraTools.petsc_load_matrix(filename)[source]

This function loads a PETSc matrix from a binary format. (Eg. what is saved using the petsc_view function).

Parameters:filename (str) – This is the name of the binary with the file stored.
Returns:matrix – The matrix that is stored in the binary file.
Return type:petsc4py matrix
proteus.LinearAlgebraTools.petsc_load_vector(filename)[source]

This function loads a PETSc vector from a binary format. (Eg. what is saved using the petsc_view function).

Parameters:filename (str) – This is the name of the binary with the file stored.
Returns:matrix – The matrix that is stored in the binary file.
Return type:petsc4py vector
proteus.LinearAlgebraTools.csr_2_petsc(size, csr)[source]

Create an petsc4py matrix from size and CSR information.

size
: tuple
A 2-tuple with the number of matrix rows and columns.
csr
: tuple
A 3-tuple with the sparse matrix csr information.

matrix : PETSc4py aij matrix

proteus.LinearAlgebraTools.superlu_sparse_2_dense(sparse_matrix, output=False)[source]

Converts a sparse superluWrapper into a dense matrix.

Parameters:
  • sparse_matrix
  • output (str) – Out file name to store the matrix.
Returns:

dense_matrix – A numpy array storing the dense matrix.

Return type:

numpy array

Notes

This function should not be used for large matrices.

proteus.LinearAlgebraTools.petsc4py_sparse_2_dense(sparse_matrix, output=False)[source]

Converts a PETSc4Py matrix to a dense numpyarray.

Parameters:
  • sparse_matrix (PETSc4py matrix) –
  • output (str) – Output file name to store the matrix.
Returns:

dense_matrix – A numpy array with the dense matrix.

Return type:

numpy array

Notes

This function is very inefficient for large matrices.

proteus.LinearAlgebraTools.superlu_2_petsc4py(sparse_superlu)[source]

Copy a sparse superlu matrix to a sparse petsc4py matrix

Parameters:sparse_superlu (proteus.superluWrappers.SparseMatrix) –
Returns:sparse_matrix
Return type:
class:p4pyPETSc.Mat
proteus.LinearAlgebraTools.petsc_create_diagonal_inv_matrix(sparse_petsc)[source]

Create an inverse diagonal petsc4py matrix from input matrix.

Parameters:sparse_petsc (p4pyPETSc.Mat) –
Returns:sparse_matrix
Return type:p4pyPETSc.Mat
proteus.LinearAlgebraTools.dense_numpy_2_petsc4py(dense_numpy, eps=1e-12)[source]

Create a sparse petsc4py matrix from a dense numpy matrix.

Note - This routine has been built mainly to support testing. It would be rare for this routine to be useful for most applications.

Parameters:
  • dense_numpy
  • eps (float) – Tolerance for non-zero values.
Returns:

sparse_matrix

Return type:

PETSc4py matrix

proteus.LinearAlgebraTools.csr_2_petsc_mpiaij(size, csr)[source]

Create an MPIaij petsc4py matrix from size and CSR information.

size
: tuple
Two entires: (num_rows, num_cols)
csr
: tuple
(row_idx, col_idx, vals)

matrix : PETSc4py MPIaij matrix

class proteus.LinearAlgebraTools.ParVec(array, blockSize, n, N, nghosts=None, subdomain2global=None, blockVecType='simple')[source]

A parallel vector built on top of daetk’s wrappers for petsc

scatter_forward_insert()[source]
scatter_reverse_add()[source]
class proteus.LinearAlgebraTools.ParVec_petsc4py(array=None, bs=None, n=None, N=None, nghosts=None, subdomain2global=None, blockVecType='simple', ghosts=None, proteus2petsc_subdomain=None, petsc2proteus_subdomain=None)[source]

Bases: petsc4py.PETSc.Vec

Parallel vector using petsc4py’s wrappers for PETSc

Parameters:
  • array (numpy_array) – A numpy array with size equal to the number of locally owned unknowns plus the number of local ghost cells.
  • bs (int) – Block size.
  • n (int) – The number of locally owned unknowns
  • N (int) – The number of unknowns in the global system
  • nghosts (int) – The number of ghost nodes for the process.
  • subdomain2global (numpy array) – Map from the process unknowns to the global uknowns.
  • blockVecType (str) –
  • ghosts (numpy array) – A numpy array with the local process uknowns that are ghost nodes.
  • proteus2petsc_subdomain (numpy array) – A numpy array that serves as a map from the proteus uknown ordering to the petsc uknown ordering
  • petsc2proteus_subdomain (numpy array) – A numpy array that serves as a map from the petsc uknown ordering to the proteus unknown ordering
scatter_forward_insert()[source]
scatter_reverse_add()[source]
save(filename)[source]

Saves to disk using a PETSc binary viewer.

class proteus.LinearAlgebraTools.ParInfo_petsc4py[source]

ARB - this class is experimental. My idea is to store the information need to constructor parallel vectors and matrices here as static class values. Then ParVec and ParMat can use these values to create parallel objects later.

print_info()[source]
class proteus.LinearAlgebraTools.ParMat_petsc4py(ghosted_csr_mat=None, par_bs=None, par_n=None, par_N=None, par_nghost=None, subdomain2global=None, blockVecType='simple', pde=None, par_nc=None, par_Nc=None, proteus_jacobian=None, nzval_proteus2petsc=None)[source]

Bases: petsc4py.PETSc.Mat

Parallel matrix based on petsc4py’s wrappers for PETSc. ghosted_csr_mat : proteus.superluWrappers.SparseMatrix

Primary CSR information for the ParMat.
par_bs
: int
The block size.
par_n
: int
The number of locally owned unknowns.
par_N
: int
The number of global unknowns.
par_nghost
: int
The number of locally owned ghost unknowns.
subdomain2global
: numpy.ndarray
A map from the local unknown to the global unknown.

blockVecType : str pde : proteus.Transport.OneLevelTransport

The Transport class defining the problem.

par_nc : int par_Nc : int proteus_jacobian : proteus.superluWrappers.SparseMatrix

Jacobian generated by Transport class’s initializeJacobian.
nzval_proteus2petsc
: numpy.ndarray
Array with index permutations for mapping between proteus and petsc degrees of freedom.
classmethod create_ParMat_from_OperatorConstructor(operator)[source]

Build a ParMat consistent with the problem from an Operator constructor matrix.

Parameters:operator (proteus.superluWrappers.SparseMatrix) – Matrix to be turned into a parallel petsc matrix.
save(filename)[source]

Saves to disk using a PETSc binary viewer.

proteus.LinearAlgebraTools.Vec(n)[source]

Build a vector of length n (using numpy)

For example:

>>> Vec(3)
array([ 0.,  0.,  0.])
proteus.LinearAlgebraTools.Mat(m, n)[source]

Build an m x n matrix (using numpy)

For example:

>>> Mat(2,3)
array([[ 0.,  0.,  0.],
      [ 0.,  0.,  0.]])
proteus.LinearAlgebraTools.SparseMatFromDict(nr, nc, aDict)[source]

Build a nr x nc sparse matrix from a dictionary representation

proteus.LinearAlgebraTools.SparseMat(nr, nc, nnz, nzval, colind, rowptr)[source]

Build a nr x nc sparse matrix from the CSR data structures

Parameters:
  • nr (int) – The number of rows.
  • nc (int) – The number of columns.
  • nnz (int) – The number of non-zero matrix entries.
  • nzval (numpy array) – Array with non-zero matrix entries.
  • colind (numpy array of 32bit integers) – CSR column array.
  • rowptr (numpy array of 32bit integers) – CSR row pointer.
Returns:

sparse_matrix – superlu sparse matrix in CSR format.

Return type:

proteus.superluWrappers.SparseMatrix

Note

For the superluWrapper, both the colind and rowptr should use 32-bit integer data types.

class proteus.LinearAlgebraTools.SparseMatShell(ghosted_csr_mat)[source]

Build a parallel matrix shell from CSR data structures.

Parameters:ghosted_csr_mat
create(A)[source]
mult(A, x, y)[source]
class proteus.LinearAlgebraTools.OperatorShell[source]

A base class for operator shells

create(A)[source]
class proteus.LinearAlgebraTools.ProductOperatorShell[source]

Bases: proteus.LinearAlgebraTools.OperatorShell

A base class for shell operators that apply multiplcation.

Operators derived from this class should have working multiplication functions.

mult(A, x, y)[source]
class proteus.LinearAlgebraTools.InvOperatorShell[source]

Bases: proteus.LinearAlgebraTools.OperatorShell

A base class for inverse operator shells

Operators derived from this class should have working apply functions.

apply(A, x, y)[source]
class proteus.LinearAlgebraTools.LSCInv_shell(Qv, B, Bt, F)[source]

Bases: proteus.LinearAlgebraTools.InvOperatorShell

Shell class for the LSC Inverse Preconditioner

This class creates a shell for the least-squares commutator (LSC) preconditioner, where \(M_{s}= (B \hat{Q^{-1}_{v}} B^{'}) (B \hat{Q^{-1}_{v}} F \hat{Q^{-1}_{v}} B^{'})^{-1} (B \hat{Q^{-1}_{v}} B^{'})\) is used to approximate the Schur complement.

Initializes the LSC inverse operator.

Parameters:
  • Qv (petsc4py matrix object) – The diagonal elements of the velocity mass matrix.
  • B (petsc4py matrix object) – The discrete divergence operator.
  • Bt (petsc4py matrix object) – The discrete gradient operator.
  • F (petsc4py matrix object) – The A-block of the linear system.
apply(A, x, y)[source]

Apply the LSC inverse operator

Parameters:
  • A (NULL) – A necessary PETSc4py placeholder for internal function operations.
  • x (PETSc4py vector) – The incoming residual vector which the operator is applied to.
Returns:

y – The result of the preconditioners action on x.

Return type:

PETSc4py vector

class proteus.LinearAlgebraTools.MatrixShell(A)[source]

Bases: proteus.LinearAlgebraTools.ProductOperatorShell

A shell class for a matrix.

Specifies a basic matrix shell.

Parameters:A (matrix) – A petsc4py matrix object
mult(A, x, y)[source]

Multiply the matrix and x.

Parameters:
  • A (matrix) – Dummy place holder for PETSc compatibility
  • x (vector) –
Returns:

y

Return type:

vector

class proteus.LinearAlgebraTools.MatrixInvShell(A)[source]

Bases: proteus.LinearAlgebraTools.InvOperatorShell

A PETSc shell class for a inverse operator.

Initializes operators and solvers for inverse operator.

Parameters:A (PETSc matrix) – This is the matrix object used to construct the inverse.
apply(A, x, y)[source]

Apply the inverse pressure mass matrix.

Parameters:
  • A (matrix) – Dummy place holder for PETSc compatibility
  • x (vector) –
Returns:

y

Return type:

vector

class proteus.LinearAlgebraTools.Sp_shell(A00, A11, A01, A10, constNullSpace=True)[source]

Bases: proteus.LinearAlgebraTools.InvOperatorShell

Shell class for the SIMPLE preconditioner which applies the following action.

\[\hat{S}^{-1} = (A_{11} - A_{01} \text{diag}(A_{00} A_{10})^{-1}\]

where \(A_{ij}\) are sub-blocks of the global saddle point system.

Parameters:
  • A00 (petsc4py matrix) – The A00 block of the global saddle point system.
  • A01 (petsc4py matrix) – The A01 block of the global saddle point system.
  • A10 (petsc4py matrix) – The A10 block of the global saddle point system.
  • A11 (petsc4py matrix) – The A11 block of the global saddle point system.
  • use_constant_null_space (bool) – Indicates whether a constant null space should be used. See note below.

Notes

For Stokes or Navier-Stokes systems, the \(S_{p}\) operator resembles a Laplcian matrix on the pressure. As with other Laplacian type operators, S_{p} has a constant null space which must be accounted for when the operator is being applied. As such, this operator’s default behavior is to use a constant null space. Should a circumstance arise in which a constant null space is not appropriate, it can be disabled with the use_constant_null_space flag.

apply(A, x, y)[source]
class proteus.LinearAlgebraTools.TwoPhase_PCDInv_shell(Qp_visc, Qp_dens, Ap_rho, Np_rho, alpha=False, delta_t=0, num_chebyshev_its=0)[source]

Bases: proteus.LinearAlgebraTools.InvOperatorShell

Shell class for the two-phase PCD preconditioner. The two-phase PCD_inverse shell applies the following operator.

\[\hat{S}^{-1} = (Q^{(1 / \mu)})^{-1} + (A_{p}^{(1 / \rho)})^{-1} (N_{p}^{(\rho)} + \dfrac{\alpha}{\Delta t} Q^{(\rho)} ) (Q^{(\rho)})^{-1}\]

where \(Q^{(1 / \mu)}\) and \(Q^{(\rho)}\) denote the pressure mass matrix scaled by the inverse dynamic viscosity and density respectively, \((A_{p}^{(1 / \rho)})^{-1}\) denotes the pressure Laplacian scaled by inverse density, and \(N_{p}^{(\rho)}\) denotes the pressure advection operator scaled by the density, and \(\alpha\) is a binary operator indicating whether the problem is temporal or steady state.

Initialize the two-phase PCD inverse operator.

Parameters:
  • Qp_visc (petsc4py matrix) – The pressure mass matrix with dynamic viscocity scaling.
  • Qp_dens (petsc4py matrix) – The pressure mass matrix with density scaling.
  • Ap_rho (petsc4py matrix) – The pressure Laplacian scaled with density scaling.
  • Np_rho (petsc4py matrix) – The pressure advection operator with inverse density scaling.
  • alpha (binary) – True if problem is temporal, False if problem is steady state.
  • delta_t (float) – Time step parameter.
  • num_chebyshev_its (int) – Number of chebyshev iteration steps to take. (0 indicates the chebyshev semi iteration is not used)
apply(A, x, y)[source]

Applies the two-phase pressure-convection-diffusion preconditioner.

Parameters:
  • A (None) – Dummy variabled needed to interface with PETSc
  • x (petsc4py vector) – Vector to which operator is applied
Returns:

y – Result of operator acting on x.

Return type:

petsc4py vector

proteus.LinearAlgebraTools.l2Norm(x)[source]

Compute the parallel \(l_2\) norm

proteus.LinearAlgebraTools.l1Norm(x)[source]

Compute the parallel \(l_1\) norm

The \(l_1\) norm of a vector \(\mathbf{x} \in \mathbb{R}^n\) is

\[\| \mathbf{x} \|_{1} = \sum_{i=0} |x_i|\]

If Python is running in parallel, then the sum is over all dimensions on all processors so that the input must not contain “ghost” entries.

This implemtation works for a distributed array with no ghost components (each component must be on a single processor).

Parameters:x – numpy array of length n
Returns:float
proteus.LinearAlgebraTools.lInfNorm(x)[source]

Compute the parallel \(l_{\infty}\) norm

The \(l_{\infty}\) norm of a vector \(\mathbf{x} \in \mathbb{R}^n\) is

\[\|x\|_{\infty} = \max_i |x_i|\]

This implemtation works for a distributed array with no ghost components (each component must be on a single processor).

Parameters:x – numpy array of length n
Returns:float
proteus.LinearAlgebraTools.wDot(x, y, h)[source]

Compute the parallel weighted dot product of vectors x and y using weight vector h.

The weighted dot product is defined for a weight vector \(\mathbf{h}\) as

\[(\mathbf{x},\mathbf{y})_h = \sum_{i} h_{i} x_{i} y_{i}\]

All weight vector components should be positive.

Parameters:x,y,h – numpy arrays for vectors and weight
Returns:the weighted dot product
proteus.LinearAlgebraTools.wl2Norm(x, h)[source]

Compute the parallel weighted l_2 norm with weight h

proteus.LinearAlgebraTools.wl1Norm(x, h)[source]

Compute the parallel weighted l_1 norm with weight h

proteus.LinearAlgebraTools.wlInfNorm(x, h)[source]

Compute the parallel weighted l_{infty} norm with weight h

proteus.LinearAlgebraTools.energyDot(x, y, A)[source]

Compute the “energy” dot product x^t A y (not parallel)

proteus.LinearAlgebraTools.energyNorm(x, A)[source]

Compute the “energy” norm x^t A x (not parallel)

proteus.LinearAlgebraTools.l2NormAvg(x)[source]

Compute the arithmetic averaged l_2 norm (root mean squared norm)

proteus.LinearAlgebraTools.rmsNorm(x)[source]

Compute the arithmetic averaged l_2 norm (root mean squared norm)

proteus.LinearAlgebraTools.l2Norm_local(x)[source]

Compute the l_2 norm for just local (processor) system (not parallel)

class proteus.LinearAlgebraTools.WeightedNorm(shape, atol, rtol)[source]

Compute the weighted norm for time step control (not currently parallel)

setWeight(y)[source]
norm(y, type)[source]