** MB-FGE Test: Axisymmetric turnaround duct.
  **************************************************************
  DISPLAY
   This  case  concerns  steady,  axisymmetric,  isothermal,
   incompressible and laminar or fully turbulent flow in the
   turnaround duct. The geometry of the case is the same, as
   was used in the 1st test case for the WUA-CFD meeting  in
   Basel (1994).
 
   The main objective of this library case is to demonstrate
   the use  of MB-FGE  technique to  calculate flows  in the
   curved ducts. For this reason default grid sizes are  not
   enough for accurate assesement of the flow.
 
   Calculations  of  the  turbulent  flow (LTURB=T) could be
   made with low-Re extensions of whether standard K-E model
   (LOWKE=T) or 2-layer K-E model (LOWKE=F).
 
   This Q1 file includes Photon Use information.
   ---------------------------------------------------------
  ENDDIS
L(PAUSE
  **************************************************************
BOOLEAN(LTURB,LOWKE); LTURB= F; LOWKE= F
  **************************************************************
  PHOTON USE
   p ; ; ; ; ;
 
   msg Computational Domain:
   mgrid 1 k 1
   mgrid 2 k 1 col 2
   mgrid 3 k 1 col 4
   msg Press Any Key to Continue...
   pause
   cl
   set vec av off
   msg Velocity Vectors:
   mvec 1 k 1 sh
   mvec 2 k 1 sh
   mvec 3 k 1 sh
   msg Press Any Key to Continue...
   pause
   cl
   msg Contours of W-velocity Resolutes:
   mcon 1 uc1 k 1 fi
   0.0001
   mcon 2 uc1 k 1 fi
   0.0001
   mcon 3 uc1 k 1 fi
   0.0001
   msg Press Any Key to Continue...
   pause
   cl
   msg Contours of Pressure:
   mcon 1 p1 k 1 fi
   0.0001
   mcon 2 p1 k 1 fi
   0.0001
   mcon 3 p1 k 1 fi
   0.0001
   msg Press E  to exit PHOTON ...
   pause
  ENDUSE
  **************************************************************
    GROUP 1. Run title and other preliminaries
IF(LTURB) THEN
+ TEXT(MB-FGE: Case 1 (WUA-CFD meeting in Basel,1994).
ELSE
+ TEXT(MB-FGE: Case 1 (WUA-CFD-94, Laminar).
ENDIF
INTEGER(NX1,NX11,NX12,NY1,NZ1,NX2,NX21,NX22,NX23,NY2,NZ2)
INTEGER(NX24,NX25,IC,NY21,NY22,NX26,IC1,JC1,IC2,JC2)
INTEGER(NX3,NY3,NZ3,IFCX,IFCY)
REAL(PI,RAD,AA,BB,RR1,CC,DD,FF,H3,ALF1,X3,Y3,XTMP,YTMP,ATMP)
REAL(HH1,HH2,RR2,RR4,ALF2,ALF3,ALF4,RR6,ALF5,ALF6,X4,Y4,X5,Y5)
REAL(X6,Y6,EE,HH6,RR7,ALF7,HH4,HH5,DHC,XE1,YE1,XE2,YE2)
REAL(DTURB,REYNO,UIN,FRIC,LMIX,TKEIN,EPSIN,DTHYD)
PI  = 3.1415;        RAD = 180./PI
AA  = 0.1158;        BB  = 0.1095;           RR1= 0.0063
ALF1= (207.5-180.)/RAD;   ALF2= 27.56/RAD;  ALF3= 38.56/RAD
X3  = AA-RR1*SIN(ALF1); Y3 = RR1*COS(ALF1);  H3 = 0.045
CC  = 0.0984;        DD  = 0.095;            FF = 0.1094
HH1 = 0.0433;        HH2 = 0.0445;           RR2= 0.0412
RR4 = 0.0405;        ALF4= 65.95/RAD;        RR6= 0.044
ALF5= 71.55/RAD;     ALF6= 108./RAD;        ALF7= 110./RAD
EE  = 0.111;         HH6 = 0.076;            RR7= 0.025
HH4 = 0.0511;        HH5 = 0.0572;           DHC= (HH6-HH5)/3.
NZ1 = 1;  NY1 = 6;   NX11= 14; NX12= 6;      NX1= NX11+NX12
NZ2 = 1;  NY21= 10;  NY22= 6;  NY2 = NY21+NY22;
NX21= 6;  NX22= 6;   NX23= 8;  NX24= 7
NX25= 6;  NX26= 14
NX2 = NX11+NX12+NX21+NX22+NX23+NX24+NX25+NX26
IFCX= 2;  IFCY= 2;   NZ3 = 1;  NY3 = IFCY*NY22
NX3= IFCX*(NX12+NX21+NX22+NX23+NX24+NX25)
  ** Initial conditions:
IF(LTURB) THEN
+ DTURB= 0.09;   REYNO= 286000.;   UIN= 1.0
+ RHO1 = 317.8;  ENUL = 1.E-4/RHO1
ELSE
+ REYNO= 1000.;  UIN = 1.0
+ RHO1 = 317.8;  ENUL= UIN*(Y3+H3)/REYNO
ENDIF
    GROUP 6. Body-fitted coordinates or grid distortion
BFC = T
  ** Define grid points and lines for the first domain:
GSET(P,P1,0.0,0.0,0.0);   GSET(P,P2, CC,0.0,0.0)
GSET(P,P3, CC, Y3,0.0);   GSET(P,P4,0.0, Y3,0.0)
GSET(L,L12,P1,P2,NX11,-1.4); GSET(L,L23,P2,P3,NY1,1.0)
GSET(L,L34,P3,P4,NX11, 1.4); GSET(L,L41,P4,P1,NY1,1.0)
GSET(P,P5,CC,0.0,0.0);   GSET(P,P6,BB,0.0,0.0)
GSET(P,P7,X3, Y3,0.0);   GSET(P,P8,CC, Y3,0.0)
GSET(L,L56,P5,P6,NX12,-1.3)
ATMP= 2.*ALF1;  YTMP= RR1*COS(ATMP);  XTMP= AA-RR1*SIN(ATMP)
GSET(L,L67,P6,P7,NY1,1.0,ARC,XTMP,YTMP,0.0)
GSET(L,L78,P7,P8,NX12, 1.2); GSET(L,L85,P8,P5,NY1,1.0)
  ** Create grid for the first domain.
GSET(D,NX1,NY1,NZ1,AA,RR1,RR1)
GSET(F,F1,P1,-,P2,-,P3,-,P4,-); GSET(M,F1,+I+J,     1,1,1)
GSET(F,F2,P5,-,P6,-,P7,-,P8,-); GSET(M,F2,+I+J,NX11+1,1,1)
GSET(C,K:NZ1+1:,F,K1,1,NX1,1,NY1,RX,0.1,0.0,0.0,INC,1.0)
DUMPC(WUA1)
  ** Define grid points and lines for the second domain:
GSET(P,P1, 0.0,    Y3,0.0);   GSET(P,P2, CC,    Y3, 0.0)
GSET(P,P23, CC,H3-DHC,0.0);   GSET(P,P3, CC,    H3, 0.0)
GSET(P,P4, 0.0,    H3,0.0);   GSET(P,P41,0.0,H3-DHC,0.0)
GSET(L, L12, P1, P2,NX11,-1.4); GSET(L,L223, P2,P23,NY21, 1.2)
GSET(L,L233,P23, P3,NY22, 1.0); GSET(L, L34, P3, P4,NX11, 1.4)
GSET(L,L441, P4,P41,NY22, 1.0); GSET(L,L411,P41, P1,NY21,-1.2)
GSET(P,P5,CC,    Y3,0.0);    GSET(P,P6,X3,    Y3,0.0)
GSET(P,S7,FF,H3-DHC,0.0);    GSET(P,S8,CC,H3-DHC,0.0)
GSET(P,P7,FF,    H3,0.0);    GSET(P,P8,CC,    H3,0.0)
GSET(L,L56,P5,P6,NX12,1.0); GSET(L,L67,P6,S7,NY21, 1.2)
GSET(L,LS7,S7,P7,NY22,1.0); GSET(L,L78,P7,P8,NX12, 1.0)
GSET(L,LS8,P8,S8,NY22,1.0); GSET(L,L85,S8,P5,NY21,-1.2)
YTMP= HH2-RR2*COS(ALF2);  XTMP= CC+RR2*SIN(ALF2)
Y4  = HH2-RR2*COS(ALF3);  X4  = CC+RR2*SIN(ALF3)
GSET(P,P9,X3,Y3,0.0);   GSET(P,P10,X4,Y4,0.0)
GSET(L,L910,P9,P10,NX21,1.0,ARC,XTMP,YTMP,0.0)
ATMP= (ALF3+ALF4)/2.
YTMP= HH2-RR4*COS(ATMP);  XTMP= CC+RR4*SIN(ATMP)
Y5  = HH2-RR4*COS(ALF4);  X5  = CC+RR4*SIN(ALF4)
GSET(P,P11,X5,Y5,0.0)
GSET(L,L111,P10,P11,NX22,1.0,ARC,XTMP,YTMP,0.0)
ATMP= (ALF5+ALF6)/2.
YTMP= HH1-RR6*COS(ATMP);  XTMP= DD+RR6*SIN(ATMP)
Y6  = HH1-RR6*COS(ALF6);  X6  = DD+RR6*SIN(ALF6)
GSET(P,P12,X6,Y6,0.0)
GSET(L,L112,P11,P12,NX23,-1.2,ARC,XTMP,YTMP,0.0)
ATMP= (PI+ALF7)/2.
YTMP= HH6-RR7-RR7*COS(ATMP);  XTMP= EE+RR7*SIN(ATMP)
GSET(P,P13,EE,HH6,0.0)
GSET(L,L123,P12,P13,NX24,1.0,ARC,XTMP,YTMP,0.0)
GSET(P,P14,FF,HH5+DHC,0.0); GSET(L,L134,P13,P14,NY21,1.2)
GSET(P,P15,FF, H3-DHC,0.0);  XTMP= FF+HH4-H3+DHC
GSET(L,L145,P14,P15,NX21+NX22+NX23+NX24,1.0,ARC,XTMP,HH4,0.0)
GSET(L,L159,P15,P9,NY21,-1.2)
GSET(P,S15,FF, H3-DHC,0.0); GSET(P,S14,FF,HH5+DHC,0.0)
GSET(L,LS45,S15,S14,NX21+NX22+NX23+NX24,1.0,ARC,XTMP,HH4,0.0)
GSET(P,P16,FF,    HH5,0.0); GSET(L,L146,S14,P16,NY22,1.0)
GSET(P,P17,FF,     H3,0.0); XTMP= FF+HH4-H3
GSET(L,L167,P16,P17,NX21+NX22+NX23+NX24,1.0,ARC,XTMP,HH4,0.0)
GSET(L,L175,P17,S15,NY22,1.0)
GSET(P,P18,EE,    HH6,0.0);   GSET(P,P19,CC,    HH6,0.0)
GSET(P,S20,CC,HH5+DHC,0.0);   GSET(P,S21,FF,HH5+DHC,0.0)
GSET(P,P20,CC,    HH5,0.0);   GSET(P,P21,FF,    HH5,0.0)
GSET(L,L189,P18,P19,NX25, 1.0); GSET(L,L190,P19,S20,NY21, 1.2)
GSET(L,LS00,S20,P20,NY22, 1.0); GSET(L,L201,P20,P21,NX25, 1.0)
GSET(L,LS11,P21,S21,NY22, 1.0); GSET(L,L218,S21,P18,NY21,-1.2)
GSET(P,T22, CC,    HH6,0.0);   GSET(P,T23,0.0,    HH6,0.0)
GSET(P,S24,0.0,HH5+DHC,0.0);   GSET(P,S25, CC,HH5+DHC,0.0)
GSET(P,T24,0.0,    HH5,0.0);   GSET(P,T25, CC,    HH5,0.0)
GSET(L,M223,T22,T23,NX26, 1.4); GSET(L,M234,T23,S24,NY21, 1.2)
GSET(L,MS44,S24,T24,NY22, 1.0); GSET(L,M245,T24,T25,NX26,-1.4)
GSET(L,MS55,T25,S25,NY22, 1.0); GSET(L,M252,S25,T22,NY21,-1.2)
  ** Create grid for the second domain.
GSET(D,NX2,NY2,NZ2,X3,H3,RR1)
GSET(F,F1,P1,-,P2,P23,P3,-,P4,P41); GSET(M,F1,+I+J,     1,1,1)
GSET(F,F2,P5,-,P6, S7,P7,-,P8, S8); GSET(M,F2,+I+J,NX11+1,1,1)
GSET(F,F3,P9,P10.P11.P12,P13,-,P14,-,P15,-)
GSET(M,F3, +I+J,NX11+NX12+1,1,1)
GSET(F,FS3,S15,-,S14,-,P16,-,P17,-)
GSET(M,FS3,+I+J,NX11+NX12+1,NY21+1,1)
IC = NX11+NX12+NX21+NX22+NX23+NX24
GSET(F,F4,P18,-,P19,S20,P20,-,P21,S21); GSET(M,F4,+I+J,IC+1,1,1)
IC = NX11+NX12+NX21+NX22+NX23+NX24+NX25
GSET(F,F5,T22,-,T23,S24,T24,-,T25,S25); GSET(M,F5,+I+J,IC+1,1,1)
  ** Define points for the embedded grid:
JC1= NY21+1;  IC1= NX11+1;  JC2= NY+1;  IC2= IC1+NX12
XE1= XC(IC1,JC1,1);  YE1= YC(IC1,JC1,1)
XE2= XC(IC2,JC1,1);  YE2= YC(IC2,JC1,1)
GSET(P,E1,XE1,YE1,0.0);   GSET(P,E2,XE2,YE2,0.0)
XE1= XC(IC1,JC2,1);  YE1= YC(IC1,JC2,1)
XE2= XC(IC2,JC2,1);  YE2= YC(IC2,JC2,1)
GSET(P,E8,XE1,YE1,0.0);   GSET(P,E7,XE2,YE2,0.0)
IC1= IC1+NX12+NX21+NX22+NX23+NX24;  IC2= IC1+NX25
XE1= XC(IC1,JC1,1);  YE1= YC(IC1,JC1,1)
XE2= XC(IC2,JC1,1);  YE2= YC(IC2,JC1,1)
GSET(P,E3,XE1,YE1,0.0);   GSET(P,E4,XE2,YE2,0.0)
XE1= XC(IC1,JC2,1); YE1= YC(IC1,JC2,1)
XE2= XC(IC2,JC2,1); YE2= YC(IC2,JC2,1)
GSET(P,E6,XE1,YE1,0.0);   GSET(P,E5,XE2,YE2,0.0)
GSET(C,K:NZ2+1:,F,K1,1,NX2,1,NY2,RX,0.1,0.0,0.0,INC,1.0)
DUMPC(WUA2)
  ** Define grid points and lines for the embedded grid:
GSET(D,NX3,NY3,NZ3,X3,H3,RR1)
GSET(L,LE12,E1,E2,IFCX*NX12,1.0); IC1= IFCX*(NX21+NX22+NX23+NX24)
XTMP= FF+HH4-H3+DHC; GSET(L,LE23,E2,E3,IC1,1.0,ARC,XTMP,HH4,0.0)
GSET(L,LE34,E3,E4,IFCX*NX25,1.0); GSET(L,LE45,E4,E5,IFCY*NY22,1.0)
GSET(L,LE56,E5,E6,IFCX*NX25,1.0)
XTMP= FF+HH4-H3; GSET(L,LE67,E6,E7,IC1,1.0,ARC,XTMP,HH4,0.0)
GSET(L,LE78,E7,E8,IFCX*NX12,1.0); GSET(L,LE81,E8,E1,IFCY*NY22,1.0)
  ** Create embedded grid:
GSET(F,F1,E1,E2.E3,E4,-,E5,E6.E7,E8,-); GSET(M,F1,+I+J,1,1,1)
GSET(C,K:NZ3+1:,F,K1,1,NX3,1,NY3,RX,0.1,0.0,0.0,INC,1.0)
DUMPC(WUA3)
  ** Assemble blocks:
NUMBLK= 3; READCO(WUA+L); GVIEW(Z); VIEW
    ** Set links:
    ** To set up LINKS you can use MBLINK or MPATCH commands
       instead of READCO(...+L).
       MBLINK(1,NORTH,2,SOUTH)
       PATCH(MBL1.2,NORTH,    1,    NX1,NY1,NY1,1,NZ1,1,LSTEP)
       PATCH(MBL2.1,SOUTH,NX1+2,2*NX1+1,  1,  1,1,NZ2,1,LSTEP)
       MBLINK(3,IN,2)
    GROUP 7. Variables stored, solved & named
STORE(VPOR); SOLVE(P1,U1,V1)
IF(LTURB) THEN
+ STORE(ENUT,LEN1)
+ IF(LOWKE) THEN
+  TURMOD(KEMODL-LOWRE)
+ ELSE
+  TURMOD(KEMODL-2L)
+ ENDIF
ENDIF
L($F150)
    GROUP 11. Initialization of variable or porosity fields
IF(LTURB) THEN
+ FRIC = 0.018;       TKEIN = 0.25*UIN*UIN*FRIC
+ LMIX = 0.09*DTURB;  EPSIN = 0.1643*TKEIN**1.5/LMIX
+ FIINIT(P1)= 1.3E-4; FIINIT(KE)= TKEIN;  FIINIT(EP)= EPSIN
ENDIF
    GROUP 13. Boundary conditions and special sources
    ** Inlet.
MPATCH(1,IN1,WEST,1,1,1,NY1,1,1,NZ1,LSTEP)
MPATCH(2,IN2,WEST,1,1,1,NY2,1,1,NZ2,LSTEP)
COVAL(IN1,P1, FIXFLU,RHO1*UIN); COVAL(IN2,P1, FIXFLU,RHO1*UIN)
COVAL(IN1,UC1,ONLYMS,     UIN); COVAL(IN1,VC1,ONLYMS,     0.0)
COVAL(IN2,UC1,ONLYMS,     UIN); COVAL(IN2,VC1,ONLYMS,     0.0)
    ** Walls.
MPATCH(1,WE,EWALL,        NX1,NX1,  1,NY1,1,NZ1,1,LSTEP)
MPATCH(2,WN,NWALL,          1,NX2,NY2,NY2,1,NZ2,1,LSTEP)
MPATCH(2,WS,SWALL,NX11+NX12+1,NX2,  1,  1,1,NZ2,1,LSTEP)
    ** Outlet.
MPATCH(2,OUT,EAST,NX2,NX2,1,NY2,1,NZ2,1,LSTEP)
COVAL(OUT,P1,  1.0E5,0.0)
COVAL(OUT,UC1,ONLYMS,0.0); COVAL(OUT,VC1,ONLYMS,0.0)
IF(LTURB) THEN
+ COVAL(IN1, KE,ONLYMS,TKEIN); COVAL(IN1, EP,ONLYMS,EPSIN)
+ COVAL(IN2, KE,ONLYMS,TKEIN); COVAL(IN2, EP,ONLYMS,EPSIN)
+ COVAL(WE, UC1, LOGLAW,0.0); COVAL(WE, VC1, LOGLAW,0.0)
+ COVAL(WE,  KE,   1.0,0.0); COVAL(WE,LTLS,   1.0,0.0)
+ COVAL(WN, UC1,LOGLAW,0.0); COVAL(WN, VC1, LOGLAW,0.0)
+ COVAL(WN,  KE,   1.0,0.0); COVAL(WN,LTLS,   1.0,0.0)
+ COVAL(WS, UC1, LOGLAW,0.0); COVAL(WS, VC1, LOGLAW,0.0)
+ COVAL(WS,  KE,   1.0,0.0); COVAL(WS,LTLS,   1.0,0.0)
+ COVAL(OUT,KE, ONLYMS,0.0); COVAL(OUT,EP, ONLYMS,0.0)
ELSE
+ COVAL(WE,UC1,1.0,0.0); COVAL(WE,VC1,1.0,0.0)
+ COVAL(WN,UC1,1.0,0.0); COVAL(WN,VC1,1.0,0.0)
+ COVAL(WS,UC1,1.0,0.0); COVAL(WS,VC1,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.1);  DTHYD= 10.*DD/NX11/UIN
RELAX(UC1,FALSDT,DTHYD); RELAX(VC1,FALSDT,DTHYD)
IF(LTURB) THEN
+ KELIN = 1
+ RELAX(KE,FALSDT,DTHYD); RELAX(EP,FALSDT,DTHYD)
ENDIF
    GROUP 19. Data communicated by satellite to GROUND
    * LSG3 = T, activates curvilinearity treatment;
    * LSG4 = T, activates nonorthogonality treatment.
LSG3= BFC;  LSG4= T
    GROUP 21. Print-out of variables
IF(LTURB) THEN
+ OUTPUT(KE, Y,N,N,Y,Y,Y); OUTPUT(EP, Y,N,N,Y,Y,Y)
ENDIF
    GROUP 22. Spot-value print-out
IXMON= NX1/2+1; IYMON= NY1/2+1; IZMON= 1