solver.f90

Go to the documentation of this file.

MODULE solve_petsc
use my_util
  TYPE solver_param
     INTEGER:: it_max
     REAL(KIND=8) :: rel_tol, abs_tol
     CHARACTER(LEN=20) :: solver, precond
     LOGICAL           :: verbose
  END TYPE solver_param

CONTAINS
  SUBROUTINE init_solver(my_par,my_ksp,matrix,communicator,solver,precond, opt_re_init)
    USE chaine_caractere
    IMPLICIT NONE
    LOGICAl, INTENT(IN),    OPTIONAL :: opt_re_init
    type(solver_param)               :: my_par
    CHARACTER(*),           OPTIONAL :: solver, precond
    LOGICAL                          :: re_init
    INTEGER  :: deb, fin
#include "petsc/finclude/petsc.h"
    mat            :: matrix
    ksp            :: my_ksp
    pc             :: prec
    petscerrorcode :: ierr
    mpi_comm       :: communicator

    IF (.NOT.present(opt_re_init)) THEN
       re_init=.false.
    ELSE
       re_init=opt_re_init
    END IF

    IF (my_par%it_max.LE.0) THEN
       my_par%it_max = 100
    END IF
    IF (my_par%rel_tol.LE.0.d0) THEN
       my_par%rel_tol = 1.d-8
    END IF
    IF (my_par%abs_tol.LE.0.d0) THEN
       my_par%abs_tol = 1.d-14
    END IF

    IF (.NOT.re_init) CALL kspcreate(communicator,my_ksp,ierr)
    !CALL KSPCreate(communicator,my_ksp,ierr)

    !CALL KSPSetOperators(my_ksp,matrix,matrix,DIFFERENT_NONZERO_PATTERN,ierr)
    !CALL KSPSetOperators(my_ksp,matrix,matrix,SAME_NONZERO_PATTERN ,ierr) !Petsc 3.4.2
    CALL kspsetoperators(my_ksp,matrix,matrix,ierr) !Petsc 3.7.2

    IF (present(solver)) THEN
       deb = start_of_string(solver)
       fin = last_of_string(solver)
       IF (solver(deb:fin)=='BCGS') THEN   
          CALL kspsettype(my_ksp, kspbcgs, ierr)
       ELSE IF (solver(deb:fin)=='GMRES') THEN 
          CALL kspsettype(my_ksp, kspgmres, ierr)
       ELSE IF (solver(deb:fin)=='FGMRES') THEN 
          CALL kspsettype(my_ksp, kspfgmres, ierr)
       ELSE IF (solver(deb:fin)=='PCR') THEN 
          CALL kspsettype(my_ksp, kspcr, ierr)
       ELSE IF (solver(deb:fin)=='CHEBYCHEV') THEN 
          CALL kspsettype(my_ksp, kspchebyshev, ierr)
       ELSE IF (solver(deb:fin)=='CG') THEN
          CALL kspsettype(my_ksp, kspcg, ierr)
       ELSE
          CALL kspsettype(my_ksp, kspfgmres, ierr)
       END IF
    ELSE
       CALL kspsettype(my_ksp, kspfgmres, ierr)
    END IF
    CALL kspsettolerances(my_ksp, my_par%rel_tol, my_par%abs_tol, &
         petsc_default_real, my_par%it_max, ierr)
    CALL kspgetpc(my_ksp, prec, ierr)
    IF (present(precond)) THEN
       deb = start_of_string(precond)
       fin = last_of_string(precond)
       IF (precond(deb:fin)=='JACOBI') THEN
          call pcsettype(prec, pcbjacobi, ierr)
       ELSE IF  (precond(deb:fin)=='HYPRE') THEN
          CALL pcsettype(prec, pchypre, ierr)
       ELSE IF  (precond(deb:fin)=='SSOR') THEN
          CALL pcsettype(prec, pcsor, ierr)
       ELSE IF  (precond(deb:fin)=='MUMPS') THEN
          CALL pcsettype(prec, pclu, ierr)
          CALL kspsettype(my_ksp, ksppreonly, ierr)
          CALL pcfactorsetmatsolverpackage(prec, matsolvermumps, ierr)
       ELSE
          CALL pcsettype(prec, pchypre, ierr)
       END IF
    ELSE
       CALL pcsettype(prec, pchypre, ierr)
    END IF
    CALL kspsetfromoptions(my_ksp, ierr)
  END SUBROUTINE init_solver

  SUBROUTINE solver(my_ksp,b,x,reinit,verbose)
    IMPLICIT NONE
    LOGICAL,  OPTIONAL :: reinit, verbose
    INTEGER            :: its
#include "petsc/finclude/petsc.h"
    ksp            :: my_ksp
    petscerrorcode :: ierr
    vec            :: x, b
    kspconvergedreason :: reason
    IF (.NOT.present(reinit)) reinit=.true.
    IF (.NOT.present(verbose)) verbose=.false.

    IF (reinit) CALL veczeroentries(x,ierr)
    CALL kspsolve(my_ksp,b,x,ierr)
    IF (verbose) THEN
       CALL kspgetiterationnumber(my_ksp, its, ierr)
       CALL kspgetconvergedreason(my_ksp, reason, ierr)
       SELECT case(reason)
       case(2) 
          WRITE(*,*) "KSP_CONVERGED_RTOL, Nb of iterations", its
       case(3) 
          WRITE(*,*) "KSP_CONVERGED_ATOL, Nb of iterations", its
       case(4) 
          WRITE(*,*) "Converged after one single iteration of the preconditioner is applied"
       case(5,6,7,8) 
          WRITE(*,*) "Converge for strange reason:", reason
       case(-2) 
          WRITE(*,*) "KSP_DIVERGED_NULL" 
       case(-3) 
          WRITE(*,*) "Not converged after it_max", its
       case(-4) 
          WRITE(*,*) "Not converged: explosion"
       case(-5,-6,-7) 
          WRITE(*,*) "Not converged for strange reasons", reason
       case(-8) 
          WRITE(*,*) "Not converged: Indefinite preconditioner"
       case(-9) 
          WRITE(*,*) "Not converged: NAN"
       case(-10) 
          WRITE(*,*) "Not converged: Indefinite matrix"
       case default
          WRITE(*,*) "Something strange happened", reason
       END SELECT
    END IF

  END SUBROUTINE solver

  SUBROUTINE create_local_petsc_matrix(communicator, LA, matrix, clean)
    USE def_type_mesh
    IMPLICIT NONE
    TYPE(petsc_csr_la)                 :: la 
    LOGICAL, OPTIONAL                  :: clean
    REAL(KIND=8),DIMENSION(:), POINTER :: aa
    INTEGER :: nnzm1, dom_np
    LOGICAL :: test_clean 

!!$  INTEGER, DIMENSION(:), POINTER :: ia, ja

#include "petsc/finclude/petsc.h"
    mpi_comm       :: communicator
    mat            :: matrix
    petscerrorcode :: ierr
    !------------------------------------------------------------------------------
    dom_np = SIZE(la%ia)-1
    nnzm1=la%ia(dom_np)-la%ia(0)-1
    ALLOCATE(aa(0:nnzm1))
    aa =0.d0


!!$  ALLOCATE(ia(0:dom_np),ja(0:nnzm1))
!!$  ia = LA%ia
!!$  ja = LA%ja
!!$  CALL MatCreateMPIAIJWithArrays(communicator,dom_np,dom_np,PETSC_DECIDE, &
!!$       PETSC_DECIDE, ia, ja, aa, matrix, ierr)
!!$DEALLOCATE(ia,ja)

    CALL matcreatempiaijwitharrays(communicator,dom_np,dom_np,petsc_decide, &
         petsc_decide, la%ia, la%ja, aa, matrix, ierr)

    DEALLOCATE(aa)
    IF (present(clean)) THEN
       test_clean=clean
    ELSE
       test_clean=.true.
    END IF
    IF (test_clean) THEN
       IF (ASSOCIATED(la%ia)) DEALLOCATE(la%ia)
       IF (ASSOCIATED(la%ja)) DEALLOCATE(la%ja)
    END IF
  END SUBROUTINE create_local_petsc_matrix

  SUBROUTINE create_local_petsc_matrix_a_detruire(communicator, aij, i_loc, matrix)
    USE def_type_mesh
    IMPLICIT NONE
    TYPE(aij_type),         INTENT(IN) :: aij
    INTEGER,     DIMENSION(2)          :: i_loc
    INTEGER,     DIMENSION(:), POINTER :: ia, ja
    REAL(KIND=8),DIMENSION(:), POINTER :: aa
    INTEGER :: nnzm1, dom_np, p, i, n
#include "petsc/finclude/petsc.h"
    mpi_comm       :: communicator
    mat            :: matrix
    petscerrorcode :: ierr
    !------------------------------------------------------------------------------
    dom_np = i_loc(2) - i_loc(1) + 1
    nnzm1 = aij%ia(i_loc(2)+1)-aij%ia(i_loc(1))-1
    ALLOCATE(ia(0:dom_np),ja(0:nnzm1))
    ia(0)=0
    DO i = 1, dom_np
       n = i_loc(1) + i -1
       ia(i) = aij%ia(n+1)-aij%ia(i_loc(1))
       DO p=aij%ia(n), aij%ia(n+1)-1
          ja(p-aij%ia(i_loc(1)))= aij%ja(p)-1
       end do
    END DO
    !------------------------------------------------------------------------------
    ALLOCATE(aa(0:nnzm1))
    aa =0
    CALL matcreatempiaijwitharrays(communicator,dom_np,dom_np,petsc_decide, &
         petsc_decide, ia, ja, aa, matrix, ierr)

    DEALLOCATE(ia,ja,aa)
  END SUBROUTINE create_local_petsc_matrix_a_detruire

  SUBROUTINE create_local_petsc_block_matrix(communicator, n_b, aij, i_loc, matrix)
    USE def_type_mesh
    IMPLICIT NONE
    TYPE(aij_type),           INTENT(IN) :: aij
    INTEGER,     DIMENSION(2)            :: i_loc
    INTEGER                              :: n_b
    INTEGER,     DIMENSION(:),   POINTER :: ia, ja
    REAL(KIND=8),DIMENSION(:),   POINTER :: aa
    INTEGER :: nnzm1, dom_np, p, i, n, ib, k
#include "petsc/finclude/petsc.h"
    mpi_comm       :: communicator
    mat            :: matrix
    petscerrorcode :: ierr

    dom_np = i_loc(2) - i_loc(1) + 1
    nnzm1 = n_b*(aij%ia(i_loc(2)+1)-aij%ia(i_loc(1))-1)
    ALLOCATE(ia(0:n_b*dom_np),ja(0:nnzm1))
    ia(0)=0
    DO k = 1, n_b
       DO i = 1, dom_np
          ib = i + (k-1)*dom_np
          n = i_loc(1) + i - 1
          ia(i) = n_b*(aij%ia(n+1)-aij%ia(i_loc(1)))
          DO p=aij%ia(n), aij%ia(n+1)-1
             ja(p-aij%ia(i_loc(1)))= aij%ja(p)-1
          end do
       END DO
    END DO
    !------------------------------------------------------------------------------
    ALLOCATE(aa(0:nnzm1))
    aa =0
    CALL matcreatempiaijwitharrays(communicator,dom_np,dom_np,petsc_decide, &
         petsc_decide, ia, ja, aa, matrix, ierr)

    DEALLOCATE(ia,ja,aa)
  END SUBROUTINE create_local_petsc_block_matrix


END MODULE solve_petsc