diff --git a/Source/ccib.f90 b/Source/ccib.f90 index eafbe5b069..9cc332a5ee 100644 --- a/Source/ccib.f90 +++ b/Source/ccib.f90 @@ -3708,7 +3708,7 @@ SUBROUTINE MESH_CC_EXCHANGE(CODE) REAL(EB), POINTER, DIMENSION(:,:,:) :: UP,UP2,VP,VP2,WP,WP2 LOGICAL, SAVE :: INITIALIZE_CC_SCALARS_FORC=.TRUE. -INTEGER :: EP,INPE,INT_NPE_LO,INT_NPE_HI,VIND,ICELL,IEDGE,IFEP,IW,IIO,JJO,KKO +INTEGER :: EP,INPE,INT_NPE_LO,INT_NPE_HI,VIND,ICELL,IEDGE,IFEP REAL(EB) :: TNOW,TINTP ! For solid phase only return. All variables exchanged currently here are gas-phase. @@ -4509,35 +4509,22 @@ SUBROUTINE MESH_CC_EXCHANGE(CODE) IF (CODE==1 .AND. M2%NICC_R(1)>0) THEN NQT2 = 4+N_TOTAL_SCALARS LL = 0 - ! Copy-cut cell scalar quantities from MESHES(NOM)%OMESH(NM) cells to MESHES(NM) (i.e. other mesh) cut-cells: - ! Use External wall cell loop: - EXTERNAL_WALL_LOOP_1 : DO IW=1,M%N_EXTERNAL_WALL_CELLS - WC=>M%WALL(IW) - EWC=>M%EXTERNAL_WALL(IW) - BC=>M%BOUNDARY_COORD(WC%BC_INDEX) - IF (.NOT.(WC%BOUNDARY_TYPE == INTERPOLATED_BOUNDARY)) CYCLE EXTERNAL_WALL_LOOP_1 - IF (EWC%NOM/=NM) CYCLE EXTERNAL_WALL_LOOP_1 - IF (M%CCVAR(BC%II,BC%JJ,BC%KK,CC_CGSC) /= CC_CUTCFE) CYCLE EXTERNAL_WALL_LOOP_1 - DO KKO=EWC%KKO_MIN,EWC%KKO_MAX - DO JJO=EWC%JJO_MIN,EWC%JJO_MAX - DO IIO=EWC%IIO_MIN,EWC%IIO_MAX - ICC = MESHES(NM)%CCVAR(IIO,JJO,KKO,CC_IDCC) - IF (ICC > 0) THEN - DO JCC=1,MESHES(NM)%CUT_CELL(ICC)%NCELL - LL = LL + 1 - MESHES(NM)%CUT_CELL(ICC)%RHOS(JCC) = M2%REAL_RECV_PKG11(NQT2*(LL-1)+1) - MESHES(NM)%CUT_CELL(ICC)%TMP(JCC) = M2%REAL_RECV_PKG11(NQT2*(LL-1)+2) - MESHES(NM)%CUT_CELL(ICC)%RSUM(JCC) = M2%REAL_RECV_PKG11(NQT2*(LL-1)+3) - MESHES(NM)%CUT_CELL(ICC)%D(JCC) = M2%REAL_RECV_PKG11(NQT2*(LL-1)+4) - DO NN=1,N_TOTAL_SCALARS - MESHES(NM)%CUT_CELL(ICC)%ZZS(NN,JCC) = M2%REAL_RECV_PKG11(NQT2*(LL-1)+4+NN) - ENDDO - ENDDO - ENDIF - ENDDO + ! Unpack from sender-validated ICC_UNKZ_CC_R list (avoids stale CCVAR for ghost cut-cells + ! rejected by the sender's validation pass). + DO ICC1=1,M2%NICC_R(1) + ICC = M2%ICC_UNKZ_CC_R(ICC1) + IF (ICC < 1) CYCLE + DO JCC=1,MESHES(NM)%CUT_CELL(ICC)%NCELL + LL = LL + 1 + MESHES(NM)%CUT_CELL(ICC)%RHOS(JCC) = M2%REAL_RECV_PKG11(NQT2*(LL-1)+1) + MESHES(NM)%CUT_CELL(ICC)%TMP(JCC) = M2%REAL_RECV_PKG11(NQT2*(LL-1)+2) + MESHES(NM)%CUT_CELL(ICC)%RSUM(JCC) = M2%REAL_RECV_PKG11(NQT2*(LL-1)+3) + MESHES(NM)%CUT_CELL(ICC)%D(JCC) = M2%REAL_RECV_PKG11(NQT2*(LL-1)+4) + DO NN=1,N_TOTAL_SCALARS + MESHES(NM)%CUT_CELL(ICC)%ZZS(NN,JCC) = M2%REAL_RECV_PKG11(NQT2*(LL-1)+4+NN) ENDDO ENDDO - ENDDO EXTERNAL_WALL_LOOP_1 + ENDDO ENDIF IF((CODE==1 .OR. CODE==4) .AND. M2%NFCC_R(2)>0) THEN @@ -4632,36 +4619,22 @@ SUBROUTINE MESH_CC_EXCHANGE(CODE) IF (CODE==4 .AND. M2%NICC_R(1)>0) THEN NQT2 = 4+N_TOTAL_SCALARS LL = 0 - ! Copy-cut cell scalar quantities from MESHES(NOM)%OMESH(NM) cells to MESHES(NM) (i.e. other mesh) cut-cells: - ! Use External wall cell loop: - EXTERNAL_WALL_LOOP_2 : DO IW=1,M%N_EXTERNAL_WALL_CELLS - WC=>M%WALL(IW) - IF (.NOT.(WC%BOUNDARY_TYPE == INTERPOLATED_BOUNDARY)) CYCLE EXTERNAL_WALL_LOOP_2 - BC=>M%BOUNDARY_COORD(WC%BC_INDEX) - EWC=>M%EXTERNAL_WALL(IW) - IF (EWC%NOM/=NM) CYCLE EXTERNAL_WALL_LOOP_2 - IF (M%CCVAR(BC%II,BC%JJ,BC%KK,CC_CGSC) /= CC_CUTCFE) & - CYCLE EXTERNAL_WALL_LOOP_2 - DO KKO=EWC%KKO_MIN,EWC%KKO_MAX - DO JJO=EWC%JJO_MIN,EWC%JJO_MAX - DO IIO=EWC%IIO_MIN,EWC%IIO_MAX - ICC = MESHES(NM)%CCVAR(IIO,JJO,KKO,CC_IDCC) - IF (ICC > 0) THEN - DO JCC=1,MESHES(NM)%CUT_CELL(ICC)%NCELL - LL = LL + 1 - MESHES(NM)%CUT_CELL(ICC)%RHO(JCC) = M2%REAL_RECV_PKG11(NQT2*(LL-1)+1) - MESHES(NM)%CUT_CELL(ICC)%TMP(JCC) = M2%REAL_RECV_PKG11(NQT2*(LL-1)+2) - MESHES(NM)%CUT_CELL(ICC)%RSUM(JCC) = M2%REAL_RECV_PKG11(NQT2*(LL-1)+3) - MESHES(NM)%CUT_CELL(ICC)%DS(JCC) = M2%REAL_RECV_PKG11(NQT2*(LL-1)+4) - DO NN=1,N_TOTAL_SCALARS - MESHES(NM)%CUT_CELL(ICC)%ZZ(NN,JCC) = M2%REAL_RECV_PKG11(NQT2*(LL-1)+4+NN) - ENDDO - ENDDO - ENDIF - ENDDO + ! Unpack from sender-validated ICC_UNKZ_CC_R list (avoids stale CCVAR for ghost cut-cells + ! rejected by the sender's validation pass). + DO ICC1=1,M2%NICC_R(1) + ICC = M2%ICC_UNKZ_CC_R(ICC1) + IF (ICC < 1) CYCLE + DO JCC=1,MESHES(NM)%CUT_CELL(ICC)%NCELL + LL = LL + 1 + MESHES(NM)%CUT_CELL(ICC)%RHO(JCC) = M2%REAL_RECV_PKG11(NQT2*(LL-1)+1) + MESHES(NM)%CUT_CELL(ICC)%TMP(JCC) = M2%REAL_RECV_PKG11(NQT2*(LL-1)+2) + MESHES(NM)%CUT_CELL(ICC)%RSUM(JCC) = M2%REAL_RECV_PKG11(NQT2*(LL-1)+3) + MESHES(NM)%CUT_CELL(ICC)%DS(JCC) = M2%REAL_RECV_PKG11(NQT2*(LL-1)+4) + DO NN=1,N_TOTAL_SCALARS + MESHES(NM)%CUT_CELL(ICC)%ZZ(NN,JCC) = M2%REAL_RECV_PKG11(NQT2*(LL-1)+4+NN) ENDDO ENDDO - ENDDO EXTERNAL_WALL_LOOP_2 + ENDDO ENDIF IF (CODE==6 .AND. M2%NFCC_R(1)>0) THEN @@ -4776,13 +4749,12 @@ SUBROUTINE MESH_CC_EXCHANGE(CODE) SUBROUTINE FILL_GCCUTCELL_SPECIES -REAL(EB):: PRFCT, VCELL, RHO_CC, TMP_CC, RSUM_CC, D_CC, ZZ_CC(1:N_TOTAL_SCALARS), VOL -TYPE (OMESH_TYPE), POINTER :: OM +REAL(EB):: VCELL, RHO_CC, TMP_CC, RSUM_CC, D_CC, ZZ_CC(1:N_TOTAL_SCALARS), VOL TYPE(CC_CUTCELL_TYPE), POINTER :: OCC -INTEGER :: NM,NOM,NN,ICC,JCC,IW,IIO,JJO,KKO +INTEGER :: NM,NOM,NN,ICC,ICC1,JCC,IW,IIO,JJO,KKO +LOGICAL :: VALID_SOURCE_CC ! Here inject OMESH cut-cell info obtained in MESH_CC_EXCHANGE into ghost-cell cc containers: -PRFCT = 0._EB; IF (PREDICTOR) PRFCT = 1._EB MESH_LOOP : DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX CALL POINT_TO_MESH(NM) EXTERNAL_WALL_LOOP : DO IW=1,N_EXTERNAL_WALL_CELLS @@ -4793,34 +4765,61 @@ SUBROUTINE FILL_GCCUTCELL_SPECIES IF (CCVAR(BC%II,BC%JJ,BC%KK,CC_CGSC) /= CC_CUTCFE) CYCLE EXTERNAL_WALL_LOOP ! Do volume average to a cell container for ghost cell II,JJ,KK: NOM = EWC%NOM - OM => MESHES(NM)%OMESH(NOM) RHO_CC=0._EB; TMP_CC=0._EB; RSUM_CC=0._EB; D_CC=0._EB; ZZ_CC(1:N_TOTAL_SCALARS)=0._EB; VOL=0._EB DO KKO=EWC%KKO_MIN,EWC%KKO_MAX DO JJO=EWC%JJO_MIN,EWC%JJO_MAX DO IIO=EWC%IIO_MIN,EWC%IIO_MAX IF (MESHES(NOM)%CELL(MESHES(NOM)%CELL_INDEX(IIO,JJO,KKO))%SOLID) CYCLE ICC = MESHES(NOM)%CCVAR(IIO,JJO,KKO,CC_IDCC) - IF (ICC > 0) THEN ! Cut-cells: + ! For remote NOM, only accept the cut-cell entry if it appears in the sender-validated + ! ICC_UNKZ_CC_R list. Otherwise the OMESH cut-cell scalars haven't been filled by + ! MESH_CC_EXCHANGE and reading them would mix in stale/zero state. + VALID_SOURCE_CC = ICC > 0 .AND. (PROCESS(NOM)==MY_RANK .OR. NM==NOM) + IF (ICC > 0 .AND. .NOT.VALID_SOURCE_CC .AND. ALLOCATED(MESHES(NM)%OMESH(NOM)%ICC_UNKZ_CC_R)) THEN + DO ICC1=1,MESHES(NM)%OMESH(NOM)%NICC_R(1) + IF (MESHES(NM)%OMESH(NOM)%ICC_UNKZ_CC_R(ICC1)==ICC) THEN + VALID_SOURCE_CC = .TRUE. + EXIT + ENDIF + ENDDO + ENDIF + IF (VALID_SOURCE_CC) THEN ! Cut-cells: OCC => MESHES(NOM)%CUT_CELL(ICC) DO JCC=1,OCC%NCELL VCELL = OCC%VOLUME(JCC) - RHO_CC = RHO_CC + (PRFCT *OCC%RHOS(JCC) + (1._EB-PRFCT)*OCC%RHO(JCC))*VCELL TMP_CC = TMP_CC + OCC%TMP(JCC)*VCELL RSUM_CC = RSUM_CC + OCC%RSUM(JCC)*VCELL - D_CC = D_CC + (PRFCT *OCC%D(JCC) + (1._EB-PRFCT)*OCC%DS(JCC))*VCELL - DO NN=1,N_TOTAL_SCALARS - ZZ_CC(NN) = ZZ_CC(NN) + (PRFCT *OCC%ZZS(NN,JCC) + (1._EB-PRFCT)*OCC%ZZ(NN,JCC))*VCELL - ENDDO + IF (PREDICTOR) THEN + RHO_CC = RHO_CC + OCC%RHOS(JCC)*VCELL + D_CC = D_CC + OCC%D(JCC)*VCELL + DO NN=1,N_TOTAL_SCALARS + ZZ_CC(NN) = ZZ_CC(NN) + OCC%ZZS(NN,JCC)*VCELL + ENDDO + ELSE + RHO_CC = RHO_CC + OCC%RHO(JCC)*VCELL + D_CC = D_CC + OCC%DS(JCC)*VCELL + DO NN=1,N_TOTAL_SCALARS + ZZ_CC(NN) = ZZ_CC(NN) + OCC%ZZ(NN,JCC)*VCELL + ENDDO + ENDIF VOL = VOL + VCELL ENDDO - ELSEIF(MESHES(NOM)%CCVAR(IIO,JJO,KKO,CC_CGSC) == CC_GASPHASE) THEN ! Regular cell: + ELSEIF(MESHES(NOM)%CCVAR(IIO,JJO,KKO,CC_CGSC) == CC_GASPHASE .OR. ICC > 0) THEN ! Regular cell or unvalidated remote cut-cell. VCELL = MESHES(NOM)%DX(IIO)*MESHES(NOM)%DY(JJO)*MESHES(NOM)%DZ(KKO) - RHO_CC = RHO_CC + (PRFCT*RHOS(BC%II,BC%JJ,BC%KK) + (1._EB-PRFCT)*RHO(BC%II,BC%JJ,BC%KK))*VCELL TMP_CC = TMP_CC + TMP(BC%II,BC%JJ,BC%KK)*VCELL - D_CC = D_CC + (PRFCT*D(BC%II,BC%JJ,BC%KK) + (1._EB-PRFCT)*DS(BC%II,BC%JJ,BC%KK))*VCELL - DO NN=1,N_TOTAL_SCALARS - ZZ_CC(NN)=ZZ_CC(NN)+(PRFCT*ZZS(BC%II,BC%JJ,BC%KK,NN)+(1._EB-PRFCT)*ZZ(BC%II,BC%JJ,BC%KK,NN))*VCELL - ENDDO + IF (PREDICTOR) THEN + RHO_CC = RHO_CC + RHOS(BC%II,BC%JJ,BC%KK)*VCELL + D_CC = D_CC + D(BC%II,BC%JJ,BC%KK)*VCELL + DO NN=1,N_TOTAL_SCALARS + ZZ_CC(NN) = ZZ_CC(NN) + ZZS(BC%II,BC%JJ,BC%KK,NN)*VCELL + ENDDO + ELSE + RHO_CC = RHO_CC + RHO(BC%II,BC%JJ,BC%KK)*VCELL + D_CC = D_CC + DS(BC%II,BC%JJ,BC%KK)*VCELL + DO NN=1,N_TOTAL_SCALARS + ZZ_CC(NN) = ZZ_CC(NN) + ZZ(BC%II,BC%JJ,BC%KK,NN)*VCELL + ENDDO + ENDIF VOL = VOL + VCELL ENDIF ENDDO @@ -19298,7 +19297,7 @@ SUBROUTINE FILL_UNKZ_GUARDCELLS TYPE (MESH_TYPE), POINTER :: M TYPE (OMESH_TYPE), POINTER :: M2,M3 TYPE (MPI_REQUEST), ALLOCATABLE, DIMENSION(:) :: REQ0,REQ0DUM -INTEGER :: N_REQ0, NICC_R, ICC, ICC1, NCELL, JCC, NICF_R, ICF +INTEGER :: N_REQ0, NICC_R, ICC, ICC1, NCELL, JCC, NICF_R, ICF, NICC_VALID, NUNKZ_VALID INTEGER, ALLOCATABLE, DIMENSION(:) :: NCC_SV INTEGER :: ISTR,IEND,JSTR,JEND,KSTR,KEND,IIO,JJO,KKO,IOR,IW,N_INT,IIOF,JJOF,KKOF,X1AXIS LOGICAL :: ALL_FLG @@ -19542,48 +19541,118 @@ SUBROUTINE FILL_UNKZ_GUARDCELLS ELSE M2 => MESHES(NOM)%OMESH(NM) M2%ICC_UNKZ_CC_S(1:3*M2%NICC_S(1)) = M3%ICC_UNKZ_CC_R(1:3*M3%NICC_R(1)) - ! Here write the ICC into ICC_UNKZ_CC_S for OMESH NM: - M => MESHES(NOM); DO ICC1=1,M2%NICC_S(1) - IIO=M2%ICC_UNKZ_CC_S(3*ICC1-2) - JJO=M2%ICC_UNKZ_CC_S(3*ICC1-1) - KKO=M2%ICC_UNKZ_CC_S(3*ICC1 ) - M2%ICC_UNKZ_CC_S(ICC1) = M%CCVAR(IIO,JJO,KKO,CC_IDCC) - ENDDO - ALLOCATE(INT1D(1:M2%NICC_S(1))); INT1D(1:M2%NICC_S(1)) = M2%ICC_UNKZ_CC_S(1:M2%NICC_S(1)) - CALL MOVE_ALLOC(FROM=INT1D,TO=M2%ICC_UNKZ_CC_S) ENDIF ENDDO ENDDO IF ((N_REQ0>0) .AND. (N_MPI_PROCESSES>1)) CALL MPI_WAITALL(N_REQ0,REQ0(1:N_REQ0),MPI_STATUSES_IGNORE,IERR) +DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX + M => MESHES(NM) + DO NOM=1,NMESHES + M3 => MESHES(NM)%OMESH(NOM) + IF (.NOT.ALLOCATED(M3%ICC_UNKZ_CC_S)) CYCLE + NICC_VALID = 0 + NUNKZ_VALID = 0 + DO ICC1=1,M3%NICC_S(1) + IIO = M3%ICC_UNKZ_CC_S(3*ICC1-2) + JJO = M3%ICC_UNKZ_CC_S(3*ICC1-1) + KKO = M3%ICC_UNKZ_CC_S(3*ICC1 ) + ICC = M%CCVAR(IIO,JJO,KKO,CC_IDCC) + IF (ICC<1 .OR. M%CCVAR(IIO,JJO,KKO,CC_CGSC)/=CC_CUTCFE) CYCLE + NICC_VALID = NICC_VALID + 1 + M3%ICC_UNKZ_CC_S(3*NICC_VALID-2:3*NICC_VALID) = (/ IIO, JJO, KKO /) + NUNKZ_VALID = NUNKZ_VALID + M%CUT_CELL(ICC)%NCELL + ENDDO + M3%NICC_S(1) = NICC_VALID + M3%NICC_S(2) = NUNKZ_VALID + ENDDO +ENDDO + +! Exchange the sender-validated cut-cell counts back to the receivers. +N_REQ0 = 0 DO NM=1,NMESHES DO NOM=LOWER_MESH_INDEX,UPPER_MESH_INDEX M2 => MESHES(NOM)%OMESH(NM) - IF (N_MPI_PROCESSES>1 .AND. NM/=NOM .AND. PROCESS(NM)/=MY_RANK .AND. M2%NICC_S(1)>0) THEN - ! Here write the ICC into ICC_UNKZ_CC_S for OMESH NM: - M => MESHES(NOM); DO ICC1=1,M2%NICC_S(1) - IIO=M2%ICC_UNKZ_CC_S(3*ICC1-2) - JJO=M2%ICC_UNKZ_CC_S(3*ICC1-1) - KKO=M2%ICC_UNKZ_CC_S(3*ICC1 ) - M2%ICC_UNKZ_CC_S(ICC1) = M%CCVAR(IIO,JJO,KKO,CC_IDCC) - ENDDO - ALLOCATE(INT1D(1:M2%NICC_S(1))); INT1D(1:M2%NICC_S(1)) = M2%ICC_UNKZ_CC_S(1:M2%NICC_S(1)) - CALL MOVE_ALLOC(FROM=INT1D,TO=M2%ICC_UNKZ_CC_S) + IF (N_MPI_PROCESSES>1 .AND. NM/=NOM .AND. PROCESS(NM)/=MY_RANK .AND. M2%NICC_R(1)>0) THEN + N_REQ0 = N_REQ0 + 1; CALL CHECK_REQ0_SIZE + CALL MPI_IRECV(M2%NICC_R(1),2,MPI_INTEGER,PROCESS(NM),NM,MPI_COMM_WORLD,REQ0(N_REQ0),IERR) ENDIF ENDDO ENDDO +DO NM=1,NMESHES + CALL MPI_BARRIER(MPI_COMM_WORLD,IERR) + IF (PROCESS(NM)/=MY_RANK) CYCLE + DO NOM=1,NMESHES + IF (NM/=NOM .AND. .NOT.MESHES(NM)%CONNECTED_MESH(NOM)) CYCLE + M3 => MESHES(NM)%OMESH(NOM) + IF (.NOT.ALLOCATED(M3%ICC_UNKZ_CC_S)) CYCLE + IF (N_MPI_PROCESSES>1 .AND. NM/=NOM .AND. PROCESS(NOM)/=MY_RANK) THEN + N_REQ0 = N_REQ0 + 1; CALL CHECK_REQ0_SIZE + CALL MPI_ISEND(M3%NICC_S(1),2,MPI_INTEGER,PROCESS(NOM),NM,MPI_COMM_WORLD,REQ0(N_REQ0),IERR) + ELSE + IF(.NOT. ALLOCATED(MESHES(NOM)%OMESH)) CYCLE + M2 => MESHES(NOM)%OMESH(NM) + M2%NICC_R(1:2) = M3%NICC_S(1:2) + ENDIF + ENDDO +ENDDO +IF ((N_REQ0>0) .AND. (N_MPI_PROCESSES>1)) CALL MPI_WAITALL(N_REQ0,REQ0(1:N_REQ0),MPI_STATUSES_IGNORE,IERR) + +! Exchange the compacted cut-cell coordinate list back to the receivers. +N_REQ0 = 0 +DO NM=1,NMESHES + DO NOM=LOWER_MESH_INDEX,UPPER_MESH_INDEX + M2 => MESHES(NOM)%OMESH(NM) + IF (N_MPI_PROCESSES>1 .AND. NM/=NOM .AND. PROCESS(NM)/=MY_RANK .AND. ALLOCATED(M2%ICC_UNKZ_CC_R)) THEN + N_REQ0 = N_REQ0 + 1; CALL CHECK_REQ0_SIZE + CALL MPI_IRECV(M2%ICC_UNKZ_CC_R(1),3*M2%NICC_R(1),MPI_INTEGER,PROCESS(NM),NM,MPI_COMM_WORLD,REQ0(N_REQ0),IERR) + ENDIF + ENDDO +ENDDO +DO NM=1,NMESHES + CALL MPI_BARRIER(MPI_COMM_WORLD,IERR) + IF (PROCESS(NM)/=MY_RANK) CYCLE + DO NOM=1,NMESHES + IF (NM/=NOM .AND. .NOT.MESHES(NM)%CONNECTED_MESH(NOM)) CYCLE + M3 => MESHES(NM)%OMESH(NOM) + IF (.NOT.ALLOCATED(M3%ICC_UNKZ_CC_S)) CYCLE + IF (N_MPI_PROCESSES>1 .AND. NM/=NOM .AND. PROCESS(NOM)/=MY_RANK) THEN + N_REQ0 = N_REQ0 + 1; CALL CHECK_REQ0_SIZE + CALL MPI_ISEND(M3%ICC_UNKZ_CC_S(1),3*M3%NICC_S(1),MPI_INTEGER,PROCESS(NOM),NM,MPI_COMM_WORLD,REQ0(N_REQ0),IERR) + ELSE + IF(.NOT. ALLOCATED(MESHES(NOM)%OMESH)) CYCLE + M2 => MESHES(NOM)%OMESH(NM) + IF (M2%NICC_R(1)>0) M2%ICC_UNKZ_CC_R(1:3*M2%NICC_R(1)) = M3%ICC_UNKZ_CC_S(1:3*M3%NICC_S(1)) + ENDIF + ENDDO +ENDDO +IF ((N_REQ0>0) .AND. (N_MPI_PROCESSES>1)) CALL MPI_WAITALL(N_REQ0,REQ0(1:N_REQ0),MPI_STATUSES_IGNORE,IERR) + +! Convert compacted coordinate lists to local ICC indices on sender and receiver. +DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX + M => MESHES(NM) + DO NOM=1,NMESHES + M3 => MESHES(NM)%OMESH(NOM) + IF (M3%NICC_S(1)<1) CYCLE + DO ICC1=1,M3%NICC_S(1) + IIO = M3%ICC_UNKZ_CC_S(3*ICC1-2) + JJO = M3%ICC_UNKZ_CC_S(3*ICC1-1) + KKO = M3%ICC_UNKZ_CC_S(3*ICC1 ) + M3%ICC_UNKZ_CC_S(ICC1) = M%CCVAR(IIO,JJO,KKO,CC_IDCC) + ENDDO + ENDDO +ENDDO DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX DO NOM=1,NMESHES M3 => MESHES(NM)%OMESH(NOM) IF (M3%NICC_R(1)<1) CYCLE - M => MESHES(NOM); DO ICC1=1,M3%NICC_R(1) - IIO=M3%ICC_UNKZ_CC_R(3*ICC1-2) - JJO=M3%ICC_UNKZ_CC_R(3*ICC1-1) - KKO=M3%ICC_UNKZ_CC_R(3*ICC1 ) + M => MESHES(NOM) + DO ICC1=1,M3%NICC_R(1) + IIO = M3%ICC_UNKZ_CC_R(3*ICC1-2) + JJO = M3%ICC_UNKZ_CC_R(3*ICC1-1) + KKO = M3%ICC_UNKZ_CC_R(3*ICC1 ) M3%ICC_UNKZ_CC_R(ICC1) = M%CCVAR(IIO,JJO,KKO,CC_IDCC) ENDDO - ALLOCATE(INT1D(1:M3%NICC_R(1))); INT1D(1:M3%NICC_R(1)) = M3%ICC_UNKZ_CC_R(1:M3%NICC_R(1)) - CALL MOVE_ALLOC(FROM=INT1D,TO=M3%ICC_UNKZ_CC_R) ENDDO ENDDO @@ -19762,49 +19831,21 @@ SUBROUTINE FILL_UNKZ_GUARDCELLS NCC_SV(:)=0 CALL POINT_TO_MESH(NM) ! Loop over cut-cells: - EXTERNAL_WALL_LOOP_2 : DO IW=1,N_EXTERNAL_WALL_CELLS - WC=>WALL(IW) - EWC=>EXTERNAL_WALL(IW) - BC=>BOUNDARY_COORD(WC%BC_INDEX) - IF (.NOT.(WC%BOUNDARY_TYPE == INTERPOLATED_BOUNDARY .OR. & - WC%BOUNDARY_TYPE == MIRROR_BOUNDARY) ) CYCLE EXTERNAL_WALL_LOOP_2 - IF ( WC%BOUNDARY_TYPE==INTERPOLATED_BOUNDARY ) THEN - IF (CCVAR(BC%II,BC%JJ,BC%KK,CC_CGSC) /= CC_CUTCFE) CYCLE EXTERNAL_WALL_LOOP_2 - NOM = EWC%NOM - DO KKO=EWC%KKO_MIN,EWC%KKO_MAX - DO JJO=EWC%JJO_MIN,EWC%JJO_MAX - DO IIO=EWC%IIO_MIN,EWC%IIO_MAX - ICC = MESHES(NOM)%CCVAR(IIO,JJO,KKO,CC_IDCC) - IF (ICC > 0) THEN - DO JCC=1,MESHES(NOM)%CUT_CELL(ICC)%NCELL - NCC_SV(NOM)=NCC_SV(NOM)+1 - MESHES(NOM)%CUT_CELL(ICC)%UNKZ(JCC) = M%OMESH(NOM)%UNKZ_CC_R(NCC_SV(NOM)) - ENDDO - ENDIF - ENDDO - ENDDO + DO NOM=1,NMESHES + M3 => M%OMESH(NOM) + IF (M3%NICC_R(1)<1) CYCLE + ! Walk the sender-validated ICC_UNKZ_CC_R list directly. This matches exactly the + ! ordering and count used by the sender to fill UNKZ_CC_R, eliminating any reliance + ! on local CCVAR for indexing into the buffer. + DO ICC1=1,M3%NICC_R(1) + ICC = M3%ICC_UNKZ_CC_R(ICC1) + IF (ICC < 1) CYCLE + DO JCC=1,MESHES(NOM)%CUT_CELL(ICC)%NCELL + NCC_SV(NOM)=NCC_SV(NOM)+1 + MESHES(NOM)%CUT_CELL(ICC)%UNKZ(JCC) = M3%UNKZ_CC_R(NCC_SV(NOM)) ENDDO - ELSEIF ( WC%BOUNDARY_TYPE==MIRROR_BOUNDARY ) THEN - IIO = BC%IIG; JJO = BC%JJG; KKO = BC%KKG - IF (CCVAR(BC%II,BC%JJ,BC%KK,CC_CGSC) /= CC_CUTCFE) CYCLE EXTERNAL_WALL_LOOP_2 - ! CYCLE if OBJECT face is in the Mirror Boundary, normal out into ghost-cell: - SELECT CASE(BC%IOR) - CASE( IAXIS); IF(FCVAR(IIO-1,JJO ,KKO ,CC_FGSC,IAXIS) == CC_SOLID) CYCLE EXTERNAL_WALL_LOOP_2 - CASE(-IAXIS); IF(FCVAR(IIO ,JJO ,KKO ,CC_FGSC,IAXIS) == CC_SOLID) CYCLE EXTERNAL_WALL_LOOP_2 - CASE( JAXIS); IF(FCVAR(IIO ,JJO-1,KKO ,CC_FGSC,JAXIS) == CC_SOLID) CYCLE EXTERNAL_WALL_LOOP_2 - CASE(-JAXIS); IF(FCVAR(IIO ,JJO ,KKO ,CC_FGSC,JAXIS) == CC_SOLID) CYCLE EXTERNAL_WALL_LOOP_2 - CASE( KAXIS); IF(FCVAR(IIO ,JJO ,KKO-1,CC_FGSC,KAXIS) == CC_SOLID) CYCLE EXTERNAL_WALL_LOOP_2 - CASE(-KAXIS); IF(FCVAR(IIO ,JJO ,KKO ,CC_FGSC,KAXIS) == CC_SOLID) CYCLE EXTERNAL_WALL_LOOP_2 - END SELECT - NOM = NM; ICC = MESHES(NOM)%CCVAR(IIO,JJO,KKO,CC_IDCC) - IF (ICC > 0) THEN - DO JCC=1,MESHES(NOM)%CUT_CELL(ICC)%NCELL - NCC_SV(NOM)=NCC_SV(NOM)+1 - MESHES(NOM)%CUT_CELL(ICC)%UNKZ(JCC) = M%OMESH(NOM)%UNKZ_CC_R(NCC_SV(NOM)) - ENDDO - ENDIF - ENDIF - ENDDO EXTERNAL_WALL_LOOP_2 + ENDDO + ENDDO ENDDO DEALLOCATE(NCC_SV) @@ -22237,10 +22278,13 @@ SUBROUTINE NUMBER_UNKH_CUTCELLS(FLAG12,NM,IPZ,NUNKH_LC) INTEGER :: ICC, JCC, I, J, K FLAG12_COND : IF (FLAG12) THEN - ! Initialize Cut-cell unknown numbers as undefined. - DO ICC=1,MESHES(NM)%N_CUTCELL_MESH - CUT_CELL(ICC)%UNKH(:) = CC_UNDEFINED - ENDDO + ! Initialize once before the pressure-zone loop assigns per-zone cut-cell unknowns. + ! Reinitializing on each IPZ pass would erase unknowns assigned to earlier zones. + IF (IPZ==0) THEN + DO ICC=1,MESHES(NM)%N_CUTCELL_MESH + CUT_CELL(ICC)%UNKH(:) = CC_UNDEFINED + ENDDO + ENDIF DO ICC=1,MESHES(NM)%N_CUTCELL_MESH CC => CUT_CELL(ICC); I=CC%IJK(IAXIS); J=CC%IJK(JAXIS); K=CC%IJK(KAXIS); IF(CELL(CELL_INDEX(I,J,K))%SOLID) CYCLE IF(ZONE_SOLVE(PRESSURE_ZONE(I,J,K))%CONNECTED_ZONE_PARENT/=IPZ ) CYCLE diff --git a/Source/geom.f90 b/Source/geom.f90 index cac28b4cd2..8671fd81a0 100644 --- a/Source/geom.f90 +++ b/Source/geom.f90 @@ -754,7 +754,8 @@ SUBROUTINE SET_CUTCELLS_3D LOGICAL, ALLOCATABLE, DIMENSION(:) :: CC_COMPUTE_MESH, CC_COMPUTE_MESH_AUX REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: GEOM_ZMAX_AUX -INTEGER :: IW,II,JJ,IIF,JJF,KKF,IIOF,JJOF,KKOF,LOHIF,IOR,CT,NCFACE_CUTCELL,NFACE_CELL,AX,SIDE,ICC,JCC,ICFC,IFC +INTEGER :: IW,II,JJ,IIF,JJF,KKF,IIOF,JJOF,KKOF,LOHIF,IOR,CT,NCFACE_CUTCELL,NFACE_CELL,AX,SIDE,ICC,JCC,ICFC,IFC, & + IBOD_DONOR,ITRI_DONOR TYPE(MESH_TYPE), POINTER :: M, M2 TYPE(EXTERNAL_WALL_TYPE), POINTER :: EWC TYPE(WALL_TYPE), POINTER :: WC @@ -767,6 +768,7 @@ SUBROUTINE SET_CUTCELLS_3D INTEGER, PARAMETER :: ADDJ(LOW_IND:HIGH_IND,IAXIS:KAXIS) = RESHAPE((/ 0,0,-1,0, 0,0/),(/2,3/)) INTEGER, PARAMETER :: ADDK(LOW_IND:HIGH_IND,IAXIS:KAXIS) = RESHAPE((/ 0,0, 0,0,-1,0/),(/2,3/)) INTEGER :: IIO,JJO,KKO,IOGC,JOGC,KOGC +LOGICAL :: FM_PENDING_BLOCK_SCAN(1:NMESHES) REAL(EB) :: TNOW @@ -1492,6 +1494,14 @@ SUBROUTINE SET_CUTCELLS_3D IF (ALLOCATED(SPCELLS_TO_BLOCK)) DEALLOCATE(SPCELLS_TO_BLOCK) ENDDO MAIN_MESH_LOOP +! Promote refinement footprints from initially blocked fine cells (e.g., slivers tagged +! BLOCKED_SMALL_CELL in GET_CARTCELL_CUTCELLS / GET_CELL_LINK_INFO): +DO NM=1,NMESHES + IF (.NOT.CC_COMPUTE_MESH(NM)) CYCLE ! Only MESHES assigned to processor and OMESHES of these. + IF (PERIODIC_TEST==105 .AND. PROCESS(NM)/=MY_RANK) CYCLE ! Don't do OMESHES for PERIODIC_TEST==105 + CALL PROMOTE_REFINEMENT_FOOTPRINTS_FROM_BLOCKED_FINE(NM) +ENDDO + CALL DEALLOCATE_BODINT_PLANE(BODINT_PLANE) CALL DEALLOCATE_BODINT_PLANE(BODINT_PLANE2) @@ -1500,6 +1510,8 @@ SUBROUTINE SET_CUTCELLS_3D DO IDIM=1,MAX_DIM +FM_PENDING_BLOCK_SCAN = .FALSE. + ! Exchange CC%NOADVANCE(JCC)>0 information among NEIGHBOURING meshes: CALL EXCHANGE_CC_NOADVANCE_INFO ! Add CC%NOADVANCE(JCC) where needed: @@ -1575,9 +1587,9 @@ SUBROUTINE SET_CUTCELLS_3D CF => M%CUT_FACE(ICF1) IIF=CF%IJK(IAXIS); JJF=CF%IJK(JAXIS); KKF=CF%IJK(KAXIS) SELECT CASE(CF%IJK(KAXIS+1)) - CASE(IAXIS); ACRT = DY(JJF)*DZ(KKF) - CASE(JAXIS); ACRT = DZ(KKF)*DX(IIF) - CASE(KAXIS); ACRT = DX(IIF)*DY(JJF) + CASE(IAXIS); ACRT = DYFACE(JJF)*DZFACE(KKF) + CASE(JAXIS); ACRT = DZFACE(KKF)*DXFACE(IIF) + CASE(KAXIS); ACRT = DXFACE(IIF)*DYFACE(JJF) END SELECT IF(SUM(CF%AREA(1:CF%NFACE))/ACRT>=CCVOL_LINK) CYCLE CC%NOADVANCE(J)=BLOCKED_CAVITY_CELL @@ -1585,38 +1597,44 @@ SUBROUTINE SET_CUTCELLS_3D ENDDO ENDDO IF (K>0) THEN + FM_PENDING_BLOCK_SCAN(NM) = .TRUE. + ENDIF + CALL DEFINE_XYZFACE_CELL(ALLOC_FLG=.FALSE.) +ENDDO MAIN_MESH_LOOP_1 + +IF (ANY(FM_PENDING_BLOCK_SCAN)) THEN + DO WHILE (ANY(FM_PENDING_BLOCK_SCAN)) + DO NM=1,NMESHES + IF (.NOT.FM_PENDING_BLOCK_SCAN(NM)) CYCLE + FM_PENDING_BLOCK_SCAN(NM) = .FALSE. + CALL PROMOTE_REFINEMENT_FOOTPRINTS_FROM_BLOCKED_FINE(NM) + ENDDO + ENDDO + MAIN_MESH_LOOP_1B : DO NM=1,NMESHES + IF (.NOT.CC_COMPUTE_MESH(NM)) CYCLE ! Only MESHES assigned to processor and OMESHES of these. + IF (PERIODIC_TEST==105 .AND. PROCESS(NM)/=MY_RANK) CYCLE ! Don't do OMESHES for PERIODIC_TEST==105 + + CALL POINT_TO_MESH(NM) + M => MESHES(NM) + CALL DEFINE_XYZFACE_CELL(ALLOC_FLG=.TRUE.) + CALL GET_CC_FACE_CELL_LIST_INFO(NM,PHASE=2) + CALL GET_CELL_LINK_INFO(NM) CALL BLOCK_SMALL_UNLINKED_CUTCELLS(NM,SUM_CCELL) IF(SUM_CCELL>0) THEN ! Rebuild incidences and cell linking information: CALL GET_CC_FACE_CELL_LIST_INFO(NM,PHASE=2) CALL GET_CELL_LINK_INFO(NM) ENDIF - ENDIF - CALL DEFINE_XYZFACE_CELL(ALLOC_FLG=.FALSE.) -ENDDO MAIN_MESH_LOOP_1 + CALL DEFINE_XYZFACE_CELL(ALLOC_FLG=.FALSE.) + ENDDO MAIN_MESH_LOOP_1B +ENDIF +CALL EXCHANGE_CC_NOADVANCE_INFO +CALL ADD_NEIGHBOR_BLOCKED_CELLS ! Call tag boundary cut-cells for blocking in refinement interfaces: CALL TAG_CC_BLOCKING_REFINEMENT ENDDO -! WRITE A file per process and mesh with the NOADVANCE cut-cells: -! DO NM=1,NMESHES -! M => MESHES(NM) -! IF(N_MPI_PROCESSES>1) THEN -! WRITE(VERBOSE_FILE,'(A,A,I0,A,I0,A)') TRIM(CHID),'_NOADVANCE_',MY_RANK,'_',NM,'.2log' -! ELSE -! WRITE(VERBOSE_FILE,'(A,A,I0,A,I0,A)') TRIM(CHID),'_NOADVANCE_',MY_RANK,'_',NM,'.1log' -! ENDIF -! OPEN(UNIT=787,FILE=TRIM(VERBOSE_FILE),STATUS='UNKNOWN') -! DO ICC=1,M%N_CUTCELL_MESH+M%N_GCCUTCELL_MESH -! CC=>M%CUT_CELL(ICC) -! DO JCC=1,CC%NCELL -! IF(CC%NOADVANCE(JCC)>0) WRITE(787,*) 'B',NM,';',ICC,JCC,CC%IJK(IAXIS:KAXIS),CC%NCELL -! ENDDO -! ENDDO -! CLOSE(787) -! ENDDO - MAIN_MESH_LOOP_3 : DO NM=1,NMESHES IF (.NOT.CC_COMPUTE_MESH(NM)) CYCLE ! Only MESHES assigned to processor and OMESHES of these. @@ -1752,6 +1770,14 @@ SUBROUTINE SET_CUTCELLS_3D MESHES(NM)%N_CC_BLOCKED = 0 IF(ALLOCATED(MESHES(NM)%XYZ_CC_BLOCKED)) DEALLOCATE(MESHES(NM)%XYZ_CC_BLOCKED) IF(ALLOCATED(MESHES(NM)%JBT_CC_BLOCKED)) DEALLOCATE(MESHES(NM)%JBT_CC_BLOCKED) + ! BODTRI_DONOR is consumed only during the setup blocking / refinement-interface passes + ! above. Free it here so it doesn't sit allocated for the lifetime of the run. + IF (ALLOCATED(MESHES(NM)%CUT_CELL)) THEN + DO ICC=1,MESHES(NM)%N_CUTCELL_MESH+MESHES(NM)%N_GCCUTCELL_MESH + IF (ALLOCATED(MESHES(NM)%CUT_CELL(ICC)%BODTRI_DONOR)) & + DEALLOCATE(MESHES(NM)%CUT_CELL(ICC)%BODTRI_DONOR) + ENDDO + ENDIF ENDDO ! Finally allocate Face and cell variables, compute area and volume factors: @@ -2226,11 +2252,294 @@ SUBROUTINE SET_CUTCELLS_3D CONTAINS +SUBROUTINE PROMOTE_REFINEMENT_FOOTPRINTS_FROM_BLOCKED_FINE(NM_FINE) + +INTEGER, INTENT(IN) :: NM_FINE +INTEGER :: NM_COARSE,IW_LOC,II_COARSE,JJ_COARSE,KK_COARSE,IOR_COARSE, & + I_FINE,J_FINE,K_FINE,IIO_LOC,JJO_LOC,KKO_LOC,IBOD_DONOR,ITRI_DONOR +LOGICAL :: FINE_BLOCKED,CELL_CHANGED +TYPE(WALL_TYPE), POINTER :: WC_COARSE +TYPE(EXTERNAL_WALL_TYPE), POINTER :: EWC_COARSE +TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC_COARSE +TYPE(MESH_TYPE), POINTER :: MC + +COARSE_MESH_LOOP : DO NM_COARSE=1,NMESHES + IF (PROCESS(NM_COARSE)/=MY_RANK) CYCLE COARSE_MESH_LOOP + MC => MESHES(NM_COARSE) + COARSE_WALL_LOOP : DO IW_LOC=1,MC%N_EXTERNAL_WALL_CELLS + WC_COARSE => MC%WALL(IW_LOC); IF (WC_COARSE%BOUNDARY_TYPE/=INTERPOLATED_BOUNDARY) CYCLE COARSE_WALL_LOOP + EWC_COARSE => MC%EXTERNAL_WALL(IW_LOC); IF (EWC_COARSE%NOM/=NM_FINE) CYCLE COARSE_WALL_LOOP + IF ((EWC_COARSE%IIO_MAX-EWC_COARSE%IIO_MIN+1) * & + (EWC_COARSE%JJO_MAX-EWC_COARSE%JJO_MIN+1) * & + (EWC_COARSE%KKO_MAX-EWC_COARSE%KKO_MIN+1) <= 1) CYCLE COARSE_WALL_LOOP + BC_COARSE => MC%BOUNDARY_COORD(WC_COARSE%BC_INDEX) + II_COARSE = BC_COARSE%IIG; JJ_COARSE = BC_COARSE%JJG; KK_COARSE = BC_COARSE%KKG; IOR_COARSE = BC_COARSE%IOR + + FINE_BLOCKED = .FALSE.; IBOD_DONOR = 0; ITRI_DONOR = 0 + DO KKO_LOC=EWC_COARSE%KKO_MIN,EWC_COARSE%KKO_MAX + DO JJO_LOC=EWC_COARSE%JJO_MIN,EWC_COARSE%JJO_MAX + DO IIO_LOC=EWC_COARSE%IIO_MIN,EWC_COARSE%IIO_MAX + CALL GET_FINE_CELL_FROM_COARSE_WALL(IOR_COARSE,IIO_LOC,JJO_LOC,KKO_LOC,I_FINE,J_FINE,K_FINE) + IF (.NOT.CELL_HAS_BLOCKED_CUTCELL(NM_FINE,I_FINE,J_FINE,K_FINE)) CYCLE + FINE_BLOCKED = .TRUE. + IF (IBOD_DONOR<1 .OR. ITRI_DONOR<1) & + CALL GET_REFINEMENT_CELL_DONOR(NM_FINE,I_FINE,J_FINE,K_FINE,IBOD_DONOR,ITRI_DONOR) + ENDDO + ENDDO + ENDDO + IF (.NOT.FINE_BLOCKED) CYCLE COARSE_WALL_LOOP + CALL TAG_CELL_BLOCKED_BY_REFINEMENT_FOOTPRINT(NM_COARSE,II_COARSE,JJ_COARSE,KK_COARSE,BLOCKED_REFI_INTER, & + IBOD_DONOR,ITRI_DONOR,CELL_CHANGED) + IF (CELL_CHANGED) FM_PENDING_BLOCK_SCAN(NM_COARSE) = .TRUE. + CALL TAG_ALL_FINE_FOOTPRINTS_FOR_COARSE_CELL(NM_COARSE,II_COARSE,JJ_COARSE,KK_COARSE,IOR_COARSE, & + BLOCKED_REFI_INTER,IBOD_DONOR,ITRI_DONOR) + ENDDO COARSE_WALL_LOOP +ENDDO COARSE_MESH_LOOP +END SUBROUTINE PROMOTE_REFINEMENT_FOOTPRINTS_FROM_BLOCKED_FINE + +SUBROUTINE TAG_ALL_FINE_FOOTPRINTS_FOR_COARSE_CELL(NM_COARSE,II_COARSE,JJ_COARSE,KK_COARSE,IOR_TRIGGER,BLOCK_TAG, & + IBOD_DONOR,ITRI_DONOR) +INTEGER, INTENT(IN) :: NM_COARSE,II_COARSE,JJ_COARSE,KK_COARSE,IOR_TRIGGER,BLOCK_TAG,IBOD_DONOR,ITRI_DONOR +INTEGER :: IW_LOC,NM_FINE +TYPE(WALL_TYPE), POINTER :: WC_COARSE +TYPE(EXTERNAL_WALL_TYPE), POINTER :: EWC_COARSE +TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC_COARSE +TYPE(MESH_TYPE), POINTER :: MC + +MC => MESHES(NM_COARSE); IF (.NOT.ALLOCATED(MC%WALL)) RETURN +DO IW_LOC=1,MC%N_EXTERNAL_WALL_CELLS + WC_COARSE => MC%WALL(IW_LOC); IF (WC_COARSE%BOUNDARY_TYPE/=INTERPOLATED_BOUNDARY) CYCLE + BC_COARSE => MC%BOUNDARY_COORD(WC_COARSE%BC_INDEX) + IF (BC_COARSE%IOR/=IOR_TRIGGER) CYCLE + IF (BC_COARSE%IIG/=II_COARSE .OR. BC_COARSE%JJG/=JJ_COARSE .OR. BC_COARSE%KKG/=KK_COARSE) CYCLE + EWC_COARSE => MC%EXTERNAL_WALL(IW_LOC) + NM_FINE = EWC_COARSE%NOM + CALL TAG_FINE_CELLS_IN_COARSE_CELL_VOLUME(NM_COARSE,II_COARSE,JJ_COARSE,KK_COARSE,NM_FINE,BLOCK_TAG, & + IBOD_DONOR,ITRI_DONOR) +ENDDO +END SUBROUTINE TAG_ALL_FINE_FOOTPRINTS_FOR_COARSE_CELL + +SUBROUTINE TAG_FINE_CELLS_IN_COARSE_CELL_VOLUME(NM_COARSE,I_COARSE,J_COARSE,K_COARSE,NM_FINE,BLOCK_TAG,IBOD_DONOR,ITRI_DONOR) +INTEGER, INTENT(IN) :: NM_COARSE,I_COARSE,J_COARSE,K_COARSE,NM_FINE,BLOCK_TAG,IBOD_DONOR,ITRI_DONOR +INTEGER :: I_FINE,J_FINE,K_FINE +REAL(EB) :: XLO,XHI,YLO,YHI,ZLO,ZHI,EPS_LOC +LOGICAL :: CELL_CHANGED +TYPE(MESH_TYPE), POINTER :: MC,MF + +IF (PROCESS(NM_FINE)/=MY_RANK) RETURN + +MC => MESHES(NM_COARSE) +MF => MESHES(NM_FINE) +XLO = MC%X(I_COARSE-1); XHI = MC%X(I_COARSE) +YLO = MC%Y(J_COARSE-1); YHI = MC%Y(J_COARSE) +ZLO = MC%Z(K_COARSE-1); ZHI = MC%Z(K_COARSE) +EPS_LOC = 10._EB*GEOMEPS +DO K_FINE=0,MF%KBP1 + IF (MF%ZC(K_FINE)ZHI+EPS_LOC) CYCLE + DO J_FINE=0,MF%JBP1 + IF (MF%YC(J_FINE)YHI+EPS_LOC) CYCLE + DO I_FINE=0,MF%IBP1 + IF (I_FINE>=1 .AND. I_FINE<=MF%IBAR .AND. & + J_FINE>=1 .AND. J_FINE<=MF%JBAR .AND. & + K_FINE>=1 .AND. K_FINE<=MF%KBAR) CYCLE + IF (MF%XC(I_FINE)XHI+EPS_LOC) CYCLE + CALL TAG_CELL_BLOCKED_BY_REFINEMENT_FOOTPRINT(NM_FINE,I_FINE,J_FINE,K_FINE,BLOCK_TAG, & + IBOD_DONOR,ITRI_DONOR,CELL_CHANGED) + IF (CELL_CHANGED) FM_PENDING_BLOCK_SCAN(NM_FINE) = .TRUE. + ENDDO + ENDDO +ENDDO +END SUBROUTINE TAG_FINE_CELLS_IN_COARSE_CELL_VOLUME + +SUBROUTINE GET_FINE_CELL_FROM_COARSE_WALL(IOR_COARSE,IIO_LOC,JJO_LOC,KKO_LOC,I_FINE,J_FINE,K_FINE) + +INTEGER, INTENT(IN) :: IOR_COARSE,IIO_LOC,JJO_LOC,KKO_LOC +INTEGER, INTENT(OUT) :: I_FINE,J_FINE,K_FINE + +I_FINE = IIO_LOC; J_FINE = JJO_LOC; K_FINE = KKO_LOC +SELECT CASE(ABS(IOR_COARSE)) +CASE(IAXIS); I_FINE = I_FINE + SIGN(1,IOR_COARSE) +CASE(JAXIS); J_FINE = J_FINE + SIGN(1,IOR_COARSE) +CASE(KAXIS); K_FINE = K_FINE + SIGN(1,IOR_COARSE) +END SELECT +END SUBROUTINE GET_FINE_CELL_FROM_COARSE_WALL + +LOGICAL FUNCTION CELL_HAS_BLOCKED_CUTCELL(NM_LOC,I_LOC,J_LOC,K_LOC) +INTEGER, INTENT(IN) :: NM_LOC,I_LOC,J_LOC,K_LOC +INTEGER :: ICC_LOC +CELL_HAS_BLOCKED_CUTCELL = .FALSE. +ICC_LOC = MESHES(NM_LOC)%CCVAR(I_LOC,J_LOC,K_LOC,CC_IDCC); IF (ICC_LOC<1) RETURN +CELL_HAS_BLOCKED_CUTCELL = ANY(MESHES(NM_LOC)%CUT_CELL(ICC_LOC)%NOADVANCE(1:MESHES(NM_LOC)%CUT_CELL(ICC_LOC)%NCELL)>0) +END FUNCTION CELL_HAS_BLOCKED_CUTCELL + +SUBROUTINE GET_REFINEMENT_CELL_DONOR(NM_LOC,I_LOC,J_LOC,K_LOC,IBOD_OUT,ITRI_OUT) + +INTEGER, INTENT(IN) :: NM_LOC,I_LOC,J_LOC,K_LOC +INTEGER, INTENT(OUT) :: IBOD_OUT,ITRI_OUT +INTEGER :: II,JJ,KK,ICC_LOC,JCC_LOC,IFC,IFACE,IFC1,JFC1,COUNT,COUNT_MAX,DUM,PASS +INTEGER, ALLOCATABLE, DIMENSION(:,:) :: BODTRI_ACC +REAL(EB), ALLOCATABLE, DIMENSION(:) :: AREA_ACC +TYPE(MESH_TYPE), POINTER :: MT + +IBOD_OUT = 0; ITRI_OUT = 0 +MT => MESHES(NM_LOC) +DO PASS=1,2 + DO KK=K_LOC-1,K_LOC+1 + DO JJ=J_LOC-1,J_LOC+1 + DO II=I_LOC-1,I_LOC+1 + ICC_LOC = MT%CCVAR(II,JJ,KK,CC_IDCC); IF (ICC_LOC<1) CYCLE + DO JCC_LOC=1,MT%CUT_CELL(ICC_LOC)%NCELL + IF (PASS==1 .AND. MT%CUT_CELL(ICC_LOC)%NOADVANCE(JCC_LOC)<=0) CYCLE + IBOD_OUT = MT%CUT_CELL(ICC_LOC)%BODTRI_DONOR(1,JCC_LOC) + ITRI_OUT = MT%CUT_CELL(ICC_LOC)%BODTRI_DONOR(2,JCC_LOC) + IF (IBOD_OUT>0 .AND. ITRI_OUT>0) RETURN + ENDDO + ENDDO + ENDDO + ENDDO +ENDDO +COUNT_MAX = 0 +DO KK=K_LOC-1,K_LOC+1 + DO JJ=J_LOC-1,J_LOC+1 + DO II=I_LOC-1,I_LOC+1 + ICC_LOC = MT%CCVAR(II,JJ,KK,CC_IDCC); IF (ICC_LOC<1) CYCLE + DO JCC_LOC=1,MT%CUT_CELL(ICC_LOC)%NCELL + COUNT_MAX = COUNT_MAX + MT%CUT_CELL(ICC_LOC)%CCELEM(1,JCC_LOC) + ENDDO + ENDDO + ENDDO +ENDDO +IF (COUNT_MAX<1) RETURN +ALLOCATE(BODTRI_ACC(1:2,COUNT_MAX+1),AREA_ACC(COUNT_MAX+1)) +BODTRI_ACC = 0; AREA_ACC = 0._EB; COUNT = 0; +DO KK=K_LOC-1,K_LOC+1 + DO JJ=J_LOC-1,J_LOC+1 + DO II=I_LOC-1,I_LOC+1 + ICC_LOC = MT%CCVAR(II,JJ,KK,CC_IDCC); IF (ICC_LOC<1) CYCLE + DO JCC_LOC=1,MT%CUT_CELL(ICC_LOC)%NCELL + DO IFC=1,MT%CUT_CELL(ICC_LOC)%CCELEM(1,JCC_LOC) + IFACE = MT%CUT_CELL(ICC_LOC)%CCELEM(IFC+1,JCC_LOC) + IF (MT%CUT_CELL(ICC_LOC)%FACE_LIST(1,IFACE)/=CC_FTYPE_CFINB) CYCLE + IFC1 = MT%CUT_CELL(ICC_LOC)%FACE_LIST(4,IFACE) + JFC1 = MT%CUT_CELL(ICC_LOC)%FACE_LIST(5,IFACE) + IF (IFC1<1 .OR. IFC1>MT%N_CUTFACE_MESH+MT%N_GCCUTFACE_MESH) CYCLE + IF (.NOT.ALLOCATED(MT%CUT_FACE(IFC1)%BODTRI)) CYCLE + IF (JFC1<1 .OR. JFC1>SIZE(MT%CUT_FACE(IFC1)%BODTRI,DIM=2)) CYCLE + CALL ACCUMULATE_BLOCKING_BODTRI(MT%CUT_FACE(IFC1)%BODTRI(1:2,JFC1),MT%CUT_FACE(IFC1)%AREA(JFC1), & + COUNT,BODTRI_ACC,AREA_ACC) + ENDDO + ENDDO + ENDDO + ENDDO +ENDDO +IF (COUNT>0) THEN + DUM = MAXLOC(AREA_ACC(1:COUNT),DIM=1) + IBOD_OUT = BODTRI_ACC(1,DUM) + ITRI_OUT = BODTRI_ACC(2,DUM) +ENDIF +DEALLOCATE(BODTRI_ACC,AREA_ACC) +END SUBROUTINE GET_REFINEMENT_CELL_DONOR + +SUBROUTINE SET_REFINEMENT_CUTCELL_DONOR(NM_LOC,ICC_LOC,JCC_LOC,IBOD_IN,ITRI_IN) +INTEGER, INTENT(IN) :: NM_LOC,ICC_LOC,JCC_LOC,IBOD_IN,ITRI_IN +INTEGER :: IBOD_LOC,ITRI_LOC + +IF (IBOD_IN>0 .AND. ITRI_IN>0) THEN + MESHES(NM_LOC)%CUT_CELL(ICC_LOC)%BODTRI_DONOR(1:2,JCC_LOC) = (/ IBOD_IN, ITRI_IN /) +ELSE + CALL GET_BLOCKING_CUTCELL_DONOR(NM_LOC,ICC_LOC,JCC_LOC,IBOD_LOC,ITRI_LOC) +ENDIF +END SUBROUTINE SET_REFINEMENT_CUTCELL_DONOR + +SUBROUTINE TAG_CELL_BLOCKED_BY_REFINEMENT_FOOTPRINT(NM_LOC,I_LOC,J_LOC,K_LOC,BLOCK_TAG,IBOD_DONOR,ITRI_DONOR,CELL_CHANGED) + +INTEGER, INTENT(IN) :: NM_LOC,I_LOC,J_LOC,K_LOC,BLOCK_TAG,IBOD_DONOR,ITRI_DONOR +LOGICAL, INTENT(OUT) :: CELL_CHANGED +INTEGER :: ICC_LOC,JCC_LOC,AX_LOC,SIDE_LOC,ICFC_LOC,CT_LOC,NCFACE_CUTCELL_LOC,NFACE_CELL_LOC,NCELL_LOC +INTEGER, ALLOCATABLE, DIMENSION(:,:) :: CCELEM_LOC,FACE_LIST_LOC +INTEGER, ALLOCATABLE, DIMENSION(:) :: NOADVANCE_LOC +REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: XYZCEN_LOC +REAL(EB), ALLOCATABLE, DIMENSION(:) :: VOLUME_LOC +TYPE(MESH_TYPE), POINTER :: MT + +CELL_CHANGED = .FALSE. +MT => MESHES(NM_LOC) + +SELECT CASE(MT%CCVAR(I_LOC,J_LOC,K_LOC,CC_CGSC)) +CASE(CC_SOLID) + RETURN +CASE(CC_CUTCFE) + ICC_LOC = MT%CCVAR(I_LOC,J_LOC,K_LOC,CC_IDCC) + IF (ICC_LOC<1) RETURN + DO JCC_LOC=1,MT%CUT_CELL(ICC_LOC)%NCELL + IF (MT%CUT_CELL(ICC_LOC)%NOADVANCE(JCC_LOC)/=NOT_BLOCKED) CYCLE + MT%CUT_CELL(ICC_LOC)%NOADVANCE(JCC_LOC) = BLOCK_TAG + CALL SET_REFINEMENT_CUTCELL_DONOR(NM_LOC,ICC_LOC,JCC_LOC,IBOD_DONOR,ITRI_DONOR) + CELL_CHANGED = .TRUE. + ENDDO +CASE(CC_GASPHASE) + CT_LOC = 6 + NCFACE_CUTCELL_LOC = CT_LOC + 1 + NCELL_LOC = 1; NFACE_CELL_LOC = CT_LOC + ALLOCATE(CCELEM_LOC(1:NCFACE_CUTCELL_LOC,1:NCELL_LOC)); CCELEM_LOC = CC_UNDEFINED + ALLOCATE(FACE_LIST_LOC(1:CC_NPARAM_CCFACE,1:NFACE_CELL_LOC)); FACE_LIST_LOC = CC_UNDEFINED + ALLOCATE(VOLUME_LOC(1:NCELL_LOC)); VOLUME_LOC(1) = MT%DX(I_LOC)*MT%DY(J_LOC)*MT%DZ(K_LOC) + ALLOCATE(XYZCEN_LOC(IAXIS:KAXIS,1:NCELL_LOC)) + XYZCEN_LOC(IAXIS:KAXIS,1) = (/ MT%XC(I_LOC), MT%YC(J_LOC), MT%ZC(K_LOC) /) + ALLOCATE(NOADVANCE_LOC(1:NCELL_LOC)); NOADVANCE_LOC(1) = BLOCK_TAG + CT_LOC = 1; CCELEM_LOC(1,1) = 0 + DO AX_LOC=IAXIS,KAXIS + DO SIDE_LOC=LOW_IND,HIGH_IND + IF (.NOT.FACE_INDEX_IN_BOUNDS(NM_LOC,I_LOC+ADDI(SIDE_LOC,AX_LOC),J_LOC+ADDJ(SIDE_LOC,AX_LOC), & + K_LOC+ADDK(SIDE_LOC,AX_LOC),AX_LOC)) CYCLE + ICFC_LOC = MT%FCVAR(I_LOC+ADDI(SIDE_LOC,AX_LOC),J_LOC+ADDJ(SIDE_LOC,AX_LOC), & + K_LOC+ADDK(SIDE_LOC,AX_LOC),CC_IDCF,AX_LOC) + IF (ICFC_LOC>0) THEN + FACE_LIST_LOC(1:CC_NPARAM_CCFACE,CT_LOC) = (/ CC_FTYPE_CFGAS, SIDE_LOC, AX_LOC,ICFC_LOC, 1, CC_UNDEFINED /) + CCELEM_LOC(1,1) = CCELEM_LOC(1,1) + 1 + CCELEM_LOC(CCELEM_LOC(1,1)+1,1) = CT_LOC + CT_LOC = CT_LOC + 1 + ELSEIF (MT%FCVAR(I_LOC+ADDI(SIDE_LOC,AX_LOC),J_LOC+ADDJ(SIDE_LOC,AX_LOC), & + K_LOC+ADDK(SIDE_LOC,AX_LOC),CC_FGSC,AX_LOC)==CC_GASPHASE) THEN + FACE_LIST_LOC(1:CC_NPARAM_CCFACE,CT_LOC) = (/ CC_FTYPE_RCGAS, SIDE_LOC, AX_LOC, 0, 0, CC_UNDEFINED /) + CCELEM_LOC(1,1) = CCELEM_LOC(1,1) + 1 + CCELEM_LOC(CCELEM_LOC(1,1)+1,1) = CT_LOC + CT_LOC = CT_LOC + 1 + ENDIF + ENDDO + ENDDO + CALL INSERT_CUT_CELL(NM_LOC,I_LOC,J_LOC,K_LOC,ICC_LOC) + MT => MESHES(NM_LOC) + CALL NEW_CELL_ALLOC(NM_LOC,ICC_LOC,NCELL_LOC,NFACE_CELL_LOC,NCFACE_CUTCELL_LOC) + MT%CUT_CELL(ICC_LOC)%NCELL = NCELL_LOC + MT%CUT_CELL(ICC_LOC)%NFACE_CELL = NFACE_CELL_LOC + CALL MOVE_ALLOC(FROM=CCELEM_LOC,TO=MT%CUT_CELL(ICC_LOC)%CCELEM) + CALL MOVE_ALLOC(FROM=FACE_LIST_LOC,TO=MT%CUT_CELL(ICC_LOC)%FACE_LIST) + CALL MOVE_ALLOC(FROM=VOLUME_LOC,TO=MT%CUT_CELL(ICC_LOC)%VOLUME) + CALL MOVE_ALLOC(FROM=XYZCEN_LOC,TO=MT%CUT_CELL(ICC_LOC)%XYZCEN) + CALL MOVE_ALLOC(FROM=NOADVANCE_LOC,TO=MT%CUT_CELL(ICC_LOC)%NOADVANCE) + CALL SET_REFINEMENT_CUTCELL_DONOR(NM_LOC,ICC_LOC,1,IBOD_DONOR,ITRI_DONOR) + CELL_CHANGED = .TRUE. +END SELECT + +END SUBROUTINE TAG_CELL_BLOCKED_BY_REFINEMENT_FOOTPRINT + + +LOGICAL FUNCTION FACE_INDEX_IN_BOUNDS(NM_LOC,I_LOC,J_LOC,K_LOC,AX_LOC) + +INTEGER, INTENT(IN) :: NM_LOC,I_LOC,J_LOC,K_LOC,AX_LOC + +FACE_INDEX_IN_BOUNDS = NM_LOC>=1 .AND. NM_LOC<=NMESHES; IF (.NOT.FACE_INDEX_IN_BOUNDS) RETURN +FACE_INDEX_IN_BOUNDS = I_LOC>=LBOUND(MESHES(NM_LOC)%FCVAR,DIM=1) .AND. I_LOC<=UBOUND(MESHES(NM_LOC)%FCVAR,DIM=1) .AND. & + J_LOC>=LBOUND(MESHES(NM_LOC)%FCVAR,DIM=2) .AND. J_LOC<=UBOUND(MESHES(NM_LOC)%FCVAR,DIM=2) .AND. & + K_LOC>=LBOUND(MESHES(NM_LOC)%FCVAR,DIM=3) .AND. K_LOC<=UBOUND(MESHES(NM_LOC)%FCVAR,DIM=3) .AND. & + AX_LOC>=LBOUND(MESHES(NM_LOC)%FCVAR,DIM=5) .AND. AX_LOC<=UBOUND(MESHES(NM_LOC)%FCVAR,DIM=5) +END FUNCTION FACE_INDEX_IN_BOUNDS + SUBROUTINE ADD_NEIGHBOR_BLOCKED_CELLS USE TRAN, ONLY: GET_IJK INTEGER :: NM2,ICELL,I2,J2,K2,BLOCK_TAG -LOGICAL :: IND_FOUND REAL(EB):: XCO,YCO,ZCO,VOL_NM,VOL_NOM,X1,Y1,Z1 TYPE(MESH_TYPE), POINTER :: M2 @@ -2252,91 +2561,67 @@ SUBROUTINE ADD_NEIGHBOR_BLOCKED_CELLS M2 => MESHES(NOM) ICELL_DO : DO ICELL=1,M2%N_CC_BLOCKED - XCO = M2%XYZ_CC_BLOCKED(IAXIS,ICELL) - YCO = M2%XYZ_CC_BLOCKED(JAXIS,ICELL) - ZCO = M2%XYZ_CC_BLOCKED(KAXIS,ICELL) + XCO = M2%XYZ_CC_BLOCKED(IAXIS,ICELL); YCO = M2%XYZ_CC_BLOCKED(JAXIS,ICELL); ZCO = M2%XYZ_CC_BLOCKED(KAXIS,ICELL) BLOCK_TAG = M2%JBT_CC_BLOCKED(2,ICELL) CALL GET_IJK(XCO,YCO,ZCO,NOM,X1,Y1,Z1,I2,J2,K2) + IF (I2<1 .OR. I2>M2%IBAR .OR. J2<1 .OR. J2>M2%JBAR .OR. K2<1 .OR. K2>M2%KBAR) CYCLE ICELL_DO VOL_NOM = M2%DX(I2)*M2%DY(J2)*M2%DZ(K2) - IF (VOL_NM > 1.5_EB * VOL_NOM) THEN ! NM is COARSE, NOM is FINE - IF(XCO < M2%XS .OR. XCO > M2%XF .OR. & - YCO < M2%YS .OR. YCO > M2%YF .OR. & - ZCO < M2%ZS .OR. ZCO > M2%ZF) CYCLE ICELL_DO + IF(XCO < M2%XS .OR. XCO > M2%XF .OR. YCO < M2%YS .OR. YCO > M2%YF .OR. ZCO < M2%ZS .OR. ZCO > M2%ZF) CYCLE ICELL_DO IF(XCO > M2%X(1) .AND. XCO < M2%X(M2%IBAR-1) .AND. & YCO > M2%Y(1) .AND. YCO < M2%Y(M2%JBAR-1) .AND. & ZCO > M2%Z(1) .AND. ZCO < M2%Z(M2%KBAR-1)) CYCLE ICELL_DO - ! Find I,J,K in NM where (XCO,YCO,ZCO) falls within cell bounds - IND_FOUND = .FALSE. - DO I=ILO_CELL-1,IHI_CELL+1 - IF (XCO < XFACE(I-1)-GEOMEPS .OR. XCO > XFACE(I)+GEOMEPS) CYCLE - DO J=JLO_CELL-1,JHI_CELL+1 - IF (YCO < YFACE(J-1)-GEOMEPS .OR. YCO > YFACE(J)+GEOMEPS) CYCLE - DO K=KLO_CELL-1,KHI_CELL+1 - IF (ZCO < ZFACE(K-1)-GEOMEPS .OR. ZCO > ZFACE(K)+GEOMEPS) CYCLE - IF (I > ILO_CELL-1 .AND. I < IHI_CELL+1 .AND. & - J > JLO_CELL-1 .AND. J < JHI_CELL+1 .AND. & - K > KLO_CELL-1 .AND. K < KHI_CELL+1) CYCLE - IND_FOUND = .TRUE. - EXIT - ENDDO - IF (IND_FOUND) EXIT - ENDDO - IF (IND_FOUND) EXIT - ENDDO - IF (.NOT.IND_FOUND) CYCLE ICELL_DO + ! Find I,J,K in NM where (XCO,YCO,ZCO) falls within cell bounds (ghost band only). + I = MINLOC(ABS(XCELL(ILO_CELL-1:IHI_CELL+1)-XCO),DIM=1) + ILO_CELL - 2 + J = MINLOC(ABS(YCELL(JLO_CELL-1:JHI_CELL+1)-YCO),DIM=1) + JLO_CELL - 2 + K = MINLOC(ABS(ZCELL(KLO_CELL-1:KHI_CELL+1)-ZCO),DIM=1) + KLO_CELL - 2 + IF (XCO < XFACE(I-1)-GEOMEPS .OR. XCO > XFACE(I)+GEOMEPS .OR. & + YCO < YFACE(J-1)-GEOMEPS .OR. YCO > YFACE(J)+GEOMEPS .OR. & + ZCO < ZFACE(K-1)-GEOMEPS .OR. ZCO > ZFACE(K)+GEOMEPS) CYCLE ICELL_DO + IF (I>ILO_CELL-1 .AND. IJLO_CELL-1 .AND. JKLO_CELL-1 .AND. K 0) THEN DO JCC = 1, M%CUT_CELL(ICC)%NCELL - IF (M%CUT_CELL(ICC)%NOADVANCE(JCC) == NOT_BLOCKED) & + IF (M%CUT_CELL(ICC)%NOADVANCE(JCC) == NOT_BLOCKED) THEN M%CUT_CELL(ICC)%NOADVANCE(JCC) = BLOCK_TAG + IF (M2%JBT_CC_BLOCKED(3,ICELL)>0 .AND. M2%JBT_CC_BLOCKED(4,ICELL)>0) & + M%CUT_CELL(ICC)%BODTRI_DONOR(1:2,JCC) = M2%JBT_CC_BLOCKED(3:4,ICELL) + ENDIF ENDDO ENDIF - + ELSEIF (VOL_NOM > 1.5_EB * VOL_NM) THEN ! NM is FINE, NOM is COARSE + CALL TAG_FINE_CELLS_IN_COARSE_CELL_VOLUME(NOM,I2,J2,K2,NM,BLOCK_TAG, & + M2%JBT_CC_BLOCKED(3,ICELL),M2%JBT_CC_BLOCKED(4,ICELL)) ELSE - ! ===================================================== ! Same refinement level (or refinement handled by EXT_WALL_LOOP) - use centroid matching - ! ===================================================== - IND_FOUND = .FALSE. - DO I=ILO_CELL-1,IHI_CELL+1 - IF (ABS(XCO-XCELL(I))=GEOMEPS .OR. ABS(YCO-YCELL(J))>=GEOMEPS .OR. & + ABS(ZCO-ZCELL(K))>=GEOMEPS) CYCLE ICELL_DO ! Here we have found the I,J,K indices of the blocked cut-cell: ICC=M%CCVAR(I,J,K,CC_IDCC) - IF(ICC>0) M%CUT_CELL(ICC)%NOADVANCE(M2%JBT_CC_BLOCKED(1,ICELL)) = BLOCK_TAG - + IF (ICC>0) THEN + IF (M2%JBT_CC_BLOCKED(1,ICELL)>0) THEN + IF (M2%JBT_CC_BLOCKED(1,ICELL)<=M%CUT_CELL(ICC)%NCELL) THEN + M%CUT_CELL(ICC)%NOADVANCE(M2%JBT_CC_BLOCKED(1,ICELL)) = BLOCK_TAG + IF (M2%JBT_CC_BLOCKED(3,ICELL)>0 .AND. M2%JBT_CC_BLOCKED(4,ICELL)>0) & + M%CUT_CELL(ICC)%BODTRI_DONOR(1:2,M2%JBT_CC_BLOCKED(1,ICELL)) = M2%JBT_CC_BLOCKED(3:4,ICELL) + ENDIF + ENDIF + ENDIF ENDIF ENDDO ICELL_DO ENDDO NEIGHBORING_MESHES_DO CALL DEFINE_XYZFACE_CELL(ALLOC_FLG=.FALSE.) ENDDO MESH_LOOP - END SUBROUTINE ADD_NEIGHBOR_BLOCKED_CELLS @@ -2470,7 +2755,6 @@ END SUBROUTINE DEFINE_XYZFACE_CELL SUBROUTINE TAG_CC_BLOCKING_REFINEMENT - LOGICAL, PARAMETER :: DO_RAY_TRACING=.TRUE. INTEGER :: DUM,II1,JJ1,KK1,IIO1,JJO1,KKO1,IIO2,JJO2,KKO2,IIG,JJG,KKG,IIOG,JJOG,KKOG @@ -2644,6 +2928,7 @@ SUBROUTINE TAG_CC_BLOCKING_REFINEMENT END SELECT IF(M%FCVAR(IIF,JJF,KKF,CC_CGSC,X1AXIS)==CC_SOLID) CYCLE EXT_WALL_LOOP ! No need to do anything. IF(M2%FCVAR(IIOF,JJOF,KKOF,CC_CGSC,X1AXIS)==CC_SOLID) THEN ! Coarse side face is solid. + CALL GET_REFINEMENT_CELL_DONOR(NOM,IIOF,JJOF,KKOF,IBOD_DONOR,ITRI_DONOR) ! Set II,JJ,KK fine cells next to this EWC for blocking. IF(M%CCVAR(II,JJ,KK,CC_CGSC)==CC_GASPHASE) THEN ! Insert cut-cell in this location, set to Block. @@ -2681,6 +2966,7 @@ SUBROUTINE TAG_CC_BLOCKING_REFINEMENT CALL MOVE_ALLOC(FROM=VOLUME ,TO=M%CUT_CELL(ICC)%VOLUME) CALL MOVE_ALLOC(FROM=XYZCEN ,TO=M%CUT_CELL(ICC)%XYZCEN) CALL MOVE_ALLOC(FROM=NOADVANCE,TO=M%CUT_CELL(ICC)%NOADVANCE) + CALL SET_REFINEMENT_CUTCELL_DONOR(NM,ICC,1,IBOD_DONOR,ITRI_DONOR) ELSEIF(M%CCVAR(II,JJ,KK,CC_CGSC)==CC_CUTCFE) THEN ! Find face in IIF,JFF,JJF,X1AXIS location and set the owner cell to block. ICC=M%CCVAR(II,JJ,KK,CC_IDCC); CC=> M%CUT_CELL(ICC) @@ -2690,6 +2976,7 @@ SUBROUTINE TAG_CC_BLOCKING_REFINEMENT IF( (CC%FACE_LIST(1,IFACE)==CC_FTYPE_RCGAS .OR. CC%FACE_LIST(1,IFACE)==CC_FTYPE_CFGAS) .AND. & CC%FACE_LIST(2,IFACE)==LOHIF .AND. CC%FACE_LIST(3,IFACE)==X1AXIS ) THEN IF(CC%NOADVANCE(JCC)==NOT_BLOCKED) CC%NOADVANCE(JCC)=BLOCKED_REFI_INTER + CALL SET_REFINEMENT_CUTCELL_DONOR(NM,ICC,JCC,IBOD_DONOR,ITRI_DONOR) CYCLE JCC_LOOP_1 ENDIF ENDDO @@ -2701,6 +2988,8 @@ SUBROUTINE TAG_CC_BLOCKING_REFINEMENT ELSEIF((EWC%IIO_MAX-EWC%IIO_MIN+1)*(EWC%JJO_MAX-EWC%JJO_MIN+1)*(EWC%KKO_MAX-EWC%KKO_MIN+1)>1) THEN IF(M%FCVAR(IIF,JJF,KKF,CC_CGSC,X1AXIS)==CC_SOLID) THEN ! Coarse side face is solid. + IF (PROCESS(NOM)/=MY_RANK) CYCLE EXT_WALL_LOOP + CALL GET_REFINEMENT_CELL_DONOR(NM,IIF,JJF,KKF,IBOD_DONOR,ITRI_DONOR) ! If needed, block ghost cells of the other mesh which has finer cells. DO KKO=EWC%KKO_MIN,EWC%KKO_MAX DO JJO=EWC%JJO_MIN,EWC%JJO_MAX @@ -2754,6 +3043,7 @@ SUBROUTINE TAG_CC_BLOCKING_REFINEMENT CALL MOVE_ALLOC(FROM=VOLUME ,TO=M2%CUT_CELL(ICC)%VOLUME) CALL MOVE_ALLOC(FROM=XYZCEN ,TO=M2%CUT_CELL(ICC)%XYZCEN) CALL MOVE_ALLOC(FROM=NOADVANCE,TO=M2%CUT_CELL(ICC)%NOADVANCE) + CALL SET_REFINEMENT_CUTCELL_DONOR(NOM,ICC,1,IBOD_DONOR,ITRI_DONOR) ELSEIF(M2%CCVAR(IOGC,JOGC,KOGC,CC_CGSC)==CC_CUTCFE) THEN ! Find face in IIF,JFF,JJF,X1AXIS location and set the owner cell to block. ICC=M2%CCVAR(IOGC,JOGC,KOGC,CC_IDCC); CC=> M2%CUT_CELL(ICC) @@ -2763,6 +3053,7 @@ SUBROUTINE TAG_CC_BLOCKING_REFINEMENT IF( (CC%FACE_LIST(1,IFACE)==CC_FTYPE_RCGAS .OR. CC%FACE_LIST(1,IFACE)==CC_FTYPE_CFGAS) .AND. & CC%FACE_LIST(2,IFACE)==LOHIF .AND. CC%FACE_LIST(3,IFACE)==X1AXIS ) THEN IF(CC%NOADVANCE(JCC)==NOT_BLOCKED) CC%NOADVANCE(JCC)=BLOCKED_REFI_INTER + CALL SET_REFINEMENT_CUTCELL_DONOR(NOM,ICC,JCC,IBOD_DONOR,ITRI_DONOR) CYCLE JCC_LOOP_3 ENDIF ENDDO @@ -2786,6 +3077,7 @@ SUBROUTINE TAG_BLOCK_CELL(NM,II1,JJ1,KK1,NOM,IIO1,JJO1,KKO1,FINE_CELL) INTEGER, INTENT(IN) :: NM,II1,JJ1,KK1,NOM,IIO1,JJO1,KKO1 LOGICAL, INTENT(IN) :: FINE_CELL +INTEGER :: JCC_LOC,IBOD_LOC,ITRI_LOC TYPE(MESH_TYPE), POINTER :: M,M2 M =>MESHES( NM) M2=>MESHES(NOM) @@ -2794,6 +3086,7 @@ SUBROUTINE TAG_BLOCK_CELL(NM,II1,JJ1,KK1,NOM,IIO1,JJO1,KKO1,FINE_CELL) ICC2 = M2%CCVAR(IIO1,JJO1,KKO1,CC_IDCC); ICC = 0 IF ( ICC2 > 0 .OR. M2%CCVAR(IIO1,JJO1,KKO1,CC_CGSC)==CC_SOLID) THEN ! There are cut-cells in omesh cartesian cell. + CALL GET_REFINEMENT_CELL_DONOR(NOM,IIO1,JJO1,KKO1,IBOD_LOC,ITRI_LOC) IF(M%CCVAR(II1,JJ1,KK1,CC_CGSC)==CC_GASPHASE) THEN ! Insert cut-cell is this location: CT = 6; @@ -2831,14 +3124,19 @@ SUBROUTINE TAG_BLOCK_CELL(NM,II1,JJ1,KK1,NOM,IIO1,JJO1,KKO1,FINE_CELL) CALL MOVE_ALLOC(FROM=VOLUME ,TO=M%CUT_CELL(ICC)%VOLUME) CALL MOVE_ALLOC(FROM=XYZCEN ,TO=M%CUT_CELL(ICC)%XYZCEN) CALL MOVE_ALLOC(FROM=NOADVANCE,TO=M%CUT_CELL(ICC)%NOADVANCE) + IF (M%CUT_CELL(ICC)%NOADVANCE(1)>0) & + CALL SET_REFINEMENT_CUTCELL_DONOR(NM,ICC,1,IBOD_LOC,ITRI_LOC) ELSEIF(M%CCVAR(II1,JJ1,KK1,CC_IDCC)>0) THEN ICC = M%CCVAR(II1,JJ1,KK1,CC_IDCC) ENDIF ! Here Test if cut-cells in II,KK,KK are blocked or not in IIO,JJO,KKO: IF(ICC>0) THEN IF(M2%CCVAR(IIO1,JJO1,KKO1,CC_CGSC)==CC_SOLID) THEN - WHERE(M%CUT_CELL(ICC)%NOADVANCE(1:M%CUT_CELL(ICC)%NCELL)==NOT_BLOCKED) & - M%CUT_CELL(ICC)%NOADVANCE(1:M%CUT_CELL(ICC)%NCELL) = BLOCKED_REFI_INTER + DO JCC_LOC=1,M%CUT_CELL(ICC)%NCELL + IF(M%CUT_CELL(ICC)%NOADVANCE(JCC_LOC)/=NOT_BLOCKED) CYCLE + M%CUT_CELL(ICC)%NOADVANCE(JCC_LOC) = BLOCKED_REFI_INTER + CALL SET_REFINEMENT_CUTCELL_DONOR(NM,ICC,JCC_LOC,IBOD_LOC,ITRI_LOC) + ENDDO ELSE; CALL TEST_CC_FOR_BLOCKING(NM,ICC,NOM,IIO1,JJO1,KKO1,ICC2) ENDIF ENDIF @@ -2846,6 +3144,7 @@ SUBROUTINE TAG_BLOCK_CELL(NM,II1,JJ1,KK1,NOM,IIO1,JJO1,KKO1,FINE_CELL) ELSE + IF (PROCESS(NOM)/=MY_RANK) RETURN ICC = M%CCVAR(II1,JJ1,KK1,CC_IDCC); ICC2 = 0 IF(ICC>0) THEN ! Set IOGC,JOGC,KOGC fine cells next to this EWC for blocking. @@ -2905,7 +3204,7 @@ SUBROUTINE TEST_CC_FOR_BLOCKING(NM,ICC,NOM,IIO1,JJO1,KKO1,ICC2) INTEGER, INTENT(IN) :: NM,ICC,NOM,IIO1,JJO1,KKO1,ICC2 INTEGER :: JCC,FC_FOUND,FC_TYPE,INBFC,INBFC_LOC,VERT_CUTFACE,NVERT,X1AXIS,X2AXIS,X3AXIS,NCROSS,DIRRAY,IFC1,JFC1,& - NVERT2,VERT_CUTFACE2,IV,IFCC,IFACE2,IFC2,JFC2 + NVERT2,VERT_CUTFACE2,IV,IFCC,IFACE2,IFC2,JFC2,IBOD_DONOR,ITRI_DONOR TYPE(MESH_TYPE), POINTER :: M,M2 TYPE(CC_CUTCELL_TYPE), POINTER :: CC,CC2 TYPE(CC_CUTFACE_TYPE), POINTER :: CF2 @@ -2922,6 +3221,7 @@ SUBROUTINE TEST_CC_FOR_BLOCKING(NM,ICC,NOM,IIO1,JJO1,KKO1,ICC2) M =>MESHES( NM) M2=>MESHES(NOM) +CALL GET_REFINEMENT_CELL_DONOR(NOM,IIO1,JJO1,KKO1,IBOD_DONOR,ITRI_DONOR) INBFC=M2%CCVAR(IIO1,JJO1,KKO1,CC_IDCF); IF(INBFC<1) RETURN ! No CC_INBOUNDARY faces in this cartesian cell. @@ -3164,7 +3464,10 @@ SUBROUTINE TEST_CC_FOR_BLOCKING(NM,ICC,NOM,IIO1,JJO1,KKO1,ICC2) ! ENDIF ENDDO INBFC_LOC_LOOP ! Here set no ADVANCE if BLOCK_CELL=T: - IF(BLOCK_CELL .AND. CC%NOADVANCE(JCC)==NOT_BLOCKED) CC%NOADVANCE(JCC) = BLOCKED_REFI_INTER + IF(BLOCK_CELL .AND. CC%NOADVANCE(JCC)==NOT_BLOCKED) THEN + CC%NOADVANCE(JCC) = BLOCKED_REFI_INTER + CALL SET_REFINEMENT_CUTCELL_DONOR(NM,ICC,JCC,IBOD_DONOR,ITRI_DONOR) + ENDIF ENDDO JCC_LOOP ! IF(NM==1 .AND. ICC<30) CLOSE(LU_CCELL) @@ -4719,7 +5022,7 @@ SUBROUTINE SET_GC_CUTCELLS_3D DO JJO=EWC%JJO_MIN,EWC%JJO_MAX DO IIO=EWC%IIO_MIN,EWC%IIO_MAX ICC = MESHES(NOM)%CCVAR(IIO,JJO,KKO,CC_IDCC) - IF (ICC > 0) THEN + IF (ICC > 0 .AND. MESHES(NOM)%CCVAR(IIO,JJO,KKO,CC_CGSC)==CC_CUTCFE) THEN N_CF = N_CF + 1 MESHES(NM)%CUT_CELL(NMICC)%NOMICC(1:2,N_CF) = (/ NOM, ICC /) NCELL = MESHES(NOM)%CUT_CELL(ICC)%NCELL @@ -4758,7 +5061,8 @@ SUBROUTINE SET_GC_CUTCELLS_3D IF(FCVAR(IIO ,JJO ,KKO ,CC_FGSC,KAXIS) == CC_SOLID) CYCLE END SELECT IF (MESHES(NM)%CCVAR(II,JJ,KK,CC_CGSC) == CC_CUTCFE) THEN - ICC = MESHES(NOM)%CCVAR(IIO,JJO,KKO,CC_IDCC); IF (ICC<1) CYCLE + ICC = MESHES(NOM)%CCVAR(IIO,JJO,KKO,CC_IDCC) + IF (ICC<1 .OR. MESHES(NOM)%CCVAR(IIO,JJO,KKO,CC_CGSC)/=CC_CUTCFE) CYCLE NMICC = MESHES(NM)%CCVAR(II,JJ,KK,CC_IDCC) NOFC = 1 ALLOCATE(MESHES(NM)%CUT_CELL(NMICC)%NOMICC(2,NOFC)); MESHES(NM)%CUT_CELL(NMICC)%NOMICC(1:2,1:NOFC) = 0 @@ -4975,7 +5279,7 @@ SUBROUTINE SNAP_GEOM_NODES ENDDO AXIS_LOOP_1 ENDDO MAIN_GEOM_LOOP_1 -! Now Mesh loop on mesh + guard planes to test against +! Now mesh loop on mesh + guard planes to test against. ! Main Loop over Meshes: MAIN_MESH_LOOP : DO NM=1,NMESHES @@ -5036,6 +5340,7 @@ SUBROUTINE SNAP_GEOM_NODES CALL DEFINE_XYZFACE_CELL(ALLOC_FLG=.FALSE.) ENDDO MAIN_MESH_LOOP + ! Deallocate SNAP_NODE in geometries: DO IG=1,N_GEOMETRY DEALLOCATE(GEOMETRY(IG)%SNAP_NODE) @@ -5051,6 +5356,127 @@ END SUBROUTINE SNAP_GEOM_NODES END SUBROUTINE SET_CUTCELLS_3D +SUBROUTINE ACCUMULATE_BLOCKING_BODTRI(BODTRI_FACE,FACE_AREA,COUNT,BODTRI_ACC,AREA_ACC) + +INTEGER, INTENT(IN) :: BODTRI_FACE(1:2) +REAL(EB), INTENT(IN) :: FACE_AREA +INTEGER, INTENT(INOUT) :: COUNT +INTEGER, INTENT(INOUT), DIMENSION(:,:) :: BODTRI_ACC +REAL(EB), INTENT(INOUT), DIMENSION(:) :: AREA_ACC +INTEGER :: DUM + +IF (SIZE(BODTRI_ACC,DIM=1)<2 .OR. SIZE(AREA_ACC)<1) RETURN +IF (BODTRI_FACE(1)<1 .OR. BODTRI_FACE(2)<1) RETURN +DUM = 1 +DO DUM=1,COUNT + IF (ALL(BODTRI_ACC(1:2,DUM)==BODTRI_FACE(1:2))) EXIT +ENDDO +IF (DUM>COUNT) THEN + BODTRI_ACC(1:2,DUM) = BODTRI_FACE(1:2) + COUNT = DUM +ENDIF +AREA_ACC(DUM) = AREA_ACC(DUM) + MAX(FACE_AREA,0._EB) + +END SUBROUTINE ACCUMULATE_BLOCKING_BODTRI + + +SUBROUTINE GET_BLOCKING_CUTCELL_DONOR(NM_LOC,ICC_LOC,JCC_LOC,IBOD_OUT,ITRI_OUT) + +INTEGER, INTENT(IN) :: NM_LOC,ICC_LOC,JCC_LOC +INTEGER, INTENT(OUT) :: IBOD_OUT,ITRI_OUT +INTEGER :: I_LOC,J_LOC,K_LOC,II,JJ,KK,ICC2,JCC2,IFC,IFACE,IFC1,JFC1,COUNT,COUNT_MAX,DUM +INTEGER, ALLOCATABLE, DIMENSION(:,:) :: BODTRI_ACC +REAL(EB), ALLOCATABLE, DIMENSION(:) :: AREA_ACC +TYPE(MESH_TYPE), POINTER :: MT + +IBOD_OUT = 0 +ITRI_OUT = 0 +MT => MESHES(NM_LOC) +IBOD_OUT = MT%CUT_CELL(ICC_LOC)%BODTRI_DONOR(1,JCC_LOC) +ITRI_OUT = MT%CUT_CELL(ICC_LOC)%BODTRI_DONOR(2,JCC_LOC) +IF (IBOD_OUT>0 .AND. ITRI_OUT>0) RETURN + +COUNT_MAX = MAX(1,MT%CUT_CELL(ICC_LOC)%CCELEM(1,JCC_LOC)) +ALLOCATE(BODTRI_ACC(1:2,COUNT_MAX+1),AREA_ACC(COUNT_MAX+1)) +BODTRI_ACC = 0; AREA_ACC = 0._EB;COUNT = 0 +DO IFC=1,MT%CUT_CELL(ICC_LOC)%CCELEM(1,JCC_LOC) + IFACE = MT%CUT_CELL(ICC_LOC)%CCELEM(IFC+1,JCC_LOC) + IF (MT%CUT_CELL(ICC_LOC)%FACE_LIST(1,IFACE)/=CC_FTYPE_CFINB) CYCLE + IFC1 = MT%CUT_CELL(ICC_LOC)%FACE_LIST(4,IFACE) + JFC1 = MT%CUT_CELL(ICC_LOC)%FACE_LIST(5,IFACE) + IF (IFC1<1 .OR. IFC1>MT%N_CUTFACE_MESH+MT%N_GCCUTFACE_MESH) CYCLE + IF (.NOT.ALLOCATED(MT%CUT_FACE(IFC1)%BODTRI)) CYCLE + IF (JFC1<1 .OR. JFC1>SIZE(MT%CUT_FACE(IFC1)%BODTRI,DIM=2)) CYCLE + CALL ACCUMULATE_BLOCKING_BODTRI(MT%CUT_FACE(IFC1)%BODTRI(1:2,JFC1),MT%CUT_FACE(IFC1)%AREA(JFC1), & + COUNT,BODTRI_ACC,AREA_ACC) +ENDDO +IF (COUNT<1) THEN + DEALLOCATE(BODTRI_ACC,AREA_ACC) + COUNT_MAX = 0 + I_LOC = MT%CUT_CELL(ICC_LOC)%IJK(IAXIS) + J_LOC = MT%CUT_CELL(ICC_LOC)%IJK(JAXIS) + K_LOC = MT%CUT_CELL(ICC_LOC)%IJK(KAXIS) + DO KK=K_LOC-1,K_LOC+1 + DO JJ=J_LOC-1,J_LOC+1 + DO II=I_LOC-1,I_LOC+1 + ICC2 = MT%CCVAR(II,JJ,KK,CC_IDCC); IF (ICC2<1) CYCLE + DO JCC2=1,MT%CUT_CELL(ICC2)%NCELL + COUNT_MAX = COUNT_MAX + MT%CUT_CELL(ICC2)%CCELEM(1,JCC2) + ENDDO + ENDDO + ENDDO + ENDDO + IF (COUNT_MAX<1) RETURN + ALLOCATE(BODTRI_ACC(1:2,COUNT_MAX+1),AREA_ACC(COUNT_MAX+1)) + BODTRI_ACC = 0 + AREA_ACC = 0._EB + COUNT = 0 + DO KK=K_LOC-1,K_LOC+1 + DO JJ=J_LOC-1,J_LOC+1 + DO II=I_LOC-1,I_LOC+1 + ICC2 = MT%CCVAR(II,JJ,KK,CC_IDCC); IF (ICC2<1) CYCLE + DO JCC2=1,MT%CUT_CELL(ICC2)%NCELL + DO IFC=1,MT%CUT_CELL(ICC2)%CCELEM(1,JCC2) + IFACE = MT%CUT_CELL(ICC2)%CCELEM(IFC+1,JCC2) + IF (MT%CUT_CELL(ICC2)%FACE_LIST(1,IFACE)/=CC_FTYPE_CFINB) CYCLE + IFC1 = MT%CUT_CELL(ICC2)%FACE_LIST(4,IFACE) + JFC1 = MT%CUT_CELL(ICC2)%FACE_LIST(5,IFACE) + IF (IFC1<1 .OR. IFC1>MT%N_CUTFACE_MESH+MT%N_GCCUTFACE_MESH) CYCLE + IF (.NOT.ALLOCATED(MT%CUT_FACE(IFC1)%BODTRI)) CYCLE + IF (JFC1<1 .OR. JFC1>SIZE(MT%CUT_FACE(IFC1)%BODTRI,DIM=2)) CYCLE + CALL ACCUMULATE_BLOCKING_BODTRI(MT%CUT_FACE(IFC1)%BODTRI(1:2,JFC1),MT%CUT_FACE(IFC1)%AREA(JFC1), & + COUNT,BODTRI_ACC,AREA_ACC) + ENDDO + ENDDO + ENDDO + ENDDO + ENDDO +ENDIF +IF (COUNT>0) THEN + DUM = MAXLOC(AREA_ACC(1:COUNT),DIM=1) + IBOD_OUT = BODTRI_ACC(1,DUM) + ITRI_OUT = BODTRI_ACC(2,DUM) + MT%CUT_CELL(ICC_LOC)%BODTRI_DONOR(1:2,JCC_LOC) = (/ IBOD_OUT, ITRI_OUT /) +ENDIF +IF (ALLOCATED(BODTRI_ACC)) DEALLOCATE(BODTRI_ACC) +IF (ALLOCATED(AREA_ACC)) DEALLOCATE(AREA_ACC) + +END SUBROUTINE GET_BLOCKING_CUTCELL_DONOR + + +LOGICAL FUNCTION VALID_GEOMETRY_FACE_DONOR(IBOD,IWSEL) + +INTEGER, INTENT(IN) :: IBOD,IWSEL + +VALID_GEOMETRY_FACE_DONOR = .FALSE. +IF (.NOT.ALLOCATED(GEOMETRY)) RETURN +IF (IBOD<1 .OR. IBOD>SIZE(GEOMETRY,DIM=1)) RETURN +IF (.NOT.ALLOCATED(GEOMETRY(IBOD)%FACES_NORMAL)) RETURN +IF (IWSEL<1 .OR. IWSEL>SIZE(GEOMETRY(IBOD)%FACES_NORMAL,DIM=2)) RETURN +VALID_GEOMETRY_FACE_DONOR = .TRUE. + +END FUNCTION VALID_GEOMETRY_FACE_DONOR + ! ----------------------- CHECK_WALL_CELL_PLANE_MATCH ---------------------------- SUBROUTINE CHECK_WALL_CELL_PLANE_MATCH @@ -5115,7 +5541,7 @@ SUBROUTINE EXCHANGE_CC_NOADVANCE_INFO USE MPI_F08 ! Local Variables: - INTEGER :: NM,NOM,N,IERR,I,J,K,ICC,JCC + INTEGER :: NM,NOM,N,IERR,I,J,K,ICC,JCC,IBOD_DONOR,ITRI_DONOR TYPE(MESH_TYPE), POINTER :: M TYPE (MPI_REQUEST), ALLOCATABLE, DIMENSION(:) :: REQ0,REQ0DUM INTEGER :: N_REQ0 @@ -5128,28 +5554,52 @@ SUBROUTINE EXCHANGE_CC_NOADVANCE_INFO ! Count cut-cells for blocking in mesh: M%N_CC_BLOCKED = 0 DO ICC=1,MESHES(NM)%N_CUTCELL_MESH - CC => CUT_CELL(ICC) + CC => M%CUT_CELL(ICC) DO JCC=1,CC%NCELL IF(CC%NOADVANCE(JCC)>0) M%N_CC_BLOCKED = M%N_CC_BLOCKED + 1 ENDDO ENDDO + ! Also count CC_SOLID interior cells in the outer 1-cell boundary band, so they get + ! projected to fine ghost cells at refinement interfaces. + DO K=1,M%KBAR + DO J=1,M%JBAR + DO I=1,M%IBAR + IF (M%CCVAR(I,J,K,CC_CGSC)/=CC_SOLID) CYCLE + IF (I>1 .AND. I1 .AND. J1 .AND. K0) THEN IF(ALLOCATED(M%XYZ_CC_BLOCKED)) DEALLOCATE(M%XYZ_CC_BLOCKED) IF(ALLOCATED(M%JBT_CC_BLOCKED)) DEALLOCATE(M%JBT_CC_BLOCKED) ALLOCATE(M%XYZ_CC_BLOCKED(3,M%N_CC_BLOCKED)) - ALLOCATE(M%JBT_CC_BLOCKED(2,M%N_CC_BLOCKED)) + ALLOCATE(M%JBT_CC_BLOCKED(4,M%N_CC_BLOCKED)) ! Fill in blocked cut-cell info: M%N_CC_BLOCKED = 0 DO ICC=1,MESHES(NM)%N_CUTCELL_MESH - CC => CUT_CELL(ICC); I = CC%IJK(IAXIS); J = CC%IJK(JAXIS); K = CC%IJK(KAXIS) + CC => M%CUT_CELL(ICC); I = CC%IJK(IAXIS); J = CC%IJK(JAXIS); K = CC%IJK(KAXIS) DO JCC=1,CC%NCELL IF(CC%NOADVANCE(JCC)>0) THEN M%N_CC_BLOCKED = M%N_CC_BLOCKED + 1 - M%XYZ_CC_BLOCKED(1:3,M%N_CC_BLOCKED) = (/XC(I),YC(J),ZC(K)/) - M%JBT_CC_BLOCKED(1:2,M%N_CC_BLOCKED) = (/JCC,CC%NOADVANCE(JCC)/) + M%XYZ_CC_BLOCKED(1:3,M%N_CC_BLOCKED) = (/M%XC(I),M%YC(J),M%ZC(K)/) + CALL GET_BLOCKING_CUTCELL_DONOR(NM,ICC,JCC,IBOD_DONOR,ITRI_DONOR) + M%JBT_CC_BLOCKED(1:4,M%N_CC_BLOCKED) = (/JCC,CC%NOADVANCE(JCC),IBOD_DONOR,ITRI_DONOR/) ENDIF ENDDO ENDDO + ! Fill in CC_SOLID boundary-band coarse cells so they project to fine ghost cells. + DO K=1,M%KBAR + DO J=1,M%JBAR + DO I=1,M%IBAR + IF (M%CCVAR(I,J,K,CC_CGSC)/=CC_SOLID) CYCLE + IF (I>1 .AND. I1 .AND. J1 .AND. K0) THEN - ALLOCATE(BODTRI(2,COUNT+1),AREA(COUNT+1)); BODTRI=0; AREA=0._EB; COUNT = 0 - DO IFC=1,M%CUT_CELL(ICC)%CCELEM(1,JCC) - IFACE = M%CUT_CELL(ICC)%CCELEM(IFC+1,JCC) - IFC1 = M%CUT_CELL(ICC)%FACE_LIST(4,IFACE) - JFC1 = M%CUT_CELL(ICC)%FACE_LIST(5,IFACE) - IF (M%CUT_CELL(ICC)%FACE_LIST(1,IFACE) /= CC_FTYPE_CFINB) CYCLE - DO DUM=1,COUNT - IF( BODTRI(1,DUM)==M%CUT_FACE(IFC1)%BODTRI(1,JFC1) .AND. & - BODTRI(2,DUM)==M%CUT_FACE(IFC1)%BODTRI(2,JFC1) ) EXIT - ENDDO - IF(DUM > COUNT) THEN ! No match in previous loop DUM=COUNT+1 - BODTRI(1:2,DUM)=M%CUT_FACE(IFC1)%BODTRI(1:2,JFC1) - COUNT = DUM - ENDIF - AREA(DUM) = AREA(DUM) + M%CUT_FACE(IFC1)%AREA(JFC1) - ENDDO - IF (COUNT>0) THEN - ! Now set IBOD, ITRI - DUM = MAXLOC(AREA,DIM=1) ! BOD,TRI and SURF_ID with max area in cc being blocked. - IBOD= BODTRI(1,DUM); ITRI= BODTRI(2,DUM) - ENDIF - DEALLOCATE(BODTRI,AREA) -ELSE - ! Look in surrounding cells: - DO KK=K-1,K+1 - DO JJ=J-1,J+1 - DO II=I-1,I+1 - ICC2=M%CCVAR(II,JJ,KK,CC_IDCC) - IF (ICC2>0) THEN - DO JCC2=1,M%CUT_CELL(ICC2)%NCELL - DO IFC=1,M%CUT_CELL(ICC2)%CCELEM(1,JCC2) - IFACE = M%CUT_CELL(ICC2)%CCELEM(IFC+1,JCC2) - IF (M%CUT_CELL(ICC2)%FACE_LIST(1,IFACE) == CC_FTYPE_CFINB) COUNT = COUNT + 1 - ENDDO - ENDDO - ENDIF - ENDDO - ENDDO - ENDDO - IF (COUNT>0) THEN - ALLOCATE(BODTRI(2,COUNT+1),AREA(COUNT+1)); BODTRI=0; AREA=0._EB; COUNT = 0 - DO KK=K-1,K+1 - DO JJ=J-1,J+1 - DO II=I-1,I+1 - ICC2=M%CCVAR(II,JJ,KK,CC_IDCC) - IF (ICC2>0) THEN - DO JCC2=1,M%CUT_CELL(ICC2)%NCELL - DO IFC=1,M%CUT_CELL(ICC2)%CCELEM(1,JCC2) - IFACE = M%CUT_CELL(ICC2)%CCELEM(IFC+1,JCC2) - IFC1 = M%CUT_CELL(ICC2)%FACE_LIST(4,IFACE) - JFC1 = M%CUT_CELL(ICC2)%FACE_LIST(5,IFACE) - IF (M%CUT_CELL(ICC2)%FACE_LIST(1,IFACE) /= CC_FTYPE_CFINB) CYCLE - DO DUM=1,COUNT - IF( BODTRI(1,DUM)==M%CUT_FACE(IFC1)%BODTRI(1,JFC1) .AND. & - BODTRI(2,DUM)==M%CUT_FACE(IFC1)%BODTRI(2,JFC1) ) EXIT - ENDDO - IF(DUM > COUNT) THEN - BODTRI(1:2,DUM)=M%CUT_FACE(IFC1)%BODTRI(1:2,JFC1) - COUNT = DUM - ENDIF - AREA(DUM) = AREA(DUM) + M%CUT_FACE(IFC1)%AREA(JFC1) - ENDDO - ENDDO - ENDIF - ENDDO - ENDDO - ENDDO - IF (COUNT>0) THEN - ! Now set IBOD, ITRI - DUM = MAXLOC(AREA,DIM=1) ! BOD,TRI and SURF_ID with max area in cc being blocked. - IBOD= BODTRI(1,DUM); ITRI= BODTRI(2,DUM) - ENDIF - DEALLOCATE(BODTRI,AREA) - ENDIF -ENDIF +CALL GET_BLOCKING_CUTCELL_DONOR(NM,ICC,JCC,IBOD,ITRI) ! For cut-cell ICC, JCC run through its boundary faces and generate new boundary EDGES, CUT-FACES and cells: BLOCK_PHASE_IF : IF(BLOCK_PHASE==1) THEN @@ -6879,7 +7247,8 @@ SUBROUTINE BLOCK_CUT_CELL(NM,ICC,JCC,BLOCK_PHASE) ! Scheme: ! 0. Add REG and CFGAS cut edges as INB cut edges for the normal faces where it corresponds: DUM=0; IF(ALLOCATED(M%CUT_FACE(IFC1)%EDGE_LIST)) DUM=SIZE(M%CUT_FACE(IFC1)%EDGE_LIST,DIM=2) - ALLOCATE(EDGE_LIST_AUX(3,DUM+M%CUT_FACE(IFCX)%CEDGES(1,JFCX))); + DUM2 = M%CUT_FACE(IFCX)%CEDGES(1,JFCX); IF(DUM2<0) DUM2 = 0 + ALLOCATE(EDGE_LIST_AUX(3,DUM+DUM2)); EDGE_LIST_AUX = CC_UNDEFINED; EDGE_LIST_REG(1,:) = CC_ETYPE_CFINB ! Initialize EDGE_LIST addition. IF(DUM>0) EDGE_LIST_AUX(1:3,1:DUM) = M%CUT_FACE(IFC1)%EDGE_LIST(1:3,1:DUM) ALLOCATE(CEDGES_AUX(SIZE(M%CUT_FACE(IFC1)%CFELEM,DIM=1),SIZE(M%CUT_FACE(IFC1)%CFELEM,DIM=2))) @@ -8182,6 +8551,9 @@ SUBROUTINE DROP_CUTFACE(NM,FTYPE,I,J,K,ILHF,X1AXIS,IFC,JFC) ICC1 = CF%CELL_LIST(2,ILH,IND(DUM)) JCC1 = CF%CELL_LIST(3,ILH,IND(DUM)) IFC1 = CF%CELL_LIST(4,ILH,IND(DUM)) + ! One side may have already been dropped to CC_SOLID (CELL_LIST entries set to + ! CC_UNDEFINED). Skip those entries safely; no cell exists on that side to reindex. + IF (ICC1<1 .OR. JCC1<1 .OR. IFC1<1) CYCLE IFACE= M%CUT_CELL(ICC1)%CCELEM(IFC1+1,JCC1) ! Dropping gas-cut cells, do not reindex local JCF for INBOUNDARY faces. These have been changed already. IF(FTYPE==CC_FTYPE_CFINB .OR. (FTYPE==CC_FTYPE_CFGAS .AND. M%CUT_CELL(ICC1)%FACE_LIST(1,IFACE)/=CC_FTYPE_CFINB)) & @@ -8432,6 +8804,7 @@ SUBROUTINE GET_CELL_LINK_INFO(NM) IFACE2 = CC%FACE_LIST(5,IFACE) IBOD = M%CUT_FACE(IFC2)%BODTRI(1,IFACE2) IWSEL = M%CUT_FACE(IFC2)%BODTRI(2,IFACE2) + IF (.NOT.VALID_GEOMETRY_FACE_DONOR(IBOD,IWSEL)) CYCLE AF = M%CUT_FACE(IFC2)%AREA( IFACE2) NRML(IAXIS:KAXIS) = NRML(IAXIS:KAXIS) + GEOMETRY(IBOD)%FACES_NORMAL(IAXIS:KAXIS,IWSEL)*AF AREA = AREA + AF @@ -8543,6 +8916,7 @@ SUBROUTINE GET_CELL_LINK_INFO(NM) IFACE2 = CC%FACE_LIST(5,IFACE) IBOD = M%CUT_FACE(IFC2)%BODTRI(1,IFACE2) IWSEL = M%CUT_FACE(IFC2)%BODTRI(2,IFACE2) + IF (.NOT.VALID_GEOMETRY_FACE_DONOR(IBOD,IWSEL)) CYCLE AF = M%CUT_FACE(IFC2)%AREA( IFACE2) NRML(IAXIS:KAXIS) = NRML(IAXIS:KAXIS) + GEOMETRY(IBOD)%FACES_NORMAL(IAXIS:KAXIS,IWSEL)*AF AREA = AREA + AF @@ -8995,7 +9369,13 @@ SUBROUTINE INIT_CFACE_CELL(NM,ICF,IFACE,CFACE_INDEX,SURF_INDEX,STAGE_FLG,IS_INB, ELSE IBOD =CF%BODTRI(1,IFACE) IWSEL=CF%BODTRI(2,IFACE) - BC%NVEC(IAXIS:KAXIS) = GEOMETRY(IBOD)%FACES_NORMAL(IAXIS:KAXIS,IWSEL) + IF (VALID_GEOMETRY_FACE_DONOR(IBOD,IWSEL)) THEN + BC%NVEC(IAXIS:KAXIS) = GEOMETRY(IBOD)%FACES_NORMAL(IAXIS:KAXIS,IWSEL) + ELSEIF (NORM2(BC%NVEC)>TWENTY_EPSILON_EB) THEN + BC%NVEC(IAXIS:KAXIS) = BC%NVEC(IAXIS:KAXIS)/NORM2(BC%NVEC) + ELSE + BC%NVEC(IAXIS:KAXIS) = 0._EB + ENDIF ENDIF X1AXIS = MAXLOC(ABS(BC%NVEC(IAXIS:KAXIS)),DIM=1) BC%IOR = INT(SIGN(1._EB,BC%NVEC(X1AXIS)))*X1AXIS @@ -9055,24 +9435,26 @@ SUBROUTINE INIT_CFACE_CELL(NM,ICF,IFACE,CFACE_INDEX,SURF_INDEX,STAGE_FLG,IS_INB, IF (SF%BACKING==EXPOSED .AND. SF%THERMAL_BC_INDEX==THERMALLY_THICK) THEN IG = CF%BODTRI(1,IFACE) TRI = CF%BODTRI(2,IFACE) - XP(IAXIS:KAXIS) = (/ BC%X, BC%Y, BC%Z /) ! CFACE centroid location. - RDIR(IAXIS:KAXIS)= - GEOMETRY(IG)%FACES_NORMAL(IAXIS:KAXIS,TRI) ! Normal into the body. - TRI_LOOP : DO IWSEL=1,GEOMETRY(IG)%N_FACES - IF (IWSEL==TRI) CYCLE - WSELEM(NOD1:NOD3) = GEOMETRY(IG)%FACES(NODS_WSEL*(IWSEL-1)+1:NODS_WSEL*IWSEL) - ! Triangles NODES coordinates: - V1(IAXIS:KAXIS) = GEOMETRY(IG)%VERTS(MAX_DIM*(WSELEM(NOD1)-1)+1:MAX_DIM*WSELEM(NOD1)) - V2(IAXIS:KAXIS) = GEOMETRY(IG)%VERTS(MAX_DIM*(WSELEM(NOD2)-1)+1:MAX_DIM*WSELEM(NOD2)) - V3(IAXIS:KAXIS) = GEOMETRY(IG)%VERTS(MAX_DIM*(WSELEM(NOD3)-1)+1:MAX_DIM*WSELEM(NOD3)) - - ! Fast triangle discard method: To do. - - ! Search for intersection point in POS(IAXIS:KAXIS): - CALL RAY_TRIANGLE_INTERSECT_PT(V1,V2,V3,XP,RDIR,IS_INTERSECT,POS) - - IF (IS_INTERSECT) EXIT TRI_LOOP - - ENDDO TRI_LOOP + IF (VALID_GEOMETRY_FACE_DONOR(IG,TRI)) THEN + XP(IAXIS:KAXIS) = (/ BC%X, BC%Y, BC%Z /) ! CFACE centroid location. + RDIR(IAXIS:KAXIS)= - GEOMETRY(IG)%FACES_NORMAL(IAXIS:KAXIS,TRI) ! Normal into the body. + TRI_LOOP : DO IWSEL=1,GEOMETRY(IG)%N_FACES + IF (IWSEL==TRI) CYCLE + WSELEM(NOD1:NOD3) = GEOMETRY(IG)%FACES(NODS_WSEL*(IWSEL-1)+1:NODS_WSEL*IWSEL) + ! Triangles NODES coordinates: + V1(IAXIS:KAXIS) = GEOMETRY(IG)%VERTS(MAX_DIM*(WSELEM(NOD1)-1)+1:MAX_DIM*WSELEM(NOD1)) + V2(IAXIS:KAXIS) = GEOMETRY(IG)%VERTS(MAX_DIM*(WSELEM(NOD2)-1)+1:MAX_DIM*WSELEM(NOD2)) + V3(IAXIS:KAXIS) = GEOMETRY(IG)%VERTS(MAX_DIM*(WSELEM(NOD3)-1)+1:MAX_DIM*WSELEM(NOD3)) + + ! Search for intersection point in POS(IAXIS:KAXIS): + CALL RAY_TRIANGLE_INTERSECT_PT(V1,V2,V3,XP,RDIR,IS_INTERSECT,POS) + + IF (IS_INTERSECT) EXIT TRI_LOOP + + ENDDO TRI_LOOP + ELSE + IS_INTERSECT = .FALSE. + ENDIF IF (IS_INTERSECT) THEN @@ -19148,14 +19530,14 @@ SUBROUTINE GET_CARTCELL_CUTFACES(NM,ISTR,IEND,JSTR,JEND,KSTR,KEND,BNDINT_FLAG) CYCLE ENDIF enddo - IF (.NOT.INLIST) THEN + IF (.NOT.INLIST .AND. SEG_CELL(6,ISEG)>0 .AND. SEG_CELL(4,ISEG)>0) THEN ! Add first triang to list: NBODTRI = NBODTRI + 1 BOD_TRI(1:2,NBODTRI) = SEG_CELL( (/ 6, 4 /) , ISEG) ENDIF ! No second triangle associated: - IF ( SEG_CELL(3,ISEG) < 2 ) CYCLE + IF ( SEG_CELL(3,ISEG) < 2 .OR. SEG_CELL(6,ISEG)<1 .OR. SEG_CELL(5,ISEG)<1 ) CYCLE ! Second triangle location INLIST = .FALSE. @@ -20357,7 +20739,7 @@ SUBROUTINE EAR_CLIP_CFACES(SIZE_CEELEM_SEG_CELL,NSEG,SEG_CELL,XYZVERT,& ENDIF IF(.NOT.FOUND_ISEG1) CYCLE - TRI = 0 + TRI = 0; BOD = 0 ! Test if triangle given by ISEG ISEG+1 DIAG is valid. ! First, drop if Body not the same: IF ( (COUNTEXT<3) .AND. (SEG_CELL2(6,ISEG)/=SEG_CELL2(6,ISEG1)) ) CYCLE @@ -20391,7 +20773,7 @@ SUBROUTINE EAR_CLIP_CFACES(SIZE_CEELEM_SEG_CELL,NSEG,SEG_CELL,XYZVERT,& ENDIF ENDIF - IF ( TRI == 0 ) THEN + IF ( BOD<1 .OR. TRI<1 ) THEN CYCLE ELSE ! Found two segments with matching triangle. @@ -21483,6 +21865,7 @@ SUBROUTINE CUT_CELL_MOVE(CUT_CELL_FROM,CUT_CELL_TO) CALL MOVE_ALLOC(FROM=CUT_CELL_FROM%DEL_RHO_D_DEL_Z_VOL ,TO=CUT_CELL_TO%DEL_RHO_D_DEL_Z_VOL) CALL MOVE_ALLOC(FROM=CUT_CELL_FROM%U_DOT_DEL_RHO_Z_VOL ,TO=CUT_CELL_TO%U_DOT_DEL_RHO_Z_VOL) CALL MOVE_ALLOC(FROM=CUT_CELL_FROM%NOADVANCE ,TO=CUT_CELL_TO%NOADVANCE) +CALL MOVE_ALLOC(FROM=CUT_CELL_FROM%BODTRI_DONOR ,TO=CUT_CELL_TO%BODTRI_DONOR) CALL MOVE_ALLOC(FROM=CUT_CELL_FROM%NOMICC ,TO=CUT_CELL_TO%NOMICC) RETURN @@ -21503,6 +21886,7 @@ SUBROUTINE CELL_DEALLOC(NM,ICC) DEALLOCATE(MESHES(NM)%CUT_CELL(ICC)%IJK_LINK,MESHES(NM)%CUT_CELL(ICC)%LINK_LEV) DEALLOCATE(MESHES(NM)%CUT_CELL(ICC)%VOLUME, MESHES(NM)%CUT_CELL(ICC)%XYZCEN) DEALLOCATE(MESHES(NM)%CUT_CELL(ICC)%NOADVANCE,MESHES(NM)%CUT_CELL(ICC)%UNKZ) +IF (ALLOCATED(MESHES(NM)%CUT_CELL(ICC)%BODTRI_DONOR)) DEALLOCATE(MESHES(NM)%CUT_CELL(ICC)%BODTRI_DONOR) RETURN @@ -21523,9 +21907,11 @@ SUBROUTINE NEW_CELL_ALLOC(NM,ICC,NCELL,NFACE_CELL,NCFACE_CUTCELL) MESHES(NM)%CUT_CELL(ICC)%IJK_LINK = CC_UNDEFINED MESHES(NM)%CUT_CELL(ICC)%LINK_LEV = 0 ! Root of link Hierarchy is zero. -ALLOCATE(MESHES(NM)%CUT_CELL(ICC)%VOLUME(1:NCELL),MESHES(NM)%CUT_CELL(ICC)%NOADVANCE(1:NCELL)) +ALLOCATE(MESHES(NM)%CUT_CELL(ICC)%VOLUME(1:NCELL),MESHES(NM)%CUT_CELL(ICC)%NOADVANCE(1:NCELL), & + MESHES(NM)%CUT_CELL(ICC)%BODTRI_DONOR(1:2,1:NCELL)) MESHES(NM)%CUT_CELL(ICC)%VOLUME = 0._EB MESHES(NM)%CUT_CELL(ICC)%NOADVANCE= NOT_BLOCKED +MESHES(NM)%CUT_CELL(ICC)%BODTRI_DONOR = 0 ALLOCATE(MESHES(NM)%CUT_CELL(ICC)%XYZCEN(IAXIS:KAXIS,1:NCELL)) MESHES(NM)%CUT_CELL(ICC)%XYZCEN = 0._EB diff --git a/Source/type.f90 b/Source/type.f90 index e0d3113834..2f902e1c2b 100644 --- a/Source/type.f90 +++ b/Source/type.f90 @@ -1389,6 +1389,7 @@ MODULE TYPES REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: XYZCEN !< Cut-cell centroid locations. (IAXIS:KAXIS,1:NCELL) INTEGER, ALLOCATABLE, DIMENSION(:) :: UNKZ !< Cut-cells unknown number for scalars. INTEGER, ALLOCATABLE, DIMENSION(:) :: NOADVANCE !< Array to define if cut-cell should be blocked. (1:NCELL) + INTEGER, ALLOCATABLE, DIMENSION(:,:) :: BODTRI_DONOR !< Donor body/triangle for blocked-cell generated inboundary faces. (1:2,1:NCELL) INTEGER :: N_NOMICC=0 !< Number of entries in NOMICC INTEGER, ALLOCATABLE, DIMENSION(:,:) :: NOMICC !< OMESH cut-cells array. (1:2,1:N_NOMICC)