PMD Example: The 1D Poisson Problem


Let us consider the following classical 1D Poisson problem:

[Poisson Equation]

PMD splits the global domain into as many subdomains as the user required processes. At this time, PMD only implements a 1D Domain decomposition ("1DD"). Each subdomain is divided into regular meshes which step size and local node number are denoted here by h and Nx respectively.
Standard second order centered differences is used to discretize the second order operator so that, in each subdomain, the previous global system now reads:

[Discretized Poisson Equation]

Let us denote by L, D and Up the lower, the main and the upper diagonals of the local operator matrix. All the user has to provide to PMD is two subroutines: MyLocalSolver() and MyFirstDerivative() (see source code examples included in the PMD package). The former, to solve the previous local linear system and the latter, to compute the normal first derivative at the subdomain interfaces. The source code, hereafter, parallel solve such a problem using few PMD routine calls:

PROGRAM Poisson
 USE PMD                            
 IMPLICIT NONE

 !... Local constants and variables
 INCLUDE 'mpif.h'                            
 TYPE(PMD_Comm_1D)               :: Comm
 TYPE(PMD_R8_SYCH)               :: Factor
 TYPE(PMD_R8_Schur)              :: Schur
 TYPE(PMD_R8_Tr_Operator)        :: Operator
 INTEGER, PARAMETER              :: Nx=1001
 REAL(KIND=R8),DIMENSION(Nx)     :: L=0, D=0, Up=0, s, x
 REAL(KIND=R8),DIMENSION(0:Nx+1) :: u, ua
 REAL(KIND=R8)                   :: h, Local_norm, Global_norm
 INTEGER                         :: rank, Nsd, Info, i, k

 !... External routines
 EXTERNAL MPI_INIT,MPI_FINALIZE,MyLocalSolver,MyFirstDerivative

 !... Initiate MPI
 CALL MPI_INIT(Info)

 !... Initiate PMD
 CALL PMD_Init_1DD(MPI_COMM_WORLD, Comm)
 CALL PMD_Comm_Info_1DD(Comm, Rank=rank, Size=Nsd)

 !... Set local operator matrix (here Tridiagonal)
 h       = 1./(Nsd*(Nx - 1))
 L(2:Nx) = -1./h**2 ; D(1:Nx) = 2./h**2 ; Up(1:Nx-1) = -1./h**2

 !... Link to PMD Operator object
 CALL PMD_Operator_1DD(Comm, L, D, Up, Operator)

 !... Build Dual-Schur complement
 CALL PMD_Schur_1DD(Comm, MyLocalSolver, MyFirstDerivative, Operator, Schur)

 !... Parallel LU-Factor the Schur matrix
 CALL PMD_Schur_Factor_1DD(Comm, Schur, Factor)

 !... Free the Schur Matrix memory allocation
 CALL PMD_Schur_Free_1DD(Comm, Schur)

 !... Set source term and analytic solution
 x(:)  = (/ ((i+rank*(Nx-1)-1)*h, i=1,Nx) /)          !... X coordinate
 ua(1:Nx) = sin(acos(-1.)*x(:)) * cos(acos(-1.)*x(:)) !... analytic solution
 !... To reach values one step outside my subdomain (e.g. ua(0) and ua(Nx+1))
 call PMD_Exchange_1DD(Comm, ua)
 !... Compute the source term
 DO i = 1, Nx
     s(i) = - ( ua(i-1) - 2. * ua(i) + ua(i+1) ) / h**2
 ENDDO

 !... Physical boundary conditions
 u(1)  = 0.  ! WEST side  (x=0)
 u(Nx) = 1.  ! EAST side  (x=1)

 !... Parallel solve using a General Direct Solver
 CALL PMD_Solve_1DD(Comm, MyLocalSolver, MyFirstDerivative, Operator, &
                    Factor, s, u)

 !... Compute the norm of the residu
 Local_norm = MAXVAL(ABS(u(1:Nx) - ua(1:Nx))) / REAL(Nx * Nsd)
 CALL MPI_ALLREDUCE(Local_norm, Global_norm, 1, MPI_DOUBLE_PRECISION, MPI_MAX, &
                    MPI_COMM_WORLD, info)

 !... Check results
 IF( rank == 0 ) PRINT *,"Gap between analytic & numerical solutions: ",Global_norm

 !... Leave PMD
 CALL PMD_End_1DD(Comm)

 !... End MPI
 CALL MPI_FINALIZE(Info)
END PROGRAM Poisson

As we notice from the previous program, PMD hides to the user the subdomain interface problem. Also, from the user point of view, the distributed problem reduces to a local one since data can be viewed as local to each process memory. Especially noteworthy is that PMD doesn't depend on the user discretization technique. It also let the user free to define his local solver and his own method to compute the local normal first derivative at the subdomain interfaces.
The PMD package includes some examples involving others linear second order operators:

1D Laplace2D Laplace1D PseudoLaplace2D PseudoLaplace1D Helmholtz2D Helmholtz2D Navier-Stokes