Skip to content

Commit 0825146

Browse files
mcgrattamcgratta
committed
FDS Source: Make ULMAT work with TUNNEL_PRECON
1 parent 1232dc2 commit 0825146

File tree

2 files changed

+72
-56
lines changed

2 files changed

+72
-56
lines changed

‎Source/pres.f90‎

Lines changed: 44 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -523,6 +523,16 @@ SUBROUTINE TUNNEL_POISSON_SOLVER
523523

524524
M => MESHES(NM)
525525

526+
IF(PRES_FLAG==ULMAT_FLAG) THEN
527+
DO K=1,M%KBAR
528+
DO J=1,M%JBAR
529+
DO I=1,M%IBAR
530+
IF(M%CELL(M%CELL_INDEX(I,J,K))%SOLID) M%PRHS(I,J,K) = 0._EB
531+
ENDDO
532+
ENDDO
533+
ENDDO
534+
ENDIF
535+
526536
TP_RDXN(I_OFFSET(NM):I_OFFSET(NM)+M%IBAR) = M%RDXN(0:M%IBAR)
527537
IF (NM>1) TP_RDXN(I_OFFSET(NM)) = 2._EB/(MESHES(NM-1)%DX(MESHES(NM-1)%IBAR)+M%DX(1))
528538
IF (NM<NMESHES) TP_RDXN(I_OFFSET(NM)+M%IBAR) = 2._EB/(MESHES(NM+1)%DX(1) +M%DX(M%IBAR))
@@ -1419,7 +1429,7 @@ SUBROUTINE ULMAT_SOLVE_ZONE(NM,IPZ)
14191429

14201430
! Local Variables:
14211431
INTEGER :: NRHS,MAXFCT,MNUM,ERROR,I,J,K,ICC,JCC,IIG,JJG,KKG,IOR,IW,IROW,NCELL,ICFACE,IFACE,JFACE,ICVL,ILH,JLH,KLH,IRC
1422-
REAL(EB):: SUM_FH(1:2),MEAN_FH,SUM_XH(1:2),MEAN_XH,DIV_FN_VOL,DIV_FN,IDX,AF,VAL,BCV
1432+
REAL(EB):: SUM_FH(1:2),MEAN_FH,SUM_XH(1:2),MEAN_XH,DIV_FN_VOL,DIV_FN,IDX,AF,VAL,BCV,DHDN
14231433
TYPE(ZONE_MESH_TYPE), POINTER :: ZM
14241434
TYPE (WALL_TYPE), POINTER :: WC
14251435
TYPE (EXTERNAL_WALL_TYPE), POINTER :: EWC
@@ -1503,6 +1513,39 @@ SUBROUTINE ULMAT_SOLVE_ZONE(NM,IPZ)
15031513
ENDIF IF_CFACE_DIRICHLET
15041514
ENDDO CFACE_LOOP
15051515

1516+
! Tunnel Preconditioner, internal wall cells normal to x take Neumann BC = - dH_BAR/dx
1517+
IF (TUNNEL_PRECONDITIONER) THEN
1518+
WALL_CELL_LOOP_0 : DO IW=N_EXTERNAL_WALL_CELLS+1,N_EXTERNAL_WALL_CELLS+N_INTERNAL_WALL_CELLS
1519+
WC => WALL(IW)
1520+
IF (WC%BOUNDARY_TYPE/=SOLID_BOUNDARY .OR. WC%CUT_FACE_INDEX>0) CYCLE WALL_CELL_LOOP_0
1521+
BC => BOUNDARY_COORD(WC%BC_INDEX);
1522+
IOR = BC%IOR; IF(ABS(IOR)/=1) CYCLE WALL_CELL_LOOP_0 ! Only in x.
1523+
! Gasphase cell indexes:
1524+
IIG = BC%IIG; JJG = BC%JJG; KKG = BC%KKG
1525+
IF(ZONE_MESH(PRESSURE_ZONE(IIG,JJG,KKG))%CONNECTED_ZONE_PARENT/=IPZ) CYCLE WALL_CELL_LOOP_0
1526+
IROW = MUNKH(IIG,JJG,KKG)
1527+
IF(CC_IBM) THEN
1528+
IF (IROW <= 0) THEN
1529+
ICC = CCVAR(IIG,JJG,KKG,CC_IDCC); IF(ICC<1) CYCLE WALL_CELL_LOOP_0
1530+
! Note: this only works with single pressure unknown per cartesian cell.
1531+
IROW = CUT_CELL(ICC)%UNKH(1)
1532+
ENDIF
1533+
ENDIF
1534+
SELECT CASE (IOR)
1535+
CASE(-1) ! -IAXIS oriented, high face of IIG cell.
1536+
DHDN = (H_BAR(I_OFFSET(NM)+IIG+1)-H_BAR(I_OFFSET(NM)+IIG))*TP_RDXN(I_OFFSET(NM)+IIG)
1537+
AF = ((1._EB-CYL_FCT)*DY(JJG) + CYL_FCT*R(IIG )) * DZ(KKG)
1538+
VAL = -DHDN*AF
1539+
CASE( 1) ! +IAXIS oriented, low face of IIG cell.
1540+
DHDN = (H_BAR(I_OFFSET(NM)+IIG)-H_BAR(I_OFFSET(NM)+IIG-1))*TP_RDXN(I_OFFSET(NM)+IIG-1)
1541+
AF = ((1._EB-CYL_FCT)*DY(JJG) + CYL_FCT*R(IIG-1)) * DZ(KKG)
1542+
VAL = DHDN*AF
1543+
END SELECT
1544+
! Add to F_H:
1545+
ZM%F_H(IROW) = ZM%F_H(IROW) - VAL
1546+
ENDDO WALL_CELL_LOOP_0
1547+
ENDIF
1548+
15061549
! Finally add External Wall cell BCs:
15071550
WALL_CELL_LOOP_1 : DO IW=1,N_EXTERNAL_WALL_CELLS
15081551
WC => WALL(IW)

‎Source/velo.f90‎

Lines changed: 28 additions & 55 deletions
Original file line numberDiff line numberDiff line change
@@ -1643,7 +1643,7 @@ SUBROUTINE VELOCITY_PREDICTOR(T,DT,DT_NEW,NM)
16431643
T_NOW=CURRENT_TIME()
16441644
ENDIF
16451645
SELECT CASE(PRES_FLAG)
1646-
CASE(GLMAT_FLAG,UGLMAT_FLAG,ULMAT_FLAG); CALL WALL_VELOCITY_NO_GRADH(NM,DT,.FALSE.)
1646+
CASE(GLMAT_FLAG,UGLMAT_FLAG,ULMAT_FLAG); CALL WALL_VELOCITY_NO_GRADH(DT,.FALSE.)
16471647
END SELECT
16481648

16491649
ENDIF FREEZE_VELOCITY_IF
@@ -1721,7 +1721,7 @@ SUBROUTINE VELOCITY_CORRECTOR(T,DT,NM)
17211721
ENDIF
17221722
SELECT CASE(PRES_FLAG)
17231723
CASE(GLMAT_FLAG,UGLMAT_FLAG,ULMAT_FLAG)
1724-
CALL WALL_VELOCITY_NO_GRADH(NM,DT,.TRUE.) ! Store U velocities on OBST surfaces.
1724+
CALL WALL_VELOCITY_NO_GRADH(DT,.TRUE.) ! Store U velocities on OBST surfaces.
17251725
END SELECT
17261726

17271727
!$OMP PARALLEL PRIVATE(I,J,K)
@@ -1765,7 +1765,7 @@ SUBROUTINE VELOCITY_CORRECTOR(T,DT,NM)
17651765
ENDIF
17661766
SELECT CASE(PRES_FLAG)
17671767
CASE(GLMAT_FLAG,UGLMAT_FLAG,ULMAT_FLAG)
1768-
CALL WALL_VELOCITY_NO_GRADH(NM,DT,.FALSE.)
1768+
CALL WALL_VELOCITY_NO_GRADH(DT,.FALSE.)
17691769
END SELECT
17701770

17711771
ENDIF FREEZE_VELOCITY_IF
@@ -3235,13 +3235,12 @@ END SUBROUTINE BAROCLINIC_CORRECTION
32353235
!> \details Ensure that the correct normal derivative of H is used on the projection. It is only used when the Poisson equation
32363236
!> for the pressure is solved .NOT. PRES_ON_WHOLE_DOMAIN (i.e. using the GLMAT solver).
32373237

3238-
SUBROUTINE WALL_VELOCITY_NO_GRADH(NM,DT,STORE_UN)
3238+
SUBROUTINE WALL_VELOCITY_NO_GRADH(DT,STORE_UN)
32393239

32403240
REAL(EB), INTENT(IN) :: DT
32413241
LOGICAL, INTENT(IN) :: STORE_UN
3242-
INTEGER, INTENT(IN) :: NM
32433242
INTEGER :: II,JJ,KK,IIG,JJG,KKG,IOR,IW,N_INTERNAL_WALL_CELLS_AUX,IC,ICG
3244-
REAL(EB) :: DHDN, VEL_N
3243+
REAL(EB) :: VEL_N
32453244
TYPE (WALL_TYPE), POINTER :: WC
32463245
REAL(EB), SAVE, ALLOCATABLE, DIMENSION(:) :: UN_WALLS
32473246
TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC
@@ -3263,10 +3262,7 @@ SUBROUTINE WALL_VELOCITY_NO_GRADH(NM,DT,STORE_UN)
32633262

32643263
WC => WALL(IW)
32653264
BC => BOUNDARY_COORD(WC%BC_INDEX)
3266-
IIG = BC%IIG
3267-
JJG = BC%JJG
3268-
KKG = BC%KKG
3269-
IOR = BC%IOR
3265+
IIG= BC%IIG; JJG= BC%JJG; KKG= BC%KKG; IOR= BC%IOR
32703266

32713267
SELECT CASE(IOR)
32723268
CASE( IAXIS)
@@ -3300,38 +3296,28 @@ SUBROUTINE WALL_VELOCITY_NO_GRADH(NM,DT,STORE_UN)
33003296
WC%BOUNDARY_TYPE==MIRROR_BOUNDARY) ) CYCLE
33013297

33023298
BC => BOUNDARY_COORD(WC%BC_INDEX)
3303-
II = BC%II
3304-
JJ = BC%JJ
3305-
KK = BC%KK
3306-
IIG = BC%IIG
3307-
JJG = BC%JJG
3308-
KKG = BC%KKG
3299+
II = BC%II; JJ = BC%JJ; KK = BC%KK; IIG = BC%IIG; JJG = BC%JJG; KKG = BC%KKG
33093300
IC = CELL_INDEX(II,JJ,KK)
33103301
ICG = CELL_INDEX(IIG,JJG,KKG)
33113302

33123303
IF (WC%BOUNDARY_TYPE==NULL_BOUNDARY .AND. .NOT.CELL(ICG)%SOLID .AND. .NOT.CELL(IC)%SOLID) CYCLE
33133304

33143305
IOR = BC%IOR
3315-
DHDN=0._EB ! Set the normal derivative of H to zero for solids.
33163306

3317-
SELECT CASE(IOR)
3318-
CASE( IAXIS)
3319-
IF (TUNNEL_PRECONDITIONER .AND. IW>N_EXTERNAL_WALL_CELLS) &
3320-
DHDN = (H_BAR(I_OFFSET(NM)+IIG)-H_BAR(I_OFFSET(NM)+IIG-1))*TP_RDXN(I_OFFSET(NM)+IIG-1)
3321-
US(IIG-1,JJG ,KKG ) = (U(IIG-1,JJG ,KKG ) - DT*( FVX(IIG-1,JJG ,KKG ) + DHDN ))
3322-
CASE(-IAXIS)
3323-
IF (TUNNEL_PRECONDITIONER .AND. IW>N_EXTERNAL_WALL_CELLS) &
3324-
DHDN = (H_BAR(I_OFFSET(NM)+IIG+1)-H_BAR(I_OFFSET(NM)+IIG))*TP_RDXN(I_OFFSET(NM)+IIG)
3325-
US(IIG ,JJG ,KKG ) = (U(IIG ,JJG ,KKG ) - DT*( FVX(IIG ,JJG ,KKG ) + DHDN ))
3326-
CASE( JAXIS)
3327-
VS(IIG ,JJG-1,KKG ) = (V(IIG ,JJG-1,KKG ) - DT*( FVY(IIG ,JJG-1,KKG ) + DHDN ))
3328-
CASE(-JAXIS)
3329-
VS(IIG ,JJG ,KKG ) = (V(IIG ,JJG ,KKG ) - DT*( FVY(IIG ,JJG ,KKG ) + DHDN ))
3330-
CASE( KAXIS)
3331-
WS(IIG ,JJG ,KKG-1) = (W(IIG ,JJG ,KKG-1) - DT*( FVZ(IIG ,JJG ,KKG-1) + DHDN ))
3332-
CASE(-KAXIS)
3333-
WS(IIG ,JJG ,KKG ) = (W(IIG ,JJG ,KKG ) - DT*( FVZ(IIG ,JJG ,KKG ) + DHDN ))
3334-
END SELECT
3307+
SELECT CASE(IOR)
3308+
CASE( IAXIS)
3309+
US(IIG-1,JJG ,KKG ) = (U(IIG-1,JJG ,KKG ) - DT*( FVX(IIG-1,JJG ,KKG ) ))
3310+
CASE(-IAXIS)
3311+
US(IIG ,JJG ,KKG ) = (U(IIG ,JJG ,KKG ) - DT*( FVX(IIG ,JJG ,KKG ) ))
3312+
CASE( JAXIS)
3313+
VS(IIG ,JJG-1,KKG ) = (V(IIG ,JJG-1,KKG ) - DT*( FVY(IIG ,JJG-1,KKG ) ))
3314+
CASE(-JAXIS)
3315+
VS(IIG ,JJG ,KKG ) = (V(IIG ,JJG ,KKG ) - DT*( FVY(IIG ,JJG ,KKG ) ))
3316+
CASE( KAXIS)
3317+
WS(IIG ,JJG ,KKG-1) = (W(IIG ,JJG ,KKG-1) - DT*( FVZ(IIG ,JJG ,KKG-1) ))
3318+
CASE(-KAXIS)
3319+
WS(IIG ,JJG ,KKG ) = (W(IIG ,JJG ,KKG ) - DT*( FVZ(IIG ,JJG ,KKG ) ))
3320+
END SELECT
33353321

33363322
ENDDO WALL_CELL_LOOP_1
33373323

@@ -3358,34 +3344,21 @@ SUBROUTINE WALL_VELOCITY_NO_GRADH(NM,DT,STORE_UN)
33583344
IF (WC%BOUNDARY_TYPE==NULL_BOUNDARY .AND. .NOT.CELL(ICG)%SOLID .AND. .NOT.CELL(IC)%SOLID) CYCLE
33593345

33603346
IOR = BC%IOR
3361-
DHDN=0._EB ! Set the normal derivative of H to zero for solids.
3362-
33633347
VEL_N = UN_WALLS(IW)
33643348

33653349
SELECT CASE(IOR)
3366-
CASE( IAXIS) ! | - Problem with this is it was modified in VELOCITY_CORRECTOR,
3367-
! V => Store the untouched U normal on internal WALLs.
3368-
IF (TUNNEL_PRECONDITIONER.AND. IW>N_EXTERNAL_WALL_CELLS) &
3369-
DHDN = (H_BAR(I_OFFSET(NM)+IIG)-H_BAR(I_OFFSET(NM)+IIG-1))*TP_RDXN(I_OFFSET(NM)+IIG-1)
3370-
U(IIG-1,JJG ,KKG ) = 0.5_EB*( VEL_N + US(IIG-1,JJG ,KKG ) - &
3371-
DT*( FVX(IIG-1,JJG ,KKG ) + DHDN ))
3350+
CASE( IAXIS)
3351+
U(IIG-1,JJG ,KKG ) = 0.5_EB*(VEL_N + US(IIG-1,JJG ,KKG ) - DT*( FVX(IIG-1,JJG ,KKG ) ))
33723352
CASE(-IAXIS)
3373-
IF (TUNNEL_PRECONDITIONER.AND. IW>N_EXTERNAL_WALL_CELLS) &
3374-
DHDN = (H_BAR(I_OFFSET(NM)+IIG+1)-H_BAR(I_OFFSET(NM)+IIG))*TP_RDXN(I_OFFSET(NM)+IIG)
3375-
U(IIG ,JJG ,KKG ) = 0.5_EB*( VEL_N + US(IIG ,JJG ,KKG ) - &
3376-
DT*( FVX(IIG ,JJG ,KKG ) + DHDN ))
3353+
U(IIG ,JJG ,KKG ) = 0.5_EB*(VEL_N + US(IIG ,JJG ,KKG ) - DT*( FVX(IIG ,JJG ,KKG ) ))
33773354
CASE( JAXIS)
3378-
V(IIG ,JJG-1,KKG ) = 0.5_EB*( VEL_N + VS(IIG ,JJG-1,KKG ) - &
3379-
DT*( FVY(IIG ,JJG-1,KKG ) + DHDN ))
3355+
V(IIG ,JJG-1,KKG ) = 0.5_EB*(VEL_N + VS(IIG ,JJG-1,KKG ) - DT*( FVY(IIG ,JJG-1,KKG ) ))
33803356
CASE(-JAXIS)
3381-
V(IIG ,JJG ,KKG ) = 0.5_EB*( VEL_N + VS(IIG ,JJG ,KKG ) - &
3382-
DT*( FVY(IIG ,JJG ,KKG ) + DHDN ))
3357+
V(IIG ,JJG ,KKG ) = 0.5_EB*(VEL_N + VS(IIG ,JJG ,KKG ) - DT*( FVY(IIG ,JJG ,KKG ) ))
33833358
CASE( KAXIS)
3384-
W(IIG ,JJG ,KKG-1) = 0.5_EB*( VEL_N + WS(IIG ,JJG ,KKG-1) - &
3385-
DT*( FVZ(IIG ,JJG ,KKG-1) + DHDN ))
3359+
W(IIG ,JJG ,KKG-1) = 0.5_EB*(VEL_N + WS(IIG ,JJG ,KKG-1) - DT*( FVZ(IIG ,JJG ,KKG-1) ))
33863360
CASE(-KAXIS)
3387-
W(IIG ,JJG ,KKG ) = 0.5_EB*( VEL_N + WS(IIG ,JJG ,KKG ) - &
3388-
DT*( FVZ(IIG ,JJG ,KKG ) + DHDN ))
3361+
W(IIG ,JJG ,KKG ) = 0.5_EB*(VEL_N + WS(IIG ,JJG ,KKG ) - DT*( FVZ(IIG ,JJG ,KKG ) ))
33893362
END SELECT
33903363

33913364
ENDDO WALL_CELL_LOOP_2

0 commit comments

Comments
 (0)