Merge pull request #14724 from marcosvanella/master

FDS Source : Address semidefinite matrices in global pressure solver.

Author: marcosvanella

Source/pres.f90

@@ -1518,7 +1518,7 @@ SUBROUTINE ULMAT_SOLVE_ZONE(NM,IPZ)
    WALL_CELL_LOOP_0 : DO IW=N_EXTERNAL_WALL_CELLS+1,N_EXTERNAL_WALL_CELLS+N_INTERNAL_WALL_CELLS
       WC => WALL(IW)
       IF (WC%BOUNDARY_TYPE/=SOLID_BOUNDARY .OR. WC%CUT_FACE_INDEX>0) CYCLE WALL_CELL_LOOP_0
-      BC  => BOUNDARY_COORD(WC%BC_INDEX); 
+      BC  => BOUNDARY_COORD(WC%BC_INDEX);
       IOR = BC%IOR; IF(ABS(IOR)/=1) CYCLE WALL_CELL_LOOP_0 ! Only in x.
       ! Gasphase cell indexes:
       IIG = BC%IIG; JJG = BC%JJG; KKG = BC%KKG
@@ -1540,7 +1540,7 @@ SUBROUTINE ULMAT_SOLVE_ZONE(NM,IPZ)
          DHDN = (H_BAR(I_OFFSET(NM)+IIG)-H_BAR(I_OFFSET(NM)+IIG-1))*TP_RDXN(I_OFFSET(NM)+IIG-1)
          AF  =  ((1._EB-CYL_FCT)*DY(JJG) + CYL_FCT*R(IIG-1)) * DZ(KKG)
          VAL =  DHDN*AF
-      END SELECT   
+      END SELECT
       ! Add to F_H:
       ZM%F_H(IROW) = ZM%F_H(IROW) - VAL
    ENDDO WALL_CELL_LOOP_0
@@ -3054,11 +3054,15 @@ MODULE GLOBMAT_SOLVER
 REAL(EB), SAVE :: CYL_FCT
 
 ! Pressure zone loops index:
-INTEGER :: IPZ
+INTEGER :: IPZ, N_ZONE_GLOBMAT
 
 ! Handle for ZONE_SOLVER array entries:
 TYPE(ZONE_SOLVE_TYPE), POINTER :: ZSL
 
+! Matrix types:
+INTEGER, PARAMETER :: SYMM_INDEFINITE       =-2
+INTEGER, PARAMETER :: SYMM_POSITIVE_DEFINITE= 2
+
 !#define SINGLE_PRECISION_PSN_SOLVE
 
 ! Timing variable:
@@ -3120,7 +3124,7 @@ SUBROUTINE GLMAT_SOLVER(T,DT)
    CALL GET_PRES_CFACE_BCS(NM,T,DT)
 ENDDO
 
-IPZ_LOOP : DO IPZ=0,N_ZONE
+IPZ_LOOP : DO IPZ=0,N_ZONE_GLOBMAT
 
    ZSL => ZONE_SOLVE(IPZ)
 
@@ -3138,7 +3142,7 @@ SUBROUTINE GLMAT_SOLVER(T,DT)
       CALL GET_FH_FROM_PRHS_AND_BCS(NM,DT,CYL_FCT,UNKH,ZSL%NUNKH_LOCAL,IPZ,ZSL%F_H)
    ENDDO MESH_LOOP_1
 
-   IF (ZSL%MTYPE==-2) THEN
+   IF (ZSL%MTYPE==SYMM_INDEFINITE) THEN
       SUM_FH = 0._EB; MEAN_FH = 0._EB
       WHOLE_DOM_IF1 : IF(.NOT.PRES_ON_WHOLE_DOMAIN) THEN
          ! Sum source F_H by Pressure Zone:
@@ -3203,6 +3207,7 @@ SUBROUTINE GLMAT_SOLVER(T,DT)
 
    CASE(HYPRE_FLAG) LIBRARY_SELECT
 #ifdef WITH_HYPRE
+   IF (ZSL%MTYPE==SYMM_INDEFINITE .AND. ZSL%UPPER_ROW==ZSL%NUNKH_TOTAL) ZSL%F_H(MAX(ZSL%NUNKH_LOCAL,1)) = 0._EB
    ! 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)
@@ -3217,7 +3222,7 @@ SUBROUTINE GLMAT_SOLVER(T,DT)
 
    END SELECT LIBRARY_SELECT
 
-   IF (ZSL%MTYPE==-2) THEN
+   IF (ZSL%MTYPE==SYMM_INDEFINITE) THEN
       SUM_XH = 0._EB; MEAN_XH = 0._EB
       WHOLE_DOM_IF2 : IF(.NOT.PRES_ON_WHOLE_DOMAIN) THEN
          ! Sum H by Pressure Zone:
@@ -3381,6 +3386,9 @@ SUBROUTINE GLMAT_SOLVER_SETUP(STAGE_FLAG)
 #endif
 #endif
 
+N_ZONE_GLOBMAT = N_ZONE
+IF(PRES_ON_WHOLE_DOMAIN) N_ZONE_GLOBMAT = 0
+
 SELECT CASE(STAGE_FLAG)
 CASE(1)
 
@@ -3923,7 +3931,7 @@ SUBROUTINE GET_H_MATRIX_LUDCMP
 INTEGER, ALLOCATABLE, DIMENSION(:,:) :: MB_FACTOR
 #endif
 #ifdef WITH_HYPRE
-INTEGER :: ONEV(1)
+INTEGER :: ONEV(1), END_ROW
 REAL(EB) :: ZEROV(1)
 #endif
 !.. All other variables
@@ -3957,7 +3965,7 @@ SUBROUTINE GET_H_MATRIX_LUDCMP
 #endif
 END SELECT
 
-IPZ_LOOP : DO IPZ=0,N_ZONE
+IPZ_LOOP : DO IPZ=0,N_ZONE_GLOBMAT
 
    ZSL => ZONE_SOLVE(IPZ); IF(ZSL%NUNKH_TOTAL==0) CYCLE
 
@@ -4092,21 +4100,49 @@ SUBROUTINE GET_H_MATRIX_LUDCMP
 
    CASE(HYPRE_FLAG) LIBRARY_SELECT
 #ifdef WITH_HYPRE
-
    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%MTYPE==SYMM_INDEFINITE) THEN
+      IF(ZSL%NUNKH_TOTAL==1) THEN ! Single unknown, zero coefficient matrix (1,1). Set coefficient to 1.
+         IF(ZSL%NUNKH_LOCAL==1) THEN
+            ZSL%NNZ_D_MAT_H(1) = 1
+            DEALLOCATE(ZSL%JD_MAT_H); ALLOCATE(ZSL%JD_MAT_H(1,1)); ZSL%JD_MAT_H(1,1) = 1
+            DEALLOCATE(ZSL%D_MAT_H ); ALLOCATE(ZSL%D_MAT_H(1,1) ); ZSL%D_MAT_H(1,1)  = 1._EB
+         ENDIF
+      ELSE ! More than one unknown
+         IF(ZSL%UPPER_ROW==ZSL%NUNKH_TOTAL) THEN
+            END_ROW = ZSL%NUNKH_LOCAL-1
+         ELSE
+            END_ROW = ZSL%NUNKH_LOCAL
+         ENDIF
+         ! Rows 1 to ZM%NUNKH-1, last column, set all to zero:
+         DO IROW=1,END_ROW
+            DO JCOL=1,ZSL%NNZ_D_MAT_H(IROW)
+               IF ( ZSL%JD_MAT_H(JCOL,IROW) /= ZSL%NUNKH_TOTAL ) CYCLE ! Make zero matrix entries in last column.
+               ZSL%D_MAT_H(JCOL,IROW) = 0._EB
+            ENDDO
+         ENDDO
+         ! Last row, all zeros except the diagonal that keeps diagonal number:
+         IF(ZSL%UPPER_ROW==ZSL%NUNKH_TOTAL) THEN
+            IROW = ZSL%NUNKH_LOCAL
+            DO JCOL=1,ZSL%NNZ_D_MAT_H(IROW)
+               IF ( ZSL%JD_MAT_H(JCOL,IROW) /= ZSL%NUNKH_TOTAL ) ZSL%D_MAT_H(JCOL,IROW) = 0._EB
+            ENDDO
+         ENDIF
+      ENDIF
+   ENDIF
    IF(ZSL%NUNKH_LOCAL > 0) THEN
-   ZSL%JD_MAT_H = ZSL%JD_MAT_H - 1
-   DO IROW=1,ZSL%NUNKH_LOCAL
+      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
+      ENDDO
    ELSE
       ZSL%HYPRE_ZSL%INDICES(1)=ZSL%LOWER_ROW-1
       ONEV(1) = 1; ZEROV(1) = 0._EB
@@ -4226,7 +4262,7 @@ SUBROUTINE GET_BCS_H_MATRIX
 REAL(EB), POINTER, DIMENSION(:,:) :: D_MAT_HP
 INTEGER :: H_MAT_IVEC
 
-IPZ_LOOP : DO IPZ=0,N_ZONE
+IPZ_LOOP : DO IPZ=0,N_ZONE_GLOBMAT
    ZSL => ZONE_SOLVE(IPZ)
 
    ! Allocate D_MAT_H:
@@ -4300,9 +4336,9 @@ SUBROUTINE GET_BCS_H_MATRIX
    ! Here all reduce with sum among MPI processes:
    IF (N_MPI_PROCESSES > 1) CALL MPI_ALLREDUCE(MPI_IN_PLACE,H_MAT_IVEC,1,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,IERR)
    IF ( H_MAT_IVEC > 0 ) THEN
-      ZSL%MTYPE = 2 ! At least one Dirichlet BC, matrix positive definite.
+      ZSL%MTYPE = SYMM_POSITIVE_DEFINITE ! At least one Dirichlet BC, matrix positive definite.
    ELSE
-      ZSL%MTYPE =-2 ! Indefinite.
+      ZSL%MTYPE = SYMM_INDEFINITE ! Indefinite.
    ENDIF
 ENDDO IPZ_LOOP
 
@@ -4332,7 +4368,7 @@ SUBROUTINE GET_H_MATRIX
 INTEGER :: WC_JD(1:2,1:2)
 LOGICAL :: FLG
 
-IPZ_LOOP : DO IPZ=0,N_ZONE
+IPZ_LOOP : DO IPZ=0,N_ZONE_GLOBMAT
 
    ZSL => ZONE_SOLVE(IPZ)
 
@@ -4553,7 +4589,7 @@ SUBROUTINE GET_MATRIXGRAPH_H_WHLDOM
 ENDDO
 
 ! Now Loop by Pressure Zone:
-IPZ_LOOP : DO IPZ=0,N_ZONE
+IPZ_LOOP : DO IPZ=0,N_ZONE_GLOBMAT
 
    ZSL => ZONE_SOLVE(IPZ)
    ZSL%NUNKH_LOCAL = SUM(ZSL%NUNKH_LOC(LOWER_MESH_INDEX:UPPER_MESH_INDEX))
@@ -4959,14 +4995,14 @@ SUBROUTINE GET_MATRIX_INDEXES_H
    ! Initialize unknown numbers for H:
    MESHES(NM)%CCVAR(:,:,:,UNKH) = IS_UNDEFINED
    ! Select the parent zone as the first in the row
-   DO IPZ=0,N_ZONE
+   DO IPZ=0,N_ZONE_GLOBMAT
       ZSL => ZONE_SOLVE(IPZ)
       IF(.NOT.PRES_ON_WHOLE_DOMAIN) &
       ZSL%CONNECTED_ZONE_PARENT = MINLOC(CONNECTED_ZONES(IPZ,:), DIM=1, MASK=CONNECTED_ZONES(IPZ,:)/=0) - 1
    ENDDO
 ENDDO
 
-PRES_ZONE_LOOP : DO IPZ=0,N_ZONE
+PRES_ZONE_LOOP : DO IPZ=0,N_ZONE_GLOBMAT
 
    ZSL => ZONE_SOLVE(IPZ)
 
@@ -5322,7 +5358,7 @@ SUBROUTINE FINISH_GLMAT_SOLVER
 ! Finalize Pardiso:
 PHASE = -1
 
-DO IPZ=0,N_ZONE
+DO IPZ=0,N_ZONE_GLOBMAT
    ZSL => ZONE_SOLVE(IPZ); IF(ZSL%NUNKH_TOTAL==0) CYCLE
    ! Finalize Cluster Sparse Solver:
    IF (UGLMAT_SOLVER_LIBRARY==MKL_CPARDISO_FLAG) THEN

Source/read.f90

@@ -9880,7 +9880,7 @@ SUBROUTINE READ_PRES
 #ifndef WITH_MKL
 #ifdef WITH_HYPRE
 IF (PRES_FLAG==ULMAT_FLAG ) ULMAT_SOLVER_LIBRARY =HYPRE_FLAG
-IF (PRES_FLAG==UGLMAT_FLAG) UGLMAT_SOLVER_LIBRARY=HYPRE_FLAG
+IF (PRES_FLAG==UGLMAT_FLAG .OR. PRES_FLAG==GLMAT_FLAG) UGLMAT_SOLVER_LIBRARY=HYPRE_FLAG
 #endif
 #endif
 #ifndef WITH_HYPRE
@@ -11065,7 +11065,7 @@ SUBROUTINE READ_OBST(QUICK_READ)
                   IF (SURF_ID6(5)      ==SURFACE(NNN)%ID) OB%SURF_INDEX(-3)   = NNN
                   IF (SURF_ID6(6)      ==SURFACE(NNN)%ID) OB%SURF_INDEX( 3)   = NNN
                ENDDO
-               
+
                DO NNN=-3,3
                   IF (SURFACE(OB%SURF_INDEX(NNN))%NODE_ID/='null') THEN
                      WRITE(MESSAGE,'(A,A,A)')  'ERROR(xxx): OBST_ID ',TRIM(ID),' cannot have a SURF with NODE_ID.'
@@ -12219,7 +12219,7 @@ SUBROUTINE READ_VENT
                ENDDO
                IF (SURFACE(VT%SURF_INDEX)%NODE_ID /='null') THEN
                   WRITE(MESSAGE,'(3A)') 'ERROR(812): VENT ',TRIM(ID),' Use SURF_ID HVAC instead of NODE_ID.'
-                  CALL SHUTDOWN(MESSAGE,PROCESS_0_ONLY=.FALSE.) ; RETURN               
+                  CALL SHUTDOWN(MESSAGE,PROCESS_0_ONLY=.FALSE.) ; RETURN
                ENDIF
 
                IF (SURF_ID=='OPEN')                            VT%TYPE_INDICATOR =  2