Merge pull request #13826 from marcosvanella/master

FDS Source : Add HYPRE PG+AMG to UGLMAT.

Author: marcosvanella

Source/cons.f90

@@ -517,8 +517,10 @@ MODULE GLOBAL_CONSTANTS
 INTEGER, PARAMETER :: UGLMAT_FLAG=2                              !< Integer pressure solver parameter UGLMAT
 INTEGER, PARAMETER :: ULMAT_FLAG=3                               !< Integer pressure solver parameter ULMAT
 INTEGER, PARAMETER :: MKL_PARDISO_FLAG=1                         !< Integer matrix solver library flag for MKL PARDISO
+INTEGER, PARAMETER :: MKL_CPARDISO_FLAG=1                        !< Integer matrix solver library flag for MKL CLUSTER PARDISO
 INTEGER, PARAMETER :: HYPRE_FLAG=2                               !< Integer matrix solver library flag for HYPRE
 INTEGER :: ULMAT_SOLVER_LIBRARY=MKL_PARDISO_FLAG                 !< Integer ULMAT library flag (defaults to MKL PARDISO)
+INTEGER :: UGLMAT_SOLVER_LIBRARY=MKL_CPARDISO_FLAG               !< Integer UGLMAT library flag (defaults to MKL CPARDISO)
 INTEGER :: PRES_FLAG = FFT_FLAG                                  !< Pressure solver
 LOGICAL :: TUNNEL_PRECONDITIONER=.FALSE.                         !< Use special pressure preconditioner for tunnels
 INTEGER :: TUNNEL_NXP                                            !< Number of x points in the entire tunnel

Source/imkl.f90

@@ -203,7 +203,7 @@ MODULE HYPRE_INTERFACE
    INTEGER :: HYPRE_IERR = 0
    INTEGER, PARAMETER :: HYPRE_SOLVER_ID = 1                 ! Preconditioned Conjugate Gradient (PCG) solver
    INTEGER, PARAMETER :: HYPRE_PRECOND_ID = 2                ! Algebraic Multi-Grid (AMG) preconditioner
-   INTEGER, PARAMETER :: HYPRE_SOLVER_MAXIT = 1000           ! Max iteratations of PCG solver
+   INTEGER, PARAMETER :: HYPRE_SOLVER_MAXIT = 1000           ! Max iterations of PCG solver
    REAL(KIND=8), PARAMETER :: HYPRE_SOLVER_TOL = 1.D-12      ! Solver tolerance
    INTEGER, PARAMETER :: HYPRE_SOLVER_SETTWONORM = 1         ! 0=use L_Infty norm (max error) for convergence, 1=use L2 norm
    INTEGER, PARAMETER :: HYPRE_SOLVER_SETPRINTLEVEL = 0      ! 0=no output, 1=minimal, 2=verbose

Source/pres.f90

@@ -3051,6 +3051,9 @@ SUBROUTINE GLMAT_SOLVER(T,DT)
 USE COMP_FUNCTIONS, ONLY: CURRENT_TIME
 USE CC_SCALARS, ONLY : GET_CUTCELL_HP,GET_PRES_CFACE_BCS,GET_FH_FROM_PRHS_AND_BCS
 USE MPI_F08
+#ifdef WITH_HYPRE
+USE HYPRE_INTERFACE
+#endif
 
 REAL(EB), INTENT(IN) :: T,DT
 
@@ -3150,6 +3153,8 @@ SUBROUTINE GLMAT_SOLVER(T,DT)
    ENDIF
    ! WRITE(LU_ERR,*) 'SUM_FH=',SUM(F_H),H_MATRIX_INDEFINITE
 
+   LIBRARY_SELECT: SELECT CASE(UGLMAT_SOLVER_LIBRARY)
+   CASE(MKL_CPARDISO_FLAG) LIBRARY_SELECT
 #ifdef WITH_MKL
    ! Dump local low an high rows assembled by this process in IPARM:
    IPARM(41) = ZSL%LOWER_ROW
@@ -3169,6 +3174,22 @@ SUBROUTINE GLMAT_SOLVER(T,DT)
 IF (ERROR /= 0 .AND. MY_RANK==0) WRITE(LU_ERR,*) 'GLMAT_SOLVER: The following ERROR was detected: ', ERROR
 #endif
 
+   CASE(HYPRE_FLAG) LIBRARY_SELECT
+#ifdef WITH_HYPRE
+   ! Solve the system using HYPRE:
+   CALL HYPRE_IJVECTORSETVALUES(ZSL%HYPRE_ZSL%F_H, ZSL%NUNKH_LOCAL, ZSL%HYPRE_ZSL%INDICES, ZSL%F_H, HYPRE_IERR)
+   CALL HYPRE_IJVECTORASSEMBLE(ZSL%HYPRE_ZSL%F_H, HYPRE_IERR)
+   CALL HYPRE_PARCSRPCGSOLVE(ZSL%HYPRE_ZSL%SOLVER, ZSL%HYPRE_ZSL%PARCSR_A_H, ZSL%HYPRE_ZSL%PAR_F_H, &
+                             ZSL%HYPRE_ZSL%PAR_X_H, HYPRE_IERR)
+   IF (CHECK_POISSON .AND. HYPRE_SOLVER_SETPRINTLEVEL>0) THEN
+      CALL HYPRE_PARCSRPCGGETNUMITERATIONS(ZSL%HYPRE_ZSL%SOLVER, ZSL%HYPRE_ZSL%NUM_ITERATIONS, HYPRE_IERR)
+      CALL HYPRE_PARCSRPCGGETFINALRELATIVE(ZSL%HYPRE_ZSL%SOLVER, ZSL%HYPRE_ZSL%FINAL_RES_NORM, HYPRE_IERR)
+   ENDIF
+   CALL HYPRE_IJVECTORGETVALUES(ZSL%HYPRE_ZSL%X_H, ZSL%NUNKH_LOCAL, ZSL%HYPRE_ZSL%INDICES, ZSL%X_H, HYPRE_IERR)
+#endif
+
+   END SELECT LIBRARY_SELECT
+
    IF (ZSL%MTYPE==-2) THEN
       SUM_XH = 0._EB; MEAN_XH = 0._EB
       WHOLE_DOM_IF2 : IF(.NOT.PRES_ON_WHOLE_DOMAIN) THEN
@@ -3321,6 +3342,18 @@ SUBROUTINE GLMAT_SOLVER_SETUP(STAGE_FLAG)
 IF (FREEZE_VELOCITY)  RETURN ! Fixed velocity soln. i.e. PERIODIC_TEST=102 => FREEZE_VELOCITY=.TRUE.
 IF (SOLID_PHASE_ONLY) RETURN
 TNOW=CURRENT_TIME()
+
+! If either MKL or HYPRE library not present stop.
+#ifndef WITH_MKL
+#ifndef WITH_HYPRE
+IF (MY_RANK==0) WRITE(LU_ERR,'(A)') &
+'Error: MKL or HYPRE Library compile flag not defined for UGLMAT pressure solver.'
+! Some error - stop flag for CALL STOP_CHECK(1).
+STOP_STATUS = SETUP_STOP
+RETURN
+#endif
+#endif
+
 SELECT CASE(STAGE_FLAG)
 CASE(1)
 
@@ -3849,9 +3882,13 @@ END SUBROUTINE COPY_CCVAR_IN_HS
 ! ------------------------------- GET_H_MATRIX_LUDCMP -------------------------------
 
 SUBROUTINE GET_H_MATRIX_LUDCMP
-#ifdef WITH_MKL
+#if defined WITH_MKL || defined WITH_HYPRE
 USE MPI_F08
 #endif
+#ifdef WITH_HYPRE
+USE HYPRE_INTERFACE
+#endif
+
 ! Local Variables:
 INTEGER :: INNZ, IROW, JCOL
 #ifdef WITH_MKL
@@ -3864,8 +3901,13 @@ SUBROUTINE GET_H_MATRIX_LUDCMP
 INTEGER, ALLOCATABLE, DIMENSION(:,:) :: MB_FACTOR
 INTEGER :: IERR
 #endif
+#ifdef WITH_HYPRE
+INTEGER :: ONEV(1)
+REAL(EB) :: ZEROV(1)
+#endif
 
 ! Define parameters:
+INNZ = 0; IROW = 0; JCOL=0
 NRHS   = 1
 MAXFCT = 1
 MNUM   = 1
@@ -3875,9 +3917,21 @@ SUBROUTINE GET_H_MATRIX_LUDCMP
 
 ERROR     = 0 ! initialize error flag
 
+SELECT CASE(UGLMAT_SOLVER_LIBRARY)
+CASE(MKL_CPARDISO_FLAG)
 #ifdef WITH_MKL
-CALL SET_CLUSTER_SOLVER_IPARM
+   CALL SET_CLUSTER_SOLVER_IPARM
 #endif
+CASE(HYPRE_FLAG)
+#ifdef WITH_HYPRE
+   IERR=0; CALL HYPRE_INITIALIZE(IERR)
+   IF (IERR==1) THEN
+      IF(MY_RANK==0) WRITE(LU_ERR,'(A)') 'Error: HYPRE pressure solver initialization error.'
+      STOP_STATUS = SETUP_STOP
+      RETURN
+   ENDIF
+#endif
+END SELECT
 
 IPZ_LOOP : DO IPZ=0,N_ZONE
 
@@ -3891,16 +3945,25 @@ SUBROUTINE GET_H_MATRIX_LUDCMP
    ! Allocate F_H ans H_H for this Process and IPZ:
    ALLOCATE( ZSL%X_H(MAX(ZSL%NUNKH_LOCAL,1)) , ZSL%F_H(MAX(ZSL%NUNKH_LOCAL,1)) ); ZSL%F_H=0._EB; ZSL%X_H=0._EB
 
+   ! Here each process defines de beginning and end rows in global numeration, for the equations
+   ! it has assembled:
+   IF(ZSL%NUNKH_LOCAL>0) THEN
+      ZSL%LOWER_ROW = ZSL%UNKH_IND(NM_START) + 1
+      ZSL%UPPER_ROW = ZSL%UNKH_IND(NM_START) + ZSL%NUNKH_LOCAL
+   ELSE
+      ZSL%LOWER_ROW = MAX(1,ZSL%UNKH_IND(NM_START))
+      ZSL%UPPER_ROW = MAX(1,ZSL%UNKH_IND(NM_START))
+   ENDIF
+
+   LIBRARY_SELECT: SELECT CASE(UGLMAT_SOLVER_LIBRARY)
+
+   CASE(MKL_CPARDISO_FLAG) LIBRARY_SELECT
+#ifdef WITH_MKL
    !--- This matrix definitoin used with MKL cluster solver -----
    ! Allocate A_H IA_H and JA_H matrices, considering all matrix coefficients:
    ALLOCATE ( ZSL%A_H(ZSL%TOT_NNZ_H) , ZSL%IA_H(MAX(ZSL%NUNKH_LOCAL,1)+1) , ZSL%JA_H(ZSL%TOT_NNZ_H) )
    ! Store upper triangular part of symmetric D_MAT_H in CSR format:
    IF(ZSL%NUNKH_LOCAL>0) THEN
-      ! Here each process defines de beginning and end rows in global numeration, for the equations
-      ! it has assembled:
-      ZSL%LOWER_ROW = ZSL%UNKH_IND(NM_START) + 1
-      ZSL%UPPER_ROW = ZSL%UNKH_IND(NM_START) + ZSL%NUNKH_LOCAL
-
       INNZ = 0
       DO IROW=1,ZSL%NUNKH_LOCAL
          ZSL%IA_H(IROW) = INNZ + 1
@@ -3913,8 +3976,6 @@ SUBROUTINE GET_H_MATRIX_LUDCMP
       ENDDO
       ZSL%IA_H(ZSL%NUNKH_LOCAL+1) = INNZ + 1
    ELSE
-      ZSL%LOWER_ROW = MAX(1,ZSL%UNKH_IND(NM_START))
-      ZSL%UPPER_ROW = MAX(1,ZSL%UNKH_IND(NM_START))
       ZSL%A_H       = 0._EB     ! Add a zero coefficient in A(ZSL%UNKH_IND(NM_START),ZSL%UNKH_IND(NM_START)).
       ZSL%IA_H(1:2) = (/1,2/)
       ZSL%JA_H(1)   = MAX(1,ZSL%UNKH_IND(NM_START))
@@ -3928,7 +3989,6 @@ SUBROUTINE GET_H_MATRIX_LUDCMP
    ALLOCATE( ZSL%A_H_FB(ZSL%TOT_NNZ_H) );   ZSL%A_H_FB(1:ZSL%TOT_NNZ_H)   = REAL(ZSL%A_H(1:ZSL%TOT_NNZ_H),FB)
 #endif
 
-#ifdef WITH_MKL
    ! Lower and uppper rows handled by this process:
    IPARM(41) = ZSL%LOWER_ROW
    IPARM(42) = ZSL%UPPER_ROW
@@ -4004,17 +4064,81 @@ SUBROUTINE GET_H_MATRIX_LUDCMP
       STOP_STATUS = SETUP_STOP
       RETURN
    ENDIF
+#endif
 
-#else
+   CASE(HYPRE_FLAG) LIBRARY_SELECT
+#ifdef WITH_HYPRE
 
-   IF (MY_RANK==0) WRITE(LU_ERR,'(A)') &
-   'Error: MKL Library compile flag was not defined for GLMAT as pressure solver.'
-   ! Some error - stop flag for CALL STOP_CHECK(1).
-   STOP_STATUS = SETUP_STOP
-   RETURN
+   IF (ALLOCATED(ZSL%HYPRE_ZSL%INDICES)) DEALLOCATE(ZSL%HYPRE_ZSL%INDICES)
+   ALLOCATE( ZSL%HYPRE_ZSL%INDICES(MAX(1,ZSL%NUNKH_LOCAL)) )
+
+   CALL HYPRE_IJMATRIXCREATE(MPI_COMM_WORLD,ZSL%LOWER_ROW-1,ZSL%UPPER_ROW-1,&
+                             ZSL%LOWER_ROW-1,ZSL%UPPER_ROW-1,ZSL%HYPRE_ZSL%A_H,HYPRE_IERR)
+   CALL HYPRE_IJMATRIXSETOBJECTTYPE(ZSL%HYPRE_ZSL%A_H,HYPRE_PARCSR,HYPRE_IERR)
+   CALL HYPRE_IJMATRIXINITIALIZE(ZSL%HYPRE_ZSL%A_H,HYPRE_IERR)
+   IF(ZSL%NUNKH_LOCAL > 0) THEN
+   ZSL%JD_MAT_H = ZSL%JD_MAT_H - 1
+   DO IROW=1,ZSL%NUNKH_LOCAL
+      ZSL%HYPRE_ZSL%INDICES(IROW)=ZSL%LOWER_ROW+IROW-2
+      CALL HYPRE_IJMATRIXSETVALUES(ZSL%HYPRE_ZSL%A_H, 1, ZSL%NNZ_D_MAT_H(IROW), ZSL%HYPRE_ZSL%INDICES(IROW), &
+            ZSL%JD_MAT_H(1:ZSL%NNZ_D_MAT_H(IROW),IROW), ZSL%D_MAT_H(1:ZSL%NNZ_D_MAT_H(IROW),IROW), HYPRE_IERR)
+   ENDDO
+   ELSE
+      ZSL%HYPRE_ZSL%INDICES(1)=ZSL%LOWER_ROW-1
+      ONEV(1) = 1; ZEROV(1) = 0._EB
+      ! Define a zero coefficient in A(ZSL%LOWER_ROW-1,ZSL%LOWER_ROW-1):
+      CALL HYPRE_IJMATRIXSETVALUES(ZSL%HYPRE_ZSL%A_H, 1, ONEV(1), ZSL%HYPRE_ZSL%INDICES(1), &
+      ZSL%HYPRE_ZSL%INDICES(1), ZEROV(1), HYPRE_IERR)
+   ENDIF
+   CALL HYPRE_IJMATRIXASSEMBLE(ZSL%HYPRE_ZSL%A_H, HYPRE_IERR)
+   CALL HYPRE_IJMATRIXGETOBJECT(ZSL%HYPRE_ZSL%A_H, ZSL%HYPRE_ZSL%PARCSR_A_H, HYPRE_IERR)
+   ! Create right hand side vector
+   CALL HYPRE_IJVECTORCREATE(MPI_COMM_WORLD, ZSL%LOWER_ROW-1, ZSL%UPPER_ROW-1, ZSL%HYPRE_ZSL%F_H, HYPRE_IERR)
+   CALL HYPRE_IJVECTORSETOBJECTTYPE(ZSL%HYPRE_ZSL%F_H, HYPRE_PARCSR, HYPRE_IERR)
+   CALL HYPRE_IJVECTORINITIALIZE(ZSL%HYPRE_ZSL%F_H, HYPRE_IERR)
+   ! Create solution vector
+   CALL HYPRE_IJVECTORCREATE(MPI_COMM_WORLD, ZSL%LOWER_ROW-1, ZSL%UPPER_ROW-1, ZSL%HYPRE_ZSL%X_H, HYPRE_IERR)
+   CALL HYPRE_IJVECTORSETOBJECTTYPE(ZSL%HYPRE_ZSL%X_H, HYPRE_PARCSR, HYPRE_IERR)
+   CALL HYPRE_IJVECTORINITIALIZE(ZSL%HYPRE_ZSL%X_H, HYPRE_IERR)
+   ! Set values
+   CALL HYPRE_IJVECTORSETVALUES(ZSL%HYPRE_ZSL%F_H, ZSL%NUNKH_LOCAL, ZSL%HYPRE_ZSL%INDICES, ZSL%F_H, HYPRE_IERR)
+   CALL HYPRE_IJVECTORSETVALUES(ZSL%HYPRE_ZSL%X_H, ZSL%NUNKH_LOCAL, ZSL%HYPRE_ZSL%INDICES, ZSL%X_H, HYPRE_IERR)
+   ! Assemble vectors
+   CALL HYPRE_IJVECTORASSEMBLE(ZSL%HYPRE_ZSL%F_H, HYPRE_IERR)
+   CALL HYPRE_IJVECTORASSEMBLE(ZSL%HYPRE_ZSL%X_H, HYPRE_IERR)
+   ! Get rhs and soln objects
+   CALL HYPRE_IJVECTORGETOBJECT(ZSL%HYPRE_ZSL%F_H, ZSL%HYPRE_ZSL%PAR_F_H, HYPRE_IERR)
+   CALL HYPRE_IJVECTORGETOBJECT(ZSL%HYPRE_ZSL%X_H, ZSL%HYPRE_ZSL%PAR_X_H, HYPRE_IERR)
+
+  ! Create solver (Parallel Compressed Sparse Row Preconditioned Conjugate Gradient)
+   CALL HYPRE_PARCSRPCGCREATE(MPI_COMM_WORLD, ZSL%HYPRE_ZSL%SOLVER, HYPRE_IERR)
+   CALL HYPRE_PARCSRPCGSETMAXITER(ZSL%HYPRE_ZSL%SOLVER, HYPRE_SOLVER_MAXIT, HYPRE_IERR)
+   CALL HYPRE_PARCSRPCGSETTOL(ZSL%HYPRE_ZSL%SOLVER, HYPRE_SOLVER_TOL, HYPRE_IERR)
+   CALL HYPRE_PARCSRPCGSETTWONORM(ZSL%HYPRE_ZSL%SOLVER, HYPRE_SOLVER_SETTWONORM, HYPRE_IERR)
+   CALL HYPRE_PARCSRPCGSETPRINTLEVEL(ZSL%HYPRE_ZSL%SOLVER, HYPRE_SOLVER_SETPRINTLEVEL, HYPRE_IERR)
+   CALL HYPRE_PARCSRPCGSETLOGGING(ZSL%HYPRE_ZSL%SOLVER, HYPRE_SOLVER_SETLOGGING, HYPRE_IERR)
 
+   ! Set up the Algebraic Multi-Grid (AMG) preconditioner and specify any parameters
+   CALL HYPRE_BOOMERAMGCREATE(ZSL%HYPRE_ZSL%PRECOND, HYPRE_IERR)
+   CALL HYPRE_BOOMERAMGSETPRINTLEVEL(ZSL%HYPRE_ZSL%PRECOND, HYPRE_PRECOND_SETPRINTLEVEL, HYPRE_IERR)
+   CALL HYPRE_BOOMERAMGSETCOARSENTYPE(ZSL%HYPRE_ZSL%PRECOND, HYPRE_PRECOND_COARSENINGTYPE, HYPRE_IERR)
+   CALL HYPRE_BOOMERAMGSETRELAXTYPE(ZSL%HYPRE_ZSL%PRECOND, HYPRE_PRECOND_SETRELAXTYPE, HYPRE_IERR)
+   CALL HYPRE_BOOMERAMGSETNUMSWEEPS(ZSL%HYPRE_ZSL%PRECOND, HYPRE_PRECOND_NUMSWEEPS, HYPRE_IERR)
+   CALL HYPRE_BOOMERAMGSETTOL(ZSL%HYPRE_ZSL%PRECOND, HYPRE_PRECOND_TOL, HYPRE_IERR)
+   CALL HYPRE_BOOMERAMGSETMAXITER(ZSL%HYPRE_ZSL%PRECOND, HYPRE_PRECOND_MAXITER, HYPRE_IERR)
+   CALL HYPRE_PARCSRPCGSETPRECOND(ZSL%HYPRE_ZSL%SOLVER, HYPRE_PRECOND_ID, ZSL%HYPRE_ZSL%PRECOND, HYPRE_IERR)
+   ! Solver setup
+   CALL HYPRE_PARCSRPCGSETUP(ZSL%HYPRE_ZSL%SOLVER, ZSL%HYPRE_ZSL%PARCSR_A_H,&
+                             ZSL%HYPRE_ZSL%PAR_F_H, ZSL%HYPRE_ZSL%PAR_X_H, HYPRE_IERR)
 #endif
 
+   END SELECT LIBRARY_SELECT
+
+   ! Deallocate MAT_H arrays:
+   IF (ALLOCATED(ZSL%NNZ_D_MAT_H)) DEALLOCATE(ZSL%NNZ_D_MAT_H)
+   IF (ALLOCATED(ZSL%D_MAT_H))     DEALLOCATE(ZSL%D_MAT_H)
+   IF (ALLOCATED(ZSL%JD_MAT_H))    DEALLOCATE(ZSL%JD_MAT_H)
+
 ENDDO IPZ_LOOP
 
 ! Set level MSG to 0 for solution:
@@ -5159,6 +5283,9 @@ END FUNCTION GRADIENT_WEIGHT
 SUBROUTINE FINISH_GLMAT_SOLVER
 
 USE MPI_F08
+#ifdef WITH_HYPRE
+USE HYPRE_INTERFACE
+#endif
 
 ! Local variables:
 INTEGER :: MAXFCT, MNUM, PHASE, NRHS, ERROR, MSGLVL
@@ -5180,13 +5307,21 @@ SUBROUTINE FINISH_GLMAT_SOLVER
 
 DO IPZ=0,N_ZONE
    ZSL => ZONE_SOLVE(IPZ); IF(ZSL%NUNKH_TOTAL==0) CYCLE
-
+   ! Finalize Cluster Sparse Solver:
+   IF (UGLMAT_SOLVER_LIBRARY==MKL_CPARDISO_FLAG) THEN
 #ifdef WITH_MKL
    CALL CLUSTER_SPARSE_SOLVER(ZSL%PT_H, MAXFCT, MNUM, ZSL%MTYPE, PHASE, ZSL%NUNKH_TOTAL, &
    ZSL%A_H, ZSL%IA_H, ZSL%JA_H, PERM, NRHS, IPARM, MSGLVL, ZSL%F_H, ZSL%X_H, MPI_COMM_WORLD, ERROR)
 #endif /* WITH_MKL */
+   ENDIF
 ENDDO
 
+IF (UGLMAT_SOLVER_LIBRARY==HYPRE_FLAG) THEN
+#ifdef WITH_HYPRE
+   CALL HYPRE_FINALIZE(HYPRE_IERR)
+#endif
+ENDIF
+
 DEALLOCATE(ZONE_SOLVE)
 
 RETURN