Merge pull request #14831 from mcgratta/master

FDS Source: Fix a few bugs in HT3D

Author: mcgratta

Source/init.f90

@@ -1587,9 +1587,8 @@ SUBROUTINE REALLOCATE_ONE_D_ARRAYS(NM,WALL_CELL,THIN_WALL_CELL)
 INTEGER, ALLOCATABLE, DIMENSION(:) :: LAYER_INDEX
 INTEGER, ALLOCATABLE, DIMENSION(:) :: N_LAYER_CELLS_OLD
 REAL(EB) :: MIN_DENSITY
-REAL(EB), DIMENSION(MAX_LAYERS) :: LAYER_DENSITY
 TYPE(MATERIAL_TYPE), POINTER :: ML
-REAL(EB), ALLOCATABLE, DIMENSION(:) :: X_S_OLD
+REAL(EB), ALLOCATABLE, DIMENSION(:) :: X_S_OLD,LAYER_DENSITY
 LOGICAL, ALLOCATABLE, DIMENSION(:) :: REMESH_LAYER
 TYPE(WALL_TYPE), POINTER :: WC
 TYPE(THIN_WALL_TYPE), POINTER :: TW
@@ -1648,7 +1647,7 @@ SUBROUTINE REALLOCATE_ONE_D_ARRAYS(NM,WALL_CELL,THIN_WALL_CELL)
 
 ONE_D%N_CELLS_INI = 0
 ONE_D%N_CELLS_MAX = 0
-LAYER_DENSITY = 0._EB
+ALLOCATE(LAYER_DENSITY(ONE_D%N_LAYERS)) ; LAYER_DENSITY = 0._EB
 
 LAYER_LOOP: DO NL=1,ONE_D%N_LAYERS
 
@@ -1754,6 +1753,7 @@ SUBROUTINE REALLOCATE_ONE_D_ARRAYS(NM,WALL_CELL,THIN_WALL_CELL)
    ENDDO MATERIAL_LOOP3
 ENDDO POINT_LOOP3
 
+DEALLOCATE(LAYER_DENSITY)
 DEALLOCATE(LAYER_INDEX)
 
 ! Reset emissivity at the surface
@@ -3828,13 +3828,13 @@ SUBROUTINE FIND_WALL_BACK_INDEX(NM,IW)
 INTEGER :: II,JJ,KK,IC,ICG,IOR,NOM,ITER,OBST_INDEX,OBST_INDEX_PREVIOUS,NN,NNN,NL,N_LAYERS_OBST,&
            N_MATL_OBST,N_LAYERS,N_MATLS,IIF,JJF,KKF
 INTEGER, DIMENSION(MAX_MATERIALS) :: MATL_INDEX_OBST,MATL_INDEX
-REAL(EB), DIMENSION(MAX_LAYERS,MAX_MATERIALS) :: MATL_MASS_FRACTION_OBST,MATL_MASS_FRACTION
-REAL(EB), DIMENSION(0:MAX_LAYERS) :: LAYER_THICKNESS,MIN_LAYER_THICKNESS
-REAL(EB), DIMENSION(MAX_LAYERS) :: LAYER_THICKNESS_OBST,HEAT_SOURCE,HEAT_SOURCE_OBST,&
+REAL(EB), DIMENSION(MAX_LAYERS_HT3D,MAX_MATERIALS) :: MATL_MASS_FRACTION_OBST,MATL_MASS_FRACTION
+REAL(EB), DIMENSION(MAX_LAYERS_HT3D) :: LAYER_THICKNESS,MIN_LAYER_THICKNESS,&
+                                   LAYER_THICKNESS_OBST,HEAT_SOURCE,HEAT_SOURCE_OBST,&
                                    STRETCH_FACTOR,STRETCH_FACTOR_OBST,CELL_SIZE,CELL_SIZE_OBST,&
                                    CELL_SIZE_FACTOR,CELL_SIZE_FACTOR_OBST
-INTEGER, DIMENSION(MAX_LAYERS) :: N_LAYER_CELLS_MAX,N_LAYER_CELLS_MAX_OBST,RAMP_IHS_INDEX,RAMP_IHS_INDEX_OBST
-LOGICAL, DIMENSION(MAX_LAYERS) :: HT3D_LAYER
+INTEGER, DIMENSION(MAX_LAYERS_HT3D) :: N_LAYER_CELLS_MAX,N_LAYER_CELLS_MAX_OBST,RAMP_IHS_INDEX,RAMP_IHS_INDEX_OBST
+LOGICAL, DIMENSION(MAX_LAYERS_HT3D) :: HT3D_LAYER
 REAL(EB) :: XXC,YYC,ZZC,THICKNESS,OLD_THICKNESS,FRONT_LINING_THICKNESS,BACK_LINING_THICKNESS,LAYER_THICKNESS_OBST_TOTAL
 CHARACTER(MESSAGE_LENGTH) :: MESSAGE
 LOGICAL :: THIN_OBSTRUCTION
@@ -3994,10 +3994,7 @@ SUBROUTINE FIND_WALL_BACK_INDEX(NM,IW)
 
       IF (OBST_INDEX/=OBST_INDEX_PREVIOUS .AND. OBST_INDEX_PREVIOUS>0 .AND. OBST_INDEX>0) THEN
          OB_PREV => OM_PREV%OBSTRUCTION(OBST_INDEX_PREVIOUS)
-         IF ( (ANY(OB%MATL_MASS_FRACTION(:)/=OB_PREV%MATL_MASS_FRACTION(:),DIM=1)) .OR. &
-              (ANY(OB%MATL_INDEX(:)        /=OB_PREV%MATL_INDEX(:)        ,DIM=1)) .OR. &
-              (OB%HEAT_SOURCE/=OB_PREV%HEAT_SOURCE) .OR. &
-              (OB%RAMP_IHS_INDEX/=OB_PREV%RAMP_IHS_INDEX) ) THEN
+         IF (OB%ORDINAL/=OB_PREV%ORDINAL) THEN
             N_LAYERS_OBST = N_LAYERS_OBST + 1
             LAYER_THICKNESS_OBST(N_LAYERS_OBST) = 0._EB
             HEAT_SOURCE_OBST(N_LAYERS_OBST)        = OB%HEAT_SOURCE    
@@ -4238,7 +4235,7 @@ SUBROUTINE FIND_THIN_WALL_BACK_INDEX(NM,ITW)
 INTEGER, INTENT(IN) :: NM,ITW
 INTEGER :: II,JJ,KK,IC,IOR,IEC,ITW2,NOM,NN,NNN,NL,N_MATLS,IIGM,IIGP,JJGM,JJGP,KKGM,KKGP,ICM,ICP
 INTEGER, DIMENSION(MAX_MATERIALS) :: MATL_INDEX
-REAL(EB), DIMENSION(MAX_LAYERS,MAX_MATERIALS) :: MATL_MASS_FRACTION
+REAL(EB), DIMENSION(MAX_LAYERS_HT3D,MAX_MATERIALS) :: MATL_MASS_FRACTION
 REAL(EB) :: XXC,YYC,ZZC
 TYPE (MESH_TYPE), POINTER :: M
 TYPE (THIN_WALL_TYPE), POINTER :: TW

Source/prec.f90

@@ -10,6 +10,7 @@ MODULE PRECISION_PARAMETERS
 INTEGER, PARAMETER :: MAX_LPC=20                     !< Maximum number of declared particle classes
 INTEGER, PARAMETER :: MAX_SPECIES=20                 !< Maximum number of declared species
 INTEGER, PARAMETER :: MAX_LAYERS=20                  !< Maximum number of solid material layers
+INTEGER, PARAMETER :: MAX_LAYERS_HT3D=500            !< Maximum number of solid material layers for an HT3D solid
 INTEGER, PARAMETER :: MAX_MATERIALS=20               !< Maximum number of solid material components
 INTEGER, PARAMETER :: MAX_MATERIALS_TOTAL=400        !< Dimension of material work array
 INTEGER, PARAMETER :: MAX_CONE_CURVES=10             !< Maximum number of cone calorimeter curves

Source/wall.f90

@@ -1593,7 +1593,7 @@ SUBROUTINE DEPOSIT_PARTICLE_MASS(LP,LPC)
 
 IF (SF%PYROLYSIS_MODEL==PYROLYSIS_PREDICTED) THEN
    ONE_D => BOUNDARY_ONE_D(LP%OD_INDEX)
-   LP%RADIUS = SUM(ONE_D%LAYER_THICKNESS(1:SF%N_LAYERS))
+   LP%RADIUS = SUM(ONE_D%LAYER_THICKNESS(1:ONE_D%N_LAYERS))
 ENDIF
 
 END SUBROUTINE DEPOSIT_PARTICLE_MASS
@@ -1838,18 +1838,18 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX,
 REAL(EB), DIMENSION(MAX_LPC) :: Q_DOT_PART_S,M_DOT_PART_S
 REAL(EB), DIMENSION(NWP_MAX) :: TMP_S,RHO_H_S
 REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: RHO_DOT,INT_WGT
-REAL(EB), DIMENSION(MAX_LAYERS) :: DX_MIN
-REAL(EB), DIMENSION(MAX_LAYERS,MAX_MATERIALS) :: RHO_ADJUSTED
+REAL(EB), DIMENSION(MAX_LAYERS_HT3D) :: DX_MIN
+REAL(EB), DIMENSION(MAX_LAYERS_HT3D,MAX_MATERIALS) :: RHO_ADJUSTED
 REAL(EB), DIMENSION(NWP_MAX) :: AAS,BBS,CCS,DDS,DDT,Q_S,Q_IR,Q_ADD,TWO_DX_KAPPA_S,DX_S,MF_FRAC,REGRID_FACTOR
 REAL(EB), DIMENSION(0:NWP_MAX+1) :: RHO_S,DELTA_TMP,RDX_S
 REAL(EB), DIMENSION(0:NWP_MAX) :: X_S_NEW,RDXN_S,R_S,R_S_NEW,DX_WGT_S
 INTEGER, DIMENSION(0:NWP_MAX+1) :: LAYER_INDEX
-INTEGER, DIMENSION(MAX_LAYERS) :: N_LAYER_CELLS_NEW
+INTEGER, DIMENSION(MAX_LAYERS_HT3D) :: N_LAYER_CELLS_NEW
 INTEGER  :: NWP_NEW,I_GRAD,IZERO,SURF_INDEX,BACKING,NWP,I,NL,N,OBST_INDEX,&
             N_CELLS,ITMP,ITER,BACK_MESH,BACK_INDEX,BACK_WALL_INDEX
 CHARACTER(MESSAGE_LENGTH) :: MESSAGE
-LOGICAL :: ISOLATED_THIN_WALL,ISOLATED_THIN_WALL_BACK,REMESH_LAYER(MAX_LAYERS),REMESH_CHECK,DIRICHLET_BACK,&
-           CELL_ZERO(MAX_LAYERS),TMP_CHECK(MAX_LAYERS)
+LOGICAL :: ISOLATED_THIN_WALL,ISOLATED_THIN_WALL_BACK,REMESH_LAYER(MAX_LAYERS_HT3D),REMESH_CHECK,DIRICHLET_BACK,&
+           CELL_ZERO(MAX_LAYERS_HT3D),TMP_CHECK(MAX_LAYERS_HT3D)
 TYPE(WALL_TYPE), POINTER :: WC,WC_BACK
 TYPE(THIN_WALL_TYPE), POINTER :: TW,TW_BACK
 TYPE(CFACE_TYPE), POINTER :: CFA,CFA_BACK
@@ -2014,7 +2014,7 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX,
 
 DO NL=1,ONE_D%N_LAYERS
    DO N=1,ONE_D%N_MATL
-      RHO_ADJUSTED(NL,N) = MATERIAL(ONE_D%MATL_INDEX(N))%RHO_S*SF%DENSITY_ADJUST_FACTOR(NL,N)
+      RHO_ADJUSTED(NL,N) = MATERIAL(ONE_D%MATL_INDEX(N))%RHO_S*SF%DENSITY_ADJUST_FACTOR(MIN(MAX_LAYERS,NL),N)
    ENDDO
 ENDDO
 
@@ -2071,7 +2071,7 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX,
          ONE_D%K_S(I) = ONE_D%K_S(I) + ONE_D%MATL_COMP(N)%RHO(I)*ML%K_S(ITMP)/RHO_ADJUSTED(LAYER_INDEX(I),N)
          ONE_D%RHO_C_S(I) = ONE_D%RHO_C_S(I) + ONE_D%MATL_COMP(N)%RHO(I)*ML%C_S(ITMP)
       ENDDO MATERIAL_LOOP0
-      IF (SF%PACKING_RATIO(LAYER_INDEX(I))>0._EB) ONE_D%K_S(I) = ONE_D%K_S(I)*SF%PACKING_RATIO(LAYER_INDEX(I))
+      IF (SF%BOUNDARY_FUEL_MODEL) ONE_D%K_S(I) = ONE_D%K_S(I)*SF%PACKING_RATIO(LAYER_INDEX(I))
 
       IF (VOLSUM > 0._EB) THEN
          ONE_D%K_S(I) = ONE_D%K_S(I)/VOLSUM
@@ -2121,7 +2121,7 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX,
    ELSEIF (PRESENT(CFACE_INDEX)) THEN
       HTCF = B1%HEAT_TRANS_COEF
    ELSEIF (PRESENT(PARTICLE_INDEX)) THEN
-      RADIUS = SF%INNER_RADIUS + SUM(ONE_D%LAYER_THICKNESS(1:SF%N_LAYERS))
+      RADIUS = SF%INNER_RADIUS + SUM(ONE_D%LAYER_THICKNESS(1:ONE_D%N_LAYERS))
       SELECT CASE(SF%GEOMETRY)
          CASE (SURF_CARTESIAN)                          ; HTC_LIMIT = 0.5_EB*RADIUS*ONE_D%RHO_C_S(1)/(      DT_BC_SUB)
          CASE (SURF_CYLINDRICAL,SURF_INNER_CYLINDRICAL) ; HTC_LIMIT = 0.5_EB*RADIUS*ONE_D%RHO_C_S(1)/(2._EB*DT_BC_SUB)
@@ -2158,7 +2158,7 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX,
          Q_CON_B = HTCB*DTMP
          Q_RAD_IN_B = EMISSIVITY_BACK*SIGMA*TMP_GAS_BACK**4
          Q_LIQUID_B = 0._EB
-         LAYER_DIVIDE = REAL(SF%N_LAYERS+1)
+         LAYER_DIVIDE = REAL(ONE_D%N_LAYERS+1)
          MF_FRAC = 1._EB
 
       CASE(INSULATED)  ! No heat transfer out the back
@@ -2774,7 +2774,7 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX,
          ONE_D%RHO_C_S(I) = ONE_D%RHO_C_S(I) + ONE_D%MATL_COMP(N)%RHO(I)*ML%C_S(ITMP)
          RHO_S(I) = RHO_S(I) + ONE_D%MATL_COMP(N)%RHO(I)
       ENDDO MATERIAL_LOOP3
-      IF (SF%PACKING_RATIO(LAYER_INDEX(I))>0._EB) ONE_D%K_S(I) = ONE_D%K_S(I)*SF%PACKING_RATIO(LAYER_INDEX(I))
+      IF (SF%BOUNDARY_FUEL_MODEL) ONE_D%K_S(I) = ONE_D%K_S(I)*SF%PACKING_RATIO(LAYER_INDEX(I))
       IF (VOLSUM > 0._EB) ONE_D%K_S(I) = ONE_D%K_S(I)/VOLSUM
       IF (ONE_D%K_S(I)<=TWO_EPSILON_EB)      ONE_D%K_S(I)      = 10000._EB
       IF (ONE_D%RHO_C_S(I)<=TWO_EPSILON_EB)  ONE_D%RHO_C_S(I)  = 0.001_EB
@@ -3387,7 +3387,7 @@ SUBROUTINE PYROLYSIS(N_MATS,MATL_INDEX,SURF_INDEX,IIG,JJG,KKG,TMP_S,TMP_F,Y_O2_F
                      MFLUX = MIN(MFLUX_MAX,MFLUX + RHO_DOT_EXTRA*DX_S(SOLID_CELL_INDEX))
                      RHO_DOT_REAC(J) =MFLUX/DX_S(SOLID_CELL_INDEX) ! kg/m3/s
                   CASE(SURF_SPHERICAL)
-                     NWP = SUM(ONE_D%N_LAYER_CELLS(1:SF%N_LAYERS))
+                     NWP = SUM(ONE_D%N_LAYER_CELLS(1:ONE_D%N_LAYERS))
                      R_S_0 = SF%INNER_RADIUS + ONE_D%X(NWP) - ONE_D%X(0)
                      R_S_1 = SF%INNER_RADIUS + ONE_D%X(NWP) - ONE_D%X(1)
                      DR = (R_S_0**3-R_S_1**3)/(3._EB*R_S_0**2)
@@ -3445,7 +3445,7 @@ SUBROUTINE PYROLYSIS(N_MATS,MATL_INDEX,SURF_INDEX,IIG,JJG,KKG,TMP_S,TMP_F,Y_O2_F
             IF (SF%BOUNDARY_FUEL_MODEL) THEN
                LENGTH_SCALE = 1._EB/(SF%SURFACE_VOLUME_RATIO(LAYER_INDEX)*SF%PACKING_RATIO(LAYER_INDEX))
             ELSE
-               LENGTH_SCALE = SF%INNER_RADIUS + SUM(ONE_D%LAYER_THICKNESS(1:SF%N_LAYERS))
+               LENGTH_SCALE = SF%INNER_RADIUS + SUM(ONE_D%LAYER_THICKNESS(1:ONE_D%N_LAYERS))
                SELECT CASE(SF%GEOMETRY)
                   CASE(SURF_SPHERICAL)
                      LENGTH_SCALE = LENGTH_SCALE/3._EB
@@ -3808,7 +3808,7 @@ SUBROUTINE HT3D_TEMPERATURE_EXCHANGE(NM)
    BC => M%BOUNDARY_COORD(TW%BC_INDEX)
    IF (ABS(BC%IOR)==M%HT_3D_SWEEP_DIRECTION) CYCLE THIN_WALL_LOOP
    ONE_D => M%BOUNDARY_ONE_D(TW%OD_INDEX)
-   NWP = SUM(ONE_D%N_LAYER_CELLS(1:SF%N_LAYERS))
+   NWP = SUM(ONE_D%N_LAYER_CELLS(1:ONE_D%N_LAYERS))
    THR_D => M%BOUNDARY_THR_D(TW%TD_INDEX)
    IF (.NOT.ALLOCATED(THR_D%NODE)) CYCLE THIN_WALL_LOOP