Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
153 changes: 0 additions & 153 deletions Source/func.f90
Original file line number Diff line number Diff line change
Expand Up @@ -1451,159 +1451,6 @@ END FUNCTION MP5
END SUBROUTINE GET_SCALAR_FACE_VALUE


SUBROUTINE GET_SCALAR_FACE_COEF(A,U,BF,I1,I2,J1,J2,K1,K2,IOR,LIMITER)

USE GLOBAL_CONSTANTS, ONLY: CENTRAL_LIMITER,GODUNOV_LIMITER,SUPERBEE_LIMITER,CHARM_LIMITER

REAL(EB), INTENT(IN) :: A(0:,0:,0:),U(-1:,-1:,-1:)
REAL(EB), INTENT(INOUT) :: BF(0:,0:,0:)
INTEGER, INTENT(IN) :: LIMITER,I1,I2,J1,J2,K1,K2,IOR
REAL(EB) :: R,B,DU_UP,DU_LOC
INTEGER :: I,J,K,IM1,JM1,KM1,IP1,JP1,KP1,IP2,JP2,KP2

SELECT CASE(LIMITER)
CASE(GODUNOV_LIMITER,CENTRAL_LIMITER)
BF=0._EB
RETURN
END SELECT

SELECT CASE(IOR)
CASE(1) ; IM1=-1 ; JM1= 0 ; KM1= 0 ; IP1=1 ; JP1=0 ; KP1=0 ; IP2=2 ; JP2=0 ; KP2=0
CASE(2) ; IM1= 0 ; JM1=-1 ; KM1= 0 ; IP1=0 ; JP1=1 ; KP1=0 ; IP2=0 ; JP2=2 ; KP2=0
CASE(3) ; IM1= 0 ; JM1= 0 ; KM1=-1 ; IP1=0 ; JP1=0 ; KP1=1 ; IP2=0 ; JP2=0 ; KP2=2
END SELECT

!$OMP PARALLEL DEFAULT(NONE) IF(I2>1) &
!$OMP& SHARED(U,A,BF,I1,I2,J1,J2,K1,K2,IP1,JP1,KP1,IM1,JM1,KM1,IP2,JP2,KP2,LIMITER) &
!$OMP& PRIVATE(I,J,K,DU_UP,DU_LOC,B,R)

SELECT CASE (LIMITER)

CASE (SUPERBEE_LIMITER)
!$OMP DO COLLAPSE(3) SCHEDULE(STATIC)
DO K=K1,K2
DO J=J1,J2
DO I=I1,I2
DU_LOC = U(I+IP1,J+JP1,K+KP1) - U(I,J,K)
IF (A(I,J,K) > 0._EB) THEN
DU_UP = U(I,J,K) - U(I+IM1,J+JM1,K+KM1)
ELSE
DU_UP = U(I+IP2,J+JP2,K+KP2) - U(I+IP1,J+JP1,K+KP1)
ENDIF
R = 0._EB ; B = 0._EB
IF (ABS(DU_LOC) > TWO_EPSILON_EB) R = DU_UP/DU_LOC
IF (R > TWO_EPSILON_EB) B = MAX(0._EB, MIN(2._EB*R,1._EB), MIN(R,2._EB))
BF(I,J,K) = MAX(0._EB, MIN(B, BF(I,J,K)))
ENDDO
ENDDO
ENDDO
!$OMP END DO

CASE (CHARM_LIMITER)
!$OMP DO COLLAPSE(3) SCHEDULE(STATIC)
DO K=K1,K2
DO J=J1,J2
DO I=I1,I2
DU_LOC = U(I+IP1,J+JP1,K+KP1) - U(I,J,K)
IF (A(I,J,K) > 0._EB) THEN
DU_UP = U(I,J,K) - U(I+IM1,J+JM1,K+KM1)
ELSE
DU_UP = U(I+IP2,J+JP2,K+KP2) - U(I+IP1,J+JP1,K+KP1)
ENDIF
R = 0._EB ; B = 0._EB
IF (ABS(DU_UP) > TWO_EPSILON_EB) R = DU_LOC/DU_UP
IF (R > TWO_EPSILON_EB) B = R*(3._EB*R+1._EB)/((R+1._EB)**2)
BF(I,J,K) = MAX(0._EB, MIN(B, BF(I,J,K)))
ENDDO
ENDDO
ENDDO
!$OMP END DO

END SELECT
!$OMP END PARALLEL

END SUBROUTINE GET_SCALAR_FACE_COEF


SUBROUTINE GET_SCALAR_FACE_VALUE_NEW(A,U,F,BF,I1,I2,J1,J2,K1,K2,IOR,LIMITER)

USE GLOBAL_CONSTANTS, ONLY: CENTRAL_LIMITER,GODUNOV_LIMITER,SUPERBEE_LIMITER,CHARM_LIMITER

REAL(EB), INTENT(IN) :: A(0:,0:,0:),U(-1:,-1:,-1:),BF(0:,0:,0:)
REAL(EB), INTENT(OUT) :: F(0:,0:,0:)
INTEGER, INTENT(IN) :: LIMITER,I1,I2,J1,J2,K1,K2,IOR
REAL(EB) :: DU_UP,DU_LOC
INTEGER :: I,J,K,IM1,JM1,KM1,IP1,JP1,KP1,IP2,JP2,KP2

SELECT CASE(IOR)
CASE(1) ; IM1=-1 ; JM1= 0 ; KM1= 0 ; IP1=1 ; JP1=0 ; KP1=0 ; IP2=2 ; JP2=0 ; KP2=0
CASE(2) ; IM1= 0 ; JM1=-1 ; KM1= 0 ; IP1=0 ; JP1=1 ; KP1=0 ; IP2=0 ; JP2=2 ; KP2=0
CASE(3) ; IM1= 0 ; JM1= 0 ; KM1=-1 ; IP1=0 ; JP1=0 ; KP1=1 ; IP2=0 ; JP2=0 ; KP2=2
END SELECT

!$OMP PARALLEL IF(I2>1)
SELECT CASE(LIMITER)
CASE(CENTRAL_LIMITER) ! central differencing
!$OMP DO SCHEDULE(STATIC)
DO K=K1,K2
DO J=J1,J2
DO I=I1,I2
F(I,J,K) = 0.5_EB*(U(I,J,K) + U(I+IP1,J+JP1,K+KP1))
ENDDO
ENDDO
ENDDO
!$OMP END DO
CASE(GODUNOV_LIMITER) ! first-order upwinding
!$OMP DO SCHEDULE(STATIC)
DO K=K1,K2
DO J=J1,J2
DO I=I1,I2
IF (A(I,J,K)>0._EB) THEN
F(I,J,K) = U(I,J,K)
ELSE
F(I,J,K) = U(I+IP1,J+JP1,K+KP1)
ENDIF
ENDDO
ENDDO
ENDDO
!$OMP END DO
CASE(SUPERBEE_LIMITER) ! SUPERBEE, Roe (1986)
!$OMP DO SCHEDULE(STATIC) PRIVATE(DU_LOC)
DO K=K1,K2
DO J=J1,J2
DO I=I1,I2
DU_LOC = U(I+IP1,J+JP1,K+KP1)-U(I,J,K)
IF (A(I,J,K)>0._EB) THEN
F(I,J,K) = U(I,J,K) + 0.5_EB*BF(I,J,K)*DU_LOC
ELSE
F(I,J,K) = U(I+IP1,J+JP1,K+KP1) - 0.5_EB*BF(I,J,K)*DU_LOC
ENDIF
ENDDO
ENDDO
ENDDO
!$OMP END DO
CASE(CHARM_LIMITER) ! CHARM
!$OMP DO SCHEDULE(STATIC) PRIVATE(DU_UP)
DO K=K1,K2
DO J=J1,J2
DO I=I1,I2
IF (A(I,J,K)>0._EB) THEN
DU_UP = U(I,J,K)-U(I+IM1,J+JM1,K+KM1)
F(I,J,K) = U(I,J,K) + 0.5_EB*BF(I,J,K)*DU_UP
ELSE
DU_UP = U(I+IP2,J+JP2,K+KP2)-U(I+IP1,J+JP1,K+KP1)
F(I,J,K) = U(I+IP1,J+JP1,K+KP1) - 0.5_EB*BF(I,J,K)*DU_UP
ENDIF
ENDDO
ENDDO
ENDDO
!$OMP END DO
END SELECT
!$OMP END PARALLEL

END SUBROUTINE GET_SCALAR_FACE_VALUE_NEW


!> \brief Random fluctuation, Theta'(t+dt) = R^2*Theta'(t) + Normal(0,sqrt(1-R^2)*SIGMA) ; R = exp(-dt/TAU)
!> \param SIGMA Standard deviation of time series
!> \param TAU Time scale or period of the time function (s)
Expand Down
Loading
Loading