Using petsc (parallel)

And here is a quick parallel code - compile as in using petsc but run like

mpirun -n <N procS> ./a.out <other options>

On my workstation, n = 2 actually results in better performance than n = 1(it has four cores but usually multicore cpus do nothing much for petsc), though not if I bump the problem up to 1024x1024

On bluecrystalp1, running a 1024x1024 problem with -pc_type sor -ksp_type bcgs (with parameters picked to make iterative solution hard), I get

1 node, 1 core : 70s 
1 node, 2 cores  : 45s
1 node, 4 cores : 32s
2 nodes 4 cores each : 14s
4 nodes, 1 core each : 23s

compared to the 49s required on my workstation with one core. Asking for 4 nodes seems to result in lengthy queuing.

using mumps instead (-pc_type lu -pc_factor mat_solver_package mumps) I get

1 node, 1 core : 77s
1 node, 2 cores: 48s
1 node, 4 cores: 31s
2 node, 1 core each: 47s
2 node, 2 cores each: 29s
2 node, 4 cores each: 26s
4 node, 2 cores each: 26s

Does this suggest that the iterative method scales better across nodes? perhaps not suprising as all it needs to do is matrix vector multiplies, which would mean the number of inteprocess comms growing like O(nCore) only.

results for BC2 & bcgs are encouraging

nodes	cores	time (s)
1	1	67

 
1 node, 4 cores  : 30.15s
2 node,  4 core each : 
4 nodes, 4 core each : 4.8s
8 nodes, 4 core each : 1.7s

!! solve a 2D n x n (a(x,y)u')' + b(x,y) u = c(x,y) 
!! with periodic boundary conditions in y, u(x=0,y) = 0 and u'(x=1,y)= d
!! a(x,y) = 1 for x < 1/2, h for x > 1/2
!! b(x,y) = 1 for x < 1/2, r for x > 1/2 
!! c(x,y) = 1 for x < 1/2, p for x > 1/2 

program main

  implicit none

#include "finclude/petsc.h"
#include "finclude/petscvec.h"
#include "finclude/petscmat.h"
#include "finclude/petscksp.h"
#include "finclude/petscpc.h"
 
  
  Mat A
  Vec x,b, res
  KSP ksp
  PetscInt ierr, ndof, n, m, i , j, lbegin, lend
  PetscScalar dfluxe(0:1), dfluxw(0:1), dfluxn(0:1), dfluxs(0:1), su(0:0), sc, &
       h, r, p , d, rnorm, bnorm, one
  PetscInt row(0:0), fcol(0:1)
  
  n = 256
  m = n/2
  ndof = n * n

  !x > 1/2 diffusion = h u'
  h = 1e-2
  !x > 1/2 dissipation = r u , set small r to have problems inverting the matrix
  !so that direct solvers win
  r = 1e-6
  !x > 1/2 source = p
  p = 0.5
  ! u' at x = 1
  d = 100
  

  ! this has to be the first petsc call 
  call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
  
  !allocate space for matrix elements
  call MatCreateMpiAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, &
       ndof, ndof, 5, PETSC_NULL_INTEGER, 1,PETSC_NULL_INTEGER, A,ierr)
  ! alocate space for rhs, residual and solution
  call MatGetVecs(A, x, b, ierr);
  call MatGetVecs(A, x, res, ierr);
  

  !establish which dofs are on this processor
  call MatGetOwnerShipRange(A, lbegin, lend, ierr)

  ! bulk
  dfluxe(0) = 1.0
  dfluxe(1) = -1.0
 
  su(0) = 1.0
  sc = 1.0

  do i = 1, n-2

     dfluxw(0) = dfluxe(0)
     dfluxw(1) = dfluxe(1)

     
     dfluxn(0) = 1.0
     dfluxn(1) = -1.0

  

     if (i > m) then
        su(0) = r
        sc = p
        dfluxe(0) = h
        dfluxe(1) = -h
        dfluxn(0) = h
        dfluxn(1) = -h
     endif
     
     dfluxs(0) = dfluxn(0)
     dfluxs(1) = dfluxn(1)
    
     do j = 1, n-2
          
        row(0) = j*n + i
        if ((lbegin .le. row(0)) .and. (lend .gt. row(0))) then
        
           ! source term
           call MatSetValues(A, 1 ,row, 1, row, su, ADD_VALUES,ierr)
           call VecSetValue( b, row(0) , sc , INSERT_VALUES, ierr)
           
           ! fluxes
           fcol(0) = row(0)
           
           fcol(1) = row(0) + 1
           call MatSetValues(A, 1 ,row, 2, fcol, dfluxe, ADD_VALUES,ierr)
           
           fcol(1) = row(0) - 1
           call MatSetValues(A, 1 ,row, 2, fcol, dfluxw, ADD_VALUES,ierr)
           
           fcol(1) = row(0) + n
           call MatSetValues(A, 1 ,row, 2, fcol, dfluxn, ADD_VALUES,ierr)
           
           fcol(1) = row(0) - n
           call MatSetValues(A, 1 ,row, 2, fcol, dfluxs, ADD_VALUES,ierr)
        end if
        
     end do

     ! periodic boundaries
     j = 0
     row(0) = j*n + i
     if ((lbegin .le. row(0)) .and. (lend .gt. row(0))) then
        ! source term
        call MatSetValues(A, 1 ,row, 1, row, su, ADD_VALUES,ierr)
        call VecSetValue( b, row(0) , sc , INSERT_VALUES, ierr)
        
        ! fluxes
        fcol(0) = row(0)
        
        fcol(1) = row(0) + 1
        call MatSetValues(A, 1 ,row, 2, fcol, dfluxe, ADD_VALUES,ierr)
        
        fcol(1) = row(0) - 1
        call MatSetValues(A, 1 ,row, 2, fcol, dfluxw, ADD_VALUES,ierr)
        
        fcol(1) = row(0) + n
        call MatSetValues(A, 1 ,row, 2, fcol, dfluxn, ADD_VALUES,ierr)
        
        fcol(1) = (n-1)*n + i
        call MatSetValues(A, 1 ,row, 2, fcol, dfluxs, ADD_VALUES,ierr)
     end if

     j = n-1
     row(0) = j*n + i
     if ((lbegin .le. row(0)) .and. (lend .gt. row(0))) then
        ! source term
        call MatSetValues(A, 1 ,row, 1, row, su, ADD_VALUES,ierr)
        call VecSetValue( b, row(0) , sc , INSERT_VALUES, ierr)
        
        ! fluxes
        fcol(0) = row(0)
        
        fcol(1) = row(0) + 1
        call MatSetValues(A, 1 ,row, 2, fcol, dfluxe, ADD_VALUES,ierr)
        
        fcol(1) = row(0) - 1
        call MatSetValues(A, 1 ,row, 2, fcol, dfluxw, ADD_VALUES,ierr)
        
        fcol(1) = row(0) - n
        call MatSetValues(A, 1 ,row, 2, fcol, dfluxn, ADD_VALUES,ierr)
        
        fcol(1) = i
        call MatSetValues(A, 1 ,row, 2, fcol, dfluxs, ADD_VALUES,ierr)
     end if
     
  end do
   
  !x-boundaries

  dfluxw(0) = h
  dfluxw(1) = -h

  do j = 0, n - 1
  
     i = 0
     row(0) = j*n + i
     if ((lbegin .le. row(0)) .and. (lend .gt. row(0))) then
        call MatSetValues(A, 1 ,row, 1, row, su, ADD_VALUES,ierr)
     end if
     
     
     i = n - 1
     row(0) = j*n + i
     if ((lbegin .le. row(0)) .and. (lend .gt. row(0))) then
        su = r
        sc = h*d
        call MatSetValues(A, 1 ,row, 1, row, su, ADD_VALUES,ierr)
        call VecSetValue( b, row(0) , sc , INSERT_VALUES, ierr)
        fcol(0) = row(0)
        fcol(1) = row(0) - 1
        call MatSetValues(A, 1 ,row, 2, fcol, dfluxw, ADD_VALUES,ierr) 
     end if
     


  end do


  !needs to be done before A can be used
  call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
  call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)

 
  !KSP context: KSP is petsc's interface to linear solvers
  call KSPCreate(PETSC_COMM_WORLD,ksp,ierr)

  !A is defines both the operator and preconditioner (but might not)
  call KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN,ierr)

  !grab all the options from the command line
  call KSPSetFromOptions(ksp,ierr)

  ! solve
  call KSPSolve(ksp,b,x,ierr)

  !check results
  call MatMult(A, x, res, ierr)
  one = -1.0
  call VecAXPY(res , one , b, ierr)
  call VecNorm(res , NORM_2, rnorm, ierr) 
  call VecNorm(b , NORM_2, bnorm, ierr)
  write(*,*) "||Ax - b||/||b|| = ", rnorm/bnorm


  !clean up
   call VecDestroy(x,ierr)
   call VecDestroy(b,ierr)
   call MatDestroy(A,ierr)

  !last petsc call
  call PetscFinalize(ierr)

end program main