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.

proteus.LinearAlgebraTools.petsc_load_matrix(filename)

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)

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.petsc_load_IS(filename)

This function loads a PETSc index-set 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 index-set that is stored in the binary file.

Return type

petsc4py IS

proteus.LinearAlgebraTools.csr_2_petsc(size, csr)

Create an petsc4py matrix from size and CSR information.

sizetuple

A 2-tuple with the number of matrix rows and columns.

csrtuple

A 3-tuple with the sparse matrix csr information.

matrix : PETSc4py aij matrix

proteus.LinearAlgebraTools.superlu_get_rank(sparse_matrix)

Returns the rank of a superluWrapper sparse matrix.

Parameters

sparse_matrix (proteus.superluWrappers.SparseMatrix) –

Returns

matrix_rank – The rank of the sparse_matrix

Return type

int

Notes

This function is a tool for debugging and should only be used for small matrices.

proteus.LinearAlgebraTools.petsc4py_get_rank(sparse_matrix)

Returns the rank of a superluWrapper sparse matrix.

Parameters

sparse_matrix (p4pyPETSc.Mat) –

Returns

matrix_rank – The rank of the sparse_matrix

Return type

int

Notes

This function is a debugging tool and should only be used for small matrices.

proteus.LinearAlgebraTools.superlu_has_pressure_null_space(sparse_matrix)

Checks whether a superluWrapper sparse matrix has a constant pressure null space.

Parameters

sparse_matrix (proteus.superluWrappers.SparseMatrix) –

Returns

does – Boolean variable indicating whether the pressure term creates a null space.

Return type

bool

Notes

Assumes interwoven dof. This function was written mainly for debugging purposes and may be slow for large matrices.

proteus.LinearAlgebraTools.petsc4py_mat_has_pressure_null_space(A)

Checks whether a PETSc4Py sparse matrix has a constant pressure null space.

Parameters

A (p4pyPETSc.Mat) –

Returns

does – Boolean variable indicating whether the pressure term creates a null space.

Return type

bool

Notes

Assumes interwoven dof. This function was written mainly for debugging purposes and may be slow for large matrices.

proteus.LinearAlgebraTools.superlu_sparse_2_dense(sparse_matrix, output=False)

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)

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)

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)

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)

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)

Create an MPIaij petsc4py matrix from size and CSR information.

sizetuple

Two entires: (num_rows, num_cols)

csrtuple

(row_idx, col_idx, vals)

matrix : PETSc4py MPIaij matrix

proteus.LinearAlgebraTools.split_PETSc_Mat(mat)

Decompose a PETSc matrix into a symmetric and skew-symmetric

matrix

mat : :class: PETSc4py Matrix

H:class: PETSc4py Matrix

Symmetric (or Hermitian) component of mat

S:class: PETSc4py Matrix

Skew-Symmetric (or skew-Hermitian) component of mat

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)

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()

scatter_reverse_add()

save(filename)

Saves to disk using a PETSc binary viewer.

class proteus.LinearAlgebraTools.ParInfo_petsc4py

Bases: object

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()

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)

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_bsint

The block size.

par_nint

The number of locally owned unknowns.

par_Nint

The number of global unknowns.

par_nghostint

The number of locally owned ghost unknowns.

subdomain2globalnumpy.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_proteus2petscnumpy.ndarray

Array with index permutations for mapping between proteus and petsc degrees of freedom.

classmethod create_ParMat_from_OperatorConstructor(operator)

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)

Saves to disk using a PETSc binary viewer.

proteus.LinearAlgebraTools.Vec(n)

Build a vector of length n (using numpy)

For example:

>>> Vec(3)
array([ 0.,  0.,  0.])

proteus.LinearAlgebraTools.Mat(m, n)

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)

Build a nr x nc sparse matrix from a dictionary representation

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

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)

Bases: object

Build a parallel matrix shell from CSR data structures.

Parameters

ghosted_csr_mat –

create(A)

mult(A, x, y)

class proteus.LinearAlgebraTools.OperatorShell

Bases: object

A base class for operator shells

create(A)

getSize()

Return the number of degrees of freedom for the operator.

class proteus.LinearAlgebraTools.ProductOperatorShell

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)

class proteus.LinearAlgebraTools.InvOperatorShell

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)

getSize()

Returns the size of InvOperatorShell.

Notes

This acts a virtual method and must be implemented for all inherited classes.

create_petsc_ksp_obj(petsc_option_prefix, matrix_operator, constant_null_space=False)

Create a PETSc4py KSP object.

Parameters

• petsc_option_prefix (str) – PETSc commandline option prefix.

• matrix_operator (mat) – PETSc matrix object for the ksp class.

• null_space (bool) – True if the KSP object has a constant null space.

Returns

ksp_obj

Return type

PETSc ksp

class proteus.LinearAlgebraTools.LSCInv_shell(Qv, B, Bt, F)

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)

Apply the LSC inverse operator

Parameters

• A (NULL) – A placeholder for internal function PETSc functions.

• x (p4pyPETSc.Vec) – Vector which LSC operator is being applied to.

Returns

y – Result of LSC acting on x.

Return type

p4pyPETSc.Vec

class proteus.LinearAlgebraTools.MatrixShell(A)

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)

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)

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)

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.SpInv_shell(A00, A11, A01, A10, constNullSpace)

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 (p4pyPETSc.Mat) – The A00 block of the global saddle point system.

• A01 (p4pyPETSc.Mat) – The A01 block of the global saddle point system.

• A10 (p4pyPETSc.Mat) – The A10 block of the global saddle point system.

• A11 (p4pyPETSc.Mat) – 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\) operator resembles a Laplcian matrix on the pressure. In cases where the global saddle point system uses pure Dirichlet boundary conditions, the \(S^{-1}\) operator has a constant null space. Since most saddle-point simulations of interest do not have pure Dirichlet conditions, the constNullSpace flag defaults to false. Having the null space set to false when the global problem uses pure Dirichlet boundary conditions will likely result in poor solver performance or failure.

apply(A, x, y)

Applies the \(S_{p}\) operator

Parameters

• A (None) – Dummy argument for PETSc interface

• x (p4pyPETSc.Vec) – Vector to which \(S\) is applied

Returns

y – Result of \(S^{-1}x\)

Return type

p4pyPETSc.Vec

class proteus.LinearAlgebraTools.TwoPhase_PCDInv_shell(Qp_visc, Qp_dens, Ap_rho, Np_rho, alpha=False, delta_t=0, num_chebyshev_its=0, strong_dirichlet_DOF=[], laplace_null_space=False, par_info=None)

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)

• strong_dirichlet_DOF (lst) – List of DOF with known, strongly enforced values.

• laplace_null_space (binary) – Indicates whether the pressure Laplace matrix has a null space or not.

• par_info (ParInfoClass) – Provides parallel info.

getSize()

Return the total number of DOF for the shell problem.

apply(A, x, y)

Applies the two-phase pressure-convection-diffusion Schur complement approximation.

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

Notes

When strong Dirichlet conditions are enforced on the pressure, the PCD operator is applied to the set of unknowns that do not have Dirichlet boundary conditions. At the end, the solution is then loaded into the original y-vector.

proteus.LinearAlgebraTools.l2Norm(x)

Compute the parallel \(l_2\) norm

proteus.LinearAlgebraTools.l1Norm(x)

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)

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)

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)

Compute the parallel weighted l_2 norm with weight h

proteus.LinearAlgebraTools.wl1Norm(x, h)

Compute the parallel weighted l_1 norm with weight h

proteus.LinearAlgebraTools.wlInfNorm(x, h)

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

proteus.LinearAlgebraTools.energyDot(x, y, A)

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

proteus.LinearAlgebraTools.energyNorm(x, A)

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

proteus.LinearAlgebraTools.l2NormAvg(x)

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

proteus.LinearAlgebraTools.rmsNorm(x)

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

proteus.LinearAlgebraTools.l2Norm_local(x)

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

class proteus.LinearAlgebraTools.WeightedNorm(shape, atol, rtol)

Bases: object

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

setWeight(y)

norm(y, type)