Merge pull request #13789 from rmcdermo/master

FDS Source: IMPORTANT add flux limter face correction to maintain iso…
firemodels · Nov 21, 2024 · 0b52ed4 · 0b52ed4
2 parents 7843940 + 898c602
commit 0b52ed4
Show file tree

Hide file tree

Showing 4 changed files with 138 additions and 5 deletions.
diff --git a/Source/cons.f90 b/Source/cons.f90
@@ -270,6 +270,7 @@ MODULE GLOBAL_CONSTANTS
 LOGICAL :: TENSOR_DIFFUSIVITY=.FALSE.               !< If true, use experimental tensor diffusivity model for spec and tmp
 LOGICAL :: OXPYRO_MODEL=.FALSE.                     !< Flag to use oxidative pyrolysis mass transfer model
 LOGICAL :: OUTPUT_WALL_QUANTITIES=.FALSE.           !< Flag to force call to WALL_MODEL
+LOGICAL :: TEST_FLUX_LIMITER_FACE_CORRECTION=.FALSE.
 
 INTEGER, ALLOCATABLE, DIMENSION(:) :: CHANGE_TIME_STEP_INDEX      !< Flag to indicate if a mesh needs to change time step
 INTEGER, ALLOCATABLE, DIMENSION(:) :: SETUP_PRESSURE_ZONES_INDEX  !< Flag to indicate if a mesh needs to keep searching for ZONEs

diff --git a/Source/init.f90 b/Source/init.f90
@@ -32,7 +32,7 @@ SUBROUTINE INITIALIZE_MESH_VARIABLES_1(DT,NM)
 USE RADCONS, ONLY: UIIDIM
 USE CONTROL_VARIABLES
 USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP
-INTEGER :: N,I,J,K,IW,IC,SURF_INDEX,IOR,IERR,IZERO,II,JJ,KK,OBST_INDEX,N_EXTERNAL_CELLS,NS
+INTEGER :: N,I,J,K,IW,IC,SURF_INDEX,IOR,IERR,IZERO,II,JJ,KK,OBST_INDEX,N_EXTERNAL_CELLS,NS,N_LOWER_SCALARS
 REAL(EB), INTENT(IN) :: DT
 INTEGER, INTENT(IN) :: NM
 REAL(EB) :: MU_N,CS,DELTA,INTEGRAL,TEMP,ZSW
@@ -187,11 +187,18 @@ SUBROUTINE INITIALIZE_MESH_VARIABLES_1(DT,NM)
 
 ! Allocate scalar face values
 
-ALLOCATE( M%FX(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO)
+! Required for cell face density correction for multicomponent mixtures
+IF (TEST_FLUX_LIMITER_FACE_CORRECTION .AND. N_TRACKED_SPECIES>2) THEN
+   N_LOWER_SCALARS=0
+ELSE
+   N_LOWER_SCALARS=1
+ENDIF
+
+ALLOCATE( M%FX(0:IBP1,0:JBP1,0:KBP1,N_LOWER_SCALARS:N_TOTAL_SCALARS),STAT=IZERO)
 CALL ChkMemErr('INIT','FX',IZERO)
-ALLOCATE( M%FY(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO)
+ALLOCATE( M%FY(0:IBP1,0:JBP1,0:KBP1,N_LOWER_SCALARS:N_TOTAL_SCALARS),STAT=IZERO)
 CALL ChkMemErr('INIT','FY',IZERO)
-ALLOCATE( M%FZ(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO)
+ALLOCATE( M%FZ(0:IBP1,0:JBP1,0:KBP1,N_LOWER_SCALARS:N_TOTAL_SCALARS),STAT=IZERO)
 CALL ChkMemErr('INIT','FZ',IZERO)
 M%FX = 0._EB
 M%FY = 0._EB

diff --git a/Source/mass.f90 b/Source/mass.f90
@@ -177,6 +177,131 @@ SUBROUTINE MASS_FINITE_DIFFERENCES(NM)
 
 ENDDO SPECIES_LOOP
 
+FACE_CORRECTION_IF: IF (TEST_FLUX_LIMITER_FACE_CORRECTION .AND. N_TRACKED_SPECIES>2) THEN
+
+   ! Repeat the above for DENSITY
+
+   CALL GET_SCALAR_FACE_VALUE(UU,RHOP,FX(:,:,:,0),1,IBM1,1,JBAR,1,KBAR,1,I_FLUX_LIMITER)
+   CALL GET_SCALAR_FACE_VALUE(VV,RHOP,FY(:,:,:,0),1,IBAR,1,JBM1,1,KBAR,2,I_FLUX_LIMITER)
+   CALL GET_SCALAR_FACE_VALUE(WW,RHOP,FZ(:,:,:,0),1,IBAR,1,JBAR,1,KBM1,3,I_FLUX_LIMITER)
+
+   !$OMP PARALLEL DO PRIVATE(IW,WC,BC,B1,II,JJ,KK,IIG,JJG,KKG,IOR,IC,U_TEMP,Z_TEMP,F_TEMP)
+   WALL_LOOP_3: DO IW=1,N_EXTERNAL_WALL_CELLS+N_INTERNAL_WALL_CELLS
+      WC=>WALL(IW)
+      IF (WC%BOUNDARY_TYPE==NULL_BOUNDARY) CYCLE WALL_LOOP_3
+      BC=>BOUNDARY_COORD(WC%BC_INDEX)
+      B1=>BOUNDARY_PROP1(WC%B1_INDEX)
+
+      II  = BC%II
+      JJ  = BC%JJ
+      KK  = BC%KK
+      IIG = BC%IIG
+      JJG = BC%JJG
+      KKG = BC%KKG
+      IOR = BC%IOR
+      IC  = CELL_INDEX(II,JJ,KK)
+
+      IF (WC%BOUNDARY_TYPE==SOLID_BOUNDARY .AND. .NOT.CELL(IC)%SOLID .AND. .NOT.CELL(IC)%EXTERIOR) THEN
+         SELECT CASE(IOR)
+            CASE( 1); FX(IIG-1,JJG,KKG,0) = 0._EB
+            CASE(-1); FX(IIG,JJG,KKG,0)   = 0._EB
+            CASE( 2); FY(IIG,JJG-1,KKG,0) = 0._EB
+            CASE(-2); FY(IIG,JJG,KKG,0)   = 0._EB
+            CASE( 3); FZ(IIG,JJG,KKG-1,0) = 0._EB
+            CASE(-3); FZ(IIG,JJG,KKG,0)   = 0._EB
+         END SELECT
+      ELSE
+         SELECT CASE(IOR)
+            CASE( 1); FX(IIG-1,JJG,KKG,0) = B1%RHO_F
+            CASE(-1); FX(IIG,JJG,KKG,0)   = B1%RHO_F
+            CASE( 2); FY(IIG,JJG-1,KKG,0) = B1%RHO_F
+            CASE(-2); FY(IIG,JJG,KKG,0)   = B1%RHO_F
+            CASE( 3); FZ(IIG,JJG,KKG-1,0) = B1%RHO_F
+            CASE(-3); FZ(IIG,JJG,KKG,0)   = B1%RHO_F
+         END SELECT
+      ENDIF
+
+      ! Overwrite first off-wall advective flux if flow is away from the wall and if the face is not also a wall cell
+
+      OFF_WALL_IF_3: IF (WC%BOUNDARY_TYPE/=INTERPOLATED_BOUNDARY .AND. WC%BOUNDARY_TYPE/=OPEN_BOUNDARY) THEN
+
+         OFF_WALL_SELECT_3: SELECT CASE(IOR)
+            CASE( 1) OFF_WALL_SELECT_3
+               !      ghost          FX/UU(II+1)
+               ! ///   II   ///  II+1  |  II+2  | ...
+               !                       ^ WALL_INDEX(II+1,+1)
+               IF ((UU(II+1,JJ,KK)>0._EB) .AND. .NOT.(CELL(CELL_INDEX(II+1,JJ,KK))%WALL_INDEX(+1)>0)) THEN
+                  Z_TEMP(0:3,1,1) = (/RHOP(II+1,JJ,KK),RHOP(II+1:II+2,JJ,KK),DUMMY/)
+                  U_TEMP(1,1,1) = UU(II+1,JJ,KK)
+                  CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER)
+                  FX(II+1,JJ,KK,0) = F_TEMP(1,1,1)
+               ENDIF
+            CASE(-1) OFF_WALL_SELECT_3
+               !            FX/UU(II-2)     ghost
+               ! ... |  II-2  |  II-1  ///   II   ///
+               !              ^ WALL_INDEX(II-1,-1)
+               IF ((UU(II-2,JJ,KK)<0._EB) .AND. .NOT.(CELL(CELL_INDEX(II-1,JJ,KK))%WALL_INDEX(-1)>0)) THEN
+                  Z_TEMP(0:3,1,1) = (/DUMMY,RHOP(II-2:II-1,JJ,KK),RHOP(II-1,JJ,KK)/)
+                  U_TEMP(1,1,1) = UU(II-2,JJ,KK)
+                  CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER)
+                  FX(II-2,JJ,KK,0) = F_TEMP(1,1,1)
+               ENDIF
+            CASE( 2) OFF_WALL_SELECT_3
+               IF ((VV(II,JJ+1,KK)>0._EB) .AND. .NOT.(CELL(CELL_INDEX(II,JJ+1,KK))%WALL_INDEX(+2)>0)) THEN
+                  Z_TEMP(1,0:3,1) = (/RHOP(II,JJ+1,KK),RHOP(II,JJ+1:JJ+2,KK),DUMMY/)
+                  U_TEMP(1,1,1) = VV(II,JJ+1,KK)
+                  CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER)
+                  FY(II,JJ+1,KK,0) = F_TEMP(1,1,1)
+               ENDIF
+            CASE(-2) OFF_WALL_SELECT_3
+               IF ((VV(II,JJ-2,KK)<0._EB) .AND. .NOT.(CELL(CELL_INDEX(II,JJ-1,KK))%WALL_INDEX(-2)>0)) THEN
+                  Z_TEMP(1,0:3,1) = (/DUMMY,RHOP(II,JJ-2:JJ-1,KK),RHOP(II,JJ-1,KK)/)
+                  U_TEMP(1,1,1) = VV(II,JJ-2,KK)
+                  CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER)
+                  FY(II,JJ-2,KK,0) = F_TEMP(1,1,1)
+               ENDIF
+            CASE( 3) OFF_WALL_SELECT_3
+               IF ((WW(II,JJ,KK+1)>0._EB) .AND. .NOT.(CELL(CELL_INDEX(II,JJ,KK+1))%WALL_INDEX(+3)>0)) THEN
+                  Z_TEMP(1,1,0:3) = (/RHOP(II,JJ,KK+1),RHOP(II,JJ,KK+1:KK+2),DUMMY/)
+                  U_TEMP(1,1,1) = WW(II,JJ,KK+1)
+                  CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER)
+                  FZ(II,JJ,KK+1,0) = F_TEMP(1,1,1)
+               ENDIF
+            CASE(-3) OFF_WALL_SELECT_3
+               IF ((WW(II,JJ,KK-2)<0._EB) .AND. .NOT.(CELL(CELL_INDEX(II,JJ,KK-1))%WALL_INDEX(-3)>0)) THEN
+                  Z_TEMP(1,1,0:3) = (/DUMMY,RHOP(II,JJ,KK-2:KK-1),RHOP(II,JJ,KK-1)/)
+                  U_TEMP(1,1,1) = WW(II,JJ,KK-2)
+                  CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER)
+                  FZ(II,JJ,KK-2,0) = F_TEMP(1,1,1)
+               ENDIF
+         END SELECT OFF_WALL_SELECT_3
+
+      ENDIF OFF_WALL_IF_3
+
+   ENDDO WALL_LOOP_3
+   !$OMP END PARALLEL DO
+
+   ! Now correct the face value of (RHO*ZZ) such that SUM(RHO*ZZ)_FACE = RHO_FACE
+
+   !$OMP PARALLEL DO PRIVATE(N) SCHEDULE(STATIC)
+   DO K=0,KBAR
+      DO J=0,JBAR
+         DO I=0,IBAR
+            N=MAXLOC(FX(I,J,K,1:N_TRACKED_SPECIES),1)
+            FX(I,J,K,N) = MAX( 0._EB, FX(I,J,K,0) - SUM(FX(I,J,K,1:(N-1))) - SUM(FX(I,J,K,(N+1):N_TRACKED_SPECIES)) )
+
+            N=MAXLOC(FY(I,J,K,1:N_TRACKED_SPECIES),1)
+            FY(I,J,K,N) = MAX( 0._EB, FY(I,J,K,0) - SUM(FY(I,J,K,1:(N-1))) - SUM(FY(I,J,K,(N+1):N_TRACKED_SPECIES)) )
+
+            N=MAXLOC(FZ(I,J,K,1:N_TRACKED_SPECIES),1)
+            FZ(I,J,K,N) = MAX( 0._EB, FZ(I,J,K,0) - SUM(FZ(I,J,K,1:(N-1))) - SUM(FZ(I,J,K,(N+1):N_TRACKED_SPECIES)) )
+         ENDDO
+      ENDDO
+   ENDDO
+   !$OMP END PARALLEL DO
+
+ENDIF FACE_CORRECTION_IF
+
 T_USED(3)=T_USED(3)+CURRENT_TIME()-TNOW
 END SUBROUTINE MASS_FINITE_DIFFERENCES
 

diff --git a/Source/read.f90 b/Source/read.f90
@@ -1745,7 +1745,7 @@ SUBROUTINE READ_MISC
                 RAMP_UX,RAMP_UY,RAMP_UZ,RAMP_VX,RAMP_VY,RAMP_VZ,RAMP_WX,RAMP_WY,RAMP_WZ,&
                 RESTART,RESTART_CHID,SC,&
                 RND_SEED,SIMULATION_MODE,SMOKE3D_16,SMOKE_ALBEDO,SOLID_PHASE_ONLY,SOOT_DENSITY,SOOT_OXIDATION,&
-                TAU_DEFAULT,TENSOR_DIFFUSIVITY,TERRAIN_IMAGE,TEXTURE_ORIGIN,&
+                TAU_DEFAULT,TENSOR_DIFFUSIVITY,TERRAIN_IMAGE,TEST_FLUX_LIMITER_FACE_CORRECTION,TEXTURE_ORIGIN,&
                 THERMOPHORETIC_DEPOSITION,THERMOPHORETIC_SETTLING,THICKEN_OBSTRUCTIONS,&
                 TMPA,TURBULENCE_MODEL,TURBULENT_DEPOSITION,UVW_FILE,&
                 VERBOSE,VISIBILITY_FACTOR,VN_MAX,VN_MIN,Y_CO2_INFTY,Y_O2_INFTY,&