diff --git a/Source/cons.f90 b/Source/cons.f90 index 7c2f174389..32355a0e76 100644 --- a/Source/cons.f90 +++ b/Source/cons.f90 @@ -583,7 +583,7 @@ MODULE GLOBAL_CONSTANTS ! Divergence Arrays -REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: DSUM,USUM,PSUM +REAL(EB), ALLOCATABLE, DIMENSION(:) :: DSUM,USUM,PSUM ! Level Set vegetation fire spread @@ -633,7 +633,7 @@ MODULE GLOBAL_CONSTANTS ! Number of initial value, pressure zone, and multiplier derived types INTEGER :: N_INIT,N_ZONE,N_MULT,N_MOVE -LOGICAL, ALLOCATABLE, DIMENSION(:,:,:) :: CONNECTED_ZONES +LOGICAL, ALLOCATABLE, DIMENSION(:,:) :: CONNECTED_ZONES INTEGER, ALLOCATABLE, DIMENSION(:,:) :: CONNECTED_ZONES_LOC REAL(EB) :: MINIMUM_ZONE_VOLUME=0._EB REAL(EB) :: PRESSURE_RELAX_TIME=1._EB diff --git a/Source/divg.f90 b/Source/divg.f90 index e7f8b8bc42..a2dc5dca8d 100644 --- a/Source/divg.f90 +++ b/Source/divg.f90 @@ -85,7 +85,7 @@ SUBROUTINE DIVERGENCE_PART_1(T,DT,NM) ! Determine if pressure ZONEs have merged -IF (N_ZONE>0) CALL MERGE_PRESSURE_ZONES(NM) +IF (N_ZONE>0) CALL MERGE_PRESSURE_ZONES ! Compute normal component of velocity at boundaries, U_NORMAL_S in the PREDICTOR step, U_NORMAL in the CORRECTOR. @@ -712,10 +712,6 @@ SUBROUTINE DIVERGENCE_PART_1(T,DT,NM) IF_PRESSURE_ZONES: IF (N_ZONE>0) THEN - USUM(1:N_ZONE,NM) = 0._EB - DSUM(1:N_ZONE,NM) = 0._EB - PSUM(1:N_ZONE,NM) = 0._EB - R_PFCT = 1._EB DO K=1,KBAR DO J=1,JBAR @@ -726,18 +722,18 @@ SUBROUTINE DIVERGENCE_PART_1(T,DT,NM) IF (IPZ<1) CYCLE IF (CELL(CELL_INDEX(I,J,K))%SOLID) CYCLE VC = DX(I)*RC(I)*VC1 - DSUM(IPZ,NM) = DSUM(IPZ,NM) + VC*DP(I,J,K) + DSUM(IPZ) = DSUM(IPZ) + VC*DP(I,J,K) IF (CC_IBM) THEN R_PFCT = 1._EB IF (CCVAR(I,J,K,CC_CGSC) == CC_SOLID) THEN CYCLE ELSEIF(CCVAR(I,J,K,CC_CGSC) == CC_CUTCFE) THEN - CALL ADD_CUTCELL_PSUM(I,J,K,PBAR_P(K,IPZ),PSUM(IPZ,NM)); CYCLE + CALL ADD_CUTCELL_PSUM(I,J,K,PBAR_P(K,IPZ),PSUM(IPZ)); CYCLE ELSEIF(CCVAR(I,J,K,CC_UNKZ) > 0) THEN - CALL ADD_LINKEDCELL_PSUM(I,J,K,VC,PBAR_P(K,IPZ),RTRM(I,J,K),PSUM(IPZ,NM)); CYCLE + CALL ADD_LINKEDCELL_PSUM(I,J,K,VC,PBAR_P(K,IPZ),RTRM(I,J,K),PSUM(IPZ)); CYCLE ENDIF ENDIF - PSUM(IPZ,NM) = PSUM(IPZ,NM) + VC*(R_PBAR(K,IPZ)*R_PFCT-RTRM(I,J,K)) + PSUM(IPZ) = PSUM(IPZ) + VC*(R_PBAR(K,IPZ)*R_PFCT-RTRM(I,J,K)) ENDDO ENDDO ENDDO @@ -752,8 +748,8 @@ SUBROUTINE DIVERGENCE_PART_1(T,DT,NM) IPZ = B1%PRESSURE_ZONE IF (IPZ<1) CYCLE WALL_LOOP4 IF (WC%BOUNDARY_TYPE/=SOLID_BOUNDARY) CYCLE WALL_LOOP4 - IF (PREDICTOR) USUM(IPZ,NM) = USUM(IPZ,NM) + B1%U_NORMAL_S*B1%AREA - IF (CORRECTOR) USUM(IPZ,NM) = USUM(IPZ,NM) + B1%U_NORMAL *B1%AREA + IF (PREDICTOR) USUM(IPZ) = USUM(IPZ) + B1%U_NORMAL_S*B1%AREA + IF (CORRECTOR) USUM(IPZ) = USUM(IPZ) + B1%U_NORMAL *B1%AREA ENDDO WALL_LOOP4 @@ -763,8 +759,8 @@ SUBROUTINE DIVERGENCE_PART_1(T,DT,NM) B1 => BOUNDARY_PROP1(CFA%B1_INDEX) IPZ = B1%PRESSURE_ZONE IF (IPZ<1) CYCLE CFACE_LOOP - IF (PREDICTOR) USUM(IPZ,NM) = USUM(IPZ,NM) + B1%U_NORMAL_S*B1%AREA - IF (CORRECTOR) USUM(IPZ,NM) = USUM(IPZ,NM) + B1%U_NORMAL *B1%AREA + IF (PREDICTOR) USUM(IPZ) = USUM(IPZ) + B1%U_NORMAL_S*B1%AREA + IF (CORRECTOR) USUM(IPZ) = USUM(IPZ) + B1%U_NORMAL *B1%AREA ENDDO CFACE_LOOP @@ -1104,15 +1100,12 @@ SUBROUTINE SPECIES_ADVECTION(N,U_DOT_DEL_RHO_Z) END SUBROUTINE SPECIES_ADVECTION -SUBROUTINE MERGE_PRESSURE_ZONES(NM) +SUBROUTINE MERGE_PRESSURE_ZONES -INTEGER, INTENT(IN) :: NM INTEGER :: IW,IPZ,IOPZ TYPE(WALL_TYPE), POINTER :: WC TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC -CONNECTED_ZONES(:,:,NM) = .FALSE. - DO IW=1,N_EXTERNAL_WALL_CELLS+N_INTERNAL_WALL_CELLS WC=>WALL(IW) IF (WC%BOUNDARY_TYPE/=NULL_BOUNDARY .AND. WC%BOUNDARY_TYPE/=OPEN_BOUNDARY .AND. WC%BOUNDARY_TYPE/=INTERPOLATED_BOUNDARY) CYCLE @@ -1121,15 +1114,15 @@ SUBROUTINE MERGE_PRESSURE_ZONES(NM) IPZ = PRESSURE_ZONE(BC%IIG,BC%JJG,BC%KKG) IOPZ = PRESSURE_ZONE(BC%II,BC%JJ,BC%KK) IF (IW>N_EXTERNAL_WALL_CELLS .AND. IPZ/=IOPZ) THEN - CONNECTED_ZONES(IOPZ,IPZ,NM) = .TRUE. - CONNECTED_ZONES(IPZ,IOPZ,NM) = .TRUE. + CONNECTED_ZONES(IOPZ,IPZ) = .TRUE. + CONNECTED_ZONES(IPZ,IOPZ) = .TRUE. ENDIF IF (WC%BOUNDARY_TYPE==OPEN_BOUNDARY) THEN - CONNECTED_ZONES(0,IPZ,NM) = .TRUE. - CONNECTED_ZONES(IPZ,0,NM) = .TRUE. + CONNECTED_ZONES(0,IPZ) = .TRUE. + CONNECTED_ZONES(IPZ,0) = .TRUE. ELSEIF (WC%BOUNDARY_TYPE==INTERPOLATED_BOUNDARY) THEN - CONNECTED_ZONES(IOPZ,IPZ,NM) = .TRUE. - CONNECTED_ZONES(IPZ,IOPZ,NM) = .TRUE. + CONNECTED_ZONES(IOPZ,IPZ) = .TRUE. + CONNECTED_ZONES(IPZ,IOPZ) = .TRUE. ENDIF ENDDO @@ -1294,36 +1287,36 @@ SUBROUTINE DIVERGENCE_PART_2(DT,NM) USUM_ADD = 0._EB DO IPZ=1,N_ZONE - SUM_P_PSUM = PBAR_P(1,IPZ)*PSUM(IPZ,NM) + SUM_P_PSUM = PBAR_P(1,IPZ)*PSUM(IPZ) OPEN_ZONE = .FALSE. - SUM_USUM = USUM(IPZ,NM) - SUM_DSUM = DSUM(IPZ,NM) - SUM_PSUM = PSUM(IPZ,NM) + SUM_USUM = USUM(IPZ) + SUM_DSUM = DSUM(IPZ) + SUM_PSUM = PSUM(IPZ) DO IOPZ=N_ZONE,0,-1 IF (IOPZ==IPZ) CYCLE - IF (CONNECTED_ZONES(IPZ,IOPZ,NM)) THEN + IF (CONNECTED_ZONES(IPZ,IOPZ)) THEN IF (IOPZ==0) THEN OPEN_ZONE = .TRUE. ELSE - SUM_P_PSUM = SUM_P_PSUM + PBAR_P(1,IOPZ)*PSUM(IOPZ,NM) - SUM_USUM = SUM_USUM + USUM(IOPZ,NM) - SUM_DSUM = SUM_DSUM + DSUM(IOPZ,NM) - SUM_PSUM = SUM_PSUM + PSUM(IOPZ,NM) + SUM_P_PSUM = SUM_P_PSUM + PBAR_P(1,IOPZ)*PSUM(IOPZ) + SUM_USUM = SUM_USUM + USUM(IOPZ) + SUM_DSUM = SUM_DSUM + DSUM(IOPZ) + SUM_PSUM = SUM_PSUM + PSUM(IOPZ) ENDIF ENDIF ENDDO IF (OPEN_ZONE) THEN P_EQ = P_0(1) - USUM_ADD(IPZ) = PSUM(IPZ,NM)*(PBAR_P(1,IPZ)-P_EQ)/PRESSURE_RELAX_TIME + DSUM(IPZ,NM) - USUM(IPZ,NM) + USUM_ADD(IPZ) = PSUM(IPZ)*(PBAR_P(1,IPZ)-P_EQ)/PRESSURE_RELAX_TIME + DSUM(IPZ) - USUM(IPZ) ELSE P_EQ = SUM_P_PSUM/SUM_PSUM - USUM_ADD(IPZ) = PSUM(IPZ,NM)*(PBAR_P(1,IPZ)-P_EQ)/PRESSURE_RELAX_TIME + DSUM(IPZ,NM) - USUM(IPZ,NM) - & - PSUM(IPZ,NM)*(SUM_DSUM-SUM_USUM)/SUM_PSUM + USUM_ADD(IPZ) = PSUM(IPZ)*(PBAR_P(1,IPZ)-P_EQ)/PRESSURE_RELAX_TIME + DSUM(IPZ) - USUM(IPZ) - & + PSUM(IPZ)*(SUM_DSUM-SUM_USUM)/SUM_PSUM ENDIF ENDDO DO IPZ=1,N_ZONE - USUM(IPZ,NM) = USUM(IPZ,NM) + USUM_ADD(IPZ) + USUM(IPZ) = USUM(IPZ) + USUM_ADD(IPZ) ENDDO ! Compute dP/dt for each pressure ZONE @@ -1336,7 +1329,7 @@ SUBROUTINE DIVERGENCE_PART_2(DT,NM) ! Compute change in background pressure DO IPZ=1,N_ZONE - IF (ABS(PSUM(IPZ,NM)) > TWO_EPSILON_EB) D_PBAR_DT_P(IPZ) = (DSUM(IPZ,NM) - USUM(IPZ,NM))/PSUM(IPZ,NM) + IF (ABS(PSUM(IPZ)) > TWO_EPSILON_EB) D_PBAR_DT_P(IPZ) = (DSUM(IPZ) - USUM(IPZ))/PSUM(IPZ) IF (CORRECTOR) P_ZONE(IPZ)%DPSTAR = D_PBAR_DT_P(IPZ) ENDDO diff --git a/Source/hvac.f90 b/Source/hvac.f90 index e9b0b86a72..83d476e9bf 100644 --- a/Source/hvac.f90 +++ b/Source/hvac.f90 @@ -1870,10 +1870,10 @@ SUBROUTINE DPSTARCALC PZ => P_ZONE(IPZ) IF (PZ%N_DUCTNODES==0) CYCLE DPSTAR(IPZ) = P_ZONE(IPZ)%DPSTAR * DT_HV - PSUM_TOT(IPZ) = PSUM(IPZ,1) + PSUM_TOT(IPZ) = PSUM(IPZ) DO IOPZ = 1,N_ZONE IF (IPZ==IOPZ) CYCLE - IF (CONNECTED_ZONES(IPZ,IOPZ,1)) PSUM_TOT(IPZ) = PSUM_TOT(IPZ) + PSUM(IOPZ,1) + IF (CONNECTED_ZONES(IPZ,IOPZ)) PSUM_TOT(IPZ) = PSUM_TOT(IPZ) + PSUM(IOPZ) ENDDO ENDDO @@ -1888,7 +1888,7 @@ SUBROUTINE DPSTARCALC DO IOPZ = 1, N_ZONE IF (IPZ==IOPZ) CYCLE IF (P_ZONE(IOPZ)%N_DUCTNODES==0) CYCLE - IF (CONNECTED_ZONES(IPZ,IOPZ,1)) THEN + IF (CONNECTED_ZONES(IPZ,IOPZ)) THEN IF (P_ZONE(IOPZ)%N_DUCTNODES==0) CYCLE DPSTAR(IOPZ) = DPSTAR(IOPZ) - DN%DIR(1) * DU%AREA * DU%VEL(OLD) * DT_HV/PSUM_TOT(IPZ) IF (DU%FIXED) DPSTAR(IOPZ) = DPSTAR(IOPZ) + DN%DIR(1) * DU%AREA * DU%VEL(NEW) * DT_HV/PSUM_TOT(IPZ) @@ -3049,7 +3049,7 @@ SUBROUTINE COLLAPSE_HVAC_BC(T) WRITE(MESSAGE,'(A,A)') 'ERROR(552): Ductnode must lie with a single pressure zone. Node: ',TRIM(DUCTNODE(NN)%ID) CALL SHUTDOWN(MESSAGE); RETURN ENDIF - IF (ANY(CONNECTED_ZONES(0,DN%ZONE_INDEX,:))) DN%ZONE_INDEX = 0 + IF (CONNECTED_ZONES(0,DN%ZONE_INDEX)) DN%ZONE_INDEX = 0 ENDIF INTERNAL_NODE_IF: IF (((DN%VENT .OR. DN%LEAKAGE) .AND. .NOT. DN%AMBIENT) .OR. & diff --git a/Source/init.f90 b/Source/init.f90 index 6d4fd000a8..b9409cb536 100644 --- a/Source/init.f90 +++ b/Source/init.f90 @@ -1108,7 +1108,7 @@ SUBROUTINE INITIALIZE_MESH_VARIABLES_2(NM) IF (N_ZONE > 0) THEN ZONE_LOOP: DO IPZ = 1,N_ZONE - PSUM(IPZ,NM) = 0._EB + PSUM(IPZ) = 0._EB DO K=1,M%KBAR DO J=1,M%JBAR DO I=1,M%IBAR @@ -1118,7 +1118,7 @@ SUBROUTINE INITIALIZE_MESH_VARIABLES_2(NM) ZZ_GET(1:N_TRACKED_SPECIES) = M%ZZ(I,J,K,1:N_TRACKED_SPECIES) CALL GET_SPECIFIC_HEAT(ZZ_GET,CP,M%TMP(I,J,K)) RTRM = M%RSUM(I,J,K)/(CP*M%PBAR(K,IPZ)) - PSUM(IPZ,NM) = PSUM(IPZ,NM) + VC*(1._EB/M%PBAR(K,IPZ)-RTRM) + PSUM(IPZ) = PSUM(IPZ) + VC*(1._EB/M%PBAR(K,IPZ)-RTRM) ENDDO ENDDO ENDDO diff --git a/Source/main.f90 b/Source/main.f90 index 9da9aa3b24..c3225d0806 100644 --- a/Source/main.f90 +++ b/Source/main.f90 @@ -52,9 +52,8 @@ PROGRAM FDS INTEGER :: LO10,NM,IZERO,ANG_INC_COUNTER REAL(EB) :: T,DT,TNOW REAL :: CPUTIME -REAL(EB), ALLOCATABLE, DIMENSION(:) :: TC_GLB,TC_LOC,DT_NEW,TI_LOC,TI_GLB,DSUM_ALL,PSUM_ALL,USUM_ALL +REAL(EB), ALLOCATABLE, DIMENSION(:) :: TC_GLB,TC_LOC,DT_NEW,TI_LOC,TI_GLB REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: TC2_GLB,TC2_LOC -LOGICAL, ALLOCATABLE, DIMENSION(:,:) :: CONNECTED_ZONES_ALL LOGICAL, ALLOCATABLE, DIMENSION(:) :: STATE_GLB,STATE_LOC INTEGER :: ITER TYPE (MESH_TYPE), POINTER :: M,M4 @@ -440,6 +439,8 @@ PROGRAM FDS ! Compute divergence just in case the flow field is not initialized to ambient +CALL INITIALIZE_DIVERGENCE_INTEGRALS + DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX CALL DIVERGENCE_PART_1(T_BEGIN,DT,NM) ENDDO @@ -672,6 +673,8 @@ PROGRAM FDS ! Boundary conditions for temperature, species, and density. Start divergence calculation. + CALL INITIALIZE_DIVERGENCE_INTEGRALS + COMPUTE_WALL_BC_LOOP_A: DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX CALL WALL_BC(T,DT,NM) IF (PARTICLE_DRAG) CALL PARTICLE_MOMENTUM_TRANSFER(DT,NM) @@ -874,6 +877,8 @@ PROGRAM FDS ! Start the computation of the divergence term. + CALL INITIALIZE_DIVERGENCE_INTEGRALS + DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX IF (N_REACTIONS>0 .OR. INIT_HRRPUV) CALL COMBUSTION_BC(NM) CALL DIVERGENCE_PART_1(T,DT,NM) @@ -1381,14 +1386,7 @@ SUBROUTINE MPI_INITIALIZATION_CHORES(TASK_NUMBER) ! Allocate a few arrays needed to exchange divergence and pressure info among meshes - IF (N_ZONE > 0) THEN - ALLOCATE(DSUM_ALL(N_ZONE),STAT=IZERO) - ALLOCATE(PSUM_ALL(N_ZONE),STAT=IZERO) - ALLOCATE(USUM_ALL(N_ZONE),STAT=IZERO) - ALLOCATE(CONNECTED_ZONES_ALL(0:N_ZONE,0:N_ZONE),STAT=IZERO) - ENDIF - - ALLOCATE(CONNECTED_ZONES(0:N_ZONE,0:N_ZONE,NMESHES),STAT=IZERO) + ALLOCATE(CONNECTED_ZONES(0:N_ZONE,0:N_ZONE),STAT=IZERO) CALL ChkMemErr('INIT','CONNECTED_ZONES',IZERO) CONNECTED_ZONES = .FALSE. @@ -1396,9 +1394,9 @@ SUBROUTINE MPI_INITIALIZATION_CHORES(TASK_NUMBER) CALL ChkMemErr('INIT','CONNECTED_ZONES_LOC',IZERO) CONNECTED_ZONES_LOC = 0 - ALLOCATE(DSUM(N_ZONE,NMESHES),STAT=IZERO) ; CALL ChkMemErr('MAIN','DSUM',IZERO) ; DSUM = 0._EB - ALLOCATE(PSUM(N_ZONE,NMESHES),STAT=IZERO) ; CALL ChkMemErr('MAIN','PSUM',IZERO) ; PSUM = 0._EB - ALLOCATE(USUM(N_ZONE,NMESHES),STAT=IZERO) ; CALL ChkMemErr('MAIN','USUM',IZERO) ; USUM = 0._EB + ALLOCATE(DSUM(N_ZONE),STAT=IZERO) ; CALL ChkMemErr('MAIN','DSUM',IZERO) ; DSUM = 0._EB + ALLOCATE(PSUM(N_ZONE),STAT=IZERO) ; CALL ChkMemErr('MAIN','PSUM',IZERO) ; PSUM = 0._EB + ALLOCATE(USUM(N_ZONE),STAT=IZERO) ; CALL ChkMemErr('MAIN','USUM',IZERO) ; USUM = 0._EB ! Determine if consumable OBST masses are to be exchanged @@ -1806,6 +1804,18 @@ SUBROUTINE END_FDS END SUBROUTINE END_FDS +!> \brief Zero out the integrals found in the divergence expression. + +SUBROUTINE INITIALIZE_DIVERGENCE_INTEGRALS + +USUM = 0._EB +PSUM = 0._EB +DSUM = 0._EB +CONNECTED_ZONES = .FALSE. + +END SUBROUTINE INITIALIZE_DIVERGENCE_INTEGRALS + + !> \brief Exchange information mesh to mesh needed for divergence integrals !> \details Sum DSUM, PSUM and USUM over all meshes controlled by the active process, then reduce over all processes @@ -1817,66 +1827,33 @@ SUBROUTINE EXCHANGE_DIVERGENCE_INFO TNOW = CURRENT_TIME() -CONNECTED_ZONES_ALL = .FALSE. - -DO IPZ=1,N_ZONE - DSUM_ALL(IPZ) = 0._EB - PSUM_ALL(IPZ) = 0._EB - USUM_ALL(IPZ) = 0._EB - DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX - DSUM_ALL(IPZ) = DSUM_ALL(IPZ) + DSUM(IPZ,NM) - PSUM_ALL(IPZ) = PSUM_ALL(IPZ) + PSUM(IPZ,NM) - USUM_ALL(IPZ) = USUM_ALL(IPZ) + USUM(IPZ,NM) - ENDDO -ENDDO -DO IPZ=0,N_ZONE - DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX - DO IOPZ=0,N_ZONE - IF (CONNECTED_ZONES(IPZ,IOPZ,NM)) CONNECTED_ZONES_ALL(IPZ,IOPZ) = .TRUE. - ENDDO - ENDDO -ENDDO +! Sum up the divergence integrals over the MPI processes IF (N_MPI_PROCESSES>1) THEN - CALL MPI_ALLREDUCE(MPI_IN_PLACE,DSUM_ALL(1),N_ZONE,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IERR) - CALL MPI_ALLREDUCE(MPI_IN_PLACE,PSUM_ALL(1),N_ZONE,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IERR) - CALL MPI_ALLREDUCE(MPI_IN_PLACE,USUM_ALL(1),N_ZONE,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IERR) - CALL MPI_ALLREDUCE(MPI_IN_PLACE,CONNECTED_ZONES_ALL(0,0),(N_ZONE+1)**2,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,IERR) + CALL MPI_ALLREDUCE(MPI_IN_PLACE,DSUM,N_ZONE,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IERR) + CALL MPI_ALLREDUCE(MPI_IN_PLACE,PSUM,N_ZONE,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IERR) + CALL MPI_ALLREDUCE(MPI_IN_PLACE,USUM,N_ZONE,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IERR) + CALL MPI_ALLREDUCE(MPI_IN_PLACE,CONNECTED_ZONES,(N_ZONE+1)**2,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,IERR) ENDIF -DO IPZ=1,N_ZONE - DO NM=1,NMESHES - DSUM(IPZ,NM) = DSUM_ALL(IPZ) - PSUM(IPZ,NM) = PSUM_ALL(IPZ) - USUM(IPZ,NM) = USUM_ALL(IPZ) - ENDDO -ENDDO -DO IPZ=0,N_ZONE - DO NM=1,NMESHES - CONNECTED_ZONES(IPZ,:,NM) = CONNECTED_ZONES_ALL(IPZ,:) - ENDDO -ENDDO - ! Connect zones to others which are not directly connected -DO NM=1,NMESHES - NEW_CONNECTION = .TRUE. - DO WHILE (NEW_CONNECTION) - NEW_CONNECTION = .FALSE. - DO IPZ=1,N_ZONE - DO IOPZ=1,N_ZONE - IF (IOPZ==IPZ) CYCLE - IF (CONNECTED_ZONES(IPZ,IOPZ,NM)) THEN - DO IOPZ2=0,N_ZONE - IF (IOPZ==IOPZ2) CYCLE - IF (CONNECTED_ZONES(IOPZ,IOPZ2,NM)) THEN - IF (.NOT.CONNECTED_ZONES(IPZ,IOPZ2,NM)) NEW_CONNECTION = .TRUE. - CONNECTED_ZONES(IPZ,IOPZ2,NM) = .TRUE. - CONNECTED_ZONES(IOPZ2,IPZ,NM) = .TRUE. - ENDIF - ENDDO - ENDIF - ENDDO +NEW_CONNECTION = .TRUE. +DO WHILE (NEW_CONNECTION) + NEW_CONNECTION = .FALSE. + DO IPZ=1,N_ZONE + DO IOPZ=1,N_ZONE + IF (IOPZ==IPZ) CYCLE + IF (CONNECTED_ZONES(IPZ,IOPZ)) THEN + DO IOPZ2=0,N_ZONE + IF (IOPZ==IOPZ2) CYCLE + IF (CONNECTED_ZONES(IOPZ,IOPZ2)) THEN + IF (.NOT.CONNECTED_ZONES(IPZ,IOPZ2)) NEW_CONNECTION = .TRUE. + CONNECTED_ZONES(IPZ,IOPZ2) = .TRUE. + CONNECTED_ZONES(IOPZ2,IPZ) = .TRUE. + ENDIF + ENDDO + ENDIF ENDDO ENDDO ENDDO diff --git a/Source/pres.f90 b/Source/pres.f90 index 1b6d8fbd95..56c673748c 100644 --- a/Source/pres.f90 +++ b/Source/pres.f90 @@ -1186,7 +1186,7 @@ SUBROUTINE ULMAT_SOLVER_SETUP(NM) CONNECTED_ZONES_LOC = 0 DO IOPZ=0,N_ZONE DO IPZ=0,N_ZONE - IF (IPZ==IOPZ .OR. CONNECTED_ZONES(IPZ,IOPZ,NM)) CONNECTED_ZONES_LOC(IPZ,IOPZ) = 1 + IF (IPZ==IOPZ .OR. CONNECTED_ZONES(IPZ,IOPZ)) CONNECTED_ZONES_LOC(IPZ,IOPZ) = 1 ENDDO ENDDO @@ -4832,7 +4832,7 @@ SUBROUTINE GET_MATRIX_INDEXES_H CONNECTED_ZONES_LOC = 0 DO IOPZ=0,N_ZONE DO IPZ=0,N_ZONE - IF (IPZ==IOPZ .OR. CONNECTED_ZONES(IPZ,IOPZ,NM)) CONNECTED_ZONES_LOC(IPZ,IOPZ) = 1 + IF (IPZ==IOPZ .OR. CONNECTED_ZONES(IPZ,IOPZ)) CONNECTED_ZONES_LOC(IPZ,IOPZ) = 1 ENDDO ENDDO ! establish sets of connected zones