PartitionCeff.f90

Go to the documentation of this file.

!!****
!!
!!  NAME
!!     Partition_Ceff
!!
!!  FUNCTION
!!     This subroutine partitions the full effective capacitance matrix
!!     into the blocks corresponding to the prescribed-prescribed, 
!!     prescribed-free and free-free DOFs.  This is accomplished by converting
!!     Ceff into its explicit form, partitioning, and converting the free-free
!!     part back to CRS format.  The prescribed-prescribed and free-prescribed
!!     blocks can stay in explicit form because they are not going to be input
!!     to BlockSolve95.
!!
!!  INPUTS
!!     ndim -- The size of one dimension of the global Ceff matrix
!!     nrows -- The number of rows assigned to this processor
!!     nnz -- The number of nonzeros in section of the Ceff matrix on this processor
!!     nstart -- The global index of the first row assigned to this processor
!!     rp1 -- The row mapping vector
!!     cval -- The collumn mapping vector
!!     aval -- The nonzero value vector
!!     ProcTemp -- The part of the global temperature vector that's assigned to this processor
!!
!!  OUTPUTS
!!     newnrows -- The number of rows assigned to this proc after BCs have been removed
!!     newnstart -- The global index of the first row assigned to this proc after BCs have been removed
!!     newndim -- The size of the global CeffFF matrix after BCs have been removed
!!
!!  USES
!!     none
!!
!!****

SUBROUTINE partition_ceff(ndim,nrows,nnz,nstart,rp1,cval,aval,newnrows,newnstart,newndim,ProcTemp,global)

  USE removebcs_global
  USE comp_row_global
  USE implicit_global
  USE precision
  USE rocstar_rocfraccomm

  IMPLICIT NONE

  include 'mpif.h'

  TYPE(rocfrac_global) :: global

  ! ... Arguments
  INTEGER :: ndim, nnz, nstart, nrows
  REAL(kind=wp), DIMENSION(nnz) :: aval
  INTEGER, DIMENSION(nnz) :: cval
  INTEGER, DIMENSION(nrows+1) :: rp1
  REAL(kind=wp), DIMENSION(LNumNp) :: proctemp

  INTEGER :: newnrows
  INTEGER :: newnstart
  INTEGER :: newndim

  ! ... local variables
  INTEGER :: numdisp  ! Local number of displacement BCs
  INTEGER, ALLOCATABLE, DIMENSION(:) :: dispbc  ! Local displacement BCs
  INTEGER, ALLOCATABLE, DIMENSION(:) :: tempintv
  
  INTEGER :: i, j, n, m, counter1, counter2, counter3, counter4, counter5

 
  INTEGER :: ncols, idof, jdof

  REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:) :: ceffexplicit, ceffff, cefffp, ceffpf, ceffpp


  ! ... Deallocate removeBCs_global variables
  IF(ALLOCATED(gtempbc)) DEALLOCATE(gtempbc)
  IF(ALLOCATED(numtempproc)) DEALLOCATE(numtempproc)

  ! ... Count the number of temperature boundary conditions in this 
  ! ... processes partition.
  numdisp = 0
  DO i = 1, global%NumNp
     IF (nodeproc(i) == myid) THEN
        ! ... If the value of the flag at node i on in this processes
        ! ... is 8, add 1 to the count of temperature BCs.
        IF(node_flag(i,1) == 8) THEN
           numdisp = numdisp + 1
        ENDIF

     ENDIF
  ENDDO
  print*,myid,' number of local disp BCs = ',numdisp

  ! Sum the local disp BCs, then broadcast to all procs
  CALL mpi_barrier(rocstar_communicator,ierr)

  ! ... MPI_REDUCE sums values on all processes to one value. (SUM(numdisp)=GNumTemp)
  CALL mpi_reduce(numdisp, gnumtemp, 1, mpi_integer, mpi_sum, 0 , rocstar_communicator, ierr)
  ! ... The number GNumTemp is then communicated to all processes, all know the global
  ! ... number of temperature boundary condtions.
  CALL mpi_bcast(gnumtemp, 1, mpi_integer, 0, rocstar_communicator, ierr)

  ! ... newndim is the size of the Global capacitance matrix after the temperature BC's are removed.
  newndim = gnumnp - gnumtemp
  print*,myid,' number of global disp BCs = ',gnumtemp

  ! ... Figure out which nodes have perscribed temperatures
  ! ... If this processes partition contains prescribed temperatures, then
  ! ... allocate an array to store the global index of that temperature BC.
  IF(numdisp > 0) THEN
     ALLOCATE(dispbc(1:numdisp))
     dispbc(:) = 0
     counter2 = 0

     DO m = 1, gnumnp
        DO i = 1, global%NumNp
           ! ... If global node m is owned by this process and there is a temperature
           ! ... BC at that node, record the node number as the next entry in array dispbc.
           IF (local2global(i) == m) THEN
              IF (nodeproc(i) == myid) THEN
                 IF(node_flag(i,1) == 8) THEN  ! Imposed constant nodal displacement
                    counter2 = counter2 + 1
                    dispbc(counter2) = m
                  ENDIF
              ENDIF
           ENDIF
        ENDDO
     ENDDO
  ENDIF
  If(numdisp >0) print*,myid,' local disp bcs = ',dispbc(:)

  ! ... Have each proc in turn broadcast its number of temperature BCs to the other procs
  ALLOCATE(numtempproc(1:nprocs))
  numtempproc(:) = 0
  numtempproc(myid+1) = numdisp
  DO i = 1, nprocs
     CALL mpi_bcast(numtempproc(i), 1, mpi_integer, i-1, rocstar_communicator, ierr)
  ENDDO

  ! ... The global row index of the first row on this process is the old index
  ! ... minus the number of temperature BCs on global nodes in the range from 
  ! ... node 1 to node nstart(the old starting row index).
  newnstart = nstart
  DO i = 1, myid
     newnstart = newnstart - numtempproc(i)
  ENDDO
  print*,myid,' number of disp bcs on other procs = ',numtempproc(:)

  ! ... Have each proc in turn broadcast the node numbers that it owns which have temperature 
  ! ... boundary conditions.  The global list of these node numbers will be stored in the array
  ! ... GTempBC.
  ALLOCATE(gtempbc(1:gnumtemp))
  gtempbc(:) = 0
  counter1 = 0
  ! ... Loop through the processes
  DO i = 1, nprocs
     ! ... If there are temperature boundary conditions on process i then allocate the temporary array
     ! ... tempintv to the number of temperature BCs on process i.
     IF (numtempproc(i) > 0) THEN
        ALLOCATE(tempintv(1:numtempproc(i)))
        tempintv(:) = 0
        ! ... It process i is this process, fill the array tempintv with the node numbers of the 
        ! ... temperature boundary conditions on this process.
        if (i-1 == myid) tempintv(:) = dispbc(:)
        ! ... The process i broadcasts the array tempintv to all other processes.
        CALL mpi_bcast(tempintv(1), numtempproc(i), mpi_integer, i-1, rocstar_communicator, ierr)
        ! ... Loop through the number of temperature BC nodes on process i and record them on 
        ! ... array GTempBC.
        DO j = 1, numtempproc(i)
           counter1 = counter1 + 1
           gtempbc(counter1) = tempintv(j)
        ENDDO
        DEALLOCATE(tempintv)
     ENDIF
  ENDDO
  print*,myid,' global disp bcs = ',gtempbc(:)

  ! ... convert Ceff to explicit form, CeffExplicit

  ! ... Allocate explicit matrix
  ALLOCATE(ceffexplicit(nrows,ndim))

  ! ... initialize CeffExplicit
  ceffexplicit = 0.0d0

  ! ... n is the CRS index for the aval and cval vectors
  n = 1

  ! ... First loop through the rows i
  DO i = 1, nrows
     ! ... loop through the columns j
     DO j = 1, ndim
        ! ... is n in the index range for row i?
        IF ((rp1(i) < n).AND.(rp1(i+1) >= n)) THEN
           ! ... if this is the correct column corresponding to index n
           IF (cval(n) == j - 1) THEN
              ceffexplicit(i,j) = aval(n)
              n = n + 1
           ENDIF
        ENDIF
        ! ... if end of CRS vectors, no need to keep looping
        IF (n > nnz) EXIT
     ENDDO
     ! ... if end of CRS vectors, no need to keep looping
     IF (n > nnz) EXIT
  ENDDO



  ! ... Partition matrix CeffExplicit into explicit forms of
  ! ... - CeffFF: Free-free DOFs
  ! ... - CeffPP: Prescribed-prescribed DOFs
  ! ... - CeffPF: Prescribed-free DOFs
  ! ... - CeffFP: Free-prescribed DOFs

  ! ... Allocate the partitioned matrices
  ALLOCATE(ceffff(nrows - numdisp, ndim - gnumtemp))
  ALLOCATE(ceffpp(numdisp,gnumtemp))
  ALLOCATE(ceffpf(numdisp, ndim - gnumtemp))
  ALLOCATE(cefffp(nrows - numdisp, gnumtemp))
  ALLOCATE(ceff_fptp(nrows - numdisp))

  ! ... Begin partitioning
  counter1 = 0
  counter2 = 0

  ! ... outer loop through rows i
  DO i = nstart, nstart + nrows - 1
     ! ... iDOF stores the DOF number of the current essential BC
     idof = 0

     ! ... reset the nmber of prescribed BSs counted in a row (counter 2)
     counter2 = 0

     ! ... check the DOFs of all local prescribed BCs to see if they are in
     ! ... this row
     DO n = 1, numdisp

        IF(i == dispbc(n)) THEN
           ! ... if i is a prescribed DOF, then increase row BC counter1
           counter1 = counter1 + 1


           ! ... if i is a prescribed DOF, then record the value of the DOF
           idof = i

           ! ... stop searching through dispBC, i can only equal one of them
           EXIT
        ENDIF
     ENDDO
     jdof = 0
     ! ... if row i corresponds to a prescribed DOF, then fill CeffPP and CeffPF
     IF(idof > 0) THEN
        ! ... loop through the columns in prescribed row i
        DO j = 1, ndim
           ! ... check the DOFs of all local prescribed BCs to see if they are in
           ! ... this column, j
           DO n = 1, gnumtemp
              IF(j == gtempbc(n)) THEN
                 ! ... if i is a prescribed DOF, then increase column BC counter2
                 counter2 = counter2 + 1
                 ! ... if i is a prescribed DOF, then record the value of the dof
                 jdof = j
                 ! ... stop searching through dispBC, j can only equal one of them
                 EXIT
              ENDIF
           ENDDO
           IF(jdof > 0) THEN
              ! ... if row j corresponds to a prescribed DOF, then fill CeffPP

              ceffpp(counter1,counter2) = ceffexplicit(idof-nstart+1,jdof)

              ! ... reset jDOF for next loop
              jdof = 0
           ELSE

              ! ... if row j corresponds to a free DOF, then fill CeffPF
              ceffpf(counter1,j-counter2) = ceffexplicit(idof-nstart+1,j)
           ENDIF
        ENDDO
     ELSE
        ! ... loop through the columns in free row i
        DO j = 1, ndim

           ! ... check the DOFs of all local prescribed BCs to see if they are in
           ! ... this column, j
           DO n = 1, gnumtemp

              IF(j == gtempbc(n)) THEN
                 ! ... if i is a prescribed DOF, then increase column BC counter2
                 counter2 = counter2 + 1
                 ! ... if i is a prescribed DOF, then record the value of the dof
                 jdof = j
                 ! ... stop searching through dispBC, j can only equal one of them
                 EXIT
              ENDIF
           ENDDO


           IF(jdof > 0) THEN
              ! ... if row j corresponds to a prescribed DOF
              cefffp(i-counter1-nstart+1,counter2) = ceffexplicit(i-nstart+1,jdof)

              ! ... reset jDOF for next loop
              jdof = 0
           ELSE
              ! ... if row j corresponds to a free DOF, then fill CeffFF
              ceffff(i-counter1-nstart+1,j-counter2) = ceffexplicit(i-nstart+1,j)
           ENDIF
        ENDDO
     ENDIF
  ENDDO

  ! ... Now convert CeffFF back to CRS format

  ! ... Count nonzeros in new matrix
  nnz_temp = 0
  counter2 = 0
  ! ... Loop throught the number of rows on this process
  DO i = 1, nrows
     ! ... Loop from the beginning to the end of each row, keeping track for the column
     ! ... index using counter2.
     DO j = rp1(i)+1, rp1(i+1)
        counter2 = counter2 + 1
        counter1 = 0
        DO m = 1, gnumtemp
           ! ... If the global index of a temperature BC matches either the row or column
           ! ... indices, then do not increment the number of non-zeros.
           IF((gtempbc(m) == i + nstart - 1).OR.(gtempbc(m) == cval(counter2)+1)) THEN
              counter1 = 1
           ENDIF
        ENDDO

        IF(counter1 == 0) nnz_temp = nnz_temp + 1
     ENDDO
  ENDDO

  ! ... Count rows in new matrix
  newnrows = 0
  counter2 = 0
  ! ... Loop throught the number of rows on this process
  DO i = 1, nrows
     counter1 = 0
     DO m = 1, gnumtemp
        ! ... If the global index of a termperature BC matches the row index, then
        ! ... do not increment the new number of rows on this process.
        IF(gtempbc(m) == i + nstart - 1) THEN
           counter1 = 1
        ENDIF
     ENDDO
     if(counter1==1) print*,myid,' row removed at ',i+nstart-1
     IF(counter1 == 0) newnrows = newnrows + 1
  ENDDO

  ! ... Allocate variables for new matrix
  ALLOCATE(rp_temp(1:newnrows+1))
  ALLOCATE(cval_temp(1:nnz_temp))
  ALLOCATE(aval_temp(1:nnz_temp))

  ! ... n is the CRS index for the vectors aval_temp and cval_temp
  n = 0
  ! ... loop through the rows in the matrix CeffFF
  DO i = 1, newnrows
     rp_temp(i) = n
     ! ... loop through the columns in the matrix CeffFF
     DO j = 1, newndim
        ! ... check for a non-zero at i,j in CeffFF
        IF (ceffff(i,j) > 0.0d0) THEN
           n = n + 1
           ! ... if there is a non-zero, it is entered into position n of 
           ! ... aval_temp and the column j(-1) is entered into positon n 
           ! ... of cval_temp
           aval_temp(n) = ceffff(i,j)
           cval_temp(n) = j - 1
        ENDIF
     ENDDO
  ENDDO
  ! ... fill in the last entry of rp_temp

  rp_temp(newnrows+1) = nnz_temp



  ! ... multiply CeffFP with the global prescribed temperature vector
  CALL prescribedload(nrows,nstart,numdisp,dispbc,cefffp,proctemp)



END SUBROUTINE partition_ceff


!!****
!!
!!  NAME
!!     PrescribedLoad
!!
!!  FUNCTION
!!     This subroutine multiplies the Free-Prescribed partition of the effective thermal capacitance
!!     matrix with the list of global prescribed temperatures.
!!
!!  INPUTS
!!     ncols -- The number of columns in Ceff_fp, same as number of prescribed temperatures
!!     nrows -- The number of rows in Ceff_fp
!!     nstart -- The global index of the first row assigned to this processor
!!     GNumTemp -- The global number of prescribed temperatures
!!     GTempBC -- Vector containing global nodes of the prescribed temperatures 
!!     NumpTemp -- The local number of prescribed temperatures
!!     TempBC -- Vector containing local nodes of the prescribed temperatures
!!     Ceff_fp -- The free-perscribed partition of the effective thermal capacitance matrix
!!     ProcTemp -- The part of the global temperature vector that's assigned to this processor 
!!
!!  OUTPUTS
!!     Ceff_fpTp -- The product of Ceff_fp and the global prescribed temperature vector
!!
!!  USES
!!     MPI
!!
!!****

SUBROUTINE prescribedload(nrows,nstart,NumTemp,TempBC,Ceff_fp,ProcTemp)

  USE precision
  USE implicit_global
  USE comp_row_global
  USE rocstar_rocfraccomm
  USE removebcs_global

  IMPLICIT NONE
  
  include 'mpif.h'
  
  ! ... Input variables
  ! ... the Ceff_fp matrix
  INTEGER :: nrows, numtemp
  REAL(kind=wp), DIMENSION(nrows-NumTemp,GNumTemp) :: ceff_fp
  ! ... the local temperature boundary condions
  INTEGER, DIMENSION(NumTemp) :: tempbc
  ! ... first row of this processor
  INTEGER :: nstart
  ! ... all temperatures on this processor
  REAL(kind=wp), DIMENSION(LNumNp) :: proctemp

  ! ... Internal variables
  ! ... the prescribed temperature vector
  REAL(kind=wp), DIMENSION(GNumTemp) :: gtempbcvalue

  ! ... communication
  REAL(kind=wp), DIMENSION(GNumNp) :: temptemp
  INTEGER, ALLOCATABLE, DIMENSION(:) :: numtempfrom
  REAL(kind=wp), ALLOCATABLE, DIMENSION(:) :: bufsnd
  INTEGER :: i, j, k, m, counter1, counter2
  REAL(kind=wp) :: tempval
  INTEGER, ALLOCATABLE, DIMENSION(:) :: req_rcv, req_snd
  INTEGER, ALLOCATABLE, DIMENSION(:,:) :: stat_rcv, stat_snd

  !
  ! Communicate how many pieces of temp to send to each other proc
  !

  IF ( nprocs > 1) THEN

     CALL mpi_barrier(rocstar_communicator,ierr)
     ALLOCATE(req_rcv(1:nprocs))
     ALLOCATE(req_snd(1:nprocs))
     ALLOCATE(stat_snd(1:mpi_status_size,1:nprocs))
     ALLOCATE(stat_rcv(1:mpi_status_size,1:nprocs))
     ALLOCATE(numtempfrom(1:nprocs))
     DO i = 1, nprocs
        req_rcv(i) = 0
        req_snd(i) = 0
        DO j = 1, mpi_status_size
           stat_rcv(j,i) = 0
           stat_snd(j,i) = 0
        ENDDO
     ENDDO
     DO i = 1, nprocs
        IF (i-1 /= myid) THEN
           CALL mpi_irecv(numtempfrom(i),1, &
                mpi_integer,i-1,10,rocstar_communicator, &
                req_rcv(i),ierr)
        ENDIF
     ENDDO
     DO i = 1, nprocs
        IF (i-1 /= myid) THEN
           !CALL MPI_ISEND(nrows,1,MPI_INTEGER, &
           !   i-1,10,ROCSTAR_COMMUNICATOR,req_snd(i),ierr)
           CALL mpi_send(nrows,1,mpi_integer, &
                i-1,10,rocstar_communicator,ierr)
        ENDIF
     ENDDO
     CALL mpi_waitall(nprocs,req_rcv,stat_rcv,ierr)
     !CALL MPI_WAITALL(nprocs,req_snd,stat_snd,ierr)
     DEALLOCATE(req_rcv)
     DEALLOCATE(req_snd)
     DEALLOCATE(stat_snd)
     DEALLOCATE(stat_rcv)


     !
     ! Communicate temperature to all the other procs
     !
     
     ! ... receive
     CALL mpi_barrier(rocstar_communicator,ierr)
     ALLOCATE(req_rcv(1:nprocs))
     ALLOCATE(req_snd(1:nprocs))
     ALLOCATE(stat_snd(1:mpi_status_size,1:nprocs))
     ALLOCATE(stat_rcv(1:mpi_status_size,1:nprocs))
     DO i = 1, nprocs
        req_rcv(i) = 0
        req_snd(i) = 0
        DO j = 1, mpi_status_size
           stat_rcv(j,i) = 0
           stat_snd(j,i) = 0
        ENDDO
     ENDDO
     ALLOCATE(frmproc(1:nprocs))
     DO i = 1, nprocs
        ! ... if processor i is not this processor (i=1=processor 0)
        IF (i-1 /= myid) THEN
         ! ... Allocate a buffer of size 2 times the number of nodes from
         ! ... processor j.  The buffer is orgnized as {node #, temp, node #, temp...}
           ALLOCATE(frmproc(i)%rcvbuf(1:2*int(numtempfrom(i))))
           CALL mpi_irecv(frmproc(i)%rcvbuf(1),2*int(numtempfrom(i)), &
                mpi_double_precision,i-1,10,rocstar_communicator, &
                req_rcv(i),ierr)
        ENDIF
     ENDDO

     ! ... send
     DO i = 1, nprocs
        IF (i-1 /= myid) THEN
           ! ... if processor i is not this processor
           ! ... allocate a buffer 2 times the number of the nodes 
           ! ... on this processor to send to processor i
           ALLOCATE(bufsnd(1:2*lnumnp))
           counter1 = 0
           DO j = 1, gnumnp
              ! ... if global node j belongs to this processor
              ! ... then the corresponding temperature is copied
              ! ... to a buffer to go to processor i
              IF ((global2local(j) /= -1) .AND. (nodeproc(global2local(j)) == myid)) THEN
                 ! ... if global node j belongs to this processor
                 ! ... then counter1 determines the index in ProcTemp where
                 ! ... the corresponding temperature is stored
                 counter1 = counter1 + 1
                 bufsnd(2*counter1-1) = j
                 bufsnd(2*counter1)   = proctemp(counter1)
              ENDIF
           ENDDO
           !CALL MPI_ISEND(bufsnd,4*INT(nrows/3),MPI_DOUBLE_PRECISION, &
           !     i-1,10,ROCSTAR_COMMUNICATOR,req_snd(i),ierr)

           ! ... send the buffer of length 2*nodes on this processor
           ! ... to processor (i-1).
           CALL mpi_send(bufsnd,2*nrows,mpi_double_precision, &
                i-1,10,rocstar_communicator,ierr)
           DEALLOCATE(bufsnd)
        ENDIF
     ENDDO
     !CALL MPI_WAITALL(NumCommProcsFrom2,req_rcv,stat_rcv,ierr)
     !CALL MPI_WAITALL(NumCommProcs2,req_snd,stat_snd,ierr)
     CALL mpi_waitall(nprocs,req_rcv,stat_rcv,ierr)
     !CALL MPI_WAITALL(nprocs,req_snd,stat_snd,ierr)
     DEALLOCATE(req_rcv)
     DEALLOCATE(req_snd)
     DEALLOCATE(stat_snd)
     DEALLOCATE(stat_rcv)
     
  ENDIF

  ! ... create a global temperature vector
  DO i = 1, gnumnp
     ! ... if global node i belongs to this processor
     ! ... then the corresponding temperature is copied
     ! ... in from the input array ProcTemp
     IF ((global2local(i) /= -1) .AND. (nodeproc(global2local(i)) == myid)) THEN
        temptemp( i ) = proctemp( i - (nstart-1))
     ELSE
     
        ! ... if global node i does not belong to this processor
        ! ... then find which processor it does belong to
        ! ... and get the temperature from the buffer received from
        ! ... that processor
        DO j = 1, nprocs
           IF (j-1 /= myid) THEN
              ! ... loop through the number of temperature values from processor j
              DO m = 1, numtempfrom(j)
                 ! ... frmproc(j)%rcvbuf(2*m-1) is where the node number
                 ! ... is stored for the temperature value stored at 
                 ! ... frmproc(j)%rcvbuf(2*m)
                 IF (int(frmproc(j)%rcvbuf(2*m-1)) == i) THEN
                    temptemp( i )   = frmproc(j)%rcvbuf(2*m)
                 ENDIF
              ENDDO
           ENDIF
        ENDDO
     ENDIF
  ENDDO
  IF (nprocs > 1) DEALLOCATE(frmproc)


  ! ... Get only the prescribed temperatures
  ! ... put in GTempBCValues
  DO i = 1, gnumnp
     ! ... loop through prescribed temperature nodes
     DO j = 1, gnumtemp
        IF (i == gtempbc(j)) THEN
           gtempbcvalue(j)= temptemp(i)
        ENDIF
     ENDDO
  ENDDO

ceff_fptp = 0.0d0

  ! ... multiply Ceff_fp with GTempBCValue
  DO i = 1, nrows-numtemp
     DO j = 1, gnumtemp
        ceff_fptp(i) = ceff_fptp(i) + ceff_fp(i,j) * gtempbcvalue(j)
     ENDDO
  ENDDO

  
END SUBROUTINE prescribedload