Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

FDS Source : Clean PRES_ON_WHOLE_DOMAIN, ccib.f90. #13942

Merged
merged 2 commits into from
Dec 24, 2024
Merged

FDS Source : Clean PRES_ON_WHOLE_DOMAIN, ccib.f90. #13942

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
2f23462
FDS Source : Clean PRES_ON_WHOLE_DOMAIN, ccib.f90.
marcosvanella Dec 23, 2024
40e3576
Merge branch 'firemodels:master' into master
marcosvanella Dec 23, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
227 changes: 36 additions & 191 deletions Source/ccib.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2469,29 +2469,9 @@ SUBROUTINE CC_PROJECT_VELOCITY(NM,DT,STORE_FLG)
END SELECT
ENDDO WALL_CELL_LOOP_1

IF(.NOT.PRES_ON_WHOLE_DOMAIN) THEN
DO K=1,KBAR
DO J=1,JBAR
DO I=0,IBAR
IF(FCVAR(I,J,K,CC_FGSC,IAXIS)==CC_SOLID) US(I,J,K) = 0._EB
ENDDO
ENDDO
ENDDO
DO K=1,KBAR
DO J=0,JBAR
DO I=1,IBAR
IF(FCVAR(I,J,K,CC_FGSC,JAXIS)==CC_SOLID) VS(I,J,K) = 0._EB
ENDDO
ENDDO
ENDDO
DO K=0,KBAR
DO J=1,JBAR
DO I=1,IBAR
IF(FCVAR(I,J,K,CC_FGSC,KAXIS)==CC_SOLID) WS(I,J,K) = 0._EB
ENDDO
ENDDO
ENDDO
ENDIF
WHERE(FCVAR(0:IBAR,1:JBAR,1:KBAR,CC_FGSC,IAXIS)==CC_SOLID) US(0:IBAR,1:JBAR,1:KBAR) = 0._EB
WHERE(FCVAR(1:IBAR,0:JBAR,1:KBAR,CC_FGSC,JAXIS)==CC_SOLID) VS(1:IBAR,0:JBAR,1:KBAR) = 0._EB
WHERE(FCVAR(1:IBAR,1:JBAR,0:KBAR,CC_FGSC,KAXIS)==CC_SOLID) WS(1:IBAR,1:JBAR,0:KBAR) = 0._EB

ELSE PRED_CORR_IF

Expand Down Expand Up @@ -2631,29 +2611,9 @@ SUBROUTINE CC_PROJECT_VELOCITY(NM,DT,STORE_FLG)

DEALLOCATE(U_STORE,V_STORE,W_STORE)

IF(.NOT.PRES_ON_WHOLE_DOMAIN) THEN
DO K=1,KBAR
DO J=1,JBAR
DO I=0,IBAR
IF(FCVAR(I,J,K,CC_FGSC,IAXIS)==CC_SOLID) U(I,J,K) = 0._EB
ENDDO
ENDDO
ENDDO
DO K=1,KBAR
DO J=0,JBAR
DO I=1,IBAR
IF(FCVAR(I,J,K,CC_FGSC,JAXIS)==CC_SOLID) V(I,J,K) = 0._EB
ENDDO
ENDDO
ENDDO
DO K=0,KBAR
DO J=1,JBAR
DO I=1,IBAR
IF(FCVAR(I,J,K,CC_FGSC,KAXIS)==CC_SOLID) W(I,J,K) = 0._EB
ENDDO
ENDDO
ENDDO
ENDIF
WHERE(FCVAR(0:IBAR,1:JBAR,1:KBAR,CC_FGSC,IAXIS)==CC_SOLID) U(0:IBAR,1:JBAR,1:KBAR) = 0._EB
WHERE(FCVAR(1:IBAR,0:JBAR,1:KBAR,CC_FGSC,JAXIS)==CC_SOLID) V(1:IBAR,0:JBAR,1:KBAR) = 0._EB
WHERE(FCVAR(1:IBAR,1:JBAR,0:KBAR,CC_FGSC,KAXIS)==CC_SOLID) W(1:IBAR,1:JBAR,0:KBAR) = 0._EB

ENDIF PRED_CORR_IF

Expand Down Expand Up @@ -6810,7 +6770,7 @@ SUBROUTINE CC_H_INTERP
! Local Variables:
REAL(EB), POINTER, DIMENSION(:,:,:) :: UP,VP,WP,HP
INTEGER :: NM, ICC, NCELL, I, J ,K
REAL(EB):: U_IBM, V_IBM, W_IBM, VCRT
REAL(EB):: VCRT
LOGICAL :: VOLFLG

! This routine interpolates H to cut cells/Cartesian cells at the end of step.
Expand Down Expand Up @@ -6866,52 +6826,12 @@ SUBROUTINE CC_H_INTERP
ENDDO ICC_LOOP

! Finally set HP to zero inside immersed solids:
IF (.NOT.PRES_ON_WHOLE_DOMAIN) THEN
DO K=0,KBP1
DO J=0,JBP1
DO I=0,IBP1
IF (MESHES(NM)%CCVAR(I,J,K,CC_CGSC) /= CC_SOLID) CYCLE
HP(I,J,K) = 0._EB
ENDDO
ENDDO
ENDDO
ENDIF

! In case of .NOT. PRES_ON_WHOLE_DOMAIN set velocities on solid faces to zero:
IF (.NOT.PRES_ON_WHOLE_DOMAIN) THEN
! Force U velocities in CC_SOLID faces to zero
U_IBM = 0._EB ! Body doesn't move.
DO K=1,KBAR
DO J=1,JBAR
DO I=0,IBAR
IF (MESHES(NM)%FCVAR(I,J,K,CC_FGSC,IAXIS) /= CC_SOLID ) CYCLE
UP(I,J,K) = U_IBM
ENDDO
ENDDO
ENDDO

! Force V velocities in CC_SOLID faces to zero
V_IBM = 0._EB ! Body doesn't move.
DO K=1,KBAR
DO J=0,JBAR
DO I=1,IBAR
IF (MESHES(NM)%FCVAR(I,J,K,CC_FGSC,JAXIS) /= CC_SOLID ) CYCLE
VP(I,J,K) = V_IBM
ENDDO
ENDDO
ENDDO
WHERE(MESHES(NM)%CCVAR(0:IBP1,0:JBP1,0:KBP1,CC_CGSC)==CC_SOLID) HP(0:IBP1,0:JBP1,0:KBP1) = 0._EB

! Force W velocities in CC_SOLID faces to zero
W_IBM = 0._EB ! Body doesn't move.
DO K=0,KBAR
DO J=1,JBAR
DO I=1,IBAR
IF (MESHES(NM)%FCVAR(I,J,K,CC_FGSC,KAXIS) /= CC_SOLID ) CYCLE
WP(I,J,K) = W_IBM
ENDDO
ENDDO
ENDDO
ENDIF
! Set velocities on solid faces to zero:
WHERE(MESHES(NM)%FCVAR(0:IBAR,1:JBAR,1:KBAR,CC_FGSC,IAXIS)==CC_SOLID) UP(0:IBAR,1:JBAR,1:KBAR) = 0._EB
WHERE(MESHES(NM)%FCVAR(1:IBAR,0:JBAR,1:KBAR,CC_FGSC,JAXIS)==CC_SOLID) VP(1:IBAR,0:JBAR,1:KBAR) = 0._EB
WHERE(MESHES(NM)%FCVAR(1:IBAR,1:JBAR,0:KBAR,CC_FGSC,KAXIS)==CC_SOLID) WP(1:IBAR,1:JBAR,0:KBAR) = 0._EB

NULLIFY(UP,VP,WP,HP)

Expand Down Expand Up @@ -7356,11 +7276,9 @@ SUBROUTINE INIT_CUTCELL_DATA(T,DT,FIRST_CALL)
LOGICAL, INTENT(IN) :: FIRST_CALL

! Local Variables:
INTEGER :: NM,I,J,K,N,ICC,JCC,X1AXIS,NFACE,ICF,IFACE
INTEGER :: NM,I,J,K,N,ICC,JCC,X1AXIS,NFACE,ICF
REAL(EB) TMP_CC,RHO_CC,AREAT,VEL_CF !,Z_CC,TMP_0_CC,P_0_CC
REAL(EB), ALLOCATABLE, DIMENSION(:) :: ZZ_CC
INTEGER :: INDADD, INDF, IFC, IFACE2, ICFC
REAL(EB):: FSCU, AREATOT
INTEGER :: IW,IROW_LOC

REAL(EB) :: TNOW
Expand Down Expand Up @@ -7599,57 +7517,8 @@ SUBROUTINE INIT_CUTCELL_DATA(T,DT,FIRST_CALL)
ENDDO
ENDDO

! Then INBOUNDARY cut-faces:
ICC_LOOP : DO ICC=1,MESHES(NM)%N_CUTCELL_MESH
I = CUT_CELL(ICC)%IJK(IAXIS)
J = CUT_CELL(ICC)%IJK(JAXIS)
K = CUT_CELL(ICC)%IJK(KAXIS)
FSCU = 0._EB

! Loop on cells neighbors and test if they are of type CC_SOLID, if so
! Add to velocity flux:
! X faces
DO INDADD=-1,1,2
INDF = I - 1 + (INDADD+1)/2
IF( FCVAR(INDF,J,K,CC_FGSC,IAXIS) /= CC_SOLID) CYCLE
FSCU = FSCU + REAL(INDADD,EB)*U(INDF,J,K)*DY(J)*DZ(K)
ENDDO
! Y faces
DO INDADD=-1,1,2
INDF = J - 1 + (INDADD+1)/2
IF( FCVAR(I,INDF,K,CC_FGSC,JAXIS) /= CC_SOLID ) CYCLE
FSCU = FSCU + REAL(INDADD,EB)*V(I,INDF,K)*DX(I)*DZ(K)
ENDDO
! Z faces
DO INDADD=-1,1,2
INDF = K - 1 + (INDADD+1)/2
IF( FCVAR(I,J,INDF,CC_FGSC,KAXIS) /= CC_SOLID ) CYCLE
FSCU = FSCU + REAL(INDADD,EB)*W(I,J,INDF)*DX(I)*DY(J)
ENDDO

! Now Define total area of INBOUNDARY cut-faces:
ICF=CCVAR(I,J,K,CC_IDCF);
ICF_COND : IF (ICF > 0) THEN
NFACE = CUT_FACE(ICF)%NFACE
AREATOT = SUM ( CUT_FACE(ICF)%AREA(1:NFACE) )
DO JCC =1,CUT_CELL(ICC)%NCELL
IFC_LOOP : DO IFC=1,CUT_CELL(ICC)%CCELEM(1,JCC)
IFACE = CUT_CELL(ICC)%CCELEM(IFC+1,JCC)
IF (CUT_CELL(ICC)%FACE_LIST(1,IFACE) == CC_FTYPE_CFINB) THEN
IFACE2 = CUT_CELL(ICC)%FACE_LIST(5,IFACE)
ICFC = CUT_FACE(ICF)%CFACE_INDEX(IFACE2)
IF(PRES_ON_WHOLE_DOMAIN) THEN
CUT_FACE(ICF)%VELS(IFACE2) = 1._EB/AREATOT*FSCU ! +ve into the solid Velocity error
CUT_FACE(ICF)%VEL( IFACE2) = 1._EB/AREATOT*FSCU
ELSE
CUT_FACE(ICF)%VELS(IFACE2) = 0._EB
CUT_FACE(ICF)%VEL( IFACE2) = 0._EB
ENDIF
ENDIF
ENDDO IFC_LOOP
ENDDO
ENDIF ICF_COND
ENDDO ICC_LOOP
! INBOUNDARY cut-faces are initialized with 0._EB velocity, that will be changed in
! CFACE_PREDICT_NORMAL_VELOCITY.

ENDIF PERIODIC_TEST_COND

Expand Down Expand Up @@ -11590,7 +11459,7 @@ SUBROUTINE CC_DENSITY_EXPLICIT(T,DT)
F_Z(:) = 0._EB
CALL GET_EXPLICIT_ADVDIFFVECTOR_SCALAR_3D(N)

! Add Advective fluxes due to PRES_ON_WHOLE_DOMAIN for F_Z:
! Add Advective fluxes for F_Z:
CALL GET_ADVDIFFVECTOR_SCALAR_3D(N)

! Here add the reaction source term M_DOT_PPP, treated explicitly:
Expand Down Expand Up @@ -12144,9 +12013,6 @@ SUBROUTINE GET_EXPLICIT_ADVDIFFVECTOR_SCALAR_3D(N)

ENDDO

! Case of PRES_ON_WHOLE_DOMAIN, we have non-zero velocities on INBOUNDARY cut-faces:
! Done on CALL GET_ADVDIFFVECTOR_SCALAR_3D(N)

! Then add (Del rho D Del Z)*dv computed on CCDIVERGENCE_PART_1:
! Loop over regular cells on CC region:
DO K=1,KBAR
Expand Down Expand Up @@ -15800,7 +15666,7 @@ SUBROUTINE CC_CHECK_DIVERGENCE(T,DT,PREDVEL)
ENDIF
ENDDO ICC_LOOP

IF(.NOT.PRES_ON_WHOLE_DOMAIN .AND. STORE_CARTESIAN_DIVERGENCE) THEN
IF(STORE_CARTESIAN_DIVERGENCE) THEN
DO K=1,KBAR
DO J=1,JBAR
DO I=1,IBAR
Expand Down Expand Up @@ -19774,16 +19640,11 @@ SUBROUTINE GET_CC_UNKH(I,J,K,IUNKH)
INTEGER :: ICC

IUNKH = CC_UNDEFINED ! This is < 0.
IF(.NOT.PRES_ON_WHOLE_DOMAIN ) THEN

! Regular gas cell, taken care of before.
! Check cut-cell:
ICC = CCVAR(I,J,K,CC_IDCC)
! If theres is a cut-cell ICC then CUT_CELL(ICC)%UNKH(1) has been populated.
IF (ICC > 0) IUNKH = CUT_CELL(ICC)%UNKH(1)

ENDIF

! Regular gas cell, taken care of before.
! Check cut-cell:
ICC = CCVAR(I,J,K,CC_IDCC)
! If theres is a cut-cell ICC then CUT_CELL(ICC)%UNKH(1) has been populated.
IF (ICC > 0) IUNKH = CUT_CELL(ICC)%UNKH(1)

RETURN
END SUBROUTINE GET_CC_UNKH
Expand All @@ -19798,12 +19659,10 @@ SUBROUTINE GET_CC_IROW(I,J,K,IPZ,IROW)

! Local variable:
INTEGER :: ICC
IROW = CC_UNDEFINED ! This is < 0.
IF(.NOT.PRES_ON_WHOLE_DOMAIN) THEN
ICC = CCVAR(I,J,K,CC_IDCC)
! If theres is a cut-cell ICC then CUT_CELL(ICC)%UNKH(1) has been populated.
IF (ICC > 0) IROW = CUT_CELL(ICC)%UNKH(1) - ZONE_SOLVE(IPZ)%UNKH_IND(NM_START)
ENDIF
IROW = CC_UNDEFINED ! This is < 0.
ICC = CCVAR(I,J,K,CC_IDCC)
! If theres is a cut-cell ICC then CUT_CELL(ICC)%UNKH(1) has been populated.
IF (ICC > 0) IROW = CUT_CELL(ICC)%UNKH(1) - ZONE_SOLVE(IPZ)%UNKH_IND(NM_START)

RETURN
END SUBROUTINE GET_CC_IROW
Expand Down Expand Up @@ -21836,14 +21695,14 @@ SUBROUTINE NUMBER_UNKH_CUTCELLS(FLAG12,NM,IPZ,NUNKH_LC)
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
IF(.NOT.PRES_ON_WHOLE_DOMAIN .AND. ZONE_SOLVE(PRESSURE_ZONE(I,J,K))%CONNECTED_ZONE_PARENT/=IPZ ) CYCLE
IF(ZONE_SOLVE(PRESSURE_ZONE(I,J,K))%CONNECTED_ZONE_PARENT/=IPZ ) CYCLE
NUNKH_LC(NM) = NUNKH_LC(NM) + 1
CUT_CELL(ICC)%UNKH(1:CC%NCELL) = NUNKH_LC(NM)
ENDDO
ELSE
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(.NOT.PRES_ON_WHOLE_DOMAIN .AND. ZONE_SOLVE(PRESSURE_ZONE(I,J,K))%CONNECTED_ZONE_PARENT/=IPZ ) CYCLE
IF(ZONE_SOLVE(PRESSURE_ZONE(I,J,K))%CONNECTED_ZONE_PARENT/=IPZ ) CYCLE
DO JCC=1,CC%NCELL
CUT_CELL(ICC)%UNKH(JCC) = CUT_CELL(ICC)%UNKH(JCC) + ZONE_SOLVE(IPZ)%UNKH_IND(NM)
ENDDO
Expand Down Expand Up @@ -21885,15 +21744,10 @@ SUBROUTINE COPY_CC_HS_TO_UNKH(NM)

ICC=CCVAR(II,JJ,KK,CC_IDCC)

IF(.NOT.PRES_ON_WHOLE_DOMAIN) THEN
IF (ICC > 0) THEN ! Cut-cells on this guard-cell Cartesian cell.
MESHES(NM)%CUT_CELL(ICC)%UNKH(1) = INT(OM%HS(IIO,JJO,KKO))
ELSE
MESHES(NM)%CCVAR(II,JJ,KK,CC_UNKH) = INT(OM%HS(IIO,JJO,KKO))
ENDIF
IF (ICC > 0) THEN ! Cut-cells on this guard-cell Cartesian cell.
MESHES(NM)%CUT_CELL(ICC)%UNKH(1) = INT(OM%HS(IIO,JJO,KKO))
ELSE
MESHES(NM)%CCVAR(II,JJ,KK,CC_UNKH) = INT(OM%HS(IIO,JJO,KKO))
IF (ICC > 0) MESHES(NM)%CUT_CELL(ICC)%UNKH(1) = MESHES(NM)%CCVAR(II,JJ,KK,CC_UNKH)
ENDIF

ENDDO EXTERNAL_WALL_LOOP
Expand All @@ -21920,21 +21774,12 @@ SUBROUTINE COPY_CC_UNKH_TO_HS(NM)
! Local Variables:
INTEGER :: I,J,K,ICC

IF(.NOT.PRES_ON_WHOLE_DOMAIN) THEN
DO ICC=1,MESHES(NM)%N_CUTCELL_MESH
I = MESHES(NM)%CUT_CELL(ICC)%IJK(IAXIS)
J = MESHES(NM)%CUT_CELL(ICC)%IJK(JAXIS)
K = MESHES(NM)%CUT_CELL(ICC)%IJK(KAXIS)
HS(I,J,K)= REAL(MESHES(NM)%CUT_CELL(ICC)%UNKH(1),EB)
ENDDO
ELSE
DO ICC=1,MESHES(NM)%N_CUTCELL_MESH
I = MESHES(NM)%CUT_CELL(ICC)%IJK(IAXIS)
J = MESHES(NM)%CUT_CELL(ICC)%IJK(JAXIS)
K = MESHES(NM)%CUT_CELL(ICC)%IJK(KAXIS)
MESHES(NM)%CUT_CELL(ICC)%UNKH(1) = INT(HS(I,J,K))
ENDDO
ENDIF
DO ICC=1,MESHES(NM)%N_CUTCELL_MESH
I = MESHES(NM)%CUT_CELL(ICC)%IJK(IAXIS)
J = MESHES(NM)%CUT_CELL(ICC)%IJK(JAXIS)
K = MESHES(NM)%CUT_CELL(ICC)%IJK(KAXIS)
HS(I,J,K)= REAL(MESHES(NM)%CUT_CELL(ICC)%UNKH(1),EB)
ENDDO

RETURN
END SUBROUTINE COPY_CC_UNKH_TO_HS
Expand Down
4 changes: 2 additions & 2 deletions Source/geom.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4745,8 +4745,8 @@ SUBROUTINE GET_GEOM_TRIBIN
DO X1AXIS=IAXIS,KAXIS

! Here reduce the X1_LOW to X1_HIGH distance to the smallest of FDS Mesh and connected meshes BBOX or Geometry:
MIN_MESHGEOM = MAX(MINMAX_MESHES( LOW_IND,X1AXIS),G%GEOM_BOX( LOW_IND,X1AXIS))
MAX_MESHGEOM = MIN(MINMAX_MESHES(HIGH_IND,X1AXIS),G%GEOM_BOX(HIGH_IND,X1AXIS))
MIN_MESHGEOM = MAX(MINMAX_MESHES( LOW_IND,X1AXIS),G%GEOM_BOX( LOW_IND,X1AXIS)-G%MEAN_LEDGE)
MAX_MESHGEOM = MIN(MINMAX_MESHES(HIGH_IND,X1AXIS),G%GEOM_BOX(HIGH_IND,X1AXIS)+G%MEAN_LEDGE)
LX1 = MAX_MESHGEOM - MIN_MESHGEOM

! Define number of bins in direction X1AXIS:
Expand Down
Loading