diff --git a/Source/ccib.f90 b/Source/ccib.f90 index 54651c46473..4a470fad4eb 100644 --- a/Source/ccib.f90 +++ b/Source/ccib.f90 @@ -8130,13 +8130,13 @@ SUBROUTINE SET_EXIMRHOZZLIM_3D(NM,N) ! Local Variables: INTEGER :: I,J,K,X1AXIS -REAL(EB), POINTER, DIMENSION(:,:,:) :: FX_ZZ=>NULL(),FY_ZZ=>NULL(),FZ_ZZ=>NULL() +REAL(EB), POINTER, DIMENSION(:,:,:,:) :: FX_ZZ=>NULL(),FY_ZZ=>NULL(),FZ_ZZ=>NULL() INTEGER :: IFACE,IW TYPE(CC_REGFACEZ_TYPE), POINTER, DIMENSION(:) :: REGFACE_Z=>NULL() -FX_ZZ=>WORK2 -FY_ZZ=>WORK3 -FZ_ZZ=>WORK4 +FX_ZZ=>SWORK1 +FY_ZZ=>SWORK2 +FZ_ZZ=>SWORK3 AXIS_DO : DO X1AXIS = IAXIS,KAXIS SELECT CASE(X1AXIS) @@ -8164,11 +8164,11 @@ SUBROUTINE SET_EXIMRHOZZLIM_3D(NM,N) ELSE SELECT CASE(X1AXIS) CASE(IAXIS) - REGFACE_Z(IFACE)%FN_ZZ(N) = FX_ZZ(I,J,K) + REGFACE_Z(IFACE)%FN_ZZ(N) = FX_ZZ(I,J,K,N) CASE(JAXIS) - REGFACE_Z(IFACE)%FN_ZZ(N) = FY_ZZ(I,J,K) + REGFACE_Z(IFACE)%FN_ZZ(N) = FY_ZZ(I,J,K,N) CASE(KAXIS) - REGFACE_Z(IFACE)%FN_ZZ(N) = FZ_ZZ(I,J,K) + REGFACE_Z(IFACE)%FN_ZZ(N) = FZ_ZZ(I,J,K,N) END SELECT ENDIF ENDDO IFACE_DO diff --git a/Source/cons.f90 b/Source/cons.f90 index ac3d51441b0..96d5cc4aded 100644 --- a/Source/cons.f90 +++ b/Source/cons.f90 @@ -270,7 +270,7 @@ MODULE GLOBAL_CONSTANTS LOGICAL :: TENSOR_DIFFUSIVITY=.FALSE. !< If true, use experimental tensor diffusivity model for spec and tmp LOGICAL :: OXPYRO_MODEL=.FALSE. !< Flag to use oxidative pyrolysis mass transfer model LOGICAL :: OUTPUT_WALL_QUANTITIES=.FALSE. !< Flag to force call to WALL_MODEL -LOGICAL :: TEST_FLUX_LIMITER_FACE_CORRECTION=.FALSE. +LOGICAL :: FLUX_LIMITER_MW_CORRECTION=.FALSE. INTEGER, ALLOCATABLE, DIMENSION(:) :: CHANGE_TIME_STEP_INDEX !< Flag to indicate if a mesh needs to change time step INTEGER, ALLOCATABLE, DIMENSION(:) :: SETUP_PRESSURE_ZONES_INDEX !< Flag to indicate if a mesh needs to keep searching for ZONEs diff --git a/Source/divg.f90 b/Source/divg.f90 index be1160bc814..da7287f941a 100644 --- a/Source/divg.f90 +++ b/Source/divg.f90 @@ -102,9 +102,9 @@ SUBROUTINE DIVERGENCE_PART_1(T,DT,NM) ! Compute species-related finite difference terms -RHO_D_DZDX => SCALAR_WORK1 -RHO_D_DZDY => SCALAR_WORK2 -RHO_D_DZDZ => SCALAR_WORK3 +RHO_D_DZDX => SWORK1 +RHO_D_DZDY => SWORK2 +RHO_D_DZDZ => SWORK3 ! Save the largest value of the material and thermal diffusion coefficients for use in Von Neumann stability constraint @@ -623,9 +623,11 @@ SUBROUTINE DIVERGENCE_PART_1(T,DT,NM) CONST_GAMMA_IF_2: IF (.NOT.CONSTANT_SPECIFIC_HEAT_RATIO) THEN + CALL SPECIES_ADVECTION_PART_1 ! Compute and store face values of (rho Z_n) + DO N=1,N_TRACKED_SPECIES - CALL SPECIES_ADVECTION(N,U_DOT_DEL_RHO_Z) ! Compute u dot grad rho Z_n + CALL SPECIES_ADVECTION_PART_2(N,U_DOT_DEL_RHO_Z) ! Compute u dot grad rho Z_n SM => SPECIES_MIXTURE(N) !$OMP PARALLEL DO PRIVATE(H_S) SCHEDULE(guided) @@ -950,40 +952,254 @@ SUBROUTINE ENTHALPY_ADVECTION(U_DOT_DEL_RHO_H_S) END SUBROUTINE ENTHALPY_ADVECTION -SUBROUTINE SPECIES_ADVECTION(N,U_DOT_DEL_RHO_Z) +SUBROUTINE SPECIES_ADVECTION_PART_1 USE MATH_FUNCTIONS, ONLY: GET_SCALAR_FACE_VALUE -INTEGER, INTENT(IN) :: N -REAL(EB), POINTER, DIMENSION(:,:,:) :: RHO_Z_P,U_DOT_DEL_RHO_Z -REAL(EB), POINTER, DIMENSION(:,:,:) :: FX_ZZ,FY_ZZ,FZ_ZZ -REAL(EB) :: UN,DU_P,DU_M,DV_P,DV_M,DW_P,DW_M,DU +USE PHYSICAL_FUNCTIONS, ONLY: GET_MOLECULAR_WEIGHT +! INTEGER, INTENT(IN) :: N +REAL(EB), POINTER, DIMENSION(:,:,:) :: RHO_Z_P,RHO_RMW +REAL(EB), POINTER, DIMENSION(:,:,:,:) :: FX_ZZ,FY_ZZ,FZ_ZZ +REAL(EB) :: ZZ_GET(1:N_TRACKED_SPECIES),MW_G +REAL(EB), PARAMETER :: DUMMY=0._EB REAL(EB), DIMENSION(0:3,0:3,0:3) :: Z_TEMP,U_TEMP,F_TEMP -INTEGER :: IC,I,J,K,IW +INTEGER :: I,J,K,IW,N !,II,JJ,KK,IOR,IC,IIG,JJG,KKG TYPE(WALL_TYPE), POINTER :: WC TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1 -FX_ZZ=>WORK2 ; FX_ZZ = 0._EB -FY_ZZ=>WORK3 ; FY_ZZ = 0._EB -FZ_ZZ=>WORK4 ; FZ_ZZ = 0._EB -RHO_Z_P=>WORK6 ; RHO_Z_P = 0._EB -U_DOT_DEL_RHO_Z=>WORK7 ; U_DOT_DEL_RHO_Z = 0._EB +FX_ZZ=>SWORK1 +FY_ZZ=>SWORK2 +FZ_ZZ=>SWORK3 +RHO_Z_P=>WORK6 -!$OMP PARALLEL DO COLLAPSE(3) -DO K=0,KBP1 - DO J=0,JBP1 - DO I=0,IBP1 - RHO_Z_P(I,J,K) = RHOP(I,J,K)*ZZP(I,J,K,N) +! Species face values + +SPECIES_LOOP: DO N=1,N_TOTAL_SCALARS + + !$OMP PARALLEL DO + DO K=0,KBP1 + DO J=0,JBP1 + DO I=0,IBP1 + RHO_Z_P(I,J,K) = RHOP(I,J,K)*ZZP(I,J,K,N) + ENDDO ENDDO ENDDO -ENDDO -!$OMP END PARALLEL DO + !$OMP END PARALLEL DO -! Compute scalar face values + ! Compute scalar face values -CALL GET_SCALAR_FACE_VALUE(UU,RHO_Z_P,FX_ZZ,1,IBM1,1,JBAR,1,KBAR,1,I_FLUX_LIMITER) -CALL GET_SCALAR_FACE_VALUE(VV,RHO_Z_P,FY_ZZ,1,IBAR,1,JBM1,1,KBAR,2,I_FLUX_LIMITER) -CALL GET_SCALAR_FACE_VALUE(WW,RHO_Z_P,FZ_ZZ,1,IBAR,1,JBAR,1,KBM1,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) + 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 + +FACE_CORRECTION_IF: IF (FLUX_LIMITER_MW_CORRECTION) THEN + + ! Repeat the above for DENSITY + + RHO_RMW=>WORK6 + + !$OMP PARALLEL DO PRIVATE(ZZ_GET,MW_G) SCHEDULE(STATIC) + DO K=0,KBP1 + DO J=0,JBP1 + DO I=0,IBP1 + 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 + + +SUBROUTINE SPECIES_ADVECTION_PART_2(N,U_DOT_DEL_RHO_Z) + +USE MATH_FUNCTIONS, ONLY: GET_SCALAR_FACE_VALUE +INTEGER, INTENT(IN) :: N +REAL(EB), POINTER, DIMENSION(:,:,:) :: U_DOT_DEL_RHO_Z +REAL(EB), POINTER, DIMENSION(:,:,:,:) :: FX_ZZ,FY_ZZ,FZ_ZZ +REAL(EB) :: UN,DU_P,DU_M,DV_P,DV_M,DW_P,DW_M,DU +INTEGER :: IC,I,J,K,IW +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 +U_DOT_DEL_RHO_Z=>WORK7 ; U_DOT_DEL_RHO_Z = 0._EB WALL_LOOP: DO IW=1,N_EXTERNAL_WALL_CELLS+N_INTERNAL_WALL_CELLS @@ -992,63 +1208,6 @@ SUBROUTINE SPECIES_ADVECTION(N,U_DOT_DEL_RHO_Z) B1 => BOUNDARY_PROP1(WC%B1_INDEX) BC => BOUNDARY_COORD(WC%BC_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) = 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) = 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) = 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) = 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) = 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) = F_TEMP(1,1,1) - ENDIF - END SELECT OFF_WALL_SELECT_2 - - ENDIF OFF_WALL_IF_2 - ! Correct flux terms at the boundary SELECT CASE(WC%BOUNDARY_TYPE) @@ -1068,7 +1227,7 @@ SUBROUTINE SPECIES_ADVECTION(N,U_DOT_DEL_RHO_Z) UN = UVW_SAVE(IW) END SELECT - DU = (B1%RHO_F*B1%ZZ_F(N) - RHO_Z_P(BC%IIG,BC%JJG,BC%KKG))*UN + DU = (B1%RHO_F*B1%ZZ_F(N) - RHOP(BC%IIG,BC%JJG,BC%KKG)*ZZP(BC%IIG,BC%JJG,BC%KKG,N))*UN U_DOT_DEL_RHO_Z(BC%IIG,BC%JJG,BC%KKG) = U_DOT_DEL_RHO_Z(BC%IIG,BC%JJG,BC%KKG) - SIGN(1._EB,REAL(BC%IOR,EB))*DU*B1%RDN ENDDO WALL_LOOP @@ -1085,19 +1244,19 @@ SUBROUTINE SPECIES_ADVECTION(N,U_DOT_DEL_RHO_Z) DV_M = 0._EB DW_P = 0._EB DW_M = 0._EB - IF (CELL(IC)%WALL_INDEX( 1)==0) DU_P = (FX_ZZ(I,J,K) - RHO_Z_P(I,J,K))*UU(I,J,K) - IF (CELL(IC)%WALL_INDEX(-1)==0) DU_M = (FX_ZZ(I-1,J,K) - RHO_Z_P(I,J,K))*UU(I-1,J,K) - IF (CELL(IC)%WALL_INDEX( 2)==0) DV_P = (FY_ZZ(I,J,K) - RHO_Z_P(I,J,K))*VV(I,J,K) - IF (CELL(IC)%WALL_INDEX(-2)==0) DV_M = (FY_ZZ(I,J-1,K) - RHO_Z_P(I,J,K))*VV(I,J-1,K) - IF (CELL(IC)%WALL_INDEX( 3)==0) DW_P = (FZ_ZZ(I,J,K) - RHO_Z_P(I,J,K))*WW(I,J,K) - IF (CELL(IC)%WALL_INDEX(-3)==0) DW_M = (FZ_ZZ(I,J,K-1) - RHO_Z_P(I,J,K))*WW(I,J,K-1) + IF (CELL(IC)%WALL_INDEX( 1)==0) DU_P = (FX_ZZ(I,J,K,N) - RHOP(I,J,K)*ZZP(I,J,K,N))*UU(I,J,K) + IF (CELL(IC)%WALL_INDEX(-1)==0) DU_M = (FX_ZZ(I-1,J,K,N) - RHOP(I,J,K)*ZZP(I,J,K,N))*UU(I-1,J,K) + IF (CELL(IC)%WALL_INDEX( 2)==0) DV_P = (FY_ZZ(I,J,K,N) - RHOP(I,J,K)*ZZP(I,J,K,N))*VV(I,J,K) + IF (CELL(IC)%WALL_INDEX(-2)==0) DV_M = (FY_ZZ(I,J-1,K,N) - RHOP(I,J,K)*ZZP(I,J,K,N))*VV(I,J-1,K) + IF (CELL(IC)%WALL_INDEX( 3)==0) DW_P = (FZ_ZZ(I,J,K,N) - RHOP(I,J,K)*ZZP(I,J,K,N))*WW(I,J,K) + IF (CELL(IC)%WALL_INDEX(-3)==0) DW_M = (FZ_ZZ(I,J,K-1,N) - RHOP(I,J,K)*ZZP(I,J,K,N))*WW(I,J,K-1) U_DOT_DEL_RHO_Z(I,J,K) = U_DOT_DEL_RHO_Z(I,J,K) + (DU_P-DU_M)*RDX(I) + (DV_P-DV_M)*RDY(J) + (DW_P-DW_M)*RDZ(K) ENDDO ENDDO ENDDO !$OMP END PARALLEL DO -END SUBROUTINE SPECIES_ADVECTION +END SUBROUTINE SPECIES_ADVECTION_PART_2 SUBROUTINE MERGE_PRESSURE_ZONES diff --git a/Source/fire.f90 b/Source/fire.f90 index 7fb2eeb5933..0729337e8b5 100644 --- a/Source/fire.f90 +++ b/Source/fire.f90 @@ -1616,7 +1616,7 @@ SUBROUTINE CONDENSATION_EVAPORATION(DT,NM) CALL POINT_TO_MESH(NM) -ZZ_INTERIM=> SCALAR_WORK1 +ZZ_INTERIM=> SWORK1 ZZ_INTERIM = ZZ RHO_INTERIM => WORK1 RHO_INTERIM = RHO diff --git a/Source/init.f90 b/Source/init.f90 index a091f434427..bf416922314 100644 --- a/Source/init.f90 +++ b/Source/init.f90 @@ -188,18 +188,15 @@ SUBROUTINE INITIALIZE_MESH_VARIABLES_1(DT,NM) ! Allocate scalar face values ! Required for cell face density correction for multicomponent mixtures -IF (TEST_FLUX_LIMITER_FACE_CORRECTION .AND. N_TRACKED_SPECIES>2) THEN +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) -ALLOCATE( M%FY(0:IBP1,0:JBP1,0:KBP1,N_LOWER_SCALARS:N_TOTAL_SCALARS),STAT=IZERO) -CALL ChkMemErr('INIT','FY',IZERO) -ALLOCATE( M%FZ(0:IBP1,0:JBP1,0:KBP1,N_LOWER_SCALARS:N_TOTAL_SCALARS),STAT=IZERO) -CALL ChkMemErr('INIT','FZ',IZERO) +ALLOCATE( M%FX(0:IBP1,0:JBP1,0:KBP1,N_LOWER_SCALARS:N_TOTAL_SCALARS),STAT=IZERO); CALL ChkMemErr('INIT','FX',IZERO) +ALLOCATE( M%FY(0:IBP1,0:JBP1,0:KBP1,N_LOWER_SCALARS:N_TOTAL_SCALARS),STAT=IZERO); CALL ChkMemErr('INIT','FY',IZERO) +ALLOCATE( M%FZ(0:IBP1,0:JBP1,0:KBP1,N_LOWER_SCALARS:N_TOTAL_SCALARS),STAT=IZERO); CALL ChkMemErr('INIT','FZ',IZERO) M%FX = 0._EB M%FY = 0._EB M%FZ = 0._EB @@ -306,16 +303,16 @@ SUBROUTINE INITIALIZE_MESH_VARIABLES_1(DT,NM) ALLOCATE(M%WORK8(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO) ; CALL ChkMemErr('INIT','WORK8',IZERO) 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%SCALAR_WORK1(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO) ; CALL ChkMemErr('INIT','SCALAR_WORK1',IZERO) -ALLOCATE(M%SCALAR_WORK2(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO) ; CALL ChkMemErr('INIT','SCALAR_WORK2',IZERO) -ALLOCATE(M%SCALAR_WORK3(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO) ; CALL ChkMemErr('INIT','SCALAR_WORK3',IZERO) -ALLOCATE(M%SCALAR_WORK4(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO) ; CALL ChkMemErr('INIT','SCALAR_WORK4',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,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%SCALAR_WORK1=0._EB -M%SCALAR_WORK2=0._EB -M%SCALAR_WORK3=0._EB -M%SCALAR_WORK4=0._EB +M%SWORK1=0._EB +M%SWORK2=0._EB +M%SWORK3=0._EB +M%SWORK4=0._EB IF (STRATIFICATION) THEN diff --git a/Source/mass.f90 b/Source/mass.f90 index 09747e2239f..35816c3dc0c 100644 --- a/Source/mass.f90 +++ b/Source/mass.f90 @@ -178,7 +178,7 @@ SUBROUTINE MASS_FINITE_DIFFERENCES(NM) ENDDO SPECIES_LOOP -FACE_CORRECTION_IF: IF (TEST_FLUX_LIMITER_FACE_CORRECTION) THEN +FACE_CORRECTION_IF: IF (FLUX_LIMITER_MW_CORRECTION) THEN ! Repeat the above for DENSITY @@ -375,7 +375,7 @@ SUBROUTINE DENSITY(T,DT,NM) UU=>WORK1 VV=>WORK2 WW=>WORK3 -DEL_RHO_D_DEL_Z__0=>SCALAR_WORK4 +DEL_RHO_D_DEL_Z__0=>SWORK4 PREDICTOR_STEP: SELECT CASE (PREDICTOR) diff --git a/Source/mesh.f90 b/Source/mesh.f90 index 85663fe0007..ad9388f7008 100644 --- a/Source/mesh.f90 +++ b/Source/mesh.f90 @@ -107,10 +107,10 @@ MODULE MESH_VARIABLES ! Work arrays - REAL(EB), ALLOCATABLE, DIMENSION(:,:,:,:) :: SCALAR_WORK1 - REAL(EB), ALLOCATABLE, DIMENSION(:,:,:,:) :: SCALAR_WORK2 - REAL(EB), ALLOCATABLE, DIMENSION(:,:,:,:) :: SCALAR_WORK3 - REAL(EB), ALLOCATABLE, DIMENSION(:,:,:,:) :: SCALAR_WORK4 + REAL(EB), ALLOCATABLE, DIMENSION(:,:,:,:) :: SWORK1 + REAL(EB), ALLOCATABLE, DIMENSION(:,:,:,:) :: SWORK2 + REAL(EB), ALLOCATABLE, DIMENSION(:,:,:,:) :: SWORK3 + REAL(EB), ALLOCATABLE, DIMENSION(:,:,:,:) :: SWORK4 REAL(EB), ALLOCATABLE, DIMENSION(:,:,:) :: WORK1,WORK2,WORK3,WORK4,WORK5,WORK6,WORK7,WORK8,WORK9 INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: IWORK1 REAL(EB), ALLOCATABLE, DIMENSION(:,:,:) :: PWORK1,PWORK2,PWORK3,PWORK4 @@ -349,7 +349,7 @@ MODULE MESH_POINTERS MU,MU_DNS,TMP,Q,KAPPA_GAS,CHI_R,QR,QR_W,UII,RSUM,D_SOURCE, & 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, & - SCALAR_WORK1,SCALAR_WORK2,SCALAR_WORK3,SCALAR_WORK4, & + SWORK1,SWORK2,SWORK3,SWORK4, & Q_REAC,AVG_DROP_DEN,AVG_DROP_TMP,AVG_DROP_RAD,AVG_DROP_AREA,M_DOT_PPP, & ADV_FX,ADV_FY,ADV_FZ,DIF_FX,DIF_FY,DIF_FZ,DIF_FXS,DIF_FYS,DIF_FZS INTEGER, POINTER, DIMENSION(:,:) :: CHEM_ACTIVE_CELLS @@ -551,10 +551,10 @@ SUBROUTINE POINT_TO_MESH(NM) DIF_FXS=>M%DIF_FXS DIF_FYS=>M%DIF_FYS DIF_FZS=>M%DIF_FZS -SCALAR_WORK1=>M%SCALAR_WORK1 -SCALAR_WORK2=>M%SCALAR_WORK2 -SCALAR_WORK3=>M%SCALAR_WORK3 -SCALAR_WORK4=>M%SCALAR_WORK4 +SWORK1=>M%SWORK1 +SWORK2=>M%SWORK2 +SWORK3=>M%SWORK3 +SWORK4=>M%SWORK4 WORK=>M%WORK LSAVE=>M%LSAVE LWORK=>M%LWORK diff --git a/Source/part.f90 b/Source/part.f90 index 9c1f60f7913..527451a667f 100644 --- a/Source/part.f90 +++ b/Source/part.f90 @@ -3634,7 +3634,7 @@ SUBROUTINE PARTICLE_MASS_ENERGY_TRANSFER(T,DT,NM) RHO_INTERIM => WORK1 ; RHO_INTERIM = RHO TMP_INTERIM => WORK2 ; TMP_INTERIM = TMP - ZZ_INTERIM => SCALAR_WORK1 ; ZZ_INTERIM = ZZ + ZZ_INTERIM => SWORK1 ; ZZ_INTERIM = ZZ ENDIF diff --git a/Source/read.f90 b/Source/read.f90 index 5a2d97b1d00..b4f3118a283 100644 --- a/Source/read.f90 +++ b/Source/read.f90 @@ -1745,7 +1745,7 @@ SUBROUTINE READ_MISC RAMP_UX,RAMP_UY,RAMP_UZ,RAMP_VX,RAMP_VY,RAMP_VZ,RAMP_WX,RAMP_WY,RAMP_WZ,& RESTART,RESTART_CHID,SC,& RND_SEED,SIMULATION_MODE,SMOKE3D_16,SMOKE_ALBEDO,SOLID_PHASE_ONLY,SOOT_DENSITY,SOOT_OXIDATION,& - TAU_DEFAULT,TENSOR_DIFFUSIVITY,TERRAIN_IMAGE,TEST_FLUX_LIMITER_FACE_CORRECTION,TEXTURE_ORIGIN,& + TAU_DEFAULT,TENSOR_DIFFUSIVITY,TERRAIN_IMAGE,FLUX_LIMITER_MW_CORRECTION,TEXTURE_ORIGIN,& THERMOPHORETIC_DEPOSITION,THERMOPHORETIC_SETTLING,THICKEN_OBSTRUCTIONS,& TMPA,TURBULENCE_MODEL,TURBULENT_DEPOSITION,UVW_FILE,& VERBOSE,VISIBILITY_FACTOR,VN_MAX,VN_MIN,Y_CO2_INFTY,Y_O2_INFTY,& @@ -2823,7 +2823,12 @@ SUBROUTINE READ_SPEC N_TRACKED_SPECIES = N_LUMPED + N_COPY + N_TOTAL_BINS + N_PREDEFINED_SMIX + N_CONDENSATION N_TOTAL_SCALARS = N_TRACKED_SPECIES + N_PASSIVE_SCALARS -! Allocate the primitive species array. +! Only allow use of flux limiter molecular weight correction if N_TRACKED_SPECIES>2 + +IF (N_TRACKED_SPECIES<=2) FLUX_LIMITER_MW_CORRECTION=.FALSE. + +! Allocate the primitive species array + ALLOCATE(SPECIES(N_SPECIES),STAT=IZERO) CALL ChkMemErr('READ','SPECIES',IZERO) ALLOCATE(Y_INDEX(N_SPECIES)) diff --git a/Source/turb.f90 b/Source/turb.f90 index b05bfc00dc7..6c982c72589 100644 --- a/Source/turb.f90 +++ b/Source/turb.f90 @@ -1651,9 +1651,9 @@ SUBROUTINE TENSOR_DIFFUSIVITY_MODEL(NM,OPT_N) ! SGS scalar flux ! CAUTION: The flux arrays must point to the same work arrays used in DIVERGENCE_PART_1 ! Note: Do not reinitialize! RHO_D_DZDX, etc., already store molecular diffusive flux - RHO_D_DZDX=>SCALAR_WORK1 - RHO_D_DZDY=>SCALAR_WORK2 - RHO_D_DZDZ=>SCALAR_WORK3 + RHO_D_DZDX=>SWORK1 + RHO_D_DZDY=>SWORK2 + RHO_D_DZDZ=>SWORK3 IF (PREDICTOR) THEN UU=>U diff --git a/Source/wall.f90 b/Source/wall.f90 index 4e959347c61..689fe356957 100644 --- a/Source/wall.f90 +++ b/Source/wall.f90 @@ -653,7 +653,7 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I ! Correction for variable molecular weights - IF (TEST_FLUX_LIMITER_FACE_CORRECTION) THEN + IF (FLUX_LIMITER_MW_CORRECTION) THEN 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 (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN @@ -692,7 +692,7 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I 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 (TEST_FLUX_LIMITER_FACE_CORRECTION) Z_TEMP(0:3,1,1) = Z_TEMP(0:3,1,1)/(/MW_OTHER_2,MW_OTHER,MW_G,MW_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) @@ -702,7 +702,7 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I 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 (TEST_FLUX_LIMITER_FACE_CORRECTION) Z_TEMP(0:3,1,1) = Z_TEMP(0:3,1,1)/(/MW_G_2,MW_G,MW_OTHER,MW_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) @@ -712,7 +712,7 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I 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 (TEST_FLUX_LIMITER_FACE_CORRECTION) Z_TEMP(1,0:3,1) = Z_TEMP(1,0:3,1)/(/MW_OTHER_2,MW_OTHER,MW_G,MW_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) @@ -722,7 +722,7 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I 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 (TEST_FLUX_LIMITER_FACE_CORRECTION) Z_TEMP(1,0:3,1) = Z_TEMP(1,0:3,1)/(/MW_G_2,MW_G,MW_OTHER,MW_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) @@ -732,7 +732,7 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I 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 (TEST_FLUX_LIMITER_FACE_CORRECTION) Z_TEMP(1,1,0:3) = Z_TEMP(1,1,0:3)/(/MW_OTHER_2,MW_OTHER,MW_G,MW_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) @@ -742,7 +742,7 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I 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 (TEST_FLUX_LIMITER_FACE_CORRECTION) Z_TEMP(1,1,0:3) = Z_TEMP(1,1,0:3)/(/MW_G_2,MW_G,MW_OTHER,MW_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) @@ -825,7 +825,7 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I ! Correction for variable molecular weights - IF (TEST_FLUX_LIMITER_FACE_CORRECTION) THEN + 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 &