From fd2ddaf7e503a42df82f5b4a0083d993a41bd4a6 Mon Sep 17 00:00:00 2001 From: mcgratta Date: Thu, 7 May 2026 12:18:02 -0400 Subject: [PATCH 1/2] FDS Source: Start over-hauling dump routine --- Source/dump.f90 | 196 +++++++++++++++++++++++++++++------------------- 1 file changed, 117 insertions(+), 79 deletions(-) diff --git a/Source/dump.f90 b/Source/dump.f90 index dc6c598ece..175cfa109a 100644 --- a/Source/dump.f90 +++ b/Source/dump.f90 @@ -31,23 +31,16 @@ MODULE DUMP TYPE (MESH_TYPE), POINTER :: M TYPE (LAGRANGIAN_PARTICLE_TYPE), POINTER :: LP -TYPE (OBSTRUCTION_TYPE), POINTER :: OB TYPE (VENTS_TYPE), POINTER :: VT TYPE (LAGRANGIAN_PARTICLE_CLASS_TYPE), POINTER :: LPC TYPE (SPECIES_TYPE), POINTER :: SS TYPE (REACTION_TYPE), POINTER :: RN TYPE (SURFACE_TYPE),POINTER :: SF TYPE (MATERIAL_TYPE),POINTER :: ML -TYPE (PROPERTY_TYPE), POINTER :: PY -TYPE (DEVICE_TYPE), POINTER :: DV, DV2 TYPE (SUBDEVICE_TYPE), POINTER :: SDV TYPE (SLICE_TYPE), POINTER :: SL TYPE (WALL_TYPE), POINTER :: WC TYPE (THIN_WALL_TYPE), POINTER :: TW -TYPE (CFACE_TYPE), POINTER :: CFA -TYPE (BOUNDARY_FILE_TYPE), POINTER :: BF -TYPE (ISOSURFACE_FILE_TYPE), POINTER :: IS -TYPE (INITIALIZATION_TYPE), POINTER :: IN PUBLIC ASSIGN_FILE_NAMES,INITIALIZE_GLOBAL_DUMPS,INITIALIZE_MESH_DUMPS,WRITE_STATUS_FILES, & TIMINGS,FLUSH_GLOBAL_BUFFERS,READ_RESTART,WRITE_DIAGNOSTICS, & @@ -544,6 +537,8 @@ SUBROUTINE INITIALIZE_GLOBAL_DUMPS(T,DT) CHARACTER(80) :: FN CHARACTER(LABEL_LENGTH) :: LAB,UNITS CHARACTER(LABEL_LENGTH+7), DIMENSION(42) :: LABEL='null' +TYPE (PROPERTY_TYPE), POINTER :: PY +TYPE (DEVICE_TYPE), POINTER :: DV TNOW=CURRENT_TIME() @@ -886,6 +881,9 @@ SUBROUTINE INITIALIZE_MESH_DUMPS(NM) LOGICAL :: OVERLAPPING_X,OVERLAPPING_Y,OVERLAPPING_Z TYPE (BOUNDARY_COORD_TYPE), POINTER :: BC TYPE (RAD_FILE_TYPE), POINTER :: RF +TYPE (OBSTRUCTION_TYPE), POINTER :: OB +TYPE (ISOSURFACE_FILE_TYPE), POINTER :: IS +TYPE (BOUNDARY_FILE_TYPE), POINTER :: BF TNOW=CURRENT_TIME() @@ -1528,6 +1526,10 @@ SUBROUTINE WRITE_SMOKEVIEW_FILE INTEGER, ALLOCATABLE, DIMENSION(:) :: RECV_USE_LEN,RECV_USE_OFF,RECV_COUNTS CHARACTER(LEN=:), ALLOCATABLE :: STR_GATHER INTEGER SHOW_BNDF(-3:3) +TYPE (PROPERTY_TYPE), POINTER :: PY +TYPE (OBSTRUCTION_TYPE), POINTER :: OB +TYPE (BOUNDARY_FILE_TYPE), POINTER :: BF +TYPE (DEVICE_TYPE), POINTER :: DV ! If this is a RESTART case but an old .smv file does not exist, shutdown with an ERROR. @@ -2716,6 +2718,10 @@ SUBROUTINE INITIALIZE_DIAGNOSTIC_FILE(DT) MU_Z,K_Z,CP_ZN,H_Z, PHI_TILDE,TMP_FLAME CHARACTER(LABEL_LENGTH) :: QUANTITY,ODE_SOLVER,OUTFORM TYPE(SPECIES_MIXTURE_TYPE),POINTER :: SM +TYPE(PROPERTY_TYPE), POINTER :: PY +TYPE(ISOSURFACE_FILE_TYPE), POINTER :: IS +TYPE(BOUNDARY_FILE_TYPE), POINTER :: BF +TYPE(DEVICE_TYPE), POINTER :: DV ! Open and initialize diagnostic output file @@ -3529,6 +3535,10 @@ SUBROUTINE DUMP_RESTART(T,DT,NM) TYPE(DUCT_TYPE), POINTER :: DU TYPE(DUCTNODE_TYPE), POINTER :: DN TYPE(STORAGE_TYPE), POINTER :: OS +TYPE(OBSTRUCTION_TYPE), POINTER :: OB +TYPE(INITIALIZATION_TYPE), POINTER :: IN +TYPE(CFACE_TYPE), POINTER :: CFA +TYPE(DEVICE_TYPE), POINTER :: DV OPEN(LU_CORE(NM),FILE=FN_CORE(NM),FORM='UNFORMATTED',STATUS='REPLACE') @@ -3709,6 +3719,10 @@ SUBROUTINE READ_RESTART(T,DT,NM) TYPE(DUCT_TYPE), POINTER :: DU TYPE(DUCTNODE_TYPE), POINTER :: DN TYPE(STORAGE_TYPE), POINTER :: OS +TYPE(OBSTRUCTION_TYPE), POINTER :: OB +TYPE(INITIALIZATION_TYPE), POINTER :: IN +TYPE(CFACE_TYPE), POINTER :: CFA +TYPE(DEVICE_TYPE), POINTER :: DV INQUIRE(FILE=FN_RESTART(NM),EXIST=EX) IF (.NOT.EX) THEN @@ -4262,6 +4276,7 @@ SUBROUTINE DUMP_ISOF(T,DT,NM) REAL(FB) :: ISO_LEVEL(1) INTEGER :: ISO_NLEVEL INTEGER :: II, JJ, KK +TYPE (ISOSURFACE_FILE_TYPE), POINTER :: IS STIME = REAL(T_BEGIN + (T-T_BEGIN)*TIME_SHRINK_FACTOR,FB) DATAFLAG = 1 @@ -4329,7 +4344,8 @@ SUBROUTINE DUMP_ISOF(T,DT,NM) DO K=0,KBP1 DO J=0,JBP1 DO I=0,IBP1 - QUANTITY(I,J,K) = GAS_PHASE_OUTPUT(T,DT,NM,I,J,K,IS%INDEX,IS%Y_INDEX,IS%Z_INDEX,0,0,IS%VELO_INDEX,0,0,0,0) + QUANTITY(I,J,K) = GAS_PHASE_OUTPUT(T,DT,NM,I,J,K,IS%INDEX,Y_INDEX=IS%Y_INDEX,Z_INDEX=IS%Z_INDEX,& + VELO_INDEX=IS%VELO_INDEX) ENDDO ENDDO ENDDO @@ -4380,7 +4396,8 @@ SUBROUTINE DUMP_ISOF(T,DT,NM) DO K=0,KBP1 DO J=0,JBP1 DO I=0,IBP1 - QUANTITY2(I,J,K) = GAS_PHASE_OUTPUT(T,DT,NM,I,J,K,IS%INDEX2,IS%Y_INDEX2,IS%Z_INDEX2,0,0,IS%VELO_INDEX2,0,0,0,0) + QUANTITY2(I,J,K) = GAS_PHASE_OUTPUT(T,DT,NM,I,J,K,IS%INDEX2,Y_INDEX=IS%Y_INDEX2,Z_INDEX=IS%Z_INDEX2,& + VELO_INDEX=IS%VELO_INDEX2) ENDDO ENDDO ENDDO @@ -4462,7 +4479,7 @@ SUBROUTINE DUMP_SMOKE3D(T,DT,NM) DO K=0,KBP1 DO J=0,JBP1 DO I=0,IBP1 - FF(I,J,K)=GAS_PHASE_OUTPUT(T,DT,NM,I,J,K,S3%QUANTITY_INDEX,S3%Y_INDEX,S3%Z_INDEX,0,0,0,0,0,0,0) + FF(I,J,K)=GAS_PHASE_OUTPUT(T,DT,NM,I,J,K,S3%QUANTITY_INDEX,Y_INDEX=S3%Y_INDEX,Z_INDEX=S3%Z_INDEX) ENDDO ENDDO ENDDO @@ -5162,7 +5179,7 @@ END SUBROUTINE GET_GEOMINFO SUBROUTINE GET_GEOMVALS(CC_INTERP2FACES,CC_CELL_CENTERED,SLICETYPE,& I1,I2,J1,J2,K1,K2,NFACES,NFACES_CUTCELLS,VALS,& IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,& - PIPE_INDEX,PROP_INDEX,REAC_INDEX,MATL_INDEX,T,DT,NM,OPT_BNDF_INDEX) + PIPE_INDEX,PROP_INDEX,REAC_INDEX,T,DT,NM,OPT_BNDF_INDEX) USE PHYSICAL_FUNCTIONS, ONLY: GET_MASS_FRACTION @@ -5170,7 +5187,7 @@ SUBROUTINE GET_GEOMVALS(CC_INTERP2FACES,CC_CELL_CENTERED,SLICETYPE,& REAL(EB), INTENT(IN) :: T,DT INTEGER, INTENT(IN) :: I1,I2,J1,J2,K1,K2,NFACES,NFACES_CUTCELLS,& - IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX,MATL_INDEX,NM + IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX,NM INTEGER, OPTIONAL,INTENT(IN) :: OPT_BNDF_INDEX CHARACTER(*), INTENT(IN) :: SLICETYPE LOGICAL, INTENT(IN) :: CC_INTERP2FACES,CC_CELL_CENTERED @@ -5237,7 +5254,7 @@ SUBROUTINE GET_GEOMVALS(CC_INTERP2FACES,CC_CELL_CENTERED,SLICETYPE,& ICF = FCVAR(SLICE,J,K,CC_IDCF,IAXIS) ! is a cut cell DO IFACECF=1,CUT_FACE(ICF)%NFACE CALL GET_GASCUTFACE_SCALAR_SLICE(VAL_CF,X1AXIS,ICF,IFACECF,CC_INTERP2FACES,CC_CELL_CENTERED,& - IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX,MATL_INDEX,T,DT,NM) + IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX,T,DT,NM) NVF=CUT_FACE(ICF)%CFELEM(1,IFACECF) DO IVCF = 1, NVF-2 IFACECUT = IFACECUT + 1 @@ -5245,7 +5262,7 @@ SUBROUTINE GET_GEOMVALS(CC_INTERP2FACES,CC_CELL_CENTERED,SLICETYPE,& ENDDO ENDDO CALL GET_SOLIDCUTFACE_SCALAR_SLICE(X1AXIS,ICF,VAL_CF,& - IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX,MATL_INDEX,T,DT,NM) + IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX,T,DT,NM) DO IFACECF=CUT_FACE(ICF)%NFACE+1,CUT_FACE(ICF)%NFACE+CUT_FACE(ICF)%NSFACE NVF=CUT_FACE(ICF)%CFELEM(1,IFACECF) DO IVCF = 1, NVF-2 @@ -5255,7 +5272,7 @@ SUBROUTINE GET_GEOMVALS(CC_INTERP2FACES,CC_CELL_CENTERED,SLICETYPE,& ENDDO ELSEIF(CELLTYPE == CC_SOLID) THEN CALL GET_SOLIDREGFACE_SCALAR_SLICE(X1AXIS,SLICE,J,K,VAL_CF,& - IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX,MATL_INDEX,T,DT,NM) + IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX,T,DT,NM) IFACE = IFACE + 1 ! is a solid or gas cell VALS(IFACE) = REAL(VAL_CF,FB) @@ -5292,7 +5309,7 @@ SUBROUTINE GET_GEOMVALS(CC_INTERP2FACES,CC_CELL_CENTERED,SLICETYPE,& ICF = FCVAR(I,SLICE,K,CC_IDCF,JAXIS) DO IFACECF=1,CUT_FACE(ICF)%NFACE CALL GET_GASCUTFACE_SCALAR_SLICE(VAL_CF,X1AXIS,ICF,IFACECF,CC_INTERP2FACES,CC_CELL_CENTERED,& - IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX,MATL_INDEX,T,DT,NM) + IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX,T,DT,NM) NVF=CUT_FACE(ICF)%CFELEM(1,IFACECF) DO IVCF = 1, NVF-2 IFACECUT = IFACECUT + 1 @@ -5300,7 +5317,7 @@ SUBROUTINE GET_GEOMVALS(CC_INTERP2FACES,CC_CELL_CENTERED,SLICETYPE,& ENDDO ENDDO CALL GET_SOLIDCUTFACE_SCALAR_SLICE(X1AXIS,ICF,VAL_CF,& - IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX,MATL_INDEX,T,DT,NM) + IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX,T,DT,NM) DO IFACECF=CUT_FACE(ICF)%NFACE+1,CUT_FACE(ICF)%NFACE+CUT_FACE(ICF)%NSFACE NVF=CUT_FACE(ICF)%CFELEM(1,IFACECF) DO IVCF = 1, NVF-2 @@ -5310,7 +5327,7 @@ SUBROUTINE GET_GEOMVALS(CC_INTERP2FACES,CC_CELL_CENTERED,SLICETYPE,& ENDDO ELSEIF(CELLTYPE == CC_SOLID) THEN CALL GET_SOLIDREGFACE_SCALAR_SLICE(X1AXIS,I,SLICE,K,VAL_CF,& - IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX,MATL_INDEX,T,DT,NM) + IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX,T,DT,NM) IFACE = IFACE + 1 ! is a solid or gas cell VALS(IFACE) = REAL(VAL_CF,FB) @@ -5346,7 +5363,7 @@ SUBROUTINE GET_GEOMVALS(CC_INTERP2FACES,CC_CELL_CENTERED,SLICETYPE,& ICF = FCVAR(I,J,SLICE,CC_IDCF,KAXIS) DO IFACECF=1,CUT_FACE(ICF)%NFACE CALL GET_GASCUTFACE_SCALAR_SLICE(VAL_CF,X1AXIS,ICF,IFACECF,CC_INTERP2FACES,CC_CELL_CENTERED,& - IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX,MATL_INDEX,T,DT,NM) + IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX,T,DT,NM) NVF=CUT_FACE(ICF)%CFELEM(1,IFACECF) DO IVCF = 1, NVF-2 IFACECUT = IFACECUT + 1 @@ -5354,7 +5371,7 @@ SUBROUTINE GET_GEOMVALS(CC_INTERP2FACES,CC_CELL_CENTERED,SLICETYPE,& ENDDO ENDDO CALL GET_SOLIDCUTFACE_SCALAR_SLICE(X1AXIS,ICF,VAL_CF,& - IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX,MATL_INDEX,T,DT,NM) + IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX,T,DT,NM) DO IFACECF=CUT_FACE(ICF)%NFACE+1,CUT_FACE(ICF)%NFACE+CUT_FACE(ICF)%NSFACE NVF=CUT_FACE(ICF)%CFELEM(1,IFACECF) DO IVCF = 1, NVF-2 @@ -5364,7 +5381,7 @@ SUBROUTINE GET_GEOMVALS(CC_INTERP2FACES,CC_CELL_CENTERED,SLICETYPE,& ENDDO ELSEIF(CELLTYPE == CC_SOLID) THEN CALL GET_SOLIDREGFACE_SCALAR_SLICE(X1AXIS,I,J,SLICE,VAL_CF,& - IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX,MATL_INDEX,T,DT,NM) + IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX,T,DT,NM) IFACE = IFACE + 1 ! is a solid or gas cell VALS(IFACE) = REAL(VAL_CF,FB) @@ -5603,12 +5620,12 @@ END SUBROUTINE DUMP_CFACES_GEOM SUBROUTINE DUMP_SLICE_GEOM_DATA(FUNIT_DATA,CC_INTERP2FACES,CC_CELL_CENTERED,SLICETYPE, & HEADER,STIME,I1,I2,J1,J2,K1,K2,DEBUG,& - IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX,MATL_INDEX,T, & + IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX,T, & DT,NM,SLICE_MIN, SLICE_MAX, OPT_BNDF_INDEX) REAL(EB), INTENT(IN) :: T,DT CHARACTER(*), INTENT(IN) :: SLICETYPE INTEGER, INTENT(IN) :: FUNIT_DATA,HEADER,I1,I2,J1,J2,K1,K2,DEBUG, & - IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX,MATL_INDEX,NM + IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX,NM INTEGER, OPTIONAL,INTENT(IN) :: OPT_BNDF_INDEX REAL(FB), INTENT(IN):: STIME LOGICAL, INTENT(IN) :: CC_INTERP2FACES,CC_CELL_CENTERED @@ -5632,7 +5649,7 @@ SUBROUTINE DUMP_SLICE_GEOM_DATA(FUNIT_DATA,CC_INTERP2FACES,CC_CELL_CENTERED,SLIC CALL GET_GEOMVALS(CC_INTERP2FACES,CC_CELL_CENTERED,SLICETYPE,& I1,I2,J1,J2,K1,K2,NFACES,NFACES_CUTCELLS,VALS,& IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,& - PIPE_INDEX,PROP_INDEX,REAC_INDEX,MATL_INDEX,T,DT,NM,OPT_BNDF_INDEX) + PIPE_INDEX,PROP_INDEX,REAC_INDEX,T,DT,NM,OPT_BNDF_INDEX) ELSE NVALS = NVERTS ALLOCATE(VALS(MAX(NVERTS,NFACES))) @@ -5647,7 +5664,7 @@ SUBROUTINE DUMP_SLICE_GEOM_DATA(FUNIT_DATA,CC_INTERP2FACES,CC_CELL_CENTERED,SLIC CALL GET_GEOMVALS(CC_INTERP2FACES,CC_CELL_CENTERED,SLICETYPE,& I1,I2,J1,J2,K1,K2,NFACES,NFACES_CUTCELLS,VALS,& IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,& - PIPE_INDEX,PROP_INDEX,REAC_INDEX,MATL_INDEX,T,DT,NM,OPT_BNDF_INDEX) + PIPE_INDEX,PROP_INDEX,REAC_INDEX,T,DT,NM,OPT_BNDF_INDEX) ! these two routines need to be moved and called only once CALL GET_GEOMINFO(SLICETYPE,I1,I2,J1,J2,K1,K2,NVERTS,NVERTS_CUTCELLS,NFACES,NFACES_CUTCELLS,VERTS,FACES,LOCATIONS) @@ -5704,12 +5721,12 @@ END SUBROUTINE DUMP_SLICE_GEOM_DATA SUBROUTINE GET_SOLIDREGFACE_SCALAR_SLICE(X1AXIS,I,J,K,VAL_CF,& - IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX,MATL_INDEX,T,DT,NM) + IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX,T,DT,NM) USE PHYSICAL_FUNCTIONS, ONLY: GET_MASS_FRACTION REAL(EB), INTENT(IN) :: T,DT -INTEGER, INTENT(IN) :: X1AXIS,I,J,K,IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX,MATL_INDEX,NM +INTEGER, INTENT(IN) :: X1AXIS,I,J,K,IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX,NM REAL(EB),INTENT(OUT):: VAL_CF ! Local Variables: @@ -5750,13 +5767,13 @@ SUBROUTINE GET_SOLIDREGFACE_SCALAR_SLICE(X1AXIS,I,J,K,VAL_CF,& IF( CCSUM > 0._EB ) CC1(LOW_IND:HIGH_IND)=CC1(LOW_IND:HIGH_IND)/CCSUM IF (CC1( LOW_IND)>TWENTY_EPSILON_EB) THEN - VAL_CF_LO = GAS_PHASE_OUTPUT(T,DT,NM,II_LO,JJ_LO,KK_LO,& - IND,Y_INDEX,Z_INDEX,0,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX,MATL_INDEX) + VAL_CF_LO = GAS_PHASE_OUTPUT(T,DT,NM,II_LO,JJ_LO,KK_LO,IND,Y_INDEX=Y_INDEX,Z_INDEX=Z_INDEX,PART_INDEX=PART_INDEX,& + VELO_INDEX=VELO_INDEX,PIPE_INDEX=PIPE_INDEX,PROP_INDEX=PROP_INDEX,REAC_INDEX=REAC_INDEX) ENDIF IF (CC1(HIGH_IND)>TWENTY_EPSILON_EB) THEN - VAL_CF_HI = GAS_PHASE_OUTPUT(T,DT,NM,II_HI,JJ_HI,KK_HI,& - IND,Y_INDEX,Z_INDEX,0,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX,MATL_INDEX) + VAL_CF_HI = GAS_PHASE_OUTPUT(T,DT,NM,II_HI,JJ_HI,KK_HI,IND,Y_INDEX=Y_INDEX,Z_INDEX=Z_INDEX,PART_INDEX=PART_INDEX,& + VELO_INDEX=VELO_INDEX,PIPE_INDEX=PIPE_INDEX,PROP_INDEX=PROP_INDEX,REAC_INDEX=REAC_INDEX) ENDIF VAL_CF = CC1(LOW_IND)*VAL_CF_LO + CC1(HIGH_IND)*VAL_CF_HI @@ -5766,12 +5783,12 @@ END SUBROUTINE GET_SOLIDREGFACE_SCALAR_SLICE SUBROUTINE GET_SOLIDCUTFACE_SCALAR_SLICE(X1AXIS,ICF,VAL_CF, & - IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX,MATL_INDEX,T,DT,NM) + IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX,T,DT,NM) USE PHYSICAL_FUNCTIONS, ONLY: GET_MASS_FRACTION REAL(EB), INTENT(IN) :: T,DT -INTEGER, INTENT(IN) :: X1AXIS,ICF,IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX,MATL_INDEX,NM +INTEGER, INTENT(IN) :: X1AXIS,ICF,IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX,NM REAL(EB),INTENT(OUT):: VAL_CF ! Local Variables: @@ -5887,8 +5904,8 @@ SUBROUTINE GET_SOLIDCUTFACE_SCALAR_SLICE(X1AXIS,ICF,VAL_CF, & ! Use closest solid Cell values for SOLID cut-face: IF (FOUND) THEN - VAL_CF = GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,& - IND,Y_INDEX,Z_INDEX,0,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX,MATL_INDEX) + VAL_CF = GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,IND,Y_INDEX=Y_INDEX,Z_INDEX=Z_INDEX,PART_INDEX=PART_INDEX,VELO_INDEX=VELO_INDEX,& + PIPE_INDEX=PIPE_INDEX,PROP_INDEX=PROP_INDEX,REAC_INDEX=REAC_INDEX) ENDIF RETURN @@ -5896,12 +5913,12 @@ END SUBROUTINE GET_SOLIDCUTFACE_SCALAR_SLICE SUBROUTINE GET_GASCUTFACE_SCALAR_SLICE(VAL_CF,X1AXIS,ICF,IFACE,CC_INTERP2FACES,CC_CELL_CENTERED,& - IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX,MATL_INDEX,T,DT,NM) + IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX,T,DT,NM) USE PHYSICAL_FUNCTIONS, ONLY: GET_MASS_FRACTION REAL(EB), INTENT(IN) :: T,DT -INTEGER, INTENT(IN) :: X1AXIS,ICF,IFACE,IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX,MATL_INDEX,NM +INTEGER, INTENT(IN) :: X1AXIS,ICF,IFACE,IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX,NM LOGICAL, INTENT(IN) :: CC_INTERP2FACES,CC_CELL_CENTERED REAL(EB),INTENT(OUT):: VAL_CF @@ -5938,9 +5955,9 @@ SUBROUTINE GET_GASCUTFACE_SCALAR_SLICE(VAL_CF,X1AXIS,ICF,IFACE,CC_INTERP2FACES,C II = CUT_CELL(ICC)%IJK(IAXIS) JJ = CUT_CELL(ICC)%IJK(JAXIS) KK = CUT_CELL(ICC)%IJK(KAXIS) - VAL_LOC(ISIDE) = GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,& - IND,Y_INDEX,Z_INDEX,0,PART_INDEX,VELO_INDEX,& - PIPE_INDEX,PROP_INDEX,REAC_INDEX,MATL_INDEX,ICC,JCC) + VAL_LOC(ISIDE) = GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,IND,Y_INDEX=Y_INDEX,Z_INDEX=Z_INDEX,PART_INDEX=PART_INDEX,& + VELO_INDEX=VELO_INDEX,PIPE_INDEX=PIPE_INDEX,PROP_INDEX=PROP_INDEX,REAC_INDEX=REAC_INDEX,& + ICC_IN=ICC,JCC_IN=JCC) END SELECT ENDDO VAL_CF = CCM1*VAL_LOC(LOW_IND) + CCP1*VAL_LOC(HIGH_IND) @@ -6149,8 +6166,8 @@ SUBROUTINE DUMP_SLCF(T,DT,NM,IFRMT) DO K=KK1,KK2 DO J=JJ1,JJ2 DO I=II1,II2 - QUANTITY(I,J,K) = GAS_PHASE_OUTPUT(T,DT,NM,I,J,K,IND,Y_INDEX,Z_INDEX,0,PART_INDEX,VELO_INDEX,0,& - PROP_INDEX,REAC_INDEX,MATL_INDEX) + QUANTITY(I,J,K) = GAS_PHASE_OUTPUT(T,DT,NM,I,J,K,IND,Y_INDEX=Y_INDEX,Z_INDEX=Z_INDEX,PART_INDEX=PART_INDEX,& + VELO_INDEX=VELO_INDEX,PROP_INDEX=PROP_INDEX,REAC_INDEX=REAC_INDEX) ENDDO ENDDO ENDDO @@ -6159,7 +6176,8 @@ SUBROUTINE DUMP_SLCF(T,DT,NM,IFRMT) DO I=II1,II2 DO J=JJ1,JJ2 KTS = K_AGL_SLICE(I,J,NTSL) - QUANTITY(I,J,K1) = GAS_PHASE_OUTPUT(T,DT,NM,I,J,KTS,IND,Y_INDEX,Z_INDEX,0,PART_INDEX,VELO_INDEX,0,0,0,0) + QUANTITY(I,J,K1) = GAS_PHASE_OUTPUT(T,DT,NM,I,J,KTS,IND,Y_INDEX=Y_INDEX,Z_INDEX=Z_INDEX,PART_INDEX=PART_INDEX,& + VELO_INDEX=VELO_INDEX) ENDDO ENDDO ENDIF @@ -6333,7 +6351,7 @@ SUBROUTINE DUMP_SLCF(T,DT,NM,IFRMT) OPEN(LU_SLCF(IQ,NM),FILE=FN_SLCF(IQ,NM),FORM='UNFORMATTED',STATUS='REPLACE') CALL DUMP_SLICE_GEOM_DATA(LU_SLCF(IQ,NM),CC_INTERP2FACES,SL%CELL_CENTERED,SL%SLICETYPE, & 1,STIME,I1,I2,J1,J2,K1,K2,DEBUG,& - IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,0,PROP_INDEX,REAC_INDEX,MATL_INDEX,T,DT,NM, & + IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,0,PROP_INDEX,REAC_INDEX,T,DT,NM, & SLICE_MIN, SLICE_MAX) SLICE_MIN_BOUND = SLICE_MIN SLICE_MAX_BOUND = SLICE_MAX @@ -6343,7 +6361,7 @@ SUBROUTINE DUMP_SLCF(T,DT,NM,IFRMT) OPEN(LU_SLCF(IQ,NM),FILE=FN_SLCF(IQ,NM),FORM='UNFORMATTED',STATUS='OLD',POSITION='APPEND') CALL DUMP_SLICE_GEOM_DATA(LU_SLCF(IQ,NM),CC_INTERP2FACES,SL%CELL_CENTERED,SL%SLICETYPE, & 0,STIME,I1,I2,J1,J2,K1,K2,DEBUG,& - IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,0,PROP_INDEX,REAC_INDEX,MATL_INDEX,T,DT,NM, & + IND,Y_INDEX,Z_INDEX,PART_INDEX,VELO_INDEX,0,PROP_INDEX,REAC_INDEX,T,DT,NM, & SLICE_MIN, SLICE_MAX) OPEN(LU_SLCF(IQ2,NM),FILE=FN_SLCF(IQ2,NM),ACTION='READ') READ(LU_SLCF(IQ2,NM),FMT=*,IOSTAT=IERROR)T_BOUND, SLICE_MIN_BOUND, SLICE_MAX_BOUND @@ -6517,6 +6535,9 @@ SUBROUTINE UPDATE_DEVICES_1(T,DT,NM) INTEGER :: N,I,J,K,IW,ICC,ICF,SURF_INDEX,LP_INDEX,IP,AXIS TYPE (BOUNDARY_PROP1_TYPE), POINTER :: B1 TYPE (BOUNDARY_COORD_TYPE), POINTER :: BC +TYPE (PROPERTY_TYPE), POINTER :: PY +TYPE (CFACE_TYPE), POINTER :: CFA +TYPE (DEVICE_TYPE), POINTER :: DV ! If any device has QUANTITY='PARTICLE FLUX N', pre-compute PARTICLE fluxes @@ -6693,14 +6714,15 @@ SUBROUTINE UPDATE_DEVICES_1(T,DT,NM) I = MIN( IBP1, MAX(0, DV%I(1)) ) J = MIN( JBP1, MAX(0, DV%J(1)) ) K = MIN( KBP1, MAX(0, DV%K(1)) ) - SDV%VALUE_1 = GAS_PHASE_OUTPUT(T,DT,NM,I,J,K,DV%QUANTITY_INDEX(1),DV%Y_INDEX,DV%Z_INDEX,DV%ELEM_INDEX,& - DV%PART_CLASS_INDEX,DV%VELO_INDEX,DV%PIPE_INDEX,DV%PROP_INDEX,DV%REAC_INDEX,& - DV%MATL_INDEX) + SDV%VALUE_1 = GAS_PHASE_OUTPUT(T,DT,NM,I,J,K,DV%QUANTITY_INDEX(1),Y_INDEX=DV%Y_INDEX,Z_INDEX=DV%Z_INDEX,& + ELEM_INDX=DV%ELEM_INDEX,PART_INDEX=DV%PART_CLASS_INDEX,VELO_INDEX=DV%VELO_INDEX,& + PIPE_INDEX=DV%PIPE_INDEX,PROP_INDEX=DV%PROP_INDEX,REAC_INDEX=DV%REAC_INDEX) IF (DV%N_QUANTITY>1) & - SDV%VALUE_2 = GAS_PHASE_OUTPUT(T,DT,NM,DV%I(2),DV%J(2),DV%K(2),DV%QUANTITY_INDEX(2),DV%Y_INDEX,DV%Z_INDEX,& - DV%ELEM_INDEX,DV%PART_CLASS_INDEX,DV%VELO_INDEX,DV%PIPE_INDEX,DV%PROP_INDEX,& - DV%REAC_INDEX,DV%MATL_INDEX) + SDV%VALUE_2 = GAS_PHASE_OUTPUT(T,DT,NM,DV%I(2),DV%J(2),DV%K(2),DV%QUANTITY_INDEX(2),Y_INDEX=DV%Y_INDEX,& + Z_INDEX=DV%Z_INDEX,ELEM_INDX=DV%ELEM_INDEX,PART_INDEX=DV%PART_CLASS_INDEX,& + VELO_INDEX=DV%VELO_INDEX,PIPE_INDEX=DV%PIPE_INDEX,PROP_INDEX=DV%PROP_INDEX,& + REAC_INDEX=DV%REAC_INDEX) CASE DEFAULT GAS_STATS_SELECT @@ -6743,9 +6765,9 @@ SUBROUTINE UPDATE_DEVICES_1(T,DT,NM) VOL = VOL/DY(J) ENDIF - VALUE = GAS_PHASE_OUTPUT(T,DT,NM,I,J,K,DV%QUANTITY_INDEX(1),DV%Y_INDEX,DV%Z_INDEX,DV%ELEM_INDEX,& - DV%PART_CLASS_INDEX,DV%VELO_INDEX,DV%PIPE_INDEX,DV%PROP_INDEX,DV%REAC_INDEX,& - DV%MATL_INDEX) + VALUE = GAS_PHASE_OUTPUT(T,DT,NM,I,J,K,DV%QUANTITY_INDEX(1),Y_INDEX=DV%Y_INDEX,Z_INDEX=DV%Z_INDEX,& + ELEM_INDX=DV%ELEM_INDEX,PART_INDEX=DV%PART_CLASS_INDEX,VELO_INDEX=DV%VELO_INDEX,& + PIPE_INDEX=DV%PIPE_INDEX,PROP_INDEX=DV%PROP_INDEX,REAC_INDEX=DV%REAC_INDEX) STATISTICS_SELECT: SELECT CASE(DV%SPATIAL_STATISTIC) CASE('MAX','MAXLOC X','MAXLOC Y','MAXLOC Z') IF (VALUE>SDV%VALUE_1) THEN @@ -7056,7 +7078,8 @@ SUBROUTINE UPDATE_DEVICES_2(T,DT) REAL(EB) :: WGT,WGT_UNBIASED INTEGER :: N,IERR,INTERVAL_INDEX REAL(EB) :: Z_INT_DENOM - +TYPE(PROPERTY_TYPE), POINTER :: PY +TYPE(DEVICE_TYPE), POINTER :: DV DEVICE_LOOP: DO N=1,N_DEVC @@ -7201,7 +7224,7 @@ SUBROUTINE UPDATE_DEVICES_2(T,DT) DV%VALUE = DV%PREVIOUS_VALUE ENDIF ELSEIF (T+DT>DV%STATISTICS_END .AND. DV%TIME_PERIOD>0._EB) THEN - CALL EXTRAPOLATE_EXTREMA + CALL EXTRAPOLATE_EXTREMA(DV) ELSE DV%VALUE = DV%TIME_MAX_VALUE(INTERVAL_INDEX) ENDIF @@ -7219,7 +7242,7 @@ SUBROUTINE UPDATE_DEVICES_2(T,DT) DV%VALUE = DV%PREVIOUS_VALUE ENDIF ELSEIF (T+DT>DV%STATISTICS_END .AND. DV%TIME_PERIOD>0._EB) THEN - CALL EXTRAPOLATE_EXTREMA + CALL EXTRAPOLATE_EXTREMA(DV) ELSE DV%VALUE = DV%TIME_MIN_VALUE(INTERVAL_INDEX) ENDIF @@ -7307,11 +7330,10 @@ END SUBROUTINE UPDATE_DEVICES_2 !> \param PIPE_INDEX Index of the pipe branch !> \param PROP_INDEX Index of the PROPerty group parameters !> \param REAC_INDEX Index of the REACtion -!> \param MATL_INDEX Index of the Material !> \param ICC_IN,JCC_IN Optional indexes of cut-cell. REAL(EB) RECURSIVE FUNCTION GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,IND,Y_INDEX,Z_INDEX,ELEM_INDX,PART_INDEX,VELO_INDEX,PIPE_INDEX,& - PROP_INDEX,REAC_INDEX,MATL_INDEX,ICC_IN,JCC_IN) RESULT(GAS_PHASE_OUTPUT_RES) + PROP_INDEX,REAC_INDEX,ICC_IN,JCC_IN) RESULT(GAS_PHASE_OUTPUT_RES) USE MEMORY_FUNCTIONS, ONLY: REALLOCATE USE MATH_FUNCTIONS, ONLY: INTERPOLATE1D,INTERPOLATE1D_UNIFORM,UPDATE_HISTOGRAM @@ -7327,9 +7349,8 @@ REAL(EB) RECURSIVE FUNCTION GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,IND,Y_INDEX,Z_INDE USE CC_SCALARS, ONLY: CC_CUTCELL_VELOCITY REAL(EB), INTENT(IN) :: T,DT -INTEGER, INTENT(IN) :: II,JJ,KK,IND,NM,Y_INDEX,Z_INDEX,ELEM_INDX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX, & - MATL_INDEX -INTEGER, INTENT(IN), OPTIONAL :: ICC_IN,JCC_IN +INTEGER, INTENT(IN) :: II,JJ,KK,IND,NM +INTEGER, INTENT(IN), OPTIONAL :: Y_INDEX,Z_INDEX,ELEM_INDX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX,ICC_IN,JCC_IN REAL(EB) :: H_TC,TMP_TC,RE_D,NUSSELT,VEL,K_G,MU_G,COSTHETA,FAC,& Q_SUM,TMP_G,UU,VV,WW,VEL2,Y_MF_INT,PATHLENGTH,EXT_COEF,MASS_EXT_COEF,ZZ_FUEL,ZZ_OX,& VELSR,WATER_VOL_FRAC,WATER_TEMP,RHS,DT_C,DT_E,T_RATIO,Y_E_LAG, H_G,H_G_SUM,CPBAR,CP,ZZ_GET(1:N_TRACKED_SPECIES),RCON,& @@ -7346,6 +7367,8 @@ REAL(EB) RECURSIVE FUNCTION GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,IND,Y_INDEX,Z_INDE TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1 TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC TYPE(BOUNDARY_ONE_D_TYPE), POINTER :: ONE_D +TYPE(PROPERTY_TYPE), POINTER :: PY +TYPE(DEVICE_TYPE), POINTER :: DV,DV2 ! Get species mass fraction if necessary @@ -7367,6 +7390,11 @@ REAL(EB) RECURSIVE FUNCTION GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,IND,Y_INDEX,Z_INDE R_Y_H2O = SPECIES(H2O_INDEX)%RCON * Y_H2O IF (Y_INDEX==H2O_INDEX) Y_SPECIES=0._EB ENDIF +IF (PRESENT(PROP_INDEX)) THEN + PY => PROPERTY(PROP_INDEX) +ELSE + PY => PROPERTY(0) +ENDIF ! Get desired output value @@ -7745,8 +7773,7 @@ REAL(EB) RECURSIVE FUNCTION GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,IND,Y_INDEX,Z_INDE S23 = 0.5_EB*(DVDZ+DWDY) GAS_PHASE_OUTPUT_RES = SQRT(2._EB*(S11**2 + S22**2 + S33**2 + 2._EB*(S12**2 + S13**2 + S23**2))) CASE(85) ! KOLMOGOROV LENGTH SCALE - SS = GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,84,Y_INDEX,Z_INDEX,ELEM_INDX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,& - REAC_INDEX,MATL_INDEX) + SS = GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,84) DISSIPATION_RATE = MU(II,JJ,KK)/RHO(II,JJ,KK)*SS**2 GAS_PHASE_OUTPUT_RES = ((MU_DNS(II,JJ,KK)/RHO(II,JJ,KK))**3/(DISSIPATION_RATE+EPS))**0.25_EB CASE(86) ! CELL REYNOLDS NUMBER @@ -7754,14 +7781,12 @@ REAL(EB) RECURSIVE FUNCTION GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,IND,Y_INDEX,Z_INDE JJJ = MAX(1,MIN(JJ,JBAR)) KKK = MAX(1,MIN(KK,KBAR)) DELTA = LES_FILTER_WIDTH(III,JJJ,KKK) - ETA = GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,85,Y_INDEX,Z_INDEX,ELEM_INDX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,& - REAC_INDEX,MATL_INDEX) + ETA = GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,85) GAS_PHASE_OUTPUT_RES = DELTA/(ETA+EPS) CASE(87) ! MOLECULAR VISCOSITY GAS_PHASE_OUTPUT_RES = MU_DNS(II,JJ,KK) CASE(88) ! DISSIPATION RATE - SS = GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,84,Y_INDEX,Z_INDEX,ELEM_INDX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,& - REAC_INDEX,MATL_INDEX) + SS = GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,84) GAS_PHASE_OUTPUT_RES = MU(II,JJ,KK)/RHO(II,JJ,KK)*SS**2 CASE(89) ! KINEMATIC VISCOSITY GAS_PHASE_OUTPUT_RES = MU(II,JJ,KK)/RHO(II,JJ,KK) @@ -7895,8 +7920,7 @@ REAL(EB) RECURSIVE FUNCTION GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,IND,Y_INDEX,Z_INDE ! McCaffrey and Heskestad, A Robust Bidirectional Low-Velocity Probe for Flame and Fire Application ! Combustion and Flame, 26, 125 - 127, (1976). IF (PY%TC) THEN - PROBE_TMP = GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,110,Y_INDEX,Z_INDEX,ELEM_INDX,PART_INDEX,VELO_INDEX,PIPE_INDEX,& - PROP_INDEX,REAC_INDEX,MATL_INDEX,ICC_IN,JCC_IN) + TMPM + PROBE_TMP = GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,110,PROP_INDEX=PROP_INDEX) + TMPM ELSE PROBE_TMP = TMP(II,JJ,KK) ENDIF @@ -8769,13 +8793,14 @@ END FUNCTION GAS_PHASE_OUTPUT !> \param OPT_PROF_INDEX Index of PROFile REAL(EB) FUNCTION SOLID_PHASE_OUTPUT(INDX,T,NM,Y_INDEX,Z_INDEX,PART_INDEX,OPT_WALL_INDEX,OPT_LP_INDEX,OPT_BNDF_INDEX,& - OPT_DEVC_INDEX,OPT_CFACE_INDEX,OPT_CUT_FACE_INDEX,OPT_NODE_INDEX,OPT_PROF_INDEX) + OPT_DEVC_INDEX,OPT_CFACE_INDEX,OPT_CUT_FACE_INDEX,OPT_NODE_INDEX,OPT_PROF_INDEX,& + OPT_PROP_INDEX) USE PHYSICAL_FUNCTIONS, ONLY: SURFACE_DENSITY,GET_MASS_FRACTION,GET_SENSIBLE_ENTHALPY,GET_SPECIFIC_HEAT,GET_CONDUCTIVITY,& GET_VISCOSITY,HEAT_TRANSFER_COEFFICIENT USE TURBULENCE, ONLY: TAU_WALL_IJ INTEGER, INTENT(IN), OPTIONAL :: OPT_WALL_INDEX,OPT_LP_INDEX,OPT_CFACE_INDEX,OPT_BNDF_INDEX,OPT_DEVC_INDEX,OPT_CUT_FACE_INDEX,& - OPT_NODE_INDEX,OPT_PROF_INDEX + OPT_NODE_INDEX,OPT_PROF_INDEX,OPT_PROP_INDEX INTEGER, INTENT(IN) :: INDX,Y_INDEX,Z_INDEX,PART_INDEX,NM REAL(EB),INTENT(IN) :: T REAL(EB) :: Q_CON,RHOSUM,VOLSUM,ZZ_GET(1:N_TRACKED_SPECIES),Y_SPECIES,DEPTH,ASH_DEPTH,UN,H_S,RHO_D_DYDN,U_CELL,V_CELL,W_CELL,& @@ -8788,6 +8813,9 @@ REAL(EB) FUNCTION SOLID_PHASE_OUTPUT(INDX,T,NM,Y_INDEX,Z_INDEX,PART_INDEX,OPT_WA TYPE(BOUNDARY_RADIA_TYPE), POINTER :: BR TYPE(BOUNDARY_ONE_D_TYPE), POINTER :: ONE_D TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC +TYPE(PROPERTY_TYPE), POINTER :: PY +TYPE(CFACE_TYPE), POINTER :: CFA +TYPE(DEVICE_TYPE), POINTER :: DV ! Assign default value to output @@ -8795,6 +8823,12 @@ REAL(EB) FUNCTION SOLID_PHASE_OUTPUT(INDX,T,NM,Y_INDEX,Z_INDEX,PART_INDEX,OPT_WA IF (PRESENT(OPT_DEVC_INDEX)) DV => DEVICE(OPT_DEVC_INDEX) +IF (PRESENT(OPT_PROP_INDEX)) THEN + PY => PROPERTY(OPT_PROP_INDEX) +ELSE + PY => PROPERTY(0) +ENDIF + IF (PRESENT(OPT_WALL_INDEX)) THEN IF (OPT_WALL_INDEX==0) RETURN @@ -9825,6 +9859,8 @@ SUBROUTINE DUMP_DEVICES(T) REAL(EB) :: STIME,DI,DD,VALUE INTEGER :: I,J,N,NN,N_OUT REAL(EB) :: DEVC_TIME,CONST,CUMSUM,COORD_FACTOR +TYPE(PROPERTY_TYPE), POINTER :: PY +TYPE(DEVICE_TYPE), POINTER :: DV ! Determine the time to write into file @@ -10282,6 +10318,7 @@ SUBROUTINE UPDATE_HRR(DT,NM) INTEGER :: I,J,K,IW,IIG,JJG,KKG,N,IND1,IND2,ICF,ICC,JCC,IP TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1 TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC +TYPE(CFACE_TYPE), POINTER :: CFA ! Compute volume integral of certain quantities like the HRR @@ -10592,6 +10629,7 @@ SUBROUTINE DUMP_BNDF(T,DT,NM) TYPE(PATCH_TYPE), POINTER :: PA REAL(FB) BNDF_TIME, BNDF_VAL_MIN, BNDF_VAL_MAX INTEGER :: CHANGE_BOUND, IERROR +TYPE(BOUNDARY_FILE_TYPE), POINTER :: BF IF (MESHES(NM)%N_PATCH==0 .AND. MESHES(NM)%N_INTERNAL_CFACE_CELLS==0) RETURN IF (.NOT. MESHES(NM)%BNDF_DUMP) RETURN @@ -10603,7 +10641,6 @@ SUBROUTINE DUMP_BNDF(T,DT,NM) FILE_LOOP: DO NF=1,N_BNDF IF (N_PATCH == 0) CYCLE FILE_LOOP BF => BOUNDARY_FILE(NF) - PY => PROPERTY(BF%PROP_INDEX) BOUND_MAX = -1.0E+33_FB BOUND_MIN = -BOUND_MAX OPEN(LU_BNDF(NF,NM),FILE=FN_BNDF(NF,NM),FORM='UNFORMATTED',STATUS='OLD',POSITION='APPEND') @@ -10643,7 +10680,7 @@ SUBROUTINE DUMP_BNDF(T,DT,NM) WALL(IW)%BOUNDARY_TYPE/=INTERPOLATED_BOUNDARY .AND. .NOT.CELL(IC)%SOLID) THEN IBK(L,N) = 1 PP(L,N) = REAL(SOLID_PHASE_OUTPUT(IND,T,NM,BF%Y_INDEX,BF%Z_INDEX,BF%PART_INDEX,OPT_WALL_INDEX=IW,& - OPT_BNDF_INDEX=NF),FB) + OPT_BNDF_INDEX=NF,OPT_PROP_INDEX=BF%PROP_INDEX),FB) ENDIF ENDDO ENDDO @@ -10682,7 +10719,7 @@ SUBROUTINE DUMP_BNDF(T,DT,NM) PPN(L,N) = PPN(L,N)/REAL(ISUM,FB) ELSE PPN(L,N) = REAL(SOLID_PHASE_OUTPUT(IND,T,NM,BF%Y_INDEX,BF%Z_INDEX,BF%PART_INDEX,OPT_WALL_INDEX=0,& - OPT_BNDF_INDEX=NF),FB) + OPT_BNDF_INDEX=NF,OPT_PROP_INDEX=BF%PROP_INDEX),FB) ENDIF ENDDO ENDDO @@ -10756,7 +10793,6 @@ SUBROUTINE DUMP_BNDF(T,DT,NM) IF (CC_IBM) THEN FILE_LOOP2 : DO NF=1,N_BNDF BF => BOUNDARY_FILE(NF) - PY => PROPERTY(BF%PROP_INDEX) IND = ABS(BF%INDEX) NC = 0 I1=0; I2=-1; J1=0; J2=-1; K1=0; K2=-1; ! Just dummy numbers, not needed for INBOUND_FACES @@ -10765,7 +10801,7 @@ SUBROUTINE DUMP_BNDF(T,DT,NM) IF (REAL(T-T_BEGIN,FB)<2._FB*EPSILON(REAL(T,FB))) THEN OPEN(LU_BNDG(NF,NM), FILE=FN_BNDG(NF,NM), FORM='UNFORMATTED',STATUS='REPLACE') CALL DUMP_SLICE_GEOM_DATA(LU_BNDG(NF,NM),.FALSE.,.TRUE.,"INBOUND_FACES",1,STIME,I1,I2,J1,J2,K1,K2,BF%DEBUG, & - IND,BF%Y_INDEX,BF%Z_INDEX,BF%PART_INDEX,0,0,BF%PROP_INDEX,0,0,T,DT,NM, & + IND,BF%Y_INDEX,BF%Z_INDEX,BF%PART_INDEX,0,0,BF%PROP_INDEX,0,T,DT,NM, & BOUND_MIN, BOUND_MAX, OPT_BNDF_INDEX=NF) BNDF_VAL_MIN = BOUND_MIN BNDF_VAL_MAX = BOUND_MAX @@ -10774,7 +10810,7 @@ SUBROUTINE DUMP_BNDF(T,DT,NM) ! data file at subsequent time steps OPEN(LU_BNDG(NF,NM), FILE=FN_BNDG(NF,NM), FORM='UNFORMATTED',STATUS='OLD',POSITION='APPEND') CALL DUMP_SLICE_GEOM_DATA(LU_BNDG(NF,NM),.FALSE.,.TRUE.,"INBOUND_FACES",0,STIME,I1,I2,J1,J2,K1,K2,BF%DEBUG, & - IND,BF%Y_INDEX,BF%Z_INDEX,BF%PART_INDEX,0,0,BF%PROP_INDEX,0,0,T,DT,NM, & + IND,BF%Y_INDEX,BF%Z_INDEX,BF%PART_INDEX,0,0,BF%PROP_INDEX,0,T,DT,NM, & BOUND_MIN, BOUND_MAX, OPT_BNDF_INDEX=NF) OPEN(LU_BNDG(NF+N_BNDF,NM),FILE=FN_BNDG(NF+N_BNDF,NM), ACTION='READ') READ(LU_BNDG(NF+N_BNDF,NM),FMT=*,IOSTAT=IERROR)BNDF_TIME, BNDF_VAL_MIN, BNDF_VAL_MAX @@ -10898,6 +10934,7 @@ SUBROUTINE TIMINGS INTEGER :: N LOGICAL :: WRITE_HEADER TYPE(CONTROL_TYPE), POINTER :: CF +TYPE(DEVICE_TYPE), POINTER :: DV ! Print out detector and control activation times @@ -11410,11 +11447,12 @@ END SUBROUTINE DUMP_MMS !> \details Given DV%N_INTERVALS values of DV%TIME_MIN(MAX)_VALUE, extrapolate the MIN(MAX) for the DV%TIME_PERIOD, !> which is typically longer than the statistical sampling duration DV%STATISTICS_END-DV%STATISTICS_START -SUBROUTINE EXTRAPOLATE_EXTREMA +SUBROUTINE EXTRAPOLATE_EXTREMA(DV) REAL(EB), ALLOCATABLE, DIMENSION(:) :: EXTREMA,YYY INTEGER :: I,J REAL(EB) :: A,B,X_AVG,Y_AVG,TT,ST2,INTERVAL_RATIO,F +TYPE(DEVICE_TYPE), POINTER :: DV ALLOCATE(EXTREMA(DV%N_INTERVALS)) ALLOCATE(YYY(DV%N_INTERVALS)) From 67cb690285adb0d2f95a73b6d0b6ffecc47e0d9a Mon Sep 17 00:00:00 2001 From: mcgratta Date: Fri, 8 May 2026 10:09:06 -0400 Subject: [PATCH 2/2] FDS Source: Fix a few bugs related to recursive --- Source/dump.f90 | 145 ++++++++++++++++++++++++++++++++++-------------- 1 file changed, 102 insertions(+), 43 deletions(-) diff --git a/Source/dump.f90 b/Source/dump.f90 index 175cfa109a..13dbbfb95d 100644 --- a/Source/dump.f90 +++ b/Source/dump.f90 @@ -14,33 +14,21 @@ MODULE DUMP USE CONTROL_VARIABLES USE OUTPUT_DATA USE PROPERTY_DATA -USE CHEMCONS, ONLY : WRITE_CVODE_SUBSTEPS, CVODE_SUBSTEP_DATA, TOTAL_SUBSTEPS_TAKEN -USE COMPLEX_GEOMETRY, ONLY : WRITE_GEOM,WRITE_GEOM_ALL,CC_FGSC,CC_IDCF,CC_IDCC,CC_UNKZ,CC_UNKF,CC_FTYPE_RCGAS,& - CC_FTYPE_CFGAS,CC_FTYPE_CFINB,CC_SOLID,CC_CGSC,CC_IDRC,CC_CUTCFE,TRIANGULATE,& - CC_VGSC,CC_GASPHASE,MAKE_UNIQUE_VERT_ARRAY,AVERAGE_FACE_VALUES +USE CHEMCONS, ONLY: WRITE_CVODE_SUBSTEPS, CVODE_SUBSTEP_DATA, TOTAL_SUBSTEPS_TAKEN +USE COMPLEX_GEOMETRY, ONLY: WRITE_GEOM,WRITE_GEOM_ALL,CC_FGSC,CC_IDCF,CC_IDCC,CC_UNKZ,CC_UNKF,CC_FTYPE_RCGAS,& + CC_FTYPE_CFGAS,CC_FTYPE_CFINB,CC_SOLID,CC_CGSC,CC_IDRC,CC_CUTCFE,TRIANGULATE,& + CC_VGSC,CC_GASPHASE,MAKE_UNIQUE_VERT_ARRAY,AVERAGE_FACE_VALUES + +USE CC_SCALARS, ONLY: ADD_Q_DOT_CUTCELLS,GET_PRES_CFACE,GET_PRES_CFACE_TEST,GET_UVWGAS_CFACE,GET_MUDNS_CFACE +USE COMP_FUNCTIONS, ONLY: SHUTDOWN -USE CC_SCALARS, ONLY : ADD_Q_DOT_CUTCELLS,GET_PRES_CFACE,GET_PRES_CFACE_TEST,GET_UVWGAS_CFACE,GET_MUDNS_CFACE -USE COMP_FUNCTIONS, ONLY : SHUTDOWN IMPLICIT NONE (TYPE,EXTERNAL) PRIVATE REAL(EB), POINTER, DIMENSION(:,:,:) :: WFX,WFY,WFZ INTEGER :: N_DEVC_FILES CHARACTER(80) :: TCFORM -LOGICAL :: EX,DRY,OPN,FROM_BNDF=.FALSE. - -TYPE (MESH_TYPE), POINTER :: M -TYPE (LAGRANGIAN_PARTICLE_TYPE), POINTER :: LP -TYPE (VENTS_TYPE), POINTER :: VT -TYPE (LAGRANGIAN_PARTICLE_CLASS_TYPE), POINTER :: LPC -TYPE (SPECIES_TYPE), POINTER :: SS -TYPE (REACTION_TYPE), POINTER :: RN -TYPE (SURFACE_TYPE),POINTER :: SF -TYPE (MATERIAL_TYPE),POINTER :: ML -TYPE (SUBDEVICE_TYPE), POINTER :: SDV -TYPE (SLICE_TYPE), POINTER :: SL -TYPE (WALL_TYPE), POINTER :: WC -TYPE (THIN_WALL_TYPE), POINTER :: TW +LOGICAL :: DRY,FROM_BNDF=.FALSE. PUBLIC ASSIGN_FILE_NAMES,INITIALIZE_GLOBAL_DUMPS,INITIALIZE_MESH_DUMPS,WRITE_STATUS_FILES, & TIMINGS,FLUSH_GLOBAL_BUFFERS,READ_RESTART,WRITE_DIAGNOSTICS, & @@ -217,6 +205,8 @@ SUBROUTINE ASSIGN_FILE_NAMES USE COMP_FUNCTIONS, ONLY: GET_FILE_NUMBER INTEGER :: NM,I,N CHARACTER(LABEL_LENGTH) :: CFORM +TYPE(MESH_TYPE), POINTER :: M +LOGICAL :: EX ! Set up file number counter @@ -532,6 +522,7 @@ SUBROUTINE INITIALIZE_GLOBAL_DUMPS(T,DT) USE HVAC_ROUTINES, ONLY: N_DUCT_QUANTITY,N_NODE_QUANTITY REAL(EB) :: TNOW REAL(EB), INTENT(IN) :: T,DT +LOGICAL :: EX INTEGER :: NN,I,N,N_OUT,N_ZONE_TMP,LU,J, N_NODE_OUT, N_DUCT_OUT, NS INTEGER, ALLOCATABLE, DIMENSION(:) :: DUCT_CELL CHARACTER(80) :: FN @@ -539,6 +530,7 @@ SUBROUTINE INITIALIZE_GLOBAL_DUMPS(T,DT) CHARACTER(LABEL_LENGTH+7), DIMENSION(42) :: LABEL='null' TYPE (PROPERTY_TYPE), POINTER :: PY TYPE (DEVICE_TYPE), POINTER :: DV +TYPE (LAGRANGIAN_PARTICLE_CLASS_TYPE), POINTER :: LPC TNOW=CURRENT_TIME() @@ -884,6 +876,10 @@ SUBROUTINE INITIALIZE_MESH_DUMPS(NM) TYPE (OBSTRUCTION_TYPE), POINTER :: OB TYPE (ISOSURFACE_FILE_TYPE), POINTER :: IS TYPE (BOUNDARY_FILE_TYPE), POINTER :: BF +TYPE (WALL_TYPE), POINTER :: WC +TYPE (SLICE_TYPE), POINTER :: SL +TYPE (LAGRANGIAN_PARTICLE_CLASS_TYPE), POINTER :: LPC +TYPE (MESH_TYPE), POINTER :: M TNOW=CURRENT_TIME() @@ -1530,6 +1526,11 @@ SUBROUTINE WRITE_SMOKEVIEW_FILE TYPE (OBSTRUCTION_TYPE), POINTER :: OB TYPE (BOUNDARY_FILE_TYPE), POINTER :: BF TYPE (DEVICE_TYPE), POINTER :: DV +TYPE (MATERIAL_TYPE),POINTER :: ML +TYPE (SURFACE_TYPE),POINTER :: SF +TYPE (LAGRANGIAN_PARTICLE_CLASS_TYPE), POINTER :: LPC +TYPE (VENTS_TYPE), POINTER :: VT +TYPE (MESH_TYPE), POINTER :: M ! If this is a RESTART case but an old .smv file does not exist, shutdown with an ERROR. @@ -2717,11 +2718,19 @@ SUBROUTINE INITIALIZE_DIAGNOSTIC_FILE(DT) REAL(EB) :: ZZ_GET(1:N_TRACKED_SPECIES),ZZ_REAC(1:N_TRACKED_SPECIES),ZZ_PROD(1:N_TRACKED_SPECIES),& MU_Z,K_Z,CP_ZN,H_Z, PHI_TILDE,TMP_FLAME CHARACTER(LABEL_LENGTH) :: QUANTITY,ODE_SOLVER,OUTFORM +LOGICAL :: EX TYPE(SPECIES_MIXTURE_TYPE),POINTER :: SM TYPE(PROPERTY_TYPE), POINTER :: PY TYPE(ISOSURFACE_FILE_TYPE), POINTER :: IS TYPE(BOUNDARY_FILE_TYPE), POINTER :: BF TYPE(DEVICE_TYPE), POINTER :: DV +TYPE(SLICE_TYPE), POINTER :: SL +TYPE(MATERIAL_TYPE),POINTER :: ML +TYPE(SURFACE_TYPE),POINTER :: SF +TYPE(REACTION_TYPE), POINTER :: RN +TYPE(SPECIES_TYPE), POINTER :: SS +TYPE(LAGRANGIAN_PARTICLE_CLASS_TYPE), POINTER :: LPC +TYPE(MESH_TYPE), POINTER :: M ! Open and initialize diagnostic output file @@ -3539,6 +3548,11 @@ SUBROUTINE DUMP_RESTART(T,DT,NM) TYPE(INITIALIZATION_TYPE), POINTER :: IN TYPE(CFACE_TYPE), POINTER :: CFA TYPE(DEVICE_TYPE), POINTER :: DV +TYPE(THIN_WALL_TYPE), POINTER :: TW +TYPE(WALL_TYPE), POINTER :: WC +TYPE(SURFACE_TYPE),POINTER :: SF +TYPE(LAGRANGIAN_PARTICLE_CLASS_TYPE), POINTER :: LPC +TYPE(LAGRANGIAN_PARTICLE_TYPE), POINTER :: LP OPEN(LU_CORE(NM),FILE=FN_CORE(NM),FORM='UNFORMATTED',STATUS='REPLACE') @@ -3723,6 +3737,11 @@ SUBROUTINE READ_RESTART(T,DT,NM) TYPE(INITIALIZATION_TYPE), POINTER :: IN TYPE(CFACE_TYPE), POINTER :: CFA TYPE(DEVICE_TYPE), POINTER :: DV +TYPE(THIN_WALL_TYPE), POINTER :: TW +TYPE(WALL_TYPE), POINTER :: WC +TYPE(SURFACE_TYPE),POINTER :: SF +TYPE(LAGRANGIAN_PARTICLE_CLASS_TYPE), POINTER :: LPC +TYPE(LAGRANGIAN_PARTICLE_TYPE), POINTER :: LP INQUIRE(FILE=FN_RESTART(NM),EXIST=EX) IF (.NOT.EX) THEN @@ -3966,6 +3985,7 @@ SUBROUTINE WRITE_DIAGNOSTICS(T,DT) CHARACTER(120) :: SIMPLE_OUTPUT,SIMPLE_OUTPUT_ERR,OUT_FORMAT CHARACTER(LABEL_LENGTH) :: DATE REAL(EB) :: TNOW,CPUTIME,STIME,DTS,MAX_CFL,MAX_VN,TROUND +TYPE(MESH_TYPE), POINTER :: M TNOW = CURRENT_TIME() @@ -4152,7 +4172,9 @@ SUBROUTINE DUMP_PART(T,NM) REAL(EB) :: PART_MIN, PART_MAX, PFACTOR REAL(FB) :: PFACTOR_FB INTEGER, PARAMETER :: PART_BOUNDFILE_VERSION=1 -TYPE (BOUNDARY_COORD_TYPE), POINTER :: BC +TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC +TYPE(LAGRANGIAN_PARTICLE_CLASS_TYPE), POINTER :: LPC +TYPE(LAGRANGIAN_PARTICLE_TYPE), POINTER :: LP ! Write the current time to the prt5 file, then start looping through the particle classes @@ -5732,10 +5754,9 @@ SUBROUTINE GET_SOLIDREGFACE_SCALAR_SLICE(X1AXIS,I,J,K,VAL_CF,& ! Local Variables: INTEGER :: II_LO,II_HI,JJ_LO,JJ_HI,KK_LO,KK_HI,SOLID_LO,SOLID_HI REAL(EB):: CC1(LOW_IND:HIGH_IND),CCSUM -REAL(EB) :: Y_SPECIES,VAL_CF_LO,VAL_CF_HI +REAL(EB) :: VAL_CF_LO,VAL_CF_HI VAL_CF = 0._EB -Y_SPECIES = 1._EB VAL_CF_LO = 0._EB VAL_CF_HI = 0._EB @@ -5794,13 +5815,11 @@ SUBROUTINE GET_SOLIDCUTFACE_SCALAR_SLICE(X1AXIS,ICF,VAL_CF, & ! Local Variables: INTEGER :: II_LO,II_HI,JJ_LO,JJ_HI,KK_LO,KK_HI,IJK(IAXIS:KAXIS),IJK2(IAXIS:KAXIS,16),ICELL,II,JJ,KK LOGICAL :: FOUND -REAL(EB):: Y_SPECIES ! Point to mesh has been called for MESHES(NM): This routine searches for a REGULAR SOLID cell in the ! vicinity of the SOLID cut-face and assigns to the latter the scalar value of the former. VAL_CF = 0._EB -Y_SPECIES = 1._EB IJK(IAXIS:KAXIS)=CUT_FACE(ICF)%IJK(IAXIS:KAXIS) @@ -5925,13 +5944,8 @@ SUBROUTINE GET_GASCUTFACE_SCALAR_SLICE(VAL_CF,X1AXIS,ICF,IFACE,CC_INTERP2FACES,C ! Local Variables: REAL(EB) :: X1F, IDX, CCM1, CCP1, VAL_LOC(LOW_IND:HIGH_IND) INTEGER :: ISIDE, ICC, JCC, LOCAL_IND, II, JJ, KK -REAL(EB) :: Y_SPECIES(LOW_IND:HIGH_IND) ! REAL(EB) :: ZZ_GET(1:N_TRACKED_SPECIES) -! Point to mesh has been called for MESHES(NM): - -Y_SPECIES(LOW_IND:HIGH_IND) = 1._EB - ! Here interpolate values from cut-cell centers: X1F= CUT_FACE(ICF)%XYZCEN(X1AXIS,IFACE) IDX= 1._EB/ ( CUT_FACE(ICF)%XCENHIGH(X1AXIS,IFACE) - & @@ -5993,6 +6007,7 @@ SUBROUTINE DUMP_SLCF(T,DT,NM,IFRMT) LOGICAL :: AGL_TERRAIN_SLICE,CC_CELL_CENTERED,CC_INTERP2FACES REAL(FB), ALLOCATABLE, DIMENSION(:) :: QQ_PACK TYPE (MESH_TYPE), POINTER :: M2 +TYPE (SLICE_TYPE), POINTER :: SL ! Return if there are no slices to process and this is not a Plot3D dump @@ -6538,6 +6553,10 @@ SUBROUTINE UPDATE_DEVICES_1(T,DT,NM) TYPE (PROPERTY_TYPE), POINTER :: PY TYPE (CFACE_TYPE), POINTER :: CFA TYPE (DEVICE_TYPE), POINTER :: DV +TYPE (WALL_TYPE), POINTER :: WC +TYPE (SUBDEVICE_TYPE), POINTER :: SDV +TYPE (LAGRANGIAN_PARTICLE_CLASS_TYPE), POINTER :: LPC +TYPE (LAGRANGIAN_PARTICLE_TYPE), POINTER :: LP ! If any device has QUANTITY='PARTICLE FLUX N', pre-compute PARTICLE fluxes @@ -6716,13 +6735,14 @@ SUBROUTINE UPDATE_DEVICES_1(T,DT,NM) K = MIN( KBP1, MAX(0, DV%K(1)) ) SDV%VALUE_1 = GAS_PHASE_OUTPUT(T,DT,NM,I,J,K,DV%QUANTITY_INDEX(1),Y_INDEX=DV%Y_INDEX,Z_INDEX=DV%Z_INDEX,& ELEM_INDX=DV%ELEM_INDEX,PART_INDEX=DV%PART_CLASS_INDEX,VELO_INDEX=DV%VELO_INDEX,& - PIPE_INDEX=DV%PIPE_INDEX,PROP_INDEX=DV%PROP_INDEX,REAC_INDEX=DV%REAC_INDEX) + PIPE_INDEX=DV%PIPE_INDEX,PROP_INDEX=DV%PROP_INDEX,REAC_INDEX=DV%REAC_INDEX,& + DEVC_INDEX=N) IF (DV%N_QUANTITY>1) & SDV%VALUE_2 = GAS_PHASE_OUTPUT(T,DT,NM,DV%I(2),DV%J(2),DV%K(2),DV%QUANTITY_INDEX(2),Y_INDEX=DV%Y_INDEX,& Z_INDEX=DV%Z_INDEX,ELEM_INDX=DV%ELEM_INDEX,PART_INDEX=DV%PART_CLASS_INDEX,& VELO_INDEX=DV%VELO_INDEX,PIPE_INDEX=DV%PIPE_INDEX,PROP_INDEX=DV%PROP_INDEX,& - REAC_INDEX=DV%REAC_INDEX) + REAC_INDEX=DV%REAC_INDEX,DEVC_INDEX=N) CASE DEFAULT GAS_STATS_SELECT @@ -6767,7 +6787,8 @@ SUBROUTINE UPDATE_DEVICES_1(T,DT,NM) VALUE = GAS_PHASE_OUTPUT(T,DT,NM,I,J,K,DV%QUANTITY_INDEX(1),Y_INDEX=DV%Y_INDEX,Z_INDEX=DV%Z_INDEX,& ELEM_INDX=DV%ELEM_INDEX,PART_INDEX=DV%PART_CLASS_INDEX,VELO_INDEX=DV%VELO_INDEX,& - PIPE_INDEX=DV%PIPE_INDEX,PROP_INDEX=DV%PROP_INDEX,REAC_INDEX=DV%REAC_INDEX) + PIPE_INDEX=DV%PIPE_INDEX,PROP_INDEX=DV%PROP_INDEX,REAC_INDEX=DV%REAC_INDEX,& + DEVC_INDEX=N) STATISTICS_SELECT: SELECT CASE(DV%SPATIAL_STATISTIC) CASE('MAX','MAXLOC X','MAXLOC Y','MAXLOC Z') IF (VALUE>SDV%VALUE_1) THEN @@ -7333,7 +7354,7 @@ END SUBROUTINE UPDATE_DEVICES_2 !> \param ICC_IN,JCC_IN Optional indexes of cut-cell. REAL(EB) RECURSIVE FUNCTION GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,IND,Y_INDEX,Z_INDEX,ELEM_INDX,PART_INDEX,VELO_INDEX,PIPE_INDEX,& - PROP_INDEX,REAC_INDEX,ICC_IN,JCC_IN) RESULT(GAS_PHASE_OUTPUT_RES) + PROP_INDEX,REAC_INDEX,DEVC_INDEX,ICC_IN,JCC_IN) RESULT(GAS_PHASE_OUTPUT_RES) USE MEMORY_FUNCTIONS, ONLY: REALLOCATE USE MATH_FUNCTIONS, ONLY: INTERPOLATE1D,INTERPOLATE1D_UNIFORM,UPDATE_HISTOGRAM @@ -7350,7 +7371,8 @@ REAL(EB) RECURSIVE FUNCTION GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,IND,Y_INDEX,Z_INDE REAL(EB), INTENT(IN) :: T,DT INTEGER, INTENT(IN) :: II,JJ,KK,IND,NM -INTEGER, INTENT(IN), OPTIONAL :: Y_INDEX,Z_INDEX,ELEM_INDX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX,ICC_IN,JCC_IN +INTEGER, INTENT(IN), OPTIONAL :: Y_INDEX,Z_INDEX,ELEM_INDX,PART_INDEX,VELO_INDEX,PIPE_INDEX,PROP_INDEX,REAC_INDEX,DEVC_INDEX,& + ICC_IN,JCC_IN REAL(EB) :: H_TC,TMP_TC,RE_D,NUSSELT,VEL,K_G,MU_G,COSTHETA,FAC,& Q_SUM,TMP_G,UU,VV,WW,VEL2,Y_MF_INT,PATHLENGTH,EXT_COEF,MASS_EXT_COEF,ZZ_FUEL,ZZ_OX,& VELSR,WATER_VOL_FRAC,WATER_TEMP,RHS,DT_C,DT_E,T_RATIO,Y_E_LAG, H_G,H_G_SUM,CPBAR,CP,ZZ_GET(1:N_TRACKED_SPECIES),RCON,& @@ -7369,6 +7391,10 @@ REAL(EB) RECURSIVE FUNCTION GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,IND,Y_INDEX,Z_INDE TYPE(BOUNDARY_ONE_D_TYPE), POINTER :: ONE_D TYPE(PROPERTY_TYPE), POINTER :: PY TYPE(DEVICE_TYPE), POINTER :: DV,DV2 +TYPE(SUBDEVICE_TYPE), POINTER :: SDV +TYPE(SURFACE_TYPE),POINTER :: SF +TYPE(LAGRANGIAN_PARTICLE_CLASS_TYPE), POINTER :: LPC +TYPE(LAGRANGIAN_PARTICLE_TYPE), POINTER :: LP ! Get species mass fraction if necessary @@ -7376,29 +7402,45 @@ REAL(EB) RECURSIVE FUNCTION GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,IND,Y_INDEX,Z_INDE R_Y_H2O = 0._EB Y_SPECIES = 1._EB -IF (Z_INDEX > 0) THEN - Y_SPECIES = ZZ(II,JJ,KK,Z_INDEX) - RCON = SPECIES_MIXTURE(Z_INDEX)%RCON -ELSEIF (Y_INDEX > 0) THEN - ZZ_GET(1:N_TRACKED_SPECIES) = ZZ(II,JJ,KK,1:N_TRACKED_SPECIES) - RCON = SPECIES(Y_INDEX)%RCON - CALL GET_MASS_FRACTION(ZZ_GET,Y_INDEX,Y_SPECIES) +IF (PRESENT(Z_INDEX)) THEN + IF (Z_INDEX > 0) THEN + Y_SPECIES = ZZ(II,JJ,KK,Z_INDEX) + RCON = SPECIES_MIXTURE(Z_INDEX)%RCON + ENDIF +ENDIF + +IF (PRESENT(Y_INDEX)) THEN + IF (Y_INDEX > 0) THEN + ZZ_GET(1:N_TRACKED_SPECIES) = ZZ(II,JJ,KK,1:N_TRACKED_SPECIES) + RCON = SPECIES(Y_INDEX)%RCON + CALL GET_MASS_FRACTION(ZZ_GET,Y_INDEX,Y_SPECIES) + ENDIF ENDIF + IF (DRY .AND. H2O_INDEX > 0) THEN ZZ_GET(1:N_TRACKED_SPECIES) = ZZ(II,JJ,KK,1:N_TRACKED_SPECIES) CALL GET_MASS_FRACTION(ZZ_GET,H2O_INDEX,Y_H2O) R_Y_H2O = SPECIES(H2O_INDEX)%RCON * Y_H2O IF (Y_INDEX==H2O_INDEX) Y_SPECIES=0._EB ENDIF + IF (PRESENT(PROP_INDEX)) THEN PY => PROPERTY(PROP_INDEX) ELSE PY => PROPERTY(0) ENDIF +IF (PRESENT(DEVC_INDEX)) THEN + DV => DEVICE(DEVC_INDEX) + IF (DV%N_SUBDEVICES>0) THEN + IF (DV%SUBDEVICE_INDEX(NM)>0) SDV => DV%SUBDEVICE(DV%SUBDEVICE_INDEX(NM)) + ENDIF +ENDIF + ! Get desired output value IND_SELECT: SELECT CASE(IND) + CASE DEFAULT ! SMOKE/WATER GAS_PHASE_OUTPUT_RES = 0._EB CASE( 1) ! DENSITY @@ -7920,7 +7962,7 @@ REAL(EB) RECURSIVE FUNCTION GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,IND,Y_INDEX,Z_INDE ! McCaffrey and Heskestad, A Robust Bidirectional Low-Velocity Probe for Flame and Fire Application ! Combustion and Flame, 26, 125 - 127, (1976). IF (PY%TC) THEN - PROBE_TMP = GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,110,PROP_INDEX=PROP_INDEX) + TMPM + PROBE_TMP = GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,110,PROP_INDEX=PROP_INDEX,DEVC_INDEX=DEVC_INDEX) + TMPM ELSE PROBE_TMP = TMP(II,JJ,KK) ENDIF @@ -8816,6 +8858,11 @@ REAL(EB) FUNCTION SOLID_PHASE_OUTPUT(INDX,T,NM,Y_INDEX,Z_INDEX,PART_INDEX,OPT_WA TYPE(PROPERTY_TYPE), POINTER :: PY TYPE(CFACE_TYPE), POINTER :: CFA TYPE(DEVICE_TYPE), POINTER :: DV +TYPE(WALL_TYPE), POINTER :: WC +TYPE(MATERIAL_TYPE),POINTER :: ML +TYPE(SURFACE_TYPE),POINTER :: SF +TYPE(LAGRANGIAN_PARTICLE_CLASS_TYPE), POINTER :: LPC +TYPE(LAGRANGIAN_PARTICLE_TYPE), POINTER :: LP ! Assign default value to output @@ -9781,6 +9828,7 @@ REAL(EB) FUNCTION PARTICLE_OUTPUT(NM,T,IND,LP_INDEX,Y_INDEX,Z_INDEX) INTEGER, INTENT(IN), OPTIONAL :: Y_INDEX,Z_INDEX REAL(EB), INTENT(IN) :: T TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC +TYPE(LAGRANGIAN_PARTICLE_TYPE), POINTER :: LP LP => LAGRANGIAN_PARTICLE(LP_INDEX) @@ -9859,6 +9907,7 @@ SUBROUTINE DUMP_DEVICES(T) REAL(EB) :: STIME,DI,DD,VALUE INTEGER :: I,J,N,NN,N_OUT REAL(EB) :: DEVC_TIME,CONST,CUMSUM,COORD_FACTOR +LOGICAL :: OPN TYPE(PROPERTY_TYPE), POINTER :: PY TYPE(DEVICE_TYPE), POINTER :: DV @@ -10076,6 +10125,9 @@ SUBROUTINE DUMP_PROF(T,NM) INTEGER, DIMENSION(0:NWP_MAX+1) :: LAYER_INDEX TYPE (PROFILE_TYPE), POINTER :: PF TYPE(BOUNDARY_ONE_D_TYPE), POINTER :: ONE_D +TYPE(WALL_TYPE), POINTER :: WC +TYPE(SURFACE_TYPE),POINTER :: SF +TYPE(LAGRANGIAN_PARTICLE_TYPE), POINTER :: LP PROF_LOOP: DO N=1,N_PROF @@ -10319,6 +10371,9 @@ SUBROUTINE UPDATE_HRR(DT,NM) TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1 TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC TYPE(CFACE_TYPE), POINTER :: CFA +TYPE(WALL_TYPE), POINTER :: WC +TYPE(LAGRANGIAN_PARTICLE_CLASS_TYPE), POINTER :: LPC +TYPE(LAGRANGIAN_PARTICLE_TYPE), POINTER :: LP ! Compute volume integral of certain quantities like the HRR @@ -10868,6 +10923,7 @@ SUBROUTINE FLUSH_GLOBAL_BUFFERS USE COMP_FUNCTIONS, ONLY : CURRENT_TIME INTEGER :: N REAL(EB) :: TNOW +LOGICAL :: OPN TNOW = CURRENT_TIME() @@ -11022,6 +11078,9 @@ SUBROUTINE COMPUTE_PARTICLE_FLUXES INTEGER :: IP,IC,IW REAL(EB) :: DROPMASS TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC +TYPE(WALL_TYPE), POINTER :: WC +TYPE(LAGRANGIAN_PARTICLE_CLASS_TYPE), POINTER :: LPC +TYPE(LAGRANGIAN_PARTICLE_TYPE), POINTER :: LP WFX => WORK4 ; WFX = 0._EB WFY => WORK5 ; WFY = 0._EB