** MB-FGE Test: Swirling flow through expanding pipe.
  **************************************************************
  DISPLAY
   This case concerns axisymmetric incompressible laminar/
   turbulent flow in the expanding pipe. The flow rotates,
   with linear swirling velocity profile at the inlet. The
   flow exhibits recirculation not only behind the step at
   the pipe wall, but in the region on the pipe axis.
 
   Calculations  of  the  turbulent  flow (LTURB=T) could be
   made with low-Re extensions of whether standard K-E model
   (LOWKE=T) or Chen-Kim K-E model (LCHEK=T).
 
   Computational domain  is built  from 2  blocks with  grid
   sizes NX1,NY1,NZ1 and NX2,NY2,NZ2 respectively. The  grid
   for the after step region can be IFCZ times finer in  the
   Z-direction as the grid at the top.
   ---------------------------------------------------------
  ENDDIS
L(PAUSE
  **************************************************************
INTEGER(IFCZ);  IFCZ= 2
BOOLEAN(LTURB,LOWKE,LCHEK); LTURB= F;  LOWKE= F;  LCHEK= F
  **************************************************************
  PHOTON USE
   p ; ; ; ; ;
 
   set vec av off
   msg Velocity vectors:
   vec i 1 sh
   msg Press any key to continue ...
   pause
   cl
   msg Pressure contours:
   con p1 i 1 fi;0.01
   pause
   msg Press any key to continue ...
   cl
   msg Swirl velocity contours:
   con uc1 i 1 fi;0.01
   pause
   msg Press any key to continue ...
   cl
   msg W-component contours:
   con w1 i 1 fi;0.01
   msg Press any key to continue ...
   pause
   cl
   msg V-component contours:
   con v1 i 1 fi;0.01
  ENDUSE
  **************************************************************
    GROUP 1. Run title and other preliminaries
IF(LTURB) THEN
+ IF(LOWRE) THEN
+  TEXT(MB-FGE: Swirl-flow (K-E low-Re).
+ ELSE
+  TEXT(MB-FGE: Swirl-flow (Chen-Kim K-E low-Re).
+ ENDIF
ELSE
+ TEXT(MB-FGE: Swirl-flow (Re=100).
ENDIF
TITLE
INTEGER(NX1,NY1,NZ1,NZST,NX2,NY2,NZ2)
REAL(REYNO,WIN,UIN,HSTEP,LSTPP,HCHAN,LCHAN,DXX,DTHYD,UCUR)
REAL(TKEIN,EPSIN)
  ** Problem definition:
IF(LTURB) THEN
+ REYNO= 1.0E5
ELSE
+ REYNO= 100.0
ENDIF
HSTEP= 1.0;  LSTPP= 2.0;  HCHAN= 2.0;  LCHAN= 20.0
WIN  = 1.0;  DXX  = 0.1;  UIN  = 3.0;  ENUL = WIN*HSTEP/REYNO
   ** First sub-domain.
NX1  = 1;    NY1  = 12;   NZ1  = 26;   NZST = 6;
   ** Second sub-domain.
NX2  = 1;    NY2  = 12;   NZ2  = (NZ1-NZST)*IFCZ
    GROUP 6. Body-fitted coordinates or grid distortion
BFC = T
  ** Define grid points and lines for the first sub-domain:
GSET(P,P1,0.0,        0.0,  0.0);GSET(P,P2,0.0,        0.0,LSTPP)
GSET(P,P3,0.0,        0.0,LCHAN);GSET(P,P4,0.0,HCHAN-HSTEP,LCHAN)
GSET(P,P5,0.0,HCHAN-HSTEP,LSTPP);GSET(P,P6,0.0,HCHAN-HSTEP,  0.0)
GSET(L,L12,P1,P2,NZST,-1.3);     GSET(L,L23,P2,P3,NZ1-NZST,1.0)
GSET(L,L34,P3,P4, NY1, 1.0);     GSET(L,L45,P4,P5,NZ1-NZST,1.0)
GSET(L,L56,P5,P6,NZST, 1.3);     GSET(L,L61,P6,P1,     NY1,1.0)
  ** Define grid points and lines for the second sub-domain:
GSET(P,P7,0.0,HCHAN-HSTEP,LSTPP);GSET(P,P8, 0.0,HCHAN-HSTEP,LCHAN)
GSET(P,P9,0.0,      HCHAN,LCHAN);GSET(P,P10,0.0,      HCHAN,LSTPP)
GSET(L,L78, P7,P8, NZ2, 1.0);   GSET(L,L89, P8, P9,NY2, 1.0)
GSET(L,L910,P9,P10,NZ2, 1.0);   GSET(L,L107,P10,P7,NY2, 1.0)
  ** Create grid for the first sub-domain.
GSET(D,NX1,NY1,NZ1,DXX,HCHAN-HSTEP,LCHAN)
GSET(F,F1,P1,P2,P3,-,P4,P5,P6,-); GSET(M,F1,+K+J,1,1,1)
GSET(C,I:NX1+1:,F,I1,1,NY1,1,NZ1,RZ,-DXX,0.0,0.0,INC,1.0)
DUMPC(MBGR1)
  ** Create grid for the second sub-domain.
GSET(D,NX2,NY2,NZ2,DXX,HSTEP,LCHAN-LSTPP)
GSET(F,F2,P7,-,P8,-,P9,-,P10,-); GSET(M,F2,+K+J,1,1,1)
GSET(C,I:NX2+1:,F,I1,1,NY2,1,NZ2,RZ,-DXX,0.0,0.0,INC,1.0)
DUMPC(MBGR2)
  ** Assemble blocks:
NUMBLK = 2; READCO(MBGR+L)
GVIEW(X); VIEW
  ** Set links:
    ** You can use MBLINK command instead of READCO(...+L).
    MBLINK(1,NORTH,2,SOUTH)
    GROUP 7. Variables stored, solved & named
STORE(VPOR); NAME(C3)= UC1; SOLVE(P1,V1,W1,UC1)
IF(LTURB) THEN
+ IF(LCHEK) THEN
+  TURMOD(KECHEN-LOWRE)
+ ELSE
+  TURMOD(KEMODL-LOWRE)
+ ENDIF
+ STORE(ENUT,LEN1)
ENDIF
L($F150)
    GROUP 11. Initialization of variable or porosity fields
INIADD= F; FIINIT(UC1)= 1.E-3
IF(LTURB) THEN
+ TKEIN= (0.05*WIN)**2; EPSIN     = TKEIN**1.5*0.1643/0.09*HSTEP
+ FIINIT(P1)= 1.3E-4;   FIINIT(KE)= TKEIN;  FIINIT(EP)= EPSIN
ENDIF
    GROUP 13. Boundary conditions and special sources
    ** Inlet.
DO II = 1,NY1
+ UCUR = (2*II-1)*UIN/NY1/2
+ MPATCH(1,INL:II:,LOW,1,NX1,II,II,1,1,1,LSTEP)
+ COVAL(INL:II:,P1,FIXFLU,WIN*RHO1);COVAL(INL:II:,UC1,ONLYMS,UCUR)
+ COVAL(INL:II:,WC1,ONLYMS,WIN);COVAL(INL:II:,VC1,ONLYMS,0.0)
+ IF(LTURB) THEN
+  COVAL(INL:II:,KE,ONLYMS,TKEIN);COVAL(INL:II:,EP,ONLYMS,EPSIN)
+ ENDIF
ENDDO
    ** Outlet.
MPATCH(1,OUT1,HIGH,1,NX1,1,NY1,NZ1,NZ1,1,LSTEP)
COVAL(OUT1,P1,1.E5,0.0)
MPATCH(2,OUT2,HIGH,1,NX2,1,NY2,NZ2,NZ2,1,LSTEP)
COVAL(OUT2,P1,1.E5,0.0)
    ** Walls.
MPATCH(1,WN1,NWALL,1,NX1,NY1,NY1,1,NZST,1,LSTEP)
MPATCH(2,WN2,NWALL,1,NX2,NY2,NY2,1, NZ2,1,LSTEP)
MPATCH(2,WL, LWALL,1,NX2,  1,NY2,1,   1,1,LSTEP)
IF(LTURB) THEN
+ COVAL(WN1,  KE,1.0,0.0); COVAL(WN2,  KE,1.0,0.0)
+ COVAL(WN1,LTLS,1.0,0.0); COVAL(WN2,LTLS,1.0,0.0)
+ COVAL( WL,  KE,1.0,0.0); COVAL( WL,LTLS,1.0,0.0)
+ COVAL(WN1,WC1,LOGLAW,0.0); COVAL(WN1,VC1,LOGLAW,0.0)
        + COVAL(WN1,UC1,LOGLAW,0.0)
+ COVAL(WN2,WC1,LOGLAW,0.0); COVAL(WN2,VC1,LOGLAW,0.0)
       + COVAL(WN2,UC1,LOGLAW,0.0)
+ COVAL(WL, WC1,LOGLAW,0.0); COVAL(WL, VC1,LOGLAW,0.0)
        + COVAL(WL, UC1,LOGLAW,0.0)
ELSE
+ COVAL(WN1,WC1,1.0,0.0); COVAL(WN1,VC1,1.0,0.0)
+ COVAL(WN1,UC1,1.0,0.0)
+ COVAL(WN2,WC1,1.0,0.0); COVAL(WN2,VC1,1.0,0.0)
+ COVAL(WN2,UC1,1.0,0.0)
+ COVAL(WL, WC1,1.0,0.0); COVAL(WL, VC1,1.0,0.0)
+ COVAL(WL, UC1,1.0,0.0)
ENDIF
    GROUP 15. Termination of sweeps
LSWEEP = 200; TSTSWP = -1
    GROUP 16. Termination of iterations
SELREF = T;   RESFAC = 1.E-3
    GROUP 17. Under-relaxation devices
DTHYD = HCHAN/(NY1*WIN)
RELAX(P1, LINRLX,  1.0); RELAX(UC1,FALSDT,DTHYD)
RELAX(VC1,FALSDT,DTHYD); RELAX(WC1,FALSDT,DTHYD)
IF(LTURB) THEN
+ KELIN= 1; RELAX(KE,FALSDT,DTHYD); RELAX(EP,FALSDT,DTHYD)
ENDIF
    GROUP 19. Data communicated by satellite to GROUND
    * LSG6= T activates special treatment of swirling
              velocity component (UC1) for NX=1 cases.
LSG6= T
    GROUP 22. Spot-value print-out
IXMON = NZ/2+1; IYMON = NY/2+1; IXMON=1