TALK=F;RUN( 1, 1)

  PHOTON USE
  P
 
 
   0.20443E+04 0.15633E+04 CR
  gr ou z 1;use patgeo;msg vectors;vec z 1 sh
  msg press  and then  to end
  pause
  ENDUSE
 
MESG(
MESG(
MESG(
MESG(     2D TURBULENT BACKWARD FACING STEP
MESG(
MESG(
MESG( This case illustrates the speed-up that MIGAL 
MESG( produces for turbulent flows with differents
MESG( convection schemes, including high order schemes.
MESG( It illustrates also the benefit one can get when
MESG( using it for the KE and EP variables. From library
MESG( case N105.
MESG(
MESG(
MESG( Enter required scheme ID :  1- QUICK for momentum (default)
MESG(                                KOREN for k and eps 
MESG(                             2- HYBRID for all variables
MESG(
INTEGER(SCH)
READVDU(SCH,INT,1)
INTEGER(SLV)
MESG( Enter required solver ID :  1 - MIGAL (default)
MESG(                             2 - SIMPLEST 
MESG(
READVDU(SLV,INT,1)

IF(SLV.EQ.1) THEN
+ TEXT(Turbulent Backward Facing Step - MIGAL
ENDIF
IF(SLV.EQ.2) THEN
+ TEXT(Turbulent Backward Facing Step -SIMPLEST
ENDIF

REAL(HIN,HSTEP,LEIN,LEUP,LERC,LEDN,RENO,VIN,KEIN,EPIN)
REAL(RLXFAC);INTEGER(NX1,NX2,NX3,NX4,NY1,NY2,NY3,NXS);CHAR(SCHM)
  ** All dimensions are based on: kg, m, sec
HIN=0.762;HSTEP=0.076;RENO=50000.0;ENUL=1.6E-5
VIN=RENO*ENUL/HSTEP
KEIN=0.5*1.E-5*VIN*VIN;EPIN=KEIN**1.5*0.1643/(0.09*HIN)
    GROUP 2. Transience; time-step specification
    GROUP 3. X-direction grid specification
NX1=10;NX2=10;NX3=30;NX4=15;NXS=NX1+NX2
LEIN=7.0*HSTEP;LEUP=5.0*HSTEP;LERC=10.0*HSTEP;LEDN=30.0*HSTEP
NREGX=4;IREGX=1;GRDPWR(X,NX1,LEIN,1.0)
IREGX=2;GRDPWR(X,NX2,LEUP,-1.3);IREGX=3;GRDPWR(X,NX3,LERC,1.2)
IREGX=4;GRDPWR(X,NX4,LEDN-LERC,1.4)
    GROUP 4. Y-direction grid specification
NY1=10;NY2=10;NY3=10
NREGY=3;IREGY=1;GRDPWR(Y,-NY1,HSTEP,1.3)
IREGY=2;GRDPWR(Y,NY2,HSTEP,1.2)
IREGY=3;GRDPWR(Y,NY3,HIN-HSTEP,1.5)
    GROUP 5. Z-direction grid specification
NZ=1;ZWLAST=0.01
    GROUP 6. Body-fitted coordinates or grid distortion
    GROUP 7. Variables stored, solved & named
SOLVE(P1,U1,V1)
SOLUTN(P1,Y,Y,Y,N,N,N);TURMOD(KEMODL);STORE(ENUT)
  ** use arithmetic averaging for vels as V2.1.3 & V2.2
     default of harmonic averaging produces unphysical U1
     velocities in near-wall cells of recirculation zone.
SOLUTN(U1,P,P,P,P,P,N);SOLUTN(V1,P,P,P,P,P,N)
  ** conjugate-gradient solver
CSG3=CNGR;RLXFAC=XULAST/(NX*VIN)
    GROUP 8. Terms (in differential equations) & devices
BOOLEAN(ISHYB)
CASE SCH
WHEN 1
+ SCHEME(QUICK,U1,V1)
+ SCHEME(KOREN,KE,EP)
+ RLXFAC=0.1*XULAST/(NX*VIN)
+ ISHYB=F
WHEN 2
+ DIFCUT=0.5
+ RLXFAC=XULAST/(NX*VIN)
+ ISHYB=T
ENDCASE
    GROUP 9. Properties of the medium (or media)
RHO1=1.2
    GROUP 10. Inter-phase-transfer processes and properties
    GROUP 11. Initialization of variable or porosity fields
FIINIT(U1)=VIN;FIINIT(P1)=1.3E-4;FIINIT(KE)=KEIN;FIINIT(EP)=EPIN
     ** Initialization of variables in blocked region
  ** since using arithmetic averaging on vels, it is necessary
     to intoduce wall friction patches around blockage separately,
     remove - from patch limits.
   CONPOR(STEP,0.0,CELL,#1,-#2,#1,-#1,#1,#1)
CONPOR(STEP,0.0,CELL,#1,#2,#1,#1,#1,#1)
    GROUP 12. Convection and diffusion adjustments
    GROUP 13. Boundary conditions and special sources
PATCH(IN1,WEST,1,1,NY1+1,NY,1,NZ,1,LSTEP)
COVAL(IN1,P1,FIXFLU,RHO1*VIN)
COVAL(IN1,U1,ONLYMS,VIN);COVAL(IN1,V1,ONLYMS,0.0)
COVAL(IN1,KE,ONLYMS,KEIN);COVAL(IN1,EP,ONLYMS,EPIN)
  ** Wall boundary conditions
PATCH(TWALL,NWALL,1,NX,NY,NY,1,NZ,1,LSTEP)
COVAL(TWALL,U1,LOGLAW,0.0)
COVAL(TWALL,KE,LOGLAW,LOGLAW);COVAL(TWALL,EP,LOGLAW,LOGLAW)
PATCH(BWALL,SWALL,#3,NX,1,1,1,NZ,1,LSTEP)
COVAL(BWALL,U1,LOGLAW,0.0)
COVAL(BWALL,KE,LOGLAW,LOGLAW);COVAL(BWALL,EP,LOGLAW,LOGLAW)
  ** Outlet boundary
PATCH(OUT,EAST,NX,NX,1,NY,1,NZ,1,LSTEP)
COVAL(OUT,P1,10.0,0.0);COVAL(OUT,KE,ONLYMS,SAME)
COVAL(OUT,EP,ONLYMS,SAME)
  ** Step wall boundary conditions
     use explicit patches & not CONPOR - arguments because
     of need to use arithmetic averaging on velocities
PATCH(STEP-WW,WWALL,NXS+1,NXS+1,1,NY1,1,1,1,1)
COVAL(STEP-WW,V1,LOGLAW,0.0);COVAL(STEP-WW,KE,LOGLAW,LOGLAW)
COVAL(STEP-WW,EP,LOGLAW,LOGLAW)
 
PATCH(STEP-SW,SWALL,1,NXS,NY1+1,NY1+1,1,1,1,1)
COVAL(STEP-SW,U1,LOGLAW,0.0);COVAL(STEP-SW,KE,LOGLAW,LOGLAW)
COVAL(STEP-SW,EP,LOGLAW,LOGLAW)
    GROUP 14. Downstream pressure for PARAB=.TRUE.
    GROUP 15. Termination of sweeps
IF(ISHYB) THEN
+ LSWEEP=600
ELSE
+ LSWEEP=2000
ENDIF
    GROUP 16. Termination of iterations
LITER(U1)=10;LITER(V1)=10;LITER(KE)=10;LITER(EP)=10
    GROUP 17. Under-relaxation devices
RELAX(P1,LINRLX,0.8);RELAX(ENUT,LINRLX,0.6)
RELAX(U1,FALSDT,RLXFAC);RELAX(V1,FALSDT,RLXFAC)
RELAX(KE,FALSDT,RLXFAC);RELAX(EP,FALSDT,RLXFAC)
    GROUP 18. Limits on variables or increments to them
    GROUP 22. Spot-value print-out
TSTSWP = -1;IYMON=NY1/2;IXMON=NX1+NX2+NX3/2
    GROUP 23. Field print-out and plot control
NPLT=1;ITABL=3

resfac=1.E-10

  PARAMETERS FOR MIGAL
  --------------------

IF(SLV.EQ.1) THEN

+ lsweep=120

   1- Enlarge false time steps

+ relax(u1,falsdt,1.E+10) 
+ relax(v1,falsdt,1.E+10) 
+ relax(ke,falsdt,rlxfac*1000) 
+ relax(ep,falsdt,rlxfac*1000)

   2- Don't underrelax enut

+ RELAX(ENUT,LINRLX,1.)

   3- Switch off multi-grid because too sharp variations for KE-EP in this case

+ spedat(MIGAL,SOLVED1, c, HYDRO)
+ spedat(MIGAL,NBGRID1, i, 1)

   4- Solve KE and EP with MIGAL

+ spedat(MIGAL,SOLVED2, c, KE)
+ spedat(MIGAL,SOLVED3, c, EP)

ENDIF