PetIBM 0.5.4
Toolbox and applications of the immersed-boundary method for distributed-memory architectures
Namespaces | Functions
Operator factories
PetIBM building blocks

Matrices & Operators used in Perot's framework (1993, 2002). More...

Collaboration diagram for Operator factories:

Namespaces

namespace  petibm::operators
 Collections of factory functions of operators.
 

Functions

PetscErrorCode petibm::operators::createR (const type::Mesh &mesh, Mat &R)
 Create a matrix of flux areas, $R$. More...
 
PetscErrorCode petibm::operators::createRInv (const type::Mesh &mesh, Mat &RInv)
 Create a matrix of inversed flux areas, $R^{-1}$. More...
 
PetscErrorCode petibm::operators::createMHead (const type::Mesh &mesh, Mat &MHead)
 Create a matrix of cell widths of velocity points, $\hat{M}$. More...
 
PetscErrorCode petibm::operators::createM (const type::Mesh &mesh, Mat &M)
 Create a matrix $M=\hat{M}R^{-1}$. More...
 
PetscErrorCode petibm::operators::createIdentity (const type::Mesh &mesh, Mat &I)
 Create an identity matrix $I$ for velocity fields. More...
 
PetscErrorCode petibm::operators::createGradient (const type::Mesh &mesh, Mat &G, const PetscBool &normalize=PETSC_TRUE)
 Create a gradient operator, $G$, for pressure field. More...
 
PetscErrorCode petibm::operators::createDivergence (const type::Mesh &mesh, const type::Boundary &bc, Mat &D, Mat &DCorrection, const PetscBool &normalize=PETSC_TRUE)
 Create a divergence operator, $D$, and corresponding boundary correction, $D_{bc}$, for velocity fields. More...
 
PetscErrorCode petibm::operators::createLaplacian (const type::Mesh &mesh, const type::Boundary &bc, Mat &L, Mat &LCorrection)
 Create a Laplacian operator, $L$, and boundary correction, $L_{bc}$ for velocity fields. More...
 
PetscErrorCode petibm::operators::createConvection (const type::Mesh &mesh, const type::Boundary &bc, Mat &H)
 Create a matrix-free Mat for convection operator, $H$. More...
 
PetscErrorCode petibm::operators::createBnHead (const Mat &Op, const PetscReal &dt, const PetscReal &coeff, const PetscInt &n, Mat &BnHead)
 Create non-normalized matrix of approximated $A^{-1}$, $B_n$. More...
 
PetscErrorCode petibm::operators::createBn (const Mat &Op, const Mat &R, const Mat &MHead, const PetscReal &dt, const PetscReal &coeff, const PetscInt &n, Mat &Bn)
 Create normalized matrix of approximated $A^{-1}$, $B_n$. More...
 
PetscErrorCode petibm::operators::createBn (const Mat &Op, const Mat &M, const PetscReal &dt, const PetscReal &coeff, const PetscInt &n, Mat &Bn)
 Create normalized matrix of approximated $A^{-1}$, $B_n$. More...
 
PetscErrorCode petibm::operators::createDelta (const type::Mesh &mesh, const type::Boundary &bc, const type::BodyPack &bodies, const delta::DeltaKernel &kernel, const PetscInt &kernelSize, Mat &Op)
 Create a Delta operator, $Delta$. More...
 

Detailed Description

Matrices & Operators used in Perot's framework (1993, 2002).

An operator is a matrix, which will be the type of PETSc Mat. It may be a sparse matrix or a matrix-free matrix. Operators in CFD are based on the discretization of domain and boundary conditions, i.e. mesh. Therefore, all operator factories require a Mesh instance as an input, and most of them require a Boundary instance as another input.

Given that PetIBM now only supports staggered grid, most of the operators only make sense to staggered Cartesian CFD methods. So the operator factories can not create the operators for arbitrary fields. For example, createLaplacian can only create Laplacian for velocity fields because in our CFD framework (Perot 1993), only the Laplacian of velocity fields is required. This factory function cannot create Laplacian operators for arbitrary fields.

For operators taking boundary conditions, we decouple the interior part and boundary corrections. For example, applying a divergence operator to velocity will become

\[(D + D_{bc})u\]

in which $D$ and $D_{bc}$ represent divergence operators of interior points and ghost points (i.e., boundary correction), respectively.

Boundary conditions in PetIBM are implemented with ghost point schemes. Most of boundary conditions can be represented with, for example,

\[ u_{i} = a_{0} u_{i-1} + a_{1} \]

where $u_{i}$ and $u_{i-1}$ represent a ghost point and its corresponding interior point, respectively. Hence when we apply ghost point schemes to operators, the $a_{0}$ part will get into the interior operator (e.g. $D$), and the $a_{1}$ part will be decoupled into the boundary correction matrix, i.e., $D_{bc}$.

Here is an example code for creating a divergence operator and its boundary correction operator, and calculating the divergence of velocity.

PetscErrorCode ierr;
Mat D, Dbc;
Vec div;
petibm::type::Solution solution;
petibm::type::Mesh mesh;
petibm::type::Boundary boundary;
petibm::type::Solution solution;
// create mesh with petibm::mesh::createMesh
// create boundary with petibm::boundary::createBoundary
// create solution with petibm::solution::createSolution, and have some
values
// create vector div with correct memory allocation
// create divergence operator
ierr = petibm::createDivergence(mesh, boundary, D, Dbc); CHKERRQ(ierr);
// calculate divergence of velocity, result will be in div
ierr = MatMult(D, solution->UGlobal, div); CHKERRQ(ierr);
ierr = MatMultAdd(Dbc, solution->UGlobal, div, div); CHKERRQ(ierr);
petibm::type::Boundary
std::shared_ptr< boundary::BoundaryBase > Boundary
Type definition of petibm::type::Boundary.
Definition: boundary.h:175
petibm::type::Mesh
std::shared_ptr< mesh::MeshBase > Mesh
Type definition of Mesh.
Definition: mesh.h:348
petibm::operators::createDivergence
PetscErrorCode createDivergence(const type::Mesh &mesh, const type::Boundary &bc, Mat &D, Mat &DCorrection, const PetscBool &normalize=PETSC_TRUE)
Create a divergence operator, , and corresponding boundary correction, , for velocity fields.
Definition: createdivergence.cpp:103
petibm::type::Solution
std::shared_ptr< solution::SolutionBase > Solution
Type definition of solution object.
Definition: solution.h:210

Function Documentation

◆ createBn() [1/2]

PetscErrorCode petibm::operators::createBn ( const Mat &  Op,
const Mat &  M,
const PetscReal &  dt,
const PetscReal &  coeff,
const PetscInt &  n,
Mat &  Bn 
)

Create normalized matrix of approximated $A^{-1}$, $B_n$.

Parameters
Op[in] Implicit operator $A$.
M[in] The $M$ matrix as defined in $A = M/dt - coeff * Op$.
dt[in] Time-step size.
coeff[in] Time-scheme coefficient for the implicit operator.
n[in] an integer; the order of Bn.
Bn[out] Matrix $B_n$.

Similar to createBnHead, except that A is now defined as

\[ A = \frac{M}{dt} - coeff \times Op \]

Definition at line 131 of file createbn.cpp.

◆ createBn() [2/2]

PetscErrorCode petibm::operators::createBn ( const Mat &  Op,
const Mat &  R,
const Mat &  MHead,
const PetscReal &  dt,
const PetscReal &  coeff,
const PetscInt &  n,
Mat &  Bn 
)

Create normalized matrix of approximated $A^{-1}$, $B_n$.

Parameters
Op[in] Implicit operator $A$.
R[in] Matrix to define $M = \hat{M}R^{-1}$.
MHead[in] Matrix to define $M = \hat{M}R^{-1}$.
dt[in] Time-step size.
coeff[in] Time-scheme coefficient for the implicit operator.
n[in] Truncation index of the Taylor series.
Bn[out] Matrix, $B_n$.

Similar to createBnHead, except that A is now defined as

\[ A = \frac{M}{dt} - coeff \times \hat{M}OpR \]

where

\[ M = \hat{M}R^{-1} \]

Definition at line 98 of file createbn.cpp.

◆ createBnHead()

PetscErrorCode petibm::operators::createBnHead ( const Mat &  Op,
const PetscReal &  dt,
const PetscReal &  coeff,
const PetscInt &  n,
Mat &  BnHead 
)

Create non-normalized matrix of approximated $A^{-1}$, $B_n$.

Parameters
Op[in] Implicit operator $A$.
dt[in] Time-step size.
coeff[in] Time-scheme coefficient for the implicit operator.
n[in] Truncation index of the Taylor series.
BnHead[out] Operator $B_n$.

$A$ is defined as *

\[A=\frac{I}{dt} - coeff \times Op\]

where $I$ is an identity matrix, $Op$ is the summation of all implicit operators, and $coeff$ is the coefficient from temporal integration. For example, for the case of Adams-Bashforth for convective terms and Crank-Nicolson for diffusion term,

\[A = \frac{I}{dt} - 0.5 \times L\]

If the implicit part is more complicated, there may be not only one coefficient for each implicit operator. Then, in order to use this function, we may need to incorporate coefficients and all implicit operators into a big implicit operator before calling this function. And then set $coeff$ as $1$. For example, if an $A$ is defined as

\[ A = \frac{I}{dt} - c_{1}\times Op_{1} - c_{2}\times Op_{2}\]

then we need to first define

\[Op = c_1\times Op_1 + c_2\times Op_2\]

so that $A$ can have the form of

\[ A = \frac{I}{dt} - 1.0\times Op \]

.

Definition at line 19 of file createbn.cpp.

◆ createConvection()

PetscErrorCode petibm::operators::createConvection ( const type::Mesh &  mesh,
const type::Boundary &  bc,
Mat &  H 
)

Create a matrix-free Mat for convection operator, $H$.

Parameters
mesh[in] Structured Cartesian mesh object.
bc[in] Data object with boundary conditions.
H[out] Operator $H$.

Though this matrix-free operator is a PETSc Mat, not all PETSc Mat functions can be applied to it because we only register MatMult and MatDestroy functions in PetIBM.

If $u$ is the composite velocity vector, then the multiplication $Hu$ will return the convection terms at velocity points.

Definition at line 372 of file createconvection.cpp.

◆ createDelta()

PetscErrorCode petibm::operators::createDelta ( const type::Mesh &  mesh,
const type::Boundary &  bc,
const type::BodyPack &  bodies,
const delta::DeltaKernel &  kernel,
const PetscInt &  kernelSize,
Mat &  Op 
)

Create a Delta operator, $Delta$.

Parameters
mesh[in] Structured Cartesian mesh object.
bc[in] Data object with boundary conditions.
bodies[in] Data object with the immersed boundaries.
kernel[in] Regularized delta kernel to use.
kernelSize[in] Size of the kernel.
Op[out] Matrix $Delta$.

In $Delta$, rows represent Lagrangian points, while columns represent velocity mesh points.

An entry in $Delta$, says $Delta_{(i,j)}$, is defined as

\[ Delta_{(i,j)} = d(x_j-\xi_i)d(y_j-\eta_i)d(z_j-\zeta_i) \]

where $(x_j,y_j,z_j)$ and $(\xi_i,\eta_i,\zeta_i)$ represent the coordinates of the $j$-th Eulerian grid point and the $i$-th Lagrangian point, respectively. $d(...)$ is discretized delta function following Roma et. al. (1999).

Definition at line 34 of file createdelta.cpp.

◆ createDivergence()

PetscErrorCode petibm::operators::createDivergence ( const type::Mesh &  mesh,
const type::Boundary &  bc,
Mat &  D,
Mat &  DCorrection,
const PetscBool &  normalize = PETSC_TRUE 
)

Create a divergence operator, $D$, and corresponding boundary correction, $D_{bc}$, for velocity fields.

Parameters
mesh[in] Structured Cartesian mesh object.
bc[in] Data object with the boundary conditions.
D[out] Divergence operator $D$.
DCorrection[out] Operator for boundary corrections, $D_{bc}$.
normalize[in] Use PETSC_TRUE to normalize operator (default is PETSC_TRUE).

PETSc matrix D should not be created before calling this function.

If normalize=PETSC_TRUE (which is default), then the divergence operator, $D$, will only contain 0 and 1.

The divergence operator calculates the divergence of velocity on pressure points. The complete divergence operator is $(D+D_{bc})$. However, the boundary correction, $D_{bc}$, represents the constant coefficients (it means the coefficients that have nothing to do with interior points, though it may still be time-dependent) in ghost point schemes, so it's safe to put the boundary correction to right-hand side in a linear system.

Definition at line 103 of file createdivergence.cpp.

◆ createGradient()

PetscErrorCode petibm::operators::createGradient ( const type::Mesh &  mesh,
Mat &  G,
const PetscBool &  normalize = PETSC_TRUE 
)

Create a gradient operator, $G$, for pressure field.

Parameters
mesh[in] Structured Cartesian mesh object.
G[out] $G$ matrix.
normalize[in] Use PETSC_TRUE to normalize the matrix (default is PETSC_TRUE).

Petsc matrix G should not be created before calling this function.

If normalize=PETSC_TRUE (which is default), then the gradient operator, $G$, will only contain 0 and 1.

This gradient operator doesn't require boundary conditions because we only need gradient of pressure at interior velocity points (see Perot 1993). And it can only be applied to pressure field. For example:

ierr = petibm::operators::createGradient(mesh, G); CHKERRQ(ierr);
ierr = MatMult(G, solution->pGlobal, grad); CHKERRQ(ierr);
petibm::operators::createGradient
PetscErrorCode createGradient(const type::Mesh &mesh, Mat &G, const PetscBool &normalize=PETSC_TRUE)
Create a gradient operator, , for pressure field.
Definition: creategradient.cpp:36

Definition at line 36 of file creategradient.cpp.

◆ createIdentity()

PetscErrorCode petibm::operators::createIdentity ( const type::Mesh &  mesh,
Mat &  I 
)

Create an identity matrix $I$ for velocity fields.

Parameters
mesh[in] Structured Cartesian mesh object.
I[out] $I$ matrix.

PETSc matrix I should not be created before calling this function.

Definition at line 210 of file creatediagmatrix.cpp.

◆ createLaplacian()

PetscErrorCode petibm::operators::createLaplacian ( const type::Mesh &  mesh,
const type::Boundary &  bc,
Mat &  L,
Mat &  LCorrection 
)

Create a Laplacian operator, $L$, and boundary correction, $L_{bc}$ for velocity fields.

Parameters
mesh[in] Structured Cartesian mesh object.
bc[in] Data object for the boundary conditions.
L[out] Laplacian operator $L$.
LCorrection[out] Operator for boundary corrections, $L_{bc}$.

PETSc matrix L should not be created before calling this function.

The complete Laplacian is $(L+L_{bc})$. The boundary correction represents the constant coefficients (the coefficients that have nothing to do with interior velocity points) in ghost point schemes. So it's often seen the boundary correction $L_{bc}$ on right-hand side of a system, while $L$ on left-hand side.

Definition at line 170 of file createlaplacian.cpp.

◆ createM()

PetscErrorCode petibm::operators::createM ( const type::Mesh &  mesh,
Mat &  M 
)

Create a matrix $M=\hat{M}R^{-1}$.

Parameters
mesh[in] Structured Cartesian mesh object.
M[out] $M$ matrix.

PETSc matrix M should not be created before calling this function.

In some CFD schemes, the combination of $\hat{M}$ and $R^{-1}$ is used frequently. So this function creates a matrix $M=\hat{M}R^{-1}$ for convenience.

Definition at line 180 of file creatediagmatrix.cpp.

◆ createMHead()

PetscErrorCode petibm::operators::createMHead ( const type::Mesh &  mesh,
Mat &  MHead 
)

Create a matrix of cell widths of velocity points, $\hat{M}$.

Parameters
mesh[in] Structured Cartesian mesh object.
MHead[out] $\hat{M}$ matrix.

PETSc matrix MHead should not be created before calling this function.

Operator $\hat{M}$ represents the cell width of velocity points. The cell width of a $u$-velocity point with index $(i, j, k)$, for example, is $dx_{i,j,k}$, and the cell width of $v$-velocity at $(i, j, k)$ is $dy_{i,j,k}$.

The multiplication of $\hat{M}$ and $R$ returns the cell volumes of velocity points, i.e., $\hat{M}R=cell\ volumes$.

Definition at line 150 of file creatediagmatrix.cpp.

◆ createR()

PetscErrorCode petibm::operators::createR ( const type::Mesh &  mesh,
Mat &  R 
)

Create a matrix of flux areas, $R$.

Parameters
mesh[in] Structured Cartesian mesh object.
R[out] $R$ matrix.

PETSc matrix R should not be created before calling this function.

Operator $R$ represents the flux areas of velocity points. For example, if $q$ is vector of velocity flux, and $u$ is velocity vector, then $Rq=u$.

See also
petibm::operators::createRInv

Definition at line 90 of file creatediagmatrix.cpp.

◆ createRInv()

PetscErrorCode petibm::operators::createRInv ( const type::Mesh &  mesh,
Mat &  RInv 
)

Create a matrix of inversed flux areas, $R^{-1}$.

Parameters
mesh[in] an instance of petibm::type::Mesh.
RInv[out] a PETSc Mat; returned matrix $R^{-1}$.

PETSc matrix RInv should not be created before calling this function.

Operator $R^{-1}$ is the inverse of $R$. And it has relationship $R^{-1}u=q$, in which $u$ and $q$ are as defined in petibm::operators::createR.

See also
petibm::operators::createR

Definition at line 120 of file creatediagmatrix.cpp.