Merge pull request #14150 from mcgratta/precon

FDS Source: Restore version of TUNNEL_PRECONDITION pre-Dec 14, 2023

Author: mcgratta

Source/cons.f90

@@ -534,9 +534,6 @@ MODULE GLOBAL_CONSTANTS
 REAL(EB), ALLOCATABLE, DIMENSION(:) :: TP_DD                     !< Diagonal of tri-diagonal matrix for tunnel pressure solver
 REAL(EB), ALLOCATABLE, DIMENSION(:) :: TP_RDXN                   !< Reciprocal of the distance between tunnel precon points
 REAL(EB), ALLOCATABLE, TARGET, DIMENSION(:) :: H_BAR             !< Pressure solution of 1-D tunnel pressure solver, predictor step
-REAL(EB), ALLOCATABLE, TARGET, DIMENSION(:) :: H_BAR_S           !< Pressure solution of 1-D tunnel pressure solver, corrector step
-REAL(EB), ALLOCATABLE, TARGET, DIMENSION(:) :: DUDT_BAR          !< Average du/dt for 1-D tunnel pressure solver, predictor step
-REAL(EB), ALLOCATABLE, TARGET, DIMENSION(:) :: DUDT_BAR_S        !< Average du/dt for 1-D tunnel pressure solver, corrector step
 INTEGER, ALLOCATABLE, DIMENSION(:) :: COUNTS_TP                  !< Counter for MPI calls used for 1-D tunnel pressure solver
 INTEGER, ALLOCATABLE, DIMENSION(:) :: DISPLS_TP                  !< Displacements for MPI calls used for 1-D tunnel pressure solver
 INTEGER, ALLOCATABLE, DIMENSION(:) :: I_OFFSET                   !< Spatial index of tunnel

Source/dump.f90

@@ -7307,8 +7307,6 @@ REAL(EB) RECURSIVE FUNCTION GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,IND,IND2,Y_INDEX,Z
    CASE( 9)  ! PRESSURE
       GAS_PHASE_OUTPUT_RES = PBAR(KK,PRESSURE_ZONE(II,JJ,KK)) + &
                              RHO(II,JJ,KK)*(0.5_EB*(H(II,JJ,KK)+HS(II,JJ,KK))-KRES(II,JJ,KK)) - P_0(KK)
-      IF (TUNNEL_PRECONDITIONER) GAS_PHASE_OUTPUT_RES = GAS_PHASE_OUTPUT_RES + &
-                                                        RHO(II,JJ,KK)*0.5_EB*(H_BAR(I_OFFSET(NM)+II)+H_BAR_S(I_OFFSET(NM)+II))
    CASE(10)  ! VELOCITY
       SELECT CASE(ABS(VELO_INDEX))
          CASE DEFAULT
@@ -7327,8 +7325,6 @@ REAL(EB) RECURSIVE FUNCTION GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,IND,IND2,Y_INDEX,Z
       GAS_PHASE_OUTPUT_RES = Q(II,JJ,KK)*0.001_EB
    CASE(12)  ! H
       GAS_PHASE_OUTPUT_RES = 0.5_EB*(HS(II,JJ,KK)+H(II,JJ,KK))
-      IF (TUNNEL_PRECONDITIONER) GAS_PHASE_OUTPUT_RES = GAS_PHASE_OUTPUT_RES + &
-                                                        0.5_EB*(H_BAR(I_OFFSET(NM)+II)+H_BAR_S(I_OFFSET(NM)+II))
    CASE(13)  ! MIXTURE FRACTION
       ! requires FUEL + AIR --> PROD (SIMPLE_CHEMISTRY, N_SIMPLE_CHEMISTRY_REACTIONS=1)
       ! f = Z_FUEL + Z_PROD/(1+S), where S is the mass stoichiometric coefficient for AIR

Source/main.f90

@@ -1419,14 +1419,11 @@ SUBROUTINE PRESSURE_ITERATION_SCHEME
       CALL PRESSURE_SOLVER_COMPUTE_RHS(T,DT,NM)
    ENDDO
 
-   ! Special case for tunnels -- filter out the lengthwise pressure gradient
-
-   IF (TUNNEL_PRECONDITIONER) CALL TUNNEL_POISSON_SOLVER
-
    ! Solve the Poission equation using either FFT or ULMAT, GLMAT, or UGLMAT
 
    SELECT CASE(PRES_FLAG)
       CASE (FFT_FLAG)
+         IF (TUNNEL_PRECONDITIONER) CALL TUNNEL_POISSON_SOLVER
          DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX
             CALL PRESSURE_SOLVER_FFT(NM)
          ENDDO
@@ -1435,6 +1432,7 @@ SUBROUTINE PRESSURE_ITERATION_SCHEME
          CALL MESH_EXCHANGE(5)
          CALL COPY_H_OMESH_TO_MESH
       CASE (ULMAT_FLAG)
+         IF (TUNNEL_PRECONDITIONER) CALL TUNNEL_POISSON_SOLVER
          DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX
             CALL ULMAT_SOLVER(NM,T,DT)
          ENDDO

Source/mesh.f90

@@ -86,6 +86,7 @@ MODULE MESH_VARIABLES
    REAL(EB), ALLOCATABLE, DIMENSION(:) :: SAVE1,SAVE2,WORK
    REAL(EB), ALLOCATABLE, DIMENSION(:,:,:) :: PRHS !< Right hand side of Poisson (pressure) equation
    REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: BXS,BXF,BYS,BYF,BZS,BZF, BXST,BXFT,BYST,BYFT,BZST,BZFT
+   REAL(EB) :: BXS_BAR,BXF_BAR
    INTEGER :: LSAVE,LWORK,LBC,MBC,NBC,LBC2,MBC2,NBC2,ITRN,JTRN,KTRN,IPS
 
    REAL(EB), ALLOCATABLE, DIMENSION(:) :: P_0         !< Ambient pressure profile, \f$ \overline{p}_0(z) \f$ (Pa)
@@ -359,6 +360,7 @@ MODULE MESH_POINTERS
 REAL(EB), POINTER, DIMENSION(:) :: SAVE1,SAVE2,WORK
 REAL(EB), POINTER, DIMENSION(:,:,:) :: PRHS
 REAL(EB), POINTER, DIMENSION(:,:) :: BXS,BXF,BYS,BYF,BZS,BZF, BXST,BXFT,BYST,BYFT,BZST,BZFT
+REAL(EB), POINTER :: BXS_BAR,BXF_BAR
 INTEGER, POINTER :: LSAVE,LWORK,LBC,MBC,NBC,LBC2,MBC2,NBC2,ITRN,JTRN,KTRN,IPS
 REAL(EB), POINTER, DIMENSION(:) :: P_0,RHO_0,TMP_0,D_PBAR_DT,D_PBAR_DT_S,U_LEAK,U_WIND,V_WIND,W_WIND
 REAL(EB), POINTER, DIMENSION(:,:) :: PBAR,PBAR_S,R_PBAR
@@ -574,6 +576,8 @@ SUBROUTINE POINT_TO_MESH(NM)
 BYFT=>M%BYFT
 BZST=>M%BZST
 BZFT=>M%BZFT
+BXS_BAR=>M%BXS_BAR
+BXF_BAR=>M%BXF_BAR
 LBC=>M%LBC
 MBC=>M%MBC
 NBC=>M%NBC