Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

FDS Source : Add HYPRE PG+AMG to UGLMAT. #13826

Merged
merged 2 commits into from
Dec 3, 2024
+235 −48
Merged

FDS Source : Add HYPRE PG+AMG to UGLMAT. #13826

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
25025bd
FDS Source : Add HYPRE PG+AMG to UGLMAT.
marcosvanella Dec 3, 2024
3991b5e
FDS Source : Initialize variables.
marcosvanella Dec 3, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions Source/cons.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
2 changes: 1 addition & 1 deletion Source/imkl.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
169 changes: 152 additions & 17 deletions Source/pres.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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

Expand All @@ -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
Expand All @@ -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))
Expand All @@ -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
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down
Loading
Loading