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:

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 Laplace

2D Laplace

1D PseudoLaplace

2D PseudoLaplace

1D Helmholtz

2D Helmholtz

2D Navier-Stokes

© CNRS - IDRIS, 13/01/2012