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

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

Equation 1.2

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

Equation 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.

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

stencil

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

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

overlap

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)

    !  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