diff --git a/Source/cons.f90 b/Source/cons.f90 index a57a166052..6dfa66c920 100644 --- 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 index b9409cb536..61d52f985c 100644 --- 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 index f53a2fe1f5..6306f6047a 100644 --- 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 index 3234e05b10..507db9b043 100644 --- 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,&