diff --git a/Source/cons.f90 b/Source/cons.f90 index 7c1b4d36c2a..a82feb807aa 100644 --- a/Source/cons.f90 +++ b/Source/cons.f90 @@ -274,7 +274,6 @@ 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 :: FLUX_LIMITER_SINGLE_COEF=.TRUE. !< Flag to base flux coefficients off of a single worst-case scalar gradient LOGICAL :: STORE_FIRE_ARRIVAL=.FALSE. !< Flag for tracking arrival of spreading fire front over a surface LOGICAL :: STORE_FIRE_RESIDENCE=.FALSE. !< Flag for tracking residence time of spreading fire front over a surface LOGICAL :: STORE_LS_SPREAD_RATE=.FALSE. !< Flag for outputting local level set spread rate magnitude diff --git a/Source/divg.f90 b/Source/divg.f90 index c146536afd5..6c4d2dd524e 100644 --- a/Source/divg.f90 +++ b/Source/divg.f90 @@ -586,7 +586,7 @@ SUBROUTINE DIVERGENCE_PART_1(T,DT,NM) CONST_GAMMA_IF_1: IF (.NOT.CONSTANT_SPECIFIC_HEAT_RATIO) THEN - CALL ENTHALPY_ADVECTION(U_DOT_DEL_RHO_H_S) ! Compute u dot grad rho h_s + CALL ENTHALPY_ADVECTION_NEW(U_DOT_DEL_RHO_H_S) ! Compute u dot grad rho h_s DO K=1,KBAR DO J=1,JBAR @@ -635,11 +635,7 @@ SUBROUTINE DIVERGENCE_PART_1(T,DT,NM) CONST_GAMMA_IF_2: IF (.NOT.CONSTANT_SPECIFIC_HEAT_RATIO) THEN - IF (FLUX_LIMITER_SINGLE_COEF) THEN - CALL SPECIES_ADVECTION_PART_1_NEW - ELSE - CALL SPECIES_ADVECTION_PART_1 ! Compute and store face values of (rho Z_n) - ENDIF + CALL SPECIES_ADVECTION_PART_1_NEW ! Compute and store face values of (rho Z_n) DO N=1,N_TRACKED_SPECIES @@ -789,14 +785,14 @@ SUBROUTINE DIVERGENCE_PART_1(T,DT,NM) END SUBROUTINE DIVERGENCE_PART_1 -SUBROUTINE ENTHALPY_ADVECTION(U_DOT_DEL_RHO_H_S) +SUBROUTINE ENTHALPY_ADVECTION_NEW(U_DOT_DEL_RHO_H_S) USE PHYSICAL_FUNCTIONS, ONLY: GET_SENSIBLE_ENTHALPY -USE MATH_FUNCTIONS, ONLY: GET_SCALAR_FACE_VALUE +USE MATH_FUNCTIONS, ONLY: GET_SCALAR_FACE_COEF,GET_SCALAR_FACE_VALUE_NEW REAL(EB), POINTER, DIMENSION(:,:,:) :: FX_H_S,FY_H_S,FZ_H_S,RHO_H_S_P,U_DOT_DEL_RHO_H_S REAL(EB) :: UN,UN_P,TMP_F_GAS,DU_P,DU_M,DV_P,DV_M,DW_P,DW_M,DU,H_S REAL(EB), ALLOCATABLE, DIMENSION(:) :: ZZ_GET -REAL(EB), DIMENSION(0:3,0:3,0:3) :: U_TEMP,F_TEMP +REAL(EB), DIMENSION(0:3,0:3,0:3) :: U_TEMP,F_TEMP,B_TEMP REAL(EB), DIMENSION(-1:3,-1:3,-1:3) :: Z_TEMP INTEGER :: IC,I,J,K,IW TYPE(WALL_TYPE), POINTER :: WC @@ -809,6 +805,10 @@ SUBROUTINE ENTHALPY_ADVECTION(U_DOT_DEL_RHO_H_S) FZ_H_S=>WORK4 ; FZ_H_S = 0._EB U_DOT_DEL_RHO_H_S=>WORK6 ; U_DOT_DEL_RHO_H_S=0._EB +BFX = 1.E10_EB +BFY = 1.E10_EB +BFZ = 1.E10_EB + ! Compute and store rho*h_s !$OMP PARALLEL PRIVATE(ZZ_GET, H_S) @@ -829,9 +829,13 @@ SUBROUTINE ENTHALPY_ADVECTION(U_DOT_DEL_RHO_H_S) ! Compute scalar face values -CALL GET_SCALAR_FACE_VALUE(UU,RHO_H_S_P,FX_H_S,1,IBM1,1,JBAR,1,KBAR,1,I_FLUX_LIMITER) -CALL GET_SCALAR_FACE_VALUE(VV,RHO_H_S_P,FY_H_S,1,IBAR,1,JBM1,1,KBAR,2,I_FLUX_LIMITER) -CALL GET_SCALAR_FACE_VALUE(WW,RHO_H_S_P,FZ_H_S,1,IBAR,1,JBAR,1,KBM1,3,I_FLUX_LIMITER) +CALL GET_SCALAR_FACE_COEF(UU,RHO_H_S_P,BFX,0,IBAR,1,JBAR,1,KBAR,1,I_FLUX_LIMITER) +CALL GET_SCALAR_FACE_COEF(VV,RHO_H_S_P,BFY,1,IBAR,0,JBAR,1,KBAR,2,I_FLUX_LIMITER) +CALL GET_SCALAR_FACE_COEF(WW,RHO_H_S_P,BFZ,1,IBAR,1,JBAR,0,KBAR,3,I_FLUX_LIMITER) + +CALL GET_SCALAR_FACE_VALUE_NEW(UU,RHO_H_S_P,FX_H_S,BFX,0,IBAR,1,JBAR,1,KBAR,1,I_FLUX_LIMITER) +CALL GET_SCALAR_FACE_VALUE_NEW(VV,RHO_H_S_P,FY_H_S,BFY,1,IBAR,0,JBAR,1,KBAR,2,I_FLUX_LIMITER) +CALL GET_SCALAR_FACE_VALUE_NEW(WW,RHO_H_S_P,FZ_H_S,BFZ,1,IBAR,1,JBAR,0,KBAR,3,I_FLUX_LIMITER) ALLOCATE(ZZ_GET(1:N_TRACKED_SPECIES)) @@ -871,7 +875,8 @@ SUBROUTINE ENTHALPY_ADVECTION(U_DOT_DEL_RHO_H_S) IF ((UU(BC%II+1,BC%JJ,BC%KK)>0._EB) .AND. .NOT.(CELL(CELL_INDEX(BC%II+1,BC%JJ,BC%KK))%WALL_INDEX(+1)>0)) THEN Z_TEMP(0:2,1,1) = (/RHO_H_S_P(BC%II+1,BC%JJ,BC%KK),RHO_H_S_P(BC%II+1:BC%II+2,BC%JJ,BC%KK)/) U_TEMP(1,1,1) = UU(BC%II+1,BC%JJ,BC%KK) - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER) + B_TEMP(1,1,1) = BFX(BC%II+1,BC%JJ,BC%KK) + CALL GET_SCALAR_FACE_VALUE_NEW(U_TEMP,Z_TEMP,F_TEMP,B_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER) FX_H_S(BC%II+1,BC%JJ,BC%KK) = F_TEMP(1,1,1) ENDIF CASE(-1) OFF_WALL_SELECT_1 @@ -881,35 +886,40 @@ SUBROUTINE ENTHALPY_ADVECTION(U_DOT_DEL_RHO_H_S) IF ((UU(BC%II-2,BC%JJ,BC%KK)<0._EB) .AND. .NOT.(CELL(CELL_INDEX(BC%II-1,BC%JJ,BC%KK))%WALL_INDEX(-1)>0)) THEN Z_TEMP(1:3,1,1) = (/RHO_H_S_P(BC%II-2:BC%II-1,BC%JJ,BC%KK),RHO_H_S_P(BC%II-1,BC%JJ,BC%KK)/) U_TEMP(1,1,1) = UU(BC%II-2,BC%JJ,BC%KK) - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER) + B_TEMP(1,1,1) = BFX(BC%II-2,BC%JJ,BC%KK) + CALL GET_SCALAR_FACE_VALUE_NEW(U_TEMP,Z_TEMP,F_TEMP,B_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER) FX_H_S(BC%II-2,BC%JJ,BC%KK) = F_TEMP(1,1,1) ENDIF CASE( 2) OFF_WALL_SELECT_1 IF ((VV(BC%II,BC%JJ+1,BC%KK)>0._EB) .AND. .NOT.(CELL(CELL_INDEX(BC%II,BC%JJ+1,BC%KK))%WALL_INDEX(+2)>0)) THEN Z_TEMP(1,0:2,1) = (/RHO_H_S_P(BC%II,BC%JJ+1,BC%KK),RHO_H_S_P(BC%II,BC%JJ+1:BC%JJ+2,BC%KK)/) U_TEMP(1,1,1) = VV(BC%II,BC%JJ+1,BC%KK) - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER) + B_TEMP(1,1,1) = BFY(BC%II,BC%JJ+1,BC%KK) + CALL GET_SCALAR_FACE_VALUE_NEW(U_TEMP,Z_TEMP,F_TEMP,B_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER) FY_H_S(BC%II,BC%JJ+1,BC%KK) = F_TEMP(1,1,1) ENDIF CASE(-2) OFF_WALL_SELECT_1 IF ((VV(BC%II,BC%JJ-2,BC%KK)<0._EB) .AND. .NOT.(CELL(CELL_INDEX(BC%II,BC%JJ-1,BC%KK))%WALL_INDEX(-2)>0)) THEN Z_TEMP(1,1:3,1) = (/RHO_H_S_P(BC%II,BC%JJ-2:BC%JJ-1,BC%KK),RHO_H_S_P(BC%II,BC%JJ-1,BC%KK)/) U_TEMP(1,1,1) = VV(BC%II,BC%JJ-2,BC%KK) - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER) + B_TEMP(1,1,1) = BFY(BC%II,BC%JJ-2,BC%KK) + CALL GET_SCALAR_FACE_VALUE_NEW(U_TEMP,Z_TEMP,F_TEMP,B_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER) FY_H_S(BC%II,BC%JJ-2,BC%KK) = F_TEMP(1,1,1) ENDIF CASE( 3) OFF_WALL_SELECT_1 IF ((WW(BC%II,BC%JJ,BC%KK+1)>0._EB) .AND. .NOT.(CELL(CELL_INDEX(BC%II,BC%JJ,BC%KK+1))%WALL_INDEX(+3)>0)) THEN Z_TEMP(1,1,0:2) = (/RHO_H_S_P(BC%II,BC%JJ,BC%KK+1),RHO_H_S_P(BC%II,BC%JJ,BC%KK+1:BC%KK+2)/) U_TEMP(1,1,1) = WW(BC%II,BC%JJ,BC%KK+1) - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER) + B_TEMP(1,1,1) = BFZ(BC%II,BC%JJ,BC%KK+1) + CALL GET_SCALAR_FACE_VALUE_NEW(U_TEMP,Z_TEMP,F_TEMP,B_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER) FZ_H_S(BC%II,BC%JJ,BC%KK+1) = F_TEMP(1,1,1) ENDIF CASE(-3) OFF_WALL_SELECT_1 IF ((WW(BC%II,BC%JJ,BC%KK-2)<0._EB) .AND. .NOT.(CELL(CELL_INDEX(BC%II,BC%JJ,BC%KK-1))%WALL_INDEX(-3)>0)) THEN Z_TEMP(1,1,1:3) = (/RHO_H_S_P(BC%II,BC%JJ,BC%KK-2:BC%KK-1),RHO_H_S_P(BC%II,BC%JJ,BC%KK-1)/) U_TEMP(1,1,1) = WW(BC%II,BC%JJ,BC%KK-2) - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER) + B_TEMP(1,1,1) = BFZ(BC%II,BC%JJ,BC%KK-2) + CALL GET_SCALAR_FACE_VALUE_NEW(U_TEMP,Z_TEMP,F_TEMP,B_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER) FZ_H_S(BC%II,BC%JJ,BC%KK-2) = F_TEMP(1,1,1) ENDIF END SELECT OFF_WALL_SELECT_1 @@ -966,117 +976,7 @@ SUBROUTINE ENTHALPY_ADVECTION(U_DOT_DEL_RHO_H_S) ENDDO !$OMP END PARALLEL DO -END SUBROUTINE ENTHALPY_ADVECTION - - -SUBROUTINE SPECIES_ADVECTION_PART_1 - -USE MATH_FUNCTIONS, ONLY: GET_SCALAR_FACE_VALUE -USE PHYSICAL_FUNCTIONS, ONLY: GET_MOLECULAR_WEIGHT -REAL(EB), POINTER, DIMENSION(:,:,:) :: RHO_Z_P -REAL(EB), POINTER, DIMENSION(:,:,:,:) :: FX_ZZ,FY_ZZ,FZ_ZZ -REAL(EB), DIMENSION(0:3,0:3,0:3) :: U_TEMP,F_TEMP -REAL(EB), DIMENSION(-1:3,-1:3,-1:3) :: Z_TEMP -INTEGER :: I,J,K,IW,N -TYPE(WALL_TYPE), POINTER :: WC -TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC -TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1 - -FX_ZZ=>SWORK1 -FY_ZZ=>SWORK2 -FZ_ZZ=>SWORK3 -RHO_Z_P=>WORK_PAD - -! Species face values - -SPECIES_LOOP: DO N=1,N_TOTAL_SCALARS - - !$OMP PARALLEL DO - DO K=-1,KBP1+1 - DO J=-1,JBP1+1 - DO I=-1,IBP1+1 - RHO_Z_P(I,J,K) = RHOP(I,J,K)*ZZP(I,J,K,N) - ENDDO - ENDDO - ENDDO - !$OMP END PARALLEL DO - - ! Compute scalar face values - - CALL GET_SCALAR_FACE_VALUE(UU,RHO_Z_P,FX_ZZ(:,:,:,N),0,IBAR,1,JBAR,1,KBAR,1,I_FLUX_LIMITER) - CALL GET_SCALAR_FACE_VALUE(VV,RHO_Z_P,FY_ZZ(:,:,:,N),1,IBAR,0,JBAR,1,KBAR,2,I_FLUX_LIMITER) - CALL GET_SCALAR_FACE_VALUE(WW,RHO_Z_P,FZ_ZZ(:,:,:,N),1,IBAR,1,JBAR,0,KBAR,3,I_FLUX_LIMITER) - - !$OMP PARALLEL DO PRIVATE(IW,WC,BC,B1,U_TEMP,Z_TEMP,F_TEMP) - WALL_LOOP_2: DO IW=1,N_EXTERNAL_WALL_CELLS+N_INTERNAL_WALL_CELLS - WC=>WALL(IW) - IF (WC%BOUNDARY_TYPE==NULL_BOUNDARY) CYCLE WALL_LOOP_2 - BC=>BOUNDARY_COORD(WC%BC_INDEX) - B1=>BOUNDARY_PROP1(WC%B1_INDEX) - - ! 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_2: IF (WC%BOUNDARY_TYPE/=INTERPOLATED_BOUNDARY .AND. WC%BOUNDARY_TYPE/=OPEN_BOUNDARY) THEN - - OFF_WALL_SELECT_2: SELECT CASE(BC%IOR) - CASE( 1) OFF_WALL_SELECT_2 - ! ghost FX/UU(II+1) - ! /// II /// II+1 | II+2 | ... - ! ^ WALL_INDEX(II+1,+1) - IF ((UU(BC%II+1,BC%JJ,BC%KK)>0._EB) .AND. .NOT.(CELL(CELL_INDEX(BC%II+1,BC%JJ,BC%KK))%WALL_INDEX(+1)>0)) THEN - Z_TEMP(0:2,1,1) = (/RHO_Z_P(BC%II+1,BC%JJ,BC%KK),RHO_Z_P(BC%II+1:BC%II+2,BC%JJ,BC%KK)/) - U_TEMP(1,1,1) = UU(BC%II+1,BC%JJ,BC%KK) - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER) - FX_ZZ(BC%II+1,BC%JJ,BC%KK,N) = F_TEMP(1,1,1) - ENDIF - CASE(-1) OFF_WALL_SELECT_2 - ! FX/UU(II-2) ghost - ! ... | II-2 | II-1 /// II /// - ! ^ WALL_INDEX(II-1,-1) - IF ((UU(BC%II-2,BC%JJ,BC%KK)<0._EB) .AND. .NOT.(CELL(CELL_INDEX(BC%II-1,BC%JJ,BC%KK))%WALL_INDEX(-1)>0)) THEN - Z_TEMP(1:3,1,1) = (/RHO_Z_P(BC%II-2:BC%II-1,BC%JJ,BC%KK),RHO_Z_P(BC%II-1,BC%JJ,BC%KK)/) - U_TEMP(1,1,1) = UU(BC%II-2,BC%JJ,BC%KK) - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER) - FX_ZZ(BC%II-2,BC%JJ,BC%KK,N) = F_TEMP(1,1,1) - ENDIF - CASE( 2) OFF_WALL_SELECT_2 - IF ((VV(BC%II,BC%JJ+1,BC%KK)>0._EB) .AND. .NOT.(CELL(CELL_INDEX(BC%II,BC%JJ+1,BC%KK))%WALL_INDEX(+2)>0)) THEN - Z_TEMP(1,0:2,1) = (/RHO_Z_P(BC%II,BC%JJ+1,BC%KK),RHO_Z_P(BC%II,BC%JJ+1:BC%JJ+2,BC%KK)/) - U_TEMP(1,1,1) = VV(BC%II,BC%JJ+1,BC%KK) - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER) - FY_ZZ(BC%II,BC%JJ+1,BC%KK,N) = F_TEMP(1,1,1) - ENDIF - CASE(-2) OFF_WALL_SELECT_2 - IF ((VV(BC%II,BC%JJ-2,BC%KK)<0._EB) .AND. .NOT.(CELL(CELL_INDEX(BC%II,BC%JJ-1,BC%KK))%WALL_INDEX(-2)>0)) THEN - Z_TEMP(1,1:3,1) = (/RHO_Z_P(BC%II,BC%JJ-2:BC%JJ-1,BC%KK),RHO_Z_P(BC%II,BC%JJ-1,BC%KK)/) - U_TEMP(1,1,1) = VV(BC%II,BC%JJ-2,BC%KK) - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER) - FY_ZZ(BC%II,BC%JJ-2,BC%KK,N) = F_TEMP(1,1,1) - ENDIF - CASE( 3) OFF_WALL_SELECT_2 - IF ((WW(BC%II,BC%JJ,BC%KK+1)>0._EB) .AND. .NOT.(CELL(CELL_INDEX(BC%II,BC%JJ,BC%KK+1))%WALL_INDEX(+3)>0)) THEN - Z_TEMP(1,1,0:2) = (/RHO_Z_P(BC%II,BC%JJ,BC%KK+1),RHO_Z_P(BC%II,BC%JJ,BC%KK+1:BC%KK+2)/) - U_TEMP(1,1,1) = WW(BC%II,BC%JJ,BC%KK+1) - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER) - FZ_ZZ(BC%II,BC%JJ,BC%KK+1,N) = F_TEMP(1,1,1) - ENDIF - CASE(-3) OFF_WALL_SELECT_2 - IF ((WW(BC%II,BC%JJ,BC%KK-2)<0._EB) .AND. .NOT.(CELL(CELL_INDEX(BC%II,BC%JJ,BC%KK-1))%WALL_INDEX(-3)>0)) THEN - Z_TEMP(1,1,1:3) = (/RHO_Z_P(BC%II,BC%JJ,BC%KK-2:BC%KK-1),RHO_Z_P(BC%II,BC%JJ,BC%KK-1)/) - U_TEMP(1,1,1) = WW(BC%II,BC%JJ,BC%KK-2) - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER) - FZ_ZZ(BC%II,BC%JJ,BC%KK-2,N) = F_TEMP(1,1,1) - ENDIF - END SELECT OFF_WALL_SELECT_2 - - ENDIF OFF_WALL_IF_2 - - ENDDO WALL_LOOP_2 - !$OMP END PARALLEL DO - -ENDDO SPECIES_LOOP - -END SUBROUTINE SPECIES_ADVECTION_PART_1 +END SUBROUTINE ENTHALPY_ADVECTION_NEW SUBROUTINE SPECIES_ADVECTION_PART_1_NEW diff --git a/Source/init.f90 b/Source/init.f90 index bfc0108a179..aa4ddfb6ead 100644 --- a/Source/init.f90 +++ b/Source/init.f90 @@ -595,18 +595,6 @@ SUBROUTINE INITIALIZE_MESH_VARIABLES_1(DT,NM) ALLOCATE(M%RSUM(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO) ; CALL ChkMemErr('INIT','RSUM',IZERO) ; M%RSUM=RSUM0 -! Allocate scalar face values - -ALLOCATE( M%FX(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO); CALL ChkMemErr('INIT','FX',IZERO) ; M%FX=0._EB -ALLOCATE( M%FY(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO); CALL ChkMemErr('INIT','FY',IZERO) ; M%FY=0._EB -ALLOCATE( M%FZ(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO); CALL ChkMemErr('INIT','FZ',IZERO) ; M%FZ=0._EB - -IF (FLUX_LIMITER_SINGLE_COEF) THEN - ALLOCATE( M%BFX(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO); CALL ChkMemErr('INIT','BFX',IZERO) ; M%BFX=0._EB - ALLOCATE( M%BFY(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO); CALL ChkMemErr('INIT','BFY',IZERO) ; M%BFY=0._EB - ALLOCATE( M%BFZ(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO); CALL ChkMemErr('INIT','BFZ',IZERO) ; M%BFZ=0._EB -ENDIF - ! Allocate storage for scalar total fluxes IF (STORE_SPECIES_FLUX) THEN @@ -831,6 +819,23 @@ SUBROUTINE INITIALIZE_MESH_VARIABLES_1(DT,NM) M%ZZS(:,:,:,N) = INITIAL_UNMIXED_FRACTION ENDDO +! Allocate scalar face values + +ALLOCATE( M%FX(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO); CALL ChkMemErr('INIT','FX',IZERO) ; M%FX=0._EB +ALLOCATE( M%FY(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO); CALL ChkMemErr('INIT','FY',IZERO) ; M%FY=0._EB +ALLOCATE( M%FZ(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO); CALL ChkMemErr('INIT','FZ',IZERO) ; M%FZ=0._EB +DO N=1,N_TRACKED_SPECIES + M%FX(:,:,:,N)=RHOA*SPECIES_MIXTURE(N)%ZZ0 + M%FY(:,:,:,N)=RHOA*SPECIES_MIXTURE(N)%ZZ0 + M%FZ(:,:,:,N)=RHOA*SPECIES_MIXTURE(N)%ZZ0 +ENDDO + +! Allocate flux limiter coefficients + +ALLOCATE( M%BFX(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO); CALL ChkMemErr('INIT','BFX',IZERO) ; M%BFX=0._EB +ALLOCATE( M%BFY(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO); CALL ChkMemErr('INIT','BFY',IZERO) ; M%BFY=0._EB +ALLOCATE( M%BFZ(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO); CALL ChkMemErr('INIT','BFZ',IZERO) ; M%BFZ=0._EB + ! Allocate and Initialize Mesh-Dependent Radiation Arrays M%QR = 0._EB diff --git a/Source/main.f90 b/Source/main.f90 index d7a282c6424..4b4968089a8 100644 --- a/Source/main.f90 +++ b/Source/main.f90 @@ -626,11 +626,7 @@ PROGRAM FDS COMPUTE_FINITE_DIFFERENCES_1: DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX CALL INSERT_ALL_PARTICLES(T,NM) IF (.NOT.SOLID_PHASE_ONLY .AND. .NOT.FREEZE_VELOCITY) CALL COMPUTE_VISCOSITY(NM,APPLY_TO_ESTIMATED_VARIABLES=.FALSE.) - IF (FLUX_LIMITER_SINGLE_COEF) THEN - CALL MASS_FINITE_DIFFERENCES_NEW(NM) - ELSE - CALL MASS_FINITE_DIFFERENCES(NM) - ENDIF + CALL MASS_FINITE_DIFFERENCES_NEW(NM) ENDDO COMPUTE_FINITE_DIFFERENCES_1 ! Estimate quantities at next time step, and decrease/increase time step if necessary based on CFL condition @@ -811,11 +807,7 @@ PROGRAM FDS COMPUTE_FINITE_DIFFERENCES_2: DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX IF (.NOT.SOLID_PHASE_ONLY .AND. .NOT.FREEZE_VELOCITY) CALL COMPUTE_VISCOSITY(NM,APPLY_TO_ESTIMATED_VARIABLES=.TRUE.) - IF (FLUX_LIMITER_SINGLE_COEF) THEN - CALL MASS_FINITE_DIFFERENCES_NEW(NM) - ELSE - CALL MASS_FINITE_DIFFERENCES(NM) - ENDIF + CALL MASS_FINITE_DIFFERENCES_NEW(NM) CALL DENSITY(T,DT,NM) IF (LEVEL_SET_MODE>0) CALL LEVEL_SET_FIRESPREAD(T,DT,NM) ENDDO COMPUTE_FINITE_DIFFERENCES_2 diff --git a/Source/mass.f90 b/Source/mass.f90 index dfb7dae2e4c..f92260f7728 100644 --- a/Source/mass.f90 +++ b/Source/mass.f90 @@ -9,7 +9,7 @@ MODULE MASS IMPLICIT NONE (TYPE,EXTERNAL) PRIVATE -PUBLIC MASS_FINITE_DIFFERENCES,MASS_FINITE_DIFFERENCES_NEW,DENSITY +PUBLIC MASS_FINITE_DIFFERENCES_NEW,DENSITY CONTAINS @@ -17,172 +17,6 @@ MODULE MASS !> \brief Compute spatial differences for mass transport equations !> \param NM Mesh index -SUBROUTINE MASS_FINITE_DIFFERENCES(NM) - -USE MATH_FUNCTIONS, ONLY: GET_SCALAR_FACE_VALUE -USE PHYSICAL_FUNCTIONS, ONLY: GET_MOLECULAR_WEIGHT - -INTEGER, INTENT(IN) :: NM -REAL(EB) :: TNOW -REAL(EB), DIMENSION(0:3,0:3,0:3) :: U_TEMP,F_TEMP -REAL(EB), DIMENSION(-1:3,-1:3,-1:3) :: Z_TEMP -REAL(EB), PARAMETER :: DUMMY=0._EB -INTEGER :: I,J,K,N,IOR,IW,IIG,JJG,KKG,II,JJ,KK,IC -REAL(EB), POINTER, DIMENSION(:,:,:) :: UU,VV,WW,RHOP -REAL(EB), POINTER, DIMENSION(:,:,:,:) :: ZZP -REAL(EB), POINTER, DIMENSION(:,:,:) :: RHO_Z_P -TYPE(WALL_TYPE), POINTER :: WC -TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC -TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1 - -IF (SOLID_PHASE_ONLY) RETURN - -TNOW=CURRENT_TIME() -CALL POINT_TO_MESH(NM) - -IF (PREDICTOR) THEN - UU => U - VV => V - WW => W - RHOP => RHO - ZZP => ZZ -ELSE - UU => US - VV => VS - WW => WS - RHOP => RHOS - ZZP => ZZS -ENDIF - -! Reset counter for CLIP_RHOMIN, CLIP_RHOMAX -! Done here so DT_RESTRICT_COUNT will persist until WRITE_DIAGNOSTICS is called - -IF (PREDICTOR) DT_RESTRICT_COUNT = 0 - -! Species face values - -SPECIES_LOOP: DO N=1,N_TOTAL_SCALARS - - RHO_Z_P=>WORK_PAD - - !$OMP PARALLEL DO - DO K=-1,KBP1+1 - DO J=-1,JBP1+1 - DO I=-1,IBP1+1 - RHO_Z_P(I,J,K) = RHOP(I,J,K)*ZZP(I,J,K,N) - ENDDO - ENDDO - ENDDO - !$OMP END PARALLEL DO - - ! Compute scalar face values - - CALL GET_SCALAR_FACE_VALUE(UU,RHO_Z_P,FX(:,:,:,N),0,IBAR,1,JBAR,1,KBAR,1,I_FLUX_LIMITER) - CALL GET_SCALAR_FACE_VALUE(VV,RHO_Z_P,FY(:,:,:,N),1,IBAR,0,JBAR,1,KBAR,2,I_FLUX_LIMITER) - CALL GET_SCALAR_FACE_VALUE(WW,RHO_Z_P,FZ(:,:,:,N),1,IBAR,1,JBAR,0,KBAR,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_2: DO IW=1,N_EXTERNAL_WALL_CELLS+N_INTERNAL_WALL_CELLS - WC=>WALL(IW) - IF (WC%BOUNDARY_TYPE==NULL_BOUNDARY) CYCLE WALL_LOOP_2 - 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,N) = 0._EB - CASE(-1); FX(IIG,JJG,KKG,N) = 0._EB - CASE( 2); FY(IIG,JJG-1,KKG,N) = 0._EB - CASE(-2); FY(IIG,JJG,KKG,N) = 0._EB - CASE( 3); FZ(IIG,JJG,KKG-1,N) = 0._EB - CASE(-3); FZ(IIG,JJG,KKG,N) = 0._EB - END SELECT - ELSE - SELECT CASE(IOR) - CASE( 1); FX(IIG-1,JJG,KKG,N) = B1%RHO_F*B1%ZZ_F(N) - CASE(-1); FX(IIG,JJG,KKG,N) = B1%RHO_F*B1%ZZ_F(N) - CASE( 2); FY(IIG,JJG-1,KKG,N) = B1%RHO_F*B1%ZZ_F(N) - CASE(-2); FY(IIG,JJG,KKG,N) = B1%RHO_F*B1%ZZ_F(N) - CASE( 3); FZ(IIG,JJG,KKG-1,N) = B1%RHO_F*B1%ZZ_F(N) - CASE(-3); FZ(IIG,JJG,KKG,N) = B1%RHO_F*B1%ZZ_F(N) - 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_2: IF (WC%BOUNDARY_TYPE/=INTERPOLATED_BOUNDARY .AND. WC%BOUNDARY_TYPE/=OPEN_BOUNDARY) THEN - - OFF_WALL_SELECT_2: SELECT CASE(IOR) - CASE( 1) OFF_WALL_SELECT_2 - ! 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) = (/RHO_Z_P(II+1,JJ,KK),RHO_Z_P(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,N) = F_TEMP(1,1,1) - ENDIF - CASE(-1) OFF_WALL_SELECT_2 - ! 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,RHO_Z_P(II-2:II-1,JJ,KK),RHO_Z_P(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,N) = F_TEMP(1,1,1) - ENDIF - CASE( 2) OFF_WALL_SELECT_2 - 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) = (/RHO_Z_P(II,JJ+1,KK),RHO_Z_P(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,N) = F_TEMP(1,1,1) - ENDIF - CASE(-2) OFF_WALL_SELECT_2 - 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,RHO_Z_P(II,JJ-2:JJ-1,KK),RHO_Z_P(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,N) = F_TEMP(1,1,1) - ENDIF - CASE( 3) OFF_WALL_SELECT_2 - 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) = (/RHO_Z_P(II,JJ,KK+1),RHO_Z_P(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,N) = F_TEMP(1,1,1) - ENDIF - CASE(-3) OFF_WALL_SELECT_2 - 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,RHO_Z_P(II,JJ,KK-2:KK-1),RHO_Z_P(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,N) = F_TEMP(1,1,1) - ENDIF - END SELECT OFF_WALL_SELECT_2 - - ENDIF OFF_WALL_IF_2 - - ENDDO WALL_LOOP_2 - !$OMP END PARALLEL DO - -ENDDO SPECIES_LOOP - -T_USED(3)=T_USED(3)+CURRENT_TIME()-TNOW -END SUBROUTINE MASS_FINITE_DIFFERENCES - - SUBROUTINE MASS_FINITE_DIFFERENCES_NEW(NM) USE MATH_FUNCTIONS, ONLY: GET_SCALAR_FACE_VALUE_NEW,GET_SCALAR_FACE_COEF @@ -294,6 +128,7 @@ SUBROUTINE MASS_FINITE_DIFFERENCES_NEW(NM) IC = CELL_INDEX(II,JJ,KK) IF (WC%BOUNDARY_TYPE==SOLID_BOUNDARY .AND. .NOT.CELL(IC)%SOLID .AND. .NOT.CELL(IC)%EXTERIOR) THEN + ! thin obstruction SELECT CASE(IOR) CASE( 1); FX(IIG-1,JJG,KKG,N) = 0._EB CASE(-1); FX(IIG,JJG,KKG,N) = 0._EB @@ -303,14 +138,16 @@ SUBROUTINE MASS_FINITE_DIFFERENCES_NEW(NM) CASE(-3); FZ(IIG,JJG,KKG,N) = 0._EB END SELECT ELSE - SELECT CASE(IOR) - CASE( 1); FX(IIG-1,JJG,KKG,N) = B1%RHO_F*B1%ZZ_F(N) - CASE(-1); FX(IIG,JJG,KKG,N) = B1%RHO_F*B1%ZZ_F(N) - CASE( 2); FY(IIG,JJG-1,KKG,N) = B1%RHO_F*B1%ZZ_F(N) - CASE(-2); FY(IIG,JJG,KKG,N) = B1%RHO_F*B1%ZZ_F(N) - CASE( 3); FZ(IIG,JJG,KKG-1,N) = B1%RHO_F*B1%ZZ_F(N) - CASE(-3); FZ(IIG,JJG,KKG,N) = B1%RHO_F*B1%ZZ_F(N) - END SELECT + IF (WC%BOUNDARY_TYPE/=INTERPOLATED_BOUNDARY) THEN + SELECT CASE(IOR) + CASE( 1); FX(IIG-1,JJG,KKG,N) = B1%RHO_F*B1%ZZ_F(N) + CASE(-1); FX(IIG,JJG,KKG,N) = B1%RHO_F*B1%ZZ_F(N) + CASE( 2); FY(IIG,JJG-1,KKG,N) = B1%RHO_F*B1%ZZ_F(N) + CASE(-2); FY(IIG,JJG,KKG,N) = B1%RHO_F*B1%ZZ_F(N) + CASE( 3); FZ(IIG,JJG,KKG-1,N) = B1%RHO_F*B1%ZZ_F(N) + CASE(-3); FZ(IIG,JJG,KKG,N) = B1%RHO_F*B1%ZZ_F(N) + END SELECT + ENDIF ENDIF ! Overwrite first off-wall advective flux if flow is away from the wall and if the face is not also a wall cell diff --git a/Source/read.f90 b/Source/read.f90 index de2085573d9..6c8070cd076 100644 --- a/Source/read.f90 +++ b/Source/read.f90 @@ -1739,7 +1739,7 @@ SUBROUTINE READ_MISC CFL_MAX,CFL_MIN,CFL_VELOCITY_NORM,CHECK_FO,CHECK_HT,CHECK_VN, & CNF_CUTOFF,CONSTANT_SPECIFIC_HEAT_RATIO,& C_SMAGORINSKY,C_VREMAN,C_WALE,DEPOSITION,EXTERNAL_FILENAME,& - FIXED_LES_FILTER_WIDTH,FLUX_LIMITER,FLUX_LIMITER_SINGLE_COEF,& + FIXED_LES_FILTER_WIDTH,FLUX_LIMITER,& FREEZE_VELOCITY,FYI,GAMMA,GRAVITATIONAL_DEPOSITION,& GRAVITATIONAL_SETTLING,GVEC,H_F_REFERENCE_TEMPERATURE,& HUMIDITY,HVAC_LOCAL_PRESSURE,HVAC_MASS_TRANSPORT_CELL_L,HVAC_PRES_RELAX,HVAC_QFAN,IBLANK_SMV,I_MAX_TEMP,& diff --git a/Source/wall.f90 b/Source/wall.f90 index 884188ed5da..ec8b4b2183d 100644 --- a/Source/wall.f90 +++ b/Source/wall.f90 @@ -270,13 +270,16 @@ END SUBROUTINE WALL_BC SUBROUTINE ASSIGN_GHOST_VALUE(WALL_INDEX,BC,B1) USE PHYSICAL_FUNCTIONS, ONLY : GET_SPECIFIC_GAS_CONSTANT +USE COMPLEX_GEOMETRY, ONLY : CC_CGSC, CC_SOLID INTEGER, INTENT(IN) :: WALL_INDEX REAL(EB) :: ARO,ZZ_GET(1:N_TRACKED_SPECIES),RHO_OTHER,RHO_ZZ_OTHER(1:N_TOTAL_SCALARS),& RSUM_TMP,RHO_OTHER_2,RHO_ZZ_OTHER_2(1:N_TOTAL_SCALARS) -INTEGER :: IIO,JJO,KKO,II2,JJ2,KK2 +INTEGER :: IIO,JJO,KKO,II2,JJ2,KK2,ICG,ICO +LOGICAL :: CC_SOLID_FLAG,SECOND_ORDER_INTERPOLATED_BOUNDARY REAL(EB), POINTER, DIMENSION(:,:,:,:) :: OM_ZZP REAL(EB), POINTER, DIMENSION(:,:,:) :: OM_RHOP TYPE(EXTERNAL_WALL_TYPE), POINTER :: EWC +TYPE(WALL_TYPE), POINTER :: WC TYPE(OMESH_TYPE), POINTER :: OM TYPE(MESH_TYPE), POINTER :: MM TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC @@ -331,13 +334,43 @@ SUBROUTINE ASSIGN_GHOST_VALUE(WALL_INDEX,BC,B1) CALL GET_SPECIFIC_GAS_CONSTANT(ZZ_GET,RSUM(BC%II,BC%JJ,BC%KK)) TMP(BC%II,BC%JJ,BC%KK) = PBAR_P(BC%KK,B1%PRESSURE_ZONE)/(RSUM(BC%II,BC%JJ,BC%KK)*RHOP(BC%II,BC%JJ,BC%KK)) -! Assign RHOP and ZZP to the second ghost cell +! Assign RHOP and ZZP to the second ghost cell (if it is not a solid) -RHOP(BC%II2,BC%JJ2,BC%KK2) = RHO_OTHER_2 -ZZP(BC%II2,BC%JJ2,BC%KK2,1:N_TOTAL_SCALARS) = MAX(0._EB,MIN(1._EB,RHO_ZZ_OTHER_2(1:N_TOTAL_SCALARS)/RHO_OTHER_2)) -ZZ_GET(1:N_TRACKED_SPECIES) = ZZP(BC%II2,BC%JJ2,BC%KK2,1:N_TRACKED_SPECIES) -CALL GET_SPECIFIC_GAS_CONSTANT(ZZ_GET,RSUM_TMP) -TMP(BC%II2,BC%JJ2,BC%KK2) = PBAR_P(BC%KK,B1%PRESSURE_ZONE)/(RSUM_TMP*RHOP(BC%II2,BC%JJ2,BC%KK2)) +SECOND_ORDER_INTERPOLATED_BOUNDARY = .FALSE. +WC => WALL(WALL_INDEX) +INTERPOLATED_IF: IF (WC%BOUNDARY_TYPE==INTERPOLATED_BOUNDARY) THEN + ! Determine if there are 4 equally sized cells spanning the interpolated boundary + CC_SOLID_FLAG = .FALSE. + IF (ABS(EWC%AREA_RATIO-1._EB)<0.01_EB) THEN + IIO = EWC%IIO_MIN + JJO = EWC%JJO_MIN + KKO = EWC%KKO_MIN + SELECT CASE(BC%IOR) + CASE( 1) ; ICG = CELL_INDEX(BC%IIG+1,BC%JJG,BC%KKG) ; ICO = MM%CELL_INDEX(IIO-1,JJO,KKO) + CASE(-1) ; ICG = CELL_INDEX(BC%IIG-1,BC%JJG,BC%KKG) ; ICO = MM%CELL_INDEX(IIO+1,JJO,KKO) + CASE( 2) ; ICG = CELL_INDEX(BC%IIG,BC%JJG+1,BC%KKG) ; ICO = MM%CELL_INDEX(IIO,JJO-1,KKO) + CASE(-2) ; ICG = CELL_INDEX(BC%IIG,BC%JJG-1,BC%KKG) ; ICO = MM%CELL_INDEX(IIO,JJO+1,KKO) + CASE( 3) ; ICG = CELL_INDEX(BC%IIG,BC%JJG,BC%KKG+1) ; ICO = MM%CELL_INDEX(IIO,JJO,KKO-1) + CASE(-3) ; ICG = CELL_INDEX(BC%IIG,BC%JJG,BC%KKG-1) ; ICO = MM%CELL_INDEX(IIO,JJO,KKO+1) + END SELECT + IF (CC_IBM) THEN ! Test if one of surrounding cells is CC_SOLID. + IF(CCVAR(BC%IIG,BC%JJG,BC%KKG,CC_CGSC)==CC_SOLID .OR. CCVAR(BC%II,BC%JJ,BC%KK,CC_CGSC)==CC_SOLID) CC_SOLID_FLAG = .TRUE. + ENDIF + IF (.NOT.CELL(ICG)%SOLID.AND..NOT.MM%CELL(ICO)%SOLID.AND..NOT.CC_SOLID_FLAG) SECOND_ORDER_INTERPOLATED_BOUNDARY = .TRUE. + ENDIF +ENDIF INTERPOLATED_IF + +IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN + RHOP(BC%II2,BC%JJ2,BC%KK2) = RHO_OTHER_2 + ZZP(BC%II2,BC%JJ2,BC%KK2,1:N_TOTAL_SCALARS) = MAX(0._EB,MIN(1._EB,RHO_ZZ_OTHER_2(1:N_TOTAL_SCALARS)/RHO_OTHER_2)) + ZZ_GET(1:N_TRACKED_SPECIES) = ZZP(BC%II2,BC%JJ2,BC%KK2,1:N_TRACKED_SPECIES) + CALL GET_SPECIFIC_GAS_CONSTANT(ZZ_GET,RSUM_TMP) + TMP(BC%II2,BC%JJ2,BC%KK2) = PBAR_P(BC%KK2,B1%PRESSURE_ZONE)/(RSUM_TMP*RHOP(BC%II2,BC%JJ2,BC%KK2)) +ELSE + RHOP(BC%II2,BC%JJ2,BC%KK2) = RHOP(BC%II,BC%JJ,BC%KK) + ZZP(BC%II2,BC%JJ2,BC%KK2,1:N_TOTAL_SCALARS) = ZZP(BC%II,BC%JJ,BC%KK,1:N_TOTAL_SCALARS) + TMP(BC%II2,BC%JJ2,BC%KK2) = TMP(BC%II,BC%JJ,BC%KK) +ENDIF END SUBROUTINE ASSIGN_GHOST_VALUE @@ -447,25 +480,21 @@ END SUBROUTINE NEAR_SURFACE_GAS_VARIABLES SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_INDEX) USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP,INTERPOLATE1D_UNIFORM,GET_SCALAR_FACE_VALUE,GET_SCALAR_FACE_VALUE_NEW,GET_SCALAR_FACE_COEF -USE COMPLEX_GEOMETRY, ONLY : CC_CGSC, CC_SOLID USE PHYSICAL_FUNCTIONS, ONLY : GET_SPECIFIC_GAS_CONSTANT,GET_VISCOSITY,GET_MOLECULAR_WEIGHT USE DEVICE_VARIABLES, ONLY : PROPERTY,PROPERTY_TYPE REAL(EB), INTENT(IN) :: T INTEGER, INTENT(IN) :: NM INTEGER, INTENT(IN), OPTIONAL :: WALL_INDEX,CFACE_INDEX,PARTICLE_INDEX -REAL(EB) :: ARO,QNET,RAMP_FACTOR,RHO_G_2,RSUM_F,PBAR_F,TSI,UN, & +REAL(EB) :: ARO,QNET,RAMP_FACTOR,RSUM_F,PBAR_F,TSI,UN, & RHO_ZZ_F(1:N_TOTAL_SCALARS),ZZ_GET(1:N_TRACKED_SPECIES),& - RHO_OTHER,RHO_OTHER_2,RHO_ZZ_OTHER(1:N_TOTAL_SCALARS),RHO_ZZ_OTHER_2,RHO_ZZ_G,RHO_ZZ_G_2, & - DDO,PBAR_G,PBAR_OTHER,DENOM,D_Z_N(0:I_MAX_TEMP),D_Z_G,D_Z_OTHER,TMP_OTHER, & - MU_DNS_G,MU_DNS_OTHER,MU_OTHER,RHO_D,RHO_D_TURB,RHO_D_DZDN,RHO_D_DZDN_OTHER,RSUM_OTHER,& + RHO_OTHER, PBAR_OTHER,D_Z_N(0:I_MAX_TEMP),D_Z_G,D_Z_OTHER,TMP_OTHER,DDO,DENOM,PBAR_G, & + MU_DNS_G,MU_DNS_OTHER,MU_OTHER,RHO_D,RHO_D_TURB,RHO_D_DZDN,RHO_D_DZDN_OTHER,RSUM_OTHER, & BBB,CCC,PPP,QQQ,RRR,UUU,YYY,WWW,HTC_OLD,RSC_LOC,DTMP,RSUM_G,MU_G -LOGICAL :: SECOND_ORDER_INTERPOLATED_BOUNDARY,ATMOSPHERIC_INTERPOLATION,CC_SOLID_FLAG -INTEGER :: IIO,JJO,KKO,N,ADCOUNT,ICG,ICO +LOGICAL :: ATMOSPHERIC_INTERPOLATION +INTEGER :: IIO,JJO,KKO,N,ADCOUNT REAL(EB), POINTER, DIMENSION(:,:,:,:) :: OM_ZZP REAL(EB), POINTER, DIMENSION(:,:,:) :: OM_MUP -REAL(EB), DIMENSION(0:3,0:3,0:3) :: U_TEMP,F_TEMP,B_TEMP -REAL(EB), DIMENSION(-1:3,-1:3,-1:3) :: Z_TEMP TYPE(EXTERNAL_WALL_TYPE), POINTER :: EWC TYPE(OMESH_TYPE), POINTER :: OM TYPE(MESH_TYPE), POINTER :: MM @@ -667,44 +696,11 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I EWC => EXTERNAL_WALL(WALL_INDEX) OM => OMESH(EWC%NOM) - IF (PREDICTOR) THEN - OM_RHOP => OM%RHOS - OM_ZZP => OM%ZZS - ELSE - OM_RHOP => OM%RHO - OM_ZZP => OM%ZZ - ENDIF - MM => MESHES(EWC%NOM) - RHO_G_2 = B1%RHO_G ! first-order - RHO_OTHER = RHOP(BC%II,BC%JJ,BC%KK) - RHO_OTHER_2 = RHOP(BC%II,BC%JJ,BC%KK) - RHO_ZZ_OTHER(1:N_TOTAL_SCALARS) = ZZP(BC%II,BC%JJ,BC%KK,1:N_TOTAL_SCALARS)*RHO_OTHER - - ! Determine if there are 4 equally sized cells spanning the interpolated boundary - - SECOND_ORDER_INTERPOLATED_BOUNDARY = .FALSE. - CC_SOLID_FLAG = .FALSE. - IF (ABS(EWC%AREA_RATIO-1._EB)<0.01_EB) THEN - IIO = EWC%IIO_MIN - JJO = EWC%JJO_MIN - KKO = EWC%KKO_MIN - SELECT CASE(BC%IOR) - CASE( 1) ; ICG = CELL_INDEX(BC%IIG+1,BC%JJG,BC%KKG) ; ICO = MM%CELL_INDEX(IIO-1,JJO,KKO) - CASE(-1) ; ICG = CELL_INDEX(BC%IIG-1,BC%JJG,BC%KKG) ; ICO = MM%CELL_INDEX(IIO+1,JJO,KKO) - CASE( 2) ; ICG = CELL_INDEX(BC%IIG,BC%JJG+1,BC%KKG) ; ICO = MM%CELL_INDEX(IIO,JJO-1,KKO) - CASE(-2) ; ICG = CELL_INDEX(BC%IIG,BC%JJG-1,BC%KKG) ; ICO = MM%CELL_INDEX(IIO,JJO+1,KKO) - CASE( 3) ; ICG = CELL_INDEX(BC%IIG,BC%JJG,BC%KKG+1) ; ICO = MM%CELL_INDEX(IIO,JJO,KKO-1) - CASE(-3) ; ICG = CELL_INDEX(BC%IIG,BC%JJG,BC%KKG-1) ; ICO = MM%CELL_INDEX(IIO,JJO,KKO+1) - END SELECT - IF (CC_IBM) THEN ! Test if one of surrounding cells is CC_SOLID. - IF(CCVAR(BC%IIG,BC%JJG,BC%KKG,CC_CGSC)==CC_SOLID .OR. CCVAR(BC%II,BC%JJ,BC%KK,CC_CGSC)==CC_SOLID) CC_SOLID_FLAG = .TRUE. - ENDIF - IF (.NOT.CELL(ICG)%SOLID.AND..NOT.MM%CELL(ICO)%SOLID.AND..NOT.CC_SOLID_FLAG) SECOND_ORDER_INTERPOLATED_BOUNDARY = .TRUE. - ENDIF - ! Density + ! interp or extrap RHO_OTHER for jump in vertical grid resolution, linear in temperature to match heat flux in divg + RHO_OTHER = RHOP(BC%II,BC%JJ,BC%KK) ATMOSPHERIC_INTERPOLATION = .FALSE. DDO = 1._EB KKO = EWC%KKO_MIN @@ -716,7 +712,6 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I ATMOSPHERIC_INTERPOLATION = .TRUE. IF (ATMOSPHERIC_INTERPOLATION) THEN - ! interp or extrap RHO_OTHER for jump in vertical grid resolution, linear in temperature to match heat flux in divg PBAR_G = PBAR_P(BC%KKG,B1%PRESSURE_ZONE) PBAR_OTHER = EVALUATE_RAMP(MM%ZC(EWC%KKO_MIN),I_RAMP_P0_Z) ZZ_GET(1:N_TRACKED_SPECIES) = MAX(0._EB,MIN(1._EB,ZZP(BC%II,BC%JJ,BC%KK,1:N_TOTAL_SCALARS))) @@ -726,199 +721,9 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I TMP(BC%II,BC%JJ,BC%KK) = PBAR_P(BC%KK,B1%PRESSURE_ZONE)/(RSUM(BC%II,BC%JJ,BC%KK)*RHOP(BC%II,BC%JJ,BC%KK)) ENDIF - SELECT CASE(BC%IOR) - CASE( 1) - IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN - RHO_G_2 = RHOP(BC%IIG+1,BC%JJG,BC%KKG) - RHO_OTHER_2 = RHOP(BC%II-1,BC%JJ,BC%KK) - ENDIF - Z_TEMP(0:3,1,1) = (/RHO_OTHER_2,RHO_OTHER,B1%RHO_G,RHO_G_2/) - U_TEMP(1,1,1) = UU(BC%II,BC%JJ,BC%KK) - IF (FLUX_LIMITER_SINGLE_COEF) THEN - B_TEMP(1,1,1) = BFX(BC%II,BC%JJ,BC%KK) - CALL GET_SCALAR_FACE_VALUE_NEW(U_TEMP,Z_TEMP,F_TEMP,B_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER) - ELSE - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER) - ENDIF - B1%RHO_F = F_TEMP(1,1,1) - CASE(-1) - IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN - RHO_G_2 = RHOP(BC%IIG-1,BC%JJG,BC%KKG) - RHO_OTHER_2 = RHOP(BC%II+1,BC%JJ,BC%KK) - ENDIF - Z_TEMP(0:3,1,1) = (/RHO_G_2,B1%RHO_G,RHO_OTHER,RHO_OTHER_2/) - U_TEMP(1,1,1) = UU(BC%II-1,BC%JJ,BC%KK) - IF (FLUX_LIMITER_SINGLE_COEF) THEN - B_TEMP(1,1,1) = BFX(BC%II-1,BC%JJ,BC%KK) - CALL GET_SCALAR_FACE_VALUE_NEW(U_TEMP,Z_TEMP,F_TEMP,B_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER) - ELSE - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER) - ENDIF - B1%RHO_F = F_TEMP(1,1,1) - CASE( 2) - IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN - RHO_G_2 = RHOP(BC%IIG,BC%JJG+1,BC%KKG) - RHO_OTHER_2 = RHOP(BC%II,BC%JJ-1,BC%KK) - ENDIF - Z_TEMP(1,0:3,1) = (/RHO_OTHER_2,RHO_OTHER,B1%RHO_G,RHO_G_2/) - U_TEMP(1,1,1) = VV(BC%II,BC%JJ,BC%KK) - IF (FLUX_LIMITER_SINGLE_COEF) THEN - B_TEMP(1,1,1) = BFY(BC%II,BC%JJ,BC%KK) - CALL GET_SCALAR_FACE_VALUE_NEW(U_TEMP,Z_TEMP,F_TEMP,B_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER) - ELSE - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER) - ENDIF - B1%RHO_F = F_TEMP(1,1,1) - CASE(-2) - IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN - RHO_G_2 = RHOP(BC%IIG,BC%JJG-1,BC%KKG) - RHO_OTHER_2 = RHOP(BC%II,BC%JJ+1,BC%KK) - ENDIF - Z_TEMP(1,0:3,1) = (/RHO_G_2,B1%RHO_G,RHO_OTHER,RHO_OTHER_2/) - U_TEMP(1,1,1) = VV(BC%II,BC%JJ-1,BC%KK) - IF (FLUX_LIMITER_SINGLE_COEF) THEN - B_TEMP(1,1,1) = BFY(BC%II,BC%JJ-1,BC%KK) - CALL GET_SCALAR_FACE_VALUE_NEW(U_TEMP,Z_TEMP,F_TEMP,B_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER) - ELSE - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER) - ENDIF - B1%RHO_F = F_TEMP(1,1,1) - CASE( 3) - IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN - RHO_G_2 = RHOP(BC%IIG,BC%JJG,BC%KKG+1) - RHO_OTHER_2 = RHOP(BC%II,BC%JJ,BC%KK-1) - ENDIF - Z_TEMP(1,1,0:3) = (/RHO_OTHER_2,RHO_OTHER,B1%RHO_G,RHO_G_2/) - U_TEMP(1,1,1) = WW(BC%II,BC%JJ,BC%KK) - IF (FLUX_LIMITER_SINGLE_COEF) THEN - B_TEMP(1,1,1) = BFZ(BC%II,BC%JJ,BC%KK) - CALL GET_SCALAR_FACE_VALUE_NEW(U_TEMP,Z_TEMP,F_TEMP,B_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER) - ELSE - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER) - ENDIF - B1%RHO_F = F_TEMP(1,1,1) - CASE(-3) - IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN - RHO_G_2 = RHOP(BC%IIG,BC%JJG,BC%KKG-1) - RHO_OTHER_2 = RHOP(BC%II,BC%JJ,BC%KK+1) - ENDIF - Z_TEMP(1,1,0:3) = (/RHO_G_2,B1%RHO_G,RHO_OTHER,RHO_OTHER_2/) - U_TEMP(1,1,1) = WW(BC%II,BC%JJ,BC%KK-1) - IF (FLUX_LIMITER_SINGLE_COEF) THEN - B_TEMP(1,1,1) = BFZ(BC%II,BC%JJ,BC%KK-1) - CALL GET_SCALAR_FACE_VALUE_NEW(U_TEMP,Z_TEMP,F_TEMP,B_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER) - ELSE - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER) - ENDIF - B1%RHO_F = F_TEMP(1,1,1) - END SELECT - - ! Species and temperature + ! Density, Species and temperature SINGLE_SPEC_IF: IF (N_TOTAL_SCALARS > 1) THEN - SPECIES_LOOP: DO N=1,N_TOTAL_SCALARS - - RHO_ZZ_G = B1%RHO_G*B1%ZZ_G(N) - RHO_ZZ_G_2 = RHO_ZZ_G ! first-order (default) - RHO_ZZ_OTHER_2 = RHO_ZZ_OTHER(N) - - SELECT CASE(BC%IOR) - CASE( 1) - IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN - RHO_ZZ_G_2 = RHOP(BC%IIG+1,BC%JJG,BC%KKG)*ZZP(BC%IIG+1,BC%JJG,BC%KKG,N) - RHO_ZZ_OTHER_2 = RHOP(BC%II-1,BC%JJ,BC%KK)*ZZP(BC%II-1,BC%JJ,BC%KK,N) - ENDIF - Z_TEMP(0:3,1,1) = (/RHO_ZZ_OTHER_2,RHO_ZZ_OTHER(N),RHO_ZZ_G,RHO_ZZ_G_2/) - U_TEMP(1,1,1) = UU(BC%II,BC%JJ,BC%KK) - IF (FLUX_LIMITER_SINGLE_COEF) THEN - B_TEMP(1,1,1) = BFX(BC%II,BC%JJ,BC%KK) - CALL GET_SCALAR_FACE_VALUE_NEW(U_TEMP,Z_TEMP,F_TEMP,B_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER) - ELSE - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER) - ENDIF - RHO_ZZ_F(N) = F_TEMP(1,1,1) - PBAR_F = PBAR_P(BC%KKG,B1%PRESSURE_ZONE) - CASE(-1) - IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN - RHO_ZZ_G_2 = RHOP(BC%IIG-1,BC%JJG,BC%KKG)*ZZP(BC%IIG-1,BC%JJG,BC%KKG,N) - RHO_ZZ_OTHER_2 = RHOP(BC%II+1,BC%JJ,BC%KK)*ZZP(BC%II+1,BC%JJ,BC%KK,N) - ENDIF - Z_TEMP(0:3,1,1) = (/RHO_ZZ_G_2,RHO_ZZ_G,RHO_ZZ_OTHER(N),RHO_ZZ_OTHER_2/) - U_TEMP(1,1,1) = UU(BC%II-1,BC%JJ,BC%KK) - IF (FLUX_LIMITER_SINGLE_COEF) THEN - B_TEMP(1,1,1) = BFX(BC%II-1,BC%JJ,BC%KK) - CALL GET_SCALAR_FACE_VALUE_NEW(U_TEMP,Z_TEMP,F_TEMP,B_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER) - ELSE - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER) - ENDIF - RHO_ZZ_F(N) = F_TEMP(1,1,1) - PBAR_F = PBAR_P(BC%KKG,B1%PRESSURE_ZONE) - CASE( 2) - IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN - RHO_ZZ_G_2 = RHOP(BC%IIG,BC%JJG+1,BC%KKG)*ZZP(BC%IIG,BC%JJG+1,BC%KKG,N) - RHO_ZZ_OTHER_2 = RHOP(BC%II,BC%JJ-1,BC%KK)*ZZP(BC%II,BC%JJ-1,BC%KK,N) - ENDIF - Z_TEMP(1,0:3,1) = (/RHO_ZZ_OTHER_2,RHO_ZZ_OTHER(N),RHO_ZZ_G,RHO_ZZ_G_2/) - U_TEMP(1,1,1) = VV(BC%II,BC%JJ,BC%KK) - IF (FLUX_LIMITER_SINGLE_COEF) THEN - B_TEMP(1,1,1) = BFY(BC%II,BC%JJ,BC%KK) - CALL GET_SCALAR_FACE_VALUE_NEW(U_TEMP,Z_TEMP,F_TEMP,B_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER) - ELSE - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER) - ENDIF - RHO_ZZ_F(N) = F_TEMP(1,1,1) - PBAR_F = PBAR_P(BC%KKG,B1%PRESSURE_ZONE) - CASE(-2) - IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN - RHO_ZZ_G_2 = RHOP(BC%IIG,BC%JJG-1,BC%KKG)*ZZP(BC%IIG,BC%JJG-1,BC%KKG,N) - RHO_ZZ_OTHER_2 = RHOP(BC%II,BC%JJ+1,BC%KK)*ZZP(BC%II,BC%JJ+1,BC%KK,N) - ENDIF - Z_TEMP(1,0:3,1) = (/RHO_ZZ_G_2,RHO_ZZ_G,RHO_ZZ_OTHER(N),RHO_ZZ_OTHER_2/) - U_TEMP(1,1,1) = VV(BC%II,BC%JJ-1,BC%KK) - IF (FLUX_LIMITER_SINGLE_COEF) THEN - B_TEMP(1,1,1) = BFY(BC%II,BC%JJ-1,BC%KK) - CALL GET_SCALAR_FACE_VALUE_NEW(U_TEMP,Z_TEMP,F_TEMP,B_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER) - ELSE - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER) - ENDIF - RHO_ZZ_F(N) = F_TEMP(1,1,1) - PBAR_F = PBAR_P(BC%KKG,B1%PRESSURE_ZONE) - CASE( 3) - IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN - RHO_ZZ_G_2 = RHOP(BC%IIG,BC%JJG,BC%KKG+1)*ZZP(BC%IIG,BC%JJG,BC%KKG+1,N) - RHO_ZZ_OTHER_2 = RHOP(BC%II,BC%JJ,BC%KK-1)*ZZP(BC%II,BC%JJ,BC%KK-1,N) - ENDIF - Z_TEMP(1,1,0:3) = (/RHO_ZZ_OTHER_2,RHO_ZZ_OTHER(N),RHO_ZZ_G,RHO_ZZ_G_2/) - U_TEMP(1,1,1) = WW(BC%II,BC%JJ,BC%KK) - IF (FLUX_LIMITER_SINGLE_COEF) THEN - B_TEMP(1,1,1) = BFZ(BC%II,BC%JJ,BC%KK) - CALL GET_SCALAR_FACE_VALUE_NEW(U_TEMP,Z_TEMP,F_TEMP,B_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER) - ELSE - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER) - ENDIF - RHO_ZZ_F(N) = F_TEMP(1,1,1) - PBAR_F = (PBAR_P(BC%KK,B1%PRESSURE_ZONE)*DZ(BC%KKG) + PBAR_P(BC%KKG,B1%PRESSURE_ZONE)*DZ(BC%KK)) / & - (DZ(BC%KK)+DZ(BC%KKG)) - CASE(-3) - IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN - RHO_ZZ_G_2 = RHOP(BC%IIG,BC%JJG,BC%KKG-1)*ZZP(BC%IIG,BC%JJG,BC%KKG-1,N) - RHO_ZZ_OTHER_2 = RHOP(BC%II,BC%JJ,BC%KK+1)*ZZP(BC%II,BC%JJ,BC%KK+1,N) - ENDIF - Z_TEMP(1,1,0:3) = (/RHO_ZZ_G_2,RHO_ZZ_G,RHO_ZZ_OTHER(N),RHO_ZZ_OTHER_2/) - U_TEMP(1,1,1) = WW(BC%II,BC%JJ,BC%KK-1) - IF (FLUX_LIMITER_SINGLE_COEF) THEN - B_TEMP(1,1,1) = BFZ(BC%II,BC%JJ,BC%KK-1) - CALL GET_SCALAR_FACE_VALUE_NEW(U_TEMP,Z_TEMP,F_TEMP,B_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER) - ELSE - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER) - ENDIF - RHO_ZZ_F(N) = F_TEMP(1,1,1) - PBAR_F = (PBAR_P(BC%KK,B1%PRESSURE_ZONE)*DZ(BC%KKG) + PBAR_P(BC%KKG,B1%PRESSURE_ZONE)*DZ(BC%KK)) / & - (DZ(BC%KK)+DZ(BC%KKG)) - END SELECT - ENDDO SPECIES_LOOP - - ! face value of temperature IF (ATMOSPHERIC_INTERPOLATION) THEN B1%TMP_F = (TMP(BC%II,BC%JJ,BC%KK)*DZ(BC%KKG) + TMP(BC%IIG,BC%JJG,BC%KKG)*DZ(BC%KK)) / (DZ(BC%KK)+DZ(BC%KKG)) @@ -926,11 +731,29 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I ZZP(BC%IIG,BC%JJG,BC%KKG,1:N_TOTAL_SCALARS)*DZ(BC%KK)) / (DZ(BC%KK)+DZ(BC%KKG)) ZZ_GET(1:N_TRACKED_SPECIES) = B1%ZZ_F(1:N_TRACKED_SPECIES) CALL GET_SPECIFIC_GAS_CONSTANT(ZZ_GET,RSUM_F) + SELECT CASE(ABS(BC%IOR)) + CASE DEFAULT + PBAR_F = PBAR_P(BC%KKG,B1%PRESSURE_ZONE) + CASE(3) + ! for very large cells with stratification the pressure variation in z can be significant + PBAR_F = (PBAR_P(BC%KK,B1%PRESSURE_ZONE)*DZ(BC%KKG) + PBAR_P(BC%KKG,B1%PRESSURE_ZONE)*DZ(BC%KK)) / & + (DZ(BC%KK)+DZ(BC%KKG)) + END SELECT B1%RHO_F = PBAR_F/(RSUM_F*B1%TMP_F) ELSE + SELECT CASE(BC%IOR) + CASE( 1); RHO_ZZ_F(1:N_TOTAL_SCALARS) = FX(BC%II ,BC%JJ,BC%KK,1:N_TOTAL_SCALARS) + CASE(-1); RHO_ZZ_F(1:N_TOTAL_SCALARS) = FX(BC%II-1,BC%JJ,BC%KK,1:N_TOTAL_SCALARS) + CASE( 2); RHO_ZZ_F(1:N_TOTAL_SCALARS) = FY(BC%II,BC%JJ ,BC%KK,1:N_TOTAL_SCALARS) + CASE(-2); RHO_ZZ_F(1:N_TOTAL_SCALARS) = FY(BC%II,BC%JJ-1,BC%KK,1:N_TOTAL_SCALARS) + CASE( 3); RHO_ZZ_F(1:N_TOTAL_SCALARS) = FZ(BC%II,BC%JJ,BC%KK ,1:N_TOTAL_SCALARS) + CASE(-3); RHO_ZZ_F(1:N_TOTAL_SCALARS) = FZ(BC%II,BC%JJ,BC%KK-1,1:N_TOTAL_SCALARS) + END SELECT + B1%RHO_F = SUM(RHO_ZZ_F(1:N_TRACKED_SPECIES)) B1%ZZ_F(1:N_TOTAL_SCALARS) = MAX(0._EB,MIN(1._EB,RHO_ZZ_F(1:N_TOTAL_SCALARS)/B1%RHO_F)) ZZ_GET(1:N_TRACKED_SPECIES) = B1%ZZ_F(1:N_TRACKED_SPECIES) CALL GET_SPECIFIC_GAS_CONSTANT(ZZ_GET,RSUM_F) + PBAR_F = PBAR_P(BC%KKG,B1%PRESSURE_ZONE) B1%TMP_F = PBAR_F/(RSUM_F*B1%RHO_F) ENDIF @@ -1493,7 +1316,7 @@ SUBROUTINE CALCULATE_ZZ_F(T,DT,WALL_INDEX,CFACE_INDEX,PARTICLE_INDEX) ENDIF IF (PREDICTOR) B1%U_NORMAL_S = -UN - IF (CORRECTOR) B1%U_NORMAL = -UN + IF (CORRECTOR) B1%U_NORMAL = -UN END SELECT METHOD_OF_MASS_TRANSFER @@ -1502,6 +1325,7 @@ SUBROUTINE CALCULATE_ZZ_F(T,DT,WALL_INDEX,CFACE_INDEX,PARTICLE_INDEX) IF (PRESENT(WALL_INDEX)) THEN IF (WALL_INDEX<=N_EXTERNAL_WALL_CELLS .AND. .NOT.CELL(IC)%SOLID .AND. .NOT.CELL(ICG)%SOLID) & ZZP(BC%II,BC%JJ,BC%KK,1:N_TRACKED_SPECIES) = 2._EB*B1%ZZ_F(1:N_TRACKED_SPECIES) - B1%ZZ_G(1:N_TRACKED_SPECIES) + ZZP(BC%II2,BC%JJ2,BC%KK2,1:N_TRACKED_SPECIES) = 2._EB*B1%ZZ_F(1:N_TRACKED_SPECIES) - B1%ZZ_G(1:N_TRACKED_SPECIES) ENDIF END SUBROUTINE CALCULATE_ZZ_F @@ -1667,7 +1491,10 @@ SUBROUTINE CALCULATE_RHO_F(BC,B1,WALL_INDEX,CFACE_INDEX) ENDIF IF (PRESENT(WALL_INDEX)) THEN - IF (WALL_INDEX<=N_EXTERNAL_WALL_CELLS .AND. BOUNDARY_TYPE/=OPEN_BOUNDARY) RHOP(BC%II,BC%JJ,BC%KK) = 2._EB*B1%RHO_F - B1%RHO_G + IF (WALL_INDEX<=N_EXTERNAL_WALL_CELLS .AND. BOUNDARY_TYPE/=OPEN_BOUNDARY) THEN + RHOP(BC%II ,BC%JJ, BC%KK ) = 2._EB*B1%RHO_F - B1%RHO_G + RHOP(BC%II2,BC%JJ2,BC%KK2) = 2._EB*B1%RHO_F - B1%RHO_G + ENDIF ENDIF END SUBROUTINE CALCULATE_RHO_F