** MB-FGE Test: Laminar/Turbulent swirling flow in a pipe with
                  linear expansion.
  **************************************************************
  DISPLAY
  ----------------------------------------------------------
   This  case  concerns  axisymmetric incompressible laminar
   (or  turbulent  LTURB=  T)  flow  in the pipe with linear
   expansion.   The  flow  exhibits  recirculation  in   the
   region on  the pipe  axis. That  kind of  flows is common
   for gas burners.
 
   The problem is  set in Y-Z  plane; the swirling  velocity
   component UC1 is solved as scalar using special treatment
   (LSG6= T). User can model the same problem as 3D  (NX>1),
   however the speed of convergence, provided the problem is
   truly axisymmetric, is rather slow.
   ---------------------------------------------------------
  ENDDIS
L(PAUSE
  **************************************************************
BOOLEAN(LUNIF,LTURB);  LTURB= F;  LUNIF= LTURB
INTEGER(IFC);  IFC= 1
  **************************************************************
  PHOTON USE
   p ; ; ; ; ;
 
   msg Computational Domain:
   mgrid 1 i 1
   mgrid 2 i 1 col 4
   msg Press Any Key to Continue...
   pause
   cl
   set vec av off
   msg Velocity Vectors:
   mvec 1 i 1 sh
   mvec 2 i 1 sh
   msg Press Any Key to Continue...
   pause
   cl
   msg Contours of Pressure:
   mcon 1 p1 i 1 fi
   0.005
   mcon 2 p1 i 1 fi
   0.005
   msg Press Any Key to Continue...
   pause
   cl
   msg Contours of Swirling Velocity:
   mcon 1 uc1 i 1 fi
   0.005
   mcon 2 uc1 i 1 fi
   0.005
   msg Press E  to exit PHOTON ...
   pause
  ENDUSE
  **************************************************************
INTEGER(NX1,NY1,NZ1,NZ11,NZ12,NZ13,NX2,NY2,NZ2)
REAL(REYNO,WIN,HCH1,LST1,LST2,HCH2,LCHAN,WCR,DTHYD,DXX,UIN,UCR)
REAL(KEIN,EPIN,HCHN)
    GROUP 1. Run title and other preliminaries
  ** Problem definition:
IF(LTURB) THEN
+ WIN= 1.0;  UIN= 1.0*WIN;   LCHAN= 10.; NZ13= 24
+ TEXT(MBFGE: Swirling flow; Lin-expan. (K-E).
ELSE
+ WIN= 1.0;  UIN= 1.75*WIN;  LCHAN= 5.0; NZ13= 12
+ TEXT(MBFGE: Swirling flow; Lin-expan. (Re=200).
ENDIF
HCH2= 0.5;  HCH1= 0.2;  LST1= 0.2;  LST2= 0.8;  DXX= 0.025
HCHN= 0.6;  NZ11= 4;    NZ12= 10;   NZ1 = NZ11+NZ12+NZ13
NY1 = 12;   NX1 = 1;    NX2 = NX1;  NY2 = 5;    NZ2= IFC*NZ13
    GROUP 6. Body-fitted coordinates or grid distortion
BFC = T; GSET(D,NX1,NY1,NZ1,DXX,HCH2,LCHAN)
    ** Create first domain:
GSET(P,P1,0.0,0.0 ,0.0      );  GSET(P,P2,0.0,0.0 ,LST1     )
GSET(P,P3,0.0,0.0 ,LST1+LST2);  GSET(P,P4,0.0,0.0 ,LCHAN    )
GSET(P,P5,0.0,HCH2,LCHAN    );  GSET(P,P6,0.0,HCH2,LST1+LST2)
GSET(P,P7,0.0,HCH1,LST1     );  GSET(P,P8,0.0,HCH1,0.0      )
GSET(L,L12,P1,P2,NZ11, 1.0);  GSET(L,L23,P2,P3,NZ12, 1.3)
GSET(L,L34,P3,P4,NZ13, 1.4);  GSET(L,L45,P4,P5,NY1,  1.0)
GSET(L,L56,P5,P6,NZ13,-1.4);  GSET(L,L67,P6,P7,NZ12,-1.3)
GSET(L,L78,P7,P8,NZ11, 1.0);  GSET(L,L81,P8,P1,NY1,  1.0)
GSET(F,F1,P1,P2.P3,P4,-,P5,P6.P7,P8,-); GSET(M,F1,+K+J,1,1,1)
GSET(C,I:NX1+1:,F,I1,RZ,-DXX*NX1,0.0,0.0,INC,1.0)
DUMPC(MBGR1)
    ** Create second domain:
GSET(D,NX2,NY2,NZ2,DXX,HCHN,LCHAN)
GSET(P,P1,0.0,HCH2,LST1+LST2);  GSET(P,P2,0.0,HCH2,    LCHAN)
GSET(P,P3,0.0,HCHN,    LCHAN);  GSET(P,P4,0.0,HCHN,LST1+LST2)
GSET(L,L12,P1,P2,NZ2, 1.4);  GSET(L,L23,P2,P3,NY2, 1.0)
GSET(L,L34,P3,P4,NZ2,-1.4);  GSET(L,L41,P4,P1,NY2, 1.0)
GSET(F,F2,P1,-,P2,-,P3,-,P4,-); GSET(M,F2,+K+J,1,1,1)
GSET(C,I:NX2+1:,F,I1,RZ,-DXX*NX2,0.0,0.0,INC,1.0)
DUMPC(MBGR2)
  ** Assemble blocks:
NUMBLK = 2;  READCO(MBGR+YL)
    ** To set up LINKS you can use MBLINK or MPATCH commands
       instead of READCO(...+YL).
    MBLINK(1,NORTH,2,SOUTH)
    GROUP 7. Variables stored, solved & named
STORE(VPOR); SOLVE(P1,V1,W1); SOLUTN(P1,Y,Y,Y,N,N,N)
NAME(C1)= UC1;  NAME(C2)= VC1;  NAME(C3)= WC1
SOLVE(UC1,VC1,WC1);      SOLUTN(UC1,Y,Y,Y,P,P,P)
SOLUTN(VC1,Y,Y,Y,P,P,P); SOLUTN(WC1,Y,Y,Y,P,P,P)
IF(NX.GT.1) THEN
+ XCYCLE= T;  SOLVE(U1)
+ TERMS(UC1,N,Y,Y,P,P,P); TERMS(U1,N,N,N,N,P,N)
ENDIF
IF(LTURB) THEN
+ TURMOD(KEMODL)
+ SOLUTN(KE,Y,Y,Y,P,P,P); SOLUTN(EP,Y,Y,Y,P,P,P)
+ KEIN= (0.05*WIN)**2;  EPIN= 0.1643*KEIN**1.5/(0.09*HCH1)
+ ENUL= 1.544E-5;  FIINIT(KE)= KEIN;  FIINIT(EP)= EPIN
ELSE
+ REYNO= 200.0;  ENUL= WIN*HCH1/REYNO
ENDIF
    GROUP 8. Terms (in differential equations) & devices
TERMS(WC1,N,Y,Y,P,P,P);  TERMS(VC1,N,Y,Y,P,P,P)
TERMS(W1, N,N,N,N,P,N);  TERMS(V1, N,N,N,N,P,N)
    GROUP 9. Properties of the medium (or media)
RHO1= 1.189
    GROUP 11. Initialization of variable or porosity fields
INIADD= F; FIINIT(WC1)= WIN
    GROUP 13. Boundary conditions and special sources
    ** Inlet.
DO II = 1,NY1
+IF(LUNIF) THEN
+ WCR = WIN
+ELSE
+ WCR = WIN*(1.0-((2*II-1)/NY1/2)**2)
+ENDIF
+ UCR = UIN*(2*II-1)/NY1/2
+ MPATCH(1,INL:II:,LOW,1,NX1,II,II,1,1,1,LSTEP)
+  COVAL(INL:II:,P1,FIXFLU,WCR*RHO1)
+  COVAL(INL:II:,UC1,ONLYMS,UCR); COVAL(INL:II:,VC1,ONLYMS,0.0)
+  COVAL(INL:II:,WC1,ONLYMS,WCR)
+IF(LTURB) THEN
+  COVAL(INL:II:,KE,ONLYMS,KEIN); COVAL(INL:II:,EP,ONLYMS,EPIN)
+ENDIF
ENDDO
    ** Walls.
MPATCH(1,WN1,NWALL,1,NX1,NY1,NY1,1,NZ11+NZ12,1,LSTEP)
MPATCH(2,WL2,LWALL,1,NX2,  1,NY2,1,        1,1,LSTEP)
MPATCH(2,WN2,NWALL,1,NX2,NY2,NY2,1,      NZ2,1,LSTEP)
    ** Outlet.
MPATCH(1,OUT1,HIGH,1,NX1,1,NY1,NZ1,NZ1,1,LSTEP)
MPATCH(2,OUT2,HIGH,1,NX2,1,NY2,NZ2,NZ2,1,LSTEP)
 COVAL(OUT1, P1,  1.E5, 0.0); COVAL(OUT2, P1,  1.E5, 0.0)
 COVAL(OUT1,UC1,ONLYMS, 0.0); COVAL(OUT2,UC1,ONLYMS, 0.0)
 COVAL(OUT1,VC1,ONLYMS,SAME); COVAL(OUT2,VC1,ONLYMS,SAME)
 COVAL(OUT1,WC1,ONLYMS,SAME); COVAL(OUT2,WC1,ONLYMS,SAME)
IF(LTURB) THEN
+ COVAL( WN1,KE, LOGLAW,LOGLAW); COVAL( WN1,EP, LOGLAW,LOGLAW)
+ COVAL( WL2,KE, LOGLAW,LOGLAW); COVAL( WL2,EP, LOGLAW,LOGLAW)
+ COVAL( WN2,KE, LOGLAW,LOGLAW); COVAL( WN2,EP, LOGLAW,LOGLAW)
+ COVAL(OUT1,KE,ONLYMS, SAME); COVAL(OUT1,EP,ONLYMS, SAME)
+ COVAL(OUT2,KE,ONLYMS, SAME); COVAL(OUT2,EP,ONLYMS, SAME)
+ COVAL(WN1,UC1,LOGLAW,0.0); COVAL(WN1,VC1,LOGLAW,0.0)
+ COVAL(WN1,WC1,LOGLAW,0.0)
+ COVAL(WL2,UC1,LOGLAW,0.0); COVAL(WL2,VC1,LOGLAW,0.0)
+ COVAL(WL2,WC1,LOGLAW,0.0)
+ COVAL(WN2,UC1,LOGLAW,0.0); COVAL(WN2,VC1,LOGLAW,0.0)
+ COVAL(WN2,WC1,LOGLAW,0.0)
ELSE
+ COVAL(WN1,UC1,1.0,0.0); COVAL(WN1,VC1,1.0,0.0)
+ COVAL(WN1,WC1,1.0,0.0)
+ COVAL(WL2,UC1,1.0,0.0); COVAL(WL2,VC1,1.0,0.0)
+ COVAL(WL2,WC1,1.0,0.0)
+ COVAL(WN2,UC1,1.0,0.0); COVAL(WN2,VC1,1.0,0.0)
+ COVAL(WN2,WC1,1.0,0.0)
ENDIF
    GROUP 15. Termination of sweeps
LSWEEP = 500; TSTSWP = -1
    GROUP 16. Termination of iterations
SELREF = T;   RESFAC = 1.E-3
    GROUP 17. Under-relaxation devices
RELAX(P1,LINRLX,0.5);  DTHYD = LCHAN/NZ/WIN
RELAX(UC1,FALSDT,DTHYD); RELAX(VC1,FALSDT,DTHYD)
RELAX(WC1,FALSDT,DTHYD); VARMIN(UC1)= 0.0
IF(LTURB) THEN
+ LSWEEP= 1000; KELIN= 1; YPLS= T
+ RELAX(KE,FALSDT,DTHYD); RELAX(EP,FALSDT,DTHYD)
ENDIF
    GROUP 19. Data communicated by satellite to GROUND
CSG3= LCRU;  LSG4= T;  LSG6= T
    GROUP 22. Spot-value print-out
IXMON= 1; IYMON= NY/2+1; IZMON= NZ/2+1