diff --git a/Source/cons.f90 b/Source/cons.f90 index 3d8df244ef..9b0214be8c 100644 --- a/Source/cons.f90 +++ b/Source/cons.f90 @@ -278,6 +278,7 @@ MODULE GLOBAL_CONSTANTS 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 LOGICAL :: TEST_NEW_CHAR_MODEL=.FALSE. !< Flag to envoke new char model +LOGICAL :: FLUX_LIMITER_MW_CORRECTION=.TRUE. !< Flag for MW correction ensure consistent equation of state at face 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/data.f90 b/Source/data.f90 index 14a33ba92f..a90cfbe240 100644 --- a/Source/data.f90 +++ b/Source/data.f90 @@ -1360,42 +1360,6 @@ SUBROUTINE DEFINE_OUTPUT_QUANTITIES OUTPUT_QUANTITY(553)%UNITS = 'm/s' OUTPUT_QUANTITY(553)%SHORT_NAME = 'V_LS' -OUTPUT_QUANTITY(560)%NAME = 'BFX' -OUTPUT_QUANTITY(560)%UNITS = '' -OUTPUT_QUANTITY(560)%SHORT_NAME = 'bfx' -OUTPUT_QUANTITY(560)%CELL_POSITION = CELL_FACE -OUTPUT_QUANTITY(560)%IOR = 1 - -OUTPUT_QUANTITY(561)%NAME = 'BFY' -OUTPUT_QUANTITY(561)%UNITS = '' -OUTPUT_QUANTITY(561)%SHORT_NAME = 'bfy' -OUTPUT_QUANTITY(561)%CELL_POSITION = CELL_FACE -OUTPUT_QUANTITY(561)%IOR = 2 - -OUTPUT_QUANTITY(562)%NAME = 'BFZ' -OUTPUT_QUANTITY(562)%UNITS = '' -OUTPUT_QUANTITY(562)%SHORT_NAME = 'bfz' -OUTPUT_QUANTITY(562)%CELL_POSITION = CELL_FACE -OUTPUT_QUANTITY(562)%IOR = 3 - -OUTPUT_QUANTITY(563)%NAME = 'BFX MINUS' -OUTPUT_QUANTITY(563)%UNITS = '' -OUTPUT_QUANTITY(563)%SHORT_NAME = 'bfx-' -OUTPUT_QUANTITY(563)%CELL_POSITION = CELL_FACE -OUTPUT_QUANTITY(563)%IOR = 1 - -OUTPUT_QUANTITY(564)%NAME = 'BFY MINUS' -OUTPUT_QUANTITY(564)%UNITS = '' -OUTPUT_QUANTITY(564)%SHORT_NAME = 'bfy-' -OUTPUT_QUANTITY(564)%CELL_POSITION = CELL_FACE -OUTPUT_QUANTITY(564)%IOR = 2 - -OUTPUT_QUANTITY(565)%NAME = 'BFZ MINUS' -OUTPUT_QUANTITY(565)%UNITS = '' -OUTPUT_QUANTITY(565)%SHORT_NAME = 'bfz-' -OUTPUT_QUANTITY(565)%CELL_POSITION = CELL_FACE -OUTPUT_QUANTITY(565)%IOR = 3 - ! Boundary Quantities (Negative indices) OUTPUT_QUANTITY(-1)%NAME = 'RADIATIVE HEAT FLUX' diff --git a/Source/divg.f90 b/Source/divg.f90 index 6c4d2dd524..9e22f3d003 100644 --- a/Source/divg.f90 +++ b/Source/divg.f90 @@ -788,11 +788,11 @@ END SUBROUTINE DIVERGENCE_PART_1 SUBROUTINE ENTHALPY_ADVECTION_NEW(U_DOT_DEL_RHO_H_S) USE PHYSICAL_FUNCTIONS, ONLY: GET_SENSIBLE_ENTHALPY -USE MATH_FUNCTIONS, ONLY: GET_SCALAR_FACE_COEF,GET_SCALAR_FACE_VALUE_NEW +USE MATH_FUNCTIONS, ONLY: GET_SCALAR_FACE_VALUE 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,B_TEMP +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 :: IC,I,J,K,IW TYPE(WALL_TYPE), POINTER :: WC @@ -805,10 +805,6 @@ SUBROUTINE ENTHALPY_ADVECTION_NEW(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,13 +825,9 @@ SUBROUTINE ENTHALPY_ADVECTION_NEW(U_DOT_DEL_RHO_H_S) ! Compute scalar face values -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) +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) ALLOCATE(ZZ_GET(1:N_TRACKED_SPECIES)) @@ -875,8 +867,7 @@ SUBROUTINE ENTHALPY_ADVECTION_NEW(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) - 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) + CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_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 @@ -886,40 +877,35 @@ SUBROUTINE ENTHALPY_ADVECTION_NEW(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) - 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) + CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_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) - 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) + CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_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) - 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) + CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_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) - 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) + CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_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) - 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) + CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_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 @@ -981,12 +967,13 @@ END SUBROUTINE ENTHALPY_ADVECTION_NEW SUBROUTINE SPECIES_ADVECTION_PART_1_NEW -USE MATH_FUNCTIONS, ONLY: GET_SCALAR_FACE_VALUE_NEW,GET_SCALAR_FACE_COEF +USE MATH_FUNCTIONS, ONLY: GET_SCALAR_FACE_VALUE USE PHYSICAL_FUNCTIONS, ONLY: GET_MOLECULAR_WEIGHT ! INTEGER, INTENT(IN) :: N -REAL(EB), POINTER, DIMENSION(:,:,:) :: RHO_Z_P +REAL(EB), POINTER, DIMENSION(:,:,:) :: RHO_Z_P,RHO_RMW REAL(EB), POINTER, DIMENSION(:,:,:,:) :: FX_ZZ,FY_ZZ,FZ_ZZ -REAL(EB), DIMENSION(0:3,0:3,0:3) :: U_TEMP,F_TEMP,B_TEMP +REAL(EB) :: ZZ_GET(1:N_TRACKED_SPECIES),MW_G +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 @@ -998,31 +985,6 @@ SUBROUTINE SPECIES_ADVECTION_PART_1_NEW FZ_ZZ=>SWORK3 RHO_Z_P=>WORK_PAD -BFX = 1.E10_EB -BFY = 1.E10_EB -BFZ = 1.E10_EB - -FACE_COEF_LOOP: DO N=1,N_TOTAL_SCALARS - - !$OMP PARALLEL DO PRIVATE(I,J,K) - 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) - END DO - END DO - END DO - !$OMP END PARALLEL DO - - ! Compute scalar face values - - CALL GET_SCALAR_FACE_COEF(UU,RHO_Z_P,BFX,0,IBAR,1,JBAR,1,KBAR,1,I_FLUX_LIMITER) - CALL GET_SCALAR_FACE_COEF(VV,RHO_Z_P,BFY,1,IBAR,0,JBAR,1,KBAR,2,I_FLUX_LIMITER) - CALL GET_SCALAR_FACE_COEF(WW,RHO_Z_P,BFZ,1,IBAR,1,JBAR,0,KBAR,3,I_FLUX_LIMITER) - -ENDDO FACE_COEF_LOOP - - ! Species face values SPECIES_LOOP: DO N=1,N_TOTAL_SCALARS @@ -1041,11 +1003,11 @@ SUBROUTINE SPECIES_ADVECTION_PART_1_NEW ! Compute scalar face values - CALL GET_SCALAR_FACE_VALUE_NEW(UU,RHO_Z_P,FX_ZZ(:,:,:,N),BFX,0,IBAR,1,JBAR,1,KBAR,1,I_FLUX_LIMITER) - CALL GET_SCALAR_FACE_VALUE_NEW(VV,RHO_Z_P,FY_ZZ(:,:,:,N),BFY,1,IBAR,0,JBAR,1,KBAR,2,I_FLUX_LIMITER) - CALL GET_SCALAR_FACE_VALUE_NEW(WW,RHO_Z_P,FZ_ZZ(:,:,:,N),BFZ,1,IBAR,1,JBAR,0,KBAR,3,I_FLUX_LIMITER) + CALL GET_SCALAR_FACE_VALUE(UU,RHO_Z_P,FX_ZZ(:,:,:,N),1,IBM1,1,JBAR,1,KBAR,1,I_FLUX_LIMITER) + CALL GET_SCALAR_FACE_VALUE(VV,RHO_Z_P,FY_ZZ(:,:,:,N),1,IBAR,1,JBM1,1,KBAR,2,I_FLUX_LIMITER) + CALL GET_SCALAR_FACE_VALUE(WW,RHO_Z_P,FZ_ZZ(:,:,:,N),1,IBAR,1,JBAR,1,KBM1,3,I_FLUX_LIMITER) - !$OMP PARALLEL DO PRIVATE(IW,WC,BC,B1,U_TEMP,Z_TEMP,F_TEMP,B_TEMP) + !$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 @@ -1064,8 +1026,7 @@ SUBROUTINE SPECIES_ADVECTION_PART_1_NEW 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) - 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) + 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 @@ -1075,40 +1036,35 @@ SUBROUTINE SPECIES_ADVECTION_PART_1_NEW 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) - 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) + 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) - 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) + 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) - 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) + 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) - 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) + 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) - 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) + 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 @@ -1120,6 +1076,126 @@ SUBROUTINE SPECIES_ADVECTION_PART_1_NEW ENDDO SPECIES_LOOP +FACE_CORRECTION_IF: IF (FLUX_LIMITER_MW_CORRECTION) THEN + + ! Repeat the above for DENSITY + + RHO_RMW=>WORK_PAD + + !$OMP PARALLEL DO PRIVATE(ZZ_GET,MW_G) SCHEDULE(STATIC) + DO K=-1,KBP1+1 + DO J=-1,JBP1+1 + DO I=-1,IBP1+1 + ZZ_GET(1:N_TRACKED_SPECIES) = ZZP(I,J,K,1:N_TRACKED_SPECIES) + CALL GET_MOLECULAR_WEIGHT(ZZ_GET,MW_G) + RHO_RMW(I,J,K) = RHOP(I,J,K)/MW_G + ENDDO + ENDDO + ENDDO + !$OMP END PARALLEL DO + + CALL GET_SCALAR_FACE_VALUE(UU,RHO_RMW,FX_ZZ(:,:,:,0),1,IBM1,1,JBAR,1,KBAR,1,I_FLUX_LIMITER) + CALL GET_SCALAR_FACE_VALUE(VV,RHO_RMW,FY_ZZ(:,:,:,0),1,IBAR,1,JBM1,1,KBAR,2,I_FLUX_LIMITER) + CALL GET_SCALAR_FACE_VALUE(WW,RHO_RMW,FZ_ZZ(:,:,:,0),1,IBAR,1,JBAR,1,KBM1,3,I_FLUX_LIMITER) + + !$OMP PARALLEL DO PRIVATE(IW,WC,BC,B1,U_TEMP,Z_TEMP,F_TEMP,ZZ_GET) + 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) + + ! 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(BC%IOR) + CASE( 1) OFF_WALL_SELECT_3 + ! 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,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(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,0) = F_TEMP(1,1,1) + ENDIF + CASE( 2) OFF_WALL_SELECT_3 + 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,0) = F_TEMP(1,1,1) + ENDIF + CASE(-2) OFF_WALL_SELECT_3 + 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,0) = F_TEMP(1,1,1) + ENDIF + CASE( 3) OFF_WALL_SELECT_3 + 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,0) = F_TEMP(1,1,1) + ENDIF + CASE(-3) OFF_WALL_SELECT_3 + 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,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 max face value of (RHO*ZZ) such that SUM(RHO*ZZ/MW)_FACE = RHO_FACE/MW_FACE + ! (necessary condition to preserve isothermal flow) + + !$OMP PARALLEL DO PRIVATE(N,MW_G) SCHEDULE(STATIC) + DO K=0,KBAR + DO J=0,JBAR + DO I=0,IBAR + N=MAXLOC(FX_ZZ(I,J,K,1:N_TRACKED_SPECIES),1) + MW_G = SPECIES_MIXTURE(N)%MW + FX_ZZ(I,J,K,N) = MW_G*MAX( 0._EB, FX_ZZ(I,J,K,0) & + - SUM(FX_ZZ(I,J,K,1:(N-1))/SPECIES_MIXTURE(1:(N-1))%MW) & + - SUM(FX_ZZ(I,J,K,(N+1):N_TRACKED_SPECIES)/SPECIES_MIXTURE((N+1):N_TRACKED_SPECIES)%MW) ) + + N=MAXLOC(FY_ZZ(I,J,K,1:N_TRACKED_SPECIES),1) + MW_G = SPECIES_MIXTURE(N)%MW + FY_ZZ(I,J,K,N) = MW_G*MAX( 0._EB, FY_ZZ(I,J,K,0) & + - SUM(FY_ZZ(I,J,K,1:(N-1))/SPECIES_MIXTURE(1:(N-1))%MW) & + - SUM(FY_ZZ(I,J,K,(N+1):N_TRACKED_SPECIES)/SPECIES_MIXTURE((N+1):N_TRACKED_SPECIES)%MW) ) + + N=MAXLOC(FZ_ZZ(I,J,K,1:N_TRACKED_SPECIES),1) + MW_G = SPECIES_MIXTURE(N)%MW + FZ_ZZ(I,J,K,N) = MW_G*MAX( 0._EB, FZ_ZZ(I,J,K,0) & + - SUM(FZ_ZZ(I,J,K,1:(N-1))/SPECIES_MIXTURE(1:(N-1))%MW) & + - SUM(FZ_ZZ(I,J,K,(N+1):N_TRACKED_SPECIES)/SPECIES_MIXTURE((N+1):N_TRACKED_SPECIES)%MW) ) + ENDDO + ENDDO + ENDDO + !$OMP END PARALLEL DO + +ENDIF FACE_CORRECTION_IF + END SUBROUTINE SPECIES_ADVECTION_PART_1_NEW diff --git a/Source/dump.f90 b/Source/dump.f90 index 7de3abe6f2..6775b9a0bc 100644 --- a/Source/dump.f90 +++ b/Source/dump.f90 @@ -8583,20 +8583,7 @@ REAL(EB) RECURSIVE FUNCTION GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,IND,IND2,Y_INDEX,Z CASE(553) ! V_LS GAS_PHASE_OUTPUT_RES = V_LS(II,JJ) - CASE(560) ! BFX - GAS_PHASE_OUTPUT_RES = BFX(II,JJ,KK) - CASE(561) ! BFY - GAS_PHASE_OUTPUT_RES = BFY(II,JJ,KK) - CASE(562) ! BFZ - GAS_PHASE_OUTPUT_RES = BFZ(II,JJ,KK) - CASE(563) ! BFX MINUS - GAS_PHASE_OUTPUT_RES = BFX(MAX(0,II-1),JJ,KK) - CASE(564) ! BFY MINUS - GAS_PHASE_OUTPUT_RES = BFY(II,MAX(0,JJ-1),KK) - CASE(565) ! BFZ MINUS - GAS_PHASE_OUTPUT_RES = BFZ(II,JJ,MAX(0,KK-1)) - - END SELECT IND_SELECT +END SELECT IND_SELECT ! Fill GAS_PHASE_OUTPUT for CUT_CELLs. ! Some variables have already been filled in fire.f90 diff --git a/Source/init.f90 b/Source/init.f90 index aa4ddfb6ea..ae3e6173ff 100644 --- a/Source/init.f90 +++ b/Source/init.f90 @@ -499,7 +499,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,IZERO,NS,IW +INTEGER :: N,I,J,K,IZERO,N_LOWER_SCALARS,NS,IW REAL(EB), INTENT(IN) :: DT INTEGER, INTENT(IN) :: NM REAL(EB) :: MU_N,CS,DELTA @@ -595,6 +595,20 @@ 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 + +! Required for cell face density correction for multicomponent mixtures + +IF (FLUX_LIMITER_MW_CORRECTION) 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) ; M%FX=0._EB +ALLOCATE( M%FY(0:IBP1,0:JBP1,0:KBP1,N_LOWER_SCALARS:N_TOTAL_SCALARS),STAT=IZERO); CALL ChkMemErr('INIT','FY',IZERO) ; M%FY=0._EB +ALLOCATE( M%FZ(0:IBP1,0:JBP1,0:KBP1,N_LOWER_SCALARS:N_TOTAL_SCALARS),STAT=IZERO); CALL ChkMemErr('INIT','FZ',IZERO) ; M%FZ=0._EB + ! Allocate storage for scalar total fluxes IF (STORE_SPECIES_FLUX) THEN @@ -682,9 +696,9 @@ SUBROUTINE INITIALIZE_MESH_VARIABLES_1(DT,NM) ALLOCATE(M%WORK9(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO) ; CALL ChkMemErr('INIT','WORK9',IZERO) ALLOCATE(M%IWORK1(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO) ; CALL ChkMemErr('INIT','IWORK1',IZERO) -ALLOCATE(M%SWORK1(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO) ; CALL ChkMemErr('INIT','SWORK1',IZERO) -ALLOCATE(M%SWORK2(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO) ; CALL ChkMemErr('INIT','SWORK2',IZERO) -ALLOCATE(M%SWORK3(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO) ; CALL ChkMemErr('INIT','SWORK3',IZERO) +ALLOCATE(M%SWORK1(0:IBP1,0:JBP1,0:KBP1,N_LOWER_SCALARS:N_TOTAL_SCALARS),STAT=IZERO) ; CALL ChkMemErr('INIT','SWORK1',IZERO) +ALLOCATE(M%SWORK2(0:IBP1,0:JBP1,0:KBP1,N_LOWER_SCALARS:N_TOTAL_SCALARS),STAT=IZERO) ; CALL ChkMemErr('INIT','SWORK2',IZERO) +ALLOCATE(M%SWORK3(0:IBP1,0:JBP1,0:KBP1,N_LOWER_SCALARS:N_TOTAL_SCALARS),STAT=IZERO) ; CALL ChkMemErr('INIT','SWORK3',IZERO) ALLOCATE(M%SWORK4(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO) ; CALL ChkMemErr('INIT','SWORK4',IZERO) M%IWORK1=0 M%SWORK1=0._EB @@ -819,23 +833,6 @@ 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/mass.f90 b/Source/mass.f90 index 3c7975c6e7..b234805383 100644 --- a/Source/mass.f90 +++ b/Source/mass.f90 @@ -19,18 +19,18 @@ MODULE MASS SUBROUTINE MASS_FINITE_DIFFERENCES_NEW(NM) -USE MATH_FUNCTIONS, ONLY: GET_SCALAR_FACE_VALUE_NEW,GET_SCALAR_FACE_COEF +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,B_TEMP +REAL(EB) :: TNOW,MW_F,MW_G,ZZ_GET(1:N_TRACKED_SPECIES) +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 +REAL(EB), POINTER, DIMENSION(:,:,:) :: RHO_Z_P,RHO_RMW TYPE(WALL_TYPE), POINTER :: WC TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1 @@ -59,40 +59,11 @@ SUBROUTINE MASS_FINITE_DIFFERENCES_NEW(NM) IF (PREDICTOR) DT_RESTRICT_COUNT = 0 -! First, compute the minimum face coefficient for each direction - -BFX = 1.E10_EB -BFY = 1.E10_EB -BFZ = 1.E10_EB - -FACE_COEF_LOOP: DO N=1,N_TOTAL_SCALARS - - RHO_Z_P=>WORK_PAD - - !$OMP PARALLEL DO PRIVATE(I,J,K) - 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) - END DO - END DO - END DO - !$OMP END PARALLEL DO - - ! Compute scalar face values - - CALL GET_SCALAR_FACE_COEF(UU,RHO_Z_P,BFX,0,IBAR,1,JBAR,1,KBAR,1,I_FLUX_LIMITER) - CALL GET_SCALAR_FACE_COEF(VV,RHO_Z_P,BFY,1,IBAR,0,JBAR,1,KBAR,2,I_FLUX_LIMITER) - CALL GET_SCALAR_FACE_COEF(WW,RHO_Z_P,BFZ,1,IBAR,1,JBAR,0,KBAR,3,I_FLUX_LIMITER) - -ENDDO FACE_COEF_LOOP - ! Species face values -SPECIES_LOOP: DO N=1,N_TOTAL_SCALARS +RHO_Z_P=>WORK_PAD - ! this is redundant and inefficient, but fine for proof of concept - RHO_Z_P=>WORK_PAD +SPECIES_LOOP: DO N=1,N_TOTAL_SCALARS !$OMP PARALLEL DO PRIVATE(I,J,K) DO K = -1, KBP1+1 @@ -107,11 +78,11 @@ SUBROUTINE MASS_FINITE_DIFFERENCES_NEW(NM) ! Compute scalar face values - CALL GET_SCALAR_FACE_VALUE_NEW(UU,RHO_Z_P,FX(:,:,:,N),BFX,0,IBAR,1,JBAR,1,KBAR,1,I_FLUX_LIMITER) - CALL GET_SCALAR_FACE_VALUE_NEW(VV,RHO_Z_P,FY(:,:,:,N),BFY,1,IBAR,0,JBAR,1,KBAR,2,I_FLUX_LIMITER) - CALL GET_SCALAR_FACE_VALUE_NEW(WW,RHO_Z_P,FZ(:,:,:,N),BFZ,1,IBAR,1,JBAR,0,KBAR,3,I_FLUX_LIMITER) + CALL GET_SCALAR_FACE_VALUE(UU,RHO_Z_P,FX(:,:,:,N),1,IBM1,1,JBAR,1,KBAR,1,I_FLUX_LIMITER) + CALL GET_SCALAR_FACE_VALUE(VV,RHO_Z_P,FY(:,:,:,N),1,IBAR,1,JBM1,1,KBAR,2,I_FLUX_LIMITER) + CALL GET_SCALAR_FACE_VALUE(WW,RHO_Z_P,FZ(:,:,:,N),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,B_TEMP) + !$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 @@ -162,8 +133,7 @@ SUBROUTINE MASS_FINITE_DIFFERENCES_NEW(NM) 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) - B_TEMP(1,1,1) = BFX(II+1,JJ,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) + 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 @@ -173,40 +143,35 @@ SUBROUTINE MASS_FINITE_DIFFERENCES_NEW(NM) 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) - B_TEMP(1,1,1) = BFX(II-2,JJ,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) + 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) - B_TEMP(1,1,1) = BFY(II,JJ+1,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) + 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) - B_TEMP(1,1,1) = BFY(II,JJ-2,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) + 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) - B_TEMP(1,1,1) = BFZ(II,JJ,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) + 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) - B_TEMP(1,1,1) = BFZ(II,JJ,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) + 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 @@ -218,6 +183,157 @@ SUBROUTINE MASS_FINITE_DIFFERENCES_NEW(NM) ENDDO SPECIES_LOOP +FACE_CORRECTION_IF: IF (FLUX_LIMITER_MW_CORRECTION) THEN + + ! Repeat the above for DENSITY + + RHO_RMW=>WORK_PAD + + !$OMP PARALLEL DO PRIVATE(ZZ_GET,MW_G) SCHEDULE(STATIC) + DO K=-1,KBP1+1 + DO J=-1,JBP1+1 + DO I=-1,IBP1+1 + ZZ_GET(1:N_TRACKED_SPECIES) = ZZP(I,J,K,1:N_TRACKED_SPECIES) + CALL GET_MOLECULAR_WEIGHT(ZZ_GET,MW_G) + RHO_RMW(I,J,K) = RHOP(I,J,K)/MW_G + ENDDO + ENDDO + ENDDO + !$OMP END PARALLEL DO + + CALL GET_SCALAR_FACE_VALUE(UU,RHO_RMW,FX(:,:,:,0),1,IBM1,1,JBAR,1,KBAR,1,I_FLUX_LIMITER) + CALL GET_SCALAR_FACE_VALUE(VV,RHO_RMW,FY(:,:,:,0),1,IBAR,1,JBM1,1,KBAR,2,I_FLUX_LIMITER) + CALL GET_SCALAR_FACE_VALUE(WW,RHO_RMW,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,ZZ_GET,MW_F) + 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 + ZZ_GET(1:N_TRACKED_SPECIES) = B1%ZZ_F(1:N_TRACKED_SPECIES) + CALL GET_MOLECULAR_WEIGHT(ZZ_GET,MW_F) + SELECT CASE(IOR) + CASE( 1); FX(IIG-1,JJG,KKG,0) = B1%RHO_F/MW_F + CASE(-1); FX(IIG,JJG,KKG,0) = B1%RHO_F/MW_F + CASE( 2); FY(IIG,JJG-1,KKG,0) = B1%RHO_F/MW_F + CASE(-2); FY(IIG,JJG,KKG,0) = B1%RHO_F/MW_F + CASE( 3); FZ(IIG,JJG,KKG-1,0) = B1%RHO_F/MW_F + CASE(-3); FZ(IIG,JJG,KKG,0) = B1%RHO_F/MW_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) = (/RHO_RMW(II+1,JJ,KK),RHO_RMW(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,RHO_RMW(II-2:II-1,JJ,KK),RHO_RMW(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) = (/RHO_RMW(II,JJ+1,KK),RHO_RMW(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,RHO_RMW(II,JJ-2:JJ-1,KK),RHO_RMW(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) = (/RHO_RMW(II,JJ,KK+1),RHO_RMW(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,RHO_RMW(II,JJ,KK-2:KK-1),RHO_RMW(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 max face value of (RHO*ZZ) such that SUM(RHO*ZZ/MW)_FACE = RHO_FACE/MW_FACE + ! (necessary condition to preserve isothermal flow) + + !$OMP PARALLEL DO PRIVATE(N,MW_G) 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) + MW_G = SPECIES_MIXTURE(N)%MW + FX(I,J,K,N) = MW_G*MAX( 0._EB, FX(I,J,K,0) & + - SUM(FX(I,J,K,1:(N-1))/SPECIES_MIXTURE(1:(N-1))%MW) & + - SUM(FX(I,J,K,(N+1):N_TRACKED_SPECIES)/SPECIES_MIXTURE((N+1):N_TRACKED_SPECIES)%MW) ) + + N=MAXLOC(FY(I,J,K,1:N_TRACKED_SPECIES),1) + MW_G = SPECIES_MIXTURE(N)%MW + FY(I,J,K,N) = MW_G*MAX( 0._EB, FY(I,J,K,0) & + - SUM(FY(I,J,K,1:(N-1))/SPECIES_MIXTURE(1:(N-1))%MW) & + - SUM(FY(I,J,K,(N+1):N_TRACKED_SPECIES)/SPECIES_MIXTURE((N+1):N_TRACKED_SPECIES)%MW) ) + + N=MAXLOC(FZ(I,J,K,1:N_TRACKED_SPECIES),1) + MW_G = SPECIES_MIXTURE(N)%MW + FZ(I,J,K,N) = MW_G*MAX( 0._EB, FZ(I,J,K,0) & + - SUM(FZ(I,J,K,1:(N-1))/SPECIES_MIXTURE(1:(N-1))%MW) & + - SUM(FZ(I,J,K,(N+1):N_TRACKED_SPECIES)/SPECIES_MIXTURE((N+1):N_TRACKED_SPECIES)%MW) ) + 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_NEW diff --git a/Source/mesh.f90 b/Source/mesh.f90 index 938cd3ee08..66f0ff6d1d 100644 --- a/Source/mesh.f90 +++ b/Source/mesh.f90 @@ -60,9 +60,6 @@ MODULE MESH_VARIABLES REAL(EB), ALLOCATABLE, DIMENSION(:,:,:) :: D_Z_MAX !< \f$ \max D_\alpha \f$ REAL(EB), ALLOCATABLE, DIMENSION(:,:,:) :: PP_RESIDUAL !< Pressure Poisson residual (debug) REAL(EB), ALLOCATABLE, DIMENSION(:,:,:) :: LES_FILTER_WIDTH !< Characteristic cell dimension (m) - REAL(EB), ALLOCATABLE, DIMENSION(:,:,:) :: BFX !< Flux limiter face B value for X-direction faces - REAL(EB), ALLOCATABLE, DIMENSION(:,:,:) :: BFY !< Flux limiter face B value for Y-direction faces - REAL(EB), ALLOCATABLE, DIMENSION(:,:,:) :: BFZ !< Flux limiter face B value for Z-direction faces REAL(EB), ALLOCATABLE, DIMENSION(:,:,:,:) :: ZZ !< Lumped species, current time step, \f$ Z_{\alpha,ijk}^n \f$ REAL(EB), ALLOCATABLE, DIMENSION(:,:,:,:) :: ZZS !< Lumped species, next time step, \f$ Z_{\alpha,ijk}^* \f$ @@ -355,7 +352,7 @@ MODULE MESH_POINTERS REAL(EB), POINTER, DIMENSION(:,:,:) :: & U,V,W,US,VS,WS,DDDT,D,DS,H,HS,KRES,FVX,FVY,FVZ,FVX_B,FVY_B,FVZ_B,FVX_D,FVY_D,FVZ_D,RHO,RHOS, & MU,MU_DNS,TMP,Q,KAPPA_GAS,CHI_R,QR,QR_W,RADIATION_EMISSION,RADIATION_ABSORPTION,UII,RSUM,D_SOURCE, & - CSD2,MTR,MSR,WEM,MIX_TIME,CHEM_SUBIT,STRAIN_RATE,D_Z_MAX,PP_RESIDUAL,LES_FILTER_WIDTH,BFX,BFY,BFZ + CSD2,MTR,MSR,WEM,MIX_TIME,CHEM_SUBIT,STRAIN_RATE,D_Z_MAX,PP_RESIDUAL,LES_FILTER_WIDTH REAL(EB), POINTER, DIMENSION(:,:,:,:) :: ZZ,ZZS,REAC_SOURCE_TERM,DEL_RHO_D_DEL_Z,FX,FY,FZ, & SWORK1,SWORK2,SWORK3,SWORK4, & Q_REAC,AVG_DROP_DEN,AVG_DROP_TMP,AVG_DROP_RAD,AVG_DROP_AREA,M_DOT_PPP, & @@ -552,9 +549,6 @@ SUBROUTINE POINT_TO_MESH(NM) FX=>M%FX FY=>M%FY FZ=>M%FZ -BFX=>M%BFX -BFY=>M%BFY -BFZ=>M%BFZ POIS_PTB=>M%POIS_PTB POIS_ERR=>M%POIS_ERR SAVE1=>M%SAVE1 diff --git a/Source/read.f90 b/Source/read.f90 index b8142c5dd2..5e5c2acf7a 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,& + FIXED_LES_FILTER_WIDTH,FLUX_LIMITER,FLUX_LIMITER_MW_CORRECTION,& 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,& @@ -3505,7 +3505,7 @@ SUBROUTINE READ_SPEC ENDDO ENDIF - +! IF (N_TRACKED_SPECIES < 3) FLUX_LIMITER_MW_CORRECTION = .FALSE. CONTAINS diff --git a/Source/wall.f90 b/Source/wall.f90 index 6d3e76b50d..b482fffa5b 100644 --- a/Source/wall.f90 +++ b/Source/wall.f90 @@ -31,7 +31,7 @@ SUBROUTINE WALL_BC(T,DT,NM) REAL(EB) :: TNOW REAL(EB), INTENT(IN) :: T,DT INTEGER, INTENT(IN) :: NM -LOGICAL :: CALL_HT_1D +LOGICAL :: CALL_HT_1D,SECOND_ORDER_INTERPOLATED_BOUNDARY REAL(EB) :: DT_BC,SLIP_COEF INTEGER :: IW,IP,ICF,ITW TYPE(WALL_TYPE), POINTER :: WC @@ -94,7 +94,7 @@ SUBROUTINE WALL_BC(T,DT,NM) WC=>WALL(IW) BC => BOUNDARY_COORD(WC%BC_INDEX) B1 => BOUNDARY_PROP1(WC%B1_INDEX) - IF (IW<=N_EXTERNAL_WALL_CELLS) CALL ASSIGN_GHOST_VALUE(IW,BC,B1) + IF (IW<=N_EXTERNAL_WALL_CELLS) CALL ASSIGN_GHOST_VALUE(IW,BC,B1,SECOND_ORDER_INTERPOLATED_BOUNDARY) IF (WC%BOUNDARY_TYPE==NULL_BOUNDARY) CYCLE WALL_CELL_LOOP_0 SF => SURFACE(WC%SURF_INDEX) CALL NEAR_SURFACE_GAS_VARIABLES(T,SF,BC,B1,WALL_INDEX=IW) @@ -135,7 +135,7 @@ SUBROUTINE WALL_BC(T,DT,NM) SF => SURFACE(WC%SURF_INDEX) IF (.NOT.SF%THERMAL_BC_INDEX==THERMALLY_THICK) THEN - CALL SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX=IW) + CALL SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX=IW,SOIB=SECOND_ORDER_INTERPOLATED_BOUNDARY) ELSEIF (CALL_HT_1D) THEN CALL SOLID_HEAT_TRANSFER(NM,T,SF%HT_DIM*DT_BC,WALL_INDEX=IW) ENDIF @@ -262,15 +262,16 @@ SUBROUTINE WALL_BC(T,DT,NM) END SUBROUTINE WALL_BC -SUBROUTINE ASSIGN_GHOST_VALUE(WALL_INDEX,BC,B1) +SUBROUTINE ASSIGN_GHOST_VALUE(WALL_INDEX,BC,B1,SECOND_ORDER_INTERPOLATED_BOUNDARY) USE PHYSICAL_FUNCTIONS, ONLY : GET_SPECIFIC_GAS_CONSTANT USE COMPLEX_GEOMETRY, ONLY : CC_CGSC, CC_SOLID INTEGER, INTENT(IN) :: WALL_INDEX +LOGICAL, INTENT(OUT) :: SECOND_ORDER_INTERPOLATED_BOUNDARY 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),PBAR_P_2 INTEGER :: IIO,JJO,KKO,II2,JJ2,KK2,ICG,ICO -LOGICAL :: CC_SOLID_FLAG,SECOND_ORDER_INTERPOLATED_BOUNDARY +LOGICAL :: CC_SOLID_FLAG REAL(EB), POINTER, DIMENSION(:,:,:,:) :: OM_ZZP REAL(EB), POINTER, DIMENSION(:,:,:) :: OM_RHOP TYPE(EXTERNAL_WALL_TYPE), POINTER :: EWC @@ -473,25 +474,29 @@ END SUBROUTINE NEAR_SURFACE_GAS_VARIABLES !> \param PARTICLE_INDEX Optional Lagrangian particle index !> \details Calculate the surface temperature TMP_F of either a WALL cell, immersed CFACE cell, or a Lagrangian particle. -SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_INDEX) +SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_INDEX,SOIB) -USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP,INTERPOLATE1D_UNIFORM,GET_SCALAR_FACE_VALUE,GET_SCALAR_FACE_VALUE_NEW,GET_SCALAR_FACE_COEF +USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP,INTERPOLATE1D_UNIFORM,GET_SCALAR_FACE_VALUE 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 +LOGICAL, INTENT(IN), OPTIONAL :: SOIB !Second Order Interpolated Boundary 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_ZZ_F(1:N_TOTAL_SCALARS),ZZ_GET(1:N_TRACKED_SPECIES),ZZ_GET_OTHER(1:N_TRACKED_SPECIES),& 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 + BBB,CCC,PPP,QQQ,RRR,UUU,YYY,WWW,HTC_OLD,RSC_LOC,DTMP,RSUM_G,MU_G,MW_G,MW_G_2,MW_OTHER,MW_OTHER_2, & + RHO_G_2,RHO_OTHER_2,RHO_ZZ_G,RHO_ZZ_G_2,RHO_ZZ_OTHER_2,RHO_ZZ_OTHER(1:N_TOTAL_SCALARS) LOGICAL :: ATMOSPHERIC_INTERPOLATION INTEGER :: IIO,JJO,KKO,N,ADCOUNT REAL(EB), POINTER, DIMENSION(:,:,:) :: OM_RHOP 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 +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 @@ -726,7 +731,199 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I ! Density, Species and temperature + ! Correction for variable molecular weights + + IF (FLUX_LIMITER_MW_CORRECTION) THEN + + 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 + + IIO = EWC%IIO_MIN + JJO = EWC%JJO_MIN + + ZZ_GET(1:N_TRACKED_SPECIES) = ZZP(BC%IIG,BC%JJG,BC%KKG,1:N_TRACKED_SPECIES); CALL GET_MOLECULAR_WEIGHT(ZZ_GET,MW_G) + ZZ_GET(1:N_TRACKED_SPECIES) = ZZP(BC%II,BC%JJ,BC%KK,1:N_TRACKED_SPECIES) ; CALL GET_MOLECULAR_WEIGHT(ZZ_GET,MW_OTHER) + IF (SOIB) THEN + SELECT CASE(BC%IOR) + CASE( 1) + ZZ_GET(1:N_TRACKED_SPECIES) = ZZP(BC%IIG+1,BC%JJG,BC%KKG,1:N_TRACKED_SPECIES) + ZZ_GET_OTHER(1:N_TRACKED_SPECIES) = OM_ZZP(IIO-1,JJO,KKO,1:N_TRACKED_SPECIES) + CASE(-1) + ZZ_GET(1:N_TRACKED_SPECIES) = ZZP(BC%IIG-1,BC%JJG,BC%KKG,1:N_TRACKED_SPECIES) + ZZ_GET_OTHER(1:N_TRACKED_SPECIES) = OM_ZZP(IIO+1,JJO,KKO,1:N_TRACKED_SPECIES) + CASE( 2) + ZZ_GET(1:N_TRACKED_SPECIES) = ZZP(BC%IIG,BC%JJG+1,BC%KKG,1:N_TRACKED_SPECIES) + ZZ_GET_OTHER(1:N_TRACKED_SPECIES) = OM_ZZP(IIO,JJO-1,KKO,1:N_TRACKED_SPECIES) + CASE(-2) + ZZ_GET(1:N_TRACKED_SPECIES) = ZZP(BC%IIG,BC%JJG-1,BC%KKG,1:N_TRACKED_SPECIES) + ZZ_GET_OTHER(1:N_TRACKED_SPECIES) = OM_ZZP(IIO,JJO+1,KKO,1:N_TRACKED_SPECIES) + CASE( 3) + ZZ_GET(1:N_TRACKED_SPECIES) = ZZP(BC%IIG,BC%JJG,BC%KKG+1,1:N_TRACKED_SPECIES) + ZZ_GET_OTHER(1:N_TRACKED_SPECIES) = OM_ZZP(IIO,JJO,KKO-1,1:N_TRACKED_SPECIES) + CASE(-3) + ZZ_GET(1:N_TRACKED_SPECIES) = ZZP(BC%IIG,BC%JJG,BC%KKG-1,1:N_TRACKED_SPECIES) + ZZ_GET_OTHER(1:N_TRACKED_SPECIES) = OM_ZZP(IIO,JJO,KKO+1,1:N_TRACKED_SPECIES) + END SELECT + CALL GET_MOLECULAR_WEIGHT(ZZ_GET,MW_G_2) + CALL GET_MOLECULAR_WEIGHT(ZZ_GET_OTHER,MW_OTHER_2) + ELSE + MW_G_2 = MW_G + MW_OTHER_2 = MW_OTHER + ENDIF + ENDIF + + SELECT CASE(BC%IOR) + CASE( 1) + IF (SOIB) THEN + RHO_G_2 = RHOP(BC%IIG+1,BC%JJG,BC%KKG) + RHO_OTHER_2 = OM_RHOP(IIO-1,JJO,KKO) + ENDIF + Z_TEMP(0:3,1,1) = (/RHO_OTHER_2,RHO_OTHER,B1%RHO_G,RHO_G_2/) + IF (FLUX_LIMITER_MW_CORRECTION) Z_TEMP(0:3,1,1) = Z_TEMP(0:3,1,1)/(/MW_OTHER_2,MW_OTHER,MW_G,MW_G_2/) + U_TEMP(1,1,1) = UU(BC%II,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) + B1%RHO_F = F_TEMP(1,1,1) + CASE(-1) + IF (SOIB) THEN + RHO_G_2 = RHOP(BC%IIG-1,BC%JJG,BC%KKG) + RHO_OTHER_2 = OM_RHOP(IIO+1,JJO,KKO) + ENDIF + Z_TEMP(0:3,1,1) = (/RHO_G_2,B1%RHO_G,RHO_OTHER,RHO_OTHER_2/) + IF (FLUX_LIMITER_MW_CORRECTION) Z_TEMP(0:3,1,1) = Z_TEMP(0:3,1,1)/(/MW_G_2,MW_G,MW_OTHER,MW_OTHER_2/) + 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) + B1%RHO_F = F_TEMP(1,1,1) + CASE( 2) + IF (SOIB) THEN + RHO_G_2 = RHOP(BC%IIG,BC%JJG+1,BC%KKG) + RHO_OTHER_2 = OM_RHOP(IIO,JJO-1,KKO) + ENDIF + Z_TEMP(1,0:3,1) = (/RHO_OTHER_2,RHO_OTHER,B1%RHO_G,RHO_G_2/) + IF (FLUX_LIMITER_MW_CORRECTION) Z_TEMP(1,0:3,1) = Z_TEMP(1,0:3,1)/(/MW_OTHER_2,MW_OTHER,MW_G,MW_G_2/) + U_TEMP(1,1,1) = VV(BC%II,BC%JJ,BC%KK) + CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER) + B1%RHO_F = F_TEMP(1,1,1) + CASE(-2) + IF (SOIB) THEN + RHO_G_2 = RHOP(BC%IIG,BC%JJG-1,BC%KKG) + RHO_OTHER_2 = OM_RHOP(IIO,JJO+1,KKO) + ENDIF + Z_TEMP(1,0:3,1) = (/RHO_G_2,B1%RHO_G,RHO_OTHER,RHO_OTHER_2/) + IF (FLUX_LIMITER_MW_CORRECTION) Z_TEMP(1,0:3,1) = Z_TEMP(1,0:3,1)/(/MW_G_2,MW_G,MW_OTHER,MW_OTHER_2/) + 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) + B1%RHO_F = F_TEMP(1,1,1) + CASE( 3) + IF (SOIB) THEN + RHO_G_2 = RHOP(BC%IIG,BC%JJG,BC%KKG+1) + RHO_OTHER_2 = OM_RHOP(IIO,JJO,KKO-1) + ENDIF + Z_TEMP(1,1,0:3) = (/RHO_OTHER_2,RHO_OTHER,B1%RHO_G,RHO_G_2/) + IF (FLUX_LIMITER_MW_CORRECTION) Z_TEMP(1,1,0:3) = Z_TEMP(1,1,0:3)/(/MW_OTHER_2,MW_OTHER,MW_G,MW_G_2/) + U_TEMP(1,1,1) = WW(BC%II,BC%JJ,BC%KK) + CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER) + B1%RHO_F = F_TEMP(1,1,1) + CASE(-3) + IF (SOIB) THEN + RHO_G_2 = RHOP(BC%IIG,BC%JJG,BC%KKG-1) + RHO_OTHER_2 = OM_RHOP(IIO,JJO,KKO+1) + ENDIF + Z_TEMP(1,1,0:3) = (/RHO_G_2,B1%RHO_G,RHO_OTHER,RHO_OTHER_2/) + IF (FLUX_LIMITER_MW_CORRECTION) Z_TEMP(1,1,0:3) = Z_TEMP(1,1,0:3)/(/MW_G_2,MW_G,MW_OTHER,MW_OTHER_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) + B1%RHO_F = F_TEMP(1,1,1) + END SELECT + + ! 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 (SOIB) 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 = OM_RHOP(IIO-1,JJO,KKO)*OM_ZZP(IIO-1,JJO,KKO,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) + CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER) + RHO_ZZ_F(N) = F_TEMP(1,1,1) + PBAR_F = PBAR_P(BC%KKG,B1%PRESSURE_ZONE) + CASE(-1) + IF (SOIB) 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 = OM_RHOP(IIO+1,JJO,KKO)*OM_ZZP(IIO+1,JJO,KKO,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) + CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER) + RHO_ZZ_F(N) = F_TEMP(1,1,1) + PBAR_F = PBAR_P(BC%KKG,B1%PRESSURE_ZONE) + CASE( 2) + IF (SOIB) 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 = OM_RHOP(IIO,JJO-1,KKO)*OM_ZZP(IIO,JJO-1,KKO,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) + CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER) + RHO_ZZ_F(N) = F_TEMP(1,1,1) + PBAR_F = PBAR_P(BC%KKG,B1%PRESSURE_ZONE) + CASE(-2) + IF (SOIB) 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 = OM_RHOP(IIO,JJO+1,KKO)*OM_ZZP(IIO,JJO+1,KKO,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) + CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER) + RHO_ZZ_F(N) = F_TEMP(1,1,1) + PBAR_F = PBAR_P(BC%KKG,B1%PRESSURE_ZONE) + CASE( 3) + IF (SOIB) 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 = OM_RHOP(IIO,JJO,KKO-1)*OM_ZZP(IIO,JJO,KKO-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) + CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER) + 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 (SOIB) 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 = OM_RHOP(IIO,JJO,KKO+1)*OM_ZZP(IIO,JJO,KKO+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) + CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER) + 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 + + ! Correction for variable molecular weights + + IF (FLUX_LIMITER_MW_CORRECTION) THEN + N=MAXLOC(RHO_ZZ_F(1:N_TRACKED_SPECIES),1) + MW_G = SPECIES_MIXTURE(N)%MW + RHO_ZZ_F(N) = MW_G*MAX( 0._EB, B1%RHO_F & + - SUM(RHO_ZZ_F(1:(N-1))/SPECIES_MIXTURE(1:(N-1))%MW) & + - SUM(RHO_ZZ_F((N+1):N_TRACKED_SPECIES)/SPECIES_MIXTURE((N+1):N_TRACKED_SPECIES)%MW) ) + B1%RHO_F = SUM(RHO_ZZ_F(1:N_TRACKED_SPECIES)) + ENDIF + + ! 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)) @@ -734,29 +931,11 @@ 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