From 98032baf65203c6b3ff5104b9b32cd98ddf4cd0b Mon Sep 17 00:00:00 2001 From: marcosvanella Date: Fri, 26 Sep 2025 17:05:06 -0400 Subject: [PATCH] FDS Source: GLOBMAT_SOLVER change matrix building arrays to handle variable size stencils. --- Source/ccib.f90 | 119 ++++++++++++++++++++--------- Source/pres.f90 | 193 ++++++++++++++++++++++++------------------------ Source/type.f90 | 12 ++- 3 files changed, 191 insertions(+), 133 deletions(-) diff --git a/Source/ccib.f90 b/Source/ccib.f90 index bd2cc423810..7f8f651fbbc 100644 --- a/Source/ccib.f90 +++ b/Source/ccib.f90 @@ -19680,22 +19680,17 @@ END SUBROUTINE GET_CUTCELL_FH ! ---------------------------- GET_H_MATRIX_CC ------------------------------------ -SUBROUTINE GET_H_MATRIX_CC(NM,NM1,IPZ,D_MAT_HP) +SUBROUTINE GET_H_MATRIX_CC(NM,NM1,IPZ) ! This routine assumes the calling subroutine has called POINT_TO_MESH for NM. INTEGER, INTENT(IN) :: NM,NM1,IPZ -REAL(EB), POINTER, DIMENSION(:,:) :: D_MAT_HP ! Local Variables: INTEGER :: X1AXIS,IFACE,ICF,I,J,K,IND(LOW_IND:HIGH_IND),IND_LOC(LOW_IND:HIGH_IND) INTEGER :: LOCROW_1,LOCROW_2,ILOC,JLOC,JCOL,IROW REAL(EB) :: AF,IDX,BIJ,KFACE(2,2) -IF (.NOT. ASSOCIATED(D_MAT_HP)) THEN - WRITE(LU_ERR,*) 'GET_H_MATRIX_CC in geom.f90: Pointer D_MAT_HP not associated.' - RETURN -ENDIF ! X direction bounds: ILO_FACE = 0 ! Low mesh boundary face index. @@ -19752,8 +19747,7 @@ SUBROUTINE GET_H_MATRIX_CC(NM,NM1,IPZ,D_MAT_HP) DO JLOC=LOW_IND,HIGH_IND ! Local col number in Kface, JD IROW=IND_LOC(ILOC) ! Process Local Unknown number. JCOL=RCF%JDH(ILOC,JLOC) ! Local position of coef in D_MAT_H - ! Add coefficient: - D_MAT_HP(JCOL,IROW) = D_MAT_HP(JCOL,IROW) + KFACE(ILOC,JLOC) + ZONE_SOLVE(IPZ)%ROW_H(IROW)%D(JCOL) = ZONE_SOLVE(IPZ)%ROW_H(IROW)%D(JCOL) + KFACE(ILOC,JLOC) ENDDO ENDDO ENDDO @@ -19801,8 +19795,7 @@ SUBROUTINE GET_H_MATRIX_CC(NM,NM1,IPZ,D_MAT_HP) DO JLOC=LOW_IND,HIGH_IND ! Local col number in Kface, JD IROW=IND_LOC(ILOC) JCOL=CF%JDH(ILOC,JLOC,IFACE) - ! Add coefficient: - D_MAT_HP(JCOL,IROW) = D_MAT_HP(JCOL,IROW) + KFACE(ILOC,JLOC) + ZONE_SOLVE(IPZ)%ROW_H(IROW)%D(JCOL) = ZONE_SOLVE(IPZ)%ROW_H(IROW)%D(JCOL) + KFACE(ILOC,JLOC) ENDDO ENDDO ENDDO @@ -19932,9 +19925,9 @@ SUBROUTINE GET_CC_MATRIXGRAPH_H(NM,NM1,IPZ,LOOP_FLAG) ! Add to global matrix arrays: DO LOCROW=LOCROW_1,LOCROW_2 DO IIND=LOW_IND,HIGH_IND - NII = ZONE_SOLVE(IPZ)%NNZ_D_MAT_H(IND_LOC(LOCROW)) + NII = ZONE_SOLVE(IPZ)%ROW_H(IND_LOC(LOCROW))%NNZ DO ILOC=1,NII - IF ( IND(IIND) == ZONE_SOLVE(IPZ)%JD_MAT_H(ILOC,IND_LOC(LOCROW)) ) THEN + IF ( IND(IIND) == ZONE_SOLVE(IPZ)%ROW_H(IND_LOC(LOCROW))%JD(ILOC) ) THEN RCF%JDH(LOCROW,IIND) = ILOC EXIT ENDIF @@ -19976,9 +19969,9 @@ SUBROUTINE GET_CC_MATRIXGRAPH_H(NM,NM1,IPZ,LOOP_FLAG) ! Add to global matrix arrays: DO LOCROW=LOCROW_1,LOCROW_2 DO IIND=LOW_IND,HIGH_IND - NII = ZONE_SOLVE(IPZ)%NNZ_D_MAT_H(IND_LOC(LOCROW)) + NII = ZONE_SOLVE(IPZ)%ROW_H(IND_LOC(LOCROW))%NNZ DO ILOC=1,NII - IF ( IND(IIND) == ZONE_SOLVE(IPZ)%JD_MAT_H(ILOC,IND_LOC(LOCROW)) ) THEN + IF ( IND(IIND) == ZONE_SOLVE(IPZ)%ROW_H(IND_LOC(LOCROW))%JD(ILOC) ) THEN CF%JDH(LOCROW,IIND,IFACE) = ILOC EXIT ENDIF @@ -19999,38 +19992,94 @@ SUBROUTINE ADD_INPLACE_NNZ_H_WHLDOM(LOCROW_1,LOCROW_2,IND,IND_LOC,IPZ) INTEGER, INTENT(IN) :: LOCROW_1,LOCROW_2,IND(LOW_IND:HIGH_IND),IND_LOC(LOW_IND:HIGH_IND),IPZ ! Local Variables: -INTEGER LOCROW, IIND, NII, ILOC, JLOC +INTEGER LOCROW, IIND, ILOC LOGICAL INLIST LOCROW_LOOP : DO LOCROW=LOCROW_1,LOCROW_2 DO IIND=LOW_IND,HIGH_IND - NII = ZONE_SOLVE(IPZ)%NNZ_D_MAT_H(IND_LOC(LOCROW)) - ! Check that column index hasn't been already counted: - INLIST = .FALSE. - DO ILOC=1,NII - IF ( IND(IIND) == ZONE_SOLVE(IPZ)%JD_MAT_H(ILOC,IND_LOC(LOCROW)) ) THEN - INLIST = .TRUE. - EXIT + ! Populate variable per-row storage ROW_H only: + ILOC = IND_LOC(LOCROW) + IF (ILOC>=1 .AND. ILOC<=SIZE(ZONE_SOLVE(IPZ)%ROW_H)) THEN + IF (.NOT. ALLOCATED(ZONE_SOLVE(IPZ)%ROW_H(ILOC)%JD)) THEN + ALLOCATE(ZONE_SOLVE(IPZ)%ROW_H(ILOC)%JD(7)) + ALLOCATE(ZONE_SOLVE(IPZ)%ROW_H(ILOC)%D(7)) + ZONE_SOLVE(IPZ)%ROW_H(ILOC)%JD(1) = IND(IIND) + ZONE_SOLVE(IPZ)%ROW_H(ILOC)%D(1) = 0._EB + ZONE_SOLVE(IPZ)%ROW_H(ILOC)%NNZ = 1 + ELSE + INLIST = ANY(ZONE_SOLVE(IPZ)%ROW_H(ILOC)%JD(1:ZONE_SOLVE(IPZ)%ROW_H(ILOC)%NNZ) == IND(IIND)) + IF (.NOT. INLIST) THEN + CALL INSERT_SORTED_INT_REAL( & + ZONE_SOLVE(IPZ)%ROW_H(ILOC)%JD, & + ZONE_SOLVE(IPZ)%ROW_H(ILOC)%D, & + ZONE_SOLVE(IPZ)%ROW_H(ILOC)%NNZ, & + IND(IIND)) + ENDIF ENDIF - ENDDO - IF ( INLIST ) CYCLE - - ! Now add in place: - NII = NII + 1 - DO ILOC=1,NII - IF ( ZONE_SOLVE(IPZ)%JD_MAT_H(ILOC,IND_LOC(LOCROW)) > IND(IIND) ) EXIT - ENDDO - DO JLOC=NII,ILOC+1,-1 - ZONE_SOLVE(IPZ)%JD_MAT_H(JLOC,IND_LOC(LOCROW)) = ZONE_SOLVE(IPZ)%JD_MAT_H(JLOC-1,IND_LOC(LOCROW)) - ENDDO - ZONE_SOLVE(IPZ)%NNZ_D_MAT_H(IND_LOC(LOCROW)) = NII - ZONE_SOLVE(IPZ)%JD_MAT_H(ILOC,IND_LOC(LOCROW)) = IND(IIND) + ENDIF ENDDO ENDDO LOCROW_LOOP RETURN END SUBROUTINE ADD_INPLACE_NNZ_H_WHLDOM +! Helper: insert NEWCOL into sorted integer array JD and mirror REAL array D, grow allocs +SUBROUTINE INSERT_SORTED_INT_REAL(JD, D, NNZ, NEWCOL) +USE PRECISION_PARAMETERS +IMPLICIT NONE (TYPE,EXTERNAL) +INTEGER, ALLOCATABLE, INTENT(INOUT) :: JD(:) +REAL(EB), ALLOCATABLE, INTENT(INOUT) :: D(:) +INTEGER, INTENT(INOUT) :: NNZ +INTEGER, INTENT(IN) :: NEWCOL + +INTEGER :: POS, CAP +INTEGER :: I +INTEGER, ALLOCATABLE :: JD_TMP(:) +REAL(EB), ALLOCATABLE :: D_TMP(:) + +! Determine current capacity from allocated size +CAP = SIZE(JD) + +! Grow capacity geometrically if needed +IF (NNZ+1 > CAP) THEN + CAP = MAX(2*MAX(CAP,1), NNZ+1) + ALLOCATE(JD_TMP(CAP)) + ALLOCATE(D_TMP(CAP)) + IF (NNZ > 0) THEN + JD_TMP(1:NNZ) = JD(1:NNZ) + D_TMP(1:NNZ) = D(1:NNZ) + ENDIF + IF (ALLOCATED(JD)) DEALLOCATE(JD) + IF (ALLOCATED(D)) DEALLOCATE(D) + CALL MOVE_ALLOC(JD_TMP, JD) + CALL MOVE_ALLOC(D_TMP, D) +ENDIF + +! Find insertion position (ascending order) +POS = 1 +DO WHILE (POS <= NNZ) + IF(JD(POS) < NEWCOL) THEN + POS = POS + 1 + ELSE + EXIT + ENDIF +ENDDO + +! Shift in-place to make room +IF (NNZ >= POS) THEN + DO I = NNZ, POS, -1 + JD(I+1) = JD(I) + D(I+1) = D(I) + ENDDO +ENDIF + +! Insert new entry +JD(POS) = NEWCOL +D(POS) = 0._EB +NNZ = NNZ + 1 + +END SUBROUTINE INSERT_SORTED_INT_REAL + ! --------------------------- GET_MMATRIX_SCALAR_3D ------------------------------- diff --git a/Source/pres.f90 b/Source/pres.f90 index e9ae0fe90fd..d64b1dff083 100644 --- a/Source/pres.f90 +++ b/Source/pres.f90 @@ -2747,7 +2747,7 @@ SUBROUTINE ULMAT_H_MATRIX_SOLVER_SETUP(NM,IPZ) ! Each MPI process builds its local set of rows. ! Matrix blocks defined on CRS distributed format. -! Total number of nonzeros for JD_MAT_H, D_MAT_H: +! Total number of nonzeros: TOT_NNZ_H = SUM( NNZ_H_MAT(1:ZM%NUNKH) ) ! Allocate Solution and RHS vectors: @@ -2766,7 +2766,7 @@ SUBROUTINE ULMAT_H_MATRIX_SOLVER_SETUP(NM,IPZ) IF (ALLOCATED(ZM%JA_H)) DEALLOCATE(ZM%JA_H) ALLOCATE ( ZM%A_H(TOT_NNZ_H) , ZM%IA_H(ZM%NUNKH+1) , ZM%JA_H(TOT_NNZ_H) ) - ! Store upper triangular part of symmetric D_MAT_H in CSR format for the ZONE_MESH (ZM%IA_H,ZM%JA_H,ZM%A_H): + ! Store upper triangular part of symmetric matrix in CSR format for the ZONE_MESH (ZM%IA_H,ZM%JA_H,ZM%A_H): INNZ = 0 DO IROW=1,ZM%NUNKH ZM%IA_H(IROW) = INNZ + 1 @@ -3943,6 +3943,10 @@ SUBROUTINE GET_H_MATRIX_LUDCMP #ifdef WITH_HYPRE INTEGER ::END_ROW, COLOR CHARACTER(FN_LENGTH) :: FN_PARCSRPCG_MATRIX +INTEGER :: NII2, K +INTEGER, ALLOCATABLE :: COLS0(:) +REAL(EB), ALLOCATABLE :: VALS(:) +INTEGER :: ROWARR(1), NCOLSARR(1) #endif !.. All other variables INTEGER :: INNZ, IROW, JCOL, IERR @@ -3982,8 +3986,17 @@ SUBROUTINE GET_H_MATRIX_LUDCMP ! Each MPI process builds its local set of rows. ! Matrix blocks defined on CRS distributed format. - ! Total number of nonzeros for ZSL%JD_MAT_H, ZSL%D_MAT_H: - ZSL%TOT_NNZ_H = 1; IF(ZSL%NUNKH_LOCAL>0) ZSL%TOT_NNZ_H = SUM( ZSL%NNZ_D_MAT_H(1:ZSL%NUNKH_LOCAL) ) + ! Total number of nonzeros: + ZSL%TOT_NNZ_H = 1 + IF (ZSL%NUNKH_LOCAL>0) THEN + IF (ALLOCATED(ZSL%ROW_H)) THEN + ZSL%TOT_NNZ_H = 0 + DO IROW=1,ZSL%NUNKH_LOCAL + ZSL%TOT_NNZ_H = ZSL%TOT_NNZ_H + ZSL%ROW_H(IROW)%NNZ + ENDDO + IF (ZSL%TOT_NNZ_H<1) ZSL%TOT_NNZ_H = 1 + ENDIF + ENDIF ! 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 @@ -4005,16 +4018,18 @@ SUBROUTINE GET_H_MATRIX_LUDCMP !--- 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: + ! Store upper triangular part of symmetric matrix in CSR format: IF(ZSL%NUNKH_LOCAL>0) THEN INNZ = 0 + ZSL%A_H = 0._EB DO IROW=1,ZSL%NUNKH_LOCAL ZSL%IA_H(IROW) = INNZ + 1 - DO JCOL=1,ZSL%NNZ_D_MAT_H(IROW) - IF ( ZSL%JD_MAT_H(JCOL,IROW) < ZSL%UNKH_IND(NM_START)+IROW ) CYCLE ! Only upper Triangular part. + I = ZSL%ROW_H(IROW)%NNZ + DO JCOL=1,I + IF ( ZSL%ROW_H(IROW)%JD(JCOL) < ZSL%UNKH_IND(NM_START)+IROW ) CYCLE ! Only upper triangular part INNZ = INNZ + 1 - ZSL%A_H(INNZ) = ZSL%D_MAT_H(JCOL,IROW) - ZSL%JA_H(INNZ) = ZSL%JD_MAT_H(JCOL,IROW) + ZSL%A_H(INNZ) = ZSL%ROW_H(IROW)%D(JCOL) + ZSL%JA_H(INNZ) = ZSL%ROW_H(IROW)%JD(JCOL) ENDDO ENDDO ZSL%IA_H(ZSL%NUNKH_LOCAL+1) = INNZ + 1 @@ -4124,41 +4139,59 @@ SUBROUTINE GET_H_MATRIX_LUDCMP 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 + ! Assemble from per-row storage + DO IROW=1,ZSL%NUNKH_LOCAL + ZSL%HYPRE_ZSL%INDICES(IROW)=ZSL%LOWER_ROW+IROW-2 + ! Determine row entries and apply indefinite adjustments if needed + NII2 = ZSL%ROW_H(IROW)%NNZ; IF (NII2 < 1) CYCLE + ALLOCATE(COLS0(NII2),VALS(NII2)) + K = 0 + IF (ZSL%MTYPE==SYMM_INDEFINITE) THEN + IF (ZSL%NUNKH_TOTAL==1 .AND. ZSL%NUNKH_LOCAL==1) THEN + DEALLOCATE(COLS0); DEALLOCATE(VALS) + ALLOCATE(COLS0(1)); ALLOCATE(VALS(1)) + COLS0(1) = ZSL%NUNKH_TOTAL-1 + VALS(1) = 1._EB + ROWARR(1) = ZSL%HYPRE_ZSL%INDICES(IROW) + NCOLSARR(1) = 1 + CALL HYPRE_IJMATRIXSETVALUES(ZSL%HYPRE_ZSL%A_H, 1, NCOLSARR, ROWARR, COLS0, VALS, HYPRE_IERR) + DEALLOCATE(COLS0); DEALLOCATE(VALS) + CYCLE ELSE - END_ROW = ZSL%NUNKH_LOCAL + IF(ZSL%UPPER_ROW==ZSL%NUNKH_TOTAL) THEN + END_ROW = ZSL%NUNKH_LOCAL-1 + ELSE + END_ROW = ZSL%NUNKH_LOCAL + ENDIF + IF (IROW <= END_ROW) THEN + DO JCOL=1,NII2 + IF (ZSL%ROW_H(IROW)%JD(JCOL) == ZSL%NUNKH_TOTAL) CYCLE + K = K + 1 + COLS0(K) = ZSL%ROW_H(IROW)%JD(JCOL) - 1 + VALS(K) = ZSL%ROW_H(IROW)%D(JCOL) + ENDDO + ELSEIF (ZSL%UPPER_ROW==ZSL%NUNKH_TOTAL .AND. IROW==ZSL%NUNKH_LOCAL) THEN + DO JCOL=1,NII2 + IF (ZSL%ROW_H(IROW)%JD(JCOL) /= ZSL%NUNKH_TOTAL) CYCLE + K = K + 1 + COLS0(K) = ZSL%ROW_H(IROW)%JD(JCOL) - 1 + VALS(K) = ZSL%ROW_H(IROW)%D(JCOL) + ENDDO + ENDIF 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 + ELSE + DO JCOL=1,NII2 + K = K + 1 + COLS0(K) = ZSL%ROW_H(IROW)%JD(JCOL) - 1 + VALS(K) = ZSL%ROW_H(IROW)%D(JCOL) 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 - - 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) + IF (K > 0) THEN + ROWARR(1) = ZSL%HYPRE_ZSL%INDICES(IROW) + NCOLSARR(1) = K + CALL HYPRE_IJMATRIXSETVALUES(ZSL%HYPRE_ZSL%A_H, 1, NCOLSARR, ROWARR, COLS0(1:K), VALS(1:K), HYPRE_IERR) + ENDIF + DEALLOCATE(COLS0); DEALLOCATE(VALS) ENDDO CALL HYPRE_IJMATRIXASSEMBLE(ZSL%HYPRE_ZSL%A_H, HYPRE_IERR) @@ -4212,10 +4245,8 @@ SUBROUTINE GET_H_MATRIX_LUDCMP 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) + ! Deallocate ROW_H: + IF(ALLOCATED(ZSL%ROW_H)) DEALLOCATE(ZSL%ROW_H) ENDDO IPZ_LOOP @@ -4277,14 +4308,11 @@ SUBROUTINE GET_BCS_H_MATRIX REAL(EB):: AF,IDX,BIJ TYPE(WALL_TYPE), POINTER :: WC TYPE (BOUNDARY_COORD_TYPE), POINTER :: BC -REAL(EB), POINTER, DIMENSION(:,:) :: D_MAT_HP INTEGER :: H_MAT_IVEC IPZ_LOOP : DO IPZ=0,N_ZONE_GLOBMAT ZSL => ZONE_SOLVE(IPZ) - ! Allocate D_MAT_H: - D_MAT_HP => ZSL%D_MAT_H NM1 = NM_START H_MAT_IVEC = 0 @@ -4335,16 +4363,16 @@ SUBROUTINE GET_BCS_H_MATRIX ! Case of unstructured projection: IF(WC%CUT_FACE_INDEX>0) CALL GET_CFACE_OPEN_BC_COEF(WC%CUT_FACE_INDEX,BOUNDARY_COORD(WC%BC_INDEX)%IOR,IDX,BIJ) - ! Find diagonal column number: + ! Find diagonal column number in per-row JD JCOL = -1 - DO JLOC = 1,ZSL%NNZ_D_MAT_H(IND_LOC(LOW_IND)) - IF (IND(LOW_IND) == ZSL%JD_MAT_H(JLOC,IND_LOC(LOW_IND))) THEN + DO JLOC = 1,ZSL%ROW_H(IND_LOC(LOW_IND))%NNZ + IF (IND(LOW_IND) == ZSL%ROW_H(IND_LOC(LOW_IND))%JD(JLOC)) THEN JCOL = JLOC EXIT ENDIF ENDDO ! Add diagonal coefficient due to DIRICHLET BC: - D_MAT_HP(JCOL,IND_LOC(LOW_IND)) = D_MAT_HP(JCOL,IND_LOC(LOW_IND)) + 2._EB*BIJ + ZSL%ROW_H(IND_LOC(LOW_IND))%D(JCOL) = ZSL%ROW_H(IND_LOC(LOW_IND))%D(JCOL) + 2._EB*BIJ ! Add to mesh dirichlet bc counter H_MAT_IVEC = H_MAT_IVEC + 1 ENDDO WALL_LOOP_1 @@ -4373,7 +4401,6 @@ SUBROUTINE GET_H_MATRIX ! Local Variables: INTEGER :: NM,NM1,NREG INTEGER :: LOW_FACE,HIGH_FACE,X1AXIS,X2AXIS,X3AXIS,IFACE -REAL(EB), POINTER, DIMENSION(:,:) :: D_MAT_HP REAL(EB), POINTER, DIMENSION(:) :: DX1,DX2,DX3 INTEGER :: I,J,K,I1,I2,I3,IIP,JJP,KKP,IIM,JJM,KKM INTEGER :: ILOC,JLOC,IROW,JCOL,IND(LOW_IND:HIGH_IND),IND_LOC(LOW_IND:HIGH_IND) @@ -4390,9 +4417,6 @@ SUBROUTINE GET_H_MATRIX ZSL => ZONE_SOLVE(IPZ) - ! Allocate D_MAT_H: - ALLOCATE( ZSL%D_MAT_H(1:NNZ_ROW_H,1:ZSL%NUNKH_LOCAL) ); ZSL%D_MAT_H = 0._EB - D_MAT_HP => ZSL%D_MAT_H NM1 = NM_START ! Main Mesh Loop: @@ -4471,9 +4495,8 @@ SUBROUTINE GET_H_MATRIX DO ILOC = LOW_IND,HIGH_IND ! Local row number in Kface DO JLOC = LOW_IND,HIGH_IND ! Local col number in Kface, JD IROW = IND_LOC(ILOC) ! Unknown number. - JCOL = RGF(IFACE)%JD(ILOC,JLOC); ! Local position of coef in D_MAT_H - ! Add coefficient: - D_MAT_HP(JCOL,IROW) = D_MAT_HP(JCOL,IROW) + KFACE(ILOC,JLOC) + JCOL = RGF(IFACE)%JD(ILOC,JLOC) + ZSL%ROW_H(IROW)%D(JCOL) = ZSL%ROW_H(IROW)%D(JCOL) + KFACE(ILOC,JLOC) ENDDO ENDDO ENDDO IFACE_LOOP_1 @@ -4534,14 +4557,13 @@ SUBROUTINE GET_H_MATRIX ILOC = LOCROW ! Local row number in Kface, only for cell IIG,JJG,KKG. DO JLOC = LOW_IND,HIGH_IND ! Local col number in Kface, JD IROW = IND_LOC(ILOC) ! Unknown number. - JCOL = WC_JD(ILOC,JLOC) ! Local position of coef in D_MAT_H - ! Add coefficient: - D_MAT_HP(JCOL,IROW) = D_MAT_HP(JCOL,IROW) + KFACE(ILOC,JLOC) + JCOL = WC_JD(ILOC,JLOC) + ZSL%ROW_H(IROW)%D(JCOL) = ZSL%ROW_H(IROW)%D(JCOL) + KFACE(ILOC,JLOC) ENDDO ENDDO WALL_LOOP_1 ! Contribution to Laplacian matrix from RC and cut-faces: - IF ( CC_IBM ) CALL GET_H_MATRIX_CC(NM,NM1,IPZ,D_MAT_HP) + IF ( CC_IBM ) CALL GET_H_MATRIX_CC(NM,NM1,IPZ) ENDDO MESH_LOOP_1 ENDDO IPZ_LOOP @@ -4561,8 +4583,7 @@ SUBROUTINE GET_MATRIXGRAPH_H_WHLDOM INTEGER :: NM INTEGER :: X1AXIS,IFACE,I,I1,J,K,IND(LOW_IND:HIGH_IND),IND_LOC(LOW_IND:HIGH_IND) INTEGER :: LOCROW,IIND,NII,ILOC -INTEGER :: NREG,IIM,JJM,KKM,IIP,JJP,KKP,LOW_FACE,HIGH_FACE,JLOC,IW,II,JJ,KK,IIG,JJG,KKG,ICF -LOGICAL :: INLIST +INTEGER :: NREG,IIM,JJM,KKM,IIP,JJP,KKP,LOW_FACE,HIGH_FACE,IW,II,JJ,KK,IIG,JJG,KKG,ICF TYPE(CC_REGFACE_TYPE), POINTER, DIMENSION(:) :: RGF TYPE(CC_RCFACE_TYPE), POINTER :: RCF TYPE(WALL_TYPE), POINTER :: WC @@ -4611,9 +4632,13 @@ SUBROUTINE GET_MATRIXGRAPH_H_WHLDOM ZSL => ZONE_SOLVE(IPZ) ZSL%NUNKH_LOCAL = SUM(ZSL%NUNKH_LOC(LOWER_MESH_INDEX:UPPER_MESH_INDEX)) - ! Allocate NNZ_D_MAT_H, JD_MAT_H: - ALLOCATE( ZSL%NNZ_D_MAT_H(1:ZSL%NUNKH_LOCAL) ); ZSL%NNZ_D_MAT_H(:) = 0 - ALLOCATE( ZSL%JD_MAT_H(1:NNZ_ROW_H,1:ZSL%NUNKH_LOCAL) ); ZSL%JD_MAT_H(:,:) = HUGE(I) + ! Allocate per-row variable storage for graph build (ROW_H): + IF (ALLOCATED(ZSL%ROW_H)) DEALLOCATE(ZSL%ROW_H) + IF (ZSL%NUNKH_LOCAL > 0) THEN + ALLOCATE(ZSL%ROW_H(1:ZSL%NUNKH_LOCAL)) + ELSE + ALLOCATE(ZSL%ROW_H(1:1)) + ENDIF ! Define NM_START: first mesh that belongs to the processor. NM_START = LOWER_MESH_INDEX @@ -4691,29 +4716,7 @@ SUBROUTINE GET_MATRIXGRAPH_H_WHLDOM IND(HIGH_IND) = CCVAR(II,JJ,KK,UNKH) ! guard-cell. IND_LOC(LOW_IND) = IND(LOW_IND) - ZSL%UNKH_IND(NM_START) ! All row indexes must refer to ind_loc. IND_LOC(HIGH_IND)= IND(HIGH_IND)- ZSL%UNKH_IND(NM_START) - ! Same code ADD_INPLACE_NNZ_H_WHLDOM(LOCROW,LOCROW,IND,IND_LOC): - DO IIND=LOW_IND,HIGH_IND - NII = ZSL%NNZ_D_MAT_H(IND_LOC(LOCROW)) - ! Check that column index hasn't been already counted: - INLIST = .FALSE. - DO ILOC=1,NII - IF ( IND(IIND) == ZSL%JD_MAT_H(ILOC,IND_LOC(LOCROW)) ) THEN - INLIST = .TRUE. - EXIT - ENDIF - ENDDO - IF (INLIST) CYCLE - ! Now add in place: - NII = NII + 1 - DO ILOC=1,NII - IF ( ZSL%JD_MAT_H(ILOC,IND_LOC(LOCROW)) > IND(IIND) ) EXIT - ENDDO - DO JLOC=NII,ILOC+1,-1 - ZSL%JD_MAT_H(JLOC,IND_LOC(LOCROW)) = ZSL%JD_MAT_H(JLOC-1,IND_LOC(LOCROW)) - ENDDO - ZSL%NNZ_D_MAT_H(IND_LOC(LOCROW)) = NII - ZSL%JD_MAT_H(ILOC,IND_LOC(LOCROW)) = IND(IIND) - ENDDO + CALL ADD_INPLACE_NNZ_H_WHLDOM(LOCROW,LOCROW,IND,IND_LOC,IPZ) ENDDO WALL_LOOP_1 ! Finally Add nonzeros corresponding to RC_FACE, CUT_FACE @@ -4768,10 +4771,10 @@ SUBROUTINE GET_MATRIXGRAPH_H_WHLDOM RGF(IFACE)%JD(1:2,1:2) = IS_UNDEFINED DO LOCROW = LOW_IND,HIGH_IND DO IIND=LOW_IND,HIGH_IND - NII = ZSL%NNZ_D_MAT_H(IND_LOC(LOCROW)) + NII = ZSL%ROW_H(IND_LOC(LOCROW))%NNZ DO ILOC=1,NII - IF ( IND(IIND) == ZSL%JD_MAT_H(ILOC,IND_LOC(LOCROW)) ) THEN - RGF(IFACE)%JD(LOCROW,IIND) = ILOC; + IF ( IND(IIND) == ZSL%ROW_H(IND_LOC(LOCROW))%JD(ILOC) ) THEN + RGF(IFACE)%JD(LOCROW,IIND) = ILOC EXIT ENDIF ENDDO @@ -4806,9 +4809,9 @@ SUBROUTINE GET_MATRIXGRAPH_H_WHLDOM IND_LOC(LOW_IND) = IND(LOW_IND) - ZSL%UNKH_IND(NM_START) ! All row indexes must refer to ind_loc. IND_LOC(HIGH_IND)= IND(HIGH_IND)- ZSL%UNKH_IND(NM_START) DO IIND=LOW_IND,HIGH_IND - NII = ZSL%NNZ_D_MAT_H(IND_LOC(LOCROW)) + NII = ZSL%ROW_H(IND_LOC(LOCROW))%NNZ DO ILOC=1,NII - IF ( IND(IIND) == ZSL%JD_MAT_H(ILOC,IND_LOC(LOCROW)) ) THEN + IF ( IND(IIND) == ZSL%ROW_H(IND_LOC(LOCROW))%JD(ILOC) ) THEN WC_JD(LOCROW,IIND) = ILOC EXIT ENDIF diff --git a/Source/type.f90 b/Source/type.f90 index 457532e8f21..bbd8fb4111d 100644 --- a/Source/type.f90 +++ b/Source/type.f90 @@ -1047,6 +1047,7 @@ MODULE TYPES INTEGER, ALLOCATABLE, DIMENSION(:) :: ICC_UNKZ_CT_S, ICC_UNKZ_CC_S, ICC_UNKZ_CT_R, ICC_UNKZ_CC_R,& ICF_UFFB_CF_S, ICF_UFFB_CF_R INTEGER, ALLOCATABLE, DIMENSION(:) :: UNKZ_CT_S, UNKZ_CC_S, UNKZ_CT_R, UNKZ_CC_R + INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: MUNKH,GSCH ! Face variables data (velocities): INTEGER, ALLOCATABLE, DIMENSION(:) :: IIO_FC_R,JJO_FC_R,KKO_FC_R,AXS_FC_R,IIO_FC_S,JJO_FC_S,KKO_FC_S,AXS_FC_S @@ -1745,6 +1746,13 @@ MODULE TYPES END TYPE ZONE_MESH_TYPE +!> \brief Per-row container for variable nonzeros in H matrix (used by GLOBMAT_SOLVER) +TYPE H_ROW_TYPE + INTEGER :: NNZ = 0 + INTEGER, ALLOCATABLE, DIMENSION(:) :: JD !< Column indices for this row + REAL(EB), ALLOCATABLE, DIMENSION(:) :: D !< Nonzero values for this row (same length as JD) +END TYPE H_ROW_TYPE + !> \brief Parameters associated with a ZONE, used in GLOBMAT_SOLVER unstructured pressure solver TYPE ZONE_SOLVE_TYPE @@ -1766,14 +1774,12 @@ MODULE TYPES INTEGER :: CONNECTED_ZONE_PARENT=0 !< Index of first zone in a connected zone list. INTEGER, ALLOCATABLE, DIMENSION(:) :: NUNKH_LOC !< Number of local H unknowns per mesh for a given pressure zone. INTEGER, ALLOCATABLE, DIMENSION(:) :: UNKH_IND !< H unknown numbering start index per mesh for a given pressure zone. - INTEGER, ALLOCATABLE, DIMENSION(:) :: NNZ_D_MAT_H!< Number of non-zeros per row of global matrix. - INTEGER, ALLOCATABLE, DIMENSION(:,:) :: JD_MAT_H !< Column index per row of non-zeros of global matrix. - REAL(EB),ALLOCATABLE, DIMENSION(:,:) :: D_MAT_H !< Nonzero values per row for global matrix. INTEGER, ALLOCATABLE, DIMENSION(:,:) :: MESH_IJK !< I,J,K positions of cell with unknown row IROW (1:3,1:NUNKH). REAL(EB),ALLOCATABLE, DIMENSION(:) :: A_H !< Matrix coefficients for pressure zone, up triang part, CSR format. INTEGER ,ALLOCATABLE, DIMENSION(:) :: IA_H,JA_H !< Matrix indexes for pressure zone, up triang part, CSR format. REAL(EB),ALLOCATABLE, DIMENSION(:) :: F_H,X_H !< RHS and Solution containers for pressure zone. REAL(FB),ALLOCATABLE, DIMENSION(:) :: A_H_FB,F_H_FB,X_H_FB!< Arrays in case of single precision solve. + TYPE(H_ROW_TYPE), ALLOCATABLE, DIMENSION(:) :: ROW_H !< Variable-nnz per-row storage for assembly. END TYPE ZONE_SOLVE_TYPE TYPE (ZONE_SOLVE_TYPE), DIMENSION(:), ALLOCATABLE, TARGET :: ZONE_SOLVE