From 941f4abd9ab0e25bac0c4ec1bdb3a054e9f44199 Mon Sep 17 00:00:00 2001 From: marcosvanella Date: Thu, 5 Feb 2026 12:22:03 -0500 Subject: [PATCH] FDS Source: Define CFACE RDN from linked-cell bounding box normal projection. --- Source/ccib.f90 | 228 +++++++++++++++++++++++++++++++++++++++--------- Source/cons.f90 | 1 + Source/geom.f90 | 2 +- 3 files changed, 189 insertions(+), 42 deletions(-) diff --git a/Source/ccib.f90 b/Source/ccib.f90 index 615cf8bb2f..aec0001097 100644 --- a/Source/ccib.f90 +++ b/Source/ccib.f90 @@ -54,6 +54,7 @@ MODULE CC_SCALARS INTEGER, ALLOCATABLE, DIMENSION(:) :: JM_MAT_Z REAL(EB), ALLOCATABLE, DIMENSION(:) :: F_Z, RZ_Z, RZ_ZS, P_0_CV, TMP_0_CV, RHO_0_CV, ZCEN_CV +REAL(EB), ALLOCATABLE, DIMENSION(:,:):: DELTA_UNKZ REAL(EB), ALLOCATABLE, DIMENSION(:,:):: F_Z0, RZ_Z0 REAL(EB), DIMENSION(0:3,0:3,0:3) :: U_TEMP,F_TEMP @@ -587,9 +588,10 @@ SUBROUTINE CHECK_CFLVN_LINKED_CELLS(NM,DT,UVWMAX,R_DX2,MUTRM) ! Local variables: INTEGER :: I,J,K,ICC,JCC,IROW,IMAX,X1AXIS,ILH,IRC,IFC,IFACE,IFC2,IFACE2,ICFA,LOWHIGH -REAL(EB):: MU_TMP,MURDN,AF,VELN,CFLMAX_TMP,VNMAX_TMP,TWOD_FCT +REAL(EB):: MU_TMP,MURDN,AF,VELN,CFLMAX_TMP,VNMAX_TMP,TWOD_FCT,MU_AVG INTEGER, ALLOCATABLE, DIMENSION(:,:) :: IJKT REAL(EB), ALLOCATABLE, DIMENSION(:) :: UVWA, MUV, MURA, VOL, DIVG +INTEGER, PARAMETER :: VN_DXN_METHOD = 2 ! 1: face-based sum, 2: face-based sum + dx/dy/dz in large CVs. IF(NUNKZ_LOC(NM)<1) RETURN @@ -776,9 +778,25 @@ SUBROUTINE CHECK_CFLVN_LINKED_CELLS(NM,DT,UVWMAX,R_DX2,MUTRM) ENDIF IF (CHECK_VN) THEN - DO IROW=UNKZ_ILC(NM)+1,UNKZ_ILC(NM)+NUNKZ_LOC(NM) - MURA(IROW) = MURA(IROW)/VOL(IROW) - ENDDO + IF (VN_DXN_METHOD==2) THEN + DO IROW=UNKZ_ILC(NM)+1,UNKZ_ILC(NM)+NUNKZ_LOC(NM) + I = IJKT(IAXIS,IROW); J = IJKT(JAXIS,IROW); K = IJKT(KAXIS,IROW) + IF (VOL(IROW) >= DEFAULT_VOLFRAC_LINK*DX(I)*DY(J)*DZ(K)) THEN ! For large linked CVs relax VN criterion. + MU_AVG = MUV(IROW)/VOL(IROW) + IF (TWO_D) THEN + MURA(IROW) = 2._EB*MU_AVG*(1._EB/DX(I)**2 + 1._EB/DZ(K)**2) + ELSE + MURA(IROW) = 2._EB*MU_AVG*(1._EB/DX(I)**2 + 1._EB/DY(J)**2 + 1._EB/DZ(K)**2) + ENDIF + ELSE + MURA(IROW) = MURA(IROW)/VOL(IROW) + ENDIF + ENDDO + ELSE + DO IROW=UNKZ_ILC(NM)+1,UNKZ_ILC(NM)+NUNKZ_LOC(NM) + MURA(IROW) = MURA(IROW)/VOL(IROW) + ENDDO + ENDIF IMAX = UNKZ_ILC(NM)+MAXLOC(MURA,DIM=1) VNMAX_TMP = DT * MURA(IMAX) IF(VNMAX_TMP > VN) THEN @@ -7133,57 +7151,77 @@ SUBROUTINE SET_CFACES_P1_RDN INTEGER :: ICF, IFACE INTEGER :: ICC, JCC, I, J, K INTEGER :: IFACE_CELL, ICF_CELL, IROW +INTEGER :: IROW_DELTA REAL(EB):: AREAI REAL(EB), ALLOCATABLE, DIMENSION(:) :: DXN_UNKZ_LOC, AREA_UNKZ_LOC, VOL_UNKZ_LOC +INTEGER, PARAMETER :: RDN_METHOD = 2 ! 1: VOL/SOLIDAREA, 2: bbox projection +REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: XYZMIN_UNKZ, XYZMAX_UNKZ +REAL(EB) :: DELTA_X, DELTA_Y, DELTA_Z, DXN, NVEC_ABS(MAX_DIM) +TYPE(CFACE_TYPE), POINTER :: CFA +TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC ! ALLOCATE local arrays ALLOCATE(DXN_UNKZ_LOC(1:NUNKZ_LOCAL)); DXN_UNKZ_LOC(:) = 0._EB ALLOCATE(AREA_UNKZ_LOC(1:NUNKZ_LOCAL)); AREA_UNKZ_LOC(:) = 0._EB ALLOCATE(VOL_UNKZ_LOC(1:NUNKZ_LOCAL)); VOL_UNKZ_LOC(:) = 0._EB +IF (RDN_METHOD==1) THEN + ! Main Loop: + MESH_LOOP_01 : DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX -! Main Loop: -MESH_LOOP_01 : DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX - - CALL POINT_TO_MESH(NM) + CALL POINT_TO_MESH(NM) - ! Do a volume weighted average of distance to wall from linked cells, if one of them is a regular cell use 1/2 the - ! distance of corner to corner sqrt(DX^2+DY^2+DZ^2). - ! 1. Regular GASPHASE cells within the cc-region: - DO K=1,KBAR - DO J=1,JBAR - DO I=1,IBAR - IF (CCVAR(I,J,K,CC_UNKZ) <= 0 ) CYCLE ! Drop if regular gas cell has not been assigned unknown number. - IROW = CCVAR(I,J,K,CC_UNKZ) - UNKZ_IND(NM_START) ! All row indexes must refer to ind_loc. - VOL_UNKZ_LOC(IROW) = VOL_UNKZ_LOC(IROW) + (DX(I)*DY(J)*DZ(K)) + ! Do a volume weighted average of distance to wall from linked cells, if one of them is a regular cell use 1/2 the + ! distance of corner to corner sqrt(DX^2+DY^2+DZ^2). + ! 1. Regular GASPHASE cells within the cc-region: + DO K=1,KBAR + DO J=1,JBAR + DO I=1,IBAR + IF (CCVAR(I,J,K,CC_UNKZ) <= 0 ) CYCLE ! Drop if regular gas cell has not been assigned unknown number. + IROW = CCVAR(I,J,K,CC_UNKZ) - UNKZ_IND(NM_START) ! All row indexes must refer to ind_loc. + VOL_UNKZ_LOC(IROW) = VOL_UNKZ_LOC(IROW) + (DX(I)*DY(J)*DZ(K)) + ENDDO ENDDO ENDDO - ENDDO - ! 2. Cut-cells: - 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 - DO JCC=1,CC%NCELL - IROW = CC%UNKZ(JCC) - UNKZ_IND(NM_START) - ! Mean INBOUNDARY cut-face distance to this cut-cell center, projected to cut-face normal: - AREAI = 0._EB - DO ICF_CELL=1,CC%CCELEM(1,JCC) - IFACE_CELL = CC%CCELEM(ICF_CELL+1,JCC) - IF (CC%FACE_LIST(1,IFACE_CELL) /= CC_FTYPE_CFINB) CYCLE - ! Indexes of INBOUNDARY cutface on CUT_FACE: - ICF = CC%FACE_LIST(4,IFACE_CELL) - IFACE = CC%FACE_LIST(5,IFACE_CELL) - ! Area sum: - AREAI = AREAI + CUT_FACE(ICF)%AREA(IFACE) + ! 2. Cut-cells: + 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 + DO JCC=1,CC%NCELL + IROW = CC%UNKZ(JCC) - UNKZ_IND(NM_START) + ! Mean INBOUNDARY cut-face distance to this cut-cell center, projected to cut-face normal: + AREAI = 0._EB + DO ICF_CELL=1,CC%CCELEM(1,JCC) + IFACE_CELL = CC%CCELEM(ICF_CELL+1,JCC) + IF (CC%FACE_LIST(1,IFACE_CELL) /= CC_FTYPE_CFINB) CYCLE + ! Indexes of INBOUNDARY cutface on CUT_FACE: + ICF = CC%FACE_LIST(4,IFACE_CELL) + IFACE = CC%FACE_LIST(5,IFACE_CELL) + ! Area sum: + AREAI = AREAI + CUT_FACE(ICF)%AREA(IFACE) + ENDDO + AREA_UNKZ_LOC(IROW) = AREA_UNKZ_LOC(IROW) + AREAI + VOL_UNKZ_LOC(IROW) = VOL_UNKZ_LOC(IROW) + CC%VOLUME(JCC) ENDDO - AREA_UNKZ_LOC(IROW) = AREA_UNKZ_LOC(IROW) + AREAI - VOL_UNKZ_LOC(IROW) = VOL_UNKZ_LOC(IROW) + CC%VOLUME(JCC) ENDDO - ENDDO -ENDDO MESH_LOOP_01 + ENDDO MESH_LOOP_01 -! Compute average DXN of all linked cells: -DXN_UNKZ_LOC = VOL_UNKZ_LOC / (AREA_UNKZ_LOC + TWENTY_EPSILON_EB) + ! Compute average DXN of all linked cells: + DXN_UNKZ_LOC = VOL_UNKZ_LOC / (AREA_UNKZ_LOC + TWO_EPSILON_EB) +ELSE + ALLOCATE(XYZMIN_UNKZ(IAXIS:KAXIS,1:NUNKZ_LOCAL)) + ALLOCATE(XYZMAX_UNKZ(IAXIS:KAXIS,1:NUNKZ_LOCAL)) + CALL GET_LINKED_CELLS_BBOX(XYZMIN_UNKZ,XYZMAX_UNKZ) + IF (ALLOCATED(DELTA_UNKZ)) THEN + IF (SIZE(DELTA_UNKZ,DIM=2) /= NUNKZ_LOCAL) DEALLOCATE(DELTA_UNKZ) + ENDIF + IF (.NOT.ALLOCATED(DELTA_UNKZ)) ALLOCATE(DELTA_UNKZ(IAXIS:KAXIS,1:NUNKZ_LOCAL)) + DO IROW_DELTA=1,NUNKZ_LOCAL + DELTA_UNKZ(IAXIS,IROW_DELTA) = MAX(0._EB,XYZMAX_UNKZ(IAXIS,IROW_DELTA)-XYZMIN_UNKZ(IAXIS,IROW_DELTA)) + DELTA_UNKZ(JAXIS,IROW_DELTA) = MAX(0._EB,XYZMAX_UNKZ(JAXIS,IROW_DELTA)-XYZMIN_UNKZ(JAXIS,IROW_DELTA)) + DELTA_UNKZ(KAXIS,IROW_DELTA) = MAX(0._EB,XYZMAX_UNKZ(KAXIS,IROW_DELTA)-XYZMIN_UNKZ(KAXIS,IROW_DELTA)) + ENDDO +ENDIF ! Finally Define B1%RDN: MESH_LOOP_02 : DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX @@ -7196,15 +7234,123 @@ SUBROUTINE SET_CFACES_P1_RDN ICC = CF%CELL_LIST(2,LOW_IND,IFACE) JCC = CF%CELL_LIST(3,LOW_IND,IFACE) IROW = CUT_CELL(ICC)%UNKZ(JCC) - UNKZ_IND(NM_START) - BOUNDARY_PROP1(CFACE(CF%CFACE_INDEX(IFACE))%B1_INDEX)%RDN = 1._EB/DXN_UNKZ_LOC(IROW) + SELECT CASE (RDN_METHOD) + CASE (1) + BOUNDARY_PROP1(CFACE(CF%CFACE_INDEX(IFACE))%B1_INDEX)%RDN = 1._EB/DXN_UNKZ_LOC(IROW) + CASE (2) + CFA => CFACE(CF%CFACE_INDEX(IFACE)) + BC => BOUNDARY_COORD(CFA%BC_INDEX) + DELTA_X = XYZMAX_UNKZ(IAXIS,IROW) - XYZMIN_UNKZ(IAXIS,IROW) + DELTA_Y = XYZMAX_UNKZ(JAXIS,IROW) - XYZMIN_UNKZ(JAXIS,IROW) + DELTA_Z = XYZMAX_UNKZ(KAXIS,IROW) - XYZMIN_UNKZ(KAXIS,IROW) + NVEC_ABS(IAXIS:KAXIS) = ABS(BC%NVEC(IAXIS:KAXIS)) + DXN = DELTA_X*NVEC_ABS(IAXIS) + DELTA_Y*NVEC_ABS(JAXIS) + DELTA_Z*NVEC_ABS(KAXIS) + BOUNDARY_PROP1(CFA%B1_INDEX)%RDN = 1._EB/(DXN + TWO_EPSILON_EB) + END SELECT ENDDO ENDDO ENDDO MESH_LOOP_02 +IF (ALLOCATED(XYZMIN_UNKZ)) DEALLOCATE(XYZMIN_UNKZ, XYZMAX_UNKZ) DEALLOCATE(DXN_UNKZ_LOC, VOL_UNKZ_LOC, AREA_UNKZ_LOC) RETURN END SUBROUTINE SET_CFACES_P1_RDN +! ---------------------- GET_LINKED_CELLS_BBOX -------------------------------- + +SUBROUTINE GET_LINKED_CELLS_BBOX(XYZMIN_UNKZ,XYZMAX_UNKZ) + +! Local Variables: +INTEGER :: I, J, K, ICC, JCC, IROW +INTEGER :: IFACE, IFC, FTYPE, X1AXIS, LOWHIGH, ILH, IFC2, IFACE2, NVFACE, IPT, IVERT +REAL(EB) :: XLO, XHI, YLO, YHI, ZLO, ZHI +REAL(EB), INTENT(OUT) :: XYZMIN_UNKZ(IAXIS:KAXIS,1:NUNKZ_LOCAL) +REAL(EB), INTENT(OUT) :: XYZMAX_UNKZ(IAXIS:KAXIS,1:NUNKZ_LOCAL) + +XYZMIN_UNKZ = HUGE(1._EB); XYZMAX_UNKZ = -HUGE(1._EB) +MESH_LOOP_01 : DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX + CALL POINT_TO_MESH(NM) + DO K=1,KBAR + DO J=1,JBAR + DO I=1,IBAR + IF (CCVAR(I,J,K,CC_UNKZ) <= 0 ) CYCLE + IROW = CCVAR(I,J,K,CC_UNKZ) - UNKZ_IND(NM_START) + XLO = X(I-1); XHI = X(I) + YLO = Y(J-1); YHI = Y(J) + ZLO = Z(K-1); ZHI = Z(K) + XYZMIN_UNKZ(IAXIS,IROW) = MIN(XYZMIN_UNKZ(IAXIS,IROW),XLO) + XYZMIN_UNKZ(JAXIS,IROW) = MIN(XYZMIN_UNKZ(JAXIS,IROW),YLO) + XYZMIN_UNKZ(KAXIS,IROW) = MIN(XYZMIN_UNKZ(KAXIS,IROW),ZLO) + XYZMAX_UNKZ(IAXIS,IROW) = MAX(XYZMAX_UNKZ(IAXIS,IROW),XHI) + XYZMAX_UNKZ(JAXIS,IROW) = MAX(XYZMAX_UNKZ(JAXIS,IROW),YHI) + XYZMAX_UNKZ(KAXIS,IROW) = MAX(XYZMAX_UNKZ(KAXIS,IROW),ZHI) + ENDDO + ENDDO + ENDDO + + 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 + DO JCC=1,CC%NCELL + IROW = CC%UNKZ(JCC) - UNKZ_IND(NM_START) + XLO = HUGE(1._EB); XHI = -HUGE(1._EB) + YLO = HUGE(1._EB); YHI = -HUGE(1._EB) + ZLO = HUGE(1._EB); ZHI = -HUGE(1._EB) + DO IFC=1,CC%CCELEM(1,JCC) + IFACE = CC%CCELEM(IFC+1,JCC) + FTYPE = CC%FACE_LIST(1,IFACE) + SELECT CASE (FTYPE) + CASE (CC_FTYPE_CFGAS,CC_FTYPE_CFINB) + IFC2 = CC%FACE_LIST(4,IFACE) + IFACE2 = CC%FACE_LIST(5,IFACE) + NVFACE = CUT_FACE(IFC2)%CFELEM(1,IFACE2) + DO IPT=1,NVFACE + IVERT = CUT_FACE(IFC2)%CFELEM(IPT+1,IFACE2) + XLO = MIN(XLO,CUT_FACE(IFC2)%XYZVERT(IAXIS,IVERT)) + XHI = MAX(XHI,CUT_FACE(IFC2)%XYZVERT(IAXIS,IVERT)) + YLO = MIN(YLO,CUT_FACE(IFC2)%XYZVERT(JAXIS,IVERT)) + YHI = MAX(YHI,CUT_FACE(IFC2)%XYZVERT(JAXIS,IVERT)) + ZLO = MIN(ZLO,CUT_FACE(IFC2)%XYZVERT(KAXIS,IVERT)) + ZHI = MAX(ZHI,CUT_FACE(IFC2)%XYZVERT(KAXIS,IVERT)) + ENDDO + CASE (CC_FTYPE_RCGAS) + LOWHIGH = CC%FACE_LIST(2,IFACE) + X1AXIS = CC%FACE_LIST(3,IFACE) + ILH = LOWHIGH - 1 + SELECT CASE (X1AXIS) + CASE (IAXIS) + XLO = MIN(XLO,X(I-1+ILH)); XHI = MAX(XHI,X(I-1+ILH)) + YLO = MIN(YLO,Y(J-1)); YHI = MAX(YHI,Y(J)) + ZLO = MIN(ZLO,Z(K-1)); ZHI = MAX(ZHI,Z(K)) + CASE (JAXIS) + XLO = MIN(XLO,X(I-1)); XHI = MAX(XHI,X(I)) + YLO = MIN(YLO,Y(J-1+ILH)); YHI = MAX(YHI,Y(J-1+ILH)) + ZLO = MIN(ZLO,Z(K-1)); ZHI = MAX(ZHI,Z(K)) + CASE (KAXIS) + XLO = MIN(XLO,X(I-1)); XHI = MAX(XHI,X(I)) + YLO = MIN(YLO,Y(J-1)); YHI = MAX(YHI,Y(J)) + ZLO = MIN(ZLO,Z(K-1+ILH)); ZHI = MAX(ZHI,Z(K-1+ILH)) + END SELECT + END SELECT + ENDDO + IF (XLO > XHI) THEN + XLO = X(I-1); XHI = X(I) + YLO = Y(J-1); YHI = Y(J) + ZLO = Z(K-1); ZHI = Z(K) + ENDIF + XYZMIN_UNKZ(IAXIS,IROW) = MIN(XYZMIN_UNKZ(IAXIS,IROW),XLO) + XYZMIN_UNKZ(JAXIS,IROW) = MIN(XYZMIN_UNKZ(JAXIS,IROW),YLO) + XYZMIN_UNKZ(KAXIS,IROW) = MIN(XYZMIN_UNKZ(KAXIS,IROW),ZLO) + XYZMAX_UNKZ(IAXIS,IROW) = MAX(XYZMAX_UNKZ(IAXIS,IROW),XHI) + XYZMAX_UNKZ(JAXIS,IROW) = MAX(XYZMAX_UNKZ(JAXIS,IROW),YHI) + XYZMAX_UNKZ(KAXIS,IROW) = MAX(XYZMAX_UNKZ(KAXIS,IROW),ZHI) + ENDDO + ENDDO +ENDDO MESH_LOOP_01 + +RETURN +END SUBROUTINE GET_LINKED_CELLS_BBOX + END SUBROUTINE CC_SET_DATA ! ------------------------------- CC_END_STEP -------------------------------- diff --git a/Source/cons.f90 b/Source/cons.f90 index a164bcc623..f3c0b24db9 100644 --- a/Source/cons.f90 +++ b/Source/cons.f90 @@ -697,6 +697,7 @@ MODULE GLOBAL_CONSTANTS ! Threshold factor for volume of cut-cells respect to volume of Cartesian cells: ! Currently used in the thermo div definition of cut-cells. +REAL(EB), PARAMETER :: DEFAULT_VOLFRAC_LINK = 0.5_EB REAL(EB) :: CCVOL_LINK=0.95_EB LOGICAL :: GET_CUTCELLS_VERBOSE=.FALSE. diff --git a/Source/geom.f90 b/Source/geom.f90 index f3292f03f1..69f76ecc30 100644 --- a/Source/geom.f90 +++ b/Source/geom.f90 @@ -23781,7 +23781,7 @@ SUBROUTINE READ_GEOM PRES_FLAG = ULMAT_FLAG ENDIF PRES_ON_WHOLE_DOMAIN = .FALSE. -IF (ABS(CCVOL_LINK-0.95_EB)