diff --git a/Source/cons.f90 b/Source/cons.f90 index 6dfa66c920e..1a929da5b5e 100644 --- a/Source/cons.f90 +++ b/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 diff --git a/Source/imkl.f90 b/Source/imkl.f90 index fc4f748f36e..edcd34fbff2 100644 --- a/Source/imkl.f90 +++ b/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 diff --git a/Source/pres.f90 b/Source/pres.f90 index b5b123bab62..f1e5b5be440 100644 --- a/Source/pres.f90 +++ b/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 diff --git a/Source/read.f90 b/Source/read.f90 index 507db9b0430..243656d1708 100644 --- a/Source/read.f90 +++ b/Source/read.f90 @@ -4431,7 +4431,7 @@ SUBROUTINE READ_COMB ENDDO COMB_LOOP 23 REWIND(LU_INPUT) ; INPUT_FILE_LINE_NUMBER = 0 -! Reread the COMB lines to check if user explicitly set DO_CHEM_LOAD_BALANCE to false. +! Reread the COMB lines to check if user explicitly set DO_CHEM_LOAD_BALANCE to false. ! For CVODE default DO_CHEM_LOAD_BALANCE is TRUE. For others, it is false. IF(TRIM(ODE_SOLVER)=='CVODE') THEN DO_CHEM_LOAD_BALANCE = .TRUE. @@ -5623,7 +5623,7 @@ SUBROUTINE PROC_REAC_2 RN%AIT_EXCLUSION_ZONE(IZ)%DEVC_INDEX,RN%AIT_EXCLUSION_ZONE(IZ)%CTRL_INDEX,IZ) ENDIF ENDDO - + IF (RN%SIMPLE_CHEMISTRY .AND. RN%HOC_COMPLETE > 2E8_EB) THEN WRITE(MESSAGE,'(A,I0,A)') 'WARNING: The heat of combustion for REACtion ',NR,' exceeds 200,000 kJ/kg.' IF (MY_RANK==0) WRITE(LU_ERR,'(A)') TRIM(MESSAGE) @@ -7083,7 +7083,7 @@ SUBROUTINE READ_MATL ML%RAMP_C_S = SPECIFIC_HEAT_RAMP ML%RAMP_K_S = CONDUCTIVITY_RAMP ML%RHO_S = DENSITY ! This is bulk density of pure material. - ML%REAC_RATE_DELTA = REAC_RATE_DELTA + ML%REAC_RATE_DELTA = REAC_RATE_DELTA ML%REFERENCE_ENTHALPY_TEMPERATURE = REFERENCE_ENTHALPY_TEMPERATURE + TMPM ML%REFERENCE_ENTHALPY = REFERENCE_ENTHALPY*1000._EB ML%RESIDUE_MATL_NAME = MATL_ID @@ -8550,7 +8550,7 @@ SUBROUTINE READ_SURF(QUICK_READ) ALLOCATE(SF%RAMP_IHS_INDEX(MAX_LAYERS)) ; SF%RAMP_IHS_INDEX = 0 ALLOCATE(SF%MINIMUM_LAYER_THICKNESS(MAX_LAYERS)) ; SF%MINIMUM_LAYER_THICKNESS = 0._EB ALLOCATE(SF%SWELL_RATIO(MAX_LAYERS)) ; SF%SWELL_RATIO = 1._EB - + COUNT_LAYERS: DO NL=1,MAX_LAYERS IF (THICKNESS(NL) < 0._EB) EXIT COUNT_LAYERS SF%N_LAYERS = SF%N_LAYERS + 1 @@ -8594,7 +8594,7 @@ SUBROUTINE READ_SURF(QUICK_READ) IF (SF%EMISSIVITY < 0._EB) SF%EMISSIVITY = EMISSIVITY_DEFAULT IF (SF%EMISSIVITY_BACK < 0._EB) SF%EMISSIVITY_BACK = EMISSIVITY_DEFAULT - + ! Define mass flux division point if the user does not specify. For all but ! surfaces with exposed backing, all pyrolyzed mass migrates to the front surface. @@ -8628,12 +8628,12 @@ SUBROUTINE READ_SURF(QUICK_READ) ENDIF ENDDO ENDDO - + WHERE (CHILD_SURF>1) CHILD_SURF = 1 SF%N_MATL = SUM(CHILD_SURF) ALLOCATE(SF%MATL_INDEX(SF%N_MATL)) ALLOCATE(SF%MATL_NAME(SF%N_MATL)) - + I=0 DO NN=1,N_MATL IF (CHILD_SURF(NN)==1) THEN @@ -8652,7 +8652,7 @@ SUBROUTINE READ_SURF(QUICK_READ) SF%TMP_GAS_BACK = TMP_GAS_BACK + TMPM IF (SF%TMP_GAS_BACK>0._EB) SF%BACKING = VOID - + ! Allocate parameters indexed by layer THERM_THICK = .FALSE. @@ -8698,7 +8698,7 @@ SUBROUTINE READ_SURF(QUICK_READ) ENDIF ! Thermal boundary conditions - + IF (SF%ADIABATIC .AND. (SF%NET_HEAT_FLUX < 1.E12_EB .OR. ABS(SF%CONVECTIVE_HEAT_FLUX)>TWO_EPSILON_EB)) THEN WRITE(MESSAGE,'(3A)') 'ERROR(336): SURF ',TRIM(SF%ID),' cannot use both ADIABATIC and NET or CONVECTIVE_HEAT_FLUX.' CALL SHUTDOWN(MESSAGE) ; RETURN @@ -9638,7 +9638,7 @@ SUBROUTINE PROC_WALL ! Compute node coordinates ! X_S_OLD provides the right size array into GET_WALL_NODE_COORDINATES. REMESH_LAYER defined as .TRUE. - + ALLOCATE(X_S_OLD(0:SF%N_CELLS_MAX)); X_S_OLD=0._EB CALL GET_WALL_NODE_COORDINATES(SF%N_CELLS_INI,SF%N_CELLS_MAX+1,SF%N_LAYERS,SF%N_LAYER_CELLS, N_LAYER_CELLS_OLD(1:SF%N_LAYERS), & SF%SMALLEST_CELL_SIZE(1:SF%N_LAYERS),SF%STRETCH_FACTOR(1:SF%N_LAYERS),REMESH_LAYER(1:SF%N_LAYERS),& @@ -9709,17 +9709,21 @@ SUBROUTINE READ_PRES SELECT CASE(TRIM(SOLVER)) - CASE('UGLMAT') + CASE('UGLMAT','UGLMAT PARDISO','UGLMAT HYPRE') PRES_METHOD = 'UGLMAT' PRES_FLAG = UGLMAT_FLAG PRES_ON_WHOLE_DOMAIN = .FALSE. IF (CHECK_POISSON) GLMAT_VERBOSE=.TRUE. + UGLMAT_SOLVER_LIBRARY = MKL_CPARDISO_FLAG + IF (TRIM(SOLVER)=='UGLMAT HYPRE') UGLMAT_SOLVER_LIBRARY = HYPRE_FLAG - CASE('GLMAT') + CASE('GLMAT','GLMAT PARDISO','GLMAT HYPRE') PRES_METHOD = 'GLMAT' PRES_FLAG = GLMAT_FLAG PRES_ON_WHOLE_DOMAIN = .TRUE. IF (CHECK_POISSON) GLMAT_VERBOSE=.TRUE. + UGLMAT_SOLVER_LIBRARY = MKL_CPARDISO_FLAG + IF (TRIM(SOLVER)=='GLMAT HYPRE') UGLMAT_SOLVER_LIBRARY = HYPRE_FLAG CASE('ULMAT','ULMAT PARDISO','ULMAT HYPRE') PRES_METHOD = 'ULMAT' @@ -9736,36 +9740,40 @@ SUBROUTINE READ_PRES CALL SHUTDOWN('ERROR(371): Pressure solver '//TRIM(SOLVER)//' not known.') ; RETURN END SELECT -ULMAT_IF: IF (PRES_FLAG==ULMAT_FLAG) THEN - ! If library is not specified and only one library is linked, choose the available library #ifndef WITH_MKL #ifdef WITH_HYPRE -ULMAT_SOLVER_LIBRARY=HYPRE_FLAG +IF (PRES_FLAG==ULMAT_FLAG ) ULMAT_SOLVER_LIBRARY =HYPRE_FLAG +IF (PRES_FLAG==UGLMAT_FLAG) UGLMAT_SOLVER_LIBRARY=HYPRE_FLAG #endif #endif #ifndef WITH_HYPRE #ifdef WITH_MKL -ULMAT_SOLVER_LIBRARY=MKL_PARDISO_FLAG +IF (PRES_FLAG==ULMAT_FLAG ) ULMAT_SOLVER_LIBRARY =MKL_PARDISO_FLAG +IF (PRES_FLAG==UGLMAT_FLAG) UGLMAT_SOLVER_LIBRARY=MKL_CPARDISO_FLAG #endif #endif ! If neither library is specified, throw an error #ifndef WITH_MKL -IF (ULMAT_SOLVER_LIBRARY==MKL_PARDISO_FLAG) THEN +IF (PRES_FLAG==ULMAT_FLAG .AND. ULMAT_SOLVER_LIBRARY==MKL_PARDISO_FLAG) THEN CALL SHUTDOWN('ERROR: MKL PARDISO selected for ULMAT solver without compiling and linking MKL library.') RETURN +ELSEIF (PRES_FLAG==UGLMAT_FLAG .AND. UGLMAT_SOLVER_LIBRARY==MKL_CPARDISO_FLAG) THEN + CALL SHUTDOWN('ERROR: MKL CLUSTER PARDISO selected for UGLMAT solver without compiling and linking MKL library.') + RETURN ENDIF #endif #ifndef WITH_HYPRE -IF (ULMAT_SOLVER_LIBRARY==HYPRE_FLAG) THEN +IF (PRES_FLAG==ULMAT_FLAG .AND. ULMAT_SOLVER_LIBRARY==HYPRE_FLAG) THEN CALL SHUTDOWN('ERROR: HYPRE selected for ULMAT solver without compiling and linking HYPRE library.') RETURN +ELSEIF (PRES_FLAG==UGLMAT_FLAG .AND. UGLMAT_SOLVER_LIBRARY==HYPRE_FLAG) THEN + CALL SHUTDOWN('ERROR: HYPRE selected for UGLMAT solver without compiling and linking HYPRE library.') + RETURN ENDIF #endif -ENDIF ULMAT_IF - ! Determine how many pressure iterations to perform per half time step. IF (VELOCITY_TOLERANCE>100._EB) THEN @@ -12922,19 +12930,19 @@ SUBROUTINE READ_INIT IN%U0 = UVW(1) IN%V0 = UVW(2) IN%W0 = UVW(3) - - + + ! Check for possible issues when N_PARTICLES given N_PARTICLE_IF: IF (N_PARTICLES > 0) THEN - + MESH_COUNT = 0 DO NM=1,NMESHES IF (IN%X1>MESHES(NM)%XF .OR. IN%X2MESHES(NM)%YF .OR. IN%Y2MESHES(NM)%ZF .OR. IN%Z2 N_PARTICLES) THEN WRITE(MESSAGE,'(3A)') 'WARNING: INIT ',TRIM(ID),' XB contains more MESHes than N_PARTICLES.' - IF (MY_RANK==0) WRITE(LU_ERR,'(A)') TRIM(MESSAGE) - ENDIF + IF (MY_RANK==0) WRITE(LU_ERR,'(A)') TRIM(MESSAGE) + ENDIF ENDIF - + ENDIF N_PARTICLE_IF ENDDO I_MULT_LOOP diff --git a/Source/type.f90 b/Source/type.f90 index 7c8fff082ad..ecbda581ed8 100644 --- a/Source/type.f90 +++ b/Source/type.f90 @@ -532,8 +532,8 @@ MODULE TYPES REAL(EB) :: ODE_REL_ERROR !< Relative error for finite rate chemistry REAL(EB) :: POLYNOMIAL_TEMP(4) !< Temperature bands for user polynomial REAL(EB) :: POLYNOMIAL_COEFF(9,3) !< Coefficients for user polynomial - REAL(EB) :: REAL_REFRACTIVE_INDEX - REAL(EB) :: COMPLEX_REFRACTIVE_INDEX + REAL(EB) :: REAL_REFRACTIVE_INDEX + REAL(EB) :: COMPLEX_REFRACTIVE_INDEX LOGICAL :: ISFUEL=.FALSE. !< Fuel species LOGICAL :: LISTED=.FALSE. !< Properties are known to FDS @@ -1732,6 +1732,11 @@ MODULE TYPES #else INTEGER, ALLOCATABLE :: PT_H(:) #endif /* WITH_MKL */ +#ifdef WITH_HYPRE + TYPE(HYPRE_ZM_TYPE) HYPRE_ZSL +#else + INTEGER :: HYPRE_ZSL +#endif /* WITH_HYPRE */ INTEGER :: NUNKH_LOCAL=0 !< SUM(NUNKH_LOC(LOWER_MESH_INDEX:UPPER_MESH_INDEX)). INTEGER :: NUNKH_TOTAL=0 !< SUM(NUNKH_TOT(LOWER_MESH_INDEX:UPPER_MESH_INDEX)). INTEGER :: TOT_NNZ_H=0 !< Total number of non-zeros owned by this process for a pres zone. diff --git a/Verification/FDS_Cases.sh b/Verification/FDS_Cases.sh index 23178e5b313..a2badee239c 100755 --- a/Verification/FDS_Cases.sh +++ b/Verification/FDS_Cases.sh @@ -394,6 +394,7 @@ $QFDS -p 5 -d Pressure_Solver dancing_eddies_embed.fds $QFDS -p 4 -d Pressure_Solver dancing_eddies_ulmat.fds $QFDS -p 4 -d Pressure_Solver dancing_eddies_ulmat_hypre.fds $QFDS -p 4 -d Pressure_Solver dancing_eddies_uglmat.fds +$QFDS -p 4 -d Pressure_Solver dancing_eddies_uglmat_hypre.fds $QFDS -p 4 -d Pressure_Solver dancing_eddies_tight.fds $QFDS -p 4 -d Pressure_Solver dancing_eddies_tight_no_precon.fds $QFDS -p 4 -d Pressure_Solver dancing_eddies_default.fds diff --git a/Verification/Pressure_Solver/dancing_eddies_uglmat_hypre.fds b/Verification/Pressure_Solver/dancing_eddies_uglmat_hypre.fds new file mode 100644 index 00000000000..6f8c26979b8 --- /dev/null +++ b/Verification/Pressure_Solver/dancing_eddies_uglmat_hypre.fds @@ -0,0 +1,36 @@ +&HEAD CHID='dancing_eddies_uglmat_hypre', TITLE='Simple 2D Tunnel test of UGLMAT HYPRE solver'/ + +&MULT ID='mesh',DX=0.075,I_UPPER=3 / +&MESH IJK=75,1,40, XB=0.0,0.075,-0.001,0.001,-0.02,0.020, MULT_ID='mesh' / + +&TIME T_END=2.0 / + +&MISC SIMULATION_MODE='DNS', STRATIFICATION=.FALSE., NOISE=.FALSE. / + +&SPEC ID='LJ AIR', VISCOSITY=1.2E-5, BACKGROUND=.TRUE. / + +&PRES CHECK_POISSON=T, SOLVER='UGLMAT HYPRE', MAX_PRESSURE_ITERATIONS=1 / + +&SURF ID='WALL', FREE_SLIP=.TRUE., DEFAULT=.TRUE. / +&SURF ID='PLATE', NO_SLIP=.TRUE. / + +&OBST XB=0.020,0.021,-0.001,0.001, 0.004,0.005, SURF_ID='PLATE' / +&OBST XB=0.100,0.102,-0.001,0.001,-0.005,0.005, SURF_ID='PLATE' / + +&SURF ID='BLOW', VEL=-0.5 / + +&VENT PBX=0.0, SURF_ID='BLOW' / +&VENT PBX=0.3, SURF_ID='OPEN' / + +&SLCF PBY=0.,QUANTITY='VELOCITY',VECTOR=.TRUE. / +&SLCF PBY=0.,QUANTITY='H' / +&SLCF PBY=0.,QUANTITY='PRESSURE' / + +&DEVC XYZ=0.2,0.,0.0, QUANTITY='PRESSURE', ID='pres' / +&DEVC XYZ=0.,0.,0., QUANTITY='MAXIMUM VELOCITY ERROR', ID='error' / +&DEVC XYZ=0.,0.,0., QUANTITY='PRESSURE ITERATIONS', ID='iter' / +&DEVC XYZ=0.,0.,0., QUANTITY='CPU TIME', ID='cpu' / + +&TAIL / + +