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).
• $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.
Advection3D.F90
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
The following PETSc routines are called
• DMDACreate1D
• DMDACreate2D
• DMDACreate3D
• DMDAGetCoordinateDA
• DMCreateGlobalVector
• VecDuplicate
• DMDAGetCorners
• DMCreateLocalVector
• DMDAGetGhostCorners
• DMDAVecGetArrayF90
• DMDAVecRestoreArrayF90
• DMGlobalToLocalBegin
• DMGlobalToLocalEnd
• VecDestroy
• DMDestroy
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
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
Advection3D.F90 (continued)
The following PETSc routines are called
• VecDuplicate
• DMDAVecGetArrayF90
• DMDAVecRestoreArrayF90
• VecCopy
• DMGlobalToLocalBegin
• DMGlobalToLocalEnd
• VecDestroy
! 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
! Release Access to Non-Ghosted Fields:
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
The time step loop is organized in the following way to overlap communications and computations:
- Update cells whose data is communicated with neighboring processes (AdvectExchangeCells3D; see AdvectionModule.F90); red cells in the above figure
- Begin ghost exchange by calling DMGlobalToLocalBegin
- Update cells whose data is not communicated with neighboring processes (AdvectNonExchangeCells3D; see AdvectionModule.F90); green cells in the above figure
- 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)
! Get Access to Fields Data:
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)
! Release Access to Fields Data:
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.F90
AdvectionModule contains AdvectExchangeCells3D and AdvectNonExchangeCells3D which call the kernel Advect3D on exchange cells and non-exchange cells, respectively.
AdvectExchangeCells3D
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
CALL Advect3D( &
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
CALL Advect3D( &
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
CALL Advect3D( &
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
CALL Advect3D( &
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
CALL Advect3D( &
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
CALL Advect3D( &
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
END SUBROUTINE AdvectExchangeCells3D
AdvectNonExchangeCells3D
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
CALL Advect3D( &
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 SUBROUTINE AdvectNonExchangeCells3D
Advect3D
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
END SUBROUTINE Advect3D
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
