This tutorial describes how to use PETSc for domain decomposition to parallelize (with MPI) a solver for the simple advection equation (1.1).

Equation 1.1

In equation (1.1), D is a scalar density and V = { Vx,Vy,Vz } is a velocity field. We assume only Vx is non-zero and constant. A simple time-explicit scheme for solving equation (1.1) is given by (1.2)

Equation 1.2

where, for example, the “upwind flux” gives (1.3)

Equation 1.3

Thus, the numerical scheme couples nearest neighbors, which means that data must be exchanged between MPI tasks during the evolution. Although the equation solved in this tutorial example is simple, the strategy outlined here can be easily generalized to more complex problems using more advanced numerical methods. The example code provided (with the driver Advection3D.F90) performs three main tasks; (1) creates a distributed logically Cartesian grid/mesh (in MeshModule.F90; Section 4), (2) creates physical fields associated with each grid cell (in FieldsModule.F90; Section 5), and (3) evolves the simple advection equation on the distributed grid (in Advection3D.F90 and AdvectionModule.F90). The main data structure used by the program (UniverseType, defined in UniverseModule.F90) holds one data structure for the mesh (MeshType, defined in MeshModule.F90) and one data structure for the physical fields (FieldsType, defined in FieldsModule.F90). The actual time stepping loop is in the main driver (Advection3D.F90).

Note: this tutorial is intended to supplement — not replace — the existing PETSc documentation, available at http://www.mcs.anl.gov/petsc/documentation. In addition, also see the useful code examples provided with the PETSc distribution

• \$PETSC_DIR/source/src/dm/examples/tutorials/
• \$PETSC_DIR/source/src/ksp/ksp/examples/tutorials/
• \$PETSC_DIR/source/src/mat/examples/tutorials/
• \$PETSC_DIR/source/src/snes/examples/tutorials/
• \$PETSC_DIR/source/src/sys/examples/tutorials/
• \$PETSC_DIR/source/src/ts/examples/tutorials/
• \$PETSC_DIR/source/src/vec/vec/examples/tutorials/

In the following section we describe the individual modules of a Fortran code that use PETSc for domain decomposition. In each section we provide links to PETSc documentation for all the PETSc routines used by the code. The user is encouraged to download, compile, and run the provided code.

Please direct any questions or comments to help@nccs.gov.

### Mesh and Data Structures

The call to CreateUniverse creates the main data structure UniverseType. PETSc (and MPI) are also initialized in CreateUniverse. It should be the first call in any program using the data structure defined by UniverseType. Notable input parameters are the number of mesh cells in each spatial dimension, the inner and outer boundaries of the computational domain, and the type of boundary conditions used. BoundaryConditions are of type DMDABoundaryType, and PETSc supports DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_GHOST, and DMDA_BOUNDARY_PERIODIC.

```  !  Program solves dD/dt+Vx*dD/dx+Vy*dD/dy+Vz*dD/dz=0 on periodic unit domain
!    with constant advection velocity (V1, V2, V3) = (1, 0, 0),
!    and initial condition D = Sin(2 * Pi * X1)

!  Create Mesh and Data Structures to Hold Physical Fields:

nZones(1:3) = 256

CALL CreateUniverse( &
U, Name = 'Advection3D',          &
nZones                            &
= nZones,                       &
InnerBoundaries                   &
= (/ Zero, Zero, Zero /),       &
OuterBoundaries                   &
= (/ One,  One,  One  /),       &
BoundaryConditions                &
= (/ DMDA_BOUNDARY_PERIODIC,    &
DMDA_BOUNDARY_PERIODIC,    &
DMDA_BOUNDARY_PERIODIC /))
```

The call to DestroyUniverse frees up memory allocated by CreateUniverse and finalizes PETSc (and MPI) to end cleanly.

```  ! Destroy Mesh and Data Structures that Hold Physical Fields:

CALL DestroyUniverse(U)
```

## UniverseModule.F90

The following PETSc routines are called
PetscInitialize
PetscFinalize
PetscGetTime

### UniverseType

UniverseType contains Mesh (MeshType) and Fields (FieldsType), and additional rudimentary variables associated with parallel execution (e.g., the integer Communicator).

```  TYPE, PUBLIC :: UniverseType
CHARACTER(80) :: &
Name
PetscInt :: &
Communicator, &
MPI_Rank, &
N_MPI_Ranks
PetscLogDouble :: &
WallTimeInitialize, &
WallTimeFinalize
TYPE(MeshType), POINTER :: &
Mesh   !
TYPE(FieldsType), POINTER :: &
Fields !
END TYPE UniverseType
```

### CreateUniverse

CreateUniverse initializes PETSc (and MPI) with the call to PetscInitialize. It also initializes the rudimentary variables and creates the data structures MeshType and FieldsType with the calls to CreateMesh and CreateFields.

```  SUBROUTINE CreateUniverse(U, Name, &
nZones, InnerBoundaries, OuterBoundaries, BoundaryConditions)

TYPE(UniverseType), POINTER :: &
U
CHARACTER(len=*), INTENT(in) :: &
Name
PetscInt, DIMENSION(:), INTENT(in) :: &
nZones
PetscReal, DIMENSION(:), INTENT(in) :: &
InnerBoundaries, &
OuterBoundaries
DMDABoundaryType, DIMENSION(:), INTENT(in) :: &
BoundaryConditions

PetscInt :: &
Error

CALL PetscInitialize(PETSC_NULL_CHARACTER, Error)

ALLOCATE(U)
U % Communicator = MPI_COMM_WORLD
CALL MPI_Comm_Rank(U % Communicator, U % MPI_Rank,    Error)
CALL MPI_Comm_Size(U % Communicator, U % N_MPI_Ranks, Error)

U % Name = Name

CALL PetscGetTime(U % WallTimeInitialize, Error)

IF(U % MPI_RANK == 0)THEN
PRINT*
PRINT*, "INFO: Creating Universe ", TRIM(U % Name)
PRINT*
PRINT*, "  INFO: Number of MPI Ranks = ", U % N_MPI_Ranks
END IF

ALLOCATE(U % Mesh)
CALL CreateMesh(U % Mesh, U % Communicator, U % MPI_Rank, &
nZones, InnerBoundaries, OuterBoundaries, BoundaryConditions)

ALLOCATE(U % Fields)
CALL CreateFields(U % Fields, U % Mesh, U % Communicator, U % MPI_Rank)

END SUBROUTINE CreateUniverse
```

The call to PetscGetTime is used to time the program (see DestroyUniverse).

### DestroyUniverse

DestroyUniverse frees up memory allocated by CreateUniverse. In particular, Fields and Mesh are deallocated by calls to DestroyFields and DestroyMesh. DestroyUniverse also finalizes PETSc (and MPI) to exit cleanly.

``` SUBROUTINE DestroyUniverse(U)

TYPE(UniverseType), POINTER :: &
U

PetscInt :: &
Error

IF(U % MPI_Rank == 0)THEN
PRINT*
PRINT*, "INFO: Destroying Universe ", TRIM(U % Name)
END IF

CALL DestroyFields(U % Fields)
DEALLOCATE(U % Fields)

CALL DestroyMesh(U % Mesh)
DEALLOCATE(U % Mesh)

CALL PetscGetTime(U % WallTimeFinalize, Error)

IF(U % MPI_Rank == 0)THEN
PRINT*
PRINT*, "  INFO: ", U % N_MPI_Ranks, " MPI Ranks Finished in ", &
U % WallTimeFinalize - U % WallTimeInitialize, " Seconds"
PRINT*
END IF

CALL PetscFinalize(Error)

DEALLOCATE(U)

END SUBROUTINE DestroyUniverse
```

## MeshModule.F90

X1,X2,X3 are integers used to access the mesh coordinates (see e.g., CreateMesh3D for usage).

```  ! Coordinate Names

PetscInt, PUBLIC, PARAMETER :: &
X1 = 0, &
X2 = 1, &
X3 = 2
```

### MeshType

MeshType contains basic information about the mesh; e.g., dimensionality (nDimensions), number of zones in each spatial dimension (nZones), the coordinates of the inner and outer boundaries (InnerBoundaries and OuterBoundaries), the type of boundary condition used in each dimension (BoundaryType). StencilWidth is the width of the computational stencil and sets the depth of the ghost zone layer to be communicated. (StencilWidth=1 in the current example, but can be larger for higher-order methods.) StencilType can be of type DMDA_STENCIL_BOX or DMDA_STENCIL_STAR to accommodate different computational stencils (see note below). The main data structures in MeshType are MeshDA and CoordinateDA, both are DM objects that are used to manage data for a structured grid in 1, 2, or 3 dimensions. MeshType also contains Positions (PositionsType), which holds mesh coordinate values.

```  TYPE, PUBLIC :: MeshType
PetscInt :: &
nDimensions, &
StencilWidth
PetscInt, DIMENSION(:), ALLOCATABLE :: &
nZones
PetscReal, DIMENSION(:), ALLOCATABLE :: &
InnerBoundaries, &
OuterBoundaries
DMDAStencilType :: &
StencilType
DMDABoundaryType, DIMENSION(:), ALLOCATABLE :: &
BoundaryType
DM :: &
MeshDA, &
CoordinateDA
TYPE(PositionsType), POINTER :: &
Positions
END TYPE MeshType
```

PETSc supports two types of computational stencils with DMDAStencilType. It determines what ghost zones will be available to a process (see Figure). With DMDA_STENCIL_STAR only the blue cells will be available. With DMDA_STENCIL_BOX both the red and blue cells are available. StencilType is set to DMDA_STENCIL_STAR in CreateMesh.

### PositionsType

PositionsType contains distributed vectors (of type Vec) to hold position coordinates. InnerEdge, Center, and OuterEdge hold the inner, center, and outer coordinates of a computational mesh cell. Vectors with the “Global” suffix (e.g., InnerEdgeGlobal) are created with DMCreateGlobalVector and do not contain ghost zones. Vectors with the “Local” suffix (e.g., InnerEdgeLocal) are created with DMCreateLocalVector and contains ghost zones (see CreatePositions in Section 4 e). The integers iBX(1:3) (“i Begin X”), iWX(1:3) (“i Width X”), and iEX(1:3) (“i End X”) are used to access specific cells on a mesh portion owned by a particular process. We have also defined the corresponding iBXG(1:3), iWXG(1:3), and iEXG(1:3) for the ghost domains. iBX and iWX, and iBXG and iWXG are set with calls to DMDAGetCorners and DMDAGetGhostedCorners (see Section 4e).

```TYPE :: PositionsType
PetscInt, DIMENSION(3) :: &
iBX,  & ! Inner Indices of Local Domain
iWX,  & ! Width of Local Domain
iEX,  & ! Outer Indices of Local Domain
iBXG, & ! Inner Indices of Local Ghosted Domain
iWXG, & ! Width of Local Ghosted Domain
iEXG    ! Outer Indices of Local Ghosted Domain
Vec :: &
InnerEdgeGlobal, &
CenterGlobal,    &
OuterEdgeGlobal, &
InnerEdgeLocal,  &
CenterLocal,     &
OuterEdgeLocal
END TYPE PositionsType
```

### CreateMesh

CreateMesh sets basic information about the mesh, and, depending on the dimensionality of the problem, calls CreateMesh1D, CreateMesh2D, or CreateMesh3D, which creates the data structures MeshDA and CoordinateDA and sets the position coordinates. After the global position coordinates are initialized (e.g., InnerEdgeGlobal), the calls DMGlobalToLocalBegin and DMGlobalToLocalEnd populate ghost values in local position coordinates (e.g., InnerEdgeLocal).

```  SUBROUTINE CreateMesh(M, Comm, Rank, &
nZones, InnerBoundaries, OuterBoundaries, BoundaryConditions)

TYPE(MeshType), POINTER :: &
M
PetscInt, INTENT(in) :: &
Comm, &
Rank
PetscInt, DIMENSION(:), INTENT(in) :: &
nZones
PetscReal, DIMENSION(:), INTENT(in) :: &
InnerBoundaries, &
OuterBoundaries
DMDABoundaryType, DIMENSION(:), INTENT(in) :: &
BoundaryConditions

PetscInt :: &
iDim, &
Error

M % nDimensions = SIZE(nZones)

ALLOCATE( M % nZones         (M % nDimensions) )
ALLOCATE( M % InnerBoundaries(M % nDimensions) )
ALLOCATE( M % OuterBoundaries(M % nDimensions) )
ALLOCATE( M % BoundaryType   (M % nDimensions) )

DO iDim = 1, M % nDimensions
M % nZones(iDim)          = nZones(iDim)
M % InnerBoundaries(iDim) = InnerBoundaries(iDim)
M % OuterBoundaries(iDim) = OuterBoundaries(iDim)
M % BoundaryType(iDim)    = BoundaryConditions(iDim)
END DO

M % StencilWidth = 1
M % StencilType  = DMDA_STENCIL_STAR

IF(Rank == 0)THEN
PRINT*
PRINT*, "  INFO: Creating Mesh"
PRINT*, "    Dimensionality = ", M % nDimensions
PRINT*, "    Stencil Width  = ", M % StencilWidth
DO iDim = 1, M % nDimensions
PRINT*, "      iDim, nZones = ", iDim, M % nZones(iDim)
PRINT*, "        Inner Boundary   = ", M % InnerBoundaries(iDim)
PRINT*, "        Outer Boundary   = ", M % OuterBoundaries(iDim)
END DO
END IF

SELECT CASE (M % nDimensions)
CASE (1)

CALL CreateMesh1D(M, Comm)

CASE (2)

CALL CreateMesh2D(M, Comm)

CASE (3)

CALL CreateMesh3D(M, Comm)

END SELECT

!  Populate Ghost Zones with Coordinates:

CALL DMGlobalToLocalBegin( &
M % CoordinateDA, &
M % Positions % InnerEdgeGlobal, INSERT_VALUES, &
M % Positions % InnerEdgeLocal, Error)
CALL DMGlobalToLocalEnd( &
M % CoordinateDA, &
M % Positions % InnerEdgeGlobal, INSERT_VALUES, &
M % Positions % InnerEdgeLocal, Error)

CALL DMGlobalToLocalBegin( &
M % CoordinateDA, &
M % Positions % OuterEdgeGlobal, INSERT_VALUES, &
M % Positions % OuterEdgeLocal, Error)
CALL DMGlobalToLocalEnd( &
M % CoordinateDA, &
M % Positions % OuterEdgeGlobal, INSERT_VALUES, &
M % Positions % OuterEdgeLocal, Error)

CALL DMGlobalToLocalBegin( &
M % CoordinateDA, &
M % Positions % CenterGlobal, INSERT_VALUES, &
M % Positions % CenterLocal, Error)
CALL DMGlobalToLocalEnd( &
M % CoordinateDA, &
M % Positions % CenterGlobal, INSERT_VALUES, &
M % Positions % CenterLocal, Error)

END SUBROUTINE CreateMesh
```

### CreateMesh3D

CreateMesh3D creates a three-dimensional mesh. The call to the PETSc routine DMDACreate3D creates the object MeshDA (of type DM) to manage the communication of three-dimensional data that is distributed across processors (see the provided code for corresponding 1D and 2D versions). The call to DMDAGetCoordinateDA initializes CoordinateDA (also of type DM) that is used to scatter between global and local coordinates and to access the coordinate vectors defined in PositionsType. The call to CreatePositions creates/allocates the position coordinate vectors and initializes indices used to access specific cells in the mesh owned by a particular process. To fill in values in a global vector (e.g., InnerEdgeLocal) we must call DMDAVecGetGetArrayF90, which returns a multi-dimensional array (e.g., InnerEdge) that shares data with the underlying vector. In the 3D case, InnerEdge is a rank-four array, where the first index denotes the coordinate (i.e., X1, X2, or X3) and the following three indices denote the logical position in the mesh. After the position coordinates have been set and access to InnerEdge is no longer needed, one must call DMDAVecResoreArrayF90.

``` SUBROUTINE CreateMesh3D(M, Comm)

TYPE(MeshType), POINTER :: &
M
PetscInt, INTENT(in) :: &
Comm

PetscInt :: &
Error, &
iX1, iX2, iX3
PetscReal, DIMENSION(3) :: &
dX
PetscReal, DIMENSION(:,:,:,:), POINTER :: &
InnerEdge, &
OuterEdge, &
Center

CALL DMDACreate3D( &
Comm, &
M % BoundaryType(1), M % BoundaryType(2), M % BoundaryType(3), &
M % StencilType, &
M % nZones(1), M % nZones(2), M % nZones(3), &
PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, &
1, M % StencilWidth, &
PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, &
M % MeshDA, Error)

CALL DMDAGetCoordinateDA( &
M % MeshDA, M % CoordinateDA, Error)

!
!  Create vectors to hold coordinate values:
!

ALLOCATE(M % Positions)
CALL CreatePositions(M % Positions, M % nDimensions, M % CoordinateDA)

!
!  Fill in coordinate values:
!

CALL DMDAVecGetArrayF90( &
M % CoordinateDA, M % Positions % InnerEdgeGlobal, InnerEdge, Error)
CALL DMDAVecGetArrayF90( &
M % CoordinateDA, M % Positions % OuterEdgeGlobal, OuterEdge, Error)
CALL DMDAVecGetArrayF90( &
M % CoordinateDA, M % Positions % CenterGlobal,    Center,    Error)

!  Equidistant mesh for simplicity
dX = ( M % OuterBoundaries - M % InnerBoundaries ) / REAL( M % nZones )

DO iX3 = M % Positions % iBX(3), M % Positions % iEX(3)
DO iX2 = M % Positions % iBX(2), M % Positions % iEX(2)
DO iX1 = M % Positions % iBX(1), M % Positions % iEX(1)

InnerEdge(X1, iX1, iX2, iX3) &
= M % InnerBoundaries(1) + dX(1) * iX1
OuterEdge(X1, iX1, iX2, iX3) &
= InnerEdge(X1, iX1, iX2, iX3) + dX(1)
Center   (X1, iX1, iX2, iX3) &
= Half * ( InnerEdge(X1, iX1, iX2, iX3) &
+ OuterEdge(X1, iX1, iX2, iX3) )

InnerEdge(X2, iX1, iX2, iX3) &
= M % InnerBoundaries(2) + dX(2) * iX2
OuterEdge(X2, iX1, iX2, iX3) &
= InnerEdge(X2, iX1, iX2, iX3) + dX(2)
Center   (X2, iX1, iX2, iX3) &
= Half * ( InnerEdge(X2, iX1, iX2, iX3) &
+ OuterEdge(X2, iX1, iX2, iX3) )

InnerEdge(X3, iX1, iX2, iX3) &
= M % InnerBoundaries(3) + dX(3) * iX3
OuterEdge(X3, iX1, iX2, iX3) &
= InnerEdge(X3, iX1, iX2, iX3) + dX(3)
Center   (X3, iX1, iX2, iX3) &
= Half * ( InnerEdge(X3, iX1, iX2, iX3) &
+ OuterEdge(X3, iX1, iX2, iX3) )

END DO
END DO
END DO

CALL DMDAVecRestoreArrayF90( &
M % CoordinateDA, M % Positions % InnerEdgeGlobal, InnerEdge, Error)
CALL DMDAVecRestoreArrayF90( &
M % CoordinateDA, M % Positions % OuterEdgeGlobal, OuterEdge, Error)
CALL DMDAVecRestoreArrayF90( &
M % CoordinateDA, M % Positions % CenterGlobal,    Center,    Error)

END SUBROUTINE CreateMesh3D
```

### CreatePositions

CreatePositions creates global and local vectors containing position coordinates (e.g., InnerEdgeLocal and InnerEdgeGlobal). Vectors created with DMCreateGlobalVector do not have ghost points accessible. Vectors created with DMCreateLocalVector have accessible ghost points. Calls to DMDAGetCorners and DMDAGetGhostedCorners returns the indices of the lower left corner of the local region and the extent (width) of the local region.

``` SUBROUTINE CreatePositions(Positions, nDimensions, CoordinateDA)

TYPE(PositionsType), POINTER :: &
Positions
PetscInt :: &
nDimensions
DM :: &
CoordinateDA

PetscInt :: &
Error, &
iDim

! Global (Non-Ghosted) Position Vectors:

CALL DMCreateGlobalVector( &
CoordinateDA,                Positions % InnerEdgeGlobal, Error)
CALL VecDuplicate( &
Positions % InnerEdgeGlobal, Positions % CenterGlobal,    Error)
CALL VecDuplicate( &
Positions % InnerEdgeGlobal, Positions % OuterEdgeGlobal, Error)

Positions % iBX(:) = 0
Positions % iWX(:) = 1
Positions % iEX(:) = 0

! Get Indices for Non-Ghosted Vectors:

SELECT CASE (nDimensions)
CASE (1)
CALL DMDAGetCorners( &
CoordinateDA, &
Positions % iBX(1), PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, &
Positions % iWX(1), PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, &
Error)

CASE (2)
CALL DMDAGetCorners( &
CoordinateDA, &
Positions % iBX(1), Positions % iBX(2), PETSC_NULL_INTEGER, &
Positions % iWX(1), Positions % iWX(2), PETSC_NULL_INTEGER, &
Error)

CASE(3)
CALL DMDAGetCorners( &
CoordinateDA, &
Positions % iBX(1), Positions % iBX(2), Positions % iBX(3), &
Positions % iWX(1), Positions % iWX(2), Positions % iWX(3), &
Error)

END SELECT

DO iDim = 1, nDimensions
Positions % iEX(iDim) &
= Positions % iBX(iDim) &
+ Positions % iWX(iDim) - 1
END DO

!  Local (Ghosted) Position Vectors:

CALL DMCreateLocalVector( &
CoordinateDA,               Positions % InnerEdgeLocal, Error)
CALL VecDuplicate( &
Positions % InnerEdgeLocal, Positions % CenterLocal,    Error)
CALL VecDuplicate( &
Positions % InnerEdgeLocal, Positions % OuterEdgeLocal, Error)

Positions % iBXG(:) = 0
Positions % iWXG(:) = 1
Positions % iEXG(:) = 0

!  Get Indices for Ghosted Vectors:

SELECT CASE (nDimensions)
CASE (1)
CALL DMDAGetGhostCorners( &
CoordinateDA, &
Positions % iBXG(1), PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, &
Positions % iWXG(1), PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, &
Error)

CASE (2)
CALL DMDAGetGhostCorners( &
CoordinateDA, &
Positions % iBXG(1), Positions % iBXG(2), PETSC_NULL_INTEGER, &
Positions % iWXG(1), Positions % iWXG(2), PETSC_NULL_INTEGER, &
Error)

CASE (3)
CALL DMDAGetGhostCorners( &
CoordinateDA, &
Positions % iBXG(1), Positions % iBXG(2), Positions % iBXG(3), &
Positions % iWXG(1), Positions % iWXG(2), Positions % iWXG(3), &
Error)

END SELECT

DO iDim = 1, nDimensions
Positions % iEXG(iDim) &
= Positions % iBXG(iDim) &
+ Positions % iWXG(iDim) - 1
END DO

END SUBROUTINE CreatePositions
```

### DestroyMesh

DestroyMesh frees memory allocated by CreateMesh

```  SUBROUTINE DestroyMesh(M)

TYPE(MeshType), POINTER :: &
M

PetscInt :: &
Error

DEALLOCATE( M % nZones )
DEALLOCATE( M % InnerBoundaries )
DEALLOCATE( M % OuterBoundaries )
DEALLOCATE( M % BoundaryType )

CALL DestroyPositions(M % Positions)
DEALLOCATE(M % Positions)

CALL DMDestroy(M % MeshDA, Error)

END SUBROUTINE DestroyMesh
```

### DestroyPositions

DestroyPositions destroys the vectors created in CreatePositions

```  SUBROUTINE DestroyPositions(Positions)

TYPE(PositionsType), POINTER :: &
Positions

PetscInt :: &
Error

CALL VecDestroy(Positions % InnerEdgeGlobal, Error)
CALL VecDestroy(Positions % CenterGlobal,    Error)
CALL VecDestroy(Positions % OuterEdgeGlobal, Error)

CALL VecDestroy(Positions % InnerEdgeLocal, Error)
CALL VecDestroy(Positions % CenterLocal,    Error)
CALL VecDestroy(Positions % OuterEdgeLocal, Error)

END SUBROUTINE DestroyPositions
```

## FieldsModule.F90

The following PETSc routines are called
DMDACreate3D
DMCreateGlobalVector
DMCreateLocalVector
VecDestroy
DMDestroy

D, V1, V2, V3 are integers used to access physical variables/fields (here representing the density and three components of velocity)

```  ! Field Names

PetscInt, PUBLIC, PARAMETER :: &
D        = 0, &
V1       = 1, &
V2       = 2, &
V3       = 3, &
N_FIELDS = 4
```

### FieldsType

FieldsType contains vectors to hold the physical fields: FieldsGlobal (non-ghosted) and FieldsLocal (ghosted). FieldsDA (of type DM) is the object used to manage data for a structured grid (e.g., exchange of data between processes). The integer DegreesOfFreedom holds the number of physical fields. It is used to create the FieldsDA object.

```  TYPE, PUBLIC :: FieldsType
PetscInt :: &
DegreesOfFreedom
Vec :: &
FieldsGlobal, &
FieldsLocal
DM :: &
FieldsDA
END TYPE FieldsType
```

### CreateFields

CreateFields initializes and creates data structures defined in FieldsType. Depending on the dimensionality of the problem, it creates the FieldsDA object by calling either DMDACreate1D, DMDACreate2D, or DMDACreate3D. It also creates the global and local vectors to hold the physical fields. Note that DegreesOfFreedom > 1 in this case, which has consequences for the shape of the array returned by DMDAVecGetArrayF90 (see initialization of FieldsGlobal in Advection3D.F90). For DegreesOfFreedom > 1 the first index in the array returned by DMDAVecGetArrayF90 is used to access specific fields. The preceding indices refer to the field’s logical position in the mesh.

```  SUBROUTINE CreateFields(F, M, Comm, Rank)

TYPE(FieldsType), POINTER :: &
F
TYPE(MeshType), POINTER :: &
M
PetscInt, INTENT(in) :: &
Comm, &
Rank

PetscInt :: &
Error

F % DegreesOfFreedom = N_FIELDS

IF(Rank == 0)THEN
PRINT*
PRINT*, "  INFO: Creating Fields"
PRINT*, "    Degrees of Freedom = ", F % DegreesOfFreedom
END IF

SELECT CASE (M % nDimensions)
CASE (1)
CALL DMDACreate1D( &
Comm,                 &
M % BoundaryType(1),  &
M % nZones(1),        &
F % DegreesOfFreedom, &
M % StencilWidth,     &
PETSC_NULL_INTEGER,   &
F % FieldsDA,         &
Error)

CASE (2)
CALL DMDACreate2D( &
Comm,                 &
M % BoundaryType(1),  &
M % BoundaryType(2),  &
M % StencilType,      &
M % nZones(1),        &
M % nZones(2),        &
PETSC_DECIDE,         &
PETSC_DECIDE,         &
F % DegreesOfFreedom, &
M % StencilWidth,     &
PETSC_NULL_INTEGER,   &
PETSC_NULL_INTEGER,   &
F % FieldsDA,         &
Error)

CASE (3)
CALL DMDACreate3D( &
Comm,                 &
M % BoundaryType(1),  &
M % BoundaryType(2),  &
M % BoundaryType(3),  &
M % StencilType,      &
M % nZones(1),        &
M % nZones(2),        &
M % nZones(3),        &
PETSC_DECIDE,         &
PETSC_DECIDE,         &
PETSC_DECIDE,         &
F % DegreesOfFreedom, &
M % StencilWidth,     &
PETSC_NULL_INTEGER,   &
PETSC_NULL_INTEGER,   &
PETSC_NULL_INTEGER,   &
F % FieldsDA,         &
Error)

END SELECT

CALL DMCreateGlobalVector( &
F % FieldsDA, F % FieldsGlobal, Error)

CALL DMCreateLocalVector( &
F % FieldsDA, F % FieldsLocal,  Error)

END SUBROUTINE CreateFields
```

### DestroyFields

DestroyFields frees memory allocated by CreateFields

```  SUBROUTINE DestroyFields(F)

TYPE(FieldsType), POINTER :: &
F

PetscInt :: &
Error

CALL VecDestroy(F % FieldsGlobal, Error)

CALL VecDestroy(F % FieldsLocal,  Error)

CALL DMDestroy(F % FieldsDA, Error)

END SUBROUTINE DestroyFields
```

The following PETSc routines are called

```!  Get Access to Mesh Coordinates:

CALL DMDAVecGetArrayF90( &
M % CoordinateDA, M % Positions % InnerEdgeGlobal, InnerEdge, Error)
CALL DMDAVecGetArrayF90( &
M % CoordinateDA, M % Positions % CenterGlobal,    Center,    Error)
CALL DMDAVecGetArrayF90( &
M % CoordinateDA, M % Positions % OuterEdgeGlobal, OuterEdge, Error)
```
```  !  Release Access to Coordinates:

CALL DMDAVecRestoreArrayF90( &
M % CoordinateDA, M % Positions % InnerEdgeGlobal, InnerEdge, Error)
CALL DMDAVecRestoreArrayF90( &
M % CoordinateDA, M % Positions % CenterGlobal,    Center,    Error)
CALL DMDAVecRestoreArrayF90( &
M % CoordinateDA, M % Positions % OuterEdgeGlobal, OuterEdge, Error)
```

### Set initial condition

To initialize FieldsGlobal, we first call DMDAVecGetArrayF90, which returns FieldsOld, a rank-four array (DegreesOfFreedom > 1) where the first dimension holds the individual fields and the next three dimensions hold the logical position in the mesh. After initialization of the data we call DMDAVecRestoreArrayF90.

```  !  Get Access to Non Ghosted Fields Data:

CALL DMDAVecGetArrayF90( &
F % FieldsDA, F % FieldsGlobal, FieldsOld, Error)

!  Set Initial Condition and Compute Time Step:

MyTimeStep = Huge

DO iX3 = iBX3, iEX3
DO iX2 = iBX2, iEX2
DO iX1 = iBX1, iEX1

FieldsOld(D,  iX1, iX2, iX3) &
= SIN( Two * Pi * Center(X1, iX1, iX2, iX3) )

FieldsOld(V1, iX1, iX2, iX3) &
= One
FieldsOld(V2, iX1, iX2, iX3) &
= Zero
FieldsOld(V3, iX1, iX2, iX3) &
= Zero

MyTimeStep &
= MIN(MyTimeStep, &
( OuterEdge(X1, iX1, iX2, iX3) &
- InnerEdge(X1, iX1, iX2, iX3) ) &
/ ( Tiny + FieldsOld(V1, iX1, iX2, iX3) ) )

MyTimeStep &
= MIN(MyTimeStep, &
( OuterEdge(X2, iX1, iX2, iX3) &
- InnerEdge(X2, iX1, iX2, iX3) ) &
/ ( Tiny + FieldsOld(V2, iX1, iX2, iX3) ) )

MyTimeStep &
= MIN(MyTimeStep, &
( OuterEdge(X3, iX1, iX2, iX3) &
- InnerEdge(X3, iX1, iX2, iX3) ) &
/ ( Tiny + FieldsOld(V3, iX1, iX2, iX3) ) )

END DO
END DO
END DO

CALL DMDAVecRestoreArrayF90( &
F % FieldsDA, F % FieldsGlobal, FieldsOld, Error)

!  Reduce Global Time Step:

CALL MPI_ALLREDUCE( &
MyTimeStep, TimeStep, 1, MPI_DOUBLE_PRECISION, MPI_MIN, &
U % Communicator, Error)

!  Populate Ghosted Fields with Initial Condition:

CALL DMGlobalToLocalBegin( &
F % FieldsDA, F % FieldsGlobal, INSERT_VALUES, F % FieldsLocal, Error)
CALL DMGlobalToLocalEnd( &
F % FieldsDA, F % FieldsGlobal, INSERT_VALUES, F % FieldsLocal, Error)

!  Keep Copy of Initial Condition:

CALL VecDuplicate(F % FieldsLocal, InitialCondition, Error)
CALL VecCopy(     F % FieldsLocal, InitialCondition, Error)
```

### Evolve initial condition

Mesh cells in a domain owned in particular process are divided into “exchange cells” (red) and “non-exchange cells” (green). A process’ exchange cells are ghost cells for neighboring processes.

The time step loop is organized in the following way to overlap communications and computations:

1. Update cells whose data is communicated with neighboring processes (AdvectExchangeCells3D; see AdvectionModule.F90); red cells in the above figure
2. Begin ghost exchange by calling DMGlobalToLocalBegin
3. Update cells whose data is not communicated with neighboring processes (AdvectNonExchangeCells3D; see AdvectionModule.F90); green cells in the above figure
4. Finish ghost exchange by calling DMGlobalToLocalEnd. Ghost cells are now populated with updated values, and ready for the next time step
```  !  Begin Time Step Loop:

Done      = .FALSE.
iCycle    = 0
MaxCycles = 100000
Time      = Zero
EndTime   = One

DO WHILE (.NOT.Done)

iCycle = iCycle + 1

IF(U % MPI_Rank == 0)THEN
PRINT*
PRINT*, "    Cycle     = ", iCycle
PRINT*, "    Time Step = ", TimeStep
PRINT*, "    Time      = ", Time
END IF

!  Copy Ghosted Data to Local Work Vector:

CALL VecCopy(F % FieldsLocal, FieldsScratch, Error)

CALL DMDAVecGetArrayF90( &
F % FieldsDA, FieldsScratch,    FieldsOld, Error)
CALL DMDAVecGetArrayF90( &
F % FieldsDA, F % FieldsGlobal, FieldsNew, Error)

!  Update Cells Whose Data is Communicated:

CALL AdvectExchangeCells3D(FieldsNew, FieldsOld, InnerEdge, OuterEdge, &
TimeStep, iBX1, iBX2, iBX3, iEX1, iEX2, iEX3, SW)

!  Begin Ghost Exchange:

CALL DMGlobalToLocalBegin( &
F % FieldsDA, F % FieldsGlobal, INSERT_VALUES, F % FieldsLocal, &
Error)

!  Update Cells Whose Data is Not Comminucated:

CALL AdvectNonExchangeCells3D(FieldsNew, FieldsOld, InnerEdge, OuterEdge, &
TimeStep, iBX1, iBX2, iBX3, iEX1, iEX2, iEX3, SW)

CALL DMDAVecRestoreArrayF90( &
F % FieldsDA, FieldsScratch,    FieldsOld, Error)
CALL DMDAVecRestoreArrayF90( &
F % FieldsDA, F % FieldsGlobal, FieldsNew, Error)

!  End Ghost Exchange:

CALL DMGlobalToLocalEnd( &
F % FieldsDA, F % FieldsGlobal, INSERT_VALUES, F % FieldsLocal, &
Error)

!  Copy New to Old:

CALL DMDAVecGetArrayF90( &
F % FieldsDA, F % FieldsLocal,  FieldsOld, Error)
CALL DMDAVecGetArrayF90( &
F % FieldsDA, F % FieldsGlobal, FieldsNew, Error)

DO iX3 = iBX3, iEX3
DO iX2 = iBX2, iEX2
DO iX1 = iBX1, iEX1
DO iField = 0, DOF - 1

FieldsOld(iField, iX1, iX2, iX3) &
= FieldsNew(iField, iX1, iX2, iX3)

END DO
END DO
END DO
END DO

CALL DMDAVecRestoreArrayF90( &
F % FieldsDA, F % FieldsLocal,  FieldsOld, Error)
CALL DMDAVecRestoreArrayF90( &
F % FieldsDA, F % FieldsGlobal, FieldsNew, Error)

Time = Time + TimeStep

IF(iCycle >= MaxCycles .OR. Time >= EndTime) Done = .TRUE.

END DO

!  Destroy Local Scratch Vector:

CALL VecDestroy(FieldsScratch, Error)
```

### Compare result with initial condition

```  !  Compare Evolved Solution with Initial Condition in the L1 Norm:

CALL DMDAVecGetArrayF90( &
F % FieldsDA, InitialCondition, FieldsOld, Error)
CALL DMDAVecGetArrayF90( &
F % FieldsDA, F % FieldsLocal,  FieldsNew, Error)

MyL1Norm = Zero
DO iX3 = iBX3, iEX3
DO iX2 = iBX2, iEX2
DO iX1 = iBX1, iEX1

MyL1Norm &
= MyL1Norm &
+ ABS(FieldsNew(D, iX1, iX2, iX3) - FieldsOld(D, iX1, iX2, iX3))

END DO
END DO
END DO

CALL MPI_REDUCE( &
MyL1Norm, L1Norm, 1, MPI_DOUBLE_PRECISION, MPI_SUM, 0, &
U % Communicator, Error)

IF(U % MPI_Rank == 0)THEN
PRINT*
PRINT*, "  INFO: L1 Norm = ", L1Norm
END IF

CALL DMDAVecRestoreArrayF90( &
F % FieldsDA, InitialCondition, FieldsOld, Error)
CALL DMDAVecRestoreArrayF90( &
F % FieldsDA, F % FieldsGlobal, FieldsNew, Error)

CALL VecDestroy(InitialCondition, Error)
```

AdvectionModule contains AdvectExchangeCells3D and AdvectNonExchangeCells3D which call the kernel Advect3D on exchange cells and non-exchange cells, respectively.

```  SUBROUTINE AdvectExchangeCells3D( &
FieldsNew, FieldsOld, InnerEdge, OuterEdge, TimeStep, &
iBX1, iBX2, iBX3, iEX1, iEX2, iEX3, SW)

PetscReal, DIMENSION(:,:,:,:), POINTER :: &
FieldsNew, &
FieldsOld, &
InnerEdge, &
OuterEdge
PetscReal, INTENT(in) :: &
TimeStep
PetscInt, INTENT(in) :: &
iBX1, iBX2, iBX3, &
iEX1, iEX2, iEX3, SW

PetscInt :: &
X1B, X1E, &
X2B, X2E, &
X3B, X3E

!  1:
X1B = iBX1
X1E = iBX1 + SW - 1
X2B = iBX2
X2E = iEX2
X3B = iBX3
X3E = iEX3
FieldsNew(D,                        &
X1B : X1E, X2B : X2E, X3B : X3E), &
FieldsOld(D,                        &
X1B - SW : X1E + SW,              &
X2B - SW : X2E + SW,              &
X3B - SW : X3E + SW),             &
FieldsNew(V1,                       &
X1B : X1E, X2B : X2E, X3B : X3E), &
OuterEdge(X1,                       &
X1B : X1E, X2B : X2E, X3B : X3E)  &
- InnerEdge(X1,                     &
X1B : X1E, X2B : X2E, X3B : X3E), &
TimeStep,                           &
X1B = X1B, X1E = X1E,               &
X2B = X2B, X2E = X2E,               &
X3B = X3B, X3E = X3E, SW = SW)

!  2:
X1B = iEX1 - SW + 1
X1E = iEX1
X2B = iBX2
X2E = iEX2
X3B = iBX3
X3E = iEX3
FieldsNew(D,                        &
X1B : X1E, X2B : X2E, X3B : X3E), &
FieldsOld(D,                        &
X1B - SW : X1E + SW,              &
X2B - SW : X2E + SW,              &
X3B - SW : X3E + SW),             &
FieldsNew(V1,                       &
X1B : X1E, X2B : X2E, X3B : X3E), &
OuterEdge(X1,                       &
X1B : X1E, X2B : X2E, X3B : X3E)  &
- InnerEdge(X1,                     &
X1B : X1E, X2B : X2E, X3B : X3E), &
TimeStep,                           &
X1B = X1B, X1E = X1E,               &
X2B = X2B, X2E = X2E,               &
X3B = X3B, X3E = X3E, SW = SW)

!  3:
X1B = iBX1 + SW
X1E = iEX1 - SW
X2B = iBX2
X2E = iBX2 + SW - 1
X3B = iBX3
X3E = iEX3
FieldsNew(D,                        &
X1B : X1E, X2B : X2E, X3B : X3E), &
FieldsOld(D,                        &
X1B - SW : X1E + SW,              &
X2B - SW : X2E + SW,              &
X3B - SW : X3E + SW),             &
FieldsNew(V1,                       &
X1B : X1E, X2B : X2E, X3B : X3E), &
OuterEdge(X1,                       &
X1B : X1E, X2B : X2E, X3B : X3E)  &
- InnerEdge(X1,                     &
X1B : X1E, X2B : X2E, X3B : X3E), &
TimeStep,                           &
X1B = X1B, X1E = X1E,               &
X2B = X2B, X2E = X2E,               &
X3B = X3B, X3E = X3E, SW = SW)

!  4:
IF(iEX2 > iBX2)THEN
X1B = iBX1 + SW
X1E = iEX1 - SW
X2B = iEX2 - SW + 1
X2E = iEX2
X3B = iBX3
X3E = iEX3
FieldsNew(D,                        &
X1B : X1E, X2B : X2E, X3B : X3E), &
FieldsOld(D,                        &
X1B - SW : X1E + SW,              &
X2B - SW : X2E + SW,              &
X3B - SW : X3E + SW),             &
FieldsNew(V1,                       &
X1B : X1E, X2B : X2E, X3B : X3E), &
OuterEdge(X1,                       &
X1B : X1E, X2B : X2E, X3B : X3E)  &
- InnerEdge(X1,                     &
X1B : X1E, X2B : X2E, X3B : X3E), &
TimeStep,                           &
X1B = X1B, X1E = X1E,               &
X2B = X2B, X2E = X2E,               &
X3B = X3B, X3E = X3E, SW = SW)
END IF

!  5:
X1B = iBX1 + SW
X1E = iEX1 - SW
X2B = iBX2 + SW
X2E = iEX2 - SW
X3B = iBX3
X3E = iBX3 + SW - 1
FieldsNew(D,                        &
X1B : X1E, X2B : X2E, X3B : X3E), &
FieldsOld(D,                        &
X1B - SW : X1E + SW,              &
X2B - SW : X2E + SW,              &
X3B - SW : X3E + SW),             &
FieldsNew(V1,                       &
X1B : X1E, X2B : X2E, X3B : X3E), &
OuterEdge(X1,                       &
X1B : X1E, X2B : X2E, X3B : X3E)  &
- InnerEdge(X1,                     &
X1B : X1E, X2B : X2E, X3B : X3E), &
TimeStep,                           &
X1B = X1B, X1E = X1E,               &
X2B = X2B, X2E = X2E,               &
X3B = X3B, X3E = X3E, SW = SW)

!  6:
IF(iEX3 > iBX3)THEN
X1B = iBX1 + SW
X1E = iEX1 - SW
X2B = iBX2 + SW
X2E = iEX2 - SW
X3B = iEX3 - SW + 1
X3E = iEX3
FieldsNew(D,                        &
X1B : X1E, X2B : X2E, X3B : X3E), &
FieldsOld(D,                        &
X1B - SW : X1E + SW,              &
X2B - SW : X2E + SW,              &
X3B - SW : X3E + SW),             &
FieldsNew(V1,                       &
X1B : X1E, X2B : X2E, X3B : X3E), &
OuterEdge(X1,                       &
X1B : X1E, X2B : X2E, X3B : X3E)  &
- InnerEdge(X1,                     &
X1B : X1E, X2B : X2E, X3B : X3E), &
TimeStep,                           &
X1B = X1B, X1E = X1E,               &
X2B = X2B, X2E = X2E,               &
X3B = X3B, X3E = X3E, SW = SW)
END IF

```

```SUBROUTINE AdvectNonExchangeCells3D( &
FieldsNew, FieldsOld, InnerEdge, OuterEdge, TimeStep, &
iBX1, iBX2, iBX3, iEX1, iEX2, iEX3, SW)

PetscReal, DIMENSION(:,:,:,:), POINTER :: &
FieldsNew, &
FieldsOld, &
InnerEdge, &
OuterEdge
PetscReal, INTENT(in) :: &
TimeStep
PetscInt, INTENT(in) :: &
iBX1, iBX2, iBX3, &
iEX1, iEX2, iEX3, SW

PetscInt :: &
X1B, X1E, &
X2B, X2E, &
X3B, X3E

X1B = iBX1 + SW
X1E = iEX1 - SW
X2B = iBX2 + SW
X2E = iEX2 - SW
X3B = iBX3 + SW
X3E = iEX3 - SW
FieldsNew(D,                        &
X1B : X1E, X2B : X2E, X3B : X3E), &
FieldsOld(D,                        &
X1B - SW : X1E + SW,              &
X2B - SW : X2E + SW,              &
X3B - SW : X3E + SW),             &
FieldsNew(V1,                       &
X1B : X1E, X2B : X2E, X3B : X3E), &
OuterEdge(X1,                       &
X1B : X1E, X2B : X2E, X3B : X3E)  &
- InnerEdge(X1,                     &
X1B : X1E, X2B : X2E, X3B : X3E), &
TimeStep,                           &
X1B = X1B, X1E = X1E,               &
X2B = X2B, X2E = X2E,               &
X3B = X3B, X3E = X3E, SW = SW)

```

```  SUBROUTINE Advect3D(D_NEW, D_OLD, V1, dX1, dt, &
X1B, X1E, X2B, X2E, X3B, X3E, SW)

PetscInt, INTENT(in) :: &
X1B, X1E, X2B, X2E, X3B, X3E, SW
PetscReal, DIMENSION(X1B : X1E, X2B : X2E, X3B : X3E), INTENT(inout) :: &
D_NEW
PetscReal, DIMENSION(X1B - SW : X1E + SW, &
X2B - SW : X2E + SW, &
X3B - SW : X3E + SW), INTENT(in) :: &
D_OLD
PetscReal, DIMENSION(X1B : X1E, X2B : X2E, X3B : X3E), INTENT(in) :: &
V1, &
dX1
PetscReal, INTENT(in) :: &
dt

PetscReal :: &
FluxOuterX1, &
FluxInnerX1
PetscInt :: &
iX1, iX2, iX3

DO iX3 = X3B, X3E
DO iX2 = X2B, X2E
DO iX1 = X1B, X1E

FluxInnerX1 &
= Half * ( ( V1(iX1, iX2, iX3) &
+ ABS(V1(iX1, iX2, iX3)) ) * D_OLD(iX1-1, iX2, iX3) &
+ ( V1(iX1, iX2, iX3) &
- ABS(V1(iX1, iX2, iX3)) ) * D_OLD(iX1, iX2, iX3) )

FluxOuterX1 &
= Half * ( ( V1(iX1, iX2, iX3) &
+ ABS(V1(iX1, iX2, iX3)) ) * D_OLD(iX1, iX2, iX3) &
+ ( V1(iX1, iX2, iX3) &
- ABS(V1(iX1, iX2, iX3)) ) * D_OLD(iX1+1, iX2, iX3) )

D_NEW(iX1, iX2, iX3) &
= D_OLD(iX1, iX2, iX3) &
- dt * ( FluxOuterX1 - FluxInnerX1 ) / dX1(iX1, iX2, iX3)

END DO
END DO
END DO

```

## Compiling

```\$ module load cray-petsc
\$ module load git
\$ git clone https://github.com/olcf/PETSc-tutorial.git
\$ cd PETSc-tutorial
\$ mkdir Executables
\$ cd Executables
\$ ftn -c ../RealKindModule.F90 ../MeshModule.F90 ../FieldsModule.F90 ../UniverseModule.F90 ../AdvectionModule.F90 ../Advection3D.F90
\$ ftn -o Advection3D RealKindModule.o MeshModule.o FieldsModule.o UniverseModule.o AdvectionModule.o Advection3D.o
```

## Running

```\$ aprun -n16 ./Advection3D

INFO: Creating Universe Advection3D

INFO: Number of MPI Ranks =            16

INFO: Creating Mesh
Dimensionality =             3
Stencil Width  =             1
iDim, nZones =             1          256
Inner Boundary   =     0.000000000000000
Outer Boundary   =     1.000000000000000
iDim, nZones =             2          256
Inner Boundary   =     0.000000000000000
Outer Boundary   =     1.000000000000000
iDim, nZones =             3          256
Inner Boundary   =     0.000000000000000
Outer Boundary   =     1.000000000000000

INFO: Creating Fields
Degrees of Freedom =             4

Cycle     =             1
Time Step =    3.9062500000000000E-003
Time      =     0.000000000000000

Cycle     =             2
Time Step =    3.9062500000000000E-003
Time      =    3.9062500000000000E-003

Cycle     =             3
Time Step =    3.9062500000000000E-003
Time      =    7.8125000000000000E-003

...

Cycle     =           256
Time Step =    3.9062500000000000E-003
Time      =    0.9960937500000000

INFO: Finished Evolution
Cycles Evolved =           256
Time           =     1.000000000000000

INFO: L1 Norm =    2.2737367544323206E-013

INFO: Destroying Universe Advection3D

INFO:            16  MPI Ranks Finished in     109.6593058109283
Seconds
```