Merge pull request #14803 from mcgratta/master

FDS Source: Second-order accuracy for LS at interpolated boundaries

Author: mcgratta

Source/main.f90

@@ -3239,14 +3239,21 @@ SUBROUTINE MESH_EXCHANGE(CODE)
             PHI_LS_P => M%PHI_LS
          ENDIF
          IF (RNODE/=SNODE) THEN
-            NQT2 = 4
+            NQT2 = 5
             PACK_REAL_SEND_PKG14: DO LL=1,M3%NIC_S
-               II1 = M3%IIO_S(LL)
-               JJ1 = M3%JJO_S(LL)
+               II1 = M3%IIO_S(LL) ; II2 = II1
+               JJ1 = M3%JJO_S(LL) ; JJ2 = JJ1
+               SELECT CASE(M3%IOR_S(LL))
+                  CASE(-1) ; II2=II1+1
+                  CASE( 1) ; II2=II1-1
+                  CASE(-2) ; JJ2=JJ1+1
+                  CASE( 2) ; JJ2=JJ1-1
+               END SELECT
                M3%REAL_SEND_PKG14(NQT2*(LL-1)+1) = PHI_LS_P(II1,JJ1)
-               M3%REAL_SEND_PKG14(NQT2*(LL-1)+2) = M%U_LS(II1,JJ1)
-               M3%REAL_SEND_PKG14(NQT2*(LL-1)+3) = M%V_LS(II1,JJ1)
-               M3%REAL_SEND_PKG14(NQT2*(LL-1)+4) = M%Z_LS(II1,JJ1)
+               M3%REAL_SEND_PKG14(NQT2*(LL-1)+2) = PHI_LS_P(II2,JJ2)
+               M3%REAL_SEND_PKG14(NQT2*(LL-1)+3) = M%U_LS(II1,JJ1)
+               M3%REAL_SEND_PKG14(NQT2*(LL-1)+4) = M%V_LS(II1,JJ1)
+               M3%REAL_SEND_PKG14(NQT2*(LL-1)+5) = M%Z_LS(II1,JJ1)
             ENDDO PACK_REAL_SEND_PKG14
          ELSE
             M2=>MESHES(NOM)%OMESH(NM)
@@ -3540,19 +3547,26 @@ SUBROUTINE MESH_EXCHANGE(CODE)
       ENDIF IF_RECEIVE_PARTICLES
 
       IF (CODE==14 .AND. M2%NIC_R>0 .AND. RNODE/=SNODE) THEN
-            NQT2 = 4
+            NQT2 = 5
             IF (PREDICTOR) THEN
                PHI_LS_P => M2%PHI1_LS
             ELSE
                PHI_LS_P => M2%PHI_LS
             ENDIF
             UNPACK_REAL_RECV_PKG14: DO LL=1,M2%NIC_R
-               II1 = M2%IIO_R(LL)
-               JJ1 = M2%JJO_R(LL)
+               II1 = M2%IIO_R(LL) ; II2 = II1
+               JJ1 = M2%JJO_R(LL) ; JJ2 = JJ1
+               SELECT CASE(M2%IOR_R(LL))
+                  CASE(-1) ; II2=II1+1
+                  CASE( 1) ; II2=II1-1
+                  CASE(-2) ; JJ2=JJ1+1
+                  CASE( 2) ; JJ2=JJ1-1
+               END SELECT
                PHI_LS_P(II1,JJ1) = M2%REAL_RECV_PKG14(NQT2*(LL-1)+1)
-               M2%U_LS(II1,JJ1)  = M2%REAL_RECV_PKG14(NQT2*(LL-1)+2)
-               M2%V_LS(II1,JJ1)  = M2%REAL_RECV_PKG14(NQT2*(LL-1)+3)
-               M2%Z_LS(II1,JJ1)  = M2%REAL_RECV_PKG14(NQT2*(LL-1)+4)
+               PHI_LS_P(II2,JJ2) = M2%REAL_RECV_PKG14(NQT2*(LL-1)+2)
+               M2%U_LS(II1,JJ1)  = M2%REAL_RECV_PKG14(NQT2*(LL-1)+3)
+               M2%V_LS(II1,JJ1)  = M2%REAL_RECV_PKG14(NQT2*(LL-1)+4)
+               M2%Z_LS(II1,JJ1)  = M2%REAL_RECV_PKG14(NQT2*(LL-1)+5)
             ENDDO UNPACK_REAL_RECV_PKG14
       ENDIF
 

Source/vege.f90

@@ -62,8 +62,10 @@ SUBROUTINE INITIALIZE_LEVEL_SET_FIRESPREAD_1(NM)
 
 ! Level set values (Phi). PHI1_LS is the first-order accurate estimate at the next time step.
 
-ALLOCATE(M%PHI_LS(0:IBP1,0:JBP1)) ; CALL ChkMemErr('VEGE:LEVEL SET','PHI_LS',IZERO)   ; PHI_LS  => M%PHI_LS  ; PHI_LS = PHI_LS_MIN
-ALLOCATE(M%PHI1_LS(0:IBP1,0:JBP1)); CALL ChkMemErr('VEGE:LEVEL SET','PHI1_LS',IZERO)  ; PHI1_LS => M%PHI1_LS ; PHI1_LS = PHI_LS_MIN
+ALLOCATE(M%PHI_LS(-1:IBAR+2,-1:JBAR+2)) ; CALL ChkMemErr('VEGE:LEVEL SET','PHI_LS',IZERO)  ; PHI_LS => M%PHI_LS   
+PHI_LS = PHI_LS_MIN
+ALLOCATE(M%PHI1_LS(-1:IBAR+2,-1:JBAR+2)); CALL ChkMemErr('VEGE:LEVEL SET','PHI1_LS',IZERO) ; PHI1_LS => M%PHI1_LS 
+PHI1_LS = PHI_LS_MIN
 
 ! Wind speed components in the center of the first gas phsae cell above the ground.
 
@@ -682,8 +684,8 @@ SUBROUTINE GET_BOUNDARY_VALUES
 SUBROUTINE FILL_BOUNDARY_VALUES
 
 USE COMPLEX_GEOMETRY, ONLY : CC_CGSC,CC_SOLID,CC_CUTCFE
-INTEGER :: IW,IIO,JJO,N_INT_CELLS,NOM,IC
-REAL(EB) :: PHI_LS_OTHER,U_LS_OTHER,V_LS_OTHER,Z_LS_OTHER
+INTEGER :: IW,IIO,JJO,IIO_2,JJO_2,N_INT_CELLS,NOM,IC
+REAL(EB) :: PHI_LS_OTHER,PHI_LS_OTHER_2,U_LS_OTHER,V_LS_OTHER,Z_LS_OTHER
 TYPE (EXTERNAL_WALL_TYPE), POINTER :: EWC
 LOGICAL :: SOLID_CELL
 
@@ -707,16 +709,27 @@ SUBROUTINE FILL_BOUNDARY_VALUES
    RETURN
 ENDIF
 
-PHI_LS_OTHER = 0._EB
+PHI_LS_OTHER   = 0._EB
+PHI_LS_OTHER_2 = 0._EB
 U_LS_OTHER = 0._EB
 V_LS_OTHER = 0._EB
 Z_LS_OTHER = 0._EB
 DO JJO=EWC%JJO_MIN,EWC%JJO_MAX
    DO IIO=EWC%IIO_MIN,EWC%IIO_MAX
+      IIO_2 = IIO
+      JJO_2 = JJO
+      SELECT CASE(IOR)
+         CASE( 1) ; IIO_2 = IIO - 1
+         CASE(-1) ; IIO_2 = IIO + 1
+         CASE( 2) ; JJO_2 = JJO - 1
+         CASE(-2) ; JJO_2 = JJO + 1
+      END SELECT
       IF (PREDICTOR) THEN
-         PHI_LS_OTHER = PHI_LS_OTHER + OMESH(NOM)%PHI_LS(IIO,JJO)
+         PHI_LS_OTHER   = PHI_LS_OTHER   + OMESH(NOM)%PHI_LS(IIO  ,JJO  )
+         PHI_LS_OTHER_2 = PHI_LS_OTHER_2 + OMESH(NOM)%PHI_LS(IIO_2,JJO_2)
       ELSE
-         PHI_LS_OTHER = PHI_LS_OTHER + OMESH(NOM)%PHI1_LS(IIO,JJO)
+         PHI_LS_OTHER   = PHI_LS_OTHER   + OMESH(NOM)%PHI1_LS(IIO  ,JJO  )
+         PHI_LS_OTHER_2 = PHI_LS_OTHER_2 + OMESH(NOM)%PHI1_LS(IIO_2,JJO_2)
       ENDIF
       U_LS_OTHER = U_LS_OTHER + OMESH(NOM)%U_LS(IIO,JJO)
       V_LS_OTHER = V_LS_OTHER + OMESH(NOM)%V_LS(IIO,JJO)
@@ -728,6 +741,10 @@ SUBROUTINE FILL_BOUNDARY_VALUES
 SELECT CASE(IOR)
    CASE(-2:2)
       PHI_LS_P(II,JJ) = PHI_LS_OTHER/REAL(N_INT_CELLS,EB)
+      IF (IOR==-2) PHI_LS_P(II,JJ-1) = PHI_LS_OTHER_2/REAL(N_INT_CELLS,EB)
+      IF (IOR== 2) PHI_LS_P(II,JJ+1) = PHI_LS_OTHER_2/REAL(N_INT_CELLS,EB)
+      IF (IOR==-1) PHI_LS_P(II-1,JJ) = PHI_LS_OTHER_2/REAL(N_INT_CELLS,EB)
+      IF (IOR== 1) PHI_LS_P(II+1,JJ) = PHI_LS_OTHER_2/REAL(N_INT_CELLS,EB)
       U_LS(II,JJ)     = U_LS_OTHER/REAL(N_INT_CELLS,EB)
       V_LS(II,JJ)     = V_LS_OTHER/REAL(N_INT_CELLS,EB)
       Z_LS(II,JJ)     = Z_LS_OTHER/REAL(N_INT_CELLS,EB)
@@ -997,7 +1014,7 @@ END SUBROUTINE LEVEL_SET_SPREAD_RATE
 
 SUBROUTINE LEVEL_SET_ADVECT_FLUX
 
-INTEGER :: I,IM1,IP1,IP2,J,JM1,JP1,JP2
+INTEGER :: I,J
 REAL(EB), DIMENSION(:) :: Z(4)
 REAL(EB), POINTER, DIMENSION(:,:) :: FLUX_LS_P,F_X,F_Y
 REAL(EB) :: DPHIDX,DPHIDY,SR_X_AVG,SR_Y_AVG
@@ -1015,28 +1032,22 @@ SUBROUTINE LEVEL_SET_ADVECT_FLUX
 
 DO J=1,JBAR
    DO I=0,IBAR
-      IM1 = I-1 ; IF (IM1<0) IM1 = I
-      IP1 = I+1
-      IP2 = I+2 ; IF (IP2>IBP1) IP2 = IP1
-      Z(1) = PHI_LS_P(IM1,J)
-      Z(2) = PHI_LS_P(I,J)
-      Z(3) = PHI_LS_P(IP1,J)
-      Z(4) = PHI_LS_P(IP2,J)
-      SR_X_AVG = 0.5_EB*(SR_X_LS(MIN(IP1,IBAR),J)+SR_X_LS(MAX(1,I),J))
+      Z(1) = PHI_LS_P(I-1,J)
+      Z(2) = PHI_LS_P(I  ,J)
+      Z(3) = PHI_LS_P(I+1,J)
+      Z(4) = PHI_LS_P(I+2,J)
+      SR_X_AVG = 0.5_EB*(SR_X_LS(MIN(I+1,IBAR),J)+SR_X_LS(MAX(1,I),J))
       F_X(I,J) = SCALAR_FACE_VALUE_LS(SR_X_AVG,Z,LIMITER_LS)
    ENDDO
 ENDDO
 
 DO J=0,JBAR
    DO I=1,IBAR
-      JM1 = J-1 ; IF (JM1<0) JM1 = J
-      JP1 = J+1
-      JP2 = J+2 ; IF (JP2>JBP1) JP2 = JP1
-      Z(1) = PHI_LS_P(I,JM1)
-      Z(2) = PHI_LS_P(I,J)
-      Z(3) = PHI_LS_P(I,JP1)
-      Z(4) = PHI_LS_P(I,JP2)
-      SR_Y_AVG = 0.5_EB*(SR_Y_LS(I,MIN(JP1,JBAR))+SR_Y_LS(I,MAX(1,J)))
+      Z(1) = PHI_LS_P(I,J-1)
+      Z(2) = PHI_LS_P(I,J  )
+      Z(3) = PHI_LS_P(I,J+1)
+      Z(4) = PHI_LS_P(I,J+2)
+      SR_Y_AVG = 0.5_EB*(SR_Y_LS(I,MIN(J+1,JBAR))+SR_Y_LS(I,MAX(1,J)))
       F_Y(I,J) = SCALAR_FACE_VALUE_LS(SR_Y_AVG,Z,LIMITER_LS)
    ENDDO
 ENDDO