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