diff --git a/Source/hvac.f90 b/Source/hvac.f90 index c57fe64db9..0546ab6a9c 100644 --- a/Source/hvac.f90 +++ b/Source/hvac.f90 @@ -1532,12 +1532,21 @@ END SUBROUTINE HVAC_CALC SUBROUTINE MATRIX_SOLVE(NNE) USE MATH_FUNCTIONS,ONLY : GAUSSJ INTEGER :: NNE,IERR,ND,NN +REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: LHS_GAUSSJ +REAL(EB), ALLOCATABLE, DIMENSION(:) :: RHS_GAUSSJ TYPE(NETWORK_TYPE), POINTER :: NE TYPE(DUCT_TYPE), POINTER :: DU TYPE(DUCTNODE_TYPE), POINTER :: DN NE =>NETWORK(NNE) -CALL GAUSSJ(LHS(1:NE%N_MATRIX,1:NE%N_MATRIX),NE%N_MATRIX,NE%N_MATRIX,RHS(1:NE%N_MATRIX),1,1,IERR) +ALLOCATE(LHS_GAUSSJ(1:NE%N_MATRIX,1:NE%N_MATRIX)) +LHS_GAUSSJ(1:NE%N_MATRIX,1:NE%N_MATRIX) = LHS(1:NE%N_MATRIX,1:NE%N_MATRIX) +ALLOCATE(RHS_GAUSSJ(1:NE%N_MATRIX)) +RHS_GAUSSJ(1:NE%N_MATRIX) = RHS(1:NE%N_MATRIX) +CALL GAUSSJ(LHS_GAUSSJ,NE%N_MATRIX,NE%N_MATRIX,RHS_GAUSSJ,1,1,IERR) +RHS(1:NE%N_MATRIX) = RHS_GAUSSJ(1:NE%N_MATRIX) +DEALLOCATE(LHS_GAUSSJ) +DEALLOCATE(RHS_GAUSSJ) DO ND = 1,NE%N_DUCTS DU=>DUCT(NE%DUCT_INDEX(ND)) IF (DU%FIXED .OR. DU%AREA < TWENTY_EPSILON_EB) CYCLE @@ -2628,7 +2637,8 @@ END SUBROUTINE DETERMINE_FIXED_ELEMENTS !> \param T Current time (s) SUBROUTINE FIND_NETWORKS(CHANGE,T) -INTEGER:: NZ,NN,ND,DUCT_COUNTER(N_DUCTS),NODE_COUNTER(N_DUCTNODES),COUNTER,COUNTER2,ZONE_COUNTER(N_ZONE) +INTEGER:: NZ,NN,ND,DUCT_COUNTER(N_DUCTS),NODE_COUNTER(N_DUCTNODES),COUNTER,COUNTER2,ZONE_COUNTER(N_ZONE),& + CONNECTED_ZONES2(0:N_ZONE,0:N_ZONE) INTEGER, DIMENSION(:), ALLOCATABLE :: NETWORK_DCOUNTER,NETWORK_NCOUNTER,RENUMBER LOGICAL, INTENT(OUT) :: CHANGE LOGICAL :: CHANGE2 @@ -2639,13 +2649,15 @@ SUBROUTINE FIND_NETWORKS(CHANGE,T) IF (.NOT. ALLOCATED(NETWORK)) CHANGE = .TRUE. IF (N_ZONE > 0) ZONE_COUNTER = 0 - +NODE_COUNTER = 0 +DUCT_COUNTER = 0 +CONNECTED_ZONES2 = CONNECTED_ZONES CALL DETERMINE_FIXED_ELEMENTS(T,CHANGE) DO NN = 1, N_DUCTNODES NZ = DUCTNODE(NN)%ZONE_INDEX IF (NZ>=1) THEN - NODE_COUNTER(NN) = NZ + NODE_COUNTER(NN) = FINDLOC(CONNECTED_ZONES2(NZ,:),1,1)-1 ZONE_COUNTER(NZ) = ZONE_COUNTER(NZ) + 1 ELSE NODE_COUNTER(NN) = NN+N_ZONE @@ -2666,101 +2678,102 @@ SUBROUTINE FIND_NETWORKS(CHANGE,T) ENDDO ENDIF +IF (ALLOCATED(NETWORK)) DEALLOCATE(NETWORK) -IF (CHANGE) THEN - IF (ALLOCATED(NETWORK)) DEALLOCATE(NETWORK) - - CHANGE2 = .TRUE. - DO WHILE (CHANGE2) - CHANGE2 = .FALSE. - DO ND = 1, N_DUCTS - DU => DUCT(ND) - IF (NODE_COUNTER(DU%NODE_INDEX(1)) /= NODE_COUNTER(DU%NODE_INDEX(2))) THEN - CHANGE2 = .TRUE. - COUNTER = MIN(NODE_COUNTER(DU%NODE_INDEX(1)),NODE_COUNTER(DU%NODE_INDEX(2))) - DUCT_COUNTER(ND) = COUNTER - NODE_COUNTER(DU%NODE_INDEX(1)) = COUNTER - NODE_COUNTER(DU%NODE_INDEX(2)) = COUNTER - ELSE - DUCT_COUNTER(ND) = NODE_COUNTER(DU%NODE_INDEX(1)) +CHANGE2 = .TRUE. +DO WHILE (CHANGE2) + CHANGE2 = .FALSE. + DO ND = 1, N_DUCTS + DU => DUCT(ND) + IF (NODE_COUNTER(DU%NODE_INDEX(1)) /= NODE_COUNTER(DU%NODE_INDEX(2))) THEN + CHANGE2 = .TRUE. + IF (NODE_COUNTER(DU%NODE_INDEX(1)) <=N_ZONE .AND. NODE_COUNTER(DU%NODE_INDEX(2)) <=N_ZONE) THEN + CONNECTED_ZONES2(NODE_COUNTER(DU%NODE_INDEX(1)),NODE_COUNTER(DU%NODE_INDEX(2)))=1 + CONNECTED_ZONES2(NODE_COUNTER(DU%NODE_INDEX(2)),NODE_COUNTER(DU%NODE_INDEX(1)))=1 ENDIF - ENDDO + COUNTER = MIN(NODE_COUNTER(DU%NODE_INDEX(1)),NODE_COUNTER(DU%NODE_INDEX(2))) + DUCT_COUNTER(ND) = COUNTER + NODE_COUNTER(DU%NODE_INDEX(1)) = COUNTER + NODE_COUNTER(DU%NODE_INDEX(2)) = COUNTER + ELSE + DUCT_COUNTER(ND) = NODE_COUNTER(DU%NODE_INDEX(1)) + ENDIF + ENDDO - IF (N_ZONE > 0) THEN - DO NZ = 1, N_ZONE - COUNTER = 1 - COUNTER2 = 1 + IF (N_ZONE > 0) THEN + DO NZ = 1, N_ZONE + COUNTER = FINDLOC(CONNECTED_ZONES2(NZ,:),1,1)-1 + DO NN = 1, P_ZONE(NZ)%N_DUCTNODES + COUNTER = MIN(COUNTER,NODE_COUNTER(P_ZONE(NZ)%NODE_INDEX(NN))) + ENDDO + IF (COUNTER < NZ) THEN + COUNTER2 = FINDLOC(CONNECTED_ZONES2(COUNTER,:),1,1)-1 + IF (COUNTER2 /= COUNTER) THEN + CHANGE2 = .TRUE. + COUNTER = COUNTER2 + ENDIF DO NN = 1, P_ZONE(NZ)%N_DUCTNODES - IF (NN==1) THEN - COUNTER = NODE_COUNTER(P_ZONE(NZ)%NODE_INDEX(NN)) - COUNTER2 = COUNTER - ELSE - IF (COUNTER /= NODE_COUNTER(P_ZONE(NZ)%NODE_INDEX(NN))) & - COUNTER2 = MAX(COUNTER2,NODE_COUNTER(P_ZONE(NZ)%NODE_INDEX(NN))) - COUNTER = MIN(COUNTER,NODE_COUNTER(P_ZONE(NZ)%NODE_INDEX(NN))) + IF (NODE_COUNTER(P_ZONE(NZ)%NODE_INDEX(NN)) /= COUNTER) THEN + CHANGE2 = .TRUE. + NODE_COUNTER(P_ZONE(NZ)%NODE_INDEX(NN)) = COUNTER ENDIF ENDDO - IF (COUNTER /= COUNTER2) THEN - CHANGE2 = .TRUE. - DO NN = 1, P_ZONE(NZ)%N_DUCTNODES - NODE_COUNTER(P_ZONE(NZ)%NODE_INDEX(NN)) = COUNTER - ENDDO - ENDIF - ENDDO - ENDIF - END DO + ENDIF + ENDDO + ENDIF +END DO - ALLOCATE(RENUMBER(MAXVAL(NODE_COUNTER))) - RENUMBER = 0 - N_NETWORKS = 0 - DO NN = 1,N_DUCTNODES - IF (RENUMBER(NODE_COUNTER(NN)) == 0) THEN - N_NETWORKS = N_NETWORKS + 1 - RENUMBER(NODE_COUNTER(NN)) = N_NETWORKS - ENDIF - ENDDO - DO NN = 1,N_DUCTNODES - NODE_COUNTER(NN) = RENUMBER(NODE_COUNTER(NN)) - ENDDO - DO ND = 1,N_DUCTS - DUCT_COUNTER(ND) = RENUMBER(DUCT_COUNTER(ND)) - ENDDO - DEALLOCATE(RENUMBER) - - ALLOCATE(NETWORK(N_NETWORKS)) - NETWORK%N_DUCTS=0 - NETWORK%N_DUCTNODES=0 - ALLOCATE(NETWORK_DCOUNTER(N_NETWORKS)) - NETWORK_DCOUNTER=0 - ALLOCATE(NETWORK_NCOUNTER(N_NETWORKS)) - NETWORK_NCOUNTER=0 - COUNTER = 0 - DO ND = 1, N_DUCTS - NETWORK(DUCT_COUNTER(ND))%N_DUCTS = NETWORK(DUCT_COUNTER(ND))%N_DUCTS + 1 - ENDDO - DO NN = 1, N_DUCTNODES - NETWORK(NODE_COUNTER(NN))%N_DUCTNODES = NETWORK(NODE_COUNTER(NN))%N_DUCTNODES + 1 - ENDDO - DO NN = 1, N_NETWORKS - ALLOCATE(NETWORK(NN)%DUCT_INDEX(NETWORK(NN)%N_DUCTS)) - ALLOCATE(NETWORK(NN)%NODE_INDEX(NETWORK(NN)%N_DUCTNODES)) - ALLOCATE(NETWORK(NN)%MATRIX_INDEX(NETWORK(NN)%N_DUCTS+NETWORK(NN)%N_DUCTNODES)) - NETWORK(NN)%MATRIX_INDEX = 0 - ENDDO - DO ND = 1, N_DUCTS - NETWORK_DCOUNTER(DUCT_COUNTER(ND)) = NETWORK_DCOUNTER(DUCT_COUNTER(ND)) + 1 - NETWORK(DUCT_COUNTER(ND))%DUCT_INDEX(NETWORK_DCOUNTER(DUCT_COUNTER(ND))) = ND - DUCT_NE(ND) = NETWORK_DCOUNTER(DUCT_COUNTER(ND)) - ENDDO - DO NN = 1, N_DUCTNODES - NETWORK_NCOUNTER(NODE_COUNTER(NN)) = NETWORK_NCOUNTER(NODE_COUNTER(NN)) + 1 - NETWORK(NODE_COUNTER(NN))%NODE_INDEX(NETWORK_NCOUNTER(NODE_COUNTER(NN))) = NN - DUCTNODE_NE(NN) = NETWORK_NCOUNTER(NODE_COUNTER(NN)) - ENDDO - DEALLOCATE(NETWORK_DCOUNTER) - DEALLOCATE(NETWORK_NCOUNTER) - CALL SETUP_SOLUTION_POINTERS -ENDIF +ALLOCATE(RENUMBER(MAXVAL(NODE_COUNTER))) +RENUMBER = 0 +N_NETWORKS = 0 +DO NN = 1,N_DUCTNODES + IF (RENUMBER(NODE_COUNTER(NN)) == 0) THEN + N_NETWORKS = N_NETWORKS + 1 + RENUMBER(NODE_COUNTER(NN)) = N_NETWORKS + ENDIF +ENDDO +DO NN = 1,N_DUCTNODES + NODE_COUNTER(NN) = RENUMBER(NODE_COUNTER(NN)) +ENDDO +DO ND = 1,N_DUCTS + DUCT_COUNTER(ND) = RENUMBER(DUCT_COUNTER(ND)) +ENDDO +DEALLOCATE(RENUMBER) + +ALLOCATE(NETWORK(N_NETWORKS)) +NETWORK%N_DUCTS=0 +NETWORK%N_DUCTNODES=0 +ALLOCATE(NETWORK_DCOUNTER(N_NETWORKS)) +NETWORK_DCOUNTER=0 +ALLOCATE(NETWORK_NCOUNTER(N_NETWORKS)) +NETWORK_NCOUNTER=0 +COUNTER = 0 +DO ND = 1, N_DUCTS + NETWORK(DUCT_COUNTER(ND))%N_DUCTS = NETWORK(DUCT_COUNTER(ND))%N_DUCTS + 1 +ENDDO +DO NN = 1, N_DUCTNODES + NETWORK(NODE_COUNTER(NN))%N_DUCTNODES = NETWORK(NODE_COUNTER(NN))%N_DUCTNODES + 1 +ENDDO +DO NN = 1, N_NETWORKS + ALLOCATE(NETWORK(NN)%DUCT_INDEX(NETWORK(NN)%N_DUCTS)) + ALLOCATE(NETWORK(NN)%NODE_INDEX(NETWORK(NN)%N_DUCTNODES)) + ALLOCATE(NETWORK(NN)%MATRIX_INDEX(NETWORK(NN)%N_DUCTS+NETWORK(NN)%N_DUCTNODES)) + NETWORK(NN)%MATRIX_INDEX = 0 +ENDDO +DO ND = 1, N_DUCTS + NETWORK_DCOUNTER(DUCT_COUNTER(ND)) = NETWORK_DCOUNTER(DUCT_COUNTER(ND)) + 1 + NETWORK(DUCT_COUNTER(ND))%DUCT_INDEX(NETWORK_DCOUNTER(DUCT_COUNTER(ND))) = ND + DUCT_NE(ND) = NETWORK_DCOUNTER(DUCT_COUNTER(ND)) +ENDDO +DO NN = 1, N_DUCTNODES + NETWORK_NCOUNTER(NODE_COUNTER(NN)) = NETWORK_NCOUNTER(NODE_COUNTER(NN)) + 1 + NETWORK(NODE_COUNTER(NN))%NODE_INDEX(NETWORK_NCOUNTER(NODE_COUNTER(NN))) = NN + DUCTNODE_NE(NN) = NETWORK_NCOUNTER(NODE_COUNTER(NN)) +ENDDO + +DEALLOCATE(NETWORK_DCOUNTER) +DEALLOCATE(NETWORK_NCOUNTER) +CALL SETUP_SOLUTION_POINTERS END SUBROUTINE FIND_NETWORKS @@ -4377,13 +4390,21 @@ SUBROUTINE MATRIX_SOLVE_QFAN(DUCTRUN_INDEX,NF) USE MATH_FUNCTIONS,ONLY : GAUSSJ INTEGER,INTENT(IN) :: DUCTRUN_INDEX,NF INTEGER :: IERR +REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: LHS_GAUSSJ +REAL(EB), ALLOCATABLE, DIMENSION(:) :: RHS_GAUSSJ TYPE(DUCTRUN_TYPE), POINTER :: DR DR =>DUCTRUN(DUCTRUN_INDEX) -CALL GAUSSJ(LHS(1:DR%N_M_DUCTS+DR%N_M_DUCTNODES,1:DR%N_M_DUCTS+DR%N_M_DUCTNODES),& - DR%N_M_DUCTS+DR%N_M_DUCTNODES,DR%N_M_DUCTS+DR%N_M_DUCTNODES,& - RHS(1:DR%N_M_DUCTS+DR%N_M_DUCTNODES),1,1,IERR) +ALLOCATE(LHS_GAUSSJ(1:DR%N_M_DUCTS+DR%N_M_DUCTNODES,1:DR%N_M_DUCTS+DR%N_M_DUCTNODES)) +LHS_GAUSSJ(1:DR%N_M_DUCTS+DR%N_M_DUCTNODES,1:DR%N_M_DUCTS+DR%N_M_DUCTNODES) = & + LHS(1:DR%N_M_DUCTS+DR%N_M_DUCTNODES,1:DR%N_M_DUCTS+DR%N_M_DUCTNODES) +ALLOCATE(RHS_GAUSSJ(1:DR%N_M_DUCTS+DR%N_M_DUCTNODES)) +RHS_GAUSSJ(1:DR%N_M_DUCTS+DR%N_M_DUCTNODES) = RHS(1:DR%N_M_DUCTS+DR%N_M_DUCTNODES) +CALL GAUSSJ(LHS_GAUSSJ,DR%N_M_DUCTS+DR%N_M_DUCTNODES,DR%N_M_DUCTS+DR%N_M_DUCTNODES,RHS_GAUSSJ,1,1,IERR) +RHS(1:DR%N_M_DUCTS+DR%N_M_DUCTNODES) = RHS_GAUSSJ(1:DR%N_M_DUCTS+DR%N_M_DUCTNODES) +DEALLOCATE(LHS_GAUSSJ) +DEALLOCATE(RHS_GAUSSJ) DR%VEL(1:DR%N_M_DUCTS,NF,NEW) =RHS(1:DR%N_M_DUCTS) DR%P(1:DR%N_M_DUCTNODES,NF,NEW) = RHS(DR%N_M_DUCTS+1:DR%N_M_DUCTS+DR%N_M_DUCTNODES)