PROGRAM poisson
USE MPI
IMPLICIT NONE
!******************************************************************************
!*  Solution of Poisson Equation on a rectangular domain
!*      x=(0,...,XMAX) y=(0,...,YMAX)
!* by means of Jacobi iteretion method (to be parallelized with MPI)
!* the source term is rho(x,y) and the psi is zero on the boundary
!*  10 nov 2005             by Andrea Bertoni
!*
!* Dividiamo il dominio in numprocs sottodomini. Ciascun processore
!* calcolerà il psi relativo al proprio dominio. I vari psi verranno
!* poi raccolti dal master(processore 0). I nodi comunicheranno
!* tra di loro per implementare l'algoritmo di Jacobi.
!******************************************************************************

!Dimensioni dominio
REAL, PARAMETER    :: XMAX= 1., YMAX= 1.
!Punti sulla griglia
INTEGER, PARAMETER :: NUMX= 100, NUMY= 100
!Max numero di iterazioni
INTEGER, PARAMETER :: NITERMAX=50000
!Errore massimo permesso(serve per terminare il ciclo)
REAL, PARAMETER    :: MAXERROR= 1e-6

REAL*8, ALLOCATABLE :: psi (:,:)    ! potenziale
REAL*8, ALLOCATABLE :: psinew (:,:)
!rho: rappresenta la densità di carica
REAL*8 :: rho(NUMX, NUMY)
!Locazione sulla griglia
REAL*8 :: x (NUMX), y(NUMY)
!Rappresentano h lungo x e lungo y
REAL*8 :: deltax, deltay
!Sommatori che tengono conto della differenza della psinew
!tra un'iterazione e l'altra
REAL*8 :: sumdiff, sumdiff2

INTEGER :: nx, ny, niter, i, j
INTEGER :: myid,numprocs,ierr
INTEGER :: STATUS(MPI_STATUS_SIZE)
!Dimensione dei sottodomini affidati ai processi
INTEGER :: fetta_dominio
INTEGER :: resto
!Per ogni processo corrispondono alla prima e ultima riga
INTEGER, ALLOCATABLE :: riga_inizio(:), riga_fine (:)



 CALL MPI_INIT( ierr )
 CALL MPI_COMM_RANK( MPI_COMM_WORLD, myid, ierr )
 CALL MPI_COMM_SIZE( MPI_COMM_WORLD, numprocs, ierr ) ! il numero di processori

 ALLOCATE(riga_inizio(0:numprocs-1),riga_fine(0:numprocs-1),stat=ierr)

 !Come viene ripartito il dominio
 fetta_dominio= NUMY / numprocs
 resto= mod(NUMY,numprocs)

 DO i=0, numprocs-1
     IF (i==0) THEN
       riga_inizio(i)=1
     ELSE
       riga_inizio(i)=riga_fine(i-1)+1
     END IF

     IF (i<resto) THEN
       riga_fine(i)=riga_inizio(i)+fetta_dominio
     ELSE
       riga_fine(i)=riga_inizio(i)+fetta_dominio-1
     END IF

 END DO

 DO i=0, resto-1
    IF (myid==i) fetta_dominio=fetta_dominio+1
 END DO



 !alloco psinew
 ALLOCATE(psinew(0:NUMX+1,0:fetta_dominio+1), stat=ierr)
 !alloco psi con dimensioni diverse in base al processore
 IF (myid==0) THEN
    ALLOCATE(psi(0:NUMX+1,0:NUMY+1), stat=ierr)
 ELSE
    ALLOCATE(psi(0:NUMX+1,0:fetta_dominio+1), stat=ierr)
 END IF

 !inizializzo psinew e psi
 psinew(:,:)=0.0
 psi(:,:)=0.0

 !calcolo h lungo x e y
 deltax = XMAX / NUMX
 deltay= YMAX / NUMY

 !definisco le posizioni dei punti
 DO nx = 1, NUMX
    x (nx)= nx * deltax
 END DO

 DO ny= 1, NUMY
    y(ny)= ny * deltay
 END DO

 !la densità viene messa a zero
 rho(:,:)=0.0

 !se mi trovo nell'intervallo [0.4,0.2] e [0.6,0.5]
 !imposto la rho


IF (myid==0) THEN

  DO nx= 1, NUMX    !inizializza la matrice di densità di carica
    DO ny= 1, NUMY  !per il processo zero
      IF (x(nx) >= 0.4 .AND. x(nx) <= 0.6) THEN
    IF (y(ny) >= 0.2 .AND. y(ny) <= 0.5) THEN
      rho(nx,ny)= -(4.*3.1415926) * (- 2.)
    END IF
      END IF
    END DO
  END DO

 ELSE

  DO nx= 1, NUMX    !inizializza la matrice di densità di carica
    DO ny= riga_inizio(myid), riga_fine(myid)
      IF (x(nx) >= 0.4 .AND. x(nx) <= 0.6) THEN
    IF (y(ny) >= 0.2 .AND. y(ny) <= 0.5) THEN
       rho(nx,(ny-riga_inizio(myid)+1))=-(4.*3.1415926)*(-2.)
    END IF
      END IF
    END DO
  END DO

 END IF

!Setto la psi uguale alla psinew del ciclo precedente e la
!aggiorno.

 DO niter = 1, NITERMAX

    DO i=1,NUMX
    DO j=1,fetta_dominio
           psi(i,j)=psinew(i,j)
    END DO
    END DO

 sumdiff = 0.

 !il processo zero spedisce la sua riga più
 !alta al processo successivo(1) e riceve
 !la più bassa del processo successivo
 IF (myid==0) THEN

    CALL MPI_RECV(psi(:,fetta_dominio+1),NUMX,MPI_DOUBLE_PRECISION,1,0, MPI_COMM_WORLD,STATUS,ierr )
    CALL MPI_SEND(psi(:,fetta_dominio),NUMX,MPI_DOUBLE_PRECISION,1,1,MPI_COMM_WORLD,ierr)

 !l'ultimo processo spedisce la sua riga
 !più bassa al processo precedente e
 !riceve la più alta del processo precedente
 ELSE IF (myid==numprocs-1) THEN

    CALL MPI_SEND(psi(:,1),NUMX,MPI_DOUBLE_PRECISION,myid-1,0,MPI_COMM_WORLD,ierr)
    CALL MPI_RECV(psi(:,0),NUMX,MPI_DOUBLE_PRECISION,myid-1,1, MPI_COMM_WORLD,STATUS,ierr )

 ELSE

 !tutti i processi devono inviare la loro riga
 !più bassa al processo precedente e ricevere
 !la riga più bassa del processo successivo
 CALL MPI_SENDRECV(psi(:,1),NUMX, MPI_DOUBLE_PRECISION, myid-1,0,  &
                &     psi(:,fetta_dominio+1),NUMX, MPI_DOUBLE_PRECISION, myid+1,0,  &
                &     MPI_COMM_WORLD, STATUS, ierr)
 !tutti i processi devono inviare la loro riga più
 !alta al processo successivo e ricevere la riga
 !più alta del processo precedente
 CALL MPI_SENDRECV(psi(:,fetta_dominio),NUMX, MPI_DOUBLE_PRECISION, myid+1,1,  &
                &     psi(:,0),NUMX, MPI_DOUBLE_PRECISION, myid-1,1,  &
                &     MPI_COMM_WORLD, STATUS, ierr)

END IF

  DO nx = 1,NUMX
    DO ny = 1,fetta_dominio
    psinew(nx,ny)= 0.25 * (-deltax*deltay*rho(nx,ny) + psi(nx-1,ny) + psi(nx,ny-1) + psi(nx,ny+1) + psi(nx+1,ny))
    sumdiff= sumdiff + ABS(psi(nx,ny) - psinew(nx,ny))
    END DO
  END DO

  !raccolgo in sumdiff2 i vari sumdiff
  CALL MPI_REDUCE(sumdiff,sumdiff2,1,MPI_DOUBLE_PRECISION,MPI_SUM,0,MPI_COMM_WORLD,ierr)

  !se sono il processo zero calcolo sumdiff2
  IF (myid==0) THEN
    sumdiff2=sumdiff2/(NUMX*NUMY)   ! per trovare la media delle differenze
    write (*,*) myid, "sumdiff2= ",sumdiff2
  END IF

  CALL MPI_BCAST(sumdiff2,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)

  IF (sumdiff2<MAXERROR) then
    write(*,*) myid,": sumddiff= ",sumdiff2
    write(*,*) "num: 8155 "
    EXIT
  END IF

END DO


  IF (myid==0) THEN
    !stampo il numero di iterazioni
    WRITE (*,*) "niter = ", NITER
    !
    DO j=fetta_dominio+1,NUMY
    CALL MPI_RECV(psi(:,j),NUMX,MPI_DOUBLE_PRECISION,MPI_ANY_SOURCE,j, MPI_COMM_WORLD,STATUS,ierr )
    END DO

    OPEN(23, FILE="psi.dat", STATUS="UNKNOWN", FORM= "FORMATTED")
    DO nx= 1, NUMX  !salva la matrice dei risultati
      DO ny= 1, NUMY
    WRITE(23,*) x(nx), y(ny), psi(nx,ny)
      END DO
      WRITE(23,*)
    END DO
    CLOSE(23)

    OPEN(24, FILE="rho.dat", STATUS="UNKNOWN", FORM= "FORMATTED")
    DO nx= 1, NUMX
      DO ny= 1, NUMY
    WRITE(24,*) x(nx), y(ny), rho(nx,ny)/(-4.*3.1415)
      END DO
      WRITE(24,*)
    END DO
    CLOSE(24)

 ELSE

 !spedisco il risultato al processo zero
    DO i=1,fetta_dominio
    CALL MPI_SEND(psi(:,i),NUMX,MPI_DOUBLE_PRECISION,0,(riga_inizio(myid)+i-1),MPI_COMM_WORLD,ierr)
    END DO
 END IF

 CALL MPI_FINALIZE(ierr)

  DEALLOCATE(psi)
  DEALLOCATE(psinew)
  DEALLOCATE(riga_inizio,riga_fine)

END PROGRAM